Prepared in cooperation with the Bureau of Reclamation

Revision and Proposed Modification of a Total Maximum 

Daily Load Model for Upper Klamath Lake, Oregon

Scientific Investigations Report 2015–5041

U.S. Department of the Interior

U.S. Geological Survey

Revision and Proposed Modification of a 

Total Maximum Daily Load Model for 

Upper Klamath Lake, Oregon

By Susan A. Wherry, Tamara M. Wood, and Chauncey W. Anderson

Prepared in cooperation with the Bureau of Reclamation

Scientific Investigations Report 2015–5041

U.S. Department of the Interior

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

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

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

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

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

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

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

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

U.S. Government.

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

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

Suggested citation:

Wherry, S.A., Wood, T.M., and Anderson, C.W., 2015, Revision and proposed modification of a total maximum daily 

load model for Upper Klamath Lake, Oregon: U.S. Geological Survey Scientific Investigations Report 2015–5041, 55 p., 

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

ISSN 2328-0328 (online)
iii

Contents

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

Comparison of Model Results.............................................................................................................1

Improvements to Modeling Cyanobacterial Dynamics ..................................................................1

Introduction ....................................................................................................................................................3

Purpose and Scope .......................................................................................................................................5

Review of Phase 1 Results ..........................................................................................................................5

Changes to Algal Submodel ...............................................................................................................6

New Information—Internal Loading and Benthic Invertebrates .................................................6

203-Year Simulation and Sensitivity Analysis...................................................................................9

Phase 2 Model Updates..............................................................................................................................10

Updated CSTR Model Calibration ....................................................................................................10

Reevaluate Phase 1 Test Statistics Using Updated Algal Settling Velocity 

Calibration, 1991–98 .....................................................................................................10

Step 1—Algal Settling Velocity and Half-Saturation Constant ................................10

Step 2—Internal Recycling Parameters.......................................................................11

Reevaluation of Phase I Test Statistics Using Extended Data, 1991–2010.......................11

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

Step 2—Algal Settling Velocity and Half-Saturation Constant.................................12

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

Step 4—Internal Recycling Parameters.......................................................................12

Model Performance ..................................................................................................................13

Long-Term Simulations ......................................................................................................................27

Lake Metabolism as a Substitute for Algal Growth and Respiration .........................................34

Calculation of Upper Klamath Lake Metabolism from Continuous Dissolved 

Oxygen, 2006–08 ...........................................................................................................36

Lake Metabolism Calculations .......................................................................................36

Multivariate Regression Models of Lake Metabolism................................................37

Incorporation of Upper Klamath Lake Metabolism Submodel into CSTR Model............40

Incorporation into CSTR Model......................................................................................40

Step 1—Algal Settling Velocity .............................................................................40

Step 2—Internal Recycling Parameters..............................................................40

Model Performance ................................................................................................41

Prediction of Upper Klamath Lake Metabolism from Meteorological 

Variables, 2006–08 ..............................................................................................44

Discussion and Suggestions for Model Improvements ........................................................................48

Model Updates....................................................................................................................................48

Continuously Stirred Tank Reactor Model Calibration ........................................................48

Metabolism as an Independent Predictor of Bloom Dynamics .........................................48

Future Model Needs...........................................................................................................................49

References Cited .........................................................................................................................................51

Appendix A. Derivation of Net Settling Velocity ....................................................................................55
iv

Figures

1. Map showing location of meteorological stations and continuous water-quality

monitoring and sampling sites, Upper Klamath Lake, Oregon..............................................4

2. Graphs showing simulated and measured water column total phosphorus

concentration, chlorophyll a concentration, and pH, Upper Klamath Lake, Oregon,

April 1991–September 2010.......................................................................................................26

3. Graphs showing 19-year average water-column total phosphorus concentration

under simulation conditions of 0-percent reduction and 40-percent reduction in

external phosphorus loads to Upper Klamath Lake, Oregon...............................................30

4. Graphs showing 19-year average sediment total phosphorus masses under

simulation conditions of 0-percent reduction and 40-percent reduction in external

phosphorus loads to Upper Klamath Lake, Oregon ..............................................................31

5. Graphs showing water-column total phosphorus concentration during the last 

19 years of a 209-year simulation (1991–2199) under conditions of 0-percent 

reduction and 40-percent reduction in external loads to Upper Klamath Lake, 

Oregon ..........................................................................................................................................32

6. Graphs showing water-column chlorophyll a concentration during the last 19 years

of a 209-year simulation (1991–2199) under conditions of 0-percent reduction and

40-percent reduction in external loads to Upper Klamath Lake, Oregon..........................33

7. Graphs showing Local Regression (LOESS) smoothed respiration as a function of

regression model of metabolism predictors, Upper Klamath Lake, Oregon, during

early, middle, and late parts of summer seasons 2006–08...................................................39

8. Graphs showing biweekly simulated and measured water-column total phosphorus

concentration, chlorophyll a concentration, and pH, Upper Klamath Lake, Oregon,

2006–08 .........................................................................................................................................41

9. Graphs showing simulated and measured water-column total phosphorus

concentration, chlorophyll a concentration, and (C) pH, Upper Klamath Lake,

Oregon, 2006–08..........................................................................................................................42

10. Graphs showing simulated and measured water-column pH in 2006, 2007, and

2008, Upper Klamath Lake, Oregon..........................................................................................43

11. Graphs showing measured values, calculated using dissolved oxygen data 

from nine, continuous, shallow-sites, U.S. Geological Survey monitors, and

simulated values of the seasonal component of the regression model of

metabolism, Upper Klamath Lake, Oregon, 2006–08.............................................................44

12. Graphs showing respiration and production values calculated from measured 

data and simulated values of the weekly component of the regression model of

metabolism, Upper Klamath Lake, Oregon, 2006–08.............................................................45

13. Graphs showing respiration and production values calculated from measured data

and simulated values from the combined seasonal and weekly components of the

regression model of metabolism, Upper Klamath Lake, Oregon, 2006–08 ........................46

14. Graphs showing comparison of the calculated and simulated net primary 

production to measured chlorophyll a values from Klamath Tribes and 

U.S. Geological Survey monitoring programs, Upper Klamath Lake, 

Oregon, 2006–08..........................................................................................................................47
v

Tables

1. Descriptions of model terms used in this report .....................................................................7

2. Descriptions, values, and sources of parameters used in the phosphorus total

maximum daily load (TMDL) model, Upper Klamath Lake, Oregon......................................8

3. pH model coefficients and statistics between the original total maximum daily load

(TMDL) model for Upper Klamath Lake, Oregon, and the recalibration considering

models with and without a Julian day component ...............................................................12

4. Descriptive statistics between measured and calculated nonalgal phosphorus

mass lost to sedimentation during March–May, Upper Klamath Lake, Oregon,

1991–2010 .....................................................................................................................................13

5. Performance statistics for total phosphorus obtained over the indicated validation

periods for model versions discussed in the report .............................................................14

6. Performance statistics for chlorophyll a obtained over the indicated validation

periods for model versions discussed in the report .............................................................16

7. Performance statistics for pH obtained over the indicated validation periods for

model versions discussed in the report..................................................................................18

8. Performance statistics for total phosphorus obtained over the indicated May–July

validation periods for model versions discussed in the report...........................................20

9. Performance statistics for chlorophyll a obtained over the indicated May–July

validation periods for model versions discussed in the report...........................................22

10. Performance statistics for pH obtained over the indicated May–July validation

periods for model versions discussed in the report .............................................................24

11. Metrics describing the long-term simulations for various model versions discussed

in the report .................................................................................................................................28

12. U.S. Geological Survey continuous monitoring sites used for metabolism analysis,

Upper Klamath Lake, Oregon, 2006–08....................................................................................37

13. Multivariate regression equations...........................................................................................38

14. Goodness-of-fit statistics for the combined (seasonal plus weekly component)

regression models for total production and community respiration, Upper Klamath

Lake, Oregon, 2006–08 ...............................................................................................................46
vi

Conversion Factors

Inch/pound to International System of Units

Multiply By To obtain

Concentration

parts per billion (ppb) 1.0 micrograms per liter (μg/L)

International System of Units to inch/pound

Multiply By To obtain

Length

meter (m) 3.281 foot (ft) 

kilometer (km) 0.6214 mile (mi)

Area

square meter (m2

) 0.0002471 acre

square meter (m2

) 10.76 square foot (ft2)

Volume

liter (L) 33.82 ounce, fluid (fl. oz)

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)

square meter per day (m2

/d) 10.764 square foot per day (ft2/d)

Mass

milligram (mg) 0.0000353 ounce, avoirdupois (oz)

gram (g) 0.03527 ounce, avoirdupois (oz)

kilogram (kg) 2.205 pound avoirdupois (lb)

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

F = (1.8 × °C) +32.

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

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

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

Supplemental Information

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

or micrograms per liter (μg/L).

Abbreviations

AIC Akaike’s Information Criterion

AFA Aphanizomenon flos-aquae 

CR community respiration

CSTR continuously stirred tank reactor

DO dissolved oxygen

GPP gross primary production

MET meteorological

NPP net primary production

NS Nash-Sutcliffe statistic

R correlation coefficient

R2 coefficient of determination

RMSE root mean square error

SOD sediment-oxygen demand

TMDL total maximum daily load

TP total phosphorus

UKL Upper Klamath Lake

USGS U.S. Geological Survey

WQ water quality

Revision and Proposed Modification of a Total Maximum 

Daily Load Model for Upper Klamath Lake, Oregon

By Susan A. Wherry, Tamara M. Wood, and Chauncey W. Anderson 

Executive Summary

This report presents Phase 2 of the review and 

development of the mass balance water-quality model, 

originally developed in 2001, that guided establishment of the 

phosphorus (P) total maximum daily load (TMDL) for Upper 

Klamath and Agency Lakes, Oregon. The purpose of Phase 2 

was to incorporate a longer (19-year) set of external phosphorus 

loading data into the lake TMDL model than had originally 

been available, and to develop a proof-of-concept method 

for modeling algal mortality and the consequent decrease in 

chlorophyll a that had not been possible with the 2001 TMDL 

model formulation. 

Using the extended 1991–2010 external phosphorus 

loading dataset, the lake TMDL model was recalibrated 

following the same procedures outlined in the Phase 1 review. 

The version of the model selected for further development 

incorporated an updated sediment initial condition, a numerical 

solution method for the chlorophyll a model, changes to light 

and phosphorus factors limiting algal growth, and a new 

pH-model regression, which removed Julian day dependence 

in order to avoid discontinuities in pH at year boundaries. 

This updated lake TMDL model was recalibrated using the 

extended dataset in order to compare calibration parameters 

to those obtained from a calibration with the original 7.5-year 

dataset. The resulting algal settling velocity calibrated from 

the extended dataset was more than twice the value calibrated 

with the original dataset, and, because the calibrated values of 

algal settling velocity and recycle rate are related (more rapid 

settling required more rapid recycling), the recycling rate also 

was larger than that determined with the original dataset. These 

changes in calibration parameters highlight the uncertainty in 

critical rates in the Upper Klamath Lake TMDL model and 

argue for their direct measurement in future data collection to 

increase confidence in the model predictions. 

Comparison of Model Results

We compared model predictions of steady-state 

water-column conditions in response to the 40-percent 

reduction in external phosphorus loads stipulated in the 2001 

TMDL, using the original model and the recalibrated version 

with revisions as described in this report. The qualitative 

prediction of reduced concentrations in response to external 

load reductions, occurring over a period on the order of 

decades, is consistent across model versions and is a reliable 

feature of the equations that describe a coupled sediment 

and water-column system. The quantitative predictions of 

the original and revised models differ in time required to 

achieve steady state (50 years and 19 years, respectively) and 

the associated steady-state concentration of water column 

total P (40 and 74 parts per billion, or ppb, respectively), with 

chlorophyll a concentrations being similar (25 and 27 ppb, 

respectively). The quantitative predictions, however, have a 

high degree of uncertainty because they depend on calibration 

parameters that vary greatly with changes in assumptions 

or the calibration dataset, and neither set of calibration 

parameters resulted in a model with good fit statistics outside 

of the calibration period. If a model that can predict water 

column concentrations in response to TMDL implementation 

with less uncertainty is required, several improvements to the 

current model can be pursued and are briefly explained. 

Improvements to Modeling Cyanobacterial 

Dynamics

After comparing recalibrated parameters and results, we 

incorporated algal mortality (bloom decline) into the algal 

submodel. The original lake TMDL algal submodel relied on 

a mass balance that included terms for growth and respiration, 

the latter a combined term that was intended to include 

mortality. Both terms were modified throughout the season 

by the product of limitation factors that were determined 

by input data averaged over the lake area and biweekly 

time step. This common modeling approach is relatively 

simple and can be run with low-resolution data (biweekly); 

however, it introduces error when used to describe a highly 

variable system such as Upper Klamath Lake instead of a 

small, homogeneous parcel of water, and does not capture 

the cumulative effects of short-term variability. Our proofof-concept alternative makes use of high-resolution data 

to better capture spatial and temporal dynamics of algal 

growth and mortality by substituting the original model terms 

with metabolic terms, gross primary production (GPP) and 

community respiration (CR), as calculated from continuously 

monitored dissolved oxygen (DO) data. 
2 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

We calculated GPP and CR during 3 years from an 

hourly mass balance of dissolved oxygen that included 

reaeration and sediment oxygen demand (SOD), and assumed 

that GPP ceased during times of darkness. To minimize 

effects of vertical transport of DO, we evaluated only the 

nine shallow sites in operation during 2006–08. Hourly 

data were accumulated each day and stoichiometric ratios 

relating oxygen to carbon, and carbon to chlorophyll a, 

were used to convert daily estimates of production and 

consumption of DO to changes in chlorophyll a concentration. 

Lake‐averaged GPP and CR were substituted for algal growth 

and respiration terms, respectively, in the model, and the 

model was recalibrated. The initial runs were promising in that 

the mid-season bloom declines were captured and statistics 

describing simulated chlorophyll a concentration (coefficient 

of determination [R2]=0.19–0.43) were comparable to those 

found over the same period for the revised and recalibrated 

model as previously described (R2=0.25–0.27). Because 

this metabolism-based approach was able to capture bloom 

declines and performed as well as the TMDL model over 

3 years of data, we developed a statistical model for predicting 

GPP and CR based on independent variables for which 

long-term records existed.

Regression analysis was used to relate GPP and CR 

to meteorological (wind and solar radiation) variables and 

water temperature. The goal was to derive relations that could 

be incorporated in the TMDL model equations to improve 

accuracy in the prediction of the state variables in response to 

implementation of TMDL constraints on external phosphorus 

loading. Wind, solar radiation, and water temperature were 

used as explanatory variables for model predictions of the 

response to changes in external loads because they are not 

expected to change due to changing phosphorus conditions. 

As a result of an initial, exploratory analysis, GPP and 

CR were determined to correlate to explanatory variables 

differently depending on season; therefore, season-specific 

regression models were developed. Each model comprised 

a seasonal component determined from biweekly data and 

a weekly component determined from daily values that 

were calculated from hourly measurements. Potentially 

explanatory variables tested in the seasonal model included 

14-day lakewide average chlorophyll a, solar radiation, 

and temperature. Potentially explanatory variables tested in 

the weekly model included wind speed, measures of water 

temperature deviation from a seasonal fit and rate of change, 

water-column stratification, and solar radiation deviation from 

a seasonal fit. The results of the combined models generally 

were good, with R2

 greater than or equal to 0.53 for CR and 

R2 greater than or equal to 0.62 for GPP. The next step in this 

analysis would be incorporation of the regression model of 

metabolism in the lake TMDL model in place of algal growth 

and respiration; however, that step would require additional 

data extraction and model calibration that was beyond the 

scope of this project. 

Based on the results of the recalibration of the 

lake TMDL model with the updated, extended dataset 

and explorations of alternative methods for modeling 

chlorophyll a, we propose five steps that could be explored to 

improve the lake TMDL model in the future: 

1. Use a daily time step for modeling. A daily time step 

is straightforward to implement and is reasonable in 

the context of the following steps that make use of 

continuously (hourly) collected data.

2. Use a spatially explicit model to interpolate between 

sampling sites and to accurately calculate (1) the 

lakewide average of variables as a function of the 

measured values at sampling sites, and (2) the 

concentrations in the lake outflow, so as to improve 

the accuracy of estimated downstream P and 

chlorophyll a loads.

3. Substitute GPP and CR for algal growth and respiration 

terms in the model. We have shown the proof-of-concept 

for using continuously measured DO to calculate the 

lake metabolism terms. Further work needs to be done to 

improve the empirical models that enable the prediction 

of GPP and CR from allowable explanatory variables 

before GPP and CR are formally substituted for the algal 

growth and respiration terms of the model.

4. Expand equations used to model internal phosphorus 

loading. Specifically, create explicit, separate terms 

for benthic invertebrate excretion, groundwater 

flux, diffusive flux, desorption at high pH, and rapid 

bacterial recycling. 

5. Improve the model treatment of net settling velocity. The 

treatment of net settling velocity as a single calibration 

parameter is not satisfactory for reasons related to 

the scaling-up process in general (appendix A), but 

particularly because the net settling velocity varies 

across time with the proportion of buoyant and nonbuoyant colonies in the water column.

 Two complementary approaches are proposed to 

improve net settling velocity. The first approach relies 

on using continuous monitor data to estimate daily 

depth-integrated algal biomass values at a monitor 

site. Estimates at individual continuous monitoring 

sites would be combined to create a lakewide average. 

Daily estimates of lakewide biomass storage in 

combination with daily values of net primary production 

(NPP) allow the daily settling flux to be calculated 

by difference. This approach can be pursued only 

in combination with the calculation of GPP and CR 

(step 3). The second approach relies on field work to 

measure rising and falling velocities and to estimate 

the fraction of rising and falling colonies as a function 
Introduction 3

of environmental variables. After measurements are 

obtained, one-dimensional water-column models that 

investigate only the interaction between buoyancy 

and stratification, and three-dimensional models that 

add horizontal transport, would be used to assist the 

scaling-up process.

The combination of the suggested improvements is 

expected to result in a modeling framework that can better 

account for spatial and temporal variations in water quality in 

the lake, while using independent variables to predict future 

algal-mortality events. Ultimately, this would provide better 

tools to represent the time frames needed to achieve desired 

outcomes based on prescribed reductions in external loading.

Introduction 

Upper Klamath Lake (UKL) in south-central Oregon 

(fig. 1) has a long history of harmful cyanobacterial 

blooms, with biomass dominated by the cyanobacterium

Aphanizomenon flos-aquae (AFA), and associated poor 

water quality, dating from at least the early 1990s (Wood 

and others, 1996; Eilers and others, 2004; Wood and others, 

2006). Among the problems associated with these blooms are 

(1) episodically low concentrations of dissolved oxygen (DO), 

high pH, and high ammonia concentrations that negatively 

affect endangered fish in the lake (Burdick and others, 2009); 

(2) the production of cyanobacterial toxins that could have 

negative health effects on wildlife and humans (VanderKooi 

and others, 2010; Eldridge and others, 2012); and (3) the 

export of nutrients, cyanobacteria, and associated toxins, 

and organic carbon to downstream reaches (Sullivan and 

others, 2011). In 2002, Oregon Department of Environmental 

Quality promulgated a Total Maximum Daily Load (TMDL) 

(Oregon Department of Environmental Quality, 2002) that 

set allowable limits on the amount of phosphorus entering the 

lake from external sources as a means to reduce the frequency, 

magnitude, and extent of algal blooms and thereby improve 

water quality. 

The TMDL was based on a lake nutrient budget 

developed by Kann and Walker (1999), and used a 

spreadsheet-based model that approximated lake average 

concentrations of total phosphorus (TP) for the combined 

Upper Klamath and Agency Lakes, assuming that the lake 

behaved as a continuously stirred tank reactor (CSTR) 

(Walker, 2001). The model boundary conditions were derived 

from external loads estimated from bi-weekly samples, and 

the model was calibrated using biweekly in-lake monitoring 

data and internal recycling data derived from mass-balance 

calculations (Walker and others, 2012). Model outputs 

included TP concentration (water-column and sedimentbound), chlorophyll a concentration, and pH.

In 2011, the U.S. Geological Survey (USGS) began a 

multi-phase technical review of the UKL TMDL model in 

order to evaluate the strengths and limitations of the model 

and to determine whether improvements could be made using 

information and data derived from studies completed since 

the development of the model. Phase 1 (Wood and others, 

2013) included a literature review of major hydrological and 

biogeochemical processes relevant to UKL, incorporation 

of recent field and laboratory studies, and a technical review 

and evaluation of the spreadsheet-based CSTR model. 

Strengths and weaknesses of the existing CSTR approach 

were described. One of the primary advantages of the existing 

model is its applicability to long-term simulations that 

are computationally achievable and that can give insights 

into large-scale lake responses to different management 

approaches. A limitation is its coarse representation of the 

lake’s water quality, both spatially and temporally, compared 

to scales that likely operate in the lake during peak growing 

season. Temporal wavelet analysis of hourly water quality 

and meteorological data in the summer of 2006 indicated 

variability at time scales as low as 5.8–7.6 days, and 

spatial wavelet analysis of phycocyanin fluorescence data 

(an indicator of cyanobacteria) collected along two lake 

transects during bloom events in 2012 showed variability at 

distances as small as 0.25 km (Eldridge and others, 2014). 

Because a whole-lake average is not necessarily the same 

as the concentration in the outflow, the model may not 

accurately predict water quality exported from the lake. 

Downstream reaches also have significant water quality 

concerns—specifically, elevated ammonia concentrations 

and anoxic conditions owing to decomposing AFA exported 

from UKL (Sullivan and others, 2011). A model that provides 

more accurate (throughout the year) and higher resolution 

boundary conditions to downstream reaches would be useful 

for water management in the Klamath River and Klamath 

Irrigation Project. 

In Phase 1, improvements were made to the growth 

limitation factors used in the original spreadsheet model, and 

the model then was re-scripted in R (R Core Team, 2014) 

with a Runge-Kutta solution method and re-calibrated. The 

lake sediment initial conditions were updated with the most 

recent information (Simon and others, 2009; Simon and 

Ingle, 2011). Modifications to the R-script model were used 

to explore the long-term behavior of model simulations. 

Further improvements were identified but could not be 

incorporated into the Phase 1 report. These improvements 

included (1) updating the model’s 7.5-year input dataset of 

external nutrient loading and lake monitoring data (Kann 

and Walker, 1999) with a 19-year dataset (Walker and others, 

2012), and (2) researching and developing algorithms to 

predict mid-season bloom declines caused by a transition 

from a primary-production-dominated to a communityrespiration-dominated (where community respiration includes 

cyanobacterial respiration and mortality) system. The 

inability to simulate mid-season bloom declines was a large 

source of error in the original model formulation (Wood and 

others, 2013). It also was suggested in the Phase 1 report that 

settling velocities could be measured rather than calibrated. 
4 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

ortac12-0723_fig01

River

Hawks

Lake

Upper

Klamath

Marsh

Lost River

Fourmile Canal

Sevenmile Canal

Modoc Point Main Canal

Swan

Lake

Round

Lake

Spring

Lake

Lake

Ewauna

Whiteline

Reservoir

Nuss

Lake

Wood

River

Marsh

John C Boyle

Reservoir

Meadow

Lake

Ball

Bay

Shoalwater Bay

Williamson

Agency Straits

Upper

Tulana

Goose

Bay

South

Marsh

Wocus

Marsh

Caledonia

Marsh

Hank's 

Marsh

Link

River

Moss

Creek

Pelican

Bay

Klamath

Lake

Aspen

Lake

Howard

Bay

Agency

Lake

Eagle

Sprague

Nannie Cr

Rock Creek

Cherry Creek

Fish

Banks

River

Varney

Creek

Clover Creek

Spencer

Creek

Wood River

Trout Creek

Anderson Creek

Lost Creek

A-Canal

Wocus

Marsh

Caledonia

Marsh

POINT ROAD

MODOC

Rattlesnake

Point

Hagelstein

County Park

Modoc Rim

Ball

Point

97

97

Bare

Island

Buck

Island

140

140

70

62

Klamath 

Falls

Plum Valley

Ridge

Modoc

Point

Link

River 

Dam

Coon

Point

Base modified from USGS and other digital sources, various scales. 

Coordinate system is UTM, Zone 10N; North American Datum of 1927.

EXPLANATION

4,110

4,143.3

Elevation, in feet above 

 Bureau of Reclamation datum

FBS Continuous water-quality

monitoring site

WMR Continuous water-quality

monitoring and sampling

site

MDN Water-quality sampling site

MDL Meteorological station

MDT

EPT

MDL

RPT

NBI

WMR

SSR

GBE

WDW

MDN

FBS

SHL

HPK

HDB

OREGON

Study

Area

122°10' 122° 121°50' 121°40'

42°

30'

42°

20'

42°

10'

0 2 4 6 8 10 MILES

0 2 4 6 8 10 KILOMETERS

KLAMATH RIVER

Figure 1. Location of meteorological stations and continuous water-quality monitoring and sampling sites, Upper Klamath Lake, 

Oregon.
Review of Phase 1 Results 5

This report further discusses the meaning of the net settling 

velocity in a system dominated by buoyant cyanobacteria, as 

well as the implications of calibration. With the availability of 

the updated monitoring data from Walker and others (2012) 

and new research on lake metabolism, the CSTR model 

calibration could be updated and new model algorithms could 

be pursued.

Although the biomass in UKL is dominated by 

cyanobacteria, the original TMDL model and its appropriate 

parameters and variables were described by the term “algal” 

instead of “cyanobacterial.” The following terms from the 

model are used in the remainder of this report: algal submodel, 

algal growth, algal settling velocity, algal phosphorus, and 

nonalgal phosphorus.

Purpose and Scope

This report presents the results of Phase 2 of the 

review of the 2001 TMDL model (Oregon Department of 

Environmental Quality, 2002) with the objective of evaluating 

modifications to the CSTR that include updated calibration 

data and alternate algorithms for important processes 

affecting water quality in the lake. The revised model is 

intended to help resource managers better (1) assess the 

efficacy of management strategies and restoration activities, 

and (2) understand topical areas needing additional data or 

research to reduce uncertainty in decision making.

Specifically, the report describes additional development 

of the R-script model from Phase 1 (Wood and others, 2013) 

and updated calibrations using newly available monitoring 

data (Walker and others, 2012), and derives a new, alternative 

algorithm to better account for removal of algal phosphorus 

from the system as the AFA bloom transitions from a 

growth-dominated to respiration-dominated regime. We 

describe the revised model calibration and verification, and 

we suggest five steps toward a new modeling framework that 

would preserve the computational speed and straightforward 

approach of the CSTR model. We anticipate that a model 

that incorporated these steps would prove to be an improved 

model, based on statistics describing its ability to simulate 

the calibration dataset. If so, such a model would provide 

more accurate predictions of water-column concentrations 

far into the future, while simulation of long-term response to 

management scenarios, such as TMDL implementation, would 

be improved. 

This report builds on Phase 1, and provides basic 

background on that effort as it pertains to specific 

modifications for the purposes of Phase 2; however, the reader 

is encouraged to consult the Phase 1 report together with this 

report in order to fully explain the history and rationale for the 

initial review, and for specific information about certain model 

parameter values or outcomes from the review. 

Review of Phase 1 Results 

In Phase 1, the verified R-script was used to make three 

changes that would improve representation of lake processes 

in the model as understood by established lake modeling 

techniques (Bowie and others, 1985; Chapra, 1997): 

1. Improve the light limitation factor by integrating the 

limitation function over the lake depth as opposed to 

calculating limitation from depth-integrated light;

2. Use Michaelis-Menten kinetics to describe phosphorus 

limitation and calibrate the half-saturation constant in 

tandem with other model parameters; and

3. Remove the assumption that chlorophyll a reaches 

equilibrium with light and nutrients over the time step. 

Following these changes, we incorporated new data collected 

in UKL since the TMDL model was originally developed, 

notably sediment phosphorus concentrations from multiple 

locations and sediment depths around the lake (Simon and 

others, 2009; Simon and Ingle, 2011), and recalibrated the 

model parameters. Once the model was recalibrated, we ran 

long-term simulations and conducted sensitivity analyses. 

These simulations allowed us to evaluate the lake’s predicted 

responses to changes in external phosphorus loading and to 

identify model parameters or areas of uncertainty that would 

benefit from additional study. Results of these analyses 

were detailed by Wood and others (2013), with additional 

suggestions for future improvement. Here, we summarize 

several of the notable Phase 1 model changes and results, 

which form the basis for additional proposed modifications 

in Phase 2.

An important aspect of the TMDL model is that 

it fundamentally relies on a coupled system whereby 

compartments representing the water column and lake 

sediments are interdependent. As stated by Wood and others 

(2013, p. 2): 

... 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 steadystate conditions.

A ramification of this model structure is that changes in 

one compartment are reflected in the other compartment at 

each biweekly time step. Thus, as calibration parameters 

change in Phases 1 and 2 resulting from incorporation of 

new datasets, there are implications for internal phosphorus 

loading, algal growth and maximum biomass, and pH. 
6 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

These implications are described for Phase 1 by Wood and 

others (2013) and for Phase 2 in section, “Phase 2 Model 

Updates” of this report, as appropriate. Nonetheless, the 

Phase 2 approach does not change the basic assumption of a 

coupled system between nutrients in the water column and 

bed sediments.

Changes to Algal Submodel 

In the TMDL model, an algal submodel uses a mass 

balance to determine biomass, with calculated lake-averaged 

terms for growth (GROWTH), respiration (RESP), and 

settling (SETTLE), as well as outflow concentrations from 

the lake (see equations 11–13 in Wood and others, 2013). 

GROWTH is a function of several factors, including a 

temperature factor, a phosphorus limiting factor (FP), and 

a light limiting factor (FL). In Phase 1, we modified the 

algorithms for FP and FL prior to recalibrating the model. 

The original formulation of FL represented light limitation 

over a 24-hour period by using depth-averaged light intensity 

in the Steele formulation and multiplying by photoperiod. 

We revised the model’s calculation of FL by integrating the 

Steele formulation over depth and time to better approximate 

average light attenuation in the water column for the entire 

lake (see equation 26 in Wood and others, 2013). FP originally 

was calculated as the percentage of nonalgal phosphorus (the 

difference between TP and calculated algal phosphorus) in 

the water column, which had the effect of underpredicting 

phosphorus limitation when nonalgal phosphorus was a 

large fraction of the TP, even if overall concentrations 

were low, and overpredicting phosphorus limitation when 

nonalgal phosphorus was a small fraction of the TP, even 

if large concentrations of phosphorus were available. We 

reformulated FP with the Michaelis-Menten kinetics approach 

for phosphorus limitation, which allows the algal growth rate 

to approach zero as the concentration of available (nonalgal) 

phosphorus approaches zero, and to approach a maximum 

value asymptotically as the concentration of the available 

phosphorus becomes large. 

After the light and phosphorus limitation factors were 

changed, the algal settling velocity, ualg, was recalibrated and 

a new parameter necessary for Michaelis-Menten kinetics, the 

half-saturation constant khalf, was introduced and calibrated 

in tandem with the recycling parameters (Wood and others, 

2013). The resulting calibrated algal setting velocity was 

41 percent lower, and the calibrated recycling rate was 

32 percent lower. The combined effect of the modifications 

to the phosphorus and light limitation factors and subsequent 

recalibration was to decrease the annual peak in predicted 

TP concentration relative to the original model by an average 

of 12 ppb. This was largely due to the lower recalibrated 

recycling rate; even so, the combination of lower algal settling 

velocity and modification to growth limitation factors resulted 

in higher peaks in chlorophyll a and pH. 

In addition to modifying the calculation of the limitation 

factors, the structure of the algal submodel was revised so that 

the assumption of maximum possible algal growth during each 

time step was removed. This resulted in chlorophyll a being 

represented as a first-order differential equation and modeled 

using the Runge-Kutta 4th-order method, as was used for the 

sediment and water-column compartment models with model 

terms and parameters defined in tables 1 and 2, respectively:

dB

dt

GROWTH RESP SETTLE Q

V

B out = − − − 



 



× (1)

The algal settling velocity was not recalibrated for this 

model, but recycling parameters were recalibrated. This 

change to the algal submodel structure resulted in a better 

model performance for the three measured variables (TP, 

chlorophyll a, pH) during the limited validation period of 

May–July and during the full validation period. 

The changes in calibration parameters among all versions 

of the model were substantial, but the performance of the 

model based on statistics calculated over the full verification 

dataset was not much changed, and was uniformly poor 

(Wood and others, 2013, tables 9–14). Furthermore, no model 

version consistently produced the best or worst results among 

the three steady-state variables—TP, 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. 

New Information—Internal Loading and 

Benthic Invertebrates

The Phase 1 review benefited from new data and insights 

into nutrient processes in the lake that were not available 

until after the implementation of the original TMDL in 2002. 

These data included revised estimates of sediment phosphorus 

concentrations, and field and laboratory measurements of the 

large densities of benthic invertebrates (worms and leeches) 

in the sediments, which previously had been unaccounted 

for and which appear to have important effects on internal 

nutrient recycling (Gardner and others, 1981; Hansen and 

others, 1998).
Review of Phase 1 Results 7

Table 1. Descriptions of model terms used in this report.

Model term Units Description

A square meters Lake area

B parts per billion Chlorophyll a concentration 

BNPP parts per billion Calculated net chlorophyll a concentration

C milligram per liter Dissolved oxygen concentration

Csat milligram per liter Saturation dissolved oxygen concentration

CR grams of oxygen per square meter per day Simulated community respiration

CRS grams of oxygen per square meter per day Seasonal component of community respiration model

CRW grams of oxygen per square meter per day Weekly component of community respiration model

D meter Lake depth

FT – Temperature adjustment factor

FP – Phosphorus limitation factor

FL – Light limitation factor

Io microeinsteins per square meter per second Light intensity at lake surface

I microeinsteins per square meter per second Site specific light intensity 

Φphoto grams of oxygen per square meter per hour Photosynthetic flux

Φresp grams of oxygen per square meter per hour Community respiration flux

Φresp,20 grams of oxygen per square meter per hour Community respiration flux at 20 degrees Celsius

ΦSOD grams of oxygen per square meter per hour Sediment oxygen demand

Φreaer grams of oxygen per square meter per hour Reaeration 

KL,20 meter per day Film transfer coefficient at 20 degrees Celsius

KR per day Maximum phosphorus recycle rate, pH-dependent

KG per day Maximum algal growth rate 

KT per day Maximum phosphorus recycle rate, temperature-dependent

Pna parts per billion Nonalgal phosphorus concentration

P parts per billion Site specific total phosphorus concentration

GPP grams of oxygen per square meter per day Simulated photosynthetic production

GPPS grams of oxygen per square meter per day Seasonal component of photosynthetic production model

GPPW grams of oxygen per square meter per day Weekly component of photosynthetic production model

GROWTH per day The algal growth component of the chlorophyll a submodel

RESP per day The respiration component of the chlorophyll a submodel

SETTLE per day The settling component of the chlorophyll a submodel

SET meters per day Space and time-dependent algal settling velocity

Qout cubic meters per day Lake outflow rate

s square meters 2-dimensional geographic location

Tlake degrees Celsius Average lake temperature 

T degrees Celsius Site specific water temperature

Tmin degrees Celsius Minimum temperature for algal growth to occur

T(min,rec) degrees Celsius Minimum temperature at which temperature-dependent recycling begins

t day Time

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

z meter Average lake depth over the time step

zfull meter Site specific full pool depth

Δt day Time step

WSPD10 meter per second Wind speed 10 meters above the water surface

WTDEV – Daily deviation in water temperature from the seasonal trend 

SSDEV – Daily deviation in shortwave solar radiation from clear sky solar radiation

WTROC degrees Celsius per day Daily rate of change in air temperature

STRAT degrees Celsius per day Degree of stratification

Appendix variables

ur meter per day Cyanobacterial rising velocity

uf meter per day Cyanobacterial falling velocity

rrise – Proportion of cyanobacterial colony rising

rfall – Proportion of cyanobacterial colony falling
8 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Table 2. Descriptions, values, and sources of parameters used in the phosphorus total maximum daily load (TMDL) model, Upper 

Klamath Lake, Oregon.

[Abbreviation: NA, not applicable]

Parameter

symbol

Parameter

description Value Units Sources of parameters

Revised chlorophyll a model recalibration

u Net settling velocity for algal phosphorus alg

pH* pH at which the total phosphorus recycle 

rate is one half the maximum

KR Maximum total phosphorus recycle rate

k Half-saturation constant for phosphorus half

limitation of algal growth

0.031

8.1

1.31

0.06

meter per day Calibrated using measured total phosphorus and 

chlorophyll a data

NA Calibrated using measured total phosphorus, 

chlorophyll a and pH data

per year Calibrated using measured total phosphorus, 

chlorophyll a and pH data

parts per million Calibrated within range provided by Bowie and 

others, 1985

Extended dataset recalibration

K Gross removal rate for nonalgal 0.016 na

phosphorus

u Net settling velocity for algal phosphorus 0.07–0.087 alg

pH* pH at which the total phosphorus recycle 8.2–8.3

rate is one half the maximum

KR Maximum total phosphorus recycle rate 1.61–2.09

k Half-saturation constant for phosphorus 0.05–0.06 half

limitation of algal growth

per day Calibrated using measured mass balance and total 

phosphorus data

meter per day Calibrated using measured total phosphorus and 

chlorophyll a data

NA Calibrated using measured total phosphorus, 

chlorophyll a and pH data

per year Calibrated using measured total phosphorus, 

chlorophyll a and pH data

parts per million Calibrated within range provided by Bowie and 

others, 1985

Calculated NPP model input recalibration

u Net settling velocity for algal phosphorus alg

pH* pH at which the total phosphorus recycle 

rate is one half the maximum

KR Maximum total phosphorus recycle rate

0.463

8.9–9

6.77–7.63

meter per day

NA

per year

Calibrated using measured total phosphorus and 

chlorophyll a data

Calibrated using measured total phosphorus, 

chlorophyll a and pH data

Calibrated using measured total phosphorus, 

chlorophyll a and pH data

Initial development of the TMDL model prior to 2002 

relied on limited data (Eilers and others, 2001) for phosphorus 

concentrations in lake sediment; however, 26 sediment cores 

subsequently were collected and analyzed during 2005–06 at 

discrete sediment depths (Simon and others, 2009; Simon and 

Ingle, 2011). The average sediment phosphorus concentration 

calculated from the 26 newer cores was about one-third 

of the concentration used in the initial model calibration; 

therefore, we performed a sensitivity analysis to determine 

the significance of the updated sediment data on model output 

and calibration parameters. The overall effect of the revised 

sediment phosphorus concentration was to increase the 

pH-dependent phosphorus recycling rate (KR) of the model.

Given the coupled water-column–lake sediment compartments 

in the model, the reduced sediment phosphorus concentration 

caused an increase in KR, resulting in an overall increase in 

phosphorus recycling in the lake, regardless of pH. 

In Phase 1, the potential contribution of benthic 

invertebrates to the lake’s nutrient cycling was investigated 

(Wood and others, 2013). Based on sampling in 2009 and 

2010, average invertebrate densities in the lake’s sediments 

were high, on the order of 20,000 individuals/m2

 (standard 

deviation =11,671). Invertebrate taxa primarily were 

Oligochaeta (worms, 65 percent), Hirudinea (leeches, 

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

Invertebrate bioturbation can significantly increase the 

diffusion and advection of solutes, including nutrients 

(Michaud and others, 2005, 2006; Mermillod-Blondin and 

Rosenberg, 2006; Mermillod-Blondin, 2011), and invertebrate 

metabolic excretion can be a significant flux from the 

sediments into the water column (Fukuhara and Yasuda, 1985; 

Fukuhara and Sakamoto, 1987; Devine and Vanni, 2002). 

In the Phase 1 report, experiments with two types of 

invertebrate common in UKL produced estimates of metabolic 

flux of 3.8 [(mg P/m2

)/d]. There is considerable spatial and 

temporal variation in the actual excretion rates, according 

to differences in biotic and abiotic factors around the lake, 

such as sediment characteristics, invertebrate taxa and 

densities, and water temperatures (which affect growth rates). 

Such variability cannot be resolved with a CSTR model. 
Review of Phase 1 Results 9

Incorporation of invertebrates into the 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. 

Scaling up from small scales to the bulk system is discussed 

further in section, “Lake Metabolism as a Substitute for Algal 

Growth and Respiration.”

Nutrient sources attributable to benthic invertebrates 

were not differentiated in the original model, but instead 

were incorporated into a lumped, residual-derived internal 

loading term. During Phase 1, the implications of bioturbation 

and excretion in the model were investigated by substituting 

the original, pH-dependent recycling mechanism with a 

temperature-dependent recycling mechanism that might 

apply to invertebrate-caused mechanisms and performing 

a sensitivity analysis. This did not significantly improve 

model performance, which indicates that neither the 

pH-dependent nor the temperature-dependent mechanism 

is uniquely suited for describing the measurements; rather, 

any rate-dependent first-order term that generally increased 

through the spring and summer (such as pH and temperature) 

would perform approximately equally well with the current 

CSTR model structure, and any such term probably would 

have poor performance statistics, particularly over the full 

validation period. 

Research is still needed to refine lake-average benthic 

invertebrate excretion rates for future use in the model; 

to resolve interactions between sediment characteristics, 

invertebrate ecology, and water-column nutrient dynamics; 

and to support a more spatially explicit lake water-quality 

model if one is developed in the future. However, the model is 

not currently configured to use both the temperature‐dependent 

and pH-dependent internal loading mechanisms. Such work 

would require separate calibrations to establish recycle rates 

that are independent from each other, which is currently 

not possible with the available data. Future laboratory work 

could provide the data necessary to include both of these 

terms if deemed useful. Incorporating benthic invertebrate 

recycling as a separate process in the model would provide 

opportunities to explore management options that could 

affect invertebrate ecology and to differentiate those effects 

from other physicochemical effects, such as high-pH induced 

phosphorus release.

203-Year Simulation and Sensitivity Analysis

Using the revised CSTR model, with variations 

corresponding to 0- and 40-percent reductions in external 

loading according to the TMDL (Oregon Department of 

Environmental Quality, 2002), Wood and others (2013) 

performed a 203-year simulation, chaining together the 

original 7-year, April 1991 to April 1998, nutrient budget 

for the lake (Kann and Walker, 1999) as input data. Model 

versions that were simulated included the original model 

and versions with and without the modifications for light and 

phosphorus limitation, temperature-dependent phosphorus 

recycling (for benthic invertebrates), and the revised solution 

(Runga-Kutta) for chlorophyll a whereby equilibrium between 

cyanobacterial biomass and light and nutrients is not necessarily 

achieved within a given 2-week time step. 

All versions of the model, when used to simulate watercolumn concentrations characterizing a steady state 203 years 

into the future, produced lower TP 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 higher 

recycling rates 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. 

The phosphorus stored in the sediment reservoir 

generally—although not always (Wood and others, 2013, 

versions O and B′)—decreased over time, and as with the 

water-column concentrations, there was a wide range in the 

steady-state value achieved. These ultimate values were related 

to the ratios of recycled phosphorus mass to sedimented 

phosphorus mass (Wood and others, 2013, fig. 13) and were 

therefore related to the values of the calibration parameters. 

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 

steady-state sediment reservoir—a lower steady-state sediment 

reservoir was associated with higher steady-state water-column 

phosphorus concentration.

The long-term simulations indicate that the sediment 

reservoir influences the long-term behavior of the modeled 

water-column concentrations in an important and possibly 

counter-intuitive way. The 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 steady-state 

conditions in the water column and sediments are determined 

by the calibration parameters. In some simulations, the size 

of the sediment phosphorus reservoir remained 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 (Wood and others, 2013).
10 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Finally, the importance of uncertainty in the values of the 

algal settling and phosphorus recycling calibration parameters 

was most apparent for long-term model simulations. 

Although short-term simulations (up to several years) using 

the different model versions produced similar results even 

though the values of these calibration parameters varied 

by as much as 40 percent, using different calibrated algal 

settling and phosphorus recycling rates resulted in relatively 

large differences in water-column concentrations when 

simulations exceeded several decades. The same was true for 

simulated sediment phosphorus concentrations. These findings 

indicate that:

1. The initial size of the sediment reservoir is not as 

important in determining steady-state conditions as 

the rates of algal and nonalgal settling and sediment 

phosphorus recycling;

2. These calibration parameters have a lag effect on 

nutrient accumulation and dynamics in the lake, where 

algal and nonalgal settling and sediment recycling in one 

year affects the water quality in the following year; and

3. The efficacy of external loading reductions is highly 

dependent on these calibration parameters, which are 

poorly understood (Wood and others, 2013).

It follows that efforts to bound the rate estimates of algal and 

nonalgal settling velocities and sediment phosphorus recycling 

with targeted field-data collection would be a high priority for 

reducing uncertainty in future lake management decisions.

Phase 2 Model Updates

Wood and others (2013) made several suggestions 

for updating the model. In particular, reducing uncertainty 

in the recycling and sedimentation terms is an important 

future research need. This can be achieved through (1) field 

measurement of specific parameters needed for the relevant 

algorithms in the model, or (2) development of new 

algorithms that would predict these terms from independent 

environmental data. The calibrated settling velocity in 

particular remains a challenge to interpret; measuring this 

parameter directly rather than determining it by calibration of 

the model (and the difficulty of doing so in a lake ecosystem 

dominated by buoyant cyanobacteria) is discussed further 

in section, “Future Model Needs.” Field work related to 

recycling by benthic invertebrates is ongoing and is expected 

to inform future modeling efforts. 

Wood and others (2013) also noted that model 

performance was limited by the inability of the model to 

simulate mid-season bloom declines, and suggested that 

the addition of a “mortality” term for the cyanobacteria 

would help capture the dynamics of bloom declines. In 

subsequent sections, we report on new equations to calculate 

cyanobacterial metabolism in the surface waters, using 

hourly measurements of DO available from the USGS lake 

monitoring program from 2006 through 2008 (Lindenberg 

and others, 2009; Kannarr and others, 2010). These estimates 

rely on independent data (water temperature, DO, and nutrient 

concentrations) to predict net primary production (NPP) in 

the lake, including cyanobacterial biomass and, importantly, 

the transition in NPP from dominance by primary production 

(autotrophy) to dominance by community respiration 

(heterotrophy), which includes aquatic respiration and 

mortality. This proof-of-concept approach is based on limited 

available data; additional data collection requirements to more 

completely develop this modeling framework are discussed in 

section, “Future Model Needs.”

Updated CSTR Model Calibration

One of the important purposes for performing the Phase 2 

evaluation of the CSTR model was inclusion of a longer 

dataset for model calibration because the original model was 

calibrated over the 1991–98 period. Since that time, data 

have been collected regularly and a 19-year set (1991–2010) 

is available for model evaluation (Walker and others, 

2012). Model terms and symbols are defined in table 1 and 

recalibrated parameters, determined through methods outlined 

below, are presented in table 2. 

Reevaluate Phase 1 Test Statistics Using 

Updated Algal Settling Velocity Calibration, 

1991–98

In Phase 1, the final modification to the model structure 

was developed for the algal submodel. The original TMDL 

model included the assumption that equilibrium between 

growth and respiration was achieved during each time step 

such that net chlorophyll a concentration equaled zero. The 

updated model (versions D and D') did not include this 

assumption and solved the first-order differential equation for 

chlorophyll a using the fourth-order Runge-Kutta numerical 

solution technique. In Phase 1, the recycling parameters were 

recalibrated for this model (version D), but algal settling 

velocity was not, and its value was assumed to be equal to 

that of the chlorophyll a equilibrium model (version B). In 

this review, we evaluated the sensitivity of algal settling 

velocity for a non-equilibrium chlorophyll a model, herein 

referred to as model version E', and the parameter calibration 

steps follow. 

Step 1—Algal Settling Velocity and Half-Saturation 

Constant 

Algal settling velocity and the half-saturation constant 

(used in the Michaelis-Menten phosphorus limitation 

algorithm) were reevaluated in tandem over the original 

7.5-year dataset following the approach outlined by Kann and 
Phase 2 Model Updates 11

Walker (1999) and used by Wood and others (2013). Algal 

settling velocity values were assumed to be between 0.001 and 

0.2 m/d, and half-saturation constant values were assumed 

to be between 0.01 and 0.06 ppm (Bowie and others, 1985); 

the best combination of parameter values was determined 

by correlating chlorophyll a as calculated by equation 1 

(using measured phosphorus concentrations) with measured 

chlorophyll a concentrations in June–July for the 1994–98 

subset of values. The calibration process minimized the sum 

of squared errors with the combination of ualg=0.031 m/d and 

khalf =0.06 ppm. These parameter values were similar to values 

previously determined in Phase 1 of the evaluation. 

In the process of recalibrating algal settling velocity, 

we also evaluated the differences in parameters when 

other subsets of values were used for calibration besides 

the previously used 1994–98 subset (Walker, 2001; Wood 

and others, 2013). We evaluated the subset from June to 

July over 1991–93 and determined that the sum of squared 

errors was minimized with the combination of ualg=0.2 m/d 

and khalf =0.02 ppm, which is different than the previously 

determined set and which indicates an algal settling velocity 

more than six times greater than the 1994–98 subset and a 

half saturation constant one-third the value of the original 

subset. Further inspection of these results showed that 1 year 

of data, 1992, was primarily responsible for the large changes 

in calibration parameters, probably due to its higher than 

average TP and chlorophyll a concentrations. When the 

settling parameters were recalibrated for June–July in all years 

excluding 1992, the sum of squared errors was minimized 

with the combination of ualg=0.048 m/d and khalf =0.06 ppm. 

This process highlights the potential sensitivity of the settling 

velocity to anomalous years and indicates the need for a 

longer calibration set or, preferably, actual measurement of 

the parameters. 

Step 2—Internal Recycling Parameters

The recycling parameters, pH* and KR, also were 

recalibrated. This calibration step evaluated the recycling 

equation in the CSTR model by running the model and 

minimizing the objective function, the product of the 

mean squared error of the three state variables (Wood and 

others, 2013, equation 24), for May–July data in 1994–98. 

The resulting recalibrated parameters were pH*=8.1 and 

KR=1.31/yr, which are similar to those parameters used 

previously for model version D'. 

Reevaluation of Phase I Test Statistics Using 

Extended Data, 1991–2010

The 1991–2010 model input dataset was developed 

using daily nutrient and hydrologic loading data (Walker and 

others, 2012), Tule Lake average daily air temperature and 

solar radiation data from California Irrigation Management 

Information System site TULELK2.A (University of 

California, 2014) and Klamath Tribes biweekly/monthly 

(depending on season) water temperature data (Kann, 2011) 

from the nine long-term sites. Following Kann and Walker 

(1999), a multiple regression equation relating average water 

temperature to the air temperature and solar radiation collected 

25 mi away at Tule Lake in northern California was used to 

determine a biweekly average water temperature value for 

UKL. The biweekly nutrient and hydrologic loading set was 

developed by accumulating the daily masses, or volumes, 

across the biweekly intervals used in the model.

A biweekly dataset of water quality in UKL for the 

longer calibration period was developed using Klamath 

Tribes measured biweekly/monthly (depending on season) 

TP, chlorophyll a, and pH data. The data were averaged 

across space and time to represent the lakewide values for 

the appropriate biweekly intervals, following the procedure 

outlined in Kann and Walker (1999). Model parameters 

were calibrated with the 1991–2010 dataset following the 

steps outlined by Kann and Walker (1999) for the version 

herein referred to as F ′, which used the updated sediment 

reservoir and growth limitation factors, and the Runge-Kutta 

algal submodel, corresponding to characteristics of model 

version E'. Recalibrated parameter values are presented in 

table 2.

Step 1—pH Regression Model Parameters

As reported previously (Wood and others, 2013, 

equation 20), we performed the initial multiple regression 

relating pH to the natural logarithm of chlorophyll a and 

Julian day (with a minimum value of 200) in Microsoft Excel®

using biweekly, lakewide average chlorophyll a concentrations 

and pH, calculated from biweekly measurements. The 

regression was evaluated during June–August 1991–2010 

and included 123 data points. The model coefficients and 

performance statistics (coefficient of determination, R2, 

and standard error) compared well to those coefficients and 

statistics originally determined over the 1991–98 period, 

with R2 slightly decreased (from 0.80 to 0.74) and standard 

error slightly increased (from 0.24 to 0.28). We developed a 

revised pH submodel regression that removed Julian day as an 

independent variable. This was done to prevent discontinuities 

in pH when Julian day changed from 365 to 1, an occurrence 

that was noted in the Phase 1 report (Wood and others, 2013). 

This new regression model was evaluated for measurements 

from all months in 1991–2010, as opposed to the shorter, 

July–August period, which included 242 data points. Again, 

the model coefficients and performance statistics compared 

well to the coefficients and statistics used in the previous 

approach, with R2

 slightly decreased from the original model 

(from 0.80 to 0.70) and standard error increased (from 0.24 to 

0.37), with the reduction in model performance due partly to 

the change in evaluation dataset. The pH model regression 

coefficients and statistics are presented in table 3.
12 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Table 3. pH model coefficients and statistics between the original total maximum daily load (TMDL) model for Upper Klamath Lake, 

Oregon, and the recalibration considering models with and without a Julian day component.

[Abbreviations: NA, not applicable; <, less than]

Quantity

Original

TMDL 

value

Recalibration value, 

calibration months July–August

Recalibration value no Julian,

calibration months January–December

1991–1998 p–value 1991–2010 p–value 1991–1998 p–value 1991–2010 p–value

pH model: intercept 7.931 7.829 <2e-16 7.729 <2e-16 7.025 <2e–16 7.161 <2e–16

pH model: chlorophyll a

coefficient

0.5338 0.5514 <2e-16 0.5271 <2e-16 0.430 <2e–16 0.4211 <2e–16

pH model: Julian day 

coefficient

–0.006 –0.006 0.0169586 –0.005 0.0110999 NA NA NA NA

Julian threshold 200 200 – 200 – NA – NA –

Coefficient of 

determination, R2

0.80 0.78 – 0.74 – 0.73 – 0.71 –

Standard error 0.24 0.25 – 0.28 – 0.37 – 0.37 –

Sample size, N – 52 – 123 – 113 – 242 –

Step 2—Algal Settling Velocity and Half-Saturation 

Constant

Algal settling velocity and the half-saturation constant 

were evaluated in tandem over the extended dataset following 

the approach by Kann and Walker (1999) and Wood and others 

(2013). Algal settling velocity values were assumed to be 

between 0.001 and 0.2 m/d and half-saturation constant values 

were assumed to be between 0.01 and 0.06 ppm (Bowie and 

others, 1985); the best value was determined by correlating 

chlorophyll a as calculated by equation 1 (using measured 

phosphorus concentrations) with measured chlorophyll a

concentrations in June–July over a dataset consisting of only 

even years, and a second dataset consisting of only odd years, 

referred to as model versions F ′

E and F ′O, respectively. For 

calibration performed using odd years, the sum of squared 

errors was minimized by the combination of ualg=0.070 m/d 

and khalf =0.06 ppm, while the calibration performed over 

even years minimized the sum of squared errors with the 

combination of ualg=0.087 m/d and khalf =0.05 ppm, indicating 

similar parameter results between calibration sets.

Step 3—Nonalgal Phosphorus Loss Rate 

The nonalgal phosphorus loss rate was evaluated by 

applying the method outlined by Wood and others (2013) to 

the longer 1991–2010 dataset. The March–May datasets were 

isolated, and negligible recycling and predominantly nonalgal 

phosphorus were assumed. A trial-by-error and least-squares 

fit process was used for calibration, with Kna=0.016/d 

determined to minimize the sum of squared errors. This value 

was only slightly larger than the value of 0.015/d determined 

in the Phase I report and the value of 0.012/d presented in 

Kann and Walker (1999). The correlation coefficient (R), 

calculated when comparing “measured” sedimented mass (as 

determined from a whole-lake mass balance calculation) to 

the simulated mass, was still low (R=0.28), but had increased 

from the correlation coefficient calculated for the shorter 

1991–98 time period (R=0.19). The regression coefficient was 

not highly significant (p=0.06). The calibrated value of Kna

resulted in a lower average value and lower standard deviation 

as compared to the “measured” set, which was consistent with 

the results from Phase 1. The comparisons between measured 

and calculated nonalgal sedimentation are presented in table 4. 

Step 4—Internal Recycling Parameters

The final calibration step was to evaluate the recycling 

equation in the CSTR model, which depended on two 

parameters, pH* and KR, by running the model and minimizing 

the objective function (Wood and others, 2013, equation 24). 

For this analysis, summation is over all May–July data in 

both odd and even years, versions F ′

O and F ′E, respectively, 

and two sets of recycling parameters were compared. For 

calibration performed using odd years, the objective function 

was minimized with the combination of pH*=8.2 and 

KR=1.61/year and for calibration performed using even years, 

the objective function was minimized with the combination of 

pH*=8.3 and KR=2.09/year. 

Two observations help with the interpretation of the 

calibration results of the various model versions. First, when 

the assumption of equilibrium was made in the algal submodel 

(versions O', A', B' and C'), the resulting pH* values were 

greater than 8.6; however, when a Runge-Kutta technique was 

used to estimate for chlorophyll a (versions D', E' and F'), 

the pH* values always were less than or equal to 8.3. This 

has important implications for the steady-state response of 
Phase 2 Model Updates 13

Table 4. Descriptive statistics between measured and calculated nonalgal phosphorus mass lost to 

sedimentation during March–May, Upper Klamath Lake, Oregon, 1991–2010.

[Description of model parameters are shown in table 2. Abbreviations: kg, kilogram; NA, not applicable]

Mass lost to sedimentation

1991–1998 1991–2010

Average

(kg)

Standard 

deviation 

(kg)

Average

(kg)

Standard 

deviation 

(kg)

Calculated from measured (mass balance) 9,008 6,608 9,415 6,001

Calculated with Kna= 0.012 d-1 (original model value) 6,781 1,813 6,800 1,567

Calculated with Kna= 0.015 d-1 (recalibrated model value) 8,476 2,267 8,499 1,959

Calculated with Kna= 0.016 d-1 (recalibrated model value) NA NA 9,066 2,089

the model (discussed in section, “Model Performance”), as 

the lower pH* threshold results in a more even distribution of 

bloom magnitude over many years, whereas the higher pH*

threshold results in a more bimodal response, in which blooms 

are large in some years and nearly nonexistent in others. 

Second, for a given initial sediment reservoir and recycling 

mechanism (but regardless of modifications to the limitation 

factors or algal submodel assumptions), a positive correlation 

was made between ualg and KR (ρ=0.71, p=0.14); for example, 

compare versions O', A', B', E', F'O and F'E. Therefore, any 

modification that increases the recycle rate also increases the 

algal settling velocity, and vice versa.

Model Performance

Recalibrating the model parameters for the updated 

chlorophyll a solution method (version E') resulted in little 

change or improvement in model performance statistics, partly 

because the parameters changed little. When the full validation 

period was considered (tables 5–7), the statistics were similar, 

although generally lower in the newly recalibrated model 

version E', for all variables. When the May–July months 

were isolated and evaluated for performance, model version 

E' performed similarly to model version D', but had slightly 

lower correlation coefficient, coefficient of determination, 

and Nash-Sutcliffe (NS) values and slightly greater root 

mean square error (RMSE) (tables 8–10). Overall, the model 

performance was not substantially changed by the parameter 

recalibration over the 1991–98 period even though the 

parameters changed. 

Recalibrating the model parameters with the extended 

dataset (versions F'O and F'E) resulted in little change or 

improvement in model performance statistics when compared 

to model version E'. The performance statistics were similar 

for all variables when respective full validation periods 

were considered (tables 5–7), although model version E'

tended to have slightly better performance for all variables. 

The performance statistics for the recalibrated models were 

best for pH, with correlation coefficients greater than 0.60 

and NS values greater than 0.34; however, statistics for TP 

and chlorophyll a did not necessarily indicate good model 

performance, with correlation coefficients less than 0.45 

and NS values less than 0.05. When May–July months were 

isolated for validation (tables 8–10), the statistics indicated 

better performance for all variables than when the full years 

were considered. However, the May–July statistics for the 

model calibrated and validated over the extended period 

(versions F'O and F'E) were not very different than statistics 

for the model calibrated and validated over the original 

7.5-year period (version E'). The performance statistics 

were best in both model versions for TP and pH and worst 

for chlorophyll a. For TP, R2 values were 0.61–0.64, and 

NS values were 0.53–0.59; for pH, R2 values were 0.60–0.66, 

and NS values were 0.48–0.59; and for chlorophyll a, R2

values were 0.39–0.45, and NS values were 0.27–0.32.

In figure 2, the simulated TP, chlorophyll a, and pH 

values were compared to measured values for 1991–2010 for 

model versions E', F'O, and F'E. The differences in the results 

of the three models appear small relative to the differences 

between measured and simulated values. All simulations 

represent a smoothed fit of the measurements, and all do a 

generally poor job of capturing the extremes in the data. In 

the graphs of TP and chlorophyll a, the downward trend is 

apparent in the simulated values that approach steady state 

(<200 ppb for TP and <150 ppb for chlorophyll a) within 

the 19-year data period. This downward trend was due to the 

combination of smaller sediment reservoir, increased recycling 

rates, and decreased pH* values that allow sediment recycling 

to begin at lower lake-average pH values and last for a longer 

time period. These three components cause the recalibrated 

model to approach steady state relatively quickly.
14 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Table 5. Performance statistics for total phosphorus obtained over the indicated validation periods for model versions discussed in the 

report.

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below are repeated from Phase 1 (Wood and others, 2013, table 9)

Original model with original calibration

O EQ JULIAN INIT REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

O′ O EQ JULIAN INIT_NEW REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

Recalibration of original model

A O EQ JULIAN INIT REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

A′ A EQ JULIAN INIT_NEW REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

Modifications to light and phosphorus limitation

B A EQ JULIAN INIT REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

B′ B EQ JULIAN INIT_NEW REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

Temperature-dependent recycling

C B EQ JULIAN INIT REC_T P_MM L_DI All months 1991–1993 and August–April 1994–1998

C′ C EQ JULIAN INIT_NEW REC_T P_MM L_DI All months 1991–1993 and August–April 1994–1998

Runge-Kutta solution of chlorophyll a

D B RK JULIAN INIT REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

D′ D RK JULIAN INIT_NEW REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias

(ppb)

RMSE

(ppb) NS

The rows below are repeated from Phase I (Wood and others, 2013, table 9)—Continued

Original model with original calibration

O 0.05 0.53 9.1 NA NA NA 0.56 0.32 83 -48.6 109.6 -0.58

O′ O 0.05 1.33 9.0 NA NA NA 0.60 0.36 83 -23.3 83.7 0.08

Recalibration of original model

A O 0.059 0.47 9.0 NA NA NA 0.59 0.35 83 -24.7 86.9 0.01

A′ A 0.059 1.24 8.9 NA NA NA 0.61 0.38 83 -10.4 75.1 0.26

Modifications to light and phosphorus limitation

B A 0.035 0.32 8.6 NA NA 0.06 0.58 0.34 83 -23.3 82.8 0.10

B′ B 0.035 1.05 8.7 NA NA 0.06 0.60 0.35 83 -23.4 83.1 0.10

Temperature-dependent recycling

C B 0.035 NA NA 0.54 12.8 0.06 0.59 0.35 83 -31.8 90.7 -0.08

C′ C 0.035 NA NA 1.62 12.8 0.06 0.60 0.36 83 -25.1 84.8 0.06

Runge-Kutta solution of chlorophyll a

D B 0.035 0.34 8.1 NA NA 0.06 0.38 0.15 83 -60.0 114.9 -0.73

D′ D 0.035 1.18 8.1 NA NA 0.06 0.39 0.15 83 -52.6 109.9 -0.59
Phase 2 Model Updates 15

Table 5. Performance statistics for total phosphorus obtained over the indicated validation periods for model versions discussed in the 

report.—Continued

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below were generated for Phase 2

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months even years and August–April odd years

F′O G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months 2006–2008

F′E E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months odd years and August–April even years

F′E G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months 2006–2008

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA All months 2006

G′2 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA All months 2007

G′3 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA All months 2008

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias

(ppb)

RMSE

(ppb) NS

The rows below were generated for Phase 2—Continued

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ 0.031 1.31 8.1 NA NA 0.06 0.34 0.11 83 -59.7 119.6 -0.88

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ 0.070 1.61 8.2 NA NA 0.06 0.38 0.15 183 7.9 90.2 0.05

F′O G′ 0.070 1.61 8.2 NA NA 0.06 0.69 0.48 33 61.7 102.7 0.04

F′E E′ 0.087 2.09 8.3 NA NA 0.05 0.39 0.16 181 2.2 92.1 -0.03

F′E G′ 0.087 2.09 8.3 NA NA 0.05 0.72 0.52 33 54.9 93.9 0.20

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ 0.463 6.77 8.9 NA NA NA 0.91 0.82 9 2.3 37.9 0.82

G′2 F′ 0.463 7.63 9 NA NA NA 0.44 0.19 9 40.5 109.5 0.04

G′3 F′ 0.463 6.76 8.9 NA NA NA 0.84 0.71 9 1.8 43.2 0.71
16 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Table 6. Performance statistics for chlorophyll a obtained over the indicated validation periods for model versions discussed in the 

report.

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below are repeated from Phase I (Wood and others, 2013, table 10)

Original model with original calibration

O EQ JULIAN INIT REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

O' O EQ JULIAN INIT_NEW REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

Recalibration of original model

A O EQ JULIAN INIT REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

A′ A EQ JULIAN INIT_NEW REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

Modifications to light and phosphorus limitation

B A EQ JULIAN INIT REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

B′ B EQ JULIAN INIT_NEW REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

Temperature-dependent recycling

C B EQ JULIAN INIT REC_T P_MM L_DI All months 1991–1993 and August–April 1994–1998

C′ C EQ JULIAN INIT_NEW REC_T P_MM L_DI All months 1991–1993 and August–April 1994–1998

Runge-Kutta solution of chlorophyll a

D B RK JULIAN INIT REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

D′ D RK JULIAN INIT_NEW REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias

(ppb)

RMSE

(ppb) NS

The rows below are repeated from Phase I (Wood and others, 2013, table 10)—Continued

Original model with original calibration

O 0.05 0.53 9.1 NA NA NA 0.42 0.17 83 -56.5 129.6 -1.77

O′ O 0.05 1.33 9.0 NA NA NA 0.43 0.18 83 -45.7 115.0 -1.18

Recalibration of original model

A O 0.059 0.47 9.0 NA NA NA 0.42 0.18 83 -41.1 112.3 -1.08

A′ A 0.059 1.24 8.9 NA NA NA 0.43 0.18 83 -34.7 103.9 -0.78

Modifications to light and phosphorus limitation

B A 0.035 0.32 8.6 NA NA 0.06 0.40 0.16 83 -51.1 127.5 -1.68

B′ B 0.035 1.05 8.7 NA NA 0.06 0.40 0.16 83 -52.0 129.6 -1.77

Temperature-dependent recycling

C B 0.035 NA NA 0.54 12.8 0.06 0.36 0.13 83 -59.2 143.3 -2.39

C′ C 0.035 NA NA 1.62 12.8 0.06 0.37 0.14 83 -55.7 138.4 -2.16

Runge-Kutta solution of chlorophyll a

D B 0.035 0.34 8.1 NA NA 0.06 0.48 0.23 83 -61.0 114.4 -1.16

D′ D 0.035 1.18 8.1 NA NA 0.06 0.48 0.23 83 -57.2 112.3 -1.08
Phase 2 Model Updates 17

Table 6. Performance statistics for chlorophyll a obtained over the indicated validation periods for model versions discussed in the 

report.—Continued

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare 

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below were generated for Phase 2

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months even years and August–April odd years

F′O G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months 2006–2008

F′E E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months odd years and August–April even years

F′E G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months 2006–2008

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA All months 2006

G′2 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA All months 2007

G′3 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA All months 2008

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias

(ppb)

RMSE

(ppb) NS

The rows below were generated for Phase 2—Continued

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ 0.031 1.31 8.1 NA NA 0.06 0.46 0.21 83 -63.9 121.9 -1.45

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ 0.070 1.61 8.2 NA NA 0.06 0.45 0.20 183 0.183 79.6 -0.02

F′O G′ 0.070 1.61 8.2 NA NA 0.06 0.52 0.27 33 22.0 77.6 0.17

F′E E′ 0.087 2.09 8.3 NA NA 0.05 0.42 0.18 181 -3.1 81.9 -0.16

F′E G′ 0.087 2.09 8.3 NA NA 0.05 0.50 0.25 33 17.0 79.2 0.14

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ 0.463 6.77 8.9 NA NA NA 0.44 0.19 9 8.2 68.0 0.13

G′2 F′ 0.463 7.63 9 NA NA NA 0.66 0.43 9 26.7 60.9 0.18

G′3 F′ 0.463 6.76 8.9 NA NA NA 0.52 0.27 9 40.1 96.6 0.08
18 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Table 7. Performance statistics for pH obtained over the indicated validation periods for model versions discussed in the report.

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below are repeated from Phase I (Wood and others, 2013, table 11)

Original model with original calibration

O EQ JULIAN INIT REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

O' O EQ JULIAN INIT_NEW REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

Recalibration of original model

A O EQ JULIAN INIT REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

A′ A EQ JULIAN INIT_NEW REC_pH P_00 L_AVG All months 1991–1993 and August–April 1994–1998

Modifications to light and phosphorus limitation

B A EQ JULIAN INIT REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

B′ B EQ JULIAN INIT_NEW REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

Temperature-dependent recycling

C B EQ JULIAN INIT REC_T P_MM L_DI All months 1991–1993 and August–April 1994–1998

C′ C EQ JULIAN INIT_NEW REC_T P_MM L_DI All months 1991–1993 and August–April 1994–1998

Runge-Kutta solution of chlorophyll a

D B RK JULIAN INIT REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

D′ D RK JULIAN INIT_NEW REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias RMSE NS

The rows below are repeated from Phase I (Wood and others, 2013, table 11)—Continued

Original model with original calibration

O 0.05 0.53 9.1 NA NA NA 0.57 0.32 83 -0.195 0.716 -0.09

O′ O 0.05 1.33 9.0 NA NA NA 0.56 0.32 83 -0.175 0.703 -0.05

Recalibration of original model

A O 0.059 0.47 9.0 NA NA NA 0.56 0.31 83 -0.155 0.700 -0.04

A′ A 0.059 1.24 8.9 NA NA NA 0.55 0.31 83 -0.143 0.694 -0.02

Modifications to light and phosphorus limitation

B A 0.035 0.32 8.6 NA NA 0.06 0.57 0.33 83 -0.161 0.704 -0.05

B′ B 0.035 1.05 8.7 NA NA 0.06 0.58 0.33 83 -0.157 0.699 -0.04

Temperature-dependent recycling

C B 0.035 NA NA 0.54 12.8 0.06 0.57 0.32 83 -0.175 0.716 -0.09

C′ C 0.035 NA NA 1.62 12.8 0.06 0.57 0.32 83 -0.168 0.712 -0.07

Runge-Kutta solution of chlorophyll a

D B 0.035 0.34 8.1 NA NA 0.06 0.64 0.41 83 -0.247 0.622 0.18

D′ D 0.035 1.18 8.1 NA NA 0.06 0.65 0.42 83 -0.235 0.612 0.21
Phase 2 Model Updates 19

Table 7. Performance statistics for pH obtained over the indicated validation periods for model versions discussed in the report.—

Continued

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare 

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below were generated for Phase 2

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months 1991–1993 and August–April 1994–1998

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months even years and August–April odd years

F′O G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months 2006–2008

F′E E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months odd years and August–April even years

F′E G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI All months 2006–2008

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA All months 2006

G′2 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA All months 2007

G′3 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA All months 2008

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias RMSE NS

The rows below were generated for Phase 2—Continued

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ 0.031 1.31 8.1 NA NA 0.06 0.55 0.30 83 -0.327 0.694 -0.02

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ 0.070 1.61 8.2 NA NA 0.06 0.63 0.40 183 -0.055 0.534 0.38

F′O G′ 0.070 1.61 8.2 NA NA 0.06 0.73 0.53 33 0.183 0.502 0.46

F′E E′ 0.087 2.09 8.3 NA NA 0.05 0.60 0.37 181 -0.038 0.552 0.34

F′E G′ 0.087 2.09 8.3 NA NA 0.05 0.72 0.53 33 0.182 0.504 0.45

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ 0.463 6.77 8.9 NA NA NA 0.70 0.50 9 0.067 0.572 0.42

G′2 F′ 0.463 7.63 9 NA NA NA 0.76 0.58 9 0.431 0.548 -0.23

G′3 F′ 0.463 6.76 8.9 NA NA NA 0.77 0.59 9 0.167 0.366 0.48
20 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Table 8. Performance statistics for total phosphorus obtained over the indicated May–July validation periods for model versions 

discussed in the report.

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below are repeated from Phase I (Wood and others, 2013, table 12)

Original model with original calibration

O EQ JULIAN INIT REC_pH P_00 L_AVG May–July 1991–1993

O′ O EQ JULIAN INIT_NEW REC_pH P_00 L_AVG May–July 1991–1993

Recalibration of original model

A O EQ JULIAN INIT REC_pH P_00 L_AVG May–July 1991–1993

A′ A EQ JULIAN INIT_NEW REC_pH P_00 L_AVG May–July 1991–1993

Modifications to light and phosphorus limitation

B A EQ JULIAN INIT REC_pH P_MM L_DI May–July 1991–1993

B′ B EQ JULIAN INIT_NEW REC_pH P_MM L_DI May–July 1991–1993

Temperature-dependent recycling

C B EQ JULIAN INIT REC_T P_MM L_DI May–July 1991–1993

C′ C EQ JULIAN INIT_NEW REC_T P_MM L_DI May–July 1991–1993

Runge-Kutta solution of chlorophyll a

D B RK JULIAN INIT REC_pH P_MM L_DI May–July 1991–1993

D′ D RK JULIAN INIT_NEW REC_pH P_MM L_DI May–July 1991–1993

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias

(ppb)

RMSE

(ppb) NS

The rows below are repeated from Phase I (Wood and others, 2013, table 12)—Continued

Original model with original calibration

O 0.05 0.53 9.1 NA NA NA 0.94 0.89 16 8.3 28.8 0.88

O′ O 0.05 1.33 9.0 NA NA NA 0.92 0.85 16 9.7 35.9 0.82

Recalibration of original model

A O 0.059 0.47 9.0 NA NA NA 0.93 0.86 16 14.9 36.5 0.81

A′ A 0.059 1.24 8.9 NA NA NA 0.90 0.81 16 13.7 43.1 0.74

Modifications to light and phosphorus limitation

B A 0.035 0.32 8.6 NA NA 0.06 0.94 0.88 16 17.4 39.3 0.78

B′ B 0.035 1.05 8.7 NA NA 0.06 0.96 0.92 16 21.6 36.6 0.81

Temperature-dependent recycling

C B 0.035 NA NA 0.54 12.8 0.06 0.95 0.91 16 25.8 44.1 0.72

C′ C 0.035 NA NA 1.62 12.8 0.06 0.95 0.90 16 26.6 45.6 0.71

Runge-Kutta solution of chlorophyll a

D B 0.035 0.34 8.1 NA NA 0.06 0.94 0.88 16 -6.8 34.4 0.83

D′ D 0.035 1.18 8.1 NA NA 0.06 0.94 0.88 16 -14.3 35.5 0.82
Phase 2 Model Updates 21

Table 8. Performance statistics for total phosphorus obtained over the indicated May–July validation periods for model versions 

discussed in the report.—Continued

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare 

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below were generated for Phase 2

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July 1991–1993

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July even years 

F′O G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July 2006–2008

F′E E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July odd years 

F′E G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July 2006–2008

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA May–July 2006

G′2 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA May–July 2007

G′3 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA May–July 2008

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias

(ppb)

RMSE

(ppb) NS

The rows below were generated for Phase 2 —Continued

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ 0.031 1.31 8.1 NA NA 0.06 0.92 0.85 16 -27.9 43.5 0.73

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ 0.070 1.61 8.2 NA NA 0.06 0.78 0.61 60 1.2 48.6 0.59

F′O G′ 0.070 1.61 8.2 NA NA 0.06 0.74 0.54 19 26.4 71.1 0.29

F′E E′ 0.087 2.09 8.3 NA NA 0.05 0.80 0.64 58 18.4 50.7 0.53

F′E G′ 0.087 2.09 8.3 NA NA 0.05 0.76 0.58 19 28.5 67.0 0.37

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ 0.463 6.77 8.9 NA NA NA 0.83 0.69 5 -0.3 42.5 0.69

G′2 F′ 0.463 7.63 9 NA NA NA 0.99 0.99 4 34.6 55.8 0.70

G′3 F′ 0.463 6.76 8.9 NA NA NA 0.95 0.90 5 -26.7 34.0 0.71
22 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Table 9. Performance statistics for chlorophyll a obtained over the indicated May–July validation periods for model versions 

discussed in the report.

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below are repeated from Phase I (Wood and others, 2013, table 13)

Original model with original calibration

O EQ JULIAN INIT REC_pH P_00 L_AVG May–July 1991–1993

O′ O EQ JULIAN INIT_NEW REC_pH P_00 L_AVG May–July 1991–1993

Recalibration of original model

A O EQ JULIAN INIT REC_pH P_00 L_AVG May–July 1991–1993

A′ A EQ JULIAN INIT_NEW REC_pH P_00 L_AVG May–July 1991–1993

Modifications to light and phosphorus limitation

B A EQ JULIAN INIT REC_pH P_MM L_DI May–July 1991–1993

B′ B EQ JULIAN INIT_NEW REC_pH P_MM L_DI May–July 1991–1993

Temperature-dependent recycling

C B EQ JULIAN INIT REC_T P_MM L_DI May–July 1991–1993

C′ C EQ JULIAN INIT_NEW REC_T P_MM L_DI May–July 1991–1993

Runge-Kutta solution of chlorophyll a

D B RK JULIAN INIT REC_pH P_MM L_DI May–July 1991–1993

D′ D RK JULIAN INIT_NEW REC_pH P_MM L_DI May–July 1991–1993

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias

(ppb)

RMSE

(ppb) NS

The rows below are repeated from Phase I (Wood and others, 2013, table 13)—Continued

Original model with original calibration

O 0.05 0.53 9.1 NA NA NA 0.74 0.55 16 -48.9 85.4 0.28

O′ O 0.05 1.33 9.0 NA NA NA 0.75 0.56 16 -48 82.3 0.33

Recalibration of original model

A O 0.059 0.47 9.0 NA NA NA 0.75 0.56 16 -39.6 77.8 0.40

A′ A 0.059 1.24 8.9 NA NA NA 0.75 0.56 16 -40.6 78.0 0.40

Modifications to light and phosphorus limitation

B A 0.035 0.32 8.6 NA NA 0.06 0.77 0.59 16 -37.6 77.1 0.41

B′ B 0.035 1.05 8.7 NA NA 0.06 0.77 0.59 16 -33.4 76.2 0.43

Temperature-dependent recycling

C B 0.035 NA NA 0.54 12.8 0.06 0.75 0.56 16 -27.5 74.0 0.46

C′ C 0.035 NA NA 1.62 12.8 0.06 0.75 0.56 16 -26.6 73.1 0.47

Runge-Kutta solution of chlorophyll a

D B 0.035 0.34 8.1 NA NA 0.06 0.81 0.65 16 11.6 61.4 0.63

D′ D 0.035 1.18 8.1 NA NA 0.06 0.80 0.65 16 7.6 61.8 0.62
Phase 2 Model Updates 23

Table 9. Performance statistics for chlorophyll a obtained over the indicated May–July validation periods for model versions 

discussed in the report.—Continued

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below were generated for Phase 2

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July 1991–1993

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July even years 

F′O G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July 2006–2008

F′E E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July odd years 

F′E G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July 2006–2008

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA May–July 2006

G′2 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA May–July 2007

G′3 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA May–July 2008

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias

(ppb)

RMSE

(ppb) NS

The rows below were generated for Phase 2—Continued

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ 0.031 1.31 8.1 NA NA 0.06 0.80 0.64 16 -1.3 64.6 0.59

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ 0.070 1.61 8.2 NA NA 0.06 0.67 0.45 60 32.2 75.1 0.32

F′O G′ 0.070 1.61 8.2 NA NA 0.06 0.44 0.19 19 18.1 83.2 0.14

F′E E′ 0.087 2.09 8.3 NA NA 0.05 0.62 0.39 58 27.6 70.3 0.27

F′E G′ 0.087 2.09 8.3 NA NA 0.05 0.42 0.18 19 15.9 84.5 0.11

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ 0.463 6.77 8.9 NA NA NA 0.96 0.92 5 -23.7 28.2 0.64

G′2 F′ 0.463 7.63 9 NA NA NA 0.94 0.89 4 -30.1 35.8 0.58

G′3 F′ 0.463 6.76 8.9 NA NA NA 0.83 0.69 5 30.7 87.3 0.55
24 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Table 10. Performance statistics for pH obtained over the indicated May–July validation periods for model versions discussed in the 

report.

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below are repeated from Phase I (Wood and others, 2013, table 14)

Original model with original calibration

O EQ JULIAN INIT REC_pH P_00 L_AVG May–July 1991–1993

O′ O EQ JULIAN INIT_NEW REC_pH P_00 L_AVG May–July 1991–1993

Recalibration of original model

A O EQ JULIAN INIT REC_pH P_00 L_AVG May–July 1991–1993

A′ A EQ JULIAN INIT_NEW REC_pH P_00 L_AVG May–July 1991–1993

Modifications to light and phosphorus limitation

B A EQ JULIAN INIT REC_pH P_MM L_DI May–July 1991–1993

B′ B EQ JULIAN INIT_NEW REC_pH P_MM L_DI May–July 1991–1993

Temperature-dependent recycling

C B EQ JULIAN INIT REC_T P_MM L_DI May–July 1991–1993

C′ C EQ JULIAN INIT_NEW REC_T P_MM L_DI May–July 1991–1993

Runge-Kutta solution of chlorophyll a

D B RK JULIAN INIT REC_pH P_MM L_DI May–July 1991–1993

D′ D RK JULIAN INIT_NEW REC_pH P_MM L_DI May–July 1991–1993

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias RMSE NS

The rows below are repeated from Phase I (Wood and others, 2013, table 14)—Continued

Original model with original calibration

O 0.05 0.53 9.1 NA NA NA 0.75 0.56 16 -0.443 0.659 0.19

O′ O 0.05 1.33 9.0 NA NA NA 0.73 0.54 16 -0.456 0.677 0.14

Recalibration of original model

A O 0.059 0.47 9.0 NA NA NA 0.75 0.56 16 -0.417 0.644 0.22

A′ A 0.059 1.24 8.9 NA NA NA 0.73 0.53 16 -0.435 0.666 0.17

Modifications to light and phosphorus limitation

B A 0.035 0.32 8.6 NA NA 0.06 0.74 0.55 16 -0.320 0.589 0.35

B′ B 0.035 1.05 8.7 NA NA 0.06 0.78 0.60 16 -0.279 0.541 0.45

Temperature-dependent recycling

C B 0.035 NA NA 0.54 12.8 0.06 0.77 0.60 16 -0.300 0.551 0.43

C′ C 0.035 NA NA 1.62 12.8 0.06 0.78 0.60 16 -0.299 0.549 0.43

Runge-Kutta solution of chlorophyll a

D B 0.035 0.34 8.1 NA NA 0.06 0.86 0.73 16 0.061 0.386 0.72

D′ D 0.035 1.18 8.1 NA NA 0.06 0.86 0.74 16 0.043 0.382 0.73
Phase 2 Model Updates 25

Table 10. Performance statistics for pH obtained over the indicated May–July validation periods for model versions discussed in the 

report.—Continued

[Model characteristics: EQ, Chlorophyll a concentration evaluated by setting dB/dt to 0; RK, Chlorophyll a concentration evaluated using Runge-Kutta 

approach; JULIAN, pH model with Julian day dependence; NO_JULIAN, pH model without Julian day dependence; INIT, Initial size of sediment reservoir 

calculated with original starting concentration 1,000 milligrams of phosphorus per kilogram of sediment; INIT_NEW, Initial size of sediment reservoir 

recalculated with starting concentration 335 milligrams of phosphorus per kilogram of 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 phosphorus; 

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 calibration parameters: Model characteristics and 

parameters are defined in the text and tables 1 and 2. Performance statistics: R, correlation coefficient; R2, coefficient of determination; N, sample size; RMSE, 

root mean square error; NS, Nash-Sutcliffe statistic. Abbreviations: m d-1, meter per day; yr-1, year; °C, degrees Celsius; ppm, parts per million]

Model 

version

Compare

to version

Model characteristics

Chlorophyll a Validation period

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation

Light 

limitation

The rows below were generated for Phase 2

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July 1991–1993

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July even years 

F′O G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July 2006–2008

F′E E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July odd years 

F′E G′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI May–July 2006–2008

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA May–July 2006

G′2 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA May–July 2007

G′3 F′ RK NO_JULIAN INIT_NEW REC_pH NA NA May–July 2008

Model 

version

Compare 

to version

Model calibration parameters Performance statistics

ualg

(m d-1)

KR

(yr-1) pH* KT(yr-1)

T(min,rec)

(°C)

khalf

(ppm) R R2 N Bias RMSE NS

The rows below were generated for Phase 2—Continued

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ 0.031 1.31 8.1 NA NA 0.06 0.86 0.74 16 0.114 0.395 0.71

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ 0.070 1.61 8.2 NA NA 0.06 0.81 0.66 60 0.177 0.463 0.59

F′O G′ 0.070 1.61 8.2 NA NA 0.06 0.67 0.45 19 0.122 0.541 0.42

F′E E′ 0.087 2.09 8.3 NA NA 0.05 0.77 0.60 58 0.232 0.518 0.48

F′E G′ 0.087 2.09 8.3 NA NA 0.05 0.68 0.46 19 0.128 0.539 0.43

Runge-Kutta solution of chlorophyll a with calculated metabolism 2006–2008

G′1 F′ 0.463 6.77 8.9 NA NA NA 0.88 0.77 5 -0.345 0.450 0.36

G′2 F′ 0.463 7.63 9 NA NA NA 0.72 0.52 4 0.232 0.455 0.35

G′3 F′ 0.463 6.76 8.9 NA NA NA 0.98 0.96 5 0.068 0.314 0.77
26 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

tac15-0975_fig02

0

100

200

300

400

500

600

700

EXPLANATION

Model version

E’

F’O

F’E

0

100

200

300

400

500

1990 1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012

7.5

7.0

8.0

8.5

9.0

9.5

10.0

10.5

11.0

A. Total phosphorus

B. Chlorophyll a

C. pH

Year

Total phosphorus, 

in parts per billion

Chlorophyll a, 

in parts per billion

pH, in standard units

Measured

Figure 2. Simulated and measured water column (A) total phosphorus (as P) concentration, (B) chlorophyll a

concentration, and (C) pH, Upper Klamath Lake, Oregon, April 1991–September 2010.
Phase 2 Model Updates 27

Long-Term Simulations

In order to better understand how long the system takes 

to achieve a steady state between the water column and 

sediment reservoir, long-term simulations were made with 

the recalibrated models (versions E', F'O, and F'E; table 11). 

The long-term simulations were achieved by chaining the 

19 years of input data (mid-April 1991–mid-April 2010) 

to produce a continuous simulation of all variables over 

209 years (eleven 19-year cycles). The recalibrated models 

were run with the chained input data for the base scenario of a 

0-percent reduction in external phosphorus load and a TMDL 

scenario of a 40-percent reduction in external phosphorus 

load. Results as represented by 19-year averages were 

consistent between models (fig. 3) and previous results (Wood 

and others, 2013; fig. 16, represented as 7-year averages). In 

all model runs with the updated sediment phosphorus initial 

condition, the model achieved steady state in about 20 years, 

and the 19-year average water-column TP concentrations 

were about 125 ppb for the base scenario and about 75 ppb 

for the TMDL scenario. When the 19-year average sediment 

TP masses were compared (fig. 4), the trends were consistent 

between models and previous results (Wood and others, 

2013; fig. 17, represented as 7-year averages); however, the 

ultimate steady-state sediment phosphorus masses varied. 

As for the water column, steady state was achieved in about 

20 years; however, the model recalibrated with the extended 

dataset (versions F'O and F'E) had lower final sediment 

phosphorus mass than the model calibrated with the original 

7.5-year dataset (version E'). This lower mass was related 

to faster sediment phosphorus depletion rates because of the 

combination of calibrated recycling parameters pH* and KR. 

The recalibrated extended set models (F'O and F'E) achieved a 

steady state of about 500 metric tons of sediment phosphorus 

with the base scenario, and a steady state of about 350 metric 

tons of sediment phosphorus with the TMDL scenario. By 

comparison, the model over the original 7.5-year time period 

(version E') achieved a steady state of about 400 metric tons 

of sediment phosphorus with the base scenario, and a steady 

state of about 250 metric tons of sediment phosphorus with the 

TMDL scenario (table 11). 

In addition to considering the amount of time required 

to achieve steady state, we evaluated the “final” steadystate 19-year time series at the end of the 209-year model 

runs. Water-column TP (fig. 5) and chlorophyll a (fig. 6) 

concentrations decreased under the TMDL scenario as 

compared to the baseline (0-percent reduction in external loads) 

scenario. The decrease in water-column TP concentration was 

39 percent, roughly the same as the reduction in external loads, 

and the decrease in chlorophyll a concentration was a little 

less but roughly similar. Across all versions of the model, the 

TMDL scenario resulted in a similar decrease in water-column 

phosphorus concentration compared to the baseline scenario 

(table 11). Note that, although chlorophyll a concentrations 

decreased in response to external load reductions, there were 

still blooms in each year (fig. 6). The blooms are smaller, and 

many are under the 100 ppb threshold that was applied in the 

Phase I review (Wood and others, 2013); however, they still 

occurred. Examination of the model parameters indicated 

that the recycling parameter pH* plays an important role in 

determining whether a bloom year will be “turned on” or 

“turned off.” When the pH* value is low, sediment phosphorus 

recycling begins at a lower lake-average pH and activates the 

feedback loop between sediment recycling and chlorophyll a. 

Because pH is a function of chlorophyll a, the recycling will 

increase rapidly and lead to increased chlorophyll a as long as 

the chlorophyll a concentration at the beginning of the season is 

high enough given growth conditions (that is, water-column TP, 

light, and temperature). When the pH* value is high, sediment 

recycling cannot begin until the lake-average pH also is 

relatively high. For example, if pH* equals 8.1, then recycling 

would begin (at a low rate) when the lakewide average pH 

reached 7.5 and the recycle rate would reach its maximum 

when the lakewide average pH reached 8.8; however, if pH* 

equals 9.1, then recycling would begin (at a low rate) when the 

lakewide average pH reached 8.5 and the maximum recycle rate 

would occur when lakewide average pH reached 9.8. Therefore, 

the models that remove the assumption of cyanobacteria 

achieving maximum growth during the 2-week time step and 

solve for chlorophyll a using a Runge-Kutta solution method 

will be less likely to have years without any bloom. 
28 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Table 11. Metrics describing the long-term simulations for various model versions discussed in the report.

[Model characteristics: 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 milligrams of phosphorus per kilogram of sediment; 

INIT_NEW, Initial size of sediment reservoir recalculated with starting concentration 335 milligrams of phosphorus per kilogram of sediment; REC_pH, Model 

relies on pH-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. %, percent; >, greater than; ppb, parts per billion]

Model 

version

Compare

to version

Model characteristics

Number of 

“final” years Chlorophyll a

model pH model Initial sediment 

reservoir

Recycle 

mechanism

Phosphorus 

limitation Light limitation

The rows below are repeated from Phase I (Wood and others, 2013, table 15)

Original model with original calibration

O′ O EQ JULIAN INIT_NEW REC_pH P_00 L_AVG 7

Recalibration of original model

A′ A EQ JULIAN INIT_NEW REC_pH P_00 L_AVG 7

Modifications to light and phosphorus limitation

B′ B EQ JULIAN INIT_NEW REC_pH P_MM L_DI 7

Runge-Kutta solution of chlorophyll a

D′ D RK JULIAN INIT_NEW REC_pH P_MM L_DI 7

Model 

version

Compare

to version

0% 40% 0% 40% 0% 40% 0% 40%

Average concentration over the final years of the 

long-term simulation (last 7 years of 203-year simulation 

or last 19 years of 209-year simulation)

Years with seasonal chlorophyll a peak >100 ppb over 

final years of the long-term simulation

Total P 

(ppb)

Sediment P 

(metric tons)

Chlorophyll a

(ppb)

Number

(out of number of “final” years)

The rows below are repeated from Phase I (Wood and others, 2013, table 15)—Continued

Original model with original calibration

O′ O 120.0 71.1 901 932 71.3 44.5 7 6

Recalibration of original model

A′ A 118.4 69.4 980 897 69.1 44.0 7 6

Modifications to light and phosphorus limitation

B′ B 121.5 73.0 871 1,028 72.0 39.8 7 4

Runge-Kutta solution of chlorophyll a

D′ D 129.0 78.4 518 355 60.3 38.4 7 6
Phase 2 Model Updates 29

Model 

version

Compare

to version

Model characteristics

Number of 

“final” years Chlorophyll a

model pH model Initial sediment 

reservoir

Recycle

mechanism

Phosphorus 

limitation

Light

limitation

The rows below were generated for Phase 2

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI 7

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI 19

F′E E′ RK NO_JULIAN INIT_NEW REC_pH P_MM L_DI 19

Model 

version

Compare

to version

0% 40% 0% 40% 0% 40% 0% 40%

Average concentration over the final years of the 

long-term simulation (last 7 years of 203-year simulation 

or last 19 years of 209-year simulation)

Years with seasonal chlorophyll a peak >100 ppb over 

final years of the long-term simulation

Total P 

(ppb)

Sediment P 

(metric tons)

Chlorophyll a

(ppb)

Number

(out of number of “final” years)

The rows below were generated for Phase 2—Continued

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 7.5-year calibration set

E′ D′ 130.1 79.6 387 251 57.0 35.3 7 4

Runge-Kutta solution of chlorophyll a with updated settling velocity and model structure, 19-year calibration set

F′O E′ 120.3 73.5 505 334 42.4 26.5 19 3

F′E E′ 120.5 73.9 521 357 44.2 28.1 19 13

Table 11. Metrics describing the long-term simulations for various model versions discussed in the report.—Continued

[Model characteristics: 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 milligrams of phosphorus per kilogram of sediment; 

INIT_NEW, Initial size of sediment reservoir recalculated with starting concentration 335 milligrams of phosphorus per kilogram of sediment; REC_pH, Model 

relies on pH-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. %, percent; >, greater than; ppb, parts per billion]
30 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

tac15-0975_fig03

50 150

0

50

100

150

200

A. 0-percent reduction

0 100 200 250 0 50 100 150 200 250

B. 40-percent reduction

EXPLANATION

19−year average water column total phosphorus, in parts per billion

Years of simulation

Model version

E’

F’O

F’E

Figure 3. 19-year average water-column total phosphorus (as P) 

concentration under simulation conditions of (A) 0-percent reduction and 

(B) 40-percent reduction in external phosphorus loads to Upper Klamath 

Lake, Oregon.
Phase 2 Model Updates 31

tac15-0975_fig04

0

200

400

600

800

1,000

0 50 100 150 200 250 0 50 100 150 200 250

Years of simulation

19-year average sediment total phosphorus, in metric tons

EXPLANATION

Model version

E’

F’O

F’E

A. 0-percent reduction B. 40-percent reduction

Figure 4. 19-year average sediment total phosphorus (as P) masses 

under simulation conditions of (A) 0-percent reduction and (B) 40-percent 

reduction in external phosphorus loads to Upper Klamath Lake, Oregon.
32 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

tac15-0975_fig05

2180 2185 2190 2195 2200 2205

Year

0

50

100

150

200

250

300

Total phosphorus, in parts per billion

A. 0-percent reduction B. 40-percent reduction

2180 2185 2190 2195 2200 2205

EXPLANATION

Model version

F’O

F’E

Figure 5. Water-column total phosphorus (as P) concentration during the last 19 years of a 209-year 

simulation (1991–2199) under conditions of (A) 0-percent reduction and (B) 40-percent reduction in external 

loads to Upper Klamath Lake, Oregon.
Phase 2 Model Updates 33

tac15-0975_fig06

Year

2180 2185 2190 2195 2200 2205

0

50

100

150

200

2180 2185 2190 2195 2200 2205

Chlorophyll a, in parts per billion

EXPLANATION

A. 0-percent reduction B. 40-percent reduction

Model version

F’O

F’E

Figure 6. Water-column chlorophyll a concentration during the last 19 years of a 209-year simulation 

(1991–2199) under conditions of (A) 0-percent reduction and (B) 40-percent reduction in external loads to 

Upper Klamath Lake, Oregon.
34 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Lake Metabolism as a Substitute for Algal Growth and Respiration

Perhaps the largest source of error in using the CSTR approach to model a large system 

that is variable in all three dimensions and time is the process by which descriptions of 

processes that are appropriate to small, homogeneous water parcels are “scaled up” and 

applied to the bulk system. This can be shown by examining the mathematical procedure for 

averaging terms of the three-dimensional mass balance equations to make them appropriate for 

a zero-dimensional CSTR model. 

As an example, consider the growth and respiration terms in the mass balance equation 

for chlorophyll a (equation 1). The growth process is described as a product of a constant 

maximum growth rate, modified by factors that describe how the growth rate is limited by 

less-than-ideal amounts of light, temperature, and nutrients. This description has been shown 

over decades to work well, and when used in process-based models to describe growth in 

discrete and approximately homogeneous segments of a model domain that are selected to 

have the appropriate temporal and spatial scales (Bowie and others, 1985; Chapra, 1997). The 

application of the small-scale description of growth to the single, large, averaged volume that 

represents the lake in the CSTR model has been done by making the following approximations 

in milligrams of chlorophyll a per liter per day:

 ( ) ( ) ( ) ( ) ( ) ( ) 1 1 ˆ D A t GL T P K F I F T F P Bdz ds dt K F T F P B ˆ ˆˆ ~ GT P D tL 0 F I dz dt

DA t D t ∆ ∆ ∗ ∆ ∆ ∫∫∫ ∫ ∫ (2)

 ( ) ( ) ˆ ˆ 1 D A t RT K F T Bdz ds dt K F T B ~ R T DA t ∆ ∆ ∫∫∫ (3)

where

s represents two-dimensional (latitude/longitude) geographic location,

1 D D

dz ∫ represents averaging over the water-column depth, 

1 A A

ds ∫ represents averaging over the area of the lake,

1 t t dt ∆ ∆ ∫ represents averaging over the 14-day time step of the model, and 

integrated functions are GROWTH and RESP as described with terms in section, “Changes 

to Algal Submodel.” The “hat” on a variable represents values with spatial and temporal 

dependence (functions of depth, area, and time) and distinguishes it from the corresponding 

CSTR variable that has been averaged so that it has no spatial dependence and takes on values 

only at 14-day intervals. For example:

 ( ) ( ) 1 ,, ˆ DA tB s z t dz ds dt B t

DA t ∆ ∆ = ′ ∫∫∫ (4)

where 

t' is a discrete time variable that is defined at 14-day intervals. 

Inaccuracy results because the spatial/temporal average of the product of multiple factors is 

not, in general, equal to the product of the spatial/temporal average of the individual factors. 

The amount of error introduced by the approximations depends on the spatial and temporal 

increments used for averaging—in this case, averaging is over 14 days and the entirety 

of UKL—and in a highly variable system, approximations (equations) 2 and 3 introduce 

more error than in a well-mixed system. Approaching the scaling-up problem in this way 

is convenient because it provides a way of parsing the growth and respiration terms into 

well-understood factors.
Phase 2 Model Updates 35

In this study, we have considered an alternative means of 

scaling up the description of growth and respiration. Consider 

that the units of equation 2 are milligrams chlorophyll a per 

liter per day. When multiplied by an appropriate stoichiometric 

coefficient δChl:C that converts milligrams chlorophyll a to 

milligrams carbon, the left-hand side of equation 2 is equal to 

lakewide average gross primary production (GPP) expressed 

in milligrams of carbon per liter per day, and when δChl:C is 

multiplied to the left had side of equation 3, the result is 

equal to the lakewide average of the autotrophic component 

of community respiration (CR), expressed in milligrams of 

carbon per liter per day. 

The use of continuous (hourly to subhourly) 

measurements of DO to estimate aquatic ecosystem 

metabolism (Odum, 1956) has been increasingly refined for 

lotic ecosystems (Chapra and Di Toro, 1991; Bales and Nardi, 

2007; Aristegi and others, 2009; Raymond and others, 2012). 

Similar methods for lentic environments have lagged behind, 

partly because of challenges in accounting for advective 

horizontal transport and physical reaeration contributing to 

spatially varying DO concentrations in lake surface waters. In 

a large productive lake, horizontal transport of DO was shown 

to lead to significant error when the estimates are made at only 

one sensor (Gelda and Effler, 2002). Antenucci and others 

(2013) used a modeling approach to determine that transport 

processes could lead to large errors if not accounted for, and 

advocated the use of non-dimensional numbers to characterize 

the likely importance of horizontal and vertical transport ahead 

of time. However, recent advances (Staehr and others, 2010) 

have improved the applications of these aquatic metabolic 

methods to lakes, offering opportunities to evaluate and 

model lake ecosystems as functions of physical and chemical 

variables that can be independently measured.

We propose the use of continuous DO measurements to 

scale up growth and respiration of the cyanobacterial bloom, 

acknowledging that the application of the technique to a 

very large, shallow lake, with strong wind-driven currents, 

substantial spatial variability, buoyant cyanobacteria, and a 

cycle of diurnal stratification and mixing does not completely 

capture the spatial variability and transport processes. As 

part of this study, an initial attempt was made to calculate the 

lakewide metabolism from continuous DO measurements 

taken during 2006–08. Efforts were made to mitigate the 

problem of vertical mixing in this initial attempt by limiting 

the analysis to sites shallower than the diurnally mixed layer. 

In doing so the estimates were weighted more heavily toward 

littoral zones where both GPP and CR are expected to be 

higher than in the pelagic zone (Van de Bogert and others, 

2007), so the estimates, particularly of GPP, are likely to be 

biased higher than the true lakewide mean. The estimates 

of CR may be higher or lower than the true lakewide mean 

because we have not included sites where the water column 

extends below the diurnally mixed layer. We have attempted 

to limit the potential bias introduced by horizontal transport by 

averaging the results at several locations.

For the purpose of substituting lake metabolism for the 

growth and respiration terms in the CSTR model, simply 

calculating lake metabolism from measurements is not enough. 

In order to simulate the response to nutrient management 

scenarios, GPP and CR must be predicted (through a submodel) 

as an explicit function of readily observable, independent 

variables. In their original form, the growth and respiration 

terms in the model are functions of the 14-day values of 

lakewide averages of chlorophyll a, phosphorus, temperature, 

and incident light. In the process of substituting lake 

metabolism terms, the suite of possible explanatory variables 

can be expanded to include other variables that are known 

or hypothesized to influence lake metabolism—in particular, 

meteorological variables such as wind and air temperature, 

and water-column stability. This would enable the model to 

take advantage of the higher frequency of data collection of 

meteorological and water quality variables that are collected 

hourly or sub-hourly to improve the time integration through the 

14-day time step. 

The complex dependencies between bloom dynamics 

and meteorological variables (and water temperature, a direct 

consequence of meteorology) at UKL have been noted in 

several previous studies (for example, Wood and others, 1996, 

2006; Morace, 2007). Air temperature and wind speed also 

determine water-column stratification, which has been related 

to aspects of bloom dynamics as well (Kann and Welch, 2005; 

Wood and Gartner, 2010). Meteorology was included in the 

algorithms of the original TMDL model through light limitation 

to growth and temperature limitation to growth and respiration, 

the form of which was broadly prescribed over the season. This 

study derived empirical functions relating (on shorter, daily 

time scales) algal growth and respiration to directly measured 

meteorological (wind and solar radiation) and water temperature 

variables, and measured indicators of the strength of watercolumn stratification based on the difference in temperature 

between the upper and lower water column. The goals were to 

derive relationships that could be incorporated into the TMDL 

model equations and to improve the accuracy and timing in 

the simulation of the state variables (chlorophyll a and TP 

concentrations, and pH) and the prediction of lake response 

to implementation of TMDL limits on external phosphorus 

loading. As noted in the TMDL model documentation (Kann 

and Walker, 1999) and in the Phase 1 report (Wood and others, 

2013), the original and new versions (through version F'E) of the 

TMDL model cannot simulate rapid declines in the bloom that 

typically occur after the second week in July, and, as a result, 

the performance of the model decreases after the early part of 

the growing season (May–July).

The approach used is an empirical investigation of the 

production and CR time series, calculated from the continuous 

monitors. Although production and CR will have broad seasonal 

dependencies on temperature, light, nutrients, and biomass 

that could be determined empirically from the 14-day dataset 

used to run the TMDL model, in this case the dependencies 

were not assumed to take the form prescribed in the TMDL 

equations. After determining the broad seasonal dependencies, 

the remaining variability on shorter time scales was investigated 

with a daily dataset as generated below.
36 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Explanatory variables were restricted to those that were 

independently observable and were not expected to change 

with nutrient management, or those that were state variables 

and thus were calculated internally by the model itself. 

Chlorophyll a and phosphorus were potential explanatory 

variables because these were state variables in the model; 

incident light, wind, and air temperature were included 

because they were meteorological variables that were not 

affected by nutrient management, and water temperature was 

allowed as a water quality variable because it is observable but 

should not change with nutrient management. In contrast, DO 

concentration was not included because it could change with 

nutrient management and is not predicted internally by the 

model. In this study, we have explored this concept by using 

our calculated GPP and CR values to develop multivariate 

regression models that could be substituted into the CSTR 

mass balance equation for chlorophyll a.

Calculation of Upper Klamath Lake Metabolism 

from Continuous Dissolved Oxygen, 2006–08

Lake Metabolism Calculations

In order to evaluate the metabolism occurring in the 

surface waters of the lake because of cyanobacterial activity, 

continuous monitor data collected over 2006–08 were used to 

quantify GPP (photosynthesis) and CR. Hourly DO and lake 

temperature data collected from nine shallow water-quality 

(WQ) sites and hourly wind speed data collected from two 

meteorological (MET) sites in the lake were used for this 

purpose (table 12). Because continuous WQ monitors are 

removed during the winter, the full daily sets covered 

May 15–October 16, 2006, May 15–November 13, 2007, 

and May 5–October 13, 2008. We focused on the shallow 

sites (less than 3.0 m full-pool depth) so that effects of vertical 

transport would be minimized. WQ sites were characterized by 

their position in the lake as “northern” and “southern” in order 

to assign proximate wind speed values, necessary for reaeration 

calculations, and a lakewide sediment oxygen demand (SOD) 

areal flux rate of 1.6 grams of oxygen per square meter per day 

[g (O2

/m2)/d] at 20 °C (Wood and others, 2013) was assumed 

for all sites. Hourly production and CR at each site then were 

evaluated by the following equation at each hour:

∆

∆ = − − +

C

t zfull Φ Φ photo resp S Φ Φ OD reaer (5)

To quantify reaeration, wind speed was used to calculate a 

film transfer coefficient and applied that transfer rate to the 

difference between saturation DO and measured DO. When 

saturation DO was less than measured DO, oxygen moved 

from the water to the air, and when saturation DO was greater 

than measured DO, oxygen moved from the air to the water. 

Saturation DO was calculated for each site from hourly water 

temperature data using the Benson and Krause equations 

(U.S. Geological Survey, 2011). The film transfer coefficient, 

in meters per day, was calculated from wind speed data using 

the following equation (Gelda and Effler, 2002):

10 10

,20 2

10 10

0.2 , 3.5 m/s

0.057 , 3.5 m/s L

WSPD WSPD

K

WSPD WSPD

 ≤ = 

 >

(6)

Because USGS MET wind speed data were measured at a 

height of 2 m above water surface, the values were translated 

to 10-m values, WSPD10, using a logarithmic approximation 

(Martin and McCutcheon, 1999). The film transfer rates were 

temperature corrected using the Van ‘t Hoff equation, with 

input parameter, θ, equal to 1.024 (Gelda and Effler, 2002), 

and applied to the DO deficit to calculate the areal reaeration 

flux for each hour:

( ) 20

,20 *1.024 * TL sat

reaer

full

K C C

z

− − Φ = (7)

The sediment-oxygen demand (SOD) areal flux also was 

temperature-corrected by the Van ‘t Hoff equation, with θ 

equal to 1.065 (Chapra, 1997), to get ΦSOD. At this point, the 

nighttime DO data (9:00 p.m.–5:00 a.m.) were extracted to 

calculate hourly respiration flux by setting photosynthetic flux 

equal to zero:

Φ Φ resp reaer full S Φ OD

C

t = − z ∆

∆ * − (8)

For each hour, respiration was reverse-temperature-corrected 

to find Φresp,20 using equation 15 from Wood and others (2013), 

with Tmin set to 4 °C instead of 14 °C. This allowed respiration 

to continue at lower temperatures than previously allowed 

in the TMDL model. Φresp,20 then was averaged over each 

continuous nighttime segment and the average value was 

applied to all hours of the prior day. For example, respiration 

averaged from June 1, 2006, 9:00 p.m., to June 2, 2006, 

5:00 a.m., was applied to all hours of June 1, 2006. With 

Φresp,20 values now temperature-corrected for each hour, hourly 

photosynthetic flux was calculated: 

Φ Φ photo full resp S Φ Φ OD reaer

C

t = z ∆∆ * + + − (9)

Daily GPP and respiration were calculated for each site 

by accumulating the hourly values, and the lakewide 

totals were determined by averaging across all sites for 

each day. Daily totals then could be converted to chlorophyll 

a concentrations by using the average lakewide depth 

and published stoichiometric ratios of oxygen to carbon 

(2.67 gO/gC) and carbon to chlorophyll a (41 gC/g 

chlorophyll a) (Bowie and others, 1985). The net chlorophyll a

concentrations (BNPP) were calculated by taking the difference 

between these daily GPP and respiration chlorophyll a

concentration totals.
Phase 2 Model Updates 37

Table 12. U.S. Geological Survey continuous monitoring sites used for metabolism analysis, Upper Klamath Lake, Oregon, 

2006–08.

[Site locations are shown in figure 1. Abbreviations: USGS, U.S. Geological Survey; WQ, water quality; MET, meterological; NA, not applicable; 

m, meter]

Site name Site name 

abbreviation Type USGS site No.

Latitude 

(north)

Longitude 

(west)

Full-pool 

measured 

depth 

(m)

Years Lake 

orientation

Fish Banks FBS WQ 422808122024400 42°28'08" 122°02'44" 2.80 2006–2008 North

North Buck Island NBI WQ 421838121513900 42°18'38" 121°51'39" 2.80 2006–2008 South

Upper Klamath Lake at 

Williamson River outlet

WMR WQ 422719121571400 42°27'19" 121°57'14" 2.50 2006–2008 North

Howard Bay HDB WQ 421933121550000 42°19'33" 121°55'00" 2.20 2006–2008 South

Hagelstein Park HPK WQ 422319121585700 42°23'19" 121°48'57" 2.60 2006 South

Skillet Handle SHL WQ 421746121522800 42°17'46" 121°52'28" 2.50 2006 South

South Shore SSR WQ 421410121492000 42°14'10" 121°49'20" 2.50 2006–2008 South

Goose Bay East GBE WQ 422749121540700 42°27'39" 121°54'08" 2.40 2006–2008 North

Williamson Delta West WDW WQ 422842121584300 42°28'48" 121°58'43" 2.20 2006–2008 North

Mid-North MDN MET 422622122004000 42°26'22" 122°00'40" NA 2006–2008 North

Mid-Lake MDL MET 422312121515900 42°23'12" 121°51'59" NA 2006–2008 South

Multivariate Regression Models of Lake Metabolism

The time series of daily calculated production, due to photosynthesis, and CR were modeled at two different temporal 

scales. We calculated the sum of a seasonal component, which represented variability at scales longer than a week, and a higher 

frequency (designated “weekly”) component, which represented variability at scales shorter than a week: GPP = GPPS + 

GPPW and CR = CRS + CRW. The seasonal component was determined by multivariate regression with the Klamath Tribes 

14-day dataset, which proceeded as follows. The seasonal component was assumed to be a function of chlorophyll a, nonalgal 

phosphorus, temperature, and solar radiation, as these variables were considered to represent long-term variability at scales 

longer than a week. 

 GPP or CR f P T I S S ,, , = (B na lake 0 ) (10)

Each of the four variable inputs were selected from a set which included the variable itself, the natural logarithm of the 

variable and the square of the variable, as shown below for chlorophyll a, which resulted in 81 unique models tested for each of 

GPPS and CRS:

 B ∈{ } B B ,ln , ( ) B2 (11)

All 3 years of data were combined for this step and an Akaike’s Information Criterion (AIC) value was computed for 

each tested model. The set of possible best models was the set of models with an AIC value within 2 of the minimum AIC 

value. A multivariate regression was run using each of these models, and the selected model was the one that resulted in the 

best combination of normally distributed residuals based on a comprehensive assessment of a variety of indicators including 

(1) visual inspection of plots of the residuals; (2) variance inflation factors indicating low collinearity of the explanatory 

variables; (3) a Durbin-Watson statistic close to 2, indicating low serial correlation in the residuals; (4) t-tests of the parameter 

estimates; and (5) parsimony when other indicators were equal. 

The weekly component was determined by multivariate regression with a daily dataset of meteorological variables and 

water temperature as follows: A time series of daily values of the seasonal component of production and CR was calculated 

from daily values of the Klamath Tribes 14-day dataset, obtained with a cubic spline interpolation, and the seasonal regression 

parameters determined in the first step. A times series of the daily values of the weekly component of production and CR was 

obtained by subtracting daily values of the seasonal component from the daily values of the total production and CR, which 

filtered out the seasonal cycle. Five possible explanatory variables were considered (table 1): 

 GPP C W W or R f = ( ) WSPD10 , , WTDEV SSDEV, , WTROC STRAT (12)
38 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

WSPD10 is the daily wind speed 10 m above the water 

surface, as determined previously. WTDEV is the daily 

deviation in water temperature from the seasonal trend, where 

the seasonal trend is the cubic-spline-interpolated water 

temperature from the 14-day dataset; it is calculated as the 

ratio of daily lakewide averaged water temperature divided by 

the interpolated daily values of the seasonal trend in lakewide 

averaged water temperature. SSDEV is the daily deviation in 

shortwave solar radiation from a clear sky calculation of solar 

radiation; it is calculated as the daily average shortwave solar 

radiation divided by the daily clear sky calculated values. The 

clear sky calculation was obtained from Environmental and 

Water Resources Institute (2005). WTROC is the daily rate of 

change in air temperature, defined as the difference between 

the water temperature of the current day and the water 

temperature 3 days prior, divided by 3. A window of 3 days 

was used to determine WTROC because the relation between 

stratification and bloom expansion has a tendency to reverse 

sign at 3 days (Eldridge and others, 2014). STRAT is a variable 

designed to be a surrogate for water-column stratification, 

and is defined as the daily cumulative sum of the difference 

in temperature measured by a datasonde located 1 m from the 

surface and a datasonde located 1 m from the bottom at the 

deepest site in the lake (MDT) monitored with continuous 

monitors. The time series in each year were divided into three 

segments—early, middle, and late season—as determined by 

a temperature rule rather than a calendar date. These segments 

based on temperature were based on a subjective examination 

of plots of chlorophyll a as a function of temperature. 

Early-season started at the beginning of the time series 

and ended when the daily lake-averaged water temperature 

exceeded 20 °C for 5 consecutive days, marking the beginning 

of mid-season. Mid-season transitioned to late-season when 

the lake-averaged water temperature was less than 20 °C for 

5 consecutive days. The 5-day time window was selected 

subjectively by examining the variability in the time series. 

The multivariate modeling procedure was performed for 

each segment. 

Calibration of the weekly component of the model, 

which was performed separately for the early, middle, and 

late seasons, considered only linear combinations of the five 

explanatory variables, and only data from 2007, which were 

selected randomly (fig. 7). The final model was determined 

from AIC values using multi-model inference (Burnham 

and Anderson, 2002). Six sets of models were created, one 

for each variable. Within each set an Akaike weight was 

calculated for each model (Burnham and Anderson, 2002). 

A model-averaged estimate of each model parameter was 

calculated as an Akaike-weighted average of the estimates 

from all models in a set, and the 80-percent confidence 

interval was calculated from a model-averaged estimate of the 

standard error of the parameter. Only those variables for which 

the 80-percent confidence interval of the model-averaged 

estimate of the parameter did not include zero were retained 

in the final model. Once the variables to be included in the 

final model were determined, the multivariate regression was 

run again with only those variables in order to determine final 

values of the model parameters. At this step, a variable was 

dropped from the model if it was found to create collinearity 

as determined by the variable infiltration factor.

The final model (table 13) combining the seasonal with 

mid- and late-season weekly components was used to re-create 

the daily time series of production and CR during 2006–08. 

The performance of the final model was assessed with a 

calculated bias (mean of simulated value minus measured 

value), a RMSE, a NS statistic, and a R2

.

Table 13. Multivariate regression equations.

[Dependent and explanatory variable definitions are provided in table 1. Abbreviations: R2, coefficient of determination; RMSE, root mean square error; 

–, no value]

Dependent 

variable

Model 

period

Number of 

observations Intercept

Coefficients of explanatory 

variables from 14-day dataset

Coefficients of explanatory 

variables from daily dataset

Regression model 

fit statistics

ln(B) I

O Tlake WSPD10 WTDEV SSDEV WTROC STRAT Adjusted 

R2 RMSE

GPPS Entire 29 –29.56 6.11 0.034 – – – – – – 0.63 5.56

CRS Entire 29 –14.21 3.71 – 0.49 – – – – – 0.55 4.45

GPPW Middle 65 – – – – –3.04 – 10.02 – – 0.23 6.06

CRW Middle 65 – – – – –2.44 14.14 – – – 0.44 4.72

GPPW Late 71 –20.09 – – – – 13.19 4.81 – 0.09 0.33 4.07

CRW Late 71 –13.03 – – – – 12.96 – – 0.13 0.16 4.54
Phase 2 Model Updates 39

tac15-0975_fig07

2 4 6 8 2 4 6 8 2 4 6 8

−5

0

5

10

15

Wind speed, in meters per second

0.7 0.9 1.1

−5

0

5

0.7 0.9 1.1 0.7 0.9 1.1

0 20 40 60 0 20 40 60 0 20 40 60

−4

−2

0

2

4

6

8

−2.0 −1.0 0.0 1.0

−5

0

5

−2.0 −1.0 0.0 1.0 −2.0 −1.0 0.0 1.0

LOESS smoothed respiration, in grams of oxygen per square meter per day

Water temperature deviation from 15−day average, as fraction

Stratification, in degrees Celsius per day

EXPLANATION

2006

2007

2008

Early Middle Late

Rate of change in air temperature, in degrees Celsius per day

Figure 7. Local Regression (LOESS) smoothed respiration as a function of regression model of metabolism predictors, 

Upper Klamath Lake, Oregon, during early, middle, and late parts of summer seasons 2006–08.
40 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Incorporation of Upper Klamath Lake Metabolism Submodel into 

CSTR Model

Incorporation into CSTR Model

The calculated net chlorophyll a concentrations, BNPP, were directly incorporated in the 

chlorophyll a component of the CSTR model by accumulating the daily values over the 2-week 

time step and substituting as below:

 dB

dt

GROWTH RESP SETTLE Q

V

B

B

t

SETTLE Q

V

out NPP out = − − − 



 



× = ∆

− + 



 



× B (13)

At this point, the settling velocity and recycling parameters were recalibrated over 2006–08 by 

following processes similar to those outlined by Kann and Walker (1999). Modified versions of 

calibration steps 1 and 2 (above) follow.

Step 1—Algal Settling Velocity

Algal settling velocity was evaluated by comparing calculated and measured chlorophyll a

concentrations when the following model equation was run with measured chlorophyll a data 

from the previous biweekly interval and measured net growth data:

 alg out t t t NPP t t

u Q BB B B t

z V −∆ −∆

  = + − + ×∆     (14)

Algal settling velocity values were assumed to be between 0.001 and 2.0 m/d (Bowie and 

others, 1985), and the best value was determined by correlating chlorophyll a, as calculated 

by equation 14, with all available measured chlorophyll a concentrations within the 2006–08 

time period. The sum of squared errors was minimized with ualg= 0.463 m/d, which is a much 

higher settling velocity (by an order of magnitude) than had previously been calibrated. For 

this model, it was unnecessary to calibrate the half-saturation constant because net growth was 

input as a measured variable and phosphorus limitation was not explicitly calculated. 

Step 2—Internal Recycling Parameters

As before, the recycling equation for the CSTR model required two calibration parameters, 

pH* and KR, which were selected based on minimization of the objective function (Wood and 

others, 2013, equation 24) over a calibration set. For this analysis, calibration was performed 

by a cross-validation leave-one-out method, where 2 full years of the 2006–08 dataset were 

used for calibration, and the third year was used for validation. In this approach, the calibration/

validation process was repeated three times to ensure that each year was included in the 

validation set at least once and included in the calibration set at least once, which resulted 

in three sets of calibrated recycling parameters (tables 5–10, G'1, G'2, G'3) . Collectively, the 

calibrated pH* value ranged from 8.9 to 9.0 and the calibrated KR value ranged from 6.77 to 

7.63/d for 2006–08,and showed little variability between results when different years were 

selected for calibration. Both the pH* and KR values were higher than those determined in the 

previous recalibration, with recycling rates as much as 6.5 times greater than those calculated 

for the CSTR model calibrated over 1991–98. The higher calibrated recycling rate is consistent 

with the higher calibrated settling velocity discussed above because faster settling requires 

faster recycling to explain the same measured dataset.
Phase 2 Model Updates 41

Model Performance

The performance statistics determined for the calibrated 

models run with calculated metabolism (tables 5–10, G'1, G'2, 

G'3) generally indicated good model performance, although 

the evaluation time period was short, so statistics were 

calculated over small sample sets. Performance statistics 

were comparable to those found for model versions F ′

O and 

F ′

E subsetted over the 2006–2008 period. The whole-year 

validation statistics generally indicated that the metabolism 

models were best at simulating TP and worst at simulating 

chlorophyll a, although version G'2 performed better for 

chlorophyll a than TP. For the whole-year validation periods, 

the metabolism model correlation coefficients generally were 

high—(0.44–0.91 for TP, 0.44–0.66 for chlorophyll a, and 

0.70–0.77 for pH (tables 5–7). When the May–July subsets 

were considered, performance statistics were improved. 

The biweekly metabolism models were able to capture the 

decreases in chlorophyll a and pH associated with the midseason bloom decline in 2006 and 2007, but not in 2008. This 

phenomenon was not achieved in any years with the previous 

CSTR models (fig. 8). 

 G’1

 G’2

 G’3

EXPLANATION

Model version

F’O

F’E

Measured

M J J A S O DN J F M M J J A S O DN A J F M M J J A S A O

2006 2007 2008

Year

0

100

200

300

400

500

600

700

0

100

200

300

400

7.5

7.0

8.0

8.5

9.0

9.5

10.0

10.5

11.0

Total phosphorus, 

in parts per billion

Chlorophyll a, 

in parts per billion pH, in standard units

A. Total phosphorus

B. Chlorophyll a

C. pH

Figure 8. Biweekly simulated and measured water-column (A) total phosphorus (as P) concentration, (B) chlorophyll 

a concentration, and (C) pH, Upper Klamath Lake, Oregon, 2006–08.
42 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

In addition to running the CSTR biweekly model with 

calculated metabolism, the model also was run with a daily 

time step using daily boundary conditions as inputs (Walker 

and others, 2012) and previously calibrated parameters. Model 

performance statistics were not calculated separately, but plots 

comparing the outputs and measured data were developed. 

The daily time step is better at capturing the mid-season 

chlorophyll a and pH decreases than the biweekly time step 

(fig. 9). Daily pH as simulated by the CSTR-metabolism 

model was compared to daily average pH as determined 

from USGS continuous monitoring data (fig. 10). Although 

the exact values are not always the same, the shapes of the 

measured and metabolism-model simulated daily curves tend 

to match well in all 3 years; the simulated curves for the three 

model versions (G ′

1

, G ′

2

, G ′

3

) lie on top of each other when 

plotted. The daily metabolism models were able to capture 

decreases in chlorophyll a and pH associated with the midseason “crash” in all years (2006–08). 

Total phosphorus, in parts per billion Chlorophyll a, in parts per billion pH, in standard units

A. Total phosphorus

B. Chlorophyll a

C. pH

0

200

400

600

800

1,000

0

100

200

300

300

7.5

7.0

8.0

8.5

9.0

9.5

10.0

10.5

 G’1

 G’2

 G’3

EXPLANATION

Model version

F’O

F’E

Measured

2006 2007 2008

Year

A M J J A S O DN J F M M J J A S O DN A J F M M J J A S O A N

Figure 9. Simulated and measured water-column (A) total phosphorus (as P) concentration, (B) chlorophyll a

concentration, and (C) pH, Upper Klamath Lake, Oregon, 2006–08.
Phase 2 Model Updates 43

tac15-0975_fig10

June July Aug. Sept. Oct.

May June July Aug. Sept. Oct.

June July Aug. Sept. Oct. Nov.

7.5

7.0

8.0

8.5

9.0

9.5

10.0

10.5

7.5

7.0

8.0

8.5

9.0

9.5

10.0

10.5

A. 2006

pH, in standard units

B. 2007

C. 2008

EXPLANATION

Continuous monitor

Model version

 F’O

 F’E

 G’1

 G’2

 G’3

Measured

Figure 10. Simulated and measured water-column pH in (A) 2006, (B) 2007, and (C) 2008, Upper 

Klamath Lake, Oregon.
44 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Prediction of Upper Klamath Lake Metabolism from 

Meteorological Variables, 2006–08 

The best regression model for the seasonal component 

of production was determined from the 14-day dataset to 

include the natural logarithm of chlorophyll a concentration 

and solar radiation, and the best regression model for the 

seasonal component of CR included the natural logarithm of 

chlorophyll a concentration and water temperature (table 13; 

equations 1 and 2). The best model for production did not 

include water temperature or nonalgal phosphorus. This is 

not an indication that production does not depend on these 

variables; rather, there was collinearity among chlorophyll a, 

solar radiation, temperature, and nonalgal phosphorus, and 

not enough independent information contained in the 14-day 

dataset to support inclusion of more than two variables in the 

regression (fig. 11; time series of seasonal overlaid on time 

series of daily, with 14-day points as symbols). 

The mid-season segment of the weekly components 

of CR and total production started on June 23, June 30, and 

June 27 in 2006, 2007, and 2008, respectively. The mid-season 

segment ended on August 29, September 3, and August 30 

in 2006, 2007, and 2008, respectively (fig. 12; time series of 

weekly components with simulated values overlaid and with 

the boundaries between early, mid-, and late-season segments 

indicated). The fit of the models to the weekly component 

generally was poor, but best for the late-season segment. 

It was not possible (based on the performance statistics) to 

determine a minimally acceptable weekly component model 

May July Sept. Nov. Jan. Mar. May July Sept. Nov. Jan. Mar. May July Sept. Nov.

0

10

-10

20

30

40

50

60

0

10

-10

20

30

40

50

60

Grams oxygen per square meter per day

2006 2007 2008

A. Respiration

B. Production

EXPLANATION

Calculated value from measurements

Simulated value, seasonal component

Biweekly simulated value, seasonal component

Figure 11. Measured values, calculated using dissolved oxygen data from nine, continuous, shallowsites, U.S. Geological Survey monitors (table 12), and simulated values of the seasonal component of the 

regression model of metabolism, Upper Klamath Lake, Oregon, 2006–08.
Phase 2 Model Updates 45

tac15-0975_fig12

Mar. May July Sept. Nov. Jan. Mar. May July Sept. Nov. Jan. Mar. May July Sept.

Grams oxygen per square meter per day

2006 2007 2008

A. Respiration

B. Production

EXPLANATION

Calculated value from measurements

Simulated value, weekly component

−10

−20

10

20

30

−10

−20

10

20

30

Early Middle Late Early Middle Late Early Middle Late

Early Middle Late Early Middle Late Early Middle Late

0

0

Figure 12. Respiration and production values calculated from measured data and simulated values of the 

weekly component of the regression model of metabolism, Upper Klamath Lake, Oregon, 2006–08.

for the early-season segment; therefore, no coefficients are 

reported in table 13 and only the seasonal component is shown 

in figure 13 (time series of CR and total production, with 

simulated values, and the boundaries between early, mid-, and 

late-season segments indicated). The coefficients in table 13

indicate that (1) mid-season GPPW and CRW are negatively 

related to wind speed, (2) CRW is positively related to water 

temperature, and (3) GPPW is positively related to shortwave 

solar radiation. Late-season relations show both GPPW and 

CRW positively related to water temperature and stratification, 

and GPPW positively related to shortwave solar radiation. 

When the seasonal and weekly models for production 

were combined (only the seasonal component for the early 

season segment), the fit statistics indicated that the combined 

model explained about 60–80 percent (R2 about 0.60–0.80) 

of the variance over the 2006–08 period, with a minimum 

of 62 percent in 2006 and a maximum of 79 percent in 

2007 (the year in which measurements were used for 

calibrating the weekly component model; table 14). When 

the seasonal and weekly models for CR were combined, the 

fit statistics indicated that the combined model explained 

about 50–70 percent (R2 about 0.50–0.70) of the variance 

over the 2006–08 period, with a minimum of 53 percent in 

2006 and a maximum of 70 percent in 2007. The combined 

production model was more successful than the CR model, as 

indicated by higher R2 values, values of the NS statistic closer 

to 1, and smaller biases. Visual inspection of the time series 

shows that the direction of change of the simulated values 

often is the same as the calculated values in the mid- and 

late-season segments, on time scales of a few days to weeks, 

even though the magnitude of the variability in the simulated 

values is not as great as in the calculated values (fig. 13). 

Model performance is similar to results from Coloso and 

others (2011), whose study attempted multivariate regression 

analysis of GPP and CR. Their R2 values for GPP and CR 

ranged from 0.22 to 0.43 and from 0.18 to 0.55, respectively, 

across 2 years and two lakes.

NPP is obtained by subtracting the CR from production 

(fig. 14). Negative NPP is indicative of declining bloom 

conditions. The daily time series of NPP calculated from 

continuous sonde data shows mid-season decreases in 

NPP, where mid-season was delineated by Julian days 

174–241 in 2006, 181–246 in 2007, and 179–243 in 2008. 
46 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

tac15-0975_fig13

May July Sept. Nov. Jan. Mar. May July Sept. Nov. Jan. Mar. May July Sept. Oct.

Grams oxygen per square meter per day

2006 2007 2008

A. Respiration

B. Production

EXPLANATION

Combined weekly and seasonal simulated values

−10

0

10

20

30

40

50

60

−10

0

10

20

30

40

50

60

Early Middle Late Early Middle Late Early Middle Late

Early Middle Late Early Middle Late Early Middle Late

Daily values of production/community respiration 

calculated from measurements

Figure 13. Respiration and production values calculated from measured data and simulated values from 

the combined seasonal and weekly components of the regression model of metabolism, Upper Klamath Lake, 

Oregon, 2006–08.

Table 14. Goodness-of-fit statistics for the combined (seasonal plus weekly component) 

regression models for total production and community respiration, Upper Klamath Lake, 

Oregon, 2006–08.

[R2, coefficient of determination; RMSE, root mean square error; NS, Nash‐Sutcliffe statistic; (g O2/ m2)/d, 

grams oxygen per square meter per day]

Year

Community respiration Production

Bias NS RMSE

[(g O2

/m2))/d]

R2

[(g O2/m2))/d] Bias NS RMSE[(g O2/m2))/d]

R2

[(g O2/m2))/d]

2006 3.30 0.14 5.62 0.53 3.24 0.37 6.02 0.62

2007 0.26 0.70 4.28 0.70 0.31 0.79 4.71 0.79

2008 2.65 0.27 5.61 0.58 1.28 0.57 5.31 0.69
Phase 2 Model Updates 47

tac15-0975_fig14

0

10

-10

20

30 A. 2006

0

100

200

300

0

10

-10

20

30

B. 2007

0

100

200

300

EXPLANATION

0

10

20

30

C. 2008

0

100

200

300

Calculated net primary production

Simulated net primary production

USGS measured chlorophyll a

Klamath Tribes measured chlorophyll a

Net primary production, in grams oxygen per square meter per day

Chlorophyll a, in parts per billion

150 200 250 300

Julian day

Figure 14. Comparison of the calculated and simulated net primary production to measured chlorophyll a values 

from Klamath Tribes and U.S. Geological Survey monitoring programs, Upper Klamath Lake, Oregon, 2006–08.
48 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Negative values of NPP in late July 2006 and 

mid-August 2008 are consistent with low chlorophyll a

concentrations near those dates, indicative of a severe bloom 

decline (Lindenberg and others 2009, Kannarr and others 

2010). The mid-season decreases in 2007 occurred about 

2 weeks earlier than in 2006 and did not reach negative 

values, consistent with an earlier and milder bloom decline 

in 2007 relative to 2006 and 2008, based on water-quality 

indicators (Kannarr and others, 2010). If the production and 

CR calculations are to be useful for modeling purposes, 

then they must be able to reproduce mid-season minima 

and negative values of NPP. The simulated values of NPP 

reproduce the timing of the mid-season minima in 2007 and 

2008, but the timing in 2006 is about 2 weeks late (fig. 14). 

Nonetheless, because simulated NPP is the difference of two 

simulated quantities, the occurrence of low and negative 

values in the time series of simulated NPP during the midseason is an important result.

Discussion and Suggestions for 

Model Improvements

Model Updates

Continuously Stirred Tank Reactor Model 

Calibration

The continuously stirred tank reactor (CSTR) model 

parameter calibrations, as performed over the original, 

shorter 1991–98 and longer, extended 1991–2010 datasets, 

provided greater insight into the sensitivity of model 

parameters to input data, calibration technique, and model 

structure. With regard to the calibration process in its current 

form, the steps were sequenced so that settling velocity is 

evaluated before recycling and the relationship between these 

parameters becomes apparent. For all model versions that 

were comparable with respect to initial sediment reservoir 

and recycling approach, but had varying light limitation and 

chlorophyll a model structures (versions O′, A', B', E', F'O

and F ′

E), there was a positive correlation between ualg (net 

settling velocity for algal phosphorus) and KR (maximum total 

phosphorus recycle rate). 

The assumptions of the chlorophyll a model also were 

important for model calibration. When equilibrium was 

assumed for the chlorophyll a model, the resulting pH* values 

always were greater than 8.6 (versions O ′, A', B' and C'); 

however, when a Runge-Kutta approach was used to simulate 

chlorophyll a, the pH* values always were less than or equal 

to 8.3. This result is important for model simulations, because 

once recycling starts, the feedback in the model leads to an 

increase in chlorophyll a (a bloom event). If the model is 

calibrated so that recycling occurs at lower pH values, then 

a bloom year is more likely to occur even if the resulting 

maximum chlorophyll a of that bloom year is low. Having 

a low pH* value makes “no bloom” years much less likely, 

regardless of recycling rate. This result is important because 

the algal submodel that uses the Runge-Kutta solution method 

is better calibrated with low pH* values; therefore, the algal 

submodel using the Runge-Kutta solution method is less likely 

to result in “no bloom” years.

Metabolism as an Independent Predictor of 

Bloom Dynamics

The lake metabolism calculations provide insight into 

the typical seasonal bloom dynamics in Upper Klamath 

Lake (UKL). Net primary production (NPP) is the difference 

between gross primary production (GPP) and community 

respiration (CR), and during most of the bloom season in this 

lake NPP is positive, indicating the system is net autotrophic. 

However, NPP is the small difference of two large numbers, 

and provides a more comprehensive insight into the nature 

of what constitutes the bloom “crash,” a rapid decline in 

chlorophyll a accompanied by low dissolved oxygen (DO) 

concentration and reduced pH that occurs in most years 

sometime between mid-July and mid-August. As long as 

GPP exceeds CR, the bloom expands, but when CR exceeds 

GPP, the bloom contracts. Even if the contributor(s) to bloom 

decline are unknown (cyanophages, bacterial infection, 

ultraviolet radiation damage, zooplankton grazing), the decline 

is manifested in the balance between GPP and CR; therefore, 

the decline can be described and potentially predicted as a 

function of external, measurable variables. Hypotheses might 

be that (1) there is a relation between bacterial abundance 

and temperature, or (2) a relation between the probability of 

cyanophage infection and the water-column mixing regime, 

or (3) a relation between ultraviolet radiation exposure 

and stratification. It might not be necessary to completely 

understand the mechanism behind the bloom decline in order 

to build an empirically based predictive model for it, although 

the results from this study were not definitive.

An assumption inherent in substituting CR for the 

growth term in the CSTR model is that autotrophic respiration 

dominates and heterotrophic respiration is a small contributor. 

This probably is a good assumption during peak bloom 

growth, but it is less clear how good an assumption it is 

when the bloom is in decline. The CR estimates in this 

study largely followed the same mid-season decline cycle 

as the GPP estimates (fig. 11). We hypothesize, therefore, 

that CR is dominated by autotrophic respiration, because 

if heterotrophic respiration dominated, we would expect to 

observe an increase in CR coincident with a decrease in GPP 

during bloom decline. Measurements of GPP and CR from 

bottles (Lindenberg and others, 2009) similarly showed that, 

with a few exceptions, CR and GPP tracked with one another. 

There was at least one indication of high heterotrophic 

CR during the bloom decline at a shallow site in Howard 

Bay, but that example also underscored the fact that such 

heterotrophy was not observed at other sites, at least down to 

the depth of the deepest light/dark bottles (1.5 meters). In this 
Discussion and Suggestions for Model Improvements 49

proof-of-concept, CR also is likely underestimated because 

the model did not include respiration below the diurnal mixed 

layer. This could be one reason that our estimates of negative 

NPP (heterotrophy) do not coincide in time perfectly with our 

other measures of bloom decline. 

Our regression results were moderately successful. In a 

previous study by Coloso and others (2011), similar attempts 

were made at multivariate regression analysis of estimates of 

GPP and CR, and the results were not better than the results 

of this study; the coefficient of determination (R2) values for 

GPP and CR ranged from 0.22 to 0.43 and from 0.18 to 0.55, 

respectively, across 2 years and two lakes. Autoregressive 

terms were included, which is an improvement that could be 

made; however, no attempt was made to calibrate in 1 year and 

validate in another, which was an appropriate way of testing 

the predictive capability of our model, for the purpose of 

substituting it into the CSTR model. We obtained substantially 

better results during our calibration year than during our 

validation years.

Incorporation of metabolism to the CSTR model was 

examined as a proof-of-concept by using GPP and CR 

calculated from measured DO data. Model performance 

statistics were comparable to those found for the recalibrated 

CSTR model over the period of 2006–08. The proof-ofconcept showed the efficacy of the metabolism approach 

in capturing mid-season bloom declines, particularly when 

running the model at the daily time step. Incorporation of 

metabolism in the CSTR model also produced a change in 

the timing of the bloom event as compared to the original 

recalibrated CSTR, which lagged by 2–3 weeks. This occurred 

because of the CSTR algal submodel reliance on a growth 

equation that was temperature-limited and that did not allow 

growth at lake temperatures of less than 14 °C. At the end of 

the season, when temperatures decreased to less than 14° C, 

chlorophyll a concentrations were approximately halved, due 

to averaging, until they decreased to less than the acceptable 

minimum threshold of 8 parts per billion. The models that 

incorporated calculated metabolism were not limited by 

temperature in this way. This demonstrates a need to better 

understand the dynamics between lake temperature and 

chlorophyll a concentrations.

Future Model Needs

After combining all the results of a technical review 

of the total maximum daily load (TMDL) model that was 

adopted in 2001, we suggest that five steps be taken toward a 

new model that preserves the benefits of the CSTR approach 

(primarily computational speed, and relative simplicity) 

while providing improved capability for simulating long-term 

response to management scenarios. Some of our findings are 

derived directly from the Phase 1 effort (Wood and others, 

2013), and others follow from the modifications examined 

in this report. These steps are presented in order of least 

difficult to most difficult to implement, given our current 

understanding: 

1. Use a daily time step. As currently scripted with the 

Runge-Kutta approach, the CSTR model runs fast 

enough to simulate many years with a daily time step 

in a reasonable amount of computation time. A daily 

time step is easy to implement and is reasonable in the 

context of the following suggestions that make use of 

data collected continuously (hourly).

2. Use a spatially explicit model as needed to interpolate 

between limited sampling sites and to accurately 

calculate (1) the lakewide average of variables as a 

function of the values at sampling sites, and (2) the 

concentrations in the outflow of the lake, so as to 

accurately calculate the load sent downstream as a 

function of the lakewide average.

3. Substitute GPP and CR for algal growth and respiration 

terms, respectively, in the model. We described the 

errors associated with the scaling-up process and the 

analytical framework for this substitution in section, 

“Lake Metabolism as a Substitute for Algal Growth and 

Respiration,” and we have demonstrated the proof-ofconcept for using continuously measured DO to calculate 

the lake metabolism terms. Improvements should be 

made, the most important being that our estimate of CR 

is low, as we limited the analysis to shallow sites within 

the diurnal mixed layer in order to avoid (for now) the 

problem of mixing from below. Much of the CR and, 

therefore, oxygen consumption in the lake occurs in the 

deepest part of the lake, in the trench along the western 

shoreline (Wood and others, 2006). The assumption that 

CR is dominated by autotrophic respiration could be 

investigated. Further work needs to be done to improve 

the empirical models that enable the prediction of GPP 

and CR from explanatory variables (meteorology and 

water quality) before these variables are substituted for 

the algal growth and respiration terms.

4. Expand equations used to model internal phosphorus 

loading. In our Phase 1 study, we determined that the 

value of the sediment phosphorus recycling rate, which 

currently is treated as a calibration parameter, exerts 

a strong effect on simulated long-term water-column 

concentrations of phosphorus, chlorophyll a, and pH, 

and on sediment phosphorus concentrations. Therefore, 

reducing the uncertainty of this parameter value in the 

model has important implications for refining long-term 

external nutrient loading strategies. This uncertainty 

can be reduced by creating explicit, separate terms 

for benthic invertebrate excretion, groundwater flux, 

diffusive flux, desorption at high pH, and rapid bacterial 

recycling for modeling internal phosphorus loading, 

and by collecting data to parameterize these terms 

where possible. Estimates are available for most of 

these terms, and work is ongoing to estimate benthic 

invertebrate excretion (Wood and others, 2013).
50 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

 To support this comprehensive description of benthic 

fluxes, sediment phosphorus would be best separated 

into two compartments, organic and inorganic, with 

benthic invertebrate excretion and bacterial recycling in 

the organic compartment (also the compartment where 

direct exchange with the water column occurs) and the 

rest in the inorganic compartment. A slow conversion 

between organic and inorganic phosphorus indicating 

sediment diagenesis likely would be needed.

5. Improve the model treatment of net settling velocity. 

As explained for the sediment phosphorus recycling 

rate, the algal settling velocity is a critical parameter for 

simulated water-column and sediment concentrations 

over long time frames, with important implications 

for external loading reduction strategies. However, the 

treatment of net settling velocity as a single calibration 

parameter is not satisfactory for the purposes of 

modeling the dramatic cyanobacterial declines observed 

in the middle of the summer season because the net 

settling velocity varies in time with the proportion of 

buoyant and non-buoyant colonies in the water column. 

Estimates of suspended-sediment concentrations and 

values of rising and falling velocity of colonies based 

on acoustic backscatter and velocity measurements 

from ADCPs showed that (1) velocities, both rising 

and falling, are high and that colonies can traverse the 

water column on sub-daily time scales, and (2) the 

distribution between buoyant and non-buoyant colonies 

varies greatly spatially, vertically within the water 

column, and temporally, particularly in response to the 

amount of thermal stability (stratification) in the water 

column, a function of wind and air temperature (Wood 

and Gartner, 2010). Rapid response to wind-induced 

changes in water-column stratification of large, buoyant, 

cyanobacterial colonies has been observed elsewhere 

(Moreno-Ostos and others, 2009). Because stratification 

varies on a daily basis and with weather patterns, 

the treatment of net settling velocity as an invariant 

calibration parameter would be inadequate.

 Two complementary approaches can better estimate 

net settling velocity. The first approach is similar to the 

approach to growth and respiration in that it avoids the 

need to quantify individual factors and relies instead on 

using measurements to estimate the term as a whole. 

This approach relies on using continuous monitoring 

data to estimate daily depth-integrated mass values at a 

monitoring site (Eldridge and others, 2014). Estimates 

at individual continuous monitoring sites could be 

combined (perhaps with the aid of a spatially explicit 

model; see step 2 to create a lakewide average). Daily 

estimates of lakewide mass storage in combination with 

daily values of NPP allow the daily settling flux to be 

calculated by difference. This approach only can be 

pursued in combination with the calculation of GPP and 

CR (step 3).

The second approach relies on field work to measure 

rising and falling velocities and to estimate the 

fraction of rising and falling colonies as a function of 

environmental variables. Some of this work has been 

done with acoustic instruments (Wood and Gartner, 

2010), but more acoustic measurements are available and 

could be analyzed for information. Sampling equipment 

could be fabricated to capture both rising and falling 

colonies over short periods of time in the water column 

and thereby measure velocities and relative amounts of 

buoyant and non-buoyant colonies.

 After measurements are obtained, this approach still 

requires modeling (both one-dimensional water-column 

models that investigate only the interaction between 

buoyancy and stratification, and three-dimensional 

models that add horizontal transport to the mix) to 

assist the scaling-up process. A modeling effort would 

be greatly facilitated by a study targeted at quantifying 

the rates of creation and destruction of the carbohydrate 

ballast that the Aphanizomenon flos‐aquae colonies use 

to regulate their buoyancy. Although there are a few 

different mechanisms that colonies use (construction of 

gas vesicles, turgor pressure collapse of vesicles, and 

ballast), it seems that the mechanism that works most 

effectively on the time scales of most interest (hours, 

because velocities indicate colonies can traverse the 

water column in hours over most of the lake) is the 

construction and destruction of carbohydrate ballast 

(Konopka and others, 1987; Whitton and Potts, 2000; 

Porat and others, 2001).

The combined results of Phases 1 and 2 of the TMDL 

model review have indicated strengths and limitations of the 

UKL TMDL model. Translating the model into the R language 

and adopting a numerical method for integrating the model 

equations has facilitated recalibration and evaluation of the 

effect of the different model terms on the resulting predictions 

of water quality. Changes were made to several model 

algorithms to make them consistent with current research on 

lake processes elsewhere. Sensitivity analysis and long-term 

simulations produced the calibration parameters that are both 

influential in determining model outcome while also having 

a high degree of uncertainty, suggesting areas for future 

research to improve understanding of the effects of changes 

in nutrient management strategies. Finally, we have suggested 

several steps to combine these findings in order to create a 

more robust model in the future and to better understand the 

likely results to changes in lake and nutrient management. 

Given the importance of water quality and nutrient status of 

Upper Klamath Lake to endangered fish and to downstream 

water quality in the Klamath River and the lower Klamath 

River Basin, finding ways to implement the suggested model 

refinements may help generate improvements in water quality 

in the Klamath River Basin.
References Cited 51

References Cited 

Antenucci, J.P., Tan, K.M., Eikaas, H.S., and Imberger, J., 

2013, The importance of transport processes and spatial 

gradients on in situ estimates of lake metabolism: 

Hydrobiologia, v. 700, p. 9–12, doi:10.1007/s10750-012-

1212-z.

Aristegi, L., Izagirre, O., and Elosegi, A., 2009, Comparison 

of several methods to calculate reaeration in streams, and 

their effects on estimation of metabolism: Hydrobiologia, 

v. 635, p. 113–124.

Bales, J.D., and Nardi, M.R., 2007, Automated routines 

for calculating whole-stream metabolism—Theoretical 

background and user's guide: U.S. Geological Survey 

Techniques and Methods, book 4, chap. C2, 33 p., http://

pubs.usgs.gov/tm/tm4c2/.

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. 

Burdick, S.M., VanderKooi, S.P., and Anderson, G.O., 2009, 

Spring and summer spatial distribution of endangered 

juvenile Lost River and shortnose suckers in relation 

to environmental variables in Upper Klamath Lake, 

Oregon—2007 Annual Report: U.S. Geological Survey 

Open-File Report 2009-1043, 57 p., http://pubs.usgs.gov/

of/2009/1043/.

Burnham, K.P., and Anderson, D.R., 2002, Model selection 

and multimodel inference—A practical informationtheoretic approach (2d ed.): New York, Springer-Verlag.

Chapra, S.C., 1997, Surface water-quality modeling: Long 

Grove, Illinois, Waveland Press.

Chapra, S.C., and Di Toro, D.M., 1991, Delta method for 

estimating primary production, respiration, and reaeration 

in streams: Journal of Environmental Engineering, v. 117, 

no. 5, p. 640–655. 

Coloso, J.J., Cole, J.J., and Pace, M.L., 2011, Difficulty in 

discerning drivers of lake ecosystem metabolism with 

high-frequency data, Ecosystems, v. 14, p. 935–948, 

doi:10.1007/s10021-011-9455-5.

Devine, J.A., and Vanni, M.J., 2002, Spatial and seasonal 

variation in nutrient excretion by benthic invertebrates 

in a utrophic reservoir: Freshwater Biology, v. 47, 

p. 1,107–1,121.

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, no. 1, p. 7, http://

dx.doi.org/10.1023/B:HYDR.0000027718.95901.ae.

Eilers, J.M., Kann, J., Cornett, J., Moser, K., St. Amand, A., 

and Gubala, C.P., 2001, Recent paleolimnology of Upper 

Klamath Lake: Bureau of Reclamation, 44 p.

Eldridge, S.L.C., Wherry, S.A., and Wood, T.M., 2014, 

Statistical analysis of the water-quality monitoring 

program, Upper Klamath Lake, Oregon, and optimization 

of the program for 2013 and beyond: U.S. Geological 

Survey Open-File Report 2014-1009, 82 p., http://dx.doi.

org/10.3133/ofr20141009.

Eldridge, S.L.C., Wood, T.M., and Echols, K.R., 2012, Spatial 

and temporal dynamics of cyanotoxins and their relation 

to other water quality variables in Upper Klamath Lake, 

Oregon, 2007–09: U.S. Geological Survey Scientific 

Investigations Report 2012-5069, 32 p. plus appendixes, 

http://pubs.usgs.gov/sir/2012/5069/.

Environmental and Water Resources Institute, 2005, The 

ASCE Standardized Reference Evapotranspiration 

Equation: ASCE-EWRI Task Committee Report, 59 p.

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.

Gelda, R.K., and Effler, S.W., 2002, Metabolic rate estimates 

for a eutrophic lake from diel dissolved oxygen signals: 

Hydrobiologia, v. 485, nos. 1–3, p. 51–66, http://dx.doi.

org/10.1023/A%3A1021327610570. 

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.

Kann, J., and Walker, W.W., 1999, Nutrient and hydrological 

loading to Upper Klamath Lake, Oregon, 1991–1998: 

Klamath Tribes, 48 p. plus appendixes.
52 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon

Kann, J., 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.

Kann, J., 2011, Upper Klamath Lake 2010 data summary 

report: Klamath Tribes, 43 p. plus appendixes.

Kannarr, K.E., Tanner, D.Q., Lindenberg, M.K., and 

Wood, T.M., 2010, Water-quality data from Upper Klamath 

and Agency Lakes, Oregon, 2007–08: U.S. Geological 

Survey Open-File Report 2010-1073, 28 p., http://pubs.

usgs.gov/of/2010/1073/.

Konopka, A., Kromkamp, J.C., and Mur, L.R., 1987, 

Regulation of gas vesicle content and buoyancy in light- or 

phosphate-limited cultures of Aphanizomenon flos-aquae

(Cyanophyta): Journal of Phycology, v. 23, no. 1, p. 70–78.

Lindenberg, M.K., Hoilman, G., and Wood, T.M., 2009, Water 

quality conditions in Upper Klamath and Agency Lakes, 

Oregon, 2006: U.S. Geological Survey Open-File Report 

2008-5201, 54 p.

Martin, J.L., and McCutcheon, S.C., 1999, Hydrodynamics 

and transport for water quality modeling: Boca Raton, Fla., 

Lewis Publishers.

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. 

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

Moreno-Ostos, E., Cruz-Pizarro, L., Basanta, A., and 

George, D.G., 2009, The influence of wind-induced 

mixing on the vertical distribution of buoyant and sinking 

phytoplankton species: Aquatic Ecology, v. 43, p. 271–284.

Porat, R., Teltsch, B., Perelman, A., and Dubinsky, Z., 

2001, Diel buoyancy changes by the cyanobacterium 

Aphanizomenon ovalisporum from a shallow reservoir: 

Journal of Plankton Research, v. 23, no. 7, p. 753–763, 

doi:10.1093/plankt/23.7.753, accessed November 23, 2010, 

at http://plankt.oxfordjournals.org/content/23/7/753.full.

Odum, H.T., 1956, Primary production in lowing waters: 

Limnology and Oceanography, v. 1, no. 2, p. 102–117.

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

http://www.deq.state.or.us/WQ/TMDLs/klamath.htm.

R Core Team, 2014, R—A language and environment for 

statistical computing: Vienna, Austria, R Foundation for 

Statistical Computing, http://www.R-project.org/.

Raymond, P.A., Zappa, C.J., Butman, D., Bott, T.L., Potter, J., 

Mulholland, P., Laursen, A.E., McDowell, W.H., and 

Newbold, D., 2012, Scaling the gas transfer velocity and 

hydraulic geometry in streams and small rivers: Limnology 

and Oceanography—Fluids and Environments, v. 2, 

p. 41–53.

Simon, N., and Ingle, S., 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, 46 p., http://pubs.usgs.gov/of/2011/1168/.

Simon, N., Lynch, D., and Gallaher, T., 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, no. 1, p. 139–153, http://dx.doi.org/10.1007/s11270-

009-0033-9.

Staehr, P.A., Bade, D., Van de Bogert, M.C., Koch, G.R., 

Williamson, C., Hanson, P. Cole, J.J., and Kratz, T., 2010, 

Lake metabolism and the diel oxygen technique—State of 

the science: Limnology and Oceanography—Methods, v. 8, 

p. 628–644.

Sullivan, A.B., Rounds, S.A., Deas, M.L., Asbill, J.R., 

Wellman, R.E., Stewart, M.A., Johnston, M.W., and 

Sogutlugil, I.E., 2011, Modeling hydrodynamics, water 

temperature, and water quality in the Klamath River 

upstream of Keno Dam, Oregon, 2006–09: U.S. Geological 

Survey Scientific Investigations Report 2011-5105, 70 p., 

http://pubs.usgs.gov/sir/2011/5105/.
References Cited 53

University of California, 2014, California weather data, 

Description of TULELK2.A: University of California 

Agriculture and Natural Resources—Statewide Integrated 

Pest Management Program Web site, accessed June 9, 2014, 

at http://169.237.140.1/calludt.cgi/WXDESCRIPTION?MA

P=&STN=TULELK2.A.

U.S. Geological Survey, 2011, Change to solubility equations 

for oxygen in water: U.S. Geological Survey Office of 

Water Quality Technical Memorandum 2011.03, 13 p.

Van de Bogert, M.C., Carpenter, S.R., Cole, J.J., and 

Pace, M.L., 2007, Assessing pelagic and benthic 

metabolism using free water measurements: Limnology and 

Oceanography—Methods, v. 5, p. 145–155.

VanderKooi, S.P., Burdick, S.M., Echols, K.R., Ottinger, 

C.A., Rosen, B.H., and Wood, T.M., 2010, Algal toxins in 

Upper Klamath Lake, Oregon—Linking water quality to 

juvenile sucker health: U.S. Geological Survey Fact Sheet 

2009-3111, 2 p., http://pubs.usgs.gov/fs/2009/3111/.

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/klamathbasin/ukldrainage/devphostmdl.pdf. 

Walker, W.W., Walker, J.D., and Kann, J., 2012, Evaluation 

oof water and nutrient balances for the Upper Klamath Lake 

Basin in water years 1992–2010: Technical Report to the 

Klamath Tribes Natural Resources Department, 50 p. plus 

appendixes. 

Whitton, B.A., and Potts, M., eds., 2000, The ecology 

of cyanobacteria—Their diversity in time and space: 

Dordrecht, The Netherlands, Kluwer Academic Publishers, 

669 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 Water-Resources Investigations Report 

96-4079, 57 p., http://pubs.usgs.gov/wri/1996/4079/report.

pdf.

Wood, T.M., and Gartner, J.W., 2010, Use of acoustic 

backscatter and vertical velocity to estimate concentration 

and dynamics of suspended solids in Upper Klamath Lake, 

south-central Oregon—Implications for Aphanizomenon 

flos-aquae: U.S. Geological Survey Scientific Investigations 

Report 2010-5203, 20 p.

Wood, T.M., Hoilman, G.R., and Lindenberg, M.K., 2006, 

Water-quality conditions in Upper Klamath Lake, 

Oregon, 2002–04: U. S. Geological Survey Scientific 

Investigations Report 2006-5209, 52 p., http://pubs.usgs.

gov/sir/2006/5209/.

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 OpenFile Report 2013-1262, 69 p. plus appendix, http://pubs.

usgs.gov/of/2013/1262/.
54 Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon
Appendix A 55

Appendix A. Derivation of Net Settling Velocity

Cyanobacterial colonies in Upper Klamath Lake, Oregon, 

probably lie along a continuum between the heaviest colonies 

that fall with a maximum velocity and the lightest colonies 

that rise with a maximum velocity, including colonies of 

neutral buoyancy. To model these dynamics, we assume that 

the cyanobacterial colonies in the lake can be divided into two 

groups—rising colonies and falling colonies. We also assume 

that rising colonies are characterized by a single, constant 

rising velocity ur and falling colonies are characterized by 

a single, constant falling velocity uf (both with a positive 

sign). At any given location in the lake and point in time, 

the chlorophyll a associated with rising colonies is given by 

( ) ( ) ,, * ,, ˆ rise r szt B szt , and the chlorophyll a associated with 

falling colonies is given by ( ) ( ) ,, * ,, ˆ fall r szt B szt , where rrise 

is the fraction of colonies that are rising (buoyant), rfall is the 

fraction of colonies that are falling (negatively buoyant), and 

rrise + rfall = 1. With these assumptions, the equation for the 

mass balance of Bˆ contains a settling term SET of the form:

( ) ˆ ˆ

, , fall rise

r f

r B r B

SET s z t u u z z

∂ ∂ = − ∂ ∂ (A1)

When this term is averaged to scale up, the result is:

SETTLE

D

u r u r r rise f fall z bottom= + { } =

1 (A2)

where the double overbar indicates averaging over area and 

time. Note that we have allowed for rising colonies at the 

bottom; in other words, for colonies that previously settled to 

re-enter the water column. In the current version of the CSTR 

model, this is approximated as:

SETTLE

D

u r u r

D

u r rise f fall z bottom~ alg

1 1

+ ( ) = =

(A3)

and ualg is the weighted average near-bottom velocity of the 

buoyant and non-buoyant colonies together, which effectively 

is a net settling velocity. To make the modeling of net settling 

velocity more manageable, the term was further simplified by 

making ualg a calibration parameter.

Publishing support provided by the U.S. Geological Survey

 Science Publishing Network, Tacoma Publishing Service Center

For more information concerning the research in this report, contact the

Director, Oregon Water Science Center

U.S. Geological Survey

2130 SW 5th Avenue

Portland, Oregon 97201

http://or.water.usgs.gov
Wherry and others—Revision and Proposed Modification of a Total Maximum Daily Load Model for Upper Klamath Lake, Oregon—Scientific Investigations Report 2015–5041

ISSN 2328-0328 (online)

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