Prepared in cooperation with the Bureau of Reclamation

Technical Evaluation of a Total Maximum Daily Load Model 

for Upper Klamath Lake, Oregon

Open-File Report 2013–1262

U.S. Department of the Interior

U.S. Geological Survey

Technical Evaluation of a Total Maximum Daily Load Model 

for Upper Klamath Lake, Oregon

By Tamara M. Wood, Susan A. Wherry, James L. Carter, James S. Kuwabara, Nancy S. Simon, and 

Stewart A. Rounds

Prepared in cooperation with the Bureau of Reclamation

Open-File Report 2013-1262

U.S. Department of the Interior

U.S. Geological Survey
ii

U.S. Department of the Interior

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

U.S. Geological Survey, Reston, Virginia 2013

For more information on the USGS—the Federal source for science about the Earth, 

its natural and living resources, natural hazards, and the environment:

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

Telephone: 1-888-ASK-USGS

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

http://www.usgs.gov/pubprod

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

Any use of trade, 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 material contained within this report.

Suggested citation:

Wood, T.M., Wherry, S.A., Carter, J.L., Kuwabara, J.S., Simon, N.S., and Rounds, S.A., 2013, Technical evaluation of 

a total maximum daily load model for Upper Klamath and Agency Lakes, Oregon: U.S. Geological Survey Open-File 

Report 2013–1262, 75 p, http://dx.doi.org/10.3133/ofr20131262

ISSN 2331-1258 (online)
iii

Contents

Executive Summary....................................................................................................................................................... 1

Rationale for this Review ........................................................................................................................................... 1

Testing the Model Using New Information ................................................................................................................. 2

Suggestions for Model Improvement ......................................................................................................................... 2

Introduction.................................................................................................................................................................... 4

Background................................................................................................................................................................ 4

Purpose and Scope ................................................................................................................................................... 8

Phosphorus in Upper Klamath Lake Bottom Sediments................................................................................................ 8

Phosphorus Content in Sediment Samples Collected Between 1969 and 2006 ........................................................ 8

Spatial Variability ......................................................................................................................................................12

Internal Loading of Phosphorus in Upper Klamath Lake...............................................................................................14

Mechanisms..............................................................................................................................................................14

Possible Role of Benthic Invertebrates in Internal Loading in Upper Klamath Lake..................................................17

Verification of the Upper Klamath Lake Phosphorus TMDL Model...............................................................................22

Water-Column Phosphorus Submodel......................................................................................................................26

Sediment-Phosphorus Submodel .............................................................................................................................26

Chlorophyll a Submodel............................................................................................................................................29

pH Submodel ............................................................................................................................................................30

Alternative Solution Method ......................................................................................................................................30

Model Calibration......................................................................................................................................................31

Step 1—pH Regression Model Parameters...........................................................................................................31

Step 2—Algal Settling Rate...................................................................................................................................31

Step 3—Nonalgal Phosphorus Loss Rate .............................................................................................................32

Step 4—Internal Recycling Parameters ................................................................................................................32

Comparison of Original to Recalibrated Model Performance.................................................................................33

Sensitivity of Calibration Parameters to Model Modifications........................................................................................41

Growth Limitation Factors .........................................................................................................................................41

Recycling Mechanism...............................................................................................................................................42

Numerical Solution of Chlorophyll a Equation...........................................................................................................46

Initial Sediment Reservoir .........................................................................................................................................47

Long-Term Model Simulation and Implications for TMDL .............................................................................................49

Summary, Conclusions, and Suggestions for Improvements to the TMDL Model ........................................................60

Verification of Model Algorithms................................................................................................................................60

Sensitivity of Calibration Parameters to Model Assumptions ....................................................................................60

Effect of Calibration Parameters on Long-Term Model Simulations..........................................................................61

Suggestions for Improvements to the TMDL Model ..................................................................................................63

References Cited..........................................................................................................................................................65

Appendix A. R-Script Code for Upper Klamath Lake Phosphorus Model .....................................................................71
iv

Figures

Figure 1. Map showing sampling sites for benthic flux and sediment core and grab sample collection in Upper 

Klamath Lake, Oregon................................................................................................................................................... 5

Figure 2. Graphs showing average concentration of phosphorus fractions as a function of depth in sediment from 26 

cores collected from Upper Klamath Lake, Oregon, 2006. .......................................................................................... 10

Figure 3. Graphs showing average concentration of carbon, nitrogen, and trace elements as a function of depth in 

sediment from 26 cores collected from Upper Klamath Lake, Oregon, 2006............................................................... 10

Figure 4. Diagrams showing varimax rotated principal component scores for all constituents in cores collected from 

Upper Klamath Lake, Oregon, integrated to 10 centimeter depth................................................................................ 13

Figure 5. Graph showing calculated whole-lake average of internal loading to Upper Klamath Lake, Oregon, April 

1991–September 1998. ............................................................................................................................................... 15

Figure 6. Diagram showing original total maximum daily load (TMDL) model dependence loops. ............................. 25

Figure 7. Diagrams showing simulated and observed water column (A) total phosphorus concentration, (B) 

chlorophyll a concentration, and (C) pH, Upper Klamath Lake, Oregon, April 1991–September 1998........................ 34

Figure 8. Diagrams showing simulated and observed water column (A) total phosphorus concentration, (B) 

chlorophyll a concentration, and (C) pH in Upper Klamath Lake, Oregon, April 1991–April 1998. .............................. 43

Figure 9. Graph showing mass of phosphorus recycled from sediments, calculated from observations of water 

column phosphorus and external loads, Upper Klamath Lake, Oregon, May–July 1991–1998, as a function of lake 

temperature. ................................................................................................................................................................ 44

Figure 10. Graphs showing simulated and observed water column (A) total phosphorus concentration, (B) chlorophyll 

a concentration, and (C) pH, Upper Klamath Lake, Oregon, April 1991–September 1998.......................................... 45

Figure 11. Graphs showing simulated and observed water column (A) total phosphorus concentration, (B) chlorophyll 

a concentration, and (C) pH, Upper Klamath Lake, Oregon, April 1991–September 1998.......................................... 47

Figure 12. Graphs showing simulated and observed water column (A) total phosphorus concentration, (B) chlorophyll 

a concentration, and (C) pH, Upper Klamath Lake, Oregon, April 1991–September 1998.......................................... 49

Figure 13. Graphs showing average (A) concentration of total phosphorus, (B) size of the sediment reservoir, (C) 

concentration of chlorophyll a during the final 7 years of a 203-year simulation for Upper Klamath Lake, Oregon, as a 

function of the average ratio of recycled mass to settled mass during the first 7 years of a 203-year simulation. (D) 

Total phosphorus concentration from graph A as a function of the sediment reservoir from graph B. ......................... 53

Figure 14. Graphs showing comparison of long-term simulated concentration of water column phosphorus, under 

conditions of (A) 0 and (B) 40 percent reductions in external loads to Upper Klamath Lake, Oregon. ........................ 54

Figure 15. Graphs showing comparison of long-term simulated concentration of sediment phosphorus reservoir in 

Upper Klamath Lake, Oregon, under conditions of (A) 0 and (B) 40 percent reductions in external loads. ................. 55

Figure 16. Graphs showing comparison of long-term simulated concentration of water column phosphorus in Upper 

Klamath Lake, Oregon, under conditions of (A) 0 and (B) 40 percent reductions in external loads............................. 56

Figure 17. Graphs showing comparison of long-term simulated concentration of sediment phosphorus reservoir in 

Upper Klamath Lake, Oregon, under conditions of (A) 0 and (B) 40 percent reductions in external loads.. ................ 56

Figure 18. Graphs showing simulated chlorophyll a concentration during the last 7 years of a 203-year simulation 

under conditions of (A) 0 percent and (B) 40 percent reductions in external loads to Upper Klamath Lake, Oregon. . 58

Figure 19. Graphs showing simulated total phosphorus concentration during the last 7 years of a 203-year simulation 

under conditions of (A) 0 percent and (B) 40 percent reductions in external loads to Upper Klamath Lake, Oregon. . 59
v

Tables

Table 1. Comparison of concentrations of total and inorganic phosphorus in sediment cores collected October 2006 

from Upper Klamath Lake, Oregon, to reported concentrations prior to 2006.............................................................. 11

Table 2. Loadings of sediment variables onto varimax rotated principal components for all constituents in cores 

collected from Upper Klamath Lake, Oregon, integrated to 10 centimeter depth......................................................... 13

Table 3. Comparison of mass-specific excretion of phosphorus by benthic invertebrates from Upper Klamath Lake, 

Oregon......................................................................................................................................................................... 20

Table 4. Values and sources of parameters used in the Upper Klamath Lake, Oregon, phosphorus TMDL model—

continued..................................................................................................................................................................... 23

Table 5. Definitions of symbols used in this report ...................................................................................................... 27

Table 6. Comparison of simulated model variables obtained from the original, semi-analytical and Runge-Kutta 

models used in the review of the original model used to establish a TMDL for Upper Klamath Lake, Oregon. ........... 31

Table 7. Comparison of pH model coefficients and statistics between the original TMDL model for Upper Klamath 

Lake, Oregon, and the model recalibration.................................................................................................................. 31

Table 8. Comparison of descriptive statistics between observed and calculated non-algal mass lost to sedimentation 

during March–May 1991–1998 .................................................................................................................................... 32

Table 9. Comparison of performance statistics for total phosphorus obtained over a validation period including all 

months in 1991–1993 and August–April in 1994–1998 from model versions discussed in the report. ........................ 35

Table 10. Comparison of performance statistics for chlorophyll a obtained over a validation period including all 

months in 1991–1993 and August–April in 1994–1998 from model versions discussed in the report. ....................... 36

Table 11. Comparison of performance statistics for pH obtained over a validation period including all months in 1991–

1993 and August–April in 1994–1998 from model versions discussed in the report.................................................... 37

Table 12. Comparison of performance statistics for total phosphorus over a validation period including May–July 

1991–1993 obtained from various model versions discussed in the report. ................................................................ 38

Table 13. Comparison of performance statistics for chlorophyll a obtained over a validation period including May–July 

1991–1993 from model versions discussed in the report............................................................................................. 39

Table 14. Comparison of performance statistics for pH obtained over a validation period including May–July 1991–

1993 from various model versions discussed in the report. ......................................................................................... 40

Table 15. Comparison of metrics describing the long-term simulations for model versions discussed in the report. .. 51
vi

Conversion Factors and Datums

Conversion Factors

Multiply By To obtain

Length

meter (m) 3.281 foot (ft) 

centimeter (cm) 2.54x inch (in.)

kilometer (km) 0.6214 mile (mi)

Area

2 square kilometer (km ) 247.1 acre

2 square kilometer (km ) 0.3861 square mile (mi2

)

2 square meter (m ) 10.7639 square foot (ft2)

Volume

3 cubic meter (m ) 35.31 cubic foot (ft3

)

Flow rate

meter per second (m/s) 3.281 foot per second (ft/s) 

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

cubic meter per day (m3/d) 35.31 cubic foot per day (ft3/d) 

Mass

gram (g) 0.03527 ounce avoirdupois (oz)

milligram (mg) 3.527×10-5 ounce avoirdupois (oz)

kilogram (kg) 2.205 pound avoirdupois (lb)

metric ton 2,205.00 pound avoirdupois (lb)

micrograms per liter (μg/L)

Concentration

parts per billion (ppb) 1.0

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as follows: °F=(1.8×°C)+32.

Concentrations of chemical constituents in water are given in milligrams per liter (mg/L), micrograms per liter (μg/L), or milligrams per gram dry weight (mg/g DW). 

Datums

Vertical coordinate information is referenced to the Bureau of Reclamation datum, which is 1.78 feet above National Geodetic 

Vertical Datum of 1929 (NGVD 29).
1 

Technical Evaluation of a Total Maximum Daily Load 

Model for Upper Klamath and Agency Lakes, Oregon

By Tamara M. Wood, Susan A. Wherry, James L. Carter, James S. Kuwabara, Nancy S. Simon, 

and Stewart A. Rounds

Executive Summary

Rationale for this Review

We reviewed a mass balance model developed in 2001 that guided establishment of the 

phosphorus total maximum daily load (TMDL)

for Upper Klamath and Agency Lakes, Oregon. 

The purpose of the review was to evaluate the 

strengths and weaknesses of the model and to determine whether improvements could be made 

using information derived from studies since the 

model was first developed. The new data have 

contributed to the understanding of processes in 

the lakes, particularly internal loading of phosphorus from sediment, and include measurements 

of diffusive fluxes of phosphorus from the bottom 

sediments, groundwater advection, desorption 

from iron oxides at high pH in a laboratory setting, and estimates of fluxes of phosphorus bound 

to iron and aluminum oxides. None of these processes in isolation, however, is large enough to 

account for the episodically high values of wholelake internal loading calculated from a mass balance, which can range from 10 to 20 milligrams 

per square meter per day for short periods. 

The possible role of benthic invertebrates in 

lake sediments in the internal loading of phosphorus in the lake has become apparent since the 

development of the TMDL model. Benthic invertebrates can increase diffusive fluxes several-fold 

through bioturbation and biodiffusion, and, if the 

invertebrates are bottom feeders, they can recycle 

phosphorus to the water column through metabolic excretion. These organisms have high densities 

(1,822–62,178 individuals per square meter) in 

Upper Klamath Lake. Conversion of the mean 

density of tubificid worms (Oligochaeta) and chironomid midges (Diptera), two of the dominant 

taxa, to an areal flux rate based on laboratory 

measurements of metabolic excretion of two 

abundant species suggested that excretion by benthic invertebrates is at least as important as any of 

the other identified processes for internal loading 

to the water column. 

Data from sediment cores collected around 

Upper Klamath Lake since the development of 

the TMDL model also contributed to this review. 

Cores were sequentially extracted to determine

the distribution of phosphorus associated with

several matrices in the sediment (freely exchangeable, metal-oxides, acid-soluble minerals, 

and residual). The concentrations of phosphorus 

in these fractions varied around the lake in patterns that reflect transport processes in the lake 

and the ultimate deposition of organic and inorganic forms of phosphorus from the water 

column. Both organic and inorganic phosphorus 

had higher concentrations in the northern part of 

the lake, in and just west of Goose Bay. At the 

time that these cores were collected, prior to restoration of the Williamson River Delta, this area 

was close to the shoreline of the lake and east of 

the Williamson River mouth. This contrasts with

erosional inputs, which, in addition to being high 

to the east of the pre-restoration Williamson River mouth, were higher in the middle of the lake 

than at the northern end. Organic forms of phos-
2 

phorus had particularly high concentrations in the 

northern bays. When these cores were used to 

calculate a new estimate of the whole-lakeaveraged concentration of total phosphorus in the 

top 10 centimeters of the lake sediments, the estimate was about one-third of the best estimate 

available when the TMDL model was developed. 

Testing the Model Using New Information

Verification of the original TMDL model 

and its calibration identified three readily implemented changes in model algorithms that 

improve the model: (1) integrating the light limitation function over depth, (2) replacing the 

phosphorus limitation function with a MichaelisMenten formulation, and (3) removing the assumption that chlorophyll a reaches instantaneous 

equilibrium with light and nutrients.

When these changes were made and the 

model was recalibrated, the calibration parameters associated with phosphorus recycling 

changed substantially, resulting in quantitative 

differences such as higher predicted total phosphorus and chlorophyll a concentrations in the 

long term (203-year) simulations that were used 

to predict the steady-state response of the lakes to 

the implementation of the TMDL. The improved 

model also predicted that, at steady state, chlorophyll a concentration would exceed 100 parts per 

billion (ppb, the units used in the original model 

documentation; equivalent to micrograms per liter) 

in 4 of 7 years, instead of 2 of 7 years. Qualitatively, however, the predictions of the improved 

model were consistent with those of the previous 

model in simulating a long-term reduction in total 

phosphorus and chlorophyll a concentrations 

relative to current conditions and in predicting 

that severe blooms (peak chlorophyll a concentration greater than (>) 100 ppb) would no longer 

happen every year. 

We used a sensitivity analysis to examine the 

result of changing the estimate of the initial sediment phosphorus reservoir to about one-third the 

value used in the original TMDL simulations, 

based on concentration data from the more recently collected sediment cores. Using the new 

estimate and updating the model calibration parameters caused the long-term predictions of the 

model to change because the calibrated recycling 

rates increased. The new recycling rates increased

predicted water column total phosphorus and 

chlorophyll a concentrations in the long term 

(203-year) simulations, a consequence of the result that, for a given a set of boundary conditions 

(percent reduction in external load) and initial 

conditions (size of sediment-phosphorus reservoir), the steady-state water-column phosphorus 

concentration generally was negatively correlated

with the steady-state sediment reservoir.

The model equations describe a coupled water column/sediment system in which the two 

compartments are inter-dependent. In such a system, as long as external inputs to the system 

continue, the sediment compartment will not 

trend independently toward depletion until it is 

gone. Rather, phosphorus concentrations in the 

water column will reach a steady state in equilibrium with the sediment reservoir over time and 

the calibration parameters (not the magnitude of 

the sediment phosphorus reservoir) will determine the steady-state conditions. In most 

simulations, the sediment reservoir decreased 

over time to a steady-state value less than the initial condition. In some simulations, however, the 

size of the sediment-phosphorus reservoir stayed 

mostly steady or increased slightly over time as 

equilibrium was approached. 

Suggestions for Model Improvement

We developed a temperature-dependent recycling term for the model that is based on the 

hypothesis that benthic invertebrate activity is an 

important driver of internal phosphorus loading 

and that higher temperature stimulates that activity. We substituted this term for the pH-dependent 

recycling term in the original model to investigate 

the behavior of the model when the feedback inherent in the pH-dependent recycling term was 

removed. Although the predictions of the modified model were consistent with previous results 

in simulating a long-term reduction in total phos-
3 

phorus and chlorophyll a, over the 203-year 

timeframe the differences in the two approaches 

became apparent: At steady state, the pHdependent recycling approach simulated large 

blooms (defined as greater than 100 ppb chlorophyll a) in only 2 of every 7 years, whereas the 

temperature-dependent approach simulated large 

blooms every year. Therefore, the form of the 

recycling term has a strong influence on model 

predictions. Because pH-dependent desorption of 

phosphorus from iron oxyhydroxides is only one 

of several mechanisms for the internal loading of 

phosphorus in the lake, adding more realistic 

complexity to the recycling term would improve 

the model. 

The calibrated settling velocity of algal 

phosphorus varied 51 percent over the implemented model modifications, indicating that the 

calibrated parameter did not represent a true 

physical or chemical process, but rather was

compensating for overly simplistic formulations 

of complex phenomena. The model likely, therefore, could be improved by reducing uncertainty 

in the sedimentation and recycling terms through 

a combination of field programs and modeling or 

empirical techniques that could better predict sedimentation and recycling as functions of other 

readily available environmental variables. For 

example, two explanatory variables for macroinvertebrate-dependent recycling might be water 

temperature and time elapsed since the start of 

the growing season. A field program could include the collection of density, body-size, and 

excretion stoichiometry data for benthic invertebrates over a few growing seasons. Because 

many quantities, such as recycling rate, sedimentphosphorus content, benthic invertebrate density, 

and sedimentation vary around the lake, a spatially explicit phosphorus transport model of the lake 

could be used to appropriately “scale up” the spatial variability from site-dependent estimates for 

use in a whole-lake model. The model could be 

further improved through the incorporation of a 

mortality term that would move algal phosphorus 

from a viable (nonsettling) to a detrital (settling) 

physiological state.

The final improvement would be to use a 

longer dataset (in terms of number of years) to 

calibrate and validate the model. Water-quality 

data from 1991 to 1998 were used to calibrate the 

model, and the behavior of the model over the 

long term depends on the calibration parameters 

derived from fitting the earliest part of the simulation to the available calibration data. The 

simulations from our model modifications were 

consistent in that all achieved a steady state within about 25–75 years, during which total 

phosphorus and chlorophyll a concentrations in 

the lakes were decreased, and most models simulated severe blooms (chlorophyll a > 100 ppb) 

happening less often after steady state is 

achieved. The consistency of this result can be 

tested in the context of potential changes to calibrated parameters based on a longer calibration 

dataset, when that becomes available. The performance of the calibrated model with a longer 

validation dataset also can be evaluated. 
4 

Introduction 

Background

Upper Klamath and Agency Lakes (fig. 1) 

are the primary components of a large, shallow, 

freshwater lake system in south-central Oregon

that has a surface area of 260 km2 and a mean 

depth of 2.8 m at full pool. The main tributaries

include the Williamson, Sprague, and Wood Rivers. Most water leaves the lake through Link Dam 

into the Link River, a 1.6-km reach that leads to 

the head of the Klamath River. Water from the 

lakes provides essential waterfowl and fish habitat and is used for irrigation, power generation, 

and to enhance and sustain downstream flows in 

the Klamath River. A restoration project completed in 2009 added 28 km2

 around the mouth of 

the Williamson River (shown as Tulana and 

Goose Bay in figure 1) to the Agency 

Lake/Upper Klamath Lake hydrologic system, 

and formed a single water body from the two 

lakes that had previously been connected only 

through Agency Straits (fig. 1). In this report, 

“Upper Klamath Lake” refers to the combined 

lake system, before or after restoration of the 

Williamson River Delta, as appropriate.

Water currents in the Upper Klamath Lake 

system are primarily wind driven. Prevailing 

north to northwest winds move the water in a 

generally clockwise direction around the lake in a 

broad, north-to-south flow on the shallow eastern 

side of the lake and a narrow, south to north return flow in the deepest part of the lake, a trench 

along the western shoreline with a maximum 

depth of about 15 m (Wood and others, 2008). 

Dense blooms dominated by the filamentous cyanobacterium Aphanizomenon flos-aquae occur in 

the lake every summer. Thermal stratification establishes in the afternoon, and the resulting water 

column stability allows the large, buoyant colonies of A. flos-aquae to accumulate in the upper 

water column. Because the lake is so shallow, the 

water column in most of the lake mixes daily 

when wind increases, usually in the afternoon, 

although in deeper parts of the lake stratification 

can sometimes persist for several days. Conditions at the sediment/water interface do not 

routinely become completely anoxic, although 

this can happen in association with a bloom decline in the deep trench along the western 

shoreline or other deep parts of the lake when 

mixing does not occur for several days (Wood 

and others, 2006), or in the bays, where particularly high concentrations of planktonic biomass 

can occur before senescence (Wood, 2001). 

Normal summer stratification, particularly during 

periods of low wind, does, however, result in low 

dissolved oxygen and high ammonia concentrations in the near-bottom water column during 

bloom declines (Kann and Welch, 2005). These 

episodes of poor water quality negatively affect 

fish in the lake, including shortnose sucker 

(Chasmistes brevirostris) and Lost River sucker 

(Deltistes luxatus), both endangered species (Perkins and others, 2000; Banish and others, 2009). 

During the rapid growth phase of the bloom, 

generally May–July, dissolved oxygen typically

is supersaturated in the near-surface water column, and photosynthetically elevated pH 

routinely exceeds 9.5 (Hoilman and others, 

2008). The probability of exceeding a specific pH 

value in this phase can be predicted from concentrations of the biomass surrogate chlorophyll a

(Kann and Smith, 1999). A. flos-aquae is a nitrogen-fixing cyanobacterium, and consequently, 

phosphorus is the primary limiting nutrient, at 

least during the growth phase of the first bloom 

(Lindenberg and others, 2009; Kann, 2011a), 

which is initiated and sustained by tributary 

phosphorus and by phosphorus stored in the bottom sediments that is transferred to the water 

column in the spring and summer through a process known as “internal loading” or “internal 

recycling” (Kann, 1998). 

A Total Maximum Daily Load (TMDL) for 

Upper Klamath and Agency Lakes was completed in 2002 to address chlorophyll a, pH, and 

dissolved-oxygen levels that did not meet waterquality standards (Oregon Department of Environmental Quality, 2002). The TMDL analysis
5 

Figure 1. Map showing sampling sites for benthic flux and sediment core and grab sample collection in Upper 

Klamath Lake, Oregon.
6 

was partly based on a dynamic mass-balance 

model that simulates the combined Upper Klamath and Agency Lakes as a continuously stirred 

tank reactor (CSTR) and predicts lake total phosphorus, chlorophyll a, and pH as a function of 

phosphorus loading (external and internal) and 

other factors (Walker, 2001). The model tracks 

phosphorus in the water column and sediment. It 

simulates lake-mean phosphorus concentrations 

(is not spatially explicit) and performs massbalance simulations on a biweekly time step. The 

model is programmed into a Microsoft® Excel®

spreadsheet.

The mass-balance model also was used to 

develop boundary conditions for naturalconditions scenarios for the Klamath River 

TMDL models (Oregon Department of Environmental Quality, 2010). The Oregon Department 

of Environmental Quality acknowledged that the 

modeled boundary conditions are one of the 

greatest sources of uncertainty in the development of the Klamath River TMDLs. As the basis 

for load allocations and regulatory actions that 

would affect point and nonpoint discharges 

downstream, operation of irrigation and hydropower systems, and attainment of beneficial uses 

in more than 320 km of the Klamath River to the 

ocean, the effects of the boundary conditions on 

the Klamath River TMDL came under scrutiny 

when reviewers of the Klamath River models 

questioned whether the Upper Klamath Lake 

TMDL model was sufficiently detailed and physically based to support these critical boundary 

conditions (Rounds and Sullivan, 2009; 2010). 

Further research, water-quality monitoring, 

and wetland restoration has occurred in and 

around the lake since the Upper Klamath Lake 

TMDL model was originally developed. The 

TMDL model was calibrated and validated using 

a whole-lake mass balance (Kann and Walker, 

2001) derived from biweekly, year-round, in-lake 

water column profiles and tributary data collected 

by the Klamath Tribes from 1991 through 1998. 

Another decade of monitoring has occurred since 

the model was completed, including a continuation of the Klamath Tribes’ dataset, which is 

documented in a series of reports (for example, 

Kann, 2011a, b). The U.S. Geological Survey 

(USGS) has operated continuous water-quality 

monitors and collected water chemistry data in 

the lake since 2002, and the results are documented in a series of reports (for example 

Eldridge and others, 2012). 

Additionally, research has led to a better 

mechanistic understanding of water-quality processes in the lake, including processes that affect 

the internal loading of phosphorus that is critical 

to the cycling of nutrients and the stimulation of 

algal growth: 

• Laboratory experiments by USGS (Fisher and 

Wood, 2004) demonstrated that the mechanism of desorption of phosphorus at high pH 

may not be the most important mechanism of 

sediment phosphorus recycling. 

• Direct measurements of diffusive fluxes 

through the sediment have been made in Upper Klamath Lake (Kuwabara and others, 

2007; 2009) and the restored Williamson 

River Delta (Kuwabara and others, 2012a). 

• Samples of benthic invertebrate communities

collected in the lake and delta have revealed 

high densities of benthic invertebrates in the 

sediment (Kuwabara and others, 2012a) that 

could contribute to internal loading through 

metabolic cycling. 

• A study of the phosphorus in the top 2 cm of 

sediment before and after the A. flos aquae

bloom in 2005 was equivocal, in that two of 

four sites showed no significant change, one 

site showed a significant increase in residual 

phosphorus, and one site showed a significant 

decrease in iron-associated phosphorus (Simon and others, 2009). 

• In 2006, 26 sediment cores were collected 

and analyzed for nutrients and metals in 10 

vertical layers, which provided a better description of the amount and forms of 

phosphorus in the bottom sediments of the 

lake than was available previously (Simon 

and Ingle, 2011). 
7 

• A hydrodynamic and heat transport model of 

the lake has provided insight into the role of 

the wind-driven circulation in transporting 

water quality constituents (Wood and others, 

2008). 

The phased restoration of deltaic wetlands 

surrounding the mouth of the Williamson River 

was a significant development. The alterations of 

the shoreline associated with this restoration were 

completed in 2009, but maturation of vegetation 

is expected to continue for many years (Elsroad 

and others, 2011), and soil accretion in the wetlands can be expected to occur on the time scale 

of decades (Graham and others, 2005). The newly 

flooded wetlands increased nutrient and DOC 

fluxes to the water column in the months after 

inundation (Wong and others, 2011) as was expected based on preliminary laboratory 

experiments with wetland soils (Aldous and others, 2005; 2007). 

Over time, however, these lands have adjusted and will continue to adjust on time scales that 

are solute specific. Nitrate and dissolved organic 

carbon (DOC) fluxes to the water column dissipated over 3 years, whereas orthophosphate, 

ammonia, iron, and manganese fluxes continue to 

be significant relative to riverine sources (Kuwabara and others, 2012b). Flux of 

orthophosphate and ammonia is declining (Wong 

and Hendrixson, 2011), but over a longer time 

scale than nitrate and DOC. The long-term role of 

the deltaic wetlands in processing water from the 

Williamson River is not certain (over the long 

term the wetland may be a sink for phosphorus 

[Graham and others, 2005]), but it is clear that 

these wetlands create spatial heterogeneity of water and sediment quality within the present-day 

shoreline. 

Because of the additional research and data 

collection since the establishment of the Upper 

Klamath Lake TMDL and the reliance of the 

Klamath River TMDL on results from the Upper 

Klamath Lake TMDL, many agencies that have a 

role in the regulation and management of Upper 

Klamath Lake have agreed that it is appropriate 

to review the status and technical basis of the 

model that informed establishment of the Upper 

Klamath Lake TMDL. As a Designated Management Agency for the TMDL process, the 

Bureau of Reclamation, the cooperator on this 

study, wants to determine whether the tools used 

to build the TMDL are consistent with the most 

up-to-date scientific information and understanding of the major hydrological and biogeochemical 

processes that occur in the lake. Further, periodic 

review of new scientific information related to a 

TMDL is an explicit part of Oregon Department 

of Environmental Quality’s adaptive management 

approach (Oregon Department of Environmental 

Quality, 2002). 

This review begins with a description of the 

most current information about phosphorus fractions in the bottom sediments of the lake, 

followed by a summary of the possible mechanisms of internal loading and how those relate 

specifically to Upper Klamath Lake. The evaluation of the TMDL model is in two parts: 

verification and sensitivity analysis. The verification step assures that (1) the equations provided 

in the documentation are error-free and consistent 

with basic physical, chemical, and biological 

principles such as mass balance and descriptions 

of algal growth, (2) the solution technique can be 

understood and is numerically sound, (3) the 

spreadsheet is consistent with the documentation 

supplied by the developers and does not contain 

coding errors, and (4) the calibration described in 

the documentation can be duplicated and the 

same calibration parameters obtained independently. 

We determined the sensitivity of the calibration parameters to the following factors: (1) the 

mathematical forms of the terms that limit algal 

growth based on the availability of light and 

phosphorus, (2) the assumption that chlorophyll a

achieves equilibrium with ambient light and nutrients during the time step, (3) the assumed 

mechanism for recycling phosphorus from the 

sediments to the water column, and (4) the estimated size of the initial sediment reservoir. The 

first two factors were determined from an exami-
8 

nation of the model code during the initial verification step of the review, the third was based on a 

recognition that recent work has identified several other mechanisms that are likely to be at least 

as important as the single, pH-dependent mechanism that was incorporated into the model, and 

the fourth was the result of recognizing that calculations based on new sediment cores produced 

an estimate of the phosphorus content in the sediments that differed substantially from the value 

used when the TMDL was developed. Finally, the 

report ends with an assessment of where the 

TMDL model is weakest and suggestions for improvements. 

Purpose and Scope

This report presents the results of a technical

review of the existing Upper Klamath Lake

TMDL mass-balance model, including an assessment of the technical competency of the 

model and an evaluation of model assumptions in 

light of information from studies that have been 

done since it was implemented. 

The review had the following objectives: 

1. Check the model code for consistency with 

the model documentation, and identify, implement, and test any easily made 

improvements to the model code. 

2. Revisit the calibration process and determine 

whether and to what extent calibration parameters would change with implementation 

of improvements identified in the first objective or when reasonable and easily 

implemented changes were made in other 

process-based model assumptions. 

3. Quantify the performance of the model 

through a validation process and determine 

the extent to which the performance of the 

model changed with the changes in calibration parameters that were determined in 

meeting the second objective. 

4. Evaluate the sensitivity of long-term model 

simulations to changes made in the process of 

meeting objectives 1 and 2 to assess the consistency of the model results that informed the 

establishment of the Upper Klamath Lake and 

Klamath River TMDLs. 

5. Suggest improvements or enhancements of 

model algorithms based on updated data and 

published research. 

Phosphorus in Upper Klamath Lake 

Bottom Sediments

Phosphorus Content in Sediment Samples Collected Between 1969 and 2006

Simon and Ingle (2011) collected 26 sediment cores in October 2006 throughout the main 

body of Upper Klamath Lake, inside of the shoreline as it existed prior to restoration of the 

Williamson River Delta, for the purpose of quantifying the amount of phosphorus (P) stored in the 

components of the bottom sediments (fig. 1). The 

cores were sampled at 1-cm intervals to a depth 

of 5 cm, and at 5-cm intervals between depths of 

5 and 30 cm. Samples were analyzed for total 

concentration of phosphorus, iron (Fe), titanium 

(Ti), aluminum (Al), and calcium (Ca) using 

sample digestion with nitric and hydrofluoric acids. Total concentrations of carbon (C) and 

nitrogen (N) were determined using flash oxidation of each sample followed by separation of the 

gaseous products using a Carlo Erba EA 1108 

Elemental Analyzer (Simon and Ingle, 2011). A 

sequential extraction was used to quantify (1) the 

loosely sorbed and interstitial water phosphorus 

(LP), (2) phosphorus bound to poorly crystalline 

iron oxides (Fe-P), and (3) phosphorus bound to 

carbonates and acid-soluble minerals (AP), which 

also includes phosphorus bound to aluminum hydroxides. The residual fraction of phosphorus 

(RP), which can include organic phosphorus and 

refractory phosphorus in minerals, was calculated 

by subtracting the phosphorus fractions recovered 

by the sequential extraction from the total P. The 

methods of collection and analysis, as well as the 

data, are reported in Simon and Ingle (2011). 

Results of averaging mass concentrations at 

each depth across all cores showed that phospho-
9 

rus fractions were surface intensified—total

phosphorus at 1 cm depth was about 1.7 times the 

value at 30 cm depth (fig. 2). The most surfaceintensified fraction was the Fe-P fraction, which 

had a surface concentration about 3.4 times the 

value at 30 cm, and the least surface-intensified 

fraction was the AP fraction, which had a concentration about 1.1 times greater at 1 cm than at 30 

cm. Compared to a core collected in 1998 (Eilers 

and others, 2004), the average total phosphorus 

concentration of the 26 cores collected in 2006 

was lower (maximum and minimum concentration over 30 cm depth was 0.5 and 0.3 mg/g DW, 

respectively, for the 2006 cores [fig. 2], and maximum and minimum concentration over 30 cm 

depth was 1.2 and 0.7 mg/g DW, respectively, for 

the 1998 core (Eilers and others, 2004, fig. 4), but 

the enrichment between the surface and 30 cm 

depth is about the same. The 2006 core collected 

closest to the site of the 1998 core (fig. 1) ranged 

in concentration from 0.45 mg/g DW at 0.5 cm to 

0.2 mg/g DW at 30 cm. The difference in absolute 

concentrations attributable to the different methods used (microwave digestion in both cases, 

with HNO3+HF for the 2006 cores and 

HNO3+HF+HCl+H3BO3 for the 1998 core) is 

minor. 

Eilers and others (2004) noted that the surface intensification of titanium and aluminum (40

percent for titanium and more than 200 percent 

for aluminum compared to the concentrations in 

the deepest part of the core) in the 1998 core indicated increased erosional inputs to Upper 

Klamath Lake in the 20th century. The 26 cores 

collected in 2006 show similar surface intensification (fig. 3). Tephra and charcoal increases in 

the upper layers of an additional sediment core 

collected in 1996 also indicated increased erosional inputs to the lake due to river 

channelization and other wetland and watershed 

disturbances over the past 100 years (Bradbury 

and others, 2004). Eilers and others (2004) also 

noted a decrease in the ratio of nitrogen to phosphorus (N:P) in the upper 18 cm of sediment, 

which they interpreted as likely due to an increase in phosphorus loading relative to nitrogen 

loading in the 20th century, although increased 

nitrogen loss processes over the same period also 

could have been partly responsible. The average 

of the 26 cores collected in 2006 shows a similar 

decrease in the N:P ratio from 26.1 at 30 cm to 

14.3 near the surface. Broadly speaking, therefore, the results of Simon and Ingle (2011) are 

consistent with the recent paleolimnological evidence presented by Eilers and others (2004). 

To investigate the possibility of a long-term

trend in the sediment phosphorus content, we 

compared results from the 2006 cores to those 

from two previous studies. The first of these studies analyzed grab samples collected from three 

sites over 16 months during 1969–70 (Wildung 

and Schmidt, 1973; Wildung and others, 1977). 

For each of the sites, the concentration (averaged 

over the top 10 cm) in cores collected nearby in 

2006 was less than the minimum concentration in 

the 12–13 samples collected in 1969–70 by as 

much as 0.13 mg/g DW at the southern end of the 

lake (Pelican Marina) and as little as 0.02 mg/g 

DW at Howard Bay (table 1). The overall variability, however, at each site where samples were 

collected in 1969–70 was large in comparison: 

0.27 mg/g DW at Pelican Marina and 0.66 mg/g 

DW at Howard Bay (table 1). The total inorganic 

phosphorus concentration in the cores collected 

in 2006 was equal to the minimum total inorganic 

phosphorus concentration in nearby samples collected from Howard Bay in 1969–70 and was less 

than the minimum total inorganic phosphorus 

concentration in samples collected from the other 

two sites by as much as 0.1 mg/g DW at the third 

site, Buck Island (table 1). Thus, the concentrations measured in 2006 were lower than those 

measured 36 years earlier at nearby sites, and at 

two of three sites most of the difference between 

the two studies was in the residual (total phosphorus minus sum of extracted inorganic 

phosphorus) fraction. 
10

Figure 2. Graphs showing average concentration of phosphorus fractions as a function of depth in sediment from 26 cores collected from Upper 

Klamath Lake, Oregon, 2006.

Figure 3. Graphs showing average concentration of carbon, nitrogen, and trace elements as a function of depth in sediment from 26 cores collected 

from Upper Klamath Lake, Oregon, 2006.
11

Table 1. Comparison of concentrations of total and inorganic phosphorus in sediment cores collected October 2006 from Upper Klamath Lake, 

Oregon, to reported concentrations prior to 2006.

[Data from 2005 and 2006 were averaged over the top 10 centimeters of sediment; samples were collected in 1969–70 with an Ekman dredge to an approximate depth of 10 cm and homogenized. This table is modified from table 1 in Simon and Ingle (2011). mg/g DW, milligrams per gram of dry weight; 

Reference 1, Wildung and Schmidt, 1973, 13–15 samples collected over 24 months at each site; analysis method for total phosphorus: high temperature ignition, extraction method for inorganic phosphorus: H2SO4; Reference 2, Simon and others (2009), samples collected twice at each site, in April and July; 

analysis method for total phosphorus: microwave digestion, HNO3 and HF, extraction method for inorganic phosphorus: sum of MgCl2 (1.0M), bicarbonatecitrate plus dithionite, and HCl (1.0M); analysis and extraction methods of Simon and Ingle (2011) same as Kuwabara and others (2007). Italicized values were 

read from a graph and are approximate.]

Prior studies Simon and Ingle (2011)

Site location

Year of

sample

collection Reference

Concentration of 

total phosphorus

(mg/g DW)

Concentration of total 

inorganic phosphorus

(mg/g DW)

Closest 

core site

Concentration 

of total 

phosphorus

(mg/g DW)

Concentration of 

total inorganic 

phosphorus

(mg/g DW)

South end of lake 1969-1970 1 0.44 - 0.71 ± 0.02 0.21 - 0.40 ± 0.01 Site 1 0.31 0.18

South of Buck Island 1969-1970 1 0.52 - 0.63 ± 0.03 0.21 - 0.31 ± 0.03 Site 2 0.41 0.11

Howard Bay 1969-1970 1 0.41 - 1.07 ± 0.05 0.13 - 0.61 ± 0.05 Site 3 0.39 0.13

Entrance to Ball Bay 2005 2 0.81 - 0.85 ± 0.11 0.13 - 0.23 ± 0.12 Site 18 0.41 0.16

South of Bare Island 2005 2 0.62 - 0.66 ± 0.15 0.12 - 0.21 ± 0.17 Site 16 0.35 0.15

Howard Bay 2005 2 0.46 - 0.63 ± 0.14 0.14 - 0.16 ± 0.07 Site 3 0.39 0.13

Middle of lake 2005 2 0.37 - 0.46 ± 0.14 0.08 - 0.09 ± 0.05 Site 9 0.19 0.10
12

To investigate further whether the difference 

between the two studies represents a true longterm trend, we compared the concentration in 

cores collected in 2006 to concentrations in cores 

collected from four nearby sites in April and July 

2005 (Simon and others, 2009). The total phosphorus concentration in cores from those sites 

was lower in 2006 than in April or July 2005 by 

as much as 0.4 mg/g DW at the entrance to Ball 

Bay, or as little as 0.07 mg/g DW at Howard Bay. 

The total inorganic phosphorus concentration was 

comparable between the 2 years (slightly higher 

in 2006 than in 2005 at three of the four sites, 

slightly lower at the fourth site [table 1]). Thus, 

the concentrations measured in April and June 

2005 were higher than concentrations measured 

in October 2006 at nearby sites by an amount 

comparable to or greater than the difference between the concentrations measured in 1969–70 

and 2006, and, again, most of the difference was 

in the residual fraction. 

The 2005 cores were collected carefully by 

divers to prevent loss of sediment at the sediment-water interface, whereas the 2006 cores 

were collected with a gravity corer; therefore, the 

collection method could account for part of the 

higher total phosphorus concentration in 2005 

compared to 2006. Samples were collected with 

an Ekman dredge in 1969–70, which likely did 

not preserve the integrity of the sediment/water 

interface as much as diver-collected cores; it is 

also unclear how the depth sampled by the Ekman dredge compares to the top 10-cm of a core. 

Therefore, although the comparison of samples 

collected 36 years apart and others collected 1 

year apart did not provide evidence of a longterm trend in sediment concentration that exceeds 

seasonal and interannual variability, differences 

in sampling methods make quantitative comparisons between these studies tentative. 

Spatial Variability

A limitation in using a CSTR model to simulate a large shallow lake is that all lake quantities 

are averaged, and appropriate averaging requires 

detailed information regarding the spatial variability in those quantities. The whole-lake 

averaged concentrations that provide the calibration and validation data for the CSTR model were 

from samples collected at 11 sites around the lake 

(Kann and Walker, 2001). Little information is 

available, however, on the spatial variability of

processes that affect the mass balance, like cyanobacterial growth and mortality, and phosphorus 

sedimentation to and recycling from the sediments. 

One factor determining the magnitude of 

phosphorus recycling from the sediments is the 

sediment phosphorus concentration. The 26 cores 

analyzed by Simon and Ingle (2011) provided 

previously unavailable information about the spatial variability of the various phosphorus fractions 

in the lake sediment. The data in each core were 

numerically integrated over 10 cm, the assumed 

depth of the active sediment layer in the TMDL 

model. The study used Principal Components 

Analysis (PCA, Haykin, 1999) with varimax rotation (Kaiser, 1958) to reduce the dimensionality 

of the dataset and aid in interpretation. With 

varimax rotation, the principal component axes 

are rotated to maximize loading of a single component onto each variable and minimize loading 

of other components onto this variable, making it

easy to associate each variable with a single 

component on which it has a high loading. The 

final set of principal components was selected on 

the basis of the Kaiser criterion, which retains 

components corresponding to eigenvalues greater 

than one (Kaiser, 1960). Principal component 

analyses were accomplished with the psych package of subroutines (“principal” function) from the 

R software (R Core Team, 2012). The PCA reduced the original 12 sediment constituents to 3 

principal components that explained 82.4 percent 

of the variance in the 10-cm integrated data (fig. 

4, table 2).
13

Figure 4. Diagrams showing varimax rotated principal component scores for all constituents in cores collected from 

Upper Klamath Lake, Oregon, integrated to 10 centimeter depth.

Table 2. Loadings of sediment variables onto varimax 

rotated principal components for all constituents in 

cores collected from Upper Klamath Lake, Oregon, 

integrated to 10 centimeter depth. 

[TP, total phosphorus; RP, residual phosphorus; LP, loosely 

bound phosphorus; AP, phosphorus bound to acid-soluble 

minerals; Fe-P, phosphorus bound to poorly crystalline iron 

oxides; Fe-cryst, poorly crystalline iron oxides]

Component

number

Variance

Explained

(percent)

Variables with 

loading > 0.5

VC1 37.8 AP, Fe, Al, Ca, Ti

VC2 20.2 C, N, LP

VC3 24.4 TP, RP, Fe-P, Fe-cryst

Iron, titanium, calcium, and aluminum are 

indicators of erosional inputs from the basin; 

these metals have similar variability in the bottom 

sediments. This variability is captured by the first 

principal component (VC1), which largely represents the spatial variability in the deposition of 

erosional material. Sites 21 and 23 east of the 

pre-restoration mouth of the Williamson River 

have high scores on this component, reflecting 

erosional inputs from the Williamson River. Under prevailing northwesterly winds, currents flow 

southeastward along the eastern shoreline and 

traverse the lake from east to west between Modoc Rim and Buck Island (Wood and others, 

2008). High scores at sites 17, 12, 8, and 9 indi-
14

cate deposition of erosional inputs along this 

pathway. High scores in Howard Bay and at thesouthern end of the lake (sites 1, 2, and 3; fig. 4) 

indicate inorganic sediment inputs to these areas 

from runoff or anthropogenic loadings. Phosphorus associated with carbonates and aluminum 

hydroxides (AP), also has high loading on VC1. 

Carbon and nitrogen have similar variability 

that is captured by component 2 (VC2), which, 

therefore, largely represents the spatial variability 

in the organic fraction of the sediments. LP also 

is associated with VC2, perhaps because mineralization of organic material results in high 

dissolved phosphorus concentrations in pore water and high concentrations of loosely sorbed 

phosphorus in the sediment matrix. The sites that 

scored highest on this component were along the 

northern shoreline of the lake and in the northern 

bays (fig. 4), indicating that these are areas characterized by high deposition of organic matter. 

Site 7 near Rattlesnake Point also has a notably 

high score on VC2. VC2 was roughly inverse to 

VC1, in that sites in the central part of the lake 

had minimum scores on VC2, but maximum 

scores on VC1. 

The third component, VC3, captures the spatial variability in total phosphorus and its largest 

fractions, Fe-P and RP. The spatial pattern of 

VC3 site scores is more similar to that of VC2 

site scores than to that of VC1 site scores (fig. 4), 

which is largely because the RP fraction, which 

contains organic phosphorus, is the largest fraction of the total phosphorus. The site scores for 

VC3 are highest in the northern and southern 

parts of Upper Klamath Lake and lowest in the 

central part; the highest scores of all are in the 

northern bays. The score in Howard Bay, which 

had a low VC2 score, is notably low in comparison. 

Internal Loading of Phosphorus in

Upper Klamath Lake

Mechanisms

Internal loading of phosphorus (loading directly from bottom sediment to the water column) 

is an important component of phosphorus cycling 

and a primary driver of massive cyanobacterial 

blooms in Upper Klamath Lake (Barbiero and 

Kann, 1994; Kann, 1998). An understanding of 

phosphorus cycling in the lake requires a good 

estimate of the magnitude of the internal loading 

and its temporal and spatial variability. The model estimated whole-lake internal loading from the

residual of the phosphorus mass balance of the 

lake on a 2-week time step for the calibration period of April 1991–September 1998 (Kann and 

Walker, 2001). Internal loading (recycling) was 

estimated to be, on average, 60 percent of the total loading to the lake (Kann and Walker, 2001). 

The time series of the whole-lake internal load 

over 14 days from the original dataset was converted to an areal flux (fig. 5). In each of the eight

growing seasons, the lake-average internal load 

exceeded 10 milligrams per square meter per day 

[(mg/m2

)/d] and in some years exceeded 20 

(mg/m2)/d. Because internal loading varies spatially (fig. 4), localized values of internal loading 

could be many times greater than indicated by the 

whole-lake average.

Several mechanisms contribute to nutrient 

fluxes from the sediments to the water column. 

Resuspension of the sediments can contribute to 

internal loading in shallow lakes (Sondergaard 

and others, 1992; Koski-Vahala and Hartikainen, 

2001; Niemisto and others, 2011). No direct connection between episodes of resuspension and 

internal loading for Upper Klamath Lake has 

been established, although resuspension can result when wind speed exceeds 4.5 m/s, a 

relatively common occurrence (Laenen and 

LeTourneau, 1996). Further, resuspended sediments from Upper Klamath Lake did not release 

more phosphorus in the long-term than sediments 

at a stable interface in laboratory experiments, 
15

Figure 5. Graph showing calculated whole-lake average of internal loading to Upper Klamath Lake, Oregon, April 

1991–September 1998. Only positive values, representing a flux of phosphorus from the sediments to the water 

column, are shown. Calculation of the average flux is based on data from Kann and Walker (2001).

and in the short-term had the capacity to sorb 

phosphorus (Fisher and Wood, 2004). Likely, resuspension in Upper Klamath Lake adds 

primarily sorbed (not bioavailable) phosphorus to 

the water column that does not directly promote 

the growth of cyanobacterial blooms or chlorophyll a production. Therefore, fluxes of dissolved 

phosphorus from the bottom of Upper Klamath 

Lake also must occur through pore water pathways by means of physical transport, either by 

molecular diffusion or advection (by groundwater 

from below or by surface water flowing horizontally and entering interstitial spaces). 

Phosphorus in the solid phase, measurements 

of which were described by Simon and Ingle 

(2011), can be made available for transport in the 

aqueous phase through various means, and the 

phosphorus fraction that is the primary source of 

dissolved phosphorus varies accordingly. The 

origin of pore water phosphorus can be chemical 

desorption (primarily from the LP fraction) or, in 

deeper, reducing sediments, dissolution of iron 

oxides (source is Fe-P fraction). Desorption can 

also occur from iron oxides in oxic conditions 

when pH exceeds 9 and hydroxide ions begin to 

displace phosphate ions from sorption sites 

(source is Fe-P fraction; Jacoby and others, 1982; 
16

Sondergaard, 1988). In the process of mineralization, bacteria convert organic sediment 

phosphorus to biologically available dissolved 

inorganic phosphorus (source is RP fraction). 

Benthic macroinvertebrates can contribute to 

benthic fluxes by enhancing diffusion through 

bioturbation, which increases the surface area in 

contact with pore water and creates macro-scale 

pathways for pore water transport (primary 

source is the LP fraction), and by excretion (primary source is the RP fraction), either to pore 

water or directly to the overlying water column.

Fluxes of solid phase phosphorus can occur 

through biotranslocation of akinetes (dormant 

cyanobacteria cells) or invertebrates that undergo 

episodic recruitment (source is the RP fraction). 

Many or all of these processes may overlap in 

time and space in shallow lakes like Upper Klamath Lake; the key to constructing a successful 

model is to have a good quantitative understanding of the most important processes. 

Direct measurements, or estimates by other 

means, of some of these processes have been 

made in Upper Klamath Lake: 

• Kuwabara and others (2007; 2009) measured 

diffusive fluxes across the sediment-water interface before, during, and after the rapid 

expansion of the cyanobacterial bloom in 

2006, using a pore water profiler that obtained pore water concentration profiles at the 

centimeter scale. Benthic fluxes of phosphorus measured in this way were consistently 

positive (from the sediments to the overlying 

water column) and ranged from 0.4 to 6.1 

(mg/m2

)/d. 

• Kuwabara and others (2009) also used heatflux modeling to estimate fluxes due to 

groundwater advection of 1.1 (mg/m2

)/d. 

• Fisher and Wood (2004) mimicked the process of desorption from iron oxides at high 

pH in a laboratory using sediment from Upper 

Klamath Lake. In the study of Fisher and 

Wood (2004), the estimated flux was 1.9 

(mg/m2)/d at pH 9.2 and 1.3 (mg/m2)/d at pH 

10.3. 

• Additional estimated fluxes of Fe-P can be 

inferred from core studies. Using data from 

cores collected by Wildung and Schmidt 

(1973) to calculate the loss of phosphorus 

from the Fe-P fraction as an areal flux during 

27 days between May 7 and June 3, 1969, 

gives an estimate of 1.8–2.5 (mg/m2

)/d (using

sediment bulk densities at nearby sites from 

Simon and Ingle, 2011). 

• Simon and others (2009) found a significant 

decrease in total phosphorus at one of their 

core sites during April–July 2005 that was associated primarily with iron oxyhydroxides. 

Averaging Fe-P over the top 10 cm of sediment and duplicate cores at that site (data 

provided in Kuwabara and others, 2007) gives 

an estimate of 0.73 (mg/m2

)/d for the loss of 

Fe-P during the 69 days between core collections. 

• Assuming that the acid-soluble phosphorus 

fraction is composed primarily of phosphorus 

bound to aluminum hydroxides (Simon and 

Ingle, 2011), which is also subject to potential 

desorption at high pH, then using data from 

the same site gives an estimate of the loss of 

Al-P of 0.45 (mg/m2)/d during the 69 days between core collections. 

• Barbiero and Kann (1994) measured the flux 

of phosphorus to the water column during cyanobacterial (A. flos aquae) recruitment at 

3.56 (mg/m2)/d, which they concluded was 

about 27 percent of the concurrent increase in 

total water column phosphorus. 

Each of these estimates indicates that the 

process involved could be episodically significant, but not likely to be the dominant process of 

internal loading, particularly during rapid, early 

spring A. flos-aquae growth, when calculations of 

internal loading can be between 10 and 20 

(mg/m2

)/d on a whole-lake basis (fig. 5). 

At the time of the development of the Upper 

Klamath Lake TMDL, the best available description of the composition of the phosphorus in the 

bottom sediments was by Wildung and Schmidt 

(1973) and Wildung and others (1977). Sediment 
17

samples used in these studies were collected approximately monthly between May 1969 and 

April 1970 from three sites at the southern end of 

Upper Klamath Lake—one near the outlet, one 

on the south side of Buck Island, and one in the 

southwest corner of Howard Bay. Wildung and 

Schmidt (1973) found that changes in total sediment phosphorus correlated with changes in total 

inorganic P. A subsequent serial extraction of 

various phosphorus fractions in a subset of the 

samples showed a decrease in phosphorus associated with hydrated iron oxides that coincided 

with a period of rapid growth of phytoplankton in 

spring 1969 at all three sites. Wildung and 

Schmidt (1973) concluded that it was primarily 

this fraction that provided bioavailable phosphorus to the water column. This conclusion was 

consistent with an identified mechanism for internal loading—desorption of Fe-bound 

phosphorus under oxic conditions at high pH 

(Jacoby and others, 1982)—and with observations that the recycling of phosphorus to the 

water column increased with pH (Kann, 1998), 

and this conclusion ultimately informed the way 

that internal loading was conceptualized and incorporated into the TMDL model. 

The conclusion, however, was largely based 

on data from one site—the Howard Bay site—

where sediment phosphorus concentrations were 

highest and most variable. For example, in sediment samples collected in April 1970, total 

phosphorus concentrations in Howard Bay, south 

of Buck Island, and near the lake outlet were 

1.03±0.05, 0.63±0.03, and 0.66±0.02 mg/g DW, 

respectively, and total inorganic phosphorus concentrations at the same three locations were 

0.57±0.05, 0.30±0.03, and 0.35±0.01 mg/g DW, 

respectively (table 3 in Wildung and Schmidt, 

1973). The range in concentration during the 16 

months of sample collection was 0.41–1.07 mg/g 

DW at Howard Bay, 0.52–0.63 mg/g DW near 

Buck Island, and 0.44–0.71 mg/g DW at the 

south end of the lake (table 1). The decrease between the two sample dates was greatest at the 

Howard Bay site, and inorganic phosphorus (as 

determined by H2SO4 extraction) accounted for 

80 percent of the reduction in total phosphorus 

(determined by high temperature ignition). 

The serial extractions that followed to distinguish which fraction of sediment phosphorus 

accounted for the loss of sediment phosphorus 

were incomplete, based on a comparison between 

the sum of the extracted fractions and the totals, 

as determined by high temperature ignition (used 

to determine total P) and H2SO4 extraction (used 

to determine total inorganic P). The residual organic fraction at the end of the serial extraction 

procedure ranged between 0.7 and 29 percent of 

the sum of the inorganic fractions. When the organic fraction was instead determined by 

difference between total phosphorus and total inorganic phosphorus, it was generally comparable 

to and often greater than the total inorganic phosphorus (Wildung and Schmidt, 1973). During 

most of the 16-month sampling period, and considering all three sites where samples were 

collected, changes in the total sediment phosphorus were correlated with changes in the both 

organic and inorganic fractions (based on visual 

inspection of figure 3 in Wildung and Schmidt, 

1973). Therefore, it is consistent with these early 

studies to conclude that fractions other than ironbound phosphorus contributed substantially to 

changes in sediment P, and that mechanisms of 

internal loading that utilize other phosphorus 

fractions as a source are as important as the single 

mechanism (desorption of iron-bound phosphorus 

under oxic conditions at high pH) that was incorporated into the TMDL model. The importance of 

the organic phosphorus fraction, in particular, as 

a source for internal loading might have been underestimated. 

Possible Role of Benthic Invertebrates in 

Internal Loading in Upper Klamath Lake

Benthic invertebrates in the bottom sediments of Upper Klamath Lake occur at densities 

that could influence diffusive fluxes through bioturbation and metabolism, both of which have the 

potential to be equally as important as the other 

physical and chemical processes discussed in the 

previous section (Andersson and others, 1988). 
18

During 2009 and 2010, we conducted a 

quantitative study of the distribution of benthic 

invertebrates in the northern part of the lake. Five 

samples were collected at five randomly selected 

sites on three dates (early June, mid-July, and 

early September) using a 0.023 m2 Ekman 

dredge. Samples were sieved using a 500 μm 

mesh and sorted at 10X magnification. Mean 

density of all invertebrates sorted was 20,625 individuals/m2

 ±11,671 (±standard deviation, 

N=150). The number of individuals per sample

ranged from 1,822 to 62,178 individuals/m2. 

Ninety-four (94) percent of all organisms comprised Oligochaeta (worms, 65 percent),

Hirudinea (leeches, 17 percent), and Chironomidae (midges, 12 percent). 

Benthic invertebrates can be sorted into functional groups based on a mechanistic assessment 

of their sediment mixing modes (Michaud and 

others, 2006). Differences in these modes likely 

influence the sediment chemistry and the solute 

distribution. For example, “upward conveyers” 

orient vertically with their heads down, and in the 

process of feeding remove sediment at depth and 

expel excrement at the sediment-water interface. 

“Gallery diffusers” dig extensive burrow systems 

and transport materials through the burrows by 

body movement. The bioturbation effects of the 

two groups are quite different. Gallery diffusers 

are effective at oxygenating the sediment to the 

maximum depth of their burrows, stimulating microbial processes. Michaud and others (2005) 

found that gallery diffusers increased oxygen uptake by a factor of two over uninhabited controls, 

and Michaud and others (2006) found that gallery 

diffusers were more effective at enhancing nutrient fluxes than another functional group—

“biodiffusers,” which move particles randomly 

over short distances, thus enhancing sediment 

diffusion (bivalves belong to this group, but are 

rare in the open waters of Upper Klamath Lake). 

When a gallery diffuser and an upward conveyer 

were compared directly, the gallery diffuser was 

more effective than the upward conveyer in increasing the depth of oxygen penetration into the 

sediment (Zhang and others, 2010), but when an 

upward conveyer and a “downward conveyer” 

(surface feeder) were compared directly, the upward conveyer more effectively enhanced solute 

fluxes (Quintana and others, 2007). 

Thus, the mode of feeding and sediment mixing must be considered on a taxon-specific basis 

to accurately assess the effect of each taxon on 

nutrient fluxes. For example, there are several 

species of tubificid worms (Oligochaeta) in Upper Klamath Lake, some of which create 

extensive galleries and move up and down 

through them, making them functionally gallery 

diffusers, whereas others simply feed head down 

without creating extensive galleries, making them 

functionally upward conveyers. The largest 

midge in Upper Klamath Lake, Chironomus plumosus, creates U-shaped tubes through which it 

brings water from the sediment surface at one end 

and expels it at the other end (functionally a mix 

of a downward conveyer and a gallery diffuser). 

At least one relatively common midge is predacious and moves about the sediment-water 

interface, acting functionally like a biodiffuser. 

With respect to the flux of phosphorus in 

particular, the chemical composition of the sediments is also important. In iron-rich sediments, 

where much of the phosphorus is in the ironbound fraction, the increased oxidation of the 

sediment to the burrow depth can decrease the 

flux of phosphorus instead of increasing it by enhancing the precipitation of iron oxyhydroxides 

to which phosphate is adsorbed or with which 

phosphate is coprecipitated. Zhang and others 

(2010) reported that in Lake Taihu, China, both 

worms and chironomids inhibited phosphorus 

release from bottom sediments, but the concentration of iron in both sediments (29.6 mg/g DW, 

top 4 cm) and in pore water (maximum about 14 

mg/L) in Zhang and others (2010) was greater 

than concentrations measured in Upper Klamath 

Lake. In contrast, sediment iron concentrations in 

Upper Klamath Lake were 16.1 mg/g DW in the 

top 4 cm, averaged over 26 cores (Simon and Ingle, 2011), and maximum iron pore water

concentrations were about 3 mg/L (Kuwabara 

and others, 2009). 
19

The influence of benthic invertebrates on nutrient fluxes also must be assessed in the context 

of the physical nature of the substrate and the hydrologic context of the environment (MermillodBlondin, 2011). In a comparison of tubificid

worms and chironomids, only the worms were 

effective at reworking fine-grained sediment deposited by stormwater and river channels, and the 

effectiveness of the tubificids was dependent on 

high organic content in the sediment (Nogaro and 

others, 2009). Mermillod-Blondin and Rosenberg 

(2006) found that bioturbation could strongly affect fluxes at the sediment-water interface in a 

diffusion-dominated (lentic) system, but had a 

smaller influence in an advection-dominated (lotic) system. Much of Upper Klamath Lake has

fine-grained sediments with moderately high organic content, and diffusion-dominated 

hydrologic processes; therefore, bioturbation 

could have the maximum possible effect under 

these conditions. There are areas of the lake, 

however, where wind-driven currents can reach 

tens of centimeters per second, and in those areas 

the effect of bioturbation on fluxes could be 

moderated.

The preceding discussion and reported 

measurements do not distinguish between the 

physical enhancement effects of bioturbation on 

diffusion and microbial degradation and the fluxes due to metabolic processing by the 

macroinvertebrates. Several studies have directly 

measured the phosphorus excretion of benthic 

invertebrates (table 3). Studies of two numerous 

benthic invertebrates in Upper Klamath Lake—

chironomids and tubificids—produced measurements of excretion that varied over a wide range, 

from 0.07 nanomoles phosphorus per milligram 

DW per hour (nmol P/mg DW/h) to 3.46 nmol 

P/mg DW/h. Significant differences among studies in the rate of phosphorus excretion are likely a 

result of differences in the species (or, more 

commonly, groups of species within a higher taxon) evaluated, temperature, study location, and 

methods used for testing and analyses. There was 

a positive dependence on temperature, but with 

regard to differences between the two types of 

invertebrates, the results were equivocal. Gardner 

and others (1981), Fukuhara and Yasuda (1985) 

and Fukuhara and Sakamoto (1987) measured 

higher excretion rates in chironomids, whereas 

Devine and Vanni (2002) measured higher excretion rates in tubificids. 

We measured excretion rates of phosphorus 

in Upper Klamath Lake resident fauna twice in 

the midge Chironomus plumosus and once in the 

tubificid worm Ilyodrilus frantzi. Animals were 

collected from Upper Klamath Lake on May 28, 

2012, and experiments were completed within 2 

days of collection using methods similar to Gardner and others (1981). Animals were kept with 

aeration in native sediment and at ambient lake 

temperature (approximately 12°C) prior to each 

experiment. All animals were individually rinsed 

twice at the beginning of each experiment to remove adhered debris. The first rinse was in 

filtered lake water and the second was in low

phosphorus artificial lake water (ALW). For each 

midge experiment, 25 fourth instar midges were 

placed separately in 5 mL of ALW in a 15 mL 

trace-contaminant-clean Teflon test vial and allowed to remain undisturbed at lake ambient 

temperature for 8 h. In contrast to a single midge, 

for the tubificid experiment, five organisms were 

used per test vial because of their much lower 

mass. 

After 8 h, test water was filtered using 0.22 

μm polycarbonate membranes and phosphorus 

was measured using U.S. Environmental Protection Agency (EPA) method 365.2 by batchautomated spectrophotometry (Aquakem® 250, 

Thermo Scientific®). Mean release rates were 

0.035 μg P/midge/h (±0.027, N=25), 0.030 μg 

P/midge/h (±0.054, N=25), and 0.001 μg 

P/tubificid/h (±0.0008, N=25) in the three experiments. Procedural blanks, which included a 1-

sec dipping of a previously rinsed midge in the 

blank ALW to account for remaining surface contamination, contained 0.001 mg P/L (±0.0012, 

N=5), which represented 2.4 percent of the mean 

phosphorus released (0.053 mg P/L) by the 8-h 

tested midges. Measurements of excretion rates
20

Table 3. Comparison of mass-specific excretion of phosphorus by benthic invertebrates from Upper Klamath 

Lake, Oregon. 

[nmole/mg DW/h = nanomole per milligram dry weight per hour, '-' = no data]

Taxon

Mass-specific excretion rate

Temperature

(oC)

Incubation 

time

(hours)

Nitrogen

(nmole/mg 

DW/h)

Phosphorus

(nmole/mg 

DW/h) Source

Chironomidae (Chironomus cf 

plumosus) 12 8 - 0.24 This study

Chironomidae (Chironomus cf 

plumosus) 12 8 - 0.21 This study

Tubificidae (Ilyodrilus frantzi) 12 8 - 0.05 This study

Chironomidae (July experiment) 5 24 - 0.57 Gardner and others (1981)

Chironomidae (September

experiment) 5 24 - 0.27 Gardner and others (1981)

Tubificidae 5 24 - 0.16 Gardner and others (1981)

Chaoborus sp. 9 4 4.34 0.99 Devine and Vanni (2002)

Chironomidae 9 4 4.48 0.28 Devine and Vanni (2002)

Tubificidae 9 4 2.34 0.47 Devine and Vanni (2002)

Chaoborus sp. 14 4 15.86 1.87 Devine and Vanni (2002)

Chironomidae 14 4 13.08 0.56 Devine and Vanni (2002)

Tubificidae 14 4 15.24 2.43 Devine and Vanni (2002)

Chaoborus sp. 15 24 - 0.78 Fukuhara and Yasuda 

(1985)

Chironomidae ( Chironomus) 15 4 - 0.19 Fukuhara and Yasuda 

(1985)

Chironomidae (Tokunagayusurika) 15 4 - 0.10 Fukuhara and Yasuda 

(1985)

Tubificidae (Limnodrilus) 15 7 - 0.07 Fukuhara and Yasuda 

(1985)

Chironomidae 20 4 4.75 0.35 Fukuhara & Sakamoto 

(1987)

Tubificidae 20 4 5.10 0.10 Fukuhara & Sakamoto 

(1987)

Tubificidae 20 0.5 - 3.46 Gardner and others (1981)

Chironomidae 20 24 - 0.81 Gardner and others (1981)

Tubificidae 20 24 - 0.19 Gardner and others (1981)

Chaoborus sp. 22 1 47.63 4.70 Devine and Vanni (2002)

Tubificidae 22 54.75 3.37 Devine and Vanni (2002)

Chironomidae 25 24.40 1.46 Devine and Vanni (2002)
21

in molar units from the two midge experiments 

(0.241 and 0.207 nmol P/mg DW/h) and the single tubificid experiment (0.054 nmol P/mg 

DW/h) are similar to those tabulated by Devine 

and Vanni (2002). 

To evaluate the contribution to internal loading of nutrient excretion by the benthos, the 

mass- and species-specific excretion rates were

converted to an areal rate. The areal rate at any 

point in time and geographic location depends on 

the density of individuals and their body size, 

species present, life stage, and temperature. Nevertheless, to get an approximation of the 

importance of excretion by benthic invertebrates 

to internal loading in the lake, we used the overall 

estimated mean densities of tubificids and chironomids derived from sampling during the 2-

year period and mean values for the rate of phosphorus excretion per taxon. An estimate of 

phosphorus release through excretion based on 

mean densities by these two taxa is 3.8 

(mg/m2

)/d. This value represents an estimate of 

the metabolic addition of phosphorus to the water 

column by tubificid worms and chironomid 

midges, but does not address the influence of bioturbation nor the contribution of phosphorus by 

other abundant taxa (leeches, for example). Although there is likely an influence of bioturbation 

on phosphorus cycling given the overall density 

and behaviors of the resident benthic invertebrates, whether this influence enhances the 

release of phosphorus is a function of many other 

biological, chemical, and mechanical processes.

For the purposes of a numerical model, it is 

necessary to adequately estimate benthic invertebrate excretion rates of specific nutrients. These 

rates are influenced by the spatial and temporal 

variation in densities of individual species or 

higher taxa throughout Upper Klamath Lake. 

Seasonal variability in taxon presence and their 

growth rates are likely dependent on temperature 

(Vannote and Sweeney, 1980). The most important influence on interannual variability in 

spring nutrient release could be year-to-year differences in the spring progression of water 

temperature, because these temperatures likely 

influence the timing of the first midge emergence

(Vannote and Sweeney, 1980). A cold spring 

would likely delay the emergence of the first 

midge cohorts relative to a warm spring. Measurements of water quality in the lake support the 

correlation of the onset of internal loading with 

spring temperatures and suggest a biological control on some of the important processes (Wood 

and others, 1996; Morace, 2007). In addition to 

temperature, year-to-year differences in food resources, predation, competition, and 

overwintering success also have a role in determining the size and distribution of Upper 

Klamath Lake benthic invertebrate populations. 

The benthos affects not only the magnitude 

of fluxes of phosphorus and nitrogen from the 

sediments to the water column, but also the nutrient ratio in those fluxes. Alves and others (2010) 

found that mass-specific nitrogen and phosphorus 

body content varied among three species of benthic invertebrates, and that the mass-specific N:P 

excretion ratio varied among those species as 

well, such that it was positively related to mass 

specific body content N:P ratio. Thus, a more 

complete understanding of relative contributions 

of the three dominant benthic invertebrate taxa in 

Upper Klamath Lake would include a comparison 

of their body stoichiometries and how those vary 

with body size, which could be incorporated into 

growth-based models of excretion. 

The possible influence of benthic invertebrates on the N:P ratio of fluxes from the 

sediments in the larger context of trophic level 

cycling must be considered. Sereda and Hudson 

(2011) noted that the N:P ratio in excretions of 

aquatic metazoans was generally below the Redfield ratio of 16:1 (a measure of the average N:P 

ration in aquatic phytoplankton [Redfield, 1934]). 

Hansen and others (1998) demonstrated that the 

effect of chironomids on nutrient sediment fluxes 

was not only an increase in fluxes, but also a 

lowering of the N:P ratio of those fluxes. In laboratory experiments, the N:P ratio of nutrient 

fluxes from unamended sediment cores was 14.3 

without the addition of chironomids, and 6.8 with 

the addition of chironomids. The addition of fresh 
22

algal material stimulated mineralization processes 

and increased nutrient fluxes from sediment cores 

with and without chironomids, but in cores with 

chironomids, the N:P ratio of those fluxes was 

lower—8.1 in chironomid-inhabited versus 9.9 in 

cores without chironomids (Hansen and others, 

1998). This implies that the N:P ratio of sediment 

fluxes into the water column of Upper Klamath 

Lake is possibly lower than it would be in the absence of the large benthic invertebrate densities. 

Thus, the role of the benthic invertebrates could 

include not only the efficient recycling of phosphorus from the sediments, but also the lowering 

of the N:P ratio to levels that favor the dominance 

of the nitrogen fixing cyanobacterium A. flos aquae over phytoplankton that are not capable of 

fixing nitrogen. 

A final consideration regarding the role of 

benthic invertebrates is that their densities vary 

spatially as well as temporally. Because the concentration of phosphorus in various sediment 

fractions varies spatially as well (fig. 4), the overlay of the density of individuals on the spatial 

distribution of phosphorus is important to determining the magnitude of their contribution to 

internal phosphorus loading in the lake. Such spatial variability cannot be resolved with a CSTR 

model. Incorporation of invertebrates into a 

CSTR TMDL model requires an intermediate 

step designed to properly scale up from individual locations, and is not likely to be as simple as 

taking a straightforward average of the quantities 

involved. 

Verification of the Upper Klamath 

Lake Phosphorus TMDL Model

The Upper Klamath Lake TMDL model 

simulates the lake as a continuously stirred tank 

reactor (CSTR) with a 2-week time step. It is 

programmed into an Excel® spreadsheet, and the 

datasets required to run the model from April 

1991 to September 1998 (hereafter referred to as 

the 7.5-year calibration dataset [but includes eight 

summers]) are supplied with the spreadsheet. The 

model requires measured hydrologic and climatologic inputs (external phosphorus loads, lake outflow, lake volume, solar radiation, air 

temperature, and water temperature), the sources 

of which are documented in Walker (2001). The 

model also requires several parameters that are 

determined through a calibration process or compiled from the literature. Those parameters also 

are documented in Walker (2001) and listed in 

table 4. 

The model has four component submodels: a 

water-column total phosphorus model, a sediment

total phosphorus model, a chlorophyll a model, 

and a pH model. These four components are interconnected and interdependent at each time step

(fig. 6) and are determined through an iterative 

process by intentionally allowing circular references and enabling the iterative calculation 

option in the spreadsheet. Required initial conditions include the initial water-column and 

sediment total phosphorus concentrations, determined from collected data. At each time step, the 

model evaluates total phosphorus in the water 

column, settled total phosphorus (algal and nonalgal components separately), buried total 

phosphorus, recycled total phosphorus from sediments, chlorophyll a concentration, and pH. The 

model can be run in two modes: the calibration 

model runs the model one time only through 7.5 

years of measured hydrologic and climatologic 

conditions; the TMDL model chains together the 

first 7 years of these same measured conditions 

(mid-April 1991 to mid-April 1998), end to end. 

When running in “TMDL mode,” the model continues to run through the same 7 years of 

measured conditions until a steady state (equilibrium) is reached. A short description of each 

submodel follows.
23

Table 4. Values and sources of parameters used in the Upper Klamath Lake, Oregon, phosphorus TMDL 

model—continued

[For definitions of variables, see table 5. P, phosphorus, chl a, chlorophyll a, -, no units]

Symbol Description

Value or value 

range Units Source

Values changed in one or more modified versions of the TMDL model

Kna Gross removal rate for non-algal 

P 

0.012–0.015 per day Calibrated using observed 

mass balance and total P 

data

khalf Half-saturation constant for P

limitation of algal growth

0.06 parts per million Calibrated within range provided by Bowie and others, 

1985

ualg Net settling velocity for algal P 0.001–0.116 meters per day Calibrated using observed 

total P and chl a data

pH∗ pH at which the P recycle rate is 

one half the maximum

8.1–9.1 - Calibrated using observed 

total P, chl a and pH data

KR Maximum P recycle rate, pHdependent

0.26–1.48 per year Calibrated using observed 

total P, chl a and pH data

KT Maximum P recycle rate, temperature-dependent

0.99 per year Calibrated using observed 

total P, chl a and pH data

T(min,rec) The minimum temperature at 

which temperature-dependent 

recycling begins

12.8 degrees Celsius Calibrated using observed 

total P, chl a and pH data

T(max,rec) The maximum temperature allowing temperature-dependent 

recycling

23.7 degrees Celsius Maximum water temperature 

in observations provided 

with model

P(sed_new) Initial sediment P concentration, 

updated based on cores collected 

in 2006

335 milligrams per 

kilogram

Simon and Ingle, 2011

Values unchanged from original TMDL model

σobs Spatial standard deviation of pH 0.25 - Calculated and reported in 

Walker, 2001

α Non-algal light extinction coefficient

1.32 per meter Calculated and reported in 

Walker, 2001

β Algal light extinction coefficient 0.0097 per (meter × parts 

per billion)

Calculated and reported in 

Walker, 2001

IK Onset of light saturation 214 microeinsteins per 

square meter per 

second

Talling, 1957

f Light reflectance fraction 0.1 - Wetzel, 1975

Tmax Temperature at which algal 

growth rate is maximum

20 degrees Celsius Bowie and others, 1985
24

Table 4. Values and sources of parameters used in the Upper Klamath Lake, Oregon, phosphorus TMDL 

model—continued

[For definitions of variables, see table 5. P, phosphorus, chl a, chlorophyll a, -, no units]

Symbol Description

Value or value 

range Units Source

Tmin Temperature at zero growth rate 14 degrees Celsius Bowie and others, 1985

KG Maximum growth rate at reference temperature, TX

1.2 per day Bowie and others, 1985

YP Algal cell P quota 0.6 milligrams P per 

milligrams chl a

Bowie and others, 1985

Kresp

Maximum respiration rate at 

reference temperature, TX

0.06 per day Bowie and others, 1985

ubur Burial velocity of P in the sediment 

1.4 millimeters per 

year

Eilers and others, 2001

zsed Active sediment thickness 10 cm Eilers and others, 2001

ρbulk Sediment bulk density 0.12 grams per cubic 

centimeter

Eilers and others, 2001

Psed Initial sediment P concentration, 

as used in original model, based 

on 1996 core

1,000 milligrams per 

kilogram

Eilers and others, 2001
25

Figure 6. Diagram showing original total maximum daily load (TMDL) model dependence loops. Arrows point from 

the dependent variable to the independent variable and each filled box is a cell in the Excel model. Bold lines indicate cells using circular references.

To better understand and work with the algorithms of the TMDL model, an R-script version 

of the model was developed. The R-script enabled expedited duplication of model calibration, 

substitution of various recycling methods for sensitivity analysis, and user-friendly modeling of 

full future time periods. The R-script version of 

the model duplicated the results of the spreadsheet model, as indicated by statistics describing 

the comparison between the simulated data from 

both models—perfect correlation (R2

=1.0), zero 

bias, zero root mean square error (RMSE), and a 

Nash-Sutcliffe statistic (NS) of 1.0 for all variables. Because the R-script model duplicated the 

results of the spreadsheet model, all of the modifications discussed in the rest of the document, as 

well as the duplication of model calibration and 

validation as presented in Walker (2001) were 

produced with the R-script model.
26

Water-Column Phosphorus Submodel

The basis of the water-column phosphorus 

submodel is the following mass balance equation, 

excluding unit conversion factors which are included in the model: 

dM

dt = QinPin + REC × S − (SED + Qout

V ) × M, (1)

where QinPin and REC × S are the external and 

internal load over the time step, respectively; 

SED × M is the sedimentation rate of algal and 

nonalgal phosphorus, and (

Qout

V ) × M is the export of phosphorus at the lake outlet. The REC

and SED terms contain calibration parameters 

that are discussed below:

REC = KRProb(pH > pH∗), (2) 

where KR and pH* are calibration parameters (table 4) and Prob(pH > pH*

) is the predicted 

probability of pH in the lake exceeding pH* as 

determined from a normal distribution with mean 

equal to pH and standard deviation σobs. The 

probabilistic relation was calibrated and verified 

with site-specific pH observations using (Walker, 

2001): Prob(pH > pH∗) = Freq(pHsite > pH*). 

The sedimentation term includes loss terms for 

algal and nonalgal phosphorus, assuming that the 

ratios of algal and nonalgal phosphorus to total

phosphorus remain constant over a time step:

SED = Knarna + ualg

z ralg, (3)

where Kna and ualg are calibration parameters. The 

solution to equation 1 over a time step takes the 

exponential form: 

M(t) = Mss + (M0 − Mss)e−�SED+Qout

V �t, (4) 

where

MSS = QinPin+REC×S

SED+Qout

V

. (5) 

All variables are defined in table 5. The solution 

proceeds as follows: At each time step, the average mass (Mavg) over the time step is calculated 

using the mass at the beginning (M0) and the 

mass at the end (M14) of the 14-day time step as 

defined by equation 4. This value of Mavg is then 

used to calculate the sedimentation and recycle 

terms in the sediment and chlorophyll a models 

(below), then the pH is calculated using the pH 

model (below). The new value of pH is then substituted back into equation (2) and the REC term 

then recalculated and substituted into equation 5. 

This process repeats until the solution converges. 

Once the solution converges, M14 becomes M0 for 

the next time step, the initial time (t0) is reset to 

zero, and the iterative procedure begins again. 

Sediment-Phosphorus Submodel

The basis of the sediment-phosphorus submodel is the mass balance equation for 

phosphorus contained in the sediment reservoir: 

dS

dt = SED × M − REC × S − BUR × S, (6)

where

BUR = uburzsed. (7) 

The solution to equation 6 over a time step takes 

the exponential form: 

S(t) = Sss + (S0 − Sss)e−(REC+BUR)t, (8)

where

SSS = SED×M

REC+BUR . (9)

Equation 8 is used to determine the average sediment mass (Savg) over the time step. The new 

value of Savg is used to calculate a new value for 

the internal loading term in equation 1, which 

leads to a new value of Mavg. This process is repeated iteratively with each new value of Mavg

and Savg until the solution converges. 
27

Table 5. Definitions of symbols used in this report—continued

[For definitions and values of parameters, see table 4. All mass, concentration, and pH variables (except pHsite) in the water 

column are lakewide averaged quantities. -, no units; P, phosphorus]

Symbol Units Description

A square meters Average lake surface area over the time step

B parts per billion Chlorophyll a concentration 

Bobs parts per billion Observed chlorophyll a concentration 

BEQ parts per billion Chlorophyll a concentration evaluated by setting dB/dt to 0

BRK parts per billion Chlorophyll a concentration evaluated using a numerical Runge-Kutta 

approach

Ec per meter Light extinction coefficient

FT - Temperature adjustment factor

FP - Phosphorus limitation factor

FPFS - Forsberg Shapiro phosphorus limitation factor

FP00 - Phosphorus limitation factor calculated as the percent of non-algal P

FPMM - Phosphorus limitation factor calculated by Michaelis-Menten kinetics

FL - Light limitation factor

FLAVG - Light limitation calculated using the average light intensity

FLDI - Average light limitation calculated by depth-integrating the light limitation equation

fp - Photoperiod (fraction of day)

I microeinsteins per 

square meter per 

second

Light intensity at a particular depth

I0 microeinsteins per 

square meter per 

second

Light intensity at lake surface

Imean microeinsteins per 

square meter per 

second

Average light intensity over the lake depth

J - Julian day

KR per day Recycle rate of P from the sediment 

KG per day Maximum algal growth rate 

ks per day Removal rate of algal P 

M kilograms Mass of total P in the water column 

Mload kilograms Mass of external total P loading to the lake 

Mload,obs kilograms Mass of external total P loading to the lake, biweekly observations 

Mrec kilograms Mass of P recycled from the sediment 
28

Table 5. Definitions of symbols used in this report—continued

[For definitions and values of parameters, see table 4. All mass, concentration, and pH variables (except pHsite) in the water 

column are lakewide averaged quantities. -, no units; P, phosphorus]

Symbol Units Description

Mrec,obs kilograms Mass of P recycled from the sediment, biweekly observations

Msed kilograms Mass of P settling from the water column 

Mout kilograms Mass of total P flowing out of the lake 

Mout, obs kilograms Mass of total P flowing out of the lake, biweekly observations 

Mss kilograms Steady state mass of water column total P calculated at each time step 

M0 kilograms Initial mass of water column total P for the time step 

M14 kilograms Mass of water column total P at the end of the 14-day time step

Mavg kilograms Average mass of total P in the water column over the time step 

Mrec_pH kilograms Mass of P recycled from the sediment, pH-dependent mechanism

Mrec T kilograms Mass of P recycled from the sediment, temperature-dependent mechanism

N - Number of observed samples used for calculation

pH - Modeled pH

pHobs - Observed pH

pHsite - Site-specific observed pH

Pin parts per billion Inflow total P concentration

Pna parts per billion Non-algal P concentration

Palg parts per billion Algal P concentration

Palg,obs parts per billion Algal P concentration, biweekly observations

PT parts per billion Total P concentration

PT,obs parts per billion Total P concentration, biweekly observations

PSS parts per billion Steady state P concentration

Qin cubic meters per day Lake inflow rate

Qout cubic meters per day Lake outflow rate

ralg - Ratio of algal P to total P

rna - Ratio of non-algal P to total P

S kilograms Mass of total P in the sediment reservoir 

Sbur kilograms Mass of P buried and removed from the active sediment reservoir 

Sss kilograms Steady state mass of sediment P calculated at each time step 

S0 kilograms Initial mass of sediment P for the time step 

S14 kilograms Mass of total P in the sediment at the end of the 14-day time step

Savg kilograms Average mass of total P in the sediment reservoir over the time step 
29

Table 5. Definitions of symbols used in this report—continued

[For definitions and values of parameters, see table 4. All mass, concentration, and pH variables (except pHsite) in the water 

column are lakewide averaged quantities. -, no units; P, phosphorus]

Symbol Units Description

Sinit kilograms Initial size of sediment reservoir calculated with original starting concentration 1,000 mg P/kg sediment

Sinit new kilograms Initial size of sediment reservoir recalculated with starting concentration 335 mg P/kg sediment

Tlake degrees Celsius Lake temperature 

t day Time

V cubic meters Mean volume of the lake over the time step

YP

′

 milligrams phosphorus 

per milligram chlorophyll a

Minimum cell quota required for photosynthesis to occur

z meter Average lake depth over the time step

∆t day Time step

Mobs kilograms Change in mass stored in the lake over the time step, biweekly observations

Π/Πmax - Photosynthesis rate at a particular depth divided by maximum photosynthesis rate

Chlorophyll a Submodel

The chlorophyll a submodel is based on a 

mass balance in which it is assumed that chlorophyll a is always at its maximum potential value 

based on phosphorus concentration and water 

temperature averaged over the time step: 

dB

dt = �GROWTH − RESP − SETTLE − Qout

V � × B = 0,

(10) 

where

GROWTH = KGFLFTFP, (11) 

RESP = KrespFT, and (12) 

SETTLE = ualg

z . (13) 

In the model, growth depends on water temperature, nutrient availability (phosphorus only), 

and light, whereas respiration depends on water 

temperature only (Walker, 2001). No mechanism 

for algal mortality that would bring on a midseason bloom decline is built into the model; therefore, the model overpredicts chlorophyll a

concentrations during and after August in most 

years, which leads to overpredicted pH and total 

phosphorus recycled from sediment. In the model, chlorophyll a concentrations decrease in the 

autumn, when growth ceases largely because of 

temperature decreases, and settling and outflow 

dominate. The mass balance (eq. 10) can be rewritten in terms of model variables with 

temperature, nutrient, and light dependencies as:

KGFLFTFPB − �KrespFT + ualg

z + QoutV � B = 0 (14) 

When lake temperature is greater than 14°C, 

the algal growth rate in the model increases linearly with lake temperature between Tmin and 

Tmax (Bowie and others, 1985), so there is a 

steady increase in growth rate as lake temperatures increase followed by a steady decrease in 

growth rate as lake temperatures decrease:
30

FT = �

Tlake−Tmin

Tmax−Tmin

� (15) 

Tmin and Tmax were determined empirically

to be 14 and 20°C, respectively. When the lake 

temperature is less than 14°C, the growth rate is 

set to zero and the chlorophyll a concentration is 

set to a baseline value of 8 ppb. When the lake 

temperature is greater than 20°C, the temperature 

adjustment factor, FT, is equal to 1. 

The light limitation factor is based on an empirically determined equation relating 

photosynthesis (Π) to light intensity (I) at a particular depth (Talling, 1957): 

Π

Πmax

=

I

IK

�1+� I

IK

�

2

�

1/2 (16) 

The light limitation factor in the model 

scales this relation to the entire water column by 

substituting the depth-averaged light intensity for 

the light intensity at a particular depth, and then 

multiplies by the photoperiod to represent integration over 24 hours: 

FLAVG = �

Imean

IK

�1+�Imean

IK �

2

�

1/2� fp (17) 

Phosphorus limitation in the model is based 

on the results of Senft (1978) and Forsberg and 

Shapiro (1981). The limitation can be expressed 

in terms of a limiting cell quota and the ratio of 

chlorophyll a concentration to total phosphorus 

concentration as (using model notation): 

FPFS = (1 − YP

′ B

PT

), (18) 

where YP

′ is the minimum cell quota (mg P per 

mg chlorophyll a) required for photosynthesis to 

occur. In the model, the average, constant cell 

quota YP is substituted for YP

′

, so that the phosphorus limitation becomes the percent of 

available (nonalgal) phosphorus in the water column: 

FP00 = �1 − YP

B

PT

� = Pna

PT

. (19) 

pH Submodel

The pH submodel is a linear regression that 

determines pH as a function of the natural logarithm of chlorophyll a concentration and the 

Julian day (J): 

pH = 7.93 + 0.534 ln (B) − 0.006 max (J, 200). (20) 

The regression quantifies the high correlation between pH and chlorophyll a during the typically 

major bloom months of June–July (R2

 = 0.87; 

Walker, 2001), the result of photosynthesis removing dissolved CO2 from the water column. 

The regression coefficients are determined empirically as part of the calibration process. 

Alternative Solution Method

The water column and sediment phosphorus 

submodels were solved in the original TMDL 

model (and the R-script version of the original 

model) using the semi-analytical approach, in 

which the variables are assumed to respond with 

exponential decay over the time step in response 

to instantaneous changes in the boundary conditions and forcing functions at the beginning of 

each time step. Vollenweider (1975) proposed 

this solution to describe situations where the response to baseline loading could be described by 

the steady-state solution and the response to one 

or more anthropogenic loadings could be described by adding exponential decay terms to the 

steady-state solution. However, Vollenweider 

(1975) used an annual or seasonal time step as 

appropriate, and the steady-state solution was not 

recalculated at the beginning of each time step.

Given the shorter 2-week time step of the TMDL

model, the advantage of simplicity that the Vollenweider-based semi-analytical approach has for 

models with long time steps is reduced. 

We developed a more efficient conventional 

numerical solution method that made it easier to 

create and modify a program in a logical and easy 

to follow sequence. A model that incorporated all 

the same submodels, but solved the water column 

and sediment phosphorus equations (which include a time derivative) using a Runge-Kutta 
31

order 4 numerical method (Burden and Faires, 

1997), was developed as an R-script for comparison to the existing model. We quantitatively 

compared the concentrations of total phosphorus, 

sediment phosphorus, chlorophyll a, and pH simulated by the semi-analytical and Runge-Kutta 

models (table 6). The Runge-Kutta model reproduced the model results reported in Walker 

(2001) well, and became the basis for all the 

modified versions of the model used in this study. 

The Runge-Kutta model was used instead of the 

semi-analytical version because of reduced computational effort (due to removal of iterative 

solution convergence), confidence in using a traditional solution method, and simplified 

evaluation of alternative recycling mechanisms.

Model Calibration

As part of the model verification, we duplicated the process of calibration described in 

Walker (2001). This process proceeds in four 

steps: (1) calibration of the regression coefficients in the pH submodel using measured 

chlorophyll a concentrations from June through 

August 1991–1998, (2) calibration of the algal 

settling velocity using observed chlorophyll a

concentrations from June through July 1994–

1998, (3) calibration of the nonalgal phosphorus 

loss rate, which primarily represents the sedimentation of particulate phosphorus, based on observed lake phosphorus concentrations from

March through May 1991–1998, when the algal 

compartment is assumed to be minimal, and (4) 

calibration of the parameters controlling phosphorus recycling using observed pH and 

concentrations of total phosphorus and chlorophyll a from May through July 1994–1998. 

Step 1—pH Regression Model Parameters

The multiple regression was carried out in 

Excel using the chlorophyll a concentration and 

pH observations from the model. These datasets 

included June–August 1991–1998 and included 

52 data points. The coefficient of determination 

R2

 and standard error determined during the duplication of the calibration compared well to 

those provided in the TMDL model documentation (table 7; Walker, 2001). The regression in 

Walker (2001) was based on at least 54 data 

points, suggesting that our minimally different 

statistics result from the slightly shorter dataset 

provided with the model than was used in the 

original calibration. 

Table 7. Comparison of pH model coefficients and 

statistics between the original TMDL model for Upper Klamath Lake, Oregon, and the model 

recalibration.

Quantity

Original 

value

Recalibration 

value p-value

pH model: intercept 7.931 7.829 1.973E-19

pH model: chl a

coefficient 0.5338 0.5514 1.595E-17

pH model: Julian 

day coefficient -0.006 -0.006 0.0169586

Coefficient of determination R2 0.80 0.78

Standard error 0.24 0.25

Step 2—Algal Settling Rate

We used trial and error to duplicate the algal 

settling rate (ualg) calibration. Parameter values 

were assumed to be between 0.001 and 0.2 m/d 

Table 6. Comparison of simulated model variables 

obtained from the original, semi-analytical and 

Runge-Kutta models used in the review of the original 

model used to establish a TMDL for Upper Klamath 

Lake, Oregon. 

[Bias is calculated as the mean value of the semianalytical result minus the Runge-Kutta result at each 

time step. R, correlation coefficient; RMSE, root mean 

square error; NS, Nash-Sutcliffe statistic. P; phosphorus; 

kg, kilogram; ppb, parts per billion]

Variable R Bias RMSE NS

Water column P 1.00 0.86 ppb 2.44 ppb 1.00

Sediment P 1.00 -2.22 kg 2.47 kg 1.00

Chlorophyll a 1.00 0.89 ppb 2.19 ppb 1.00

pH 1.00 0.00 0.01 1.00
32

(Bowie and others, 1985). The best value was 

determined by correlating chlorophyll a as calculated by equation (14) (using observed 

phosphorus concentrations) with observed chlorophyll a concentrations in June–July 1994–1998. 

The sum of squared errors was minimized at 

ualg = 0.059 m/d, in close agreement with the 

original value (0.05 m/d) provided in the model 

documentation. Using a value of ualg = 0.059

m/d, the correlation coefficient R between the 

simulated and observed chlorophyll a values was 

0.79 (p<0.0001). 

Step 3—Nonalgal Phosphorus Loss Rate

The original model documentation does not 

specify the method used to calibrate the nonalgal

phosphorus loss rate Kna, but Kann and Walker 

(2001) discuss the general approach. The calibration of the nonalgal phosphorus loss rate uses the

biweekly mass balance of Kann and Walker 

(2001) during March–May, when recycling from 

the sediments is assumed to be nonexistent and 

total phosphorus in the water column is assumed 

to be predominantly nonalgal. Under these assumptions, the mass lost to sedimentation can be

calculated by rearranging equation (1). This provides a set of “observations” of mass lost to 

sedimentation, which was compared to a set of 

calculated values of mass lost to sedimentation. 

The calculated values were determined by assuming a value of Kna and calculating the mass lost

to sedimentation from lake total phosphorus concentration and volume at each time step: 

Msed = KnaVPTΔt. (21) 

We used a least-squares fit as a calibration 

criterion in the absence of other information and 

determined by trial and error that the value of 

Kna that resulted in the lowest sum of squared 

errors was 0.015/d, compared to the value of 

0.012/d that is supplied in Walker (2001). Using 

either value, the correlation coefficient between 

the “observed” and calculated values of mass lost

to sedimentation is the same (and low: R=0.19), 

and the regression coefficient is not significant 

(p=0.35), indicating that for practical purposes 

the regression model does not have predictive 

capability. This parameter is difficult to calibrate, 

likely because the model does not distinguish between different forms of nonalgal P. Particleassociated nonalgal phosphorus has a different 

loss rate from dissolved nonalgal P. Both values 

of Kna result in lower average values and smaller

standard deviations of Msed than characterize the 

observations (table 8). Even so, the model assumes that the loss rate determined in this way is 

applicable for all seasons and all nonalgal phosphorus concentrations throughout the year (Kann 

and Walker, 2001). 

Table 8. Comparison of descriptive statistics between observed and calculated non-algal mass lost 

to sedimentation during March–May 1991–1998

Mass lost to sedimentation

Average 

(kilograms) 

Standard 

Deviation 

(kilograms)

Observed (from mass balance) 9008 6608

Calculated with Kna=0.012 d-1 

(original model value) 6781 1813

Calculated with Kna=0.015 d-1 

(recalibrated model value) 8476 2267

Step 4—Internal Recycling Parameters

The equation for the mass of phosphorus 

transferred from the sediments to the water column (“recycled,” Mrec_pH) at any time step is an 

empirically derived function of the bottom sediment phosphorus concentration and the water 

column pH: 

Mrec_pH = SavgKRProb(pH > pH∗)∆t , (22) 

where

Prob(pH > pH∗) = 1 − CNDF �

pH∗−pH

σobs

� (23) 

and CNDF is the cumulative standard normal distribution of the expression in the square brackets. 

Equation 22 requires two calibration parameters: pH∗ and KR. The parameter pH∗ is the 

calibrated threshold at which the calibrated recycle rate KR is halved. At pH values much higher 

than pH*

, Prob(pH > pH∗) approaches 1, and at 
33

pH values much less than pH*, Prob(pH > pH∗)

approaches zero. The best values of these calibration parameters were those that minimized the 

following objective function:

O(PT, B, pH) =

min �

∑�PT,obs−PT�

2

∙∑(Bobs−B)2∙∑(pHobs−pH)2

N3 �, (24) 

where the summation is over all May–July 1994–

1998 data. The calibration period was restricted 

to May–July, when bloom growth is most rapid, 

and to avoid later summer months, when the 

model overpredicts chlorophyll a. This approach 

relies on the model structure and calibrated parameters from steps 1–3 being accurate. The 

parameters supplied with the original model were 

pH∗ = 9.1 and KR=0.53/yr. The values obtained 

in the process of duplicating the calibration with 

the datasets supplied with the model were 

pH∗ = 9.0 and KR=0.47/yr. While duplicating the 

calibration, however, it was noted that if the objective function (eq. 24) was modified to

incorporate May–July data for 1991–1998, instead of just for 1994–1998, the parameters that 

minimized the objective function were pH∗ = 9.1

and KR = 0.52/yr, much closer to those supplied 

with the original model. One explanation, therefore, for the differences between the parameters 

obtained when the calibration was duplicated and 

the original values is that the original calibration 

may have been performed using more years of 

data than was documented in Walker (2001). 

Comparison of Original to Recalibrated Model Performance

The recalibrated model differed from the 

original model in two ways: The calibrated recycling rate was 11 percent lower (0.47 versus 0.53 

per year), resulting in lower peaks of phosphorus 

in the water column during the time of year when 

recycling was important (fig. 7A). Because chlorophyll a is in equilibrium with available 

phosphorus at each time step, the lower watercolumn phosphorus translated to lower peaks in 

chlorophyll a (fig. 7B). The recalibrated algal settling velocity was 15 percent higher (0.059 versus

0.05 m/d) than the original, which acted to reinforce the lower peaks in chlorophyll a; lower 

peaks in chlorophyll a resulted in lower pH (fig. 

7C), which reinforced the lower recycling rate.

The original model was validated using data 

from all months in 1991–1993 and August–April 

in 1994–1998, which comprised the data set not 

used in the calibration. Although the documentation states a validation period, Walker (2001) 

presents no performance statistics for the final 

calibrated model. To fill this gap, we calculated 

several standard model statistics from simulations 

produced with the original model and the parameter values supplied with it. We compared these 

statistics to the same statistics calculated from 

simulations produced with the recalibrated model, which used parameter values determined by 

duplicating the entire calibration process and 

which differed in some cases from those supplied 

with the original model. This process was a 

means of quantifying the sensitivity of the model 

to the calibration parameters. The statistics considered are the correlation coefficient (R), the 

coefficient of determination (R2

), the mean of 

simulated minus observed values (Bias), the Nash 

Sutcliffe statistic (NS; Nash and Sutcliffe, 1970), 

and the root mean square error (RMSE). 

For all variables, the recalibrated model resulted in a slight improvement in model statistics: 

smaller bias, smaller RMSE, and increased NS 

(compare model versions O and A, tables 9–11). 

Overall, both the original and recalibrated model 

simulated total phosphorus concentration more 

successfully than they simulated chlorophyll a

concentration. Both models showed a negative 

bias for all variables, indicating that the models 

tended to underpredict the observed values during 

the validation period. Neither model performed 

well based on standard measures: R2

 varied from 

0.35 for total phosphorus to 0.18 for chlorophyll 

a in the recalibrated model; the NS statistic was 

negative or near zero for all variables in both 

models, indicating that the mean squared error 

was greater than or close to the variance in the 

dataset.
34

Figure 7. Diagrams showing simulated and observed water column (A) total phosphorus concentration, (B) chlorophyll a concentration, and (C) pH, Upper Klamath Lake, Oregon, April 1991–September 1998. Model version O is 

the original model with the original calibration parameters, and model version A is the recalibrated version of the 

original model.
35

Table 9. Comparison of performance statistics for total phosphorus obtained over a validation period including all months in 1991–1993 and August–

April in 1994–1998 from model versions discussed in the report.

[The statistics are calculated over a validation period including all months in 1991–1993 and August–April in 1994–1998 (N=83). R=correlation coefficient; 

R2

=coefficient of determination; RMSE=root mean square error; NS=Nash-Sutcliffe statistic; EQ=Chlorophyll a concentration evaluated by setting dB/dt to 0; 

RK=Chlorophyll a concentration evaluated using Runge-Kutta approach; INIT=Initial size of sediment reservoir calculated with original starting concentration 

1,000 mg P/kg sediment; INIT_NEW=Initial size of sediment reservoir recalculated with starting concentration 335 mg P/kg sediment; REC_pH=Model relies on 

pH-dependent recycle mechanism; REC_T=model relies on temperature-dependent recycle mechanism; P_00=Phosphorus limitation factor calculated as the percent of non-algal P; P_MM=Phosphorus limitation factor calculated by Michaelis-Menten kinetics; L_AVG=Light limitation factor calculated using the average 

light intensity; L_DI=Average light limitation factor calculated by depth-integrating the light limitation equation. Model characteristics and parameters are defined 

in the text.]

Model Characteristics Model Calibration Parameters Performance Statistics

Model 

version

Compare 

to 

version

Chlorophyll a

model

Initial 

sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light

limitation

ualg

(m/d)

KR

(per yr) pH*

KT

(per yr)

T*

(°C)

khalf

(ppm) R R2

Bias 

(ppb)

RMSE 

(ppb) NS

Original model with original calibration

O EQ INIT REC_pH P_00 L_AVG 0.05 0.53 9.1 NA NA NA 0.56 0.32 -48.6 109.6 -0.58

O' O EQ INIT_NEW REC_pH P_00 L_AVG 0.05 1.33 9.0 NA NA NA 0.60 0.36 -23.3 83.7 0.08

Recalibration of original model

A O EQ INIT REC_pH P_00 L_AVG 0.059 0.47 9.0 NA NA NA 0.59 0.35 -24.7 86.9 0.01

A' A EQ INIT_NEW REC_pH P_00 L_AVG 0.059 1.24 8.9 NA NA NA 0.61 0.38 -10.4 75.1 0.26

Modifications to light and phosphorus limitation

B A EQ INIT REC_pH P_MM L_DI 0.035 0.32 8.6 NA NA 0.06 0.58 0.34 -23.3 82.8 0.10

B' B EQ INIT_NEW REC_pH P_MM L_DI 0.035 1.05 8.7 NA NA 0.06 0.60 0.35 -23.4 83.1 0.10

Temperature-dependent recycling

C B EQ INIT REC_T P_MM L_DI 0.035 NA NA 0.54 12.8 0.06 0.59 0.35 -31.8 90.7 -0.08

C' C EQ INIT_NEW REC_T P_MM L_DI 0.035 NA NA 1.62 12.8 0.06 0.60 0.36 -25.1 84.8 0.06

Runge-Kutta solution of chlorophyll a

D B RK INIT REC_pH P_MM L_DI 0.035 0.34 8.1 NA NA 0.06 0.38 0.15 -60.0 114.9 -0.73

D' D RK INIT_NEW REC_pH P_MM L_DI 0.035 1.18 8.1 NA NA 0.06 0.39 0.15 -52.6 109.9 -0.59
36

Table 10. Comparison of performance statistics for chlorophyll a obtained over a validation period including all months in 1991–1993 and August–

April in 1994–1998 from model versions discussed in the report.

[The statistics are calculated over a validation period including all months in 1991-1993 and August-April in 1994-1998 (N=83). R=correlation coefficient; 

R2

=coefficient of determination; RMSE=root mean square error; NS=Nash-Sutcliffe statistic; EQ=Chlorophyll a concentration evaluated by setting dB/dt to 0; 

RK=Chlorophyll a concentration evaluated using Runge-Kutta approach; INIT=Initial size of sediment reservoir calculated with original starting concentration 

1,000 mg P/kg sediment; INIT_NEW=Initial size of sediment reservoir recalculated with starting concentration 335 mg P/kg sediment; REC_pH=Model relies on 

pH-dependent recycle mechanism; REC_T=model relies on temperature-dependent recycle mechanism; P_00=Phosphorus limitation factor calculated as the percent of non-algal P; P_MM=Phosphorus limitation factor calculated by Michaelis-Menten kinetics; L_AVG=Light limitation factor calculated using the average 

light intensity; L_DI=Average light limitation factor calculated by depth-integrating the light limitation equation. Model characteristics and parameters are defined 

in the text.] 

 Model Characteristics Model Calibration Parameters Performance Statistics

Model 

version

Compare 

to 

version

Chlorophyll a

model

Initial

sediment 

reservoir

Recycle 

mechanismPhosphorus 

limitation

Light 

limitation

ualg

(m/d)

KR

(per yr) pH*

KT

(per yr)

T*

(°C)

khalf

(ppm) R R2

Bias 

(ppb)

RMSE 

(ppb) NS

Original model with original calibration

O EQ INIT REC_pH P_00 L_AVG 0.05 0.53 9.1 NA NA NA 0.42 0.17 -56.5 129.6 -1.77

O' O EQ INIT_NEW REC_pH P_00 L_AVG 0.05 1.33 9.0 NA NA NA 0.43 0.18 -45.7 115.0 -1.18

Recalibration of original model

A A' EQ INIT REC_pH P_00 L_AVG 0.059 0.47 9.0 NA NA NA 0.42 0.18 -41.1 112.3 -1.08

A' A EQ INIT_NEW REC_pH P_00 L_AVG 0.059 1.24 8.9 NA NA NA 0.43 0.18 -34.7 103.9 -0.78

Modifications to light and phosphorus limitation

B A EQ INIT REC_pH P_MM L_DI 0.035 0.32 8.6 NA NA 0.06 0.40 0.16 -51.1 127.5 -1.68

B' B EQ INIT_NEW REC_pH P_MM L_DI 0.035 1.05 8.7 NA NA 0.06 0.40 0.16 -52.0 129.6 -1.77

Temperature-dependent recycling

C B EQ INIT REC_T P_MM L_DI 0.035 NA NA 0.54 12.8 0.06 0.36 0.13 -59.2 143.3 -2.39

C' C EQ INIT_NEW REC_T P_MM L_DI 0.035 NA NA 1.62 12.8 0.06 0.37 0.14 -55.7 138.4 -2.16

Runge-Kutta solution of chlorophyll a

D B RK INIT REC_pH P_MM L_DI 0.035 0.34 8.1 NA NA 0.06 0.48 0.23 -61.0 114.4 -1.16

D' D RK INIT_NEW REC_pH P_MM L_DI 0.035 1.18 8.1 NA NA 0.06 0.48 0.23 -57.2 112.3 -1.08
37

Table 11. Comparison of performance statistics for pH obtained over a validation period including all months in 1991–1993 and August–April in 1994–

1998 from model versions discussed in the report.

[The statistics are calculated over a validation period including all months, 1991–1993 and August–April, 1994–1998 (N=83). R=correlation coefficient; 

R2

=coefficient of determination; RMSE=root mean square error; NS=Nash-Sutcliffe statistic; EQ=Chlorophyll a concentration evaluated by setting dB/dt to 0; 

RK=Chlorophyll a concentration evaluated using Runge-Kutta approach; INIT=Initial size of sediment reservoir calculated with original starting concentration 

1,000 mg P/kg sediment; INIT_NEW=Initial size of sediment reservoir recalculated with starting concentration 335 mg P/kg sediment; REC_pH=Model relies on 

pH-dependent recycle mechanism; REC_T=model relies on temperature-dependent recycle mechanism; P_00=Phosphorus limitation factor calculated as the percent of non-algal P; P_MM=Phosphorus limitation factor calculated by Michaelis-Menten kinetics; L_AVG=Light limitation factor calculated using the average 

light intensity; L_DI=Average light limitation factor calculated by depth-integrating the light limitation equation. Model characteristics and parameters are defined 

in the text.]

Model characteristics Model calibration parameters Performance statistics

Model 

version

Compare 

to 

version

Chlorophyll a

model

Initial

sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

ualg

(m/d)

KR

(per yr) pH*

KT

(per yr)

T*

(°C)

khalf

(ppm) R R2 Bias RMSE NS

Original model with original calibration

O EQ INIT REC_pH P_00 L_AVG 0.05 0.53 9.1 NA NA NA 0.57 0.32 -0.195 0.716 -0.09

O' O EQ INIT_NEW REC_pH P_00 L_AVG 0.05 1.33 9.0 NA NA NA 0.56 0.32 -0.175 0.703 -0.05

Recalibration of original model

A A' EQ INIT REC_pH P_00 L_AVG 0.059 0.47 9.0 NA NA NA 0.56 0.31 -0.155 0.700 -0.04

A' A EQ INIT_NEW REC_pH P_00 L_AVG 0.059 1.24 8.9 NA NA NA 0.55 0.31 -0.143 0.694 -0.02

Modifications to light and phosphorus limitation

B A EQ INIT REC_pH P_MM L_DI 0.035 0.32 8.6 NA NA 0.06 0.57 0.33 -0.161 0.704 -0.05

B' B EQ INIT_NEW REC_pH P_MM L_DI 0.035 1.05 8.7 NA NA 0.06 0.58 0.33 -0.157 0.699 -0.04

Temperature-dependent recycling

C B EQ INIT REC_T P_MM L_DI 0.035 NA NA 0.54 12.8 0.06 0.57 0.32 -0.175 0.716 -0.09

C' C EQ INIT_NEW REC_T P_MM L_DI 0.035 NA NA 1.62 12.8 0.06 0.57 0.32 -0.168 0.712 -0.07

Runge-Kutta solution of chlorophyll a

D B RK INIT REC_pH P_MM L_DI 0.035 0.34 8.1 NA NA 0.06 0.64 0.41 -0.247 0.622 0.18

D' D RK INIT_NEW REC_pH P_MM L_DI 0.035 1.18 8.1 NA NA 0.06 0.65 0.42 -0.235 0.612 0.21

Several parameters in the model were calibrated using data 

from only May, June, and July, when the bloom is in a rapid 

growth phase; the performance of the model was tested over 

those months alone. Tables 12–14 compare model statistics for 

May–July 1991–1993 (years not used in the calibration). Model 

performance was better during this limited validation period than 

during the longer validation period included in the statistics presented in tables 9–11 for all variables and both models, as 

indicated by larger R2 values, smaller bias, and positive NS statistics. R2

 was 0.89 or 0.86 for total phosphorus, 0.55 or 0.56 for 
38

chlorophyll a (version O and A, respectively), and 0.56 (both 

versions) for pH. The NS statistic was 0.88 and 0.81 for total 

phosphorus, 0.28 or 0.40 for chlorophyll a, and 0.19 or 0.22 for 

pH (version O and A, respectively). The bias was positive for 

total phosphorus and negative for chlorophyll a and pH (both 

versions), indicating that over this limited validation period that 

included only spring through mid-summer months, the model 

overpredicted total phosphorus, but still underpredicted chlorophyll a and pH. Based on these statistics, the recalibrated model

(version A) simulated total phosphorus concentration slightly 

worse than, and pH and chlorophyll a concentration slightly better than, the original model and calibration (version O), but the 

differences in the performance statistics were small. 

Table 12. Comparison of performance statistics for total phosphorus over a validation period including May–July 1991–1993 obtained from various 

model versions discussed in the report.

[The statistics are calculated over a validation period including May–July 1991–1993 (N=16). R=correlation coefficient; R2=coefficient of determination; 

RMSE=root mean square error; NS=Nash-Sutcliffe statistic; EQ=Chlorophyll a concentration evaluated by setting dB/dt to 0; RK=Chlorophyll a concentration 

evaluated using Runge-Kutta approach; INIT=Initial size of sediment reservoir calculated with original starting concentration 1,000 mg P/kg sediment; 

INIT_NEW=Initial size of sediment reservoir recalculated with starting concentration 335 mg P/kg sediment; REC_pH=Model relies on pH-dependent recycle 

mechanism; REC_T=model relies on temperature-dependent recycle mechanism; P_00=Phosphorus limitation factor calculated as the percent of non-algal P; 

P_MM=Phosphorus limitation factor calculated by Michaelis-Menten kinetics; L_AVG=Light limitation factor calculated using the average light intensity; 

L_DI=Average light limitation factor calculated by depth-integrating the light limitation equation. Model characteristics and parameters are defined in the text.] 

Model characteristics Model calibration parameters Performance statistics

Model 

version

Compare 

to 

version

Chlorophyll a

model

Initial

sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

ualg

(m/d)

KR

(per yr) pH*

KT

(per yr)

T*

(°C)

khalf

(ppm) R R2

Bias 

(ppb)

RMSE 

(ppb) NS

Original model with original calibration

O EQ INIT REC_pH P_00 L_AVG 0.05 0.53 9.1 NA NA NA 0.94 0.89 8.3 28.8 0.88

O' O EQ INIT_NEW REC_pH P_00 L_AVG 0.05 1.33 9.0 NA NA NA 0.92 0.85 9.7 35.9 0.82

Recalibration of original model

A O EQ INIT REC_pH P_00 L_AVG 0.059 0.47 9.0 NA NA NA 0.93 0.86 14.9 36.5 0.81

A' A EQ INIT_NEW REC_pH P_00 L_AVG 0.059 1.24 8.9 NA NA NA 0.90 0.81 13.7 43.1 0.74

Modifications to light and phosphorus limitation

B A EQ INIT REC_pH P_MM L_DI 0.035 0.32 8.6 NA NA 0.06 0.94 0.88 17.4 39.3 0.78

B' B EQ INIT_NEW REC_pH P_MM L_DI 0.035 1.05 8.7 NA NA 0.06 0.96 0.92 21.6 36.6 0.81

Temperature-dependent recycling

C B EQ INIT REC_T P_MM L_DI 0.035 NA NA 0.54 12.8 0.06 0.95 0.91 25.8 44.1 0.72

C' C EQ INIT_NEW REC_T P_MM L_DI 0.035 NA NA 1.62 12.8 0.06 0.95 0.90 26.6 45.6 0.71

Runge-Kutta solution of chlorophyll a

D B RK INIT REC_pH P_MM L_DI 0.035 0.34 8.1 NA NA 0.06 0.94 0.88 -6.8 34.4 0.83

D' D RK INIT_NEW REC_pH P_MM L_DI 0.035 1.18 8.1 NA NA 0.06 0.94 0.88 -14.3 35.5 0.82
39

Table 13. Comparison of performance statistics for chlorophyll a obtained over a validation period including May–July 1991–1993 from model versions 

discussed in the report.

[The statistics are calculated over a validation period including May–July 1991–1993 (N=16). R=correlation coefficient; R2=coefficient of determination; 

RMSE=root mean square error; NS=Nash-Sutcliffe statistic; EQ=Chlorophyll a concentration evaluated by setting dB/dt to 0; RK=Chlorophyll a concentration 

evaluated using Runge-Kutta approach; INIT=Initial size of sediment reservoir calculated with original starting concentration 1,000 mg P/kg sediment; 

INIT_NEW=Initial size of sediment reservoir recalculated with starting concentration 335 mg P/kg sediment; REC_pH=Model relies on pH-dependent recycle 

mechanism; REC_T=model relies on temperature-dependent recycle mechanism; P_00=Phosphorus limitation factor calculated as the percent of non-algal P; 

P_MM=Phosphorus limitation factor calculated by Michaelis-Menten kinetics; L_AVG=Light limitation factor calculated using the average light intensity; 

L_DI=Average light limitation factor calculated by depth-integrating the light limitation equation. Model characteristics and parameters are defined in the text.R, 

correlation coefficient; R2

, coefficient of determination; RMSE, root mean square error; NS, Nash-Sutcliffe statistic. Model characteristics and parameters are defined in the text.]

 Model characteristics Model calibration parameters Performance statistics

Model 

version

Compare 

to 

version

Chlorophyll a

model

Initial

sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

ualg

(m/d)

KR

(per yr) pH*

KT

(per yr)

T*

(°C)

khalf

(ppm) R R2

Bias 

(ppb)

RMSE 

(ppb) NS

Original model with original calibration

O EQ INIT REC_pH P_00 L_AVG 0.05 0.53 9.1 NA NA NA 0.74 0.55 -48.9 85.4 0.28

O' O EQ INIT_NEW REC_pH P_00 L_AVG 0.05 1.33 9.0 NA NA NA 0.75 0.56 -48 82.3 0.33

Recalibration of original model

A A' EQ INIT REC_pH P_00 L_AVG 0.059 0.47 9.0 NA NA NA 0.75 0.56 -39.6 77.8 0.40

A' A EQ INIT_NEW REC_pH P_00 L_AVG 0.059 1.24 8.9 NA NA NA 0.75 0.56 -40.6 78.0 0.40

Modifications to light and phosphorus limitation

B A EQ INIT REC_pH P_MM L_DI 0.035 0.32 8.6 NA NA 0.06 0.77 0.59 -37.6 77.1 0.41

B' B EQ INIT_NEW REC_pH P_MM L_DI 0.035 1.05 8.7 NA NA 0.06 0.77 0.59 -33.4 76.2 0.43

Temperature-dependent recycling

C B EQ INIT REC_T P_MM L_DI 0.035 NA NA 0.54 12.8 0.06 0.75 0.56 -27.5 74.0 0.46

C' C EQ INIT_NEW REC_T P_MM L_DI 0.035 NA NA 1.62 12.8 0.06 0.75 0.56 -26.6 73.1 0.47

Runge-Kutta solution of chlorophyll a

D B RK INIT REC_pH P_MM L_DI 0.035 0.34 8.1 NA NA 0.06 0.81 0.65 11.6 61.4 0.63

D' D RK INIT_NEW REC_pH P_MM L_DI 0.035 1.18 8.1 NA NA 0.06 0.80 0.65 7.6 61.8 0.62
40

Table 14. Comparison of performance statistics for pH obtained over a validation period including May–July 1991–1993 from various model versions discussed in the report. 

[The statistics are calculated over a validation period including May-July, 1991–1993 (N=16). R=correlation coefficient; R2=coefficient of determination; 

RMSE=root mean square error; NS=Nash-Sutcliffe statistic; EQ=Chlorophyll a concentration evaluated by setting dB/dt to 0; RK=Chlorophyll a concentration evaluated using Runge-Kutta approach; INIT=Initial size of sediment reservoir calculated with original starting concentration 1000 mg P/kg sediment; INIT_NEW=Initial 

size of sediment reservoir recalculated with starting concentration 335 mg P/kg sediment; REC_pH=Model relies on pH-dependent recycle mechanism; 

REC_T=model relies on temperature-dependent recycle mechanism; P_00=Phosphorus limitation factor calculated as the percent of non-algal P; P_MM=Phosphorus 

limitation factor calculated by Michaelis-Menten kinetics; L_AVG=Light limitation factor calculated using the average light intensity; L_DI=Average light limitation 

factor calculated by depth-integrating the light limitation equation. Model characteristics and parameters are defined in the text.]

 Model characteristics Model calibration parameters Performance statistics

Model 

version

Compare 

to 

version

Chlorophyll a

model

Initial

Sediment 

Reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

ualg

(m/d)

KR

(per yr) pH*

KT

(per yr)

T*

(°C)

khalf

(ppm) R R2 Bias RMSE NS

Original model with original calibration

O EQ INIT REC_pH P_00 L_AVG 0.05 0.53 9.1 NA NA NA 0.75 0.56 -0.443 0.659 0.19

O' O EQ INIT_NEW REC_pH P_00 L_AVG 0.05 1.33 9.0 NA NA NA 0.73 0.54 -0.456 0.677 0.14

Recalibration of original model

A A' EQ INIT REC_pH P_00 L_AVG 0.059 0.47 9.0 NA NA NA 0.75 0.56 -0.417 0.644 0.22

A' A EQ INIT_NEW REC_pH P_00 L_AVG 0.059 1.24 8.9 NA NA NA 0.73 0.53 -0.435 0.666 0.17

Modifications to light and phosphorus limitation

B A EQ INIT REC_pH P_MM L_DI 0.035 0.32 8.6 NA NA 0.06 0.74 0.55 -0.320 0.589 0.35

B' B EQ INIT_NEW REC_pH P_MM L_DI 0.035 1.05 8.7 NA NA 0.06 0.78 0.60 -0.279 0.541 0.45

Temperature-dependent recycling

C B EQ INIT REC_T P_MM L_DI 0.035 NA NA 0.54 12.8 0.06 0.77 0.60 -0.300 0.551 0.43

C' C EQ INIT_NEW REC_T P_MM L_DI 0.035 NA NA 1.62 12.8 0.06 0.78 0.60 -0.299 0.549 0.43

Runge-Kutta solution of chlorophyll a

D B RK INIT REC_pH P_MM L_DI 0.035 0.34 8.1 NA NA 0.06 0.86 0.73 0.061 0.386 0.72

D' D RK INIT_NEW REC_pH P_MM L_DI 0.035 1.18 8.1 NA NA 0.06 0.86 0.74 0.043 0.382 0.73
41

The primary reason for poorer model statistics when months outside of May–July were 

included in the calculation is the inability of the 

model to replicate algal decline in August–

October. Simulated values of chlorophyll a and 

total phosphorus (of which algal phosphorus is 

the largest component) were much higher (as 

much as 4 times) than observed values during 

those months (fig. 7). A positive feedback inherent in the model contributes to the discrepancy—

because there is no modeled midseason bloom 

decline, higher-than-observed chlorophyll a values result in higher-than-observed pH values, 

according to the pH submodel (eq. 20). These

higher-than-observed pH values induce more recycling (eq. 2), which elevates the amount of 

bioavailable phosphorus in the water column, 

which feeds back into the growth term in the 

chlorophyll a model (eq. 11). 

Sensitivity of Calibration Parameters 

to Model Modifications

In the process of reviewing the model equations, we identified a few easily implemented 

improvements. These included modifications to 

the mathematical forms of limitation factors 

(light and phosphorus) and elimination of the assumption that chlorophyll a concentration is in 

instantaneous equilibrium with light, temperature, 

and nutrients. 

Updated information regarding the size and 

forms of the phosphorus reservoir in the bed sediments prompted investigation of two additional 

modifications. The first modification was a 

change in the initial size of the sediment phosphorus reservoir. The second was a change in the 

form of the phosphorus recycling term to simulate a mechanism of recycling that would depend 

more on temperature than on pH and likely would 

draw on the sizable organic phosphorus fraction 

as a source. Each of the considered modifications 

to the model required that the calibrations of the 

algal settling rate and (or) internal recycling be 

revisited. One of the modifications (temperaturedependent recycling) replaced the calibration parameters in the internal recycling component of 

the model with new parameters. We investigated 

the sensitivity of the calibration parameters to the 

various model modifications and compared the 

simulated results produced by the original and 

recalibrated model versions. 

Growth Limitation Factors

Phosphorus and light availability limit algal 

growth in the original model (section “Chlorophyll a Submodel”). The phosphorus limitation is 

calculated as the percentage of available (nonalgal) phophorus in the water column, and light 

limitation is calculated with an application of Talling's (1957) equation, substituting average light 

intensity in the water column for light at a particular depth. As part of the validation process, we 

replaced these equations with other known and 

accepted equations (Bowie and others, 1985) using a Michaelis-Menten kinetics approach for 

phosphorus limitation and a depth-integrated Talling equation approach for light limitation. 

In the original model, using the percentage of

nonalgal phosphorus as FP (eq. 19) can result in 

little phosphorus limitation (even at very low 

concentrations) when nonalgal phosphorus is a 

large fraction of the total phosphorus, and much 

phosphorus limitation (even at very high concentrations) when nonalgal phosphorus is a small 

fraction of the total. As an alternative, MichaelisMenten kinetics is a well-known and accepted 

functionality between a reaction rate and the concentration of a substrate that allows the reaction 

rate (in this case, algal growth rate) to go to zero

as the concentration of the substrate (in this case, 

nonalgal phosphorus concentration) goes to zero, 

and to approach a maximum value asymptotically 

as the concentration of the substrate becomes 

very large. The equation for FP presented in 

Bowie and others (1985), maintains MichaelisMenten kinetics for all values: 

FPMM = Pna

khalf+Pna

. (25) 
42

Equation 25 was substituted for equation 19 in a 

modified model, denoted version B in tables 9–

14. 

The light-limitation function in the original 

model provides the light-limitation factor for 

whatever depth is associated with the average 

light intensity (eq. 17). Alternatively, integrating 

the Talling function over depth provides the average light limitation over the entire lake depth, z: 

FLDI = fp

Ecz ln

⎣

⎢

⎢

⎡ I0

IK

+�1+�I0

IK

�

2

�

1

2

I0

IK

e−Ecz+�1+�I0

IK

e−Ecz�

2

�

1

2

⎦

⎥

⎥

⎤

(26) 

Equation 26 was substituted for equation 17 in 

version B of the model (tables 9–14). 

Because FP and FL are used in the chlorophyll a submodel, we performed a recalibration 

of the settling parameter, ualg (table 4), as well as 

calibration of the new parameter, the halfsaturation constant khalf, followed by a recalibration of the recycling parameters (steps 2 and 4 of 

the calibration process). Calibration of ualg and

the new parameter khalf was accomplished with

a trial-and-error approach by evaluating combinations of these two parameters within the range of 

0.001–0.2 m/d and 0.01–0.06 ppm, respectively, 

based on literature values (Bowie and others, 

1985). The effect of the modified limitation factors was to decrease the calibrated settling 

velocity ualg from 0.059 to 0.035 m/d, decrease 

the calibrated recycle rate KR from 0.47 to 

0.32/yr, and decrease the pH at which recycling is 

one-half the maximum (pH∗) from 9.0 to 8.6 (tables 9–14). These changes in the calibration 

parameters are not the result of mechanistic, process-based changes, but rather are a consequence 

of a change in the balance of terms in equation 14 

that was used to calibrate ualg in step 2 of the 

calibration process, which in turn changed the 

calibration of KR and pH∗ in step 4 of the calibration. 

The combined effect of the modifications to 

the phosphorus and light limitation factors was to 

decrease the annual peak in total phosphorus concentration relative to the original, recalibrated

model because the recycling rate was smaller; 

even so, the combination of lower algal settling 

and modification to growth limitation factors resulted in somewhat higher peaks in chlorophyll a

and pH (fig. 8). 

Model performance was not substantially 

changed by the modification of the phosphorus 

and light limitation expressions, either when statistics were calculated over the full validation 

period, when performance was particularly poor 

based on low R2

 values, high bias, and low NS 

values (tables 9–11), or over the limited May–

July validation period, when performance improved (tables 12–14).

Recycling Mechanism

Several mechanisms likely contribute to the 

internal loading of phosphorus from the sediments, including diffusion, desorption, diffusion 

enhanced by bioturbation, and metabolic cycling. 

It is impossible to know which of these mechanisms dominates at any given point in time or 

location in the lake, but given the high density of 

benthic macroinvertebrates, mechanisms involving benthic macroinvertebrates probably are 

important at some times and locations. Furthermore, these mechanisms are more likely to 

depend on temperature than on pH and would 

likely draw on the sizable organic phosphorus 

fraction as a source. A well-developed mathematical description of recycling by the benthic 

macroinvertebrate community is outside the 

scope of this study, but one of the dominant features of such a description would include a 

dependence on temperature. To investigate the 

sensitivity of the calibration of recycling parameters to the form of the recycling term, we 

developed a temperature-dependent recycling 

mechanism analogous to the temperature dependency of cyanobacterial growth in the model. We 

substituted the temperature-dependent recycling 

mechanism for the pH-dependent recycling 

mechanism in model version B, which, like the 

original model, incorporates an internal recycling
43

 Figure 8. Diagrams showing simulated and observed water column (A) total phosphorus concentration, (B) chlorophyll a concentration, and (C) pH in Upper Klamath Lake, Oregon, April 1991–April 1998. Model version A is the 

original model (recalibrated) and model version B incorporates modifications to the phosphorus and light limitations 

to algal growth.mechanism (eq. 22) that depends on pH, the percentage of the lake (represented by percentage of 

sampling sites) exceeding a threshold pH denoted pH* (the pH at which the recycle rate is one-half the maximum 

recycle rate, KR), and the average mass of phosphorus in the sediment reservoir. Thus, equation 22 is a mathematical representation of pH-dependent desorption of phosphorus from poorly crystalline iron oxides in the sediments. 

The form of equation 22 derives from a 

demonstrated dependence of internal recycling on 

pH during May–July (Walker, 2001). The dependence of internal recycling on temperature 

can be similarly demonstrated. The calibrated 

sedimentation terms (nonalgal phosphorus loss 

rate Kna and the algal settling velocity ualg from 

calibration steps 2 and 3) were substituted into a 

discretized version of the water column phosphorus mass balance (eq. 1), and the “observed” 

recycling at each time step was calculated as:

Mrec,obs = ∆Mobs + Mout,obs + �KnaPT,obsV∆t +

ualgPalg,obsA∆t� − Mload,obs (27) 
44

When these observations of Mrec were plotted against temperature, recycling was seen to 

increase with temperature, particularly above approximately 13°C (fig. 9). A relatively simple 

form for temperature-dependent recycling in 

which recycling is proportional to temperature 

above a threshold value describes the relation evident in figure 9, and also bears a functional 

resemblance to equation 22: 

Mrec_T =

�

Tlake−Tmin,rec

Tmax,rec−Tmin,rec

KTSavg∆t, when Tlake ≥ Tmin,rec

0, when Tlake < Tmin,rec

(28) 

in which the rate constant KT and the minimum 

temperature at which recycling can occur Tmin,rec

are calibration parameters. The calibration of 

these parameters proceeded, in analogy to the calibration of KR and pH∗, by trial and error. The 

values KT = 0.54/yr and Tmin,rec = 12.8°C were

the values that minimized the objective function 

given in equation 24 (model version C in tables 

9–14). 

Figure 9. Graph showing mass of phosphorus recycled from sediments, calculated from observations of water column phosphorus and external loads, Upper Klamath Lake, Oregon, May–July 1991–1998, as a function of the lake 

temperature. 

Replacing the pH-dependent recycling 

mechanism with a temperature-dependent recycling mechanism resulted in an increase in most 

of the peak phosphorus and chlorophyll a concentrations during the 7.5-year calibration period

(fig. 10). The model modified with this simple 

temperature-dependent recycling mechanism has 

generally slightly worse performance statistics 

overall than the preceding model (model version 

C compared to model version B in tables 9–14), 

particularly for chlorophyll a. Over the full validation period, R2

 changed from 0.34 to 0.35 and 

NS changed from 0.101 to -0.079, and over the 

limited validation period R2

 changed from 0.88 to 

0.91 and NS changed from 0.78 to 0.72 for total 

phosphorus (version B to version C). Over the 
45

full validation period, R2 changed from 0.16 to 

0.13 and NS changed from -1.68 to -2.39, and 

over the limited validation period R2

 changed 

from 0.59 to 0.56 and NS changed from 0.41 to 

0.46 for chlorophyll a (version B to version C). 

Given the major conceptual change that this modification represented, these changes are small, 

and indicate that neither the pH-dependent nor 

the temperature-dependent mechanism is uniquely suited to describing the observations; rather, 

any rate-dependent first-order term that generally 

increased through the spring and summer like pH 

and temperature would perform approximately 

equally well, and any such term probably would 

have generally poor performance statistics, particularly over the full validation period. 

Figure 10. Graphs showing simulated and observed water column (A) total phosphorus concentration, (B) chlorophyll a concentration, and (C) pH, Upper Klamath Lake, Oregon, April 1991–September 1998. Model version B

incorporates modifications to the phosphorus and light limitations to algal growth, and model version C has been 

modified to replace a pH-dependent recycling mechanism with a temperature-dependent recycling mechanism.
46

Numerical Solution of Chlorophyll a Equation

In this sensitivity analysis, we examined the 

implications of the assumption of equilibration of 

chlorophyll a concentration to temperature, light, 

and nutrients at each time step, which is, in effect, an assumption that the cyanobacterial bloom 

achieves maximum possible growth during each 

time step. The same chlorophyll a mass balance 

(eq. 10) was used, but the assumption that the 

time derivative was zero was not made. A 4th 

order Runge-Kutta numerical method was used to 

solve the differential equation in an analogous 

manner to the solution of the water column and 

sediment phosphorus equations. 

Model version B, which incorporated modifications to the phosphorus and light limitation 

factors was revisited using the Runge-Kutta solution for chlorophyll a. This version of the model 

was denoted version D (tables 9–14). The calibration of the internal recycling parameters (step 4) 

was repeated, with the result that the calibration 

parameter pH∗ was lower by one-half pH unit, 

indicating accelerated recycling as a function of 

pH, and that the newly calibrated recycle rate KR

was nearly unchanged (0.34/yr versus 0.32/yr, 

comparing version D to version B; tables 9–14). 

As a result, the modified model (version D) simulated higher baseline concentrations of water 

column phosphorus than version B; however, because chlorophyll a was no longer required to 

achieve its maximum value in the 2-week time 

interval, peaks in chlorophyll a were uniformly 

reduced over the 7.5-year calibration period 

(fig. 11). 

The model modified with the Runge-Kutta 

solution for chlorophyll a was the best performing model of all the models tested for all three 

variables during the limited validation period 

comprising May–July: highest correlation coefficients (0.88, 0.65, and 0.73 for total phosphorus, 

chlorophyll a, and pH, respectively), lowest 

RMSE, and NS closer to 1 (0.83, 0.63, and 0.72 

for total phosphorus, chlorophyll a, and pH, respectively (tables 12–14). During the full 

validation period, this model gave the best results 

for chlorophyll a (R2

=0.23) and pH (R2=0.41; 

tables 10 and 11), but it did not simulate total 

phosphorus as well as the other models that were 

tested (R2=0.15; table 9). The winter spikes in pH 

(fig. 11C) resulted from discontinuity in the pH 

regression model (eq. 20) at the end of the calendar year and can be eliminated in future model 

revisions. 
47

Figure 11. Graphs showing simulated and observed water column (A) total phosphorus concentration, (B) chlorophyll a concentration, and (C) pH, Upper Klamath Lake, Oregon, April 1991–September 1998. Model version B

incorporates modifications to the phosphorus and light limitations to algal growth, and model version D incorporates 

a nonequilibrium approach to the solution of the chlorophyll a mass balance.

Initial Sediment Reservoir

One of two initial conditions of the model is

the initial mass of total phosphorus in the sediment reservoir (the other being the initial mass of 

total phosphorus in the water column). The initial 

mass in the reservoir is calculated from

Sinit = ρbulkzsedAPsed. (29) 

The sediment bulk density ρbulk used was 0.12 

g/cm3 and the active sediment thickness zsed was

10 cm, based on Eilers and others (2001). These

values are consistent with more recent data (Simon and Ingle, 2011); however, the initial 

sediment phosphorus concentration used in the 

original model was 1,000 mg/kg, based on the 

single core analyzed by Eilers and others (2001). 

Because more abundant core data have become

available (Simon and Ingle, 2011), the average 
48

concentration of phosphorus in the sediments, 

Psed, was recalculated based on the new information. The new average sediment concentration, 

Psed_new, based on all 26 cores at depths to 10 cm 

was 335 mg/kg, only about one-third of that reported by Eilers and others (2001), making it 

necessary to investigate sensitivity of the model 

calibration to initial sediment phosphorus concentration. 

We ran each of the model versions again using the new initial condition for sediment 

phosphorus. Recalibration of algal settling velocity was not necessary, because the balance of 

terms used in that step did not change with this 

change to initial conditions. The calibration of the 

recycling parameters (step 4) was repeated, however, and the value of the calibration parameter 

pH∗ changed little—between 0 and 0.2 pH units 

(tables 9–14)—but the newly calibrated recycle 

rate KR (or KT in version C′) more than doubled 

in some cases (versions O′ and A′) or increased

by over threefold (versions B′, C′, and D′), indicating accelerated recycling at all values of pH. 

Even so, the effect of the lower initial sediment 

reservoir was to reduce peaks in all of total phosphorus, pH, and chlorophyll a during the 7.5-year 

calibration period relative to original initial condition (fig. 12, model version A as an example). 

Statistics indicated that the lower initial sediment 

reservoir generally did not greatly affect model 

performance, regardless of whether statistics 

were calculated over the full (tables 9–11) or limited (tables 12–14) validation period. Model 

performance increased marginally over the full 

validation period, particularly in the simulation of 

total phosphorus, as evidenced by small increases 

in R2 and NS, and decreases in RMSE (tables 9–

11). Model performance decreased marginally 

over the limited validation period (tables 12–14) 

in the simulation of total phosphorus, and changes between un-primed and primed versions of the 

model were small and mixed in the simulation of 

chlorophyll a and pH.
49

Figure 12. Graphs showing simulated and observed water column (A) total phosphorus concentration, (B) chlorophyll a concentration, and (C) pH, Upper Klamath Lake, Oregon, April 1991–September 1998. Model version A is 

the original model (recalibrated), and model version A′ is the original model recalibrated with a lower initial condition 

for the sediment phosphorus reservoir.

Long-Term Model Simulation and 

Implications for TMDL

The original model was used to address 

questions that preceded the implementation of the 

TMDL because of its ability to predict the response of water quality in Upper Klamath Lake 

to changes in management of the land and point 

sources of phosphorus in the Upper Klamath 

Lake drainage basin. The CSTR is well-suited to 

addressing these questions because it is not computationally intensive and can provide 

simulations of conditions many decades into the 

future, consistent with the time frame for response to management changes. The 2-week time 

step and the input data requirements of the model 

were well-matched to the biweekly sampling 

strategy and the collected constituents of the 
50

long-term monitoring program that provided 

boundary conditions, observations for calibration 

of parameters, and mass balance terms for validation. The simplicity of the algorithms also allows

insight into the factors that determine long-term 

system behavior that would be difficult to generate with a more complicated model.

The model predicts the response to a reduction in external loads of phosphorus to the lake. 

The original model included a “TMDL” option, 

which looped through the 7 years of input data

(mid-April 1991–mid-April 1998) until a steady 

state (equilibrium) was achieved without storing 

interim values. The steady state was defined as 

the point at which the simulated water column 

and sediment phosphorus over the 7-year period 

were essentially unchanged with additional cycles 

through the input data. This method provided the 

end point, but did not indicate how long it took to 

achieve steady state (how many 7-year cycles) or 

how concentrations evolved during the transition 

period. 

As part of the current validation process, we 

modified the R-script model to chain the 7 years 

of model input data to produce a continuous simulation of all the variables over 203 years (29 

seven-year cycles), which brought the simulations to a steady state. These simulations were 

obtained with all versions of the model discussed 

previously. A second set of simulations was obtained with all versions of the model, but after 

first reducing the external load boundary condition by 40 percent, as recommended in the Upper 

Klamath Lake TMDL (Oregon Department of 

Environmental Quality, 2002). Thus, the simulations were arranged into four groups of five 

simulations each, defined by the four combinations of original or updated initial sediment 

reservoir and 0 or 40 percent reduction in the external load boundary condition. 

Within a group of simulations, the behavior 

of each model can be summarized with the set of 

metrics compiled in table 15. The ratio of the recycling term to the sedimentation term in the 

water column phosphorus balance 

(Mrec/Msed , eq. 1) averaged over the first 7 

years of simulation was positively correlated with 

the average water column phosphorus concentration over the last 7 years of the simulation (fig. 

13A) and negatively correlated with the sediment 

reservoir averaged over the last 7 years of the 

simulation (fig. 13B). Given a set of boundary 

and initial conditions, more recycling relative to 

settling resulted in a steady state characterized by 

lower sediment phosphorus and higher water column phosphorus concentration as compared to 

less recycling relative to settling; thus, those two 

quantities were negatively correlated (fig. 13D). 

The ratio Mrec/Msed also was correlated with 

chlorophyll a concentration within a group of 

simulations (fig. 13C), with the exception of 

model versions D and D′; those two versions of 

the model lack an assumption that chlorophyll a

(a surrogate for cyanobacterial biomass) achieves 

equilibrium with the available light and nutrients 

within a time step. 
51

Table 15. Comparison of metrics describing the long-term simulations for model versions discussed in the report.

[Shading indicates four groups of model simulations defined by four combinations of external loading and initial sediment reservoir. Mrec=mass of total phosphorus recycled from sediments to water column; Msed=mass of total phosphorus settled from water column to sediments; EQ=Chlorophyll a concentration evaluated 

by setting dB/dt to 0; RK=Chlorophyll a concentration evaluated using Runge-Kutta approach; INIT=Initial size of sediment reservoir calculated with original 

starting concentration 1,000 mg P/kg sediment; INIT_NEW=Initial size of sediment reservoir recalculated with starting concentration 335 mg P/kg sediment; 

REC_pH=Model relies on pH-dependent recycle mechanism; REC_T=model relies on temperature-dependent recycle mechanism; P_00=Phosphorus limitation 

factor calculated as the percent of non-algal P; P_MM=Phosphorus limitation factor calculated by Michaelis-Menten kinetics; L_AVG=Light limitation factor 

calculated using the average light intensity; L_DI=Average light limitation factor calculated by depth-integrating the light limitation equation.]

Reduction in External Loading Boundary Condition Applied to the Model

0% 40% 0% 40% 0% 40% 0% 40% 0% 40%

Model Characteristics

Ratio of recycling to 

sedimentation, 

Mrec/Msed, average 

over first 7 years of 

203-year simulation

Average concentration over last 7 years of 203-

year simulation

Years with seasonal

 chlorophyll a peak 

>100 ppb over last 7 

years of 203-year 

simulation (out of 7)

Model 

version

Compare 

to 

version

Chlorophyll a

model

Initial 

sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light

limitation Total P (ppb)

Sediment P

(metric tons)

Chlorophyll a

(ppb)

Original model with original calibration

O EQ INIT REC_pH P_00 L_AVG 0.62 0.33 105.9 39.7 2,302 3,338 62.6 25.1 7 2

O' O EQ INIT_NEW REC_pH P_00 L_AVG 0.60 0.52 120.0 71.1 901 932 71.3 44.5 7 6

Recalibration of original model

A O EQ INIT REC_pH P_00 L_AVG 0.57 0.42 102.8 42.6 2,441 3,018 60.9 27.1 7 3 

A' A EQ INIT_NEW REC_pH P_00 L_AVG 0.57 0.53 118.4 69.4 980 897 69.1 44.0 7 6 

Modifications to light and phosphorus limitation

B A EQ INIT REC_pH P_MM L_DI 0.70 0.62 105.9 58.9 2,224 2,067 62.5 32.0 7 5 

B' B EQ INIT_NEW REC_pH P_MM L_DI 0.69 0.43 121.5 73.0 871 1,028 72.0 39.8 7 4 

Temperature-dependent recycling

C B EQ INIT REC_T P_MM L_DI 0.69 0.79 106.3 65.0 2,117 1,299 63.6 39.5 7 7 

C' C EQ INIT_NEW REC_T P_MM L_DI 0.68 0.76 120.1 73.3 867 534 72.6 46.1 7 7 

Runge-Kutta solution of chlorophyll a

D B RK INIT REC_pH P_MM L_DI 1.09 1.18 116.8 70.2 1,569 1,061 54.4 33.9 7 4 

D' D RK INIT_NEW REC_pH P_MM L_DI 1.06 1.13 129.0 78.4 518 355 60.3 38.4 7 6 

Table continues on next page. 
52

Table 15. Comparison of metrics describing the long-term simulations for model versions discussed in the report—continued

[Shading indicates four groups of model simulations defined by four combinations of external loading and initial sediment reservoir.] 

Reduction in External Loading Boundary Condition Applied to the Model

0% 40% 0% 40% 0% 40% 0% 40%

Average concentration over last 7 years of 203-year 

simulation

Years with seasonal

 chlorophyll a peak 

>100 ppb over last 7 

years of 200-year 

simulation (out of 7)

Total P (ppb)

Sediment P

(metric tons)

Chlorophyll a 

(ppb)

Steady state minimum, maximum, and range over all model versions

Minimum, original phosphorus reservoir 102.8 39.7 1,569 1,061 54.4 25.1 7 2

Maximum, original phosphorus reservoir 116.8 70.2 2,441 3,338 63.6 39.5 7 7

Range, original phosphorus reservoir: 14.1 30.5 872 2276 9.3 14.4 - -

Range, original phosphorus reservoir, as percent of average: 12.8% 55.5% 43.5% 103.5% 15.7% 44.6% - -

Minimum, updated phosphorus reservoir: 118.4 69.4 518 355 60.3 38.4 7 4

Maximum, updated phosphorus reservoir: 129.0 78.4 980 1,028 72.6 46.1 7 7 

Range, updated phosphorus reservoir: 10.7 8.9 462 673 12.3 7.7 - -

Range, updated phosphorus reservoir, as percent of average: 8.6% 12.1% 61.7% 97.3% 18.5% 18.3% - -
53

Figure 13. Graphs showing average (A) concentration of total phosphorus, (B) size of the sediment reservoir, (C) 

concentration of chlorophyll a during the final 7 years of a 203-year simulation for Upper Klamath Lake, Oregon, as 

a function of the average ratio of recycled mass to settled mass during the first 7 years of a 203-year simulation. 

(D) Total phosphorus concentration from graph A as a function of the sediment reservoir from graph B. 

All versions of the model produced simulations in which the steady-state phosphorus 

concentration was less than the initial concentration (fig. 14A, B). However, the steady-state 

concentration varied over a range of 14.0 ppb 

(12.8 percent of the average), from 102.8 to 

116.8 ppb when no reduction in external loads 

was applied as a boundary condition (table 15). 

When a reduction in external loads of 40 percent

was applied, the steady-state phosphorus concentration was lower for all model versions than 

when no reduction in external loads was simulat-
54

ed, but varied over a range of 30.5 ppb (55.5 percent of the average), from 39.7 to 70.2 ppb (table 

15). When a 40 percent reduction in external 

loads was simulated, the range in steady-state

concentration was larger and the effect of the difference in settling velocity ualg was apparent in 

the grouping of versions O and A (ualg=0.05 and 

0.059 m/d, respectively) and the grouping of versions B, C and D (ualg=0.035 m/d; therefore, 

lower settling relative to recycling than in versions O and A) (fig. 14B). 

Figure 14. Graphs showing comparison of long-term simulated concentration of water column phosphorus, under 

conditions of (A) 0 and (B) 40 percent reductions in external loads to Upper Klamath Lake, Oregon. Version 

O=original model, original calibration; version A=original model, recalibrated; version B=modifications to light and 

phosphorus limitations to algal growth; version C=version B with replacement of pH-dependent recycling with temperature-dependent recycling; version D=version B without assumption of maximum chlorophyll a growth during 

each time step.

Mechanistic connections between the assumptions that distinguish the model versions and 

the steady-state outcomes are not clear, however. 

For example, versions B, C, and D share the corrections made to phosphorus and light limitations 

to cyanobacterial growth, making those limitation 

factors generally lower (data not shown). That led 

to a smaller value for ualg through the balance 

used to calibrate the algal settling velocity (eq. 

14), even though there was no direct connection 

between how light and phosphorus limitation 

were calculated and the settling velocity. Thus, 

reasonable differences in model assumptions led 

to quite different values of calibration parameters. 

Over the 7.5 years for which observations were 

available, the differences between the calibrated 

model versions were barely perceptible (figs. 8, 

10–12), and they all simulated the 7.5 years of 

available observations with roughly the same 

success (tables 9–14). However, when those same 

models were used to simulate conditions 203 

years into the future, the differences in calibration 

parameters manifested themselves in potentially 

meaningful quantitative differences in steadystate conditions. 

Like water column phosphorus concentration, the size of the sediment reservoir of 

phosphorus achieved equilibrium during the 203-

year simulations (fig. 15). When no reduction in 

external loading was applied as a boundary condition, the steady-state values were always less 

than the initial size of the reservoir, but the range 
55

was large (872 metric tons, or 43.5 percent of the 

average). Water column phosphorus concentration and the size of the sediment reservoir were 

inversely correlated, and the model versions that 

simulated the highest steady-state phosphorus 

concentration also simulated the smallest steadystate sediment phosphorus reservoir (fig. 13D). 

When a 40 percent reduction in external loads 

was applied as a boundary condition, the range in 

steady-state phosphorus reservoir from 1,061 to 

3,338 metric tons (table 15, fig. 15B) encompassed the range in the steady-state phosphorus 

reservoir that resulted when no reduction in external loads was applied (fig. 15A) and at 2,276 

metric tons was more than twice as large (table 

15). Under conditions of reduced external loading, model version O, which had the smallest

Mrec/Msed ratio (table 15) and simulated the 

lowest overall steady-state water column phosphorus concentration (fig. 14B) simulated a 

steady-state phosphorus reservoir that was larger 

than the initial condition (fig. 15B). 

Figure 15. Graphs showing comparison of long-term simulated concentration of sediment phosphorus reservoir in 

Upper Klamath Lake, Oregon, under conditions of (A) 0 and (B) 40 percent reductions in external loads. Version 

O=original model, original calibration; version A=original model, recalibrated; version B=modifications to light and 

phosphorus limitations to algal growth; version C=version B with replacement of pH-dependent recycling with temperature-dependent recycling; version D=version B without assumption of maximum chlorophyll a growth during 

each time step.

When the updated estimate of the initial size 

of the sediment reservoir was used to start the 

simulations, the resulting range in steady-state

phosphorus concentrations was smaller than 

when the original size of the sediment reservoir 

was used (fig. 16A, B)—8.6 percent and 12.1 percent of the average when no reductions and when 

40 percent reductions in external loads were applied as a boundary condition, respectively. The 

steady-state values overall were, however, higher—118.4–129.0 ppb and 69.4–78.4 ppb when no 

reductions and 40 percent reductions in external 

loads were applied, respectively (fig. 13A, table 

15). The updated phosphorus reservoir was 

smaller at the beginning of the simulations but 

also decreased less before achieving a steady 

state (fig. 17A, B). When a 40 percent reduction 

in external loads was applied, model version C′, 

which was characterized by the smallest Mrec/

Msed in this group of models (table 15), produced a small increase in the phosphorus 

reservoir to achieve steady-state conditions. 
56

Figure 16. Graphs showing comparison of long-term simulated concentration of water column phosphorus in Upper 

Klamath Lake, Oregon, under conditions of (A) 0 and (B) 40 percent reductions in external loads. Versions of the 

model marked with a prime symbol (O′, for example) were calibrated using a lower initial condition of the sediment 

phosphorus reservoir than the corresponding versions without the prime symbol. Version O′=original model, original 

calibration; version A′=original model, recalibrated; version B′=modifications to light and phosphorus limitations to 

algal growth; version C′=version B′ with replacement of pH-dependent recycling with temperature-dependent recycling; version D′=version B′ without assumption of maximum chlorophyll a growth during each time step. 

Figure 17. Graphs showing comparison of long-term simulated concentration of sediment phosphorus reservoir in 

Upper Klamath Lake, Oregon, under conditions of (A) 0 and (B) 40 percent reductions in external loads. Versions of 

the model marked with a prime symbol (O′, for example) were calibrated using a lower initial condition of the sediment phosphorus reservoir than the corresponding versions without the prime symbol. Version O′=original model, 

original calibration; version A′=original model, recalibrated; version B′=modifications to light and phosphorus limitations to algal growth; version C′=version B′ with replacement of pH-dependent recycling with temperaturedependent recycling; version D′=version B′ without assumption of maximum chlorophyll a growth during each time 

step. 
57

The number of years in the final 7 years of 

the 203-year simulations in which maximum

chlorophyll a concentration exceeded 100 ppb is 

another metric that can be used to characterize 

the long-term simulation (table 15; fig. 18). With

no reduction in external loads, simulated chlorophyll a exceeded 100 ppb in each of 7 years 

regardless of the model version or the size of the 

initial sediment reservoir (fig. 18A, C). Consistent 

with Walker (2001), when a 40 percent reduction 

in external loads was applied as a boundary condition (fig 18B), model version O showed 

simulated maximum chlorophyll a exceeding 100 

ppb in only 2 of 7 years. Other model versions 

simulated maximum chlorophyll a exceeding 100 

ppb in 3–7 of 7 years (fig. 18B). When the updated estimate of the initial sediment reservoir was 

used, the number of years in which maximum

chlorophyll a exceeded 100 ppb was greater for 

each version of the model than when the larger, 

original initial sediment phosphorus reservoir was 

used (table 15, fig. 18D). Total water column 

phosphorus concentration simulated by each version of the model shows year-to-year variability 

in maximum phosphorus concentration that corresponds directly to the year-to-year variability in 

maximum chlorophyll a (fig. 19).

Model versions C and C′ simulated maximum chlorophyll a exceeding 100 ppb in 7 of 7 

years, demonstrating an important difference between these model versions (in which a 

temperature-dependent recycling mechanism was 

substituted for a pH-dependent one) and the others. Because of the proportionality between pH 

and chlorophyll a, the pH-dependent mechanism 

has built-in positive feedback that the temperature-dependent mechanism does not. In model 

versions O, O′, A, A′, B, B′, D and D′, if cyanobacterial growth does not exceed a threshold, pH 

remains below a critical value for recycling 

(roughly represented by pH*) and recycling is 

constrained, whereas when the threshold is exceeded, recycling is accelerated. This results in 

uneven changes in peak chlorophyll a from year 

to year when the external loading is reduced, with 

very large blooms in some years and very small 

blooms in others, depending on whether the 

threshold is exceeded (fig. 18). Total phosphorus 

concentration shows analogous year-to-year variability (fig. 19). Because the temperaturedependent mechanism has no built-in feedback, 

changes occur more evenly across years, leading 

to maximum chlorophyll a being smaller in all 

years, but not less than 100 ppb in any year. 

The 203-year simulations provided insight 

into the length of time required to achieve steady 

state for each initial condition. Times were longer 

for the original sediment reservoir initial condition than for the updated sediment reservoir 

initial condition: about 50–70 years at the higher 

initial condition (figs. 14A, B and 15A, B) and 

about 25–30 years at the lower initial condition 

(figs. 16A, B and 17A, B). 
58

Figure 18. Graphs showing simulated chlorophyll a concentration during the last 7 years of a 203-year simulation under conditions of (A and C) 0 

percent and (B and D) 40 percent reductions in external loads to Upper Klamath Lake, Oregon. The x-axis is labeled with the date in the 7-year cycle 

of observations that was merged end-to-end to create 203 years of boundary condition data to run the model. Versions of the model marked with a 

prime symbol (O′, for example) were calibrated using a lower initial condition of the sediment phosphorus reservoir than the corresponding versions 

without the prime symbol. Versions O and O′=original model, original calibration; versions A and A′=original model, recalibrated; versions B and 

B′=modifications to light and phosphorus limitations to algal growth; versions C and C′=version B with replacement of pH-dependent recycling with 

temperature-dependent recycling; versions D and D′=version B without assumption of maximum chlorophyll a growth during each time step. 
59

Figure 19. Graphs showing simulated total phosphorus concentration during the last 7 years of a 203-year simulation under conditions of (A and C) 

0 percent and (B and D) 40 percent reductions in external loads to Upper Klamath Lake, Oregon. The x-axis is labeled with the date in the 7-year 

cycle of observations that was merged end-to-end to create 203 years of boundary condition data to run the model. Versions of the model marked 

with a prime symbol (O′, for example) were calibrated using a lower initial condition of the sediment phosphorus reservoir than the corresponding 

versions without the prime symbol. Versions O and O′=original model, original calibration; versions A and A′=original model, recalibrated; versions B

and B′=modifications to light and phosphorus limitations to algal growth; versions C and C′=version B with replacement of pH-dependent recycling 

with temperature-dependent recycling; versions D and D′=version B without assumption of maximum chlorophyll a growth during each time step.
60

Summary, Conclusions, and Suggestions for Improvements to the TMDL 

Model

Verification of Model Algorithms

Verification of the calibration of the original 

TMDL model revealed a few, minor discrepancies 

in calibration parameters; however, these could be 

explained by our interpretation of ambiguous descriptions of the calibration dataset in the model 

documentation. For the most part, the calibration 

of the model was duplicated closely with the original calibration dataset, which included 

phosphorus loading and water-quality data from 

1991 through 1998. Because the semi-analytical 

approach in the original model was difficult to follow and modify, we developed an R-script version 

of the model that used a conventional numerical 

solution technique, and used that code for all the 

simulations reported herein. The R-script code is 

supplied as an appendix to this report. 

The verification process identified three 

changes in model algorithms that would be helpful 

in improving the representation of lake processes

in the model: 

1. The light limitation function would better 

represent the available light for photosynthesis if it were integrated over depth, 

rather than substituting a depth-averaged 

light value into the light-limitation function defined for a specific depth. 

2. The phosphorus limitation function would 

be more defensible if it were modified to 

incorporate a Michaelis-Menten formulation, rather than using a simple percentage 

of nonalgal phosphorus in the water column. 

3. The assumption that chlorophyll a reaches 

equilibrium with light and nutrients over 

the time step does not accurately reflect 

the response time of the algal population to 

these factors, and the equation for chlorophyll a would best be solved numerically, 

as are the equations for water column and 

sediment phosphorus. Implementation of 

change 3 would be facilitated by eliminating the discontinuity in the pH regression 

model at the end of the calendar year, 

which was not done in this study. 

Sensitivity of Calibration Parameters to Model 

Assumptions

We developed a simple temperaturedependent recycling term based on a hypothesis 

that benthic invertebrate activity is an important 

driver of internal phosphorus loading and that 

such activity is stimulated by higher temperature. 

This term used a recycling rate analogous to the 

pH-dependent recycling rate and a threshold temperature that also was treated as a calibration 

parameter. The temperature-dependent recycling 

term, therefore, required two calibration parameters, as did the pH-dependent term. As a result, 

when this temperature-dependent recycling mechanism was substituted for the pH-dependent 

mechanism in the original model, the total number 

of calibration parameters did not change. 

We incorporated each of these three modifications to the original model—changes to the light 

and phosphorus limitation to growth, substitution 

of a temperature-dependent recycling term for the 

pH-dependent term, and elimination of the assumption of maximum cyanobacterial growth 

during each time step—into new versions of the 

model. The algal settling velocity and recycling 

rate, or the algal settling velocity alone, as appropriate, were recalibrated using the same 

calibration dataset that was used in the calibration 

of the original model. The three new versions of 

the model, along with the original model and original calibration, as well as a recalibration of the 

original model done as part of this study, provided 

five calibrated versions of the model that could be 

compared by their performance over the validation 

dataset. 

The calibrated settling velocity supplied with 

the original model was 0.05 m/d. The value of this 

parameter in the four versions of the model creat-
61

ed as part of this study varied from 0.059 m/d 

(+18 percent) to 0.035 m/d (-30 percent). The calibrated recycling rate supplied with the original 

model was 0.53/yr. The value of this parameter in 

three versions of the model created as part of this 

study that included the pH-dependent recycling 

mechanism varied from 0.32/yr (-39.6 percent) to 

0.47/yr (-11.3 percent). The value of the recycling 

rate in the version of the model that incorporated 

temperature-dependent recycling instead of pHdependent recycling was 0.54/yr (+1.9 percent). 

These represent substantial changes to calibration 

parameters that are not easily related mechanistically to the modifications that define the different 

versions of the model. For example, the 30 percent

decrease in the settling velocity and 40 percent

decrease in the recycling rate occurred when the 

light and phosphorus limitation factors were modified as indicated. Rather, these changes in 

calibration parameters are the consequence of repeating the calibration process after making 

seemingly unrelated changes to the model assumptions. 

A final sensitivity analysis examined the consequences of changing the initialization of the 

sediment phosphorus reservoir. The most recent 

estimates, based on a set of 26 cores that provided 

a more complete coverage of the lake than was 

available when the TMDL was adopted, indicate 

that the previous initialization of the sediment 

phosphorus reservoir was too high by as much as 

66 percent. The recycling rate in each of the five 

model versions was recalibrated using this new 

estimate. The calibrated recycling rate in each version of the model increased to compensate for the 

lower sediment concentration. (Because the recycled mass is proportional to both the size of the 

reservoir and the recycling rate, to achieve approximately the same balance in equation 1, the 

rate must increase if the phosphorus content decreases.) The pH-dependent recycling rate in this 

second set of five versions of the model varied 

from 1.05 to 1.33/yr, and the temperaturedependent recycling rate was 1.62/ yr. 

Although the changes in calibration parameters among all 10 versions of the model were 

substantial, the performance of the model based 

on statistics calculated over the full verification 

data set was not much changed, and was uniformly poor. Further, no model version consistently 

produced the best or worse results among the three 

state variables—total phosphorus, chlorophyll a, 

and pH. When only May–July was considered in 

calculating statistics, performance was better. 

Generally, the differences among the various 

model versions were small, but the version of the 

model that did not assume maximum growth during the time step performed better overall in 

simulating chlorophyll a and pH. All versions of 

the model, therefore, performed better during 

months of rapid bloom expansion than during the 

other months of the year. 

Effect of Calibration Parameters on Long-Term 

Model Simulations

Although the differences in model performance statistics among the various versions tested 

were small, when the models were used to examine the long-term steady-state conditions, the 

differences in the value of calibration parameters 

manifested themselves as sometimes large differences in concentration over time frames exceeding 

several decades. All versions of the model, when 

used to simulate water column concentrations 

characterizing a steady state 203 years into the 

future, produced lower total phosphorus concentrations over time, even when no reductions in the 

external load boundary condition were applied. 

When 40 percent reductions in the external load 

boundary condition were applied, the steady-state

concentrations were always lower than when no 

reduction was applied. Within that broadly consistent result, however, there were large variations 

in the steady-state concentration that could be related to the ratio of the recycling term to the 

sedimentation term in the mass balance, with more 

recycling relative to sedimentation resulting in 

higher concentrations. The values of the recycling 

and sedimentation terms are directly dependent on 

the value of the calibration parameters—in particular the algal settling velocity and the recycling 

rate.
62

The phosphorus stored in the sediment reservoir generally—although not always—decreased 

over time as well, and as with the water column 

concentrations, there was a wide range in the 

steady-state value that could be related to the ratio 

of the recycling term to the sedimentation term in 

the mass balance and the values of the calibration 

parameters that largely determine that ratio. The 

higher the ratio, indicating that there is more recycling relative to sedimentation, the lower the 

steady-state phosphorus reservoir. It follows that 

the steady-state water column phosphorus concentration was negatively correlated with the steadystate sediment reservoir—a lower steady-state

sediment reservoir was associated with higher 

steady-state water column phosphorus concentration. 

This last point became apparent when versions of the model that were calibrated to an 

updated, lower initial sediment reservoir were 

used to produce 203-year simulations. Models calibrated to the lower initial sediment reservoir 

simulated higher steady-state water column phosphorus concentrations than the counterpart 

calibrated to the higher sediment reservoir. For 

example, for models calibrated to a higher sediment reservoir, the steady-state phosphorus 

concentration ranged from about 103 to 117 ppb 

when no reduction in the external load boundary 

condition was applied, and from about 40 to 70 

ppb when a 40 percent reduction in the external 

load boundary condition was applied. For models 

calibrated to the lower sediment reservoir, the 

steady-state phosphorus concentration ranged 

from about 118 to 129 ppb when no reduction in 

the external load boundary condition was applied 

and from about 69 to 78 ppb when a 40 percent 

reduction in the external load boundary condition 

was applied. Therefore, the sediment reservoir influences the long-term behavior of the modeled 

water column concentrations in an important and 

possibly counterintuitive way. These model results 

indicate that sediment phosphorus does not act as 

a supply that moves uniformly toward depletion 

over time. Rather, the model algorithms require 

the lake water column to reach a steady state in 

equilibrium with the sediments, and the steadystate conditions in the water column and sediments are determined by the calibration 

parameters. In some simulations, the size of the 

sediment phosphorus reservoir stayed mostly 

steady or increased slightly over time as that equilibrium state was approached. A corollary is that 

the initial size of the sediment reservoir is not as 

important in determining steady-state conditions 

as the rate of the lake processes that add phosphorus to the sediment through settling and recycle 

phosphorus from sediments to the water column. 

Another metric that described steady-state

conditions was the number of years in the final 7-

year cycle in which peak chlorophyll a exceeded 

100 ppb. This number was 7 out of 7 for all 10 

model versions when no reduction in external load 

boundary conditions was applied—in other words, 

without a reduction in external loads, even though 

phosphorus and chlorophyll a concentrations decreased over time, the peak chlorophyll a

concentration still exceeded 100 ppb every year. 

When a 40 percent reduction in external load was 

applied, however, the number of years in which 

peak chlorophyll a exceeded 100 ppb decreased to 

between 2 and 6 years out of 7, depending on the 

model version, and for only those versions that 

incorporated a pH-dependent recycle mechanism. 

Because pH is proportional to the logarithm 

of chlorophyll a, and recycling is in turn proportional to pH when a pH-dependent recycle 

mechanism is used, reductions in chlorophyll a are 

manifested as reductions in recycling, which lead 

to further reductions in chlorophyll a, and so on. 

This leads to a bimodal distribution of chlorophyll 

a, corresponding to very large blooms in some 

years and very small blooms in other years. Any 

recycling mechanism that does not have the same 

feedback would lead to a unimodal distribution. 

This was demonstrated in this study with a simple

temperature-dependent recycling mechanism. 

Versions of the model in which a temperaturedependent recycle mechanism was substituted for 

a pH-dependent mechanism simulated chlorophyll 

a concentration maxima at steady state that were 

lower during large bloom years, but still exceeded 
63

100 ppb in every year because recycling was more 

even across all years. Because new evidence indicates that pH-dependent desorption of phosphorus 

from iron oxyhydroxides is one of several mechanisms for the internal loading of phosphorus in the 

lake, and given that the form of the recycling term 

has a strong influence on the model outcome, efforts to improve the predictive power of the model 

likely will be most successful if they include a refinement of this recycling term. 

Suggestions for Improvements to the TMDL 

Model

The goal of the sensitivity analysis was not to 

create a more “correct” model or to determine the 

“correct” calibration parameters, but rather to 

demonstrate sensitivity of the calibration parameters to model assumptions and initial conditions. If 

the comparison is limited to model versions that 

were calibrated with the updated initial phosphorus reservoir, calibrated values of the recycle rate 

over all the model modifications tested ranged

from 1.05 to 1.62 per year. The settling velocity of 

algal phosphorus ranged from 0.035 to 0.059 m/d. 

The nonalgal phosphorus loss rate was not revisited as part of the sensitivity analysis because it was 

calculated based on observations of wintertime 

total phosphorus, but the linear model used to determine the nonalgal phosphorus loss rate did not 

have predictive power and, therefore, confidence 

in the parameter value was low. Such a wide range 

in calibrated values indicates that the parameter 

calibration is not converging on a value that represents a true physical or chemical process. Rather, 

it indicates that the parameter calibration is compensating for imprecise formulations of those 

processes and that some of the streamlined equations in the model are not adequate placeholders 

for the more complex phenomena that occur in the 

lake. 

The model could be improved by reducing 

uncertainty in the recycling and sedimentation 

terms. One approach would be to undertake a field 

program to directly measure or estimate by other 

means the parameters required for the recycling 

and sedimentation terms in the model, and in the 

process to either determine appropriate constant 

values (as are used in the original model for algal 

settling rate and nonalgal loss rate) or develop 

new algorithms that can be implemented within 

the model to predict these terms as functions of 

other environmental variables that are available as 

input data. In the case of macro-invertebratedependent recycling, for example, two important 

explanatory variables might be water temperature 

and time elapsed since the start of the growing 

season, which would be related to the timing of 

increased invertebrate growth and activity in the 

spring and the progression of species cohorts. 

The contribution of benthic invertebrates to 

nutrient cycling in Upper Klamath Lake has begun 

to be quantified; however, much remains to be 

learned regarding their role. For example, leeches 

(functionally biodiffusers) represent a significant 

fraction of the benthic assemblage at all spatial 

and temporal scales thus far examined, but we 

have yet to consider their roles in phosphorus release. Additionally, the rate at which nutrients are 

released by a species possibly changes at different 

life stages. Given the synchronicity of many invertebrate life stages, this effect, if present, could be 

important relative to nutrient recycling both spatially and temporally. Further, trophic interactions 

likely affect A. flos-aquae, Microcystis, tubificid

worm, leech, midge, and fish population dynamics 

in ways that influence phosphorus cycling in Upper Klamath Lake. Additionally, the influence of 

zooplankton on nutrient cycling has yet to be studied in the lake. Better defining the interactions 

between climate, hydrodynamics, chemistry, and 

biology would greatly enhance our understanding 

of the ecology of Upper Klamath Lake. It is not 

within the scope of this study to provide a detailed 

plan for making such measurements/estimates; 

however, some tools and techniques that likely 

would be useful can be identified. 

Quantitative assessment of the contribution of 

metabolic cycling and excretion by benthic invertebrates to the recycling of phosphorus from the 

sediments through a season would require collection of density, body-size, and body-stoichiometry

data for the most dominant taxa at representative 
64

sites in the lake on at least a biweekly time scale

(the time step of the model) and over a few growing seasons to capture intra-annual variability.

Inverse modeling techniques combined with pore 

water profiles can be used to determine bioturbation and bio-irrigation parameters (Meile and 

others, 2001) that define the enhancement of diffusion processes by bioturbation and bioirrigation. Such techniques could be used in combination with silica profiles (silica acts like a 

conservative tracer) that already have been collected in the lake and wetlands with a pore water 

profiler (Kuwabara and others, 2007). A calibrated 

early diagenesis sediment model in combination 

with pore water profiles of reacting solutes, such 

as redox sensitive species and pH (Berg and others, 2003), could be used to infer degradation rates 

of settled organic material. Because recycling can 

be expected to vary spatially around the lake, as 

does sediment phosphorus, benthic invertebrate 

density, and the sedimentation of fresh organic 

material, a spatially explicit phosphorus transport 

model of the lake could be used to appropriately 

“scale up” the spatial variability from sitedependent estimates for use in a CSTR model. 

A fundamental problem regarding sedimentation is that the model uses sedimentation to 

account for the process of removal of algal phosphorus from the water column, but no mortality 

term in the equations moves algal phosphorus 

from a viable (largely non-settling) to a detrital 

(settling) physiological state. Uniform settling of 

the large buoyant colonies that dominate the cyanobacterial population in Upper Klamath Lake 

through the season is unlikely. Rather, more settling likely occurs during bloom declines. 

Experimentation with sediment traps in Upper 

Klamath Lake (Eldridge and others, 2013) determined that settling algal material was partially 

degraded, not bright green and healthy, which 

supports the idea that settling increases when the 

bloom is declining. Large, net negative internal 

loading based on mass balance calculations also 

indicates that large quantities of suspended material settle during periods of algal decline (Kann 

and Walker, 2001). This is the reason that the performance of the model degrades when it is validated with data from outside May–July, when the 

bloom is usually growing rapidly. A rapid decline 

in the bloom that the model cannot simulate occurs at some point after the second week in July in 

most years, roughly coincident with the seasonal 

maximum temperature (Lindenberg and others, 

2009). Consequently, the model consistently 

overpredicts total phosphorus (by overpredicting 

the algal phosphorus component) and particularly 

chlorophyll a concentrations after the bloom decline, as noted by Walker (2001). 

To address this problem, the model could be 

improved by incorporating a mortality term. As 

noted by Walker (2001), the reason for the seasonal rapid bloom decline was not known when 

the model was developed, so no term was developed to describe the phenomenon. Although the 

cause still remains largely unknown, 20 years of 

biweekly profiles collected by the Klamath Tribes 

and 8 years of continuous monitor data collected 

by the USGS provide a description of the seasonal 

bloom dynamics over a wide range of conditions. 

An adequate empirical description of the midseason bloom decline possibly could be developed 

for incorporation into the TMDL model with the 

use of those data. 

The previously noted difficulty with calibration of the nonalgal phosphorus loss rate is caused 

by the inability of the model to distinguish different forms of nonalgal phosphorus that have 

substantially different loss rates. For example, the 

principal components analysis of bed sediment 

cores provided evidence that particle-associated 

inorganic phosphorus settles out of the water column downstream of the Williamson River mouth. 

The model could be improved by tracking particle-associated and dissolved nonalgal phosphorus 

separately. This would require information about 

the distribution of phosphorus between the two 

fractions at the tributary boundaries. 

Finally, the TMDL model parameters, and 

consequently the character of the simulations (particularly the long-term 203-year response), are 

likely to be dependent on the dataset used during 
65

calibration. Several simulations of the long-term 

behavior of the model have been presented in this 

report, and these simulations generally indicate 

that a steady state is achieved within about 50 

years. All the calibrated models were consistent in 

that they resulted in a final steady-state in which 

total phosphorus concentration in the lakes was

reduced, and most models simulated severe 

blooms (>100 ppb) happening less often by the 

time steady state was reached. However, the calibration process essentially fits the early part of the 

simulation to the available data, but the long-term 

response is an extrapolation that depends on the 

form of the model and its fitting parameters. In 

this case, the water-quality data from 1991 

through 1998 that were used to calibrate the model 

showed decreasing trends in the average inflow 

and lake concentrations by about 2.8 percent per 

year and 8.2 percent per year, respectively (Walker, 2001). Those decreasing trends are 

incorporated into the fitting parameters. An updated mass balance dataset was not available for 

updating the model, but when a longer term dataset becomes available, the potential changes to 

calibrated parameters based on a longer calibration dataset and the performance of the model 

over a longer validation dataset can be evaluated.

References Cited

Aldous, A.R., Craft, C.B., Stevens, C.J., Barry, 

M.J., and Bach, L.B., 2007, Soil phosphorus release from a restoration wetland, Upper 

Klamath Lake, Oregon: Wetlands, v. 27, no. 4, 

p. 1025–1035. 

Aldous, A., McCormich, P., Ferguson, C., Graham, S., and Craft, C., 2005, Hydrologic regime 

controls soil phosphorus fluxes in restoration 

and undisturbed wetlands: Restoration Ecology, 

v. 13, no. 2, p. 341–347. 

Alves, J.M., Caliman, A., Guariento, R.D., 

Figueiredo-Barros, M.P., Carneiro, L.S., Farjalla, V.F., Bozelli, R.L., and Esteves, F.A., 

2010, Stoichiometry of benthic invertebrate nutrient recycling—Interspecific variation and the 

role of body mass: Aquatic Ecology, v. 44, p. 

421–430. 

Andersson, G., Graneli, W., and Stenson, J., 1988, 

The influence of animals on phosphorus cycling 

in lake ecosystems: Hydrobiologia, v. 170, p. 

267–284. 

Banish, N.P., Adams, B.J., and Shively, R.S., 

2009, Distribution and habitat associations of 

radio-tagged adult Lost River suckers and 

shortnose suckers in Upper Klamath Lake, Oregon: Transactions of the American Fisheries 

Society, v. 138, p. 153–168. 

Barbiero, R.P., and Kann, J., 1994, The importance of benthic recruitment to the 

population development of Aphanizomenon 

flos-aquae and internal loading in a shallow 

lake: Journal of Plankton Research, v. 16, p. 

1581–1588. 

Berg, P., Rysgaard, S., and Thamdrup, B., 2003, 

Dynamic modeling of early diagenesis and nutrient cycling—A case study in an arctic marine 

sediment: American Journal of Science, v. 303, 

p. 905–955. 

Bowie, G.L., Mills, W.B., Porcella, D.B., Campbell, C.L., Pagenkopf, J.R., Rupp, G.L., 

Johnson, K.M., Chan, P.W.H., Gherini, S.A., 

and Chamberlin, C.E., 1985, Rates, constants, 

and kinetics formulations in surface water quality modeling: U.S. Environmental Protection 

Agency, EPA/600/3-85/040, 455 p. 

Bradbury, J.P., Colman, S.M., and Reynolds, R.L., 

2004, The history of recent limnological changes and human impact on Upper Klamath Lake, 

Oregon: Journal of Paleolimnology, v. 31, p. 

151–165. 

Burden, R.L., and Faires, J.D., 1997, Numerical 

analysis: Pacific Grove, California, Brooks/Cole 

Publishing Co., 811 p. 

Devine, J.A., and Vanni, M.J., 2002, Spatial and 

seasonal variation in nutrient excretion by benthic invertebrates in a eutrophic reservoir: 

Freshwater Biology, v. 47, p. 1107–1121. 

Eilers, J.M., Kann, J., Cornett, J., Moser, K., St. 

Amand, A., and Gubala, C., 2001, Recent paleo-
66

limnology of Upper Klamath Lake, Oregon: 

Roseburg, Oreg., JC Headwaters, Inc., 44 p. 

Eilers, J.M., Kann, J., Cornett, J., Moser, K., and 

St. Amand, A., 2004, Paleolimnological evidence of change in a shallow, hypereutrophic 

lake—Upper Klamath Lake, Oregon, USA: Hydrobiologia, v. 520, p. 7–18. 

Eldridge, S.L., Wood, T.M., Echols, K.R., and 

Topping, B.R., 2013, Microcystins, nutrient dynamics, and other environmental factors during 

blooms of non-microcystin-producing Aphanizomenon flos-aquae in Upper Klamath Lake, 

Oregon, 2009: Lake and Reservoir Management, v. 29, p. 68-81. DOI: 

10.1080/10402381.2013.775199. 

Eldridge, D.B., Caldwell Eldridge, S.L., Schenk, 

L.S., and Wood, T.M., 2012, Water quality data 

from Upper Klamath and Agency Lakes, Oregon, 2009–10: U.S. Geological Survey OpenFile Report 2012-1142, 32 p. (Also available at 

http://pubs.usgs.gov/of/2012/1142.) 

Elsroad, A., Rudd, N., and Hendrixson, H., 2011, 

Williamson River Delta Preserve vegetation 

monitoring—Tulana third-year post-breaching results: Klamath Falls, Oregon, The Nature 

Conservancy, 14 p., accessed August 26, 2013, at 

http://conserveonline.org/library/williamsonriver-delta-preserve-vegetation-2/@@view.html.

Fisher, L.H., and Wood, T.M., 2004, Effect of water-column pH on sediment-phosphorus release 

rates in Upper Klamath Lake, Oregon, 2001: 

U.S. Geological Survey Water-Resources Investigations Report 03-4271, 25 p. 

Forsberg, B.R., and Shapiro, J., 1981, The effects 

of artificial destratification on algal populations, 

in Stefan, H.G., ed., Proceedings of the Symposium on Surface Water Impoundments, June 2–

5, 1980, Minneapolis, Minnesota: New York, 

American Society of Civil Engineers, p. 851–

864. 

Fukuhara, H., and Sakamoto, M., 1987, Enhancement of inorganic nitrogen and phosphate 

release from lake sediment by tubificid worms 

and chironomid larvae: Oikos, v. 38, p. 312–

320. 

Fukuhara, H., and Yasuda, K., 1985, Phosphorus 

excretion by some zoobenthos in a eutrophic 

freshwater lake and its temperature dependency:

Japanese Journal of Limnology, v. 46, p. 287–

296. 

Gardner, W.S., Nalepa, T.F., Quigley, M.A., and 

Malczyk, J.M., 1981, Release of phosphorus by 

certain benthic invertebrates: Canadian Journal 

of Fisheries and Aquatic Science, v. 38, p. 978–

981. 

Graham, S.A., Craft, C.B., McCormick, P.V., and 

Aldous, A., 2005, Forms and accumulation of 

soil P in natural and recently restored peatlands—Upper Klamath Lake, Oregon, USA: 

Wetlands, v. 25, no. 3, p. 594–606. 

Hansen, K., Mouridsen, S., and Kristensen, E., 

1998, The impact of Chironomus plumosus larvae on organic matter decay and nutrient (N,P) 

exchange in a shallow eutrophic lake sediment 

following a phytoplankton sedimentation: Hydrobiologia, v. 364, p. 65–74. 

Haykin, S., 1999, Neural Networks: A Comprehensive Foundation 2nd Edition, Prentice Hall, 

Upper Saddle River, New Jersey. 

Hoilman, G.R., Lindenberg, M.K., and Wood, 

T.M., 2008, Water quality conditions in Upper 

Klamath and Agency Lakes, Oregon, 2005: U.S. 

Geological Survey Scientific Investigations Report 2008-5026, 44 p. 

Jacoby, J.M., Lynch, D.D., Welch, E.B., and Perkins, M.A., 1982, Internal phosphorus loading 

in a shallow eutrophic lake: Water Research, v. 

16, p. 911–919. 

Kaiser, H.F., 1958, The varimax criterion for analytic rotation in factor analysis: Psychometrika, 

vol. 23, no. 3, p. 187-200. 

Kaiser, H.F., 1960, The application of electronic 

computers to factor analysis: Educational and 

Psychological Measurement, vol. 20, 141-151. 

Kann, J., 1998, Ecology and water quality dynamics of a shallow hypereutrophic lake dominated 

by cyanobacteria (Aphanizomenon flos-aquae): 

Chapel Hill, University of North Carolina, PhD 

dissertation, 110 p. 
67

Kann, J., 2011a, Upper Klamath Lake 2010 data 

summary report: Ashland, Oregon, Aquatic 

Ecosystem Sciences LLC, 54 p.

Kann, J., 2011b, Upper Klamath Lake tributary 

loading—2010 data summary report: Ashland, 

Oregon, Aquatic Ecosystem Sciences, LLC, 31 

p. 

Kann, J., and Smith, V.H., 1999, Estimating the 

probability of exceeding elevated pH values 

critical to fish populations in a hypereutrophic 

lake: Canadian Journal of Fisheries and Aquatic 

Sciences, v. 56, p. 2,262–2,270. 

Kann, Jacob, and Walker, W.W., 2001, Nutrient 

and hydrologic loading to Upper Klamath Lake:

Klamath Tribes Natural Resource Department,

U.S. Bureau of Reclamation Cooperative Studies, 48 p. 

Kann, Jacob, and Welch, J.B., 2005, Wind control 

on water quality in shallow, hypereutrophic Upper Klamath Lake, Oregon: Lake and Reservoir 

Management, v. 21, no. 2, p. 149–158. 

Koski-Vahala, J., and Hartikainen, H., 2001, Assessment of the risk of phosphorus loading due 

to resuspended sediment: Journal of Environmental Quality, v. 30, p. 960–966. 

Kuwabara, J.S., Lynch, D.D., Topping, B.R., 

Murphy, F., Carter, J.L., Simon, N.S., Parchaso, 

F., Wood, T.M., Lindenberg, M.K., Wiese, K., 

and Avanzino, R.J., 2007, Quantifying the benthic source of nutrients to the water column of 

Upper Klamath Lake, Oregon: U.S. Geological 

Survey Open-File Report 2007–1276, 39 p. (Also available at 

http://pubs.usgs.gov/of/2007/1276/) 

Kuwabara, J.S., Topping, B.R., Carter, J.L., 

Wood, T.M., Parchaso, F., Cameron, J.M., Asbill, J.R., Carlson, R.A., Fend, S.V., 2012a, 

Time scales of change in chemical and biological parameters after engineered levee breaches 

adjacent to Upper Klamath and Agency Lakes, 

Oregon: U.S. Geological Survey Open-File Report 2012–1057, 26 p. (Also available at 

http://pubs.usgs.gov/of/2012/1057/) 

Kuwabara, J.S., Topping, B.R., Carter, J.L., 

Wood, T.M., Cameron, J.M., Asbill-Case, J.R., 

Carlson, R.A., 2012b, Changes in benthic nutrient sources within a wetland after hydrologic 

reconnection: Environmental Toxicology and 

Chemistry, v. 31, no. 9, p. 1995–2013. 

Kuwabara, J.S., Topping, B.R., Lynch, D.D., 

Carter, J.L., and Essaid, H.I., 2009, Benthic nutrient sources to hypereutrophic Upper Klamath 

Lake, Oregon, USA: Environmental Toxicology 

and Chemistry, v. 28, p. 516–524. 

Laenen, A., and LeTourneau, A.P., 1996, Upper 

Klamath basin nutrient-loading study—Estimate 

of wind-induced resuspension of bed sediment 

during periods of low lake elevation: U.S. Geological Survey Open File Report 95-414, 11 p. 

(Also available at 

http://pubs.er.usgs.gov/publication/ofr95414.) 

Lindenberg, M.K., Hoilman, Gene, and Wood, 

T.M., 2009, Water quality conditions in Upper 

Klamath and Agency Lakes, Oregon, 2006: U.S. 

Geological Survey Scientific Investigations Report 2008-5201, 54 p. (Also available at 

http://pubs.usgs.gov/sir/2008/5201/.) 

Meile, C., Koretsky, C.M., and Van Cappellen, P., 

2001, Quantifying bioirrigation in aquatic sediments—An inverse modeling approach: 

Limnology and Oceanography, v. 46, p. 164–

177. 

Mermillod-Blondin, F., 2011, The functional significance of bioturbation and biodeposition on 

biogeochemical processes at the water-sediment 

interface in freshwater and marine ecosystems: 

Journal of the North American Benthological 

Society, v. 30, no. 3, p. 770–778. 

Mermillod-Blondin, F., and Rosenberg, R., 2006, 

Ecosystem engineering: the impact of bioturbation on biogeochemical processes in marine and 

freshwater benthic habitats: Aquatic Science, v. 

68, p. 434–442. 

Michaud, E., Desrosiers, G., Mermillod-Blondin, 

F., Sundby, B., and Stora, G., 2005, The functional group approach to bioturbation—The 

effects of biodiffusers and gallery-diffusers of 

the Macoma balthica community on sediment 

oxygen uptake: Journal of Experimental Marine 

Biology and Ecology, v. 326, p. 77–88. 
68

Michaud, E., Desrosiers, G., Mermillod-Blondin, 

F., Sundby, B., and Stora, G., 2006, The functional group approach to bioturbation—II. The 

effects of the Macoma balthica community on 

fluxes of nutrients and dissolved organic carbon 

across the sediment-water interface: Journal of 

Experimental Marine Biology and Ecology, v. 

337, p. 178–189. 

Morace, J.L., 2007, Relation between selected water-quality variables, climatic factors, and lake 

levels in Upper Klamath and Agency Lakes, Oregon, 1990–2006: U.S. Geological Survey 

Scientific Investigations Report 2007–5117, 54 p. 

(Also available at 

http://pubs.usgs.gov/sir/2007/5117/.) 

Nash, J.E. and Sutcliffe, J.V., 1970, River flow 

forecasting through conceptual models part I —

A discussion of principles: Journal of Hydrology, v. 10, no. 3, p. 282–290.

Niemisto, J., Holmroos, H., and Horppila, J., 2011, 

Water pH and sediment resuspension regulating 

internal phosphorus loading in a shallow lake—

Field experiment on diurnal variation: Journal of 

Limnology, v. 70, no. 1, p. 3–10.

Nogaro, G., Mermillod-Blondin, F., Valett, M.H., 

Francois-Carcaillet, F., Gaudet, J.-P., Lafont, 

M., and Gibert, J., 2009, Ecosystem engineering 

at the sediment-water interface—Bioturbation 

and consumer-substrate interaction: Oecologia, 

v. 161, p. 125–138. 

Oregon Department of Environmental Quality, 

2002, Upper Klamath Lake drainage Total Maximum Daily Load (TMDL) and Water Quality

Management Plan (WQMP): Oregon Department 

of Environmental Quality, 188 p., accessed August 

8, 2013, at 

http://www.deq.state.or.us/WQ/TMDLs/klamath.

htm. 

Oregon Department of Environmental Quality, 

2010, Upper Klamath and Lost River subbasins 

total maximum daily load (TMDL) and water

quality management plan (WQMP): Oregon Department of Environmental Quality, 231 p., plus 

app., accessed August 8, 2013, at 

http://www.deq.state.or.us/WQ/TMDLs/klamath.

htm.) 

Perkins, D.L., Kann, J., and Scoppettone, G.G., 

2000, The role of poor water quality and fish 

kills in the decline of endangered Lost River and 

shortnose suckers in Upper Klamath Lake, U. S. 

Geological Survey, Biological Resources Division Report submitted to the U. S. Bureau of 

Reclamation, Contract 4-AA-29-12160: Klamath Falls, Oregon, U.S. Geological Survey, 

Klamath Falls Project Office.

Quintana, C.O., Tang, Min, and Kristensen, Erik, 

2007, Simultaneous study of particle reworking, 

irrigation transport and reaction rates in sediment bioturbated by the polychaetes 

Heteromastus and Marenzelleria: Journal of 

Experimental Marine Biology and Ecology, v. 

352, p. 392–406. 

R Core Team, 2012, R–A language and environment for statistical computing: Vienna, Austria, 

R Foundation for Statistical Computing, ISBN 

3-900051-07-0, accessed August 20, 2013, at 

http://www.R-project.org/. 

Redfield, A.C., 1934, On the proportions of organic derivations in sea water and their relation to 

the composition of plankton, in Daniel, R.J., ed., 

James Johnstone Memorial Volume: University 

Press of Liverpool, p. 177–192. 

Rounds, S.A., and Sullivan, A.B., 2009, Review of

Klamath River total maximum daily load models 

from Link River Dam to Keno Dam, Oregon: 

U.S. Geological Survey Administrative Report, 

37 p., accessed August 20, 2013, at 

http://www.usbr.gov/mp/kbao/waterquality/klamath

_river_model_review_final.pdf. 

Rounds, S.A., and Sullivan, A.B., 2013, Review of

revised Klamath River Total Maximum Daily

Load models from Link River Dam to Keno

Dam, Oregon: U.S. Geological Survey OpenFile Report 2013-1136, 31 p. (Also available at 

http://pubs.usgs.gov/of/2013/1136/.) 

Senft, W.H., 1978, Dependence of light-saturated 

rates of algal photosynthesis on intra-cellular concentrations of phosphorus: Limnology and 

Oceanography, v. 23, p. 709–718. 
69

Sereda, J.F., and Hudson, J.J., 2011, Empirical models for predicting the excretion of nutrients (N and 

P) by aquatic metazoans—Taxonomic differences 

in rates and element ratios: Freshwater Biology, v. 

56, p. 250–263.

Simon, N.S., Lynch, D., and Gallaher, T.N., 2009, 

Phosphorus fractionation in sediment cores collected in 2005 before and after onset of an 

Aphanizomenon flos-aquae bloom in Upper 

Klamath Lake, OR, USA: Water, Air, and Soil 

Pollution, v. 204, p. 139–153, accessed August 

20, 2013, at 

http://springerlink.com/content/06580113p6q415v

7/fulltext.pdf.

Simon, N.S., and Ingle, S.N., 2011, Physical and 

chemical characteristics including total and geochemical forms of phosphorus in sediment from 

the top 30 centimeters of cores collected in October 2006 at 26 sites in Upper Klamath Lake, 

Oregon: U.S. Geological Survey Open File Report 2011-1168, 49 p. (Available at 

http://pubs.usgs.gov/of/2011/1168/.) 

Sondergaard, M., 1988, Seasonal variations in the 

loosely sorbed phosphorus fraction of the sediment of a shallow and hypereutrophic lake: 

Environmental Geology and Water Sciences, v. 

11, p. 115–121. 

Sondergaard, M., Kristensen, P., and Jeppesen, E., 

1992, Phosphorus release from resuspended sediment in the shallow and wind exposed Lake 

Arreso, Denmark: Hydrobiologia, v. 228, p. 91–

99. 

Talling, J.F., 1957, Photosynthetic characteristics 

of some freshwater plankton diatoms in relation 

to underwater radiation: New Phytologist, v. 56, 

no. 1, p. 29–50. 

Vannote, R.L., and Sweeney, B.W., 1980, Geographic analysis of thermal equilibria—A 

conceptual model for evaluating the effect of 

natural and modified thermal regimes on aquatic 

insect communities: The American Naturalist, v. 

115, no. 5, p. 667–693. 

Vollenweider, R.A., 1975, Input-output models 

with special reference to the phosphorus loading 

concept in limnology: Schweizererische 

Zeitschrift für Hydrologie, v. 37, p. 53–84. 

Walker, W.W., Jr., 2001, Development of a phosphorus TMDL for Upper Klamath Lake, 

Oregon: prepared for Oregon Department of 

Environmental Quality, Bend, Oregon, 80 p., 

accessed August 26, 2013, at 

http://www.deq.state.or.us/wq/tmdls/docs/klama

thbasin/ukldrainage/devphostmdl.pdf. 

Wetzel, R.G., 1975, Limnology: W.B. Saunders, 

743 p. 

Wildung, R.E., and Schmidt, R.L., 1973, Phosphorus release from lake sediments, prepared for 

Office of Research and Monitoring, U.S. Environmental Protection Agency, Washington, D.C.: 

U.S. Environmental Protection Agency report 

EPA-R3-73-024, 185 p. 

Wildung, R.E., Schmidt, R.L., and Routson, R.C., 

1977, The phosphorus status of eutrophic lake 

sediments as related to changes in limnological 

conditions—Phosphorus mineral components: 

Journal of Environmental Quality, v. 6, no. 1, p. 

100–104, accessed August 9, 2013, at 

http://dx.doi.org/10.2134/jeq1977.0047242500060

0010022x. 

Wong, S.W., Barry, M.J., Aldous, A.R., Rudd, N.T., 

Hendrixson, H.A., and Doehring, C.M., 2011, Nutrient release from a recently flooded delta 

wetland—Comparison of field measurements to 

laboratory results: Wetlands, v. 31, p. 433–443. 

Wong, S., and Hendrixson, H., 2011, Water quality 

conditions in the Williamson River delta, Oregon—Three years post-restoration: Klamath Falls, 

Oregon, The Nature Conservancy Klamath Basin 

Field Office, 2010 Annual Data Report, 38 p. 

Wood, T.M., Fuhrer, G.J., and Morace, J.L., 1996, 

Relation between selected water-quality variables 

and lake level in Upper Klamath and Agency 

Lakes, Oregon: U.S. Geological Survey WaterResources Investigations Report 96-4079, 57 p. 

(Also available at 

http://pubs.er.usgs.gov/publication/wri964079) 

Wood, T.M., 2001, Sediment oxygen demand in 

Upper Klamath and Agency Lakes, Oregon, 1999: 

U.S. Geological Survey Water-Resources Investigations Report 01-4080, 13 p. 
70

Wood, T.M., Hoilman, G.R., and Lindenberg, M.K., 

2006, Water-quality conditions in Upper Klamath 

Lake, Oregon, 2002–2004: U.S. Geological Survey Scientific Investigations Report 2006-5209, 

52 p. (Also available at 

http://pubs.usgs.gov/sir/2006/5209/.) 

Wood, T.M., Cheng, R.T., Gartner, J.W., Hoilman,

G.R., Lindenberg, M.K., and Wellman, R.E., 2008, 

Modeling hydrodynamics and heat transport in Upper Klamath Lake, Oregon, and implications for

water quality: U.S. Geological Survey Scientific Investigations Report 2008–5076, 48 p. (Also available

at http://pubs.usgs.gov/sir/2008/5076/)

Zhang, L., Gu, X., Fan, C., Shang, J., Shen, Q., 

Wang, Z., and Shen, J., 2010, Impact of different 

benthic animals on phosphorus dynamics across 

the sediment-water interface: Journal of Environmental Sciences, v. 22, no. 11, p. 1674–1682.
71

Appendix A. R-Script Code for Upper Klamath Lake Phosphorus Model

#Load observed, bi-weekly input data

setwd("E:/users/swherry/UKL_TMDL_Project/TMDL_model/TMDL_model_R/")

input<-read.csv("input_data.csv",header=TRUE)

#Convert input Date column to date class

input$Date<-as.Date(input$Date,"%m/%d/%Y")

#Define iterative Chlorophyll a model termination criteria

iteration<-50

difference<-10^-16

#Define model parameters

kg<-0.012 #Gross removal (sedimentation) rate of non-algal P from water column (1/day)

ph_sigma<-0.25 #pH spatial standard deviation

ph2_int<-7.93101869324701 #pH vs chlorophyll a intercept

ph2_chla<-0.533758530498359 #pH vs chlorophyll a slope

ph2_julian<--0.00596604468236423 #pH vs Julian day slope

jdmin<-200 #Minimum Julian day for model

ph_min<-7.8 #Minimum predicted pH

beta<-0.0097 #Algal light extinction coefficient (1/m/ppb)

alpha<-1.2 #Background light extinction (1/m)

ihalf<-214 #Onset of light saturation (uE/m2-sec)

gmax<-1.2 #Max algal growth rate (1/day)

rmax<-0.06 #Max algal respiration rate (1/day)

topt<-20 #Optimum temp for algal growth (deg-C)

tzero<-14 #Temperature when zero algal growth (deg-C)

yp<-0.6 #Algal cell P quota (mg P/mg chlorophyll a)

chla_min<-8 #Minimum predicted chlorophyll a (ppb)

zsed<-10 #Active sediment thickness (cm)

sed_dens<-0.12 #Sediment bulk density (g/cm3)

uburial<-1.4 #Burial velocity (mm/yr)

latitude<-42.5 #Latitude

reflect<-0.1 #Light reflection fraction

par<-0.714285714285714 #PAR/langley conversion

wc_p0<-74 #Initial water column total P concentration (ppb)

deltat<-14 #Time increment (days)

area<-271 #Total lake area (km2)

v0<-622 #Initial lake volume (hm3)

khalf<-0.06 #Phosphorus half saturation constant for Michaelis-Menten kinetics

#Request inputs from the model user

reduc<-as.numeric(readline("How much load reduction (%) is there? "))/100 

#Watershed load Reduction (%)

ualgae<-as.numeric(readline("What is the algal settling velocity (m/d)? ")) 

#Algal settling rate (m/d)

ph_threshold<-as.numeric(readline("What is pH* (pH at half max recycle)? ")) 

#pH at half the max recycle rate

ks_ph<-as.numeric(readline("What is the maximum total P recycle rate, KR? ")) 

#Maximum total P recycle rate (1/yr)

sed_p0<-as.numeric(readline("What is the initial sediment total P concentration (Original=1000; Recalc=335)? 

")) #Initial sediment total P concentration (mg/kg)

#Basic lake characteristics calculations

tmdl_load<-input$Load*(1-reduc) #Phosphorus load for given TMDL reduction

load_tot<-(tmdl_load+input$Precip_load) #Total P load from inflow and precip

v_0<-input$End_Storage-input$Change_Storage #Initial lake volumes for each time step

v_mean<-input$End_Storage-(input$Change_Storage/2) #Mean lake volume for each time step

z_0<-v_0/input$Area #Initial lake depth for each time step

z_mean<-v_mean/input$Area #Mean lake depth for each time step

q_out<-input$Outflow/deltat #Lake outflow for each time step

#Calculated parameters from defined parameters

sed_mass<-sed_dens*zsed*10 #Sediment soil mass (kg/m2)
72

sinit<-sed_mass*area*sed_p0 #Initial sediment total P reservoir (kg)

minit<-wc_p0*v_mean[1] #Initial water column total P (kg)

step<-deltat/365.25 #Time step (year)

kb<-uburial/10/zsed #Burial rate (1/yr)

#Define the unknown variables the model will be solving for

Mo<-vector(mode="numeric",length=nrow(input)) #Initial water column total P mass for each time 

step

Mt<-vector(mode="numeric",length=nrow(input)) #Final water column total P mass for each time step 

(14 days later)

Mavg<-vector(mode="numeric",length=nrow(input)) #Average water column total P mass over the time 

step

Cavg<-vector(mode="numeric",length=nrow(input)) #Average water column total P concentration over the 

time step

C_NA<-vector(mode="numeric",length=nrow(input)) #Non-algal water column total P concentration over 

the time step

C_alg<-vector(mode="numeric",length=nrow(input)) #Algal water column total P concentration

Sed_alg<-vector(mode="numeric",length=nrow(input)) #Mass of algal P settling out of the water column

Sed_NA<-vector(mode="numeric",length=nrow(input)) #Mass of non-algal P settling out of the water column

Mout<-vector(mode="numeric",length=nrow(input)) #Mass of total P flowing out of the lake

Mrec<-vector(mode="numeric",length=nrow(input)) #Mass of total P recycled into the water column from 

the sediment

delM<-vector(mode="numeric",length=nrow(input)) #Change in mass of total P over the time step

So<-vector(mode="numeric",length=nrow(input)) #Initial sediment total P mass for each time step

St<-vector(mode="numeric",length=nrow(input)) #Final sediment total P mass for each timstep (14 

days later)

Savg<-vector(mode="numeric",length=nrow(input)) #Average sediment total P mass over the time step

Sb<-vector(mode="numeric",length=nrow(input)) #Mass of buried total P in the sediment

Chla<-vector(mode="numeric",length=nrow(input)) #Chlorophyll-a concentration 

pH<-vector(mode="numeric",length=nrow(input)) #pH

W_check<-vector(mode="numeric",length=nrow(input))

S_check<-vector(mode="numeric",length=nrow(input))

#Request nutrient and light limitation information from model user

fp_method<-as.numeric(readline("Which phosphorus limitation approach should be used (1=%Non-alg P; 2=MichaelisMenten)? "))

fl_method<-as.numeric(readline("Which light limitation approach should be used (1=Average; 2=Depth-integrated)? 

"))

#Define a function that runs the chlorophyll a model with a water column total P mass input

chla_model<-function(m) {

n<-1 

while(n<iteration) {

if(n==1) {

chla<-chla_min

chla_calc<-chla_min} else {chla<-mean(c(chla,chla_calc[n-1]))}

#Calculate the average water column total P concentrations (total and non-algal)

cavg<-m/v_mean[i]

na<-cavg-yp*chla

#Calculate the light extinction coefficient

Ec<-alpha+beta*chla

#Calculate the amount of daylight, in hours, based on Julian day and geography

if(input$Julian[i]<8) d<-365+input$Julian[i]-7 else d<-input$Julian[i]

declin<-0.4102*sin((2*pi/365)*(d-80.25))

day_length<-7.72*acos(-tan((2*pi/360)*latitude)*tan(declin))

#Calculate light intensity at the lake surface

Io<-input$SolarRad[i]*par*(24/day_length)*(1-reflect)

#Determine average light intensity over the water column

Im<-(Io-Io*exp(-z_mean[i]*Ec))/(z_mean[i]*Ec)

#Calculate the temperature adjustment factor
73

ft<-max(0,(input$WtrTemp[i]-tzero)/(topt-tzero))

#Determine the growth rate

gt<-max(0,gmax*min(1,ft))

#Determine the combination respiration, settling and outflow term

r<-max(0.01,rmax*ft+(ualgae/z_mean[i])+(q_out[i]/v_mean[i]))

#Determine the light and nutrient limitation factors based on the user-specified methods

fl<-

ifelse(fl_method==1,(Im/ihalf)/sqrt(1+(Im/ihalf)^2)*(day_length/24),((day_length/24)/(Ec*z_mean[i]))*log(((Io/i

half)+sqrt(1+((Io/ihalf)^2)))/((Io*exp(-z_mean[i]*Ec)/ihalf)+sqrt(1+((Io*exp(-z_mean[i]*Ec)/ihalf)^2)))))

fp<-ifelse(fp_method==1,min(1,na/cavg),max(0,min(1,(na/1000)/((na/1000)+khalf))))

#Calculate the chlorophyll a concentration such that its minimum is never less than 8 ppb

#and its maximum never exceeds the cell quota with respect to available phosphorus

chla_calc[n]<-min(max(chla_min,ifelse(gt==0,chla_min,gt*fl*fp*chla/r)),cavg/yp)

#If chlorophyll concentration hasn't changed by more than the threshold since the last iteration, 

#finish the iterative process otherwise continue

if (abs(chla-chla_calc[n])<difference) n<-iteration else n<-n+1

} 

chla_calc[length(chla_calc)]

} 

#Begin the model loop, the outer loop runs through each time step "i"

for(i in 1:nrow(input)) {

#Set the initial water column and sediment total P masses, the calculated values for the first

#time step i=1, and end values from the previous timstep for i>1

if(i==1) { 

m_o<-minit

s_o<-sinit} else {

m_o<-Mt[i-1] 

s_o<-St[i-1]}

#Determine the corresponding chlorophyll a concentration for m_o

chla0<-chla_model(m_o)

#pH Model (Dependent on chlorophyll a concentration)

ph<-max(ph_min,ph2_int+ph2_chla*log(chla0)+ph2_julian*max(jdmin,input$Julian[i]))

#Calculate the recycled mass of total P from the sediment at t=0 

m_rec0<-((1-pnorm((ph_threshold-ph)/ph_sigma))*ks_ph*s_o*step)/deltat

#Calculate the average water column non-algal P concentration at t=0

c_na<-m_o/v_mean[i]-chla0*yp

#Calculate the mass of algal P settling out of the water column

sed_na<-c_na*v_mean[i]*kg

#Calculate the mass of non-algal P settling out of the water column 

sed_alg<-chla0*yp*input$Area[i]*ualgae

#Calculate the mass of total P flowing out of the lake

m_out<-((m_o/v_mean[i])*input$Outflow[i])/deltat

#Sediment Model

#Calculate the mass of buried total P in the sediment

s_b<-(kb*s_o*step)/deltat

#Calculate the first set of Runge-Kutta factors, k11 and k12, or the first

#approximated change in mass

k11<-(load_tot[i]/deltat+m_rec0-sed_alg-sed_na-m_out)*deltat

k12<-(sed_alg+sed_na-m_rec0-s_b)*deltat
74

#Calculate the variables modified by k11 and k12

#Determine the corresponding chlorophyll a concentration for m_o

chla1<-chla_model(m_o+0.5*k11)

ph1<-max(ph_min,ph2_int+ph2_chla*log(chla1)+ph2_julian*max(jdmin,input$Julian[i]))

m_rec1<-((1-pnorm((ph_threshold-ph1)/ph_sigma))*ks_ph*(s_o+0.5*k12)*step)/deltat

sed_na1<-kg*((m_o+0.5*k11)-yp*chla1*v_mean[i])

sed_alg1<-chla1*yp*input$Area[i]*ualgae

m_out1<-((m_o+0.5*k11)/v_mean[i])*input$Outflow[i]/deltat

s_b1<-kb*(s_o+0.5*k12)*step/deltat

#Calculate the second set of Runge-Kutta factors, k21 and k22, or the second

#approximated change in mass

k21<-(load_tot[i]/deltat+m_rec1-sed_alg1-sed_na1-m_out1)*deltat

k22<-(sed_alg1+sed_na1-m_rec1-s_b1)*deltat

#Calculate the variables modified by k21 and k22

#Determine the corresponding chlorophyll a concentration for m_o

chla2<-chla_model(m_o+0.5*k21)

ph2<-max(ph_min,ph2_int+ph2_chla*log(chla2)+ph2_julian*max(jdmin,input$Julian[i]))

m_rec2<-((1-pnorm((ph_threshold-ph2)/ph_sigma))*ks_ph*(s_o+0.5*k22)*step)/deltat

sed_na2<-kg*((m_o+0.5*k21)-yp*chla2*v_mean[i])

sed_alg2<-chla2*yp*input$Area[i]*ualgae

m_out2<-((m_o+0.5*k21)/v_mean[i])*input$Outflow[i]/deltat

s_b2<-kb*(s_o+0.5*k22)*step/deltat

#Calculate the third set of Runge-Kutta factors, k31 and k32, or the third

#approximated change in mass

k31<-(load_tot[i]/deltat+m_rec2-sed_alg2-sed_na2-m_out2)*deltat

k32<-(sed_alg2+sed_na2-m_rec2-s_b2)*deltat

#Calculate the variables modified by k31 and k32

#Determine the corresponding chlorophyll a concentration for m_o

chla3<-chla_model(m_o+k31)

ph3<-max(ph_min,ph2_int+ph2_chla*log(chla3)+ph2_julian*max(jdmin,input$Julian[i]))

m_rec3<-((1-pnorm((ph_threshold-ph3)/ph_sigma))*ks_ph*(s_o+k32)*step)/deltat

sed_na3<-kg*((m_o+k31)-yp*chla3*v_mean[i])

sed_alg3<-chla3*yp*input$Area[i]*ualgae

m_out3<-((m_o+k31)/v_mean[i])*input$Outflow[i]/deltat

s_b3<-kb*(s_o+k32)*step/deltat

#Calculate the first set of Runge-Kutta factors, k41 and k42, or the fourth

#approximated change in mass

k41<-(load_tot[i]/deltat+m_rec3-sed_alg3-sed_na3-m_out3)*deltat

k42<-(sed_alg3+sed_na3-m_rec3-s_b3)*deltat

#Calculate the final water column total P mass 14 days later based on the four approximations;

#the second and third approximations are given the greatest weight

m_t<-m_o+(1/6)*(k11+2*k21+2*k31+k41)

#Calculate the average water column total P mass over the time step

m_avg<-mean(c(m_o,m_t)) 

#Calculate the average water column total P concentration over the time step

c_avg<-m_avg/v_mean[i]

#Calculate the final sediment total P mass 14 days later

s_t<-s_o+(1/6)*(k12+2*k22+2*k32+k42)

#Calculate the change in mass of total P over the time step

delm<-m_t-m_o 

#Sediment Model

#Calculate the change in mass of total P over the time step

dels<-s_t-s_o

#Calculate the average sediment total P mass over the time step

s_avg<-mean(c(s_o,s_t))
75

#Calculate the average mass of total P flowing out of the lake over the time step

m_out_avg<-(c_avg*input$Outflow[i])/deltat

#Calculate the average mass of buried total P in the sediment over the time step

s_b_avg<-(kb*s_avg*step)/deltat

#Determine average chlorophyll a concentration as based on average water column total P mass over the time step

chla_avg<-chla_model(m_avg)

#Determine the average pH over the time step

ph_avg<-max(ph_min,ph2_int+ph2_chla*log(chla_avg)+ph2_julian*max(jdmin,input$Julian[i]))

#Calculate the average total P recycled mass over the time step

m_rec_avg<-((1-pnorm((ph_threshold-ph_avg)/ph_sigma))*ks_ph*s_avg*step)/deltat

#Calculate the average water column non-algal P concentration over the time step

c_na_avg<-c_avg-chla_avg*yp

#Calculate the average mass of algal P settling out of the water column over the time step

sed_na_avg<-c_na_avg*v_mean[i]*kg

#Calculate the average mass of non-algal P settling out of the water column over the time step

sed_alg_avg<-chla_avg*yp*input$Area[i]*ualgae

#Check the water column mass balance to see if it equals or is close to zero

WC_check<-load_tot[i]+((m_rec_avg-sed_alg_avg-sed_na_avg-m_out_avg)*deltat)-delm 

#Check the sediment mass balance to see if it equals or is close to zero

Sed_check<-((sed_alg_avg+sed_na_avg-m_rec_avg-s_b_avg)*deltat)-dels

#Assign the variables determined through iteration for each time step

Mo[i]<-m_o

Mt[i]<-m_t

Mavg[i]<-m_avg

Cavg[i]<-c_avg

C_NA[i]<-c_na_avg

C_alg[i]<-c_avg-c_na_avg

Sed_alg[i]<-sed_alg_avg*deltat

Sed_NA[i]<-sed_na_avg*deltat

Mout[i]<-m_out_avg*deltat

Mrec[i]<-m_rec_avg*deltat

So[i]<-s_o

St[i]<-s_t

Savg[i]<-s_avg

Sb[i]<-s_b_avg*deltat

Chla[i]<-chla_avg

pH[i]<-ph_avg

W_check[i]<-WC_check

S_check[i]<-Sed_check

} 

model_output<-

data.frame(Date=input$Date,Mo,Mt,Mavg,Cavg,C_NA,C_alg,Sed_alg,Sed_NA,Mout,Mrec,So,St,Savg,Sb,Chla,pH)
Wood and others— Technical Evaluation of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon—Open-File Report 2013-1262

ISSN 2331-1258 (online)

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