Modeling Hydrodynamics and Heat Transport in 

Upper Klamath Lake, Oregon, and Implications for 

Water Quality

Prepared in cooperation with the Bureau of Reclamation

Scientific Investigations Report 2008–5076

U.S. Department of the Interior

U.S. Geological Survey
Cover: Photograph of Upper Klamath Lake, Oregon, with pelican. 

(Photograph taken by Lee Simons, U.S. Geological Survey, 2007.)
Modeling Hydrodynamics and Heat Transport 

in Upper Klamath Lake, Oregon, and 

Implications for Water Quality

By Tamara M. Wood, Ralph T. Cheng, Jeffrey W. Gartner, Gene R. Hoilman, Mary K. 

Lindenberg, and Roy E. Wellman

Prepared in cooperation with the Bureau of Reclamation

Scientific Investigations Report 2008–5076

U.S. Department of the Interior

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

DIRK KEMPTHORNE, Secretary

U.S. Geological Survey

Mark D. Myers, Director

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

For product and ordering information:

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

Telephone: 1-888-ASK-USGS

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

natural hazards, and the environment:

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

Telephone: 1-888-ASK-USGS

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

U.S. Government.

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

reproduce any copyrighted materials contained within this report.

Suggested citation:

Wood, T.M., Cheng, R.T., Gartner, J.W., Hoilman, G.R., Lindenberg, M.K., and Wellman, R.E., 2008, Modeling 

hydrodynamics and heat transport in Upper Klamath Lake, Oregon, and implications for water quality: U.S. Geological 

Survey Scientific Investigations Report 2008–5076, 48 p.
iii

Contents

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

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

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

Model Description .........................................................................................................................................4

Governing Equations ............................................................................................................................4

Orthogonal Unstructured Grids .........................................................................................................5

Numerical Approximation ..................................................................................................................6

Turbulence Closure...............................................................................................................................9

Surface and Bottom Boundary Conditions.......................................................................................9

Wind Interpolation..............................................................................................................................10

Heat Transport.....................................................................................................................................10

Datasets.........................................................................................................................................................12

Tributary Inflows and Outflows.........................................................................................................12

Inflows .........................................................................................................................................12

Meteorological Data...........................................................................................................................15

Calibration Data...................................................................................................................................15

Surface Elevation.......................................................................................................................15

Velocity ........................................................................................................................................15

Water Temperature ...................................................................................................................15

Other Datasets.....................................................................................................................................15

Simulation Results ......................................................................................................................................17

Hydrodynamics....................................................................................................................................17

Heat Transport.....................................................................................................................................27

Implications for Water Quality...................................................................................................................37

Acknowledgments.......................................................................................................................................42

Summary........................................................................................................................................................42

References Cited..........................................................................................................................................44

Appendix A. Calculation of the Rate of Change in Dissolved Oxygen Over 24 Hours from 

Light/Dark Bottle Incubations and the Extinction Coefficient.................................................47
iv

Figures

Figure 1. Maps showing location of monitoring sites in Upper Klamath Lake, Oregon ...... 2

Figure 2. Diagram showing an example of an unstructured orthogonal grid .................. 6

Figure 3. Diagram showing unstructured orthogonal grid used in model simulations of 

Upper Klamath Lake, Oregon............................................................... 7

Figure 4. Graphs showing inflows to the Upper Klamath Lake/Agency Lake hydrologic 

system in Oregon in (A) 2005 and (B) 2006................................................ 12

Figure 5. Graphs showing outflows from the Upper Klamath Lake/Agency Lake 

hydrologic system in Oregon in (A) 2005 and (B) 2006 ................................. 14

Figure 6. Graphs showing hourly values of lake elevation at Upper Klamath Lake, 

Oregon, measured at two gages, 2005 ................................................... 19

Figure 7. Graphs showing wind speed and direction at site MDL in Upper Klamath Lake, 

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

Figure 8. Diagram showing simulated depth-averaged currents under prevailing 

(northwest) wind conditions, Upper Klamath Lake, Oregon ........................... 21

Figure 9. Diagram showing simulated depth-averaged currents under reversed 

(southwest) wind conditions, Upper Klamath Lake, Oregon ........................... 21

Figure 10. Graphs showing observed and simulated current speeds at Acoustic Doppler 

Current Profiler (ADCP) sites in Upper Klamath Lake, Oregon, 2005 .................. 22

Figure 11. Graphs showing observed and simulated current direction at Acoustic Doppler 

Current Profiler (ADCP) sites in Upper Klamath Lake, Oregon, 2005 .................. 22

Figure 12. Graphs showing goodness-of-fit velocity statistics at site ADCP1 in Upper 

Klamath Lake, Oregon, 2005 and 2006 ................................................... 23

Figure 13. Graphs showing observed and simulated east-west and north-south current 

components at site ADCP3 in Upper Klamath Lake, Oregon, 2005..................... 23

Figure 14. Graphs showing hourly values of lake elevation at Upper Klamath Lake, 

Oregon, measured at two gages, 2006 ................................................... 25

Figure 15. Graphs showing observed and simulated current speed and direction at 

two Acoustic Doppler Current Profiler (ADCP) sites in Upper Klamath Lake, 

Oregon, 2006 ................................................................................. 26

Figure 16. Graphs showing observed and simulated temperatures at 10 shallow sites in 

Upper Klamath Lake, Oregon, 2005 ...................................................... 28

Figure 17. Graphs showing goodness-of-fit temperature statistics at 13 sites in Upper 

Klamath Lake, Oregon, 2005 and 2006 ................................................... 29

Figure 18. Graphs showing air temperature at site MDL in 2005 and 2006, and air 

temperature at site WMR-MET in 2005, Upper Klamath Lake, Oregon ............... 31

Figure 19. Graphs showing observed and simulated temperatures at three deep sites in 

Upper Klamath Lake, Oregon, 2005 ...................................................... 31

Figure 20. Graphs showing daily maximum and minimum difference in temperature 

between near-surface and near-bottom of the water column at three deep sites 

in Upper Klamath Lake, Oregon, 2005 ................................................... 32

Figure 21. Graphs showing observed and simulated temperatures at four deep sites in 

Upper Klamath Lake, Oregon, 2006 ...................................................... 33

Figure 22. Graphs showing observed and simulated temperatures at eight shallow sites 

in Upper Klamath Lake, Oregon, 2006 ................................................... 34

Figure 23. Graphs showing observed and simulated temperatures at five nearshore sites 

in Upper Klamath Lake, Oregon, 2006 ................................................... 35
v

Figures—Continued

Figure 24. Graphs showing daily maximum and minimum difference in temperature 

between near-surface and near-bottom of the water column at three deep 

sites in Upper Klamath Lake, Oregon, 2006 ............................................. 36

Figure 25. Diagram showing the location of polygons in the numerical grid where the depth 

was greater than 4.5 meters and where tracer concentrations were defined for 

numerical experiments ..................................................................... 37

Figure 26. Diagram showing concentration of tracers (A) T1 and (B) T2 in the surface layer 

of Upper Klamath Lake, Oregon, at day 5 of numerical experiments starting on 

August 1, 2005 .............................................................................. 38

Figure 27. Graphs showing time series of the simulated concentration of tracers T1 and 

T2 during 10 days of numerical experiments starting on August 1, 2005, at sites 

MDN and MDL in Upper Klamath Lake, Oregon ....................................... 39

Figure 28. Graphs showing the 24-hour change in dissolved oxygen concentration that 

would result from the oxygen production and consumption rates measured in 

light and dark bottle incubations at sites (A) MDN, (B) MDT, and (C) RPT, in Upper 

Klamath Lake, Oregon, between June and October, 2006, as a function of the 

chlorophyll a concentration measured at the same site on the same date ......... 40

Figure 29. Graphs showing time series of dissolved oxygen concentration data collected 

at sites MDN, MDT, and MDL in Upper Klamath Lake, Oregon, 2005.................. 40

Figure 30. Graphs showing average hourly vertical velocities at site ADCP1 (A) near 

surface and (B) near bottom, as measured by an Acoustic Doppler Current 

Profiler (ADCP), Upper Klamath Lake, Oregon .......................................... 41

Tables

Table 1. List of symbols used in this report ......................................................... 8

Table 2. Summary of the dates of availability of meteorological, velocity and water 

temperature data in Upper Klamath Lake between May and October of 2005 and 

2006; summary of the simulation periods in 2005 and 2006, and the dates over 

which error statistics were calculated in 2005 and 2006 .............................. 16

Table 3. List of inputs to the hydrodynamic and heat transport model of Upper Klamath 

Lake, Oregon ................................................................................. 18

Table 4. Goodness-of-fit velocity statistics for the UnTRIM model of Upper Klamath Lake, 

Oregon, 2005 ................................................................................. 24

Table 5. Goodness-of-fit velocity statistics for the UnTRIM model of Upper Klamath Lake, 

Oregon, 2006 ................................................................................. 27

Table 6. Goodness-of-fit temperature statistics for the UnTRIM model of Upper Klamath 

Lake, Oregon, 2005 ........................................................................ 29

Table 7. Goodness-of-fit temperature statistics for the UnTRIM model of Upper Klamath 

Lake, Oregon, 2006. ........................................................................ 30
vi

Conversion Factors and Datums

Conversion Factors

Multiply By To obtain

centimeter (cm) 0.3937 inch (in.)

cubic meter (m3

) 1.308 cubic yard (yd3) 

cubic meter per second (m3

/s) 70.07 acre-foot per day (acre-ft/d) 

cubic meter per second (m3

/s) 35.31 cubic foot per second (ft3/s)

cubic meter per second (m3

/s) 22.83 million gallons per day (Mgal/d) 

gram (g) 0.03527 ounce, avoirdupois (oz)

joule (J) 0.0000002 kilowatthour (kWh)

kilogram (kg) 2.205 pound avoirdupois (lb)

kilogram per cubic meter (kg/m3

) 0.06242 pound per cubic foot (lb/ft3) 

kilometer (km) 0.6214 mile (mi)

liter (L) 1.057 quart (qt)

meter (m) 3.2808 feet (ft) 

millibar (mb) 0.02953 inches of mercury at 0°C

square kilometer (km2

) 0.3861 square mile (mi2)

square meter (m2

) 10.76 square foot (ft2) 

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

°F=(1.8×°C)+32.

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

°C=(°F–32)/1.8.

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

Datums

Vertical coordinate information is referenced to the National Geodetic Vertical Datum of 1929 

(NGVD 29). ”Elevation,” as used in this report, refers to distance above the vertical datum.

Horizontal coordinate information is referenced to the North American Datum of 1927 (NAD 27).
Modeling Hydrodynamics and Heat Transport in Upper 

Klamath Lake, Oregon, and Implications for Water Quality

By Tamara M. Wood, Ralph T. Cheng, Jeffrey W. Gartner, Gene R. Hoilman, Mary K. Lindenberg, and 

Roy E. Wellman

Abstract 

The three-dimensional numerical model UnTRIM was 

used to model hydrodynamics and heat transport in Upper 

Klamath Lake, Oregon, between mid-June and mid-September 

in 2005 and between mid-May and mid-October in 2006. 

Data from as many as six meteorological stations were used 

to generate a spatially interpolated wind field to use as a 

forcing function. Solar radiation, air temperature, and relative 

humidity data all were available at one or more sites. In 

general, because the available data for all inflows and outflows 

did not adequately close the water budget as calculated from 

lake elevation and stage-capacity information, a residual 

inflow or outflow was used to assure closure of the water 

budget. 

Data used for calibration in 2005 included lake elevation 

at 3 water-level gages around the lake, water currents at 

5 Acoustic Doppler Current Profiler (ADCP) sites, and 

temperature at 16 water-quality monitoring locations. The 

calibrated model accurately simulated the fluctuations of the 

surface of the lake caused by daily wind patterns. The use of 

a spatially variable surface wind interpolated from two sites 

on the lake and four sites on the shoreline generally resulted 

in more accurate simulation of the currents than the use of a 

spatially invariant surface wind as observed at only one site 

on the lake. The simulation of currents was most accurate at 

the deepest site (ADCP1, where the velocities were highest) 

using a spatially variable surface wind; the mean error (ME) 

and root mean square error (RMSE) for the depth-averaged 

speed over a 37-day simulation from July 26 to August 31, 

2005, were 0.50 centimeter per second (cm/s) and 3.08 cm/s, 

respectively. Simulated currents at the remaining sites were 

less accurate and, in general, underestimated the measured 

currents. The maximum errors in simulated currents were 

at a site near the southern end of the trench at the mouth of 

Howard Bay (ADCP7), where the ME and RMSE in the 

depth-averaged speed were 3.02 and 4.38 cm/s, respectively. 

The range in ME of the temperature simulations over the same 

period was –0.94 to 0.73 degrees Celsius (°C), and the RMSE 

ranged from 0.43 to 1.12°C. The model adequately simulated 

periods of stratification in the deep trench when complete 

mixing did not occur for several days at a time.

The model was validated using boundary conditions and 

forcing functions from 2006 without changing any calibration 

parameters. A spatially variable wind was used. Data for the 

model validation periods in 2006 included lake elevation at 4 

gages around the lake, currents collected at 2 ADCP sites, and 

temperature collected at 21 water-quality monitoring locations. 

Errors generally were larger than in 2005. ME and RMSE 

in the simulated velocity at ADCP1 were 2.30 cm/s and 

3.88 cm/s, respectively, for the same 37-day simulation over 

which errors were computed for 2005. The ME in temperature 

over the same period ranged from –0.56 to 1.5°C and the 

RMSE ranged from 0.41 to 1.86°C. 

Numerical experiments with conservative tracers were 

used to demonstrate the prevailing clockwise circulation 

patterns in the lake, and to show the influence of water from 

the deep trench located along the western shoreline of the lake 

on fish habitat in the northern part of the lake. Because water 

exiting the trench is split into two pathways, the numerical 

experiments indicate that bottom water from the trench has 

a stronger influence on water quality in the northern part of 

the lake, and surface water from the trench has a stronger 

influence on the southern part of the lake. This may be part of 

the explanation for why episodes of low dissolved oxygen tend 

to be more severe in the northern than in the southern part of 

the lake. 

Introduction

Upper Klamath Lake (UKL) is a large (surface area 

232 km2) and shallow (mean depth 2.8 m at full pool) lake 

located in southern Oregon, in the semiarid landscape in the 

rain shadow east of the Cascade Range, at a full-pool elevation 

of 1,262.9 m above the National Geodetic Vertical Datum 

of 1929 (fig. 1). Agency Lake is connected to the northern 

end of Upper Klamath Lake through a narrow channel 

(Agency Straits), and adds about 38 km2

 of surface area. 
2 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality Figure 1. Location of monitoring sites in Upper Klamath Lake, Oregon. UpperKlamath Upper Klamath Lake Lake Aspen Lake Aspen Lake

River

Williamson

Agency

Straits

Agency

Straits

River

Williamson

Howard

Bay

Sevenmile

Canal

Ball

Point

Shoalwater Bay

Wocus

Marsh

Caledonia

Marsh

Pelican

Bay

Link

River Lake

Ewauna

Agency

Lake

Agency

Lake

Moss

Creek

Wood

River

Rock

Creek

Sprague

River

Howard

Bay

Ball

Bay

Shoalwater

Bay

Wocus

Marsh

Caledonia

Marsh

Pelican

Bay

Fourmile

Canal

Link

River Lake

Ewauna

Moss

Creek

Wood

River

Rock

Creek

Sprague

River

Ball

Point

A Canal

Sevenmile

Canal

Fourmile

Canal

A Canal

97

97

Buck

Island

Bare

Island

Bare

Island

140

161

70

66

97

97

140

161

70

66

Buck

Island

ModocPoint

Klamath 

Falls

Klamath 

Falls

Plum Valley

Eagle Ridge

Rattlesnake

Point

Modoc

Point

ModocPoint

Modoc

Point

Rattlesnake

Point

Link River 

Dam

Plum Valley

Eagle Ridge

Link River 

Dam

593

592 591

593

592 591

ADCP7

ADCP6

ADCP5 ADCP3

ADCP1

ADCP9

ADCP1B

ADCP1A

MDL

EHB

RPT

MPT

EPT

FBS

WMR

MDT

MDN

BLB

NBI

HDB

SHB

WMR-MET

HDBMET

BLB-MET

SSHRMET

WMR-MET

HDBMET

BLB-MET

SSR

WDW

SHL

HPK

GBE WMR

RPT

NBI

MRM

MDL

HDB

FBS

EBB

EHB

SET-U

SET-L

MDN-U

MDN-L

MDT-U

MDT-L

EPT-U

EPT-L

SSHRMET

11507500 11507500

11502500

11507000

11505900

11505800

AGKO

11502500

11507000

11505900

11505800

AGKO

2005 2006

122°00’ 121°45'

42°30’

122°00’ 121°45'

42°30’

42°15’

42°15’

0 5 MILES0 5 KILOMETERS

Studyarea

OREGON

EXPLANATION

Depth of Klamath and Agency Lakes,

in feet

ADCP location–2005-2006

0

50

USGS water level gage

Meteorological site

11505800

USGS surface-water gaging site 11502500

MPT

EBB

MDL

MDL

BLB MET

ADCP6

Continuous water-quality monitoring 

site

2005200620052006

Continuous water-quality monitoring 

site, and meteorological site

Graham, Matthews and Associates 

gage

591
Introduction 3

Paleolimnological evidence indicates that the lake has been 

highly productive for at least the last 1,000 years (Eilers and 

others 2004). This can be attributed to, among other things, 

the fact that much of the basin drains volcanic soils with high 

phosphorus content, and, because most of the lake is shallow, 

photosynthetically active radiation penetrates most of its 

volume on a daily basis, providing energy that is converted to 

biomass through the photosynthesis performed by algae and 

cyanobacteria. 

Starting about 150 years ago, changes in land use in the 

UKL basin led to increased sedimentation rates and nutrient 

loads, and to a decrease in the nitrogen to phosphorus ratio 

in the nutrient loads to the lake. Simultaneous with these 

changes, the diverse taxa that characterized the assemblage 

of phytoplankton prior to 150 years ago gradually came to be 

dominated by a single species of buoyant cyanobacterium, 

Aphanizomenon flos aquae (AFA) (Phinney and Peek, 1960; 

Miller and Tash, 1967; Kann, 1998; Eilers and others, 2004; 

Eilers and others, 2001; Bradbury, Colman, and Dean, 2004; 

Bradbury, Colman, and Reynolds, 2004; Colman and others, 

2004; Colman, Bradbury, and Rosenbaum, 2004). The 

shallowness of the lake also makes the transfer of nutrients 

from the sediments to the water column and back again to 

the sediments particularly efficient, although the specific 

mechanism for the transfer is still debated (Kann, 1998; 

Walker, 2001; Fisher and Wood, 2004). Annual blooms of 

AFA are perpetuated primarily by this internal recycling 

of nutrients stored in the sediments, rather than by external 

sources (Kann and Walker, 1999; Walker, 2001; Kann, 1998).

The annual cyanobacterial blooms have important 

implications for water quality. Highly supersaturated 

dissolved oxygen concentrations with large diel swings, and 

high pH values caused by the photosynthetic removal of 

carbon dioxide from the water column, are associated with 

the rapidly expanding phase of the bloom. A rapid bloom 

decline is accompanied by undersaturated dissolved oxygen 

concentrations as photosynthetic production of oxygen slows 

dramatically and ongoing respiratory demands continue, and 

decay processes associated with cell senescence consume 

oxygen rapidly. Severe low dissolved oxygen events (LDOEs) 

associated with bloom declines are detrimental to the survival 

of the shortnose and Lost River suckers, which are listed as 

endangered under the Endangered Species Act by the U.S. 

Fish and Wildlife Service (2001). As determined from an 

analysis of three large fish die-offs during 1995–97 (Perkins 

and others, 2000), severe LDOEs in which dissolved oxygen 

concentrations of less than 4 mg/L occur throughout the water 

column for a large part of the day and persist for several 

days at a time may cause fish death. These conditions may 

also cause fish death indirectly by forcing fish into crowded 

conditions and facilitating the spread of disease among 

animals already weakened by exposure to hypoxia, high 

pH, and the high un-ionized ammonia concentrations that 

sometimes accompany conditions of high pH. 

Continuous water-quality monitors with well-documented 

data quality were first installed by the U.S. Geological 

Survey (USGS) for long-term deployment in 2002. Thus, 

reconstructing the details of the LDOEs of 1995–97 on time 

scales of a few days to a week is not possible, but information 

gathered since 2002 indicates that there is a rapid decline 

in the AFA bloom in most years around the end of July or 

beginning of August (Wood and others, 2006; Hoilman 

and others, 2008). The longest dataset of water-quality 

measurements in UKL is a 17-year record of biweekly profiles 

of conventional water-quality variables, as well as depthaveraged chlorophyll a and nutrient concentrations (Kris 

Fisher, Klamath Tribes, oral commun., 2007). This dataset 

does not provide enough temporal resolution to determine 

precise dates of bloom peaks and troughs, but the general 

pattern appears consistent over this long record. The effect 

of the late July to early August decline on dissolved oxygen 

concentrations over a large area in the northern part of the 

lake, which is the preferred habitat for adult suckers, varies 

in severity from year to year. In 2003, for example, dissolved 

oxygen concentrations at a centrally located site in the 

northern part of the lake were less than 3 mg/L continuously 

for 8 days, and the spatial extent of the LDOE was nearly 

40 km2

 (Wood and others, 2006). This event culminated in 

a smaller fish die-off than those of the mid-1990s (Adams 

and others, 2003). A die-off of the severity of the mid-1990s 

has not occurred since the water-quality monitors were 

installed in 2002; therefore, it is not an annual occurrence. 

The impact of large fish die-offs on the endangered sucker 

populations is, nonetheless, devastating; thus, it is important to 

understand what combination of conditions results in the most 

severe events before a management strategy can be devised 

to modulate the algal bloom-decline cycle. Understanding 

the causes of the precipitous bloom decline is critical in 

this regard because the loss of photosynthetic production in 

addition to the added oxygen demand generated by senescing 

cells is the root cause of LDOEs. It also has been observed that 

the events are more severe in the northern sucker habitat area 

of the lake than in areas of similar depth in the rest of the lake, 

and observations of currents have indicated that circulation 

patterns play an important role. 

The surface area to volume ratio of a shallow lake 

contributes to greater primary production per volume and 

enhanced nutrient cycling between the water column and 

sediments in comparison to a deep lake (Scheffer, 1998). 

Thus, the hypereutrophication of UKL has similarities to 

processes observed in other shallow lakes that experience 

massive cyanobacterial blooms. At the same time, each lake is 

unique in some way, and UKL is no exception. A particularly 

important feature of UKL that has implications for both the 

hydrodynamics and water quality of the lake is its bathymetry. 

Although most of the lake is shallow, a relatively deep 

trench runs along the western shoreline (fig. 1). This trench 

is apparent in the bathymetry as far south as Buck Island. 
4 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

0 ,

where

 and are components of velocity in the 

horizontal directions and ,

respectively, and

 is the component of velocity in the

vertical direction, .

u v u

x y z

u v

x y

w

z

∂ ∂ ∂

+ + =

∂ ∂ ∂

(1)

After integrating over depth and applying the kinematic 

condition at the free surface, the mass conservation equation 

becomes

0

( )

∂η ∂ ∂ + + =

∂ ∂ ∂

η

 ,

where

is the free surface, and

 and are the horizontal components of,

velocity integrated over the depth of

the water column from the bottom 

 = - to the free surface at = 

U V

t x y

x,y,t

U V

z h z η.

(2)

 and . h h U udz V vdz η η

− − = = ∫ ∫ (3)

The momentum conservation equations in three spatial 

dimensions are

0

2 2

2 2

1

 ,

a

h z

u u u u p u v w fv g

t x y z x x

q u u u g K K

x x x y z z

∂ ∂ ∂ ∂ ∂ ∂η + + + − = − −

∂ ∂ ∂ ∂ ρ ∂ ∂

∂ρ ∂ ∂ ∂ ∂ ∂ ′     − − + + +     ∂ ∂   ∂ ∂ ∂ ∂   (4)

0

2 2

2 2

1

 ,

a

h z

v v v v p u v w fu g

t x y z y y

q v v v g K K

y y x y z z

∂ ∂ ∂ ∂ ∂ ∂η + + + + = − −

∂ ∂ ∂ ∂ ρ ∂ ∂

∂ρ ∂ ∂ ∂ ∂ ∂ ′     − − + + +     ∂ ∂   ∂ ∂ ∂ ∂   (5)

and

2 2

2 2

 .

h

z

w w w w q w w u v w K

t x y z z x y

w K

z z

∂ ∂ ∂ ∂ ∂ ∂ ∂  

+ + + = − + +   ∂ ∂ ∂ ∂ ∂   ∂ ∂

∂ ∂   +   ∂ ∂  

(6)

It runs across the mouth of Howard Bay, along the western 

shoreline and to the west of Bare Island, and then turns west 

around Eagle Point, runs across the entrance to Shoalwater 

Bay and along Ball Point, before turning north at the entrance 

to Ball Bay and fading away. The currents in the trench are 

the strongest in the lake and are aligned with the bathymetry; 

under prevailing wind conditions, water flows northward 

through the trench (Gartner and others, 2007).

Recognizing that a detailed understanding of the 

movement of water around the lake was essential to a 

complete understanding of the water quality of the lake, 

the Bureau of Reclamation and the U.S. Geological Survey 

(USGS) entered into a cooperative agreement in 2005 to 

develop a hydrodynamic model of the lake. This report is the 

first comprehensive documentation of the UKL model and 

its application, although singular aspects of the model have 

appeared in proceedings papers (Cheng and others, 2005; 

Wood and Cheng, 2006). 

Purpose and Scope

The purpose of this report is to provide a comprehensive 

accounting of the boundary, forcing, and calibration 

data available during the summer field season in 2005 

(approximately mid-June through mid-September) and 2006 

(approximately mid-May through mid-October), to provide 

a detailed description of the source and boundary terms 

that have been added to the model in order to accurately 

simulate heat transport, and to document the calibration and 

validation of the numerical model. The report also discusses 

the implications of lake circulation for water quality, with 

particular emphasis on how circulation might affect dissolved 

oxygen concentrations in the northern part of the lake. The 

most important numerical features of the model are briefly 

summarized below, but the reader is referred to the references 

cited in this report for the details of the numerical methods 

used in the computational core of the UnTRIM model.

Model Description

Governing Equations

The governing equations for three-dimensional, 

baroclinic circulation and the transport of scalar variables 

in a lake include the conservation equations of mass and 

momentum, a kinematic free-surface equation (derived from 

mass conservation), an equation of state relating density to 

temperature, and a conservation equation for each scalar 

variable. The lake is assumed to be sufficiently large for a 

Coriolis acceleration term to be included in the momentum 

equations. To simplify the governing equations, the water 

is assumed to be incompressible and it is assumed that the 

Boussinesq approximation applies. In Cartesian coordinates, 

the mass conservation equation is
Model Description 5

The Coriolis parameter f is assumed constant over the domain. 

Kz and Kh and are the vertical and horizontal viscosity 

coefficients, respectively, ρ0 is a constant reference density of 

water (1,000 kg/m3), g is gravitational acceleration (9.81 m s-2), 

and ρ' is the vertically integrated and normalized deviation 

from reference density. The above equations have been derived 

by decomposing the pressure P into a hydrostatic component, 

dependent on the free-surface elevation η and the spatially 

varying density ρ, and a nonhydrostatic component q as:

( )

0 ( [ ] )

( , , , ) , ,

( , , ) ( , , ) ( , , , ) ,

a p x y z t p x y t

g x y t z g x y t q x y z t

=

+ ρ η − + ρ + ′ (7)

where 

0

0

.

 ,

where

is the variable of integration, and 

is atmospheric pressure

η ρ − ρ ρ = ζ ′ ρ

ζ

∫z

a

d

p

(8)

The equation for conservative transport of a scalar variable 

(including heat and solute) with concentration C is:

( ) ( ) ( )

 ,

where

is a settling velocity for the scalar, and 

and are the vertical and horizontal

diffusivity coefficients, r

∂ ∂ ∂ ∂ −     ∂ ∂   + + + =   ∂ ∂ ∂ ∂ ∂ ∂  

∂ ∂ ∂ ∂     +   +   ∂ ∂ ∂ ∂    

s

h

h z

s

z h

C uC vC w w C C D

t x y z x x

C C D D

y y z z

w

D D

espectfully.

(9)

The final equation is the equation of state, which in the 

general case relates density to temperature and salinity. In this 

freshwater application the equation of state is an empirical 

relation that relates density in kg m-3 only to water temperature 

Tw in °C (Gill, 1982):

2

3 2 4 3 6 4

9 5

999.842594 6.793952 10 9.095290

10 1.001685 10 1.120083 10

6.536332 10 .

w

w w w

w

T

T T T

T

−

− − −

−

ρ = + × −

× + × − ×

+ × (10)

Equations (2) through (10) are the governing equations 

solved within the computational core of the UnTRIM 

model. If the water density can be approximated as constant 

(barotropic flows), the scalar transport equation (9) is 

uncoupled from the momentum equations (4), (5), and (6). 

When the hydrostatic approximation is made, the horizontal 

gradients of q drop out of equations (4) and (5), and equation 

(6) is replaced by the hydrostatic approximation. 

The governing system of equations can be solved 

efficiently by a semi-implicit finite-difference method that 

is computationally fast, accurate, and stable over a regular 

computational mesh as discussed by Casulli and Cheng 

(1992) and Casulli and Cattani (1994). A weak CourantFriedrich-Lewy (CFL) stability condition is imposed on the 

computational time step due to the explicit treatment of the 

horizontal diffusion in the momentum equations. If density 

stratification in the lake is considered, the transport equations 

are coupled with the momentum equations through the density 

gradient terms. In this case, the density gradients in the 

momentum equations and the mass conservation of the scalar 

variables (including heat and solute) are solved explicitly. 

The resulting numerical scheme is subject to an additional 

weak CFL stability condition on the computational time step 

due to the explicit treatment of the transport equation and the 

baroclinic pressure terms in the momentum equations.

Orthogonal Unstructured Grids 

The governing equations are solved in physical space 

without invoking any coordinate transformation in the 

horizontal or vertical directions. The stability properties of the 

governing partial differential equations are controlled by using 

a semi-implicit finite-difference scheme (Casulli, 1990; Casulli 

and Cheng, 1992); the resulting numerical algorithm is robust 

and computationally efficient. Traditional finite-difference 

schemes resort to refining the rectangular finite-difference 

mesh when a complicated domain is encountered in order to 

resolve the flow distributions in narrow and confined regions. 

The resulting fine-resolution grids in broad and open regions 

are unnecessary, and the computational mesh consumes a large 

portion of computing resources, which is not computationally 

efficient. In the UnTRIM model, the semi-implicit finitedifference method is applied over an unstructured grid (Casulli 

and Zanolli, 1998; Casulli and Walters 2000) in which fine 

grid resolutions are used in complex regions, and relatively 

coarse grids are used in broad and open areas. 

When a computational mesh is created, the horizontal 

domain is covered by a set of nonoverlapping convex 

polygons. Each side of a polygon is either a boundary line 

or a side of an adjacent polygon. The center of each polygon 

is defined such that the segment joining the centers of two 

adjacent polygons intersects the side shared by the two 

polygons and is orthogonal to it (fig. 2). The center of a 
6 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

polygon does not necessarily coincide with its geometrical 

center, except in special cases such as rectangular finitedifference grids and grids of uniform equilateral triangles. 

This type of computational mesh is called an unstructured 

orthogonal grid (Casulli and Zanolli, 1998; Casulli and 

Walters, 2000). In an unstructured grid representation within 

the UnTRIM model, each polygon has either three or four 

sides. Each polygon, side, and vertex in the grid is assigned 

a unique number. The x–y coordinates of the vertices of each 

polygon must be specified, as well as the information that 

relates each polygon to its vertices and sides. 

In the vertical direction, the water body to be modeled 

is divided into layers by horizontal planes. The polygons 

between the horizontal planes become a stack of prisms whose 

thickness is the prescribed layer thickness. The water-surface 

elevation is assumed to be uniform within each polygon and is 

defined at the center of the polygon. The velocity component 

normal to each face of a prism is assumed to be uniform over 

the face. The velocity is defined at each vertex in the middle 

of each layer, and the spatial distribution of velocity through 

the prism is obtained by interpolating between the velocities 

at the vertices. Finally, the water depth below each polygon is 

specified and assumed uniform on each side.

The unstructured orthogonal grid for UKL (fig. 3) was 

created and optimized with the JANET software (Lippert 

and Sellerhoff, 2006). This software enables semiautomated 

creation of an initial grid, followed by manual fine-tuning 

of individual polygons that are identified by tests within 

the software to be outside of a prescribed tolerance for 

orthogonality. The resulting grid for UKL has 8,389 polygons, 

4,555 vertices, and 28 layers of 0.5 m thickness. Because most 

of the lake is shallow, however, the number of active prisms 

in each layer varies from more than 8,000 in each of the top 

3 layers to less than 100 in each of the 11 deepest layers of 

the grid. The length of the sides of the polygons varies from 

approximately 100 m in areas along the trench where the 

bathymetry changes rapidly to approximately 400 m in the 

open areas of the lake. 

Numerical Approximation 

A semi-implicit finite-difference scheme is used to 

obtain an efficient numerical algorithm for which stability is 

independent of free-surface gravity waves, wind stress, vertical 

viscosity and bottom friction, for both the hydrostatic (Casulli 

and Walters, 2000; Casulli, 1990) and nonhydrostatic case 

(Casulli and Zanolli, 2002; Casulli, 1999a, Casulli, 1999b). 

The momentum equations (4)–(6) are finite-differenced in the 

direction normal to each vertical face of each computational 

prism (along oa, ob, and oc in fig. 2). The momentum 

equations relate the gradient of water-surface elevation 

between adjoining prisms to the face velocity on the common 

face between these prisms. The vertical mixing and the 

bottom friction are discretized implicitly in time for numerical 

stability, whereas an explicit finite-difference operator is used 

to solve the wind stress and the advection and horizontal 

dispersion terms. This operator can take several particular 

forms; an Eulerian–Lagrangian scheme is used in UnTRIM 

(Casulli and Cheng, 1992). For stability, the implicitness 

factor θ must be in the range 1/2 ≤ θ ≤1 (Casulli and Cattani, 

1994). In the vertical direction, a simple finite-difference 

discretization that does not require uniform layers is adopted. 

The vertical thickness of the top and bottom layers can vary 

spatially and the thickness of the top layer also can vary with 

time. The vertical thickness of the top layer is allowed to 

become zero, in order to accommodate the drying of cells.

The free-surface equation (2) is discretized semiimplicitly (Casulli and Cattani, 1994; Casulli and Walters, 

2000), and only the face velocities are needed to complete the 

finite volume balance of total volume below each polygon. 

The surface elevation at the center of each polygon is 

determined by substituting the finite-differenced momentum 

equations on all faces of the polygon into the continuity 

equation. The resulting matrix equation governs the watersurface elevation over the entire domain. This matrix equation 

is strongly diagonally dominant, symmetric, and positive 

definite; thus its unique solution can be efficiently determined 

by preconditioned conjugate gradient iterations until the 

residual norm becomes smaller than a given tolerance (Golub 

and van Loan, 1996). Once the free surface of the entire 

domain for the next time level has been calculated, the normal 

velocities on the faces of the prisms are calculated by back 

Figure 2. An example of an unstructured orthogonal grid.

OR19_0141_UKL_fig02

c

o a

b
Model Description 7

substitution. The details of the 

finite-difference equations are not 

reproduced here and readers are 

referred to Casulli and Walters 

(2000).

The transport equation (9) 

for scalar variables is solved 

explicitly by using the velocity 

field obtained from the previous 

time step. An explicit finitevolume discretization that both 

conserves mass and guarantees 

that the solution will be bounded 

by the maximum and minimum of 

the initial and boundary values is 

used. In addition, a “flux limiter” 

is used in the discretization of 

the equation in order to limit 

numerical diffusion and preserve 

high accuracy at grid resolutions 

that are well-matched to the 

domain (Casulli and Zanolli, 

2005). The flux limiter used in the 

simulations discussed here is the 

Superbee function (Roe, 1986).

In summary, the numerical 

solution scheme used in the 

computational core of UnTRIM to 

solve the governing equations is 

designed to achieve computational 

efficiency while retaining the 

accuracy of the solution. For 

the UKL simulations discussed 

in this report, the hydrostatic 

approximation was made. The 

simulations included the transport 

of three scalar quantities: 

turbulent kinetic energy, the rate 

of dissipation of turbulent kinetic 

energy, and heat. The equations 

describing the transport of each 

of these scalars are discussed 

below. Using the numerical 

solution scheme described above, 

the UnTRIM model required 

an average of 1.3 seconds of 

computational time per 2-minute 

simulation time step to solve the 

governing equations on the UKL 

grid, using a Dell 620 Latitude 

1.66 GHz notebook computer with 

an Intel T2300 dual processor. A 

list of symbols used in this report 

is provided in table 1.

Figure 3. Unstructured orthogonal grid used in model simulations of Upper Klamath 

Lake, Oregon.

OR19_0141_UKL_fig03

0 1 2 3 4 5 MILES

0 1 2 3 4 5 KILOMETERS
8 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

Symbol Name Units

a Constant in equation for surface drag coefficient (m s-1)-0.5

B Production of turbulent kinetic energy by buoyancy m2

 s-3

C Concentration of scalar variable arbitrary

C10 10 m wind surface drag coefficient –

CB Empirical coefficient mb °C -1

cp Specific heat capacity of water at 20°C J g-1 oC-1

c0 Constant used to relate turbulent kinetic energy to rate 

of turbulent kinetic energy dissipation

–

Dz Vertical turbulent diffusivity coefficient m2

 s-1

Dh Horizontal turbulent diffusivity coefficient m2

 s-1

es Saturation vapor pressure of air at water surface 

temperature

mb

ea Vapor pressure of air mb

f Coriolis parameter s-1

g Gravitational acceleration m s-2

HB Back radiation from water surface W m-2

HE Evaporative cooling W m-2

HLW Incoming longwave radiation W m-2

HS Sensible heat loss W m-2

HSW Incoming shortwave radiation W m-2

h Depth of the water column m

Kz Vertical turbulent viscosity coefficient m2

 s-1

Kh Horizontal turbulent viscosity coefficient m2

 s-1

k Turbulent kinetic energy m2

 s-2

ke Extinction coefficient m-1

Lw Latent heat of evaporation J kg-1

l Length scale of turbulent eddies –

M Fraction of the sky covered by clouds –

P Production of turbulent kinetic energy by mean shear m2

 s-3

p Pressure kg m-1 s-2

pa Atmospheric pressure kg m-1 s-2

psl Atmospheric pressure at sea level kg m-1 s-2

q Non-hydrostatic component of pressure m2

 s-2

RSW Shortwave radiation W m-2

rB Bottom friction coefficient –

rh Relative humidity –

rLW Fraction of incident longwave radiation reflected at the 

water surface

–

rSW Fraction of incident shortwave radiation reflected at the 

water surface

–

SD Stability function used in calculation of turbulent 

diffusivity for constituent transport

–

Symbol Name Units

SK Stability function used in calculation of turbulent 

diffusivity for momentum transport

–

Ta Air temperature °C

Tw Water temperature °C

t Time s

U,V Depth-integrated components of velocity in the 

horizontal plane

m s-1

u,v Components of velocity in the horizontal plane m s-1

r

uB

Near-bottom velocity m s-1

r

ui

Wind vector at i m above the surface m s-1

r

uT

Current vector at the surface of the water column m s-1

r

u*

Friction velocity m s-1

w Veritical component of velocity m s-1

ws Settling velocity m s-1

x,y Cartesian coordinates in the horizontal plane m

z Vertical cartesian coordinate m

zB Location of bottom computational grid point m

z0 Bottom roughness height m

α0 Proportionality constant °K-2

η Free surface elevation m

ε Rate of dissipation of turbulent kinetic energy m2

 s-3

a ε Emmisivity of air –

w ε Emmisivity of water –

κ von Karman constant –

ρ Density of water kg m-3

ρ ́ Depth-integrated deviation of water density from 

reference density, normalized by reference density

m

ρ0 Reference density of water; density at 20°C kg m-3

ρa Density of air kg m-3

σ Stefan-Boltzmann constant W m-2 °K-4

σk Turbulent Schmidt number for turbulent kinetic energy –

σε Turbulent Schmidt number for dissipation of turbulent 

kinetic energy

–

B τ

r Stress vector at the sediment/water interface N m-2

T τ

r Stress vector at the water surface N m-2

ζ Variable of integration –

Table 1. List of symbols used in this report.

[Abbreviations: m, meter; s, second; J, joules; g, gram; °C, degrees Celsius; mb, millibar; W, Watt; kg, kilogram; N, Newton; °K, degrees Kelvin; –, no units]
Model Description 9

1/2 (2 ) . D k lS z = D (14)

The length scale of the turbulent eddies, l, is related to ε by 0 3 3/2 1 l c k ( ) − = ε μ , where 0 cμ is a constant that depends on

the stability functions used. The stability functions SK and 

SD are algebraic functions of k and ε that are derived by 

converting the system of differential equations for the secondorder moments of turbulent quantities to a system of algebraic 

equations, thus eliminating a great deal of computational 

effort. The stability functions used in this application were 

those derived by Kantha and Clayson (1994). The solution 

of the equations (11) and (12) for k and ε is accomplished 

sequentially at each time step. The advective and diffusive 

transport is solved as described above for a scalar quantity. 

Then P, B, and ε are calculated from the velocity and density 

gradients, and then the source and sink terms of equations (11) 

and (12) are used to update k and ε. The implementation of 

this turbulence closure in UnTRIM is described by Celebioglu 

and Piasecki (2006a and 2006b).

Surface and Bottom Boundary Conditions

The shear stress at the surface of the water column, T τ

r , is 

specified in terms of a surface drag coefficient C10:

10 10 10

10

( ) 

where

is the density of the air,

is the wind vector at 10 m above the water

surface, and

is the current vector at the surface of the

water column.

, T a T T

a

T

C u u u u

u

u

τ = ρ − −

ρ

r r r r r

r

r

(15)

The surface drag coefficient is in general an increasing 

function of the wind speed, and various formulations have 

been proposed and evaluated (Garratt, 1977; Smith and Banke, 

1975). At very low wind speeds (less than about 3 m s-1) 

experimental evidence shows that the wind surface drag 

coefficient increases with decreasing wind speed (Wuest and 

Lorke, 2003). In this application, an increasing C10 at low wind 

speeds did not bring the simulated currents closer to observed 

currents, and the following function for C10 was used: 

C a u 10 = 10

r for 10-m wind speeds less than 15 m s-1

and 10 C a = 15 for greater wind speeds. The value of a was 

treated as a calibration parameter, although the final choice of 

a = 0.0005 was the value proposed by Wu (1969).

The shear stress at the sediment/water interface B τ

r is 

defined in terms of a bottom friction coefficient rB and the 

near-bottom current speed B u

r as: 

Turbulence Closure

The vertical turbulent dispersion of momentum and 

diffusion of heat are calculated within the model as spatial 

and temporal functions of two localized flow characteristics—

turbulent kinetic energy k and the rate of dissipation of 

turbulent kinetic energy ε. Equations for the transport of these 

two turbulence characteristics can be written as:

 

( ) ( ) ( )

 ,

h h

z

k

k uk vk wk k k D D

t x y z x x y y

K k P B

z z

∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂    

+ + + =   +   ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂    

∂ ∂  

+   + + − ε

∂ σ ∂  

(11)

( )

( 1 3 2 )

( ) ( )

 

where

 and are the production of turbulent kinetic

energy by velocity shear and

 buoyancy, respectively, an

, ε

∂ε ∂ ε ∂ ε ∂ ε ∂ ∂ε ∂ ∂ε    

+ + + = +     ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂    

∂ ∂ε ε  

+   + + − ε

∂ σ ∂  

h h

z

u v w

D D

t x y z x x y y

K

c P c B c

z z k

P B

, .

d 

 and are the turbulance Schmidt numbers for

 and respectively

σ σ

ε ε k k

(12)

The implementation of this two-equation turbulence 

closure model within UnTRIM is accomplished by using the 

generic length-scale approach of Umlauf and Burchard (2003) 

and Warner and others (2005), which provides a convenient 

and flexible means for implementing a wide variety of twoequation turbulence closure models. The choice of the generic 

coefficients—c1, c2, and c3—that apply to the source and 

dissipation terms for the length-scale parameter ε determines 

which turbulence closure model is implemented. The k – ε

model (Rodi, 1993) was used in the simulations presented 

in this report. In accordance with Warner and others (2005), 

the choice of c1

 = 1.44, c2 = 1.92, and c3 = –0.52 was used to 

implement the k – ε model in the case of stable stratification 

(B < 0). In the case of unstable stratification (B > 0), the weak 

thermal stratification in UKL generally resulted in too much 

mixing when buoyancy was considered to be a source of k

and ε, so the buoyancy term B was set to zero in equations 

(11) and (12) when unstable stratification occurred. Equations 

(11) and (12) are solved within the model and the turbulent 

diffusivities, Kz

 and Dz for momentum and constituent 

transport, respectively, are calculated at each time step as:

1/2 (2 ) K k lS z = K (13)

and
10 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

temperature profiles in the lower atmosphere are available, 

the thickness of the surface flow is adjusted to the height at 

which the kinetic energy in the wind is equal to the work done 

against buoyant restoring forces. Surface winds flow over 

features of the terrain that are lower than this thickness, and 

around features that are higher than this thickness. Because 

temperature soundings in the near-surface around UKL are 

not part of the routine data collection around the lake, the full 

functionality of this model was not used. Instead, a stable 

and invariant adiabatic temperature profile in the atmospheric 

boundary layer was assumed.

Heat Transport

Source/sink terms that describe the transfer of heat 

across the air/water interface were added to the equation for 

conservative transport of a scalar variable (eq. 9). Surface 

terms include incoming shortwave and longwave radiation, 

outgoing longwave radiation, evaporative cooling, and 

conduction/convection. Computationally, all of these terms 

except for incoming shortwave radiation are included in the 

surface boundary condition. Incoming shortwave radiation is 

treated as an internal source/sink term that allows the radiation 

to be absorbed through a finite distance in the upper layers 

of the model water column rather than only at the air/water 

interface. The three-dimensional equation describing the 

transport of heat is:

 ( )

0

( ) ( )

1 ,

w w w w w

h

w w SW

h z

p

T uT vT wT T

D

t x y z x x

T T R

D D

y y z z c z

∂ ∂ ∂ ∂ ∂   ∂

+ + + =   ∂ ∂ ∂ ∂ ∂ ∂  

∂   ∂ ∂   ∂ ∂

+   + +   ∂ ∂ ∂ ∂ ρ ∂    

(19)

where Tw is the water temperature, and

0

1 SW

p

R

c z

∂

ρ ∂

(20)

is the rate of temperature change at any depth z due to local 

absorption of shortwave radiation. The shortwave radiation at 

any depth z in the water column is given by:

( ) (1 )

2

 

where

is the shortwave radiation incident at the water

surface in units of Watts/m ,

is the fraction of that is reflected at the

 water surface (the albedo),

is the ext

, − = − e k z

SW SW SW

SW

SW SW

e

R z H r e

H

r H

k

o C .

inction coefficient, and

is the specific heat of water at 20 p c

(21)

0

 * *

where

is the friction velocity. *

, B

B B B

r u u u u

u

ρ

τ

= =

r

r r r r

r

(16)

Within the bottom logarithmic profile layer, the speed of the 

currents at distance z' away from the interface is given by 

(Schlichting, 1955):

*

0

0

' ( ') ln 

where

is the von Karmen constant (0.41) and

is the roughness height of the bed. 

, u z

u z

z

z

  =  

κ  

κ

r r (17)

Combining these two expressions allows the solution of rB in 

terms of the distance zB of the lowest computational grid point 

away from the bottom, where B u

r is defined:

2

2

0

ln . B

B

z r z

−

   

= κ          

(18)

Wind Interpolation

The effort to develop an UnTRIM model of UKL began 

in 2003, and these efforts were reasonably successful at 

reproducing the hydrodynamics of the lake (Cheng and others, 

2005). The early simulations showed that the currents in the 

lake were highly responsive to wind forcing at the surface, 

and it was determined that the biggest improvement in model 

performance would be obtained by generating a spatially 

accurate surface wind, rather than using a uniform surface 

wind over the entire lake. By 2005, additional monitoring 

sites had been established around the lake, as described 

below. Observations from these monitoring sites showed that 

the surrounding land topography influences the wind field 

over the lake, as evidenced, for example, by the fact that the 

prevailing wind direction over the northern part of the lake is 

westerly, but is northwesterly over the southern part of the lake 

(Hoilman and others, 2008). 

An objective analysis model of the atmospheric boundary 

layer was used to interpolate the observations from the various 

wind-monitoring locations to a uniform grid (500-m spacing, 

both east-west and north-south) over the lake (Ludwig and 

others, 1991). The procedure for doing this involves starting 

with an initial “guess” interpolation and then modifying 

the initial guess iteratively to satisfy mass continuity. When 
Model Description 11

The internal source/sink term is applied to the water 

temperature outside of the computational core of the model, 

with the result that at each time step, temperature over the 

entire domain is first updated to an intermediate value by 

physical transport; then, the change in heat resulting from 

the absorption of shortwave radiation is applied to the 

intermediate temperature. 

The rest of the surface heat transfer terms are treated 

computationally as a boundary condition because they apply 

directly at the water surface. The boundary condition at the 

water surface is:

0

( )

 , w LW B E S

z

p

T H H H H

D

z c

∂ − − − = − ∂ ρ

(22)

where HLW , HB, HE , and HS are the incoming longwave 

radiation, reflected longwave radiation, evaporative heat flux, 

and sensible heat flux, respectively, all with units of Watts/m2

. 

This boundary condition is updated at every time step and the 

new value is supplied to the computational core of UnTRIM.

Of the surface heat transfer terms, only the incoming 

shortwave radiation HSW is directly measured. The rest of the 

surface terms are calculated within the model at each time step 

by using standard formulations (McCutcheon, 1989; Martin 

and McCutcheon, 1999). Atmospheric long wave radiation is 

calculated from the Stefan-Boltzmann law, modified for the 

emissivity of the air:

4 (1 )( 273.16)

,

LW a LW a

a

LW

a

H r T

r

T

= ε σ − +

ε

σ

 

where

is the emissivity of the air,

is the Stefan-Boltzmann constant

is the fraction of longwave radiation that is

reflected at the water surface, and

is air temperatu

,

re in o C.

(23)

The emissivity of the air is calculated as:

2 2

0

0

(1 0.17 )( 273.16) a M Ta

M

ε = α + +

α

 

where

is a proportionality constant, and

is the fraction of the sky covered by clouds.

, (24)

Back radiation from the water surface also is calculated 

from the Stefan-Boltzmann law:

4 ( 273.16) 

where

is the emissivity of the water, taken to be a

constant.

= ε σ + ,

ε

B w w

w

H T (25)

Evaporative heat loss is calculated as:

0 ( ) , H L F e e E w s a = ρ − (26)

where

2

2

0.26(0.5 0.54 ) 

where

as the function of wind speed measured at height

 2 m ( ) is a modified version of the original

 Penman equation (Janssen, 2006),

is the latent heat of evaporation and is

= + , r

r

w

F u

F

u

L taken to be

constant, and

is the vapor pressure in the air, both in millibars. s e

(27)

The saturation vapor pressure at the water-surface temperature 

is given by (Dingman, 2002): 

17.3 6.11exp 237.3 , w

s

w

T

e

T

  =     +

(28)

and the vapor pressure in the air is calculated from relative 

humidity rh as:

17.3 6.11exp 237.3 100 . a

a

a

T rh

e

T

   =      +   (29)

Sensible heat loss is often calculated from evaporative 

heat loss using the Bowen ratio:

 

where

is an empirical coefficient,

is local atmospheric pressure, and

is atmospheric pressure at sea level.

,   − =     −

a w a

S E B

sl s a

B

a

sl

p T T H H C

p e e

C

p

p

(30)

The sensible heat loss term is small compared to the other 

terms in eq. (22). In the UKL model the inclusion of this term 

tended to reduce the accuracy of the simulation, particularly 

during periods of atmospheric cooling when HS was negative, 

so sensible heat loss was not included in the model simulations 

presented in this report. 
12 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

Datasets

Tributary Inflows and Outflows

Inflows

The inflows to the UKL/Agency Lake hydrologic system 

include the Williamson River, the Wood River, the Fourmile 

and Sevenmile Canals, precipitation, and ground water. The 

biggest tributary is the Williamson River, which flows into 

UKL north of Eagle Ridge on the opposite shore (fig. 1). The 

Williamson River streamflow gage is located at river mile 

10.3; daily values from the gage were used as input to the 

model (USGS site 11502500, fig. 1; Herrett and others, 2006; 

U.S. Geological Survey, 2007). Daily values of inflow from 

the next three largest tributaries—the Wood River, Sevenmile 

Canal, and Fourmile Canal—were provided by Graham 

Matthews and Associates (GMA) for 2005 (Cort Pryor, 

Graham Matthews and Associates, oral commun., 2006). 

These three tributaries empty into Agency Lake from the 

Wood River Valley (fig. 1). In 2005, the discharge at Fourmile 

Canal was estimated from August 13 through the end of the 

water year because equipment at the site failed. In 2006, daily 

values of inflow at Sevenmile Canal were provided by GMA 

(Cort Pryor, Graham Matthews and Associates, oral commun., 

2006), but daily values were unavailable at either the Wood 

River or the Fourmile Canal. The discharge from these two 

sites was set to 3 m3/s and 5 m3/s, respectively, which was the 

approximate base flow based on the 2005 data. The remainder 

of the discharge from those two tributaries, which probably 

was substantial only during the early part of the summer, 

became part of the calculation of ungaged inputs during 2006 

(fig. 4). 

Daily precipitation data were collected at the Bureau 

of Reclamation (2008) Pacific Northwest Cooperative 

Agricultural Weather Network (AgriMet) site located near 

Agency Lake (AGKO, fig. 1). Because the UnTRIM model 

of UKL does not accommodate precipitation as a distributed 

input, one-third of the daily precipitation was added to 

the Wood River inflow, and two-thirds was added to the 

Williamson River inflow. 

The ungaged inflows to the lake, including ground-water 

inflows, were calculated by using the daily lake stage (USGS 

site 11507001), which is a weighted average of the stage 

measured at three sites around the lake at Rattlesnake Point, 

Rocky Point, and Klamath Falls (fig. 1), and the volume of 

the lake as a function of lake stage, which was calculated 

from the same bathymetry used to create the model grid. 

DISCHARGE OR PRECIPITATION, IN CUBIC METERS PER SECOND 

15 1 1 15 1 15 1 15 30

MAY JUNE JULY AUGUST SEPTEMBER

1

2005

0

20

40

60

80

100

120

Williamson River near Chiloquin

Wood River at Dike Road

Sevenmile Canal at Dike Road

Fourmile Canal at Lower Weir

Precipitation

Ungaged inflow

A.

Figure 4. Inflows to the Upper Klamath Lake/Agency Lake hydrologic system in Oregon in (A) 2005 and (B) 2006.
Datasets 13

This bathymetry is based on a survey of approximately 

94,200 soundings done in May and June of 1996. These 

soundings were oriented in east/west transects, providing a 

200-m north/south spacing between the transects, and were 

thinned to a minimum of 50-m spacing along the transects 

(M. Neuman, Bureau of Reclamation, oral commun., 2007). If 

the change in volume of the lake at each time step is denoted 

ΔV, then the ungaged inflow, R, at each time step is calculated 

as:

R outflows inflows V = ∑ ∑− + ∆ . (31)

The ungaged inflow is sometimes large compared to 

inflow from tributaries to the lake, particularly in the spring 

and early summer, and at other times is negative, indicating 

ungaged outflow rather than inflow (fig. 5). It is also, at times, 

large when compared to the estimate of 8.5 m3/s of groundwater inflow to the lake (Hubbard, 1970). Within the model, 

when the ungaged inflow is positive it is added to the Wood 

River Valley inflows. When the ungaged inflow is negative 

(indicating ungaged outflow) it is added to the Link River 

discharge at the outlet of the lake. It is apparent that the gaged 

inflows to and outflows from the lake are not sufficient to 

close the water balance with high accuracy. The method used 

for calculating the ungaged inflow and applying it to the 

model assures, however, that a mass balance of water in the 

model is achieved. 

The outflow from the lake at the Link River is gaged by 

the USGS (site ID 11507500; fig. 1). Daily values were used 

as input to the model. This gaged discharge does not include 

the portion of the flow that, prior to 2005, was diverted for 

power generation through a canal that ultimately discharges 

back into the Link River downstream of the gage. In 2005, no 

water was diverted through the power canal, so the discharge 

at the gage is an accurate measure of the total flow over the 

Link River Dam. In 2006, there may have been diversions 

through the power canal, but the discharge measurements 

through the canal were unavailable and, therefore, were 

incorporated into the model as an ungaged outflow (fig. 5). 

The A canal (fig. 1) diverts water for irrigation from 

the southern end of the lake above the Link River dam. The 

flow through the A canal is gaged by the Klamath Irrigation 

District and the daily values were obtained from the Bureau 

of Reclamation (Paul Gregory, Bureau of Reclamation, oral 

commun., 2007).

Daily evaporation data were collected at the AgriMet 

site located near Agency Lake (AGKO, fig. 1) (Bureau of 

Reclamation, 2008). Because the UnTRIM model of UKL 

does not accommodate evaporation as a distributed outflow, 

the daily evaporation was added to the Link River outflow.

OR19_0141_UKL_fig_4b DISCHARGE OR PRECIPITATION, IN CUBIC METERS PER SECOND 

15 1 15 1 15 1 15 1 15 30

MAY JUNE JULY AUGUST SEPTEMBER

1

2006

0

20

40

60

80

100

120

Williamson River near Chiloquin

Wood River at Dike Road

Sevenmile Canal at Dike Road

Fourmile Canal at Lower Weir

Precipitation

Ungaged inflow

B.

Figure 4.—Continued.
14 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

OR19_0141_UKL_fig_5a DISCHARGE, IN CUBIC METERS PER SECOND 

15 1 15 1 15 1 15 1 15 30

MAY JUNE JULY AUGUST SEPTEMBER

1

2005

0

20

40

60

80

100

120

Link River

Evaporation

A Canal

Ungaged outflow

A.

DISCHARGE , IN CUBIC METERS PER SECOND 

MAY JUNE JULY AUGUST SEPTEMBER

2006

0

20

40

60

80

100

120

Link River

Evaporation

A Canal

Ungaged outflow

1 15 1 15 1 15 1 15 1 15 30

B.

Figure 5. Outflows from the Upper Klamath Lake/Agency Lake hydrologic system in Oregon in (A) 2005 and (B) 2006.
Datasets 15

Meteorological Data

Wind data were collected during 2005 and 2006 at two 

sites located on the lake and four sites located around the 

shoreline (fig. 1). After July 26, 2005, reliable data were 

available from at least five of the six sites (Hoilman and 

others, 2008). Prior to July 26, however, wind data were 

available from fewer sites. For this reason, three types of 

wind forcing covering different periods were used to run 

model simulations in 2005 (table 2). From June 8 (when the 

raft-mounted wind sensor at MDL was deployed in the lake 

at the beginning of the field season) to August 31 (when the 

wind sensors at MDL failed), the model was forced with a 

spatially uniform wind as measured at site MDL (denoted the 

UW-MDL simulation). From July 26 to October 12 (when 

the rafts on the lake were removed at the end of the field 

season), the model was forced with a spatially variable wind, 

interpolated as described above between sites on and around 

the lake (denoted the VW simulation). This resulted in an 

overlap period in midsummer between July 26 and August 

31 for which these two simulations could be compared. To 

provide an additional comparison, during the same overlap 

period the model was forced with a spatially uniform wind as 

measured at site MDN (denoted the UW-MDN simulation). 

During 2006, wind data were available from six sites and a 

spatially variable wind forcing was used between May 20 

and October 17. All meteorological data were collected as 

15-minute averages in 2005 and 10-minute averages in 2006.

Air temperature and relative humidity were measured 

at five sites—WMR-MET, BLB-MET, HDB-MET, SSHRMET, and MDL. Shortwave solar radiation was measured at 

sites WMR-MET and SSHR-MET. The measurements of air 

temperature and relative humidity at site MDL were used in 

these simulations with the exception of the period following 

September 4, 2005, when the sensors at site MDL failed. After 

September 4, 2005, the air-temperature data measured at the 

land-based station WMR-MET were used. The measurements 

of solar radiation from site WMR-MET were used in the 

simulations, with the exception of a few days prior to June 15, 

2005, when the data from site WMR-MET were not yet 

available; solar data from the AgriMet site at Agency Lake 

(Bureau of Reclamation, 2008) were used for those days. 

Calibration Data

Surface Elevation

USGS water-level gage data collected at three sites 

around the lake were used for water-level calibration. The sites 

were Upper Klamath Lake at Rocky Point (site ID 11505800), 

Upper Klamath Lake at Rattlesnake Point (site ID 11505900), 

and Upper Klamath Lake at Klamath Falls (site ID 11507000, 

fig. 1).

Velocity

Data used for calibration of water velocity came from five 

acoustic Doppler current profilers (ADCPs) deployed in the 

lake in 2005 and two ADCPs deployed in 2006 (fig. 1). These 

instruments and the velocity measurements collected have 

been described elsewhere (Gartner and others, 2007). Halfhourly velocity data were collected from the ADCPs between 

June 21 and September 12 in 2005, and between May 24 

and September 25 in 2006. On July 5, 2006, the ADCP 

located at site ADCP1 was dragged by a barge from its initial 

deployment location to a location approximately 280 m away, 

where the meter continued to collect data.

Water Temperature

The data used for water temperature calibration in this 

report have been described elsewhere (Hoilman and others, 

2008). These data were measured at 13 sites (1 m from the 

top and 1 m from the bottom at 3 sites) in 2005 and 12 sites 

(1 m from the top and 1 m from the bottom at 4 sites) in 

2006 (fig. 1). At 3 of these sites in 2005 and 4 of these sites 

in 2006, a water-quality monitor was placed at 1 m from the 

surface and 1 m from the bottom, resulting in a total of 16 

water-quality monitors in both years. At sites where only 

one monitor was used, it was placed at 1 m from the bottom 

or at the middle of the water column if the total depth was 

2 m or less. Hourly data were collected from the 16 waterquality monitors between June 3 and October 11, 2005, and 

between May 15 and October 16, 2006. Five additional sites 

were added within 100 m of the shoreline in 2006, starting on 

June 14 (fig. 1, table 2).

Other Datasets

The extinction coefficient used in the model was 

determined from weekly vertical profiles of photosynthetically 

active radiation (PAR) at six water-quality monitoring sites. 

Light intensity was measured from the water surface to 2.5 m 

depth (or the lake bottom) in 0.5 m increments with a LiCor 

LI-193 spherical quantum sensor. These measurements 

were used to estimate the extinction coefficient by fitting an 

exponential relation based on Beer’s Law to the vertical profile 

of the light meter readings. The methods are described in more 

detail by Hoilman and others (2008). A mean value of the 

extinction coefficient was used in the simulations presented 

in this report, which was obtained as a spatial and temporal 

average of the measured values for July and August of 2005 at 

all six sites. 

The vertical profiles of light intensity were collected to 

support a weekly program of light and dark bottle incubations 

in 2005 and 2006. Calculations based on some of these 
16 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality Table 2. Summary of the dates of availability of meteorological, velocity and water temperature data in Upper Klamath Lake between May and October of 2005 and 2006; summary of the simulation periods in 2005 and 2006, and the dates over which error statistics were calculated in 2005 and 2006. May June July August September October Meteorological stations BLB-MET HDB-MET WMR-MET SSHR-MET MDN MDL Acoustic Doppler Current Profilers ALL Water-temperature monitors ALL Simulations UW-MDL UW-MDN VW Error calculations Midsummer Meteorological stations BLB-MET HDB-MET WMR-MET SSHR-MET MDN MDL Acoustic Doppler Current Profilers ADCP1a ADCP1b ADCP9 Water temperature monitors ALL Water temperature monitors: Nearshore ALL Simulations VW Error calculations

Full seasonMidsummer

Air temperature sensor failure 

Instrument dragged from original location

2005 2006
Simulation Results 17

data that were collected in 2006 are presented below in 

the discussion of the role of currents in determining water 

quality in the lake. The values of production and consumption 

obtained from the bottle incubations were converted to a 

change in dissolved oxygen concentration over a 24-hour 

period as detailed in Appendix A. Further details of datacollection methods are provided by Hoilman and others 

(2008).

The single value for the fraction of the sky covered by 

clouds that was used in the model simulations is the long term 

July-August average at the Medford, Oregon, airport (Oregon 

Climate Service, 2008). Because Medford is located about 

130 km to the west of Klamath Falls and at about 830 m lower 

elevation, the daily values of cloud cover collected there are 

not expected to be accurate for Klamath Falls. The months of 

July and August in Klamath Falls are generally quite clear, 

but weather patterns in the spring and fall can result in periods 

of several cloudy days at a time. However, cloud cover is 

used only in the calculation of incoming longwave radiation 

(eq. 24), and the errors associated with using a constant value 

are expected to be small.

Simulation Results 

Hydrodynamics

The hydrodynamic model was calibrated for 2005 by 

adjusting parameters a and z0, which determine the surface 

and bottom friction coefficient, respectively (table 3). The 

same values for a and z0 were used to run the model with the 

three different types of wind (UW-MDL, UW-MDN, and 

VW). In the UW-MDL simulation, a spatially uniform wind 

as measured at site MDL, a raft site in the southern part of 

the lake, was used. In the UW-MDN simulation, a spatially 

uniform wind as measured at site MDN, a raft site in the 

northern part of the lake, was used. In the VW simulation, a 

spatially variable wind was used. The spatially variable wind 

was determined from an interpolation of six sites, including 

the two raft sites and four others—WMR-MET, BLB-MET, 

HDB-MET, and SSHR-MET—located on the shoreline 

(fig. 1), with the exception that the data from site SSHR-MET 

were not included in the interpolation until August 18 when 

the sensors at that site started collecting data, and the data 

from site MDL were not included in the interpolation between 

September 1 and September 8 because the sensors at that site 

failed. The availability of data from the various meteorological 

sites and the periods covered by the three simulations were 

discussed previously and are summarized in table 2. In all 

cases the model was run for 5 days or more prior to the start of 

any error calculations in order to allow oscillations resulting 

from the startup of the model to dissipate. 

The gaged and simulated lake elevation at two sites, 

UKL at Rocky Point, near the northern end of the lake, and 

UKL at Klamath Falls, near the southern end, are shown in 

fig. 6. These two gages show that the lake surface rises at the 

southern end and falls at the northern end on a daily basis, in 

response to the diel wind fluctuation. The model reproduces 

the timing of the oscillations well, although the simulated 

amplitude is larger than observed at the gages, particularly at 

the Klamath Falls gage, as can be seen in detail from July 15 

to July 30. The VW simulation results in simulated oscillations 

closer to the observed amplitude than the UW-MDL 

simulation (fig. 6). 

Prevailing winds over UKL blow over the northern 

part of the lake from the west to slightly northwest and are 

constrained by the surrounding topography to a northwest 

wind over the southern two-thirds of the lake, where the wind 

sensor at site MDL was located (Hoilman and others, 2008, 

fig. 7). On a typical day, wind speed picks up in the early to 

midafternoon and then dies down in the late evening or early 

morning. Typically, early morning is the calmest time of day. 

The response of the currents to moderate, prevailing winds is 

a relatively well-defined clockwise circulation consisting of 

broad, shallow flow in the direction of the wind on the eastern 

side of the lake and passing to the east of Bare Island, and 

a narrow, deep flow opposing the wind through the trench 

along the western shoreline and passing to the west of Bare 

Island. Usually, a large fraction of the northward flow in the 

trench veers eastward both south and north of Bare Island, 

forming a gyre that circulates water between Bare Island and 

Rattlesnake Point. The rest of the northward flow continues 

into the northern part of the lake, completing a clockwise 

loop that circulates water between Rattlesnake Point and the 

northernmost part of UKL (fig. 8). A reversal in the wind 

direction that is sustained for several hours causes the currents 

in the trench to reverse, resulting in a stall in the circulation 

pattern. If the wind reversal lasts long enough (at least a day), 

the current loop between Bare Island and Rattlesnake Point 

can reverse as well (fig. 9). 

All of the ADCP data collected during 2003–06 are 

described in Gartner and others (2007), but it is useful to 

describe the basic circulation pattern in the lake in order to 

provide context for the ADCP data used in calibrating the 

model. Site ADCP1 was located in the deepest part of the 

trench (fig. 1). The current speeds are higher in the trench 

than in the shallower parts of the lake; site ADCP1 was 

strategically located to capture some of the highest velocities 

in the lake (fig. 10). The direction of the currents at site 

ADCP1 is constrained by the bathymetry of the trench to 

approximately 310 degrees under prevailing wind conditions; 

during a wind reversal the currents change direction by 

180 degrees (fig. 11). The other ADCPs also were in areas that 

have pronounced prevailing current directions, although not as 

pronounced as at site ADCP1. Site ADCP7 was located near 
18 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

the southern terminus of the trench at the mouth of Howard 

Bay. Current speeds there were the second highest of the five 

locations (fig. 10), and current directions were predominantly 

about 30 degrees east of due north, in alignment with the 

direction of the trench at that site. Sites ADCP3 and ADCP5 

were located in shallower water in the northern third of the 

lake: site ADCP5 in a shallower section of the trench to the 

west of Eagle Point, and ADCP3 near the terminus of the 

trench north of Ball Bay (fig. 1). The current directions at both 

of these sites are not as tightly constrained by the bathymetry 

but still show a prevailing direction that indicates currents 

moving in a clockwise direction around the northern third of 

the lake (fig. 11). Current speeds at these sites are much lower 

than in the deeper sections of the trench (fig. 10). Site ADCP6 

was located in the broad, shallow flow that moves over the 

eastern part of the lake; the current direction there was the 

most variable of all the sites, but it indicates flow primarily to 

the southeast (fig. 11). Current speeds at site ADCP6 were the 

lowest of those measured at the five sites (fig. 10). 

Error statistics for the three simulations (UW-MDL, 

UW-MDN, and VW) for the midsummer overlap period 

between July 26 and August 31 are provided in table 4 and 

are presented visually for ADCP1 in fig. 12. Error statistics 

were calculated at each ADCP site for both the depth-averaged 

speed, which is calculated from all the bins in the water 

column, and for individual east-west and north-south velocity 

Table 3. List of inputs to the hydrodynamic and heat transport model of Upper Klamath Lake, Oregon.

[How determined: m indicates a measured quantity; c indicates a calibration parameter; l indicates a value determined 

from literature; p indicates a physical constant. Abbreviations: m, meter; s, second; J, joules; g, gram; °C, degrees 

Celsius; mb, millibar; W, Watt; kg, kilogram; N, Newton; °K, degrees Kelvin; –, no units]

Symbol Name Value Units How 

determined

a Constant in equation for surface drag coefficient 0.00056 (m s-1)-0.5 c

CB Bowen ratio empirical coefficient 0.61 mb °C -1 3

1

cp Specific heat capacity of water at 20°C 4.186 J g-1 oC-1 p

c0 Constant used to relate turbulent kinetic energy to 

rate of turbulent kinetic energy dissipation

0.5544 – 1l

Dh Horizontal turbulent diffusivity coefficient 0.01 m2

 s-1 2l

g Gravitational acceleration 9.81 m s-2 p

HSW Incoming shortwave radiation Variable W m-2 m

Kh Horizontal turbulent viscosity coefficient 0.10 m2

 s-1 2l

ke Extinction coefficient 2.9 m-1 m

Lw Latent heat of evaporation 2.5 × 106 J kg-1 p

M Fraction of the sky covered by clouds 0.23 – m

pa Atmospheric pressure 871.4 mb m

psl Atmospheric pressure at sea level 1,013.25 mb p

rh Relative humidity Variable – m

rLW Fraction of incident longwave radiation reflected 

at the water surface

0.02 – c

rSW Fraction of incident shortwave radiation reflected 

at the water surface

0.05 – c

Ta Air temperature Variable °C m

r

ui

Wind vector at i m above the surface Variable m s-1 m

z0 Bottom roughness height 0.0005 m c

α0 Proportionality constant 0.937 × 10-5 °K-2 4

l

w ε Emmisivity of water 0.97 – p

κ von Karman constant 0.41 – p

ρ0 Reference of density of water 1,000 kg m-3 p

ρa Density of air (dry air at 20°C and Pa) 1.036 kg m-3 p

σ Stefan-Boltzmann constant 5.67 ×10-8 W m-2 °K-4 p

1

Warner and others (2005).

2

Bowie and others, 1985 and Lerman (1971).

3

Bowen (1926).

4

Martin and McCutcheon (1998), Swinbank (1963), and Wunderlich (1972).
Simulation Results 19

OR19_0141_UKL_fig_06

AREA COVERED

IN GRAPH ABOVE

15 15

10 14

1 1 15

30 1 5

1 15 1

JUNE JULY

AUGUST

AUGUST

JULY

SEPTEMBER

1 14

OCTOBER

2005

1,261.6

1,261.8

1,262.0

1,262.2

1,262.4

1,262.6

1,262.8

1,263.0

ELEVATION OF WATER SURFACE, IN METERS NGVD29

Klamath Falls

Rocky Point

Klamath Falls—simulated

Rocky Point—simulated

1,262.0

1,262.2

1,262.4

1,262.6

UNIFORM WIND FORCING VARIABLE WIND FORCING

Figure 6. Hourly values of lake elevation at Upper Klamath Lake, Oregon, measured at two gages, 2005. 
20 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

OR19_0141_UKL_fig_07

0

5

10

15

0

5

10

15

15 1 15 1 15 1 15 1 15 1 14

MAY JUNE JULY AUGUST SEPTEMBER OCTOBER 

2006

15 1 15 1 15 1 15 1 15 1 14

MAY JUNE JULY AUGUST SEPTEMBER OCTOBER 

2005

SPEED, IN METERS

PER SECOND

100

200

300

400

500

DIRECTION SPEED, IN METERS

PER SECOND

100

200

300

400

500

DIRECTION

Figure 7. Wind speed and direction at site MDL in Upper Klamath Lake, Oregon. The wind direction ranges from 100 

to 460 degrees because 360 degrees have been added to values of wind direction less than 100 degrees to avoid rapid 

numerical changes when the wind direction moves from east to west of north, and vice versa. The wind direction is 

referenced clockwise from north, and is the direction from which the wind is blowing.

components at two points in the water column—the bin closest 

to 1 m from the bottom and the bin closest to 1 m from the 

surface. The mean error (ME) in the simulated currents at 

all of the observation sites is provided as a measure of the 

overall bias of the simulation at each site. Errors are calculated 

as observed value minus simulated value, so a positive ME 

indicates simulated values less than observed values. The root 

mean squared error (RMSE) is provided as a measure of the 

overall goodness-of-fit of the simulation to the observations. 

The mean of the half-hourly observations and simulations also 

is provided for comparison. During the calibration process, 

emphasis was placed on correctly simulating the highest 

velocities through the trench, and as a result the ME and 

RMSE at site ADCP1 are small in comparison to the depthaveraged speed and the individual velocity components at that 

site (table 4, row 1 and rows 6–9, and fig. 10). In addition, 

both the ME and RMSE at ADCP1 are smaller for the VW 

simulation than for either of the uniform wind simulations, 

particularly in comparison to the UW-MDN simulation, which 

used only the wind measured in the northern part of the lake. 

The second-best fit appears to be at site ADCP6 (table 4, 

row 4), where velocities are low. Comparison of the individual 

velocity components table 4, rows 18–21) indicates that, while 

the speed at this site is captured well by the VW simulation, 

the direction of the simulated currents is rotated to the west of 

the observations (fig. 11). At site ADCP6, the VW simulation 

produces a better result than the UW-MDL simulation, even 

though site ADCP6 was coincident with the meteorological 

site at MDL.

The simulated mean current speeds at the remaining 

sites are biased low in the VW simulation, as indicated by a 

positive ME, particularly at sites ADCP5 and ADCP7 (table 4, 

rows 3 and 5) and to a lesser extent at site ADCP3 (table 4, 

row 2). The ME of the average current speed at sites ADCP5 

and ADCP7 is less for the UW-MDL simulation than for the 

VW simulation (table 4, rows 3 and 5); at site ADCP5 the 

RMSE also is less for the UW-MDL simulation than for the 

VW simulation (table 4, row 3). Site ADCP5 is located in 

the northern part of the lake; yet, if a uniform wind forcing is 

used, the currents there are better simulated by using winds 

measured in the central part of the lake than in the northern 

part of the lake closer to the site. This underscores that fact 

that success in accurately simulating the currents anywhere in 

the lake is dependent on having an accurately measured wind 

over as much of the lake as possible, rather than an accurate 

wind measured nearby at a single point of interest. In the case 

that wind data are available at only a single point, site MDL is 

a better site for collecting wind to force the model because the 
Simulation Results 21

Figure 8. Simulated depth-averaged currents under prevailing 

(northwest) wind conditions, Upper Klamath Lake, Oregon. 

Figure 9. Simulated depth-averaged currents under reversed 

(southwest) wind conditions, Upper Klamath Lake, Oregon. 

wind data collected there probably represents a large surface 

area in the broad, central part of the lake where much of the 

momentum transfer from the wind to the water occurs; site 

MDN is not as good a choice because the wind data collected 

there represent a smaller surface area of the lake.

The ME statistics for the east-west and north-south 

velocity components (table 4, rows 6–25) are mixed with 

regard to whether the VW or the UW-MDL scenario 

produces the best simulation of the observations. The RMSE 

statistics, however, show consistent improvement with the 

VW simulation, with the exception of the near-bottom, 

north-south velocity component at site ADCP6 (table 4, row 

20). Also apparent in the statistics (particularly the RMSE) 

for the individual components is the fact that the errors are 

consistently larger at the surface than at the bottom. 

At site ADCP3, for example, the simulated north-south 

component of the surface currents has a tendency to be too 

large in the early afternoon, whereas the simulated east-west 

component has a tendency to be too small in the late afternoon 

(fig. 13 and table 4). The VW simulation reproduces the 

directional shift between the bottom and the surface better 

than the UW-MDL simulation (fig. 13). 

At sites ADCP1, ADCP6, and ADCP7, the comparison 

between observed and simulated currents worsens from 

September 1 to September 8 (fig. 10). Because there is a data 

gap in the wind observations at site MDL during this time, 

a likely explanation is that the interpolated wind field is still 

dependent on observations at site MDL to provide accuracy in 

the wind forcing over the lake.

OR19_0141_UKL_fig09

Date: August 26, 2005 00:00

0 1 2 3 4 5 MILES

0 1 2 3 4 5 KILOMETERS

EXPLANATION

Length equivalent to

10 centimeters per second

L_fig08

0 1 2 3 4 5 MILES

0 1 2 3 4 5 KILOMETERS

EXPLANATION

Length equivalent to

10 centimeters per second

Date: July 29, 2005 04:00
22 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

Figure 10. Observed 

and simulated current 

speeds at Acoustic 

Doppler Current Profiler 

(ADCP) sites in Upper 

Klamath Lake, Oregon, 

2005. The error bar 

placed on graphs of 

speed indicates the 

root mean square error 

(RMSE) between the 

observed and simulated 

currents over the length 

of the simulation.

Figure 11. Observed 

and simulated current 

direction at Acoustic 

Doppler Current Profiler 

(ADCP) sites in Upper 

Klamath Lake, Oregon, 

2005. The current 

direction ranges from 

100 to 460 degrees 

because 360 degrees 

have been added 

to values of current 

direction less than 

100 degrees to avoid 

rapid changes when the 

wind direction moves 

from east to west of 

north, and vice versa. 

The current direction is 

referenced clockwise 

from north, and is the 

direction to which the 

current is moving.

R19_0141_TWood_fig_10

Observations

Simulations

10

0

20

30

40

50

10

0

20

30

40

50

10

0

20

30

40

50

10

0

20

30

40

50

10

0

20

30

40

50

15 1 15 1 15 1 14

JUNE JULY AUGUST SEPTEMBER

UNIFORM WIND FORCING VARIABLE WIND FORCING

2005

CURRENT SPEED, IN CENTIMETERS PER SECOND

ADCP3

ADCP5

ADCP6

ADCP7

ADCP1 RMSE RMSE

RMSE RMSE

RMSE RMSE

RMSE RMSE

RMSE RMSE

CURRENT DIRECTION, DEGREES CLOCKWISE FROM NORTH

Observations

Simulations

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

15 1 15 1 15 1 14

JUNE JULY AUGUST SEPTEMBER

UNIFORM WIND FORCING VARIABLE WIND FORCING

2005

ADCP3

ADCP5

ADCP6

ADCP7

ADCP1
Simulation Results 23

The model was validated using 2006 boundary 

conditions and forcing functions from May 20 to October 

15. The values of the two calibration parameters, a and 

z0, were the same as established during calibration for the 

2005 season. This provided an opportunity to validate the 

performance of the model over nearly 5 months. The endof-summer lake elevation in 2006 was lower than that in 

2005 and revealed some inaccuracies in the representation 

of the bathymetry and shoreline in the model, as the 

simulated gages were without water after mid-September 

(fig. 14).

The observed and simulated current speed and 

direction at the two 2006 ADCP sites are shown in fig. 15. 

The visual comparison suggests that the underestimation 

of current speeds by the model simulation at site ADCP1 

is greater in 2006 than in 2005, and this is confirmed by 

the ME statistics and the means (table 5, rows 1 and 2). 

In order to check whether the 2006 errors were larger 

-20

0

-10

20

10

-20

0

-10

20

10

-20

0

-10

20

10

-20

0

-10

20

10

-20

0

-10

20

10

-20

0

-10

20

10

CURRENT SPEED, IN CENTIMETERS PER SECOND

1 2 3 4 5 6 7 8 9 10

AUGUST 2005

OBSERVED EAST COMPONENT

OBSERVED NORTH COMPONENT

SIMULATED EAST COMPONENT- VW

SIMULATED NORTH COMPONENT- VW

SIMULATED EAST COMPONENT- UW-MDL

SIMULATED NORTH COMPONENT- UW-MDL

Bottom

Surface

Figure 13. Observed and simulated east-west and north-south current components at site ADCP3 in Upper Klamath 

Lake, Oregon, 2005.

Figure 12. Goodness-of-fit velocity statistics at site ADCP1 in 

Upper Klamath Lake, Oregon, 2005 and 2006.

OR19_0141_UKL_fig_12ab

MEAN ERROR (ME) IN CENTIMETERS PER SECOND

2005 VW

2005 UW-MDL

2005 UW-MDN

2006 VW

-3

-4

1

0

-1

-2

2

4

3

5

ADCP1 ADCP1

-u(b)

ADCP1

-u(s)

ADCP1

-v(b)

ADCP1

-v(s)

ROOT MEAN SQUARED ERROR (RMSE), 

IN CENTIMETERS PER SECOND

0

4

3

2

1

5

7

6

8

ADCP1 ADCP1

-u(b)

ADCP1

-u(s)

ADCP1

-v(b)

ADCP1

-v(s)
24 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

The simulation of the currents at the second 

2006 ADCP site was less successful. This site, 

ADCP9, was located in the northern part of 

the lake, farther from the trench than either site 

ADCP3 or ADCP5 in 2005. The mean depthaveraged speed at this site during August was 

3.42 cm/s (table 5), which was less than the speed 

at either of the northern sites in 2005 (table 4). 

The RMSE of 2.49 for the July 26 to August 31 

period (table 5, row 3) is comparable to the errors 

at the northern sites in 2005 (table 4, rows 2 and 

3), but a visual examination of the data shows 

that the simulation does not reproduce well the 

fluctuations in the observations, particularly on the 

time scale of several days (fig. 15). Site ADCP9 

was located in a place where the currents were 

particularly complicated but not strong. The mean 

direction of the currents was toward 280–290 

degrees clockwise from north at the surface, 

rotated about 20 degrees counterclockwise from 

the bottom (Gartner and others, 2007), so at this 

site more than at the northern ADCP sites in 2005, 

the currents directly opposed the prevailing wind 

direction observed at site MDN, where the strongest 

winds were approximately westerly, coming from 

270–280 degrees clockwise from north (Hoilman 

and others, 2008). Thus, errors in the wind 

forcing at the site, as well as inaccuracies in the 

characterization of surface or bottom friction, will 

result in more noticeable discrepancies between 

observations and simulations at this site. 

When putting the error statistics into 

perspective, it is important to understand that 

although ADCP technology is generally reliable 

and accurate for measuring water-current 

profiles, there are some inherent limitations and 

assumptions associated with those measurements, 

including two primary assumptions. First, since 

the instrument determines beam velocity from 

Doppler (or frequency) shift in the acoustic beams 

from backscattered signals, it is actually measuring 

the motion of particles in the water, not the actual 

water motion, which may be different. Second, 

the instrument transforms the beam velocities into 

earth coordinates, which requires the assumption 

that flow is essentially homogeneous in the 

four beams. There are error checks to discard 

obviously bad measurements but the assumption of 

homogeneous flow may impart some error. Because 

the four acoustic beams are oriented 20 degrees 

from vertical, the validity of that assumption is a 

function of distance from the transducer. Further, the 

instrument has a limited accuracy for each acoustic 

measurement. That accuracy (or “single ping” 

because of the longer simulation period in 2006, errors also were 

calculated over the same period (July 26 through August 31) that was 

used for the error calculations in 2005. Errors over this shorter period 

were comparable to those over the longer simulation, which seems to 

suggest that errors are not growing with the length of the simulation 

period. Nonetheless, the 2006 errors were larger for the July 26 through 

August 31 period than for the same period in 2005. For example, the 

RMSE at site ADCP1 was 3.88 for July 26 to August 31, 2006, as 

compared to 3.08 for the same period in 2005. It is not obvious why 

the model underestimates the currents in 2006, but it indicates that the 

calibration obtained for the 2005 data is not necessarily optimized for 

multiple years.

Table 4. Goodness-of-fit velocity statistics for the UnTRIM model of Upper 

Klamath Lake, Oregon, 2005. 

[Errors are calculated half-hourly as observed value-simulated value. Site modifier: -u 

indicates east/west component of velocity; -v indicates north/south component of velocity; 

(b) indicates 1 meter from bottom, (s) indicates 1 meter from surface. Mean and Root 

Mean Square: VW indicates that the model was run using spatially variable wind forcing; 

UW indicates that the model was run using spatially uniform wind forcing as observed 

at the site indicated (MDN or MDL). Abbreviations: cm/s, centimeter per second; °C, 

degrees Celsius]

Row Site

Mean of 

observations

(°C)

Mean error

(cm/s) 

Root mean 

square error 

(cm/s)

VW UWMDL

UWMDN VW UWMDLUWMDN

Depth-averaged speed

1 ADCP1 11.23 0.50 -0.73 3.73 3.08 3.35 5.65

2 ADCP3 4.08 .80 -.48 1.10 2.04 2.20 2.41

3 ADCP5 5.81 2.58 1.89 1.95 3.77 3.52 3.80

4 ADCP6 4.44 -.07 2.54 2.08 1.77 3.25 3.14

5 ADCP7 8.38 3.02 1.66 5.01 4.38 4.60 7.11

Velocity components

6 ADCP1-u(b) -6.15 0.54 1.74 -1.15 2.80 3.70 3.93

7 ADCP1-u(s) -6.05 -.32 -.25 -3.82 4.19 5.23 7.11

8 ADCP1-v(b) 6.44 -1.02 -2.73 .96 3.32 4.98 4.10

9 ADCP1-v(s) 5.96 2.00 4.20 4.06 4.75 7.02 6.84

10 ADCP3-u(b) -1.03 .29 1.16 .44 2.21 2.66 2.55

11 ADCP3-u(s) -.40 .79 1.02 .79 4.73 5.34 4.61

12 ADCP3-v(b) 3.87 .63 -.84 1.47 2.15 2.78 2.67

13 ADCP3-v(s) 1.98 1.13 1.82 1.25 3.76 5.37 3.62

14 ADCP5-u(b) -3.51 -1.84 -.65 -1.34 3.22 3.54 3.49

15 ADCP5-u(s) -3.74 -1.26 -1.36 -1.85 3.78 4.63 4.45

16 ADCP5-v(b) -2.30 -.99 -1.49 -.63 2.66 3.14 2.70

17 ADCP5-v(s) -4.52 -.67 2.16 -.69 3.44 5.04 3.32

18 ADCP6-u(b) 1.50 .57 1.49 2.47 2.19 2.70 3.46

19 ADCP6-u(s) 6.07 5.05 5.82 6.85 8.10 8.62 9.97

20 ADCP6-v(b) -1.28 3.57 .28 .14 4.44 2.90 2.54

21 ADCP6-v(s) -8.55 -.48 -3.48 -5.44 5.61 6.96 8.96

22 ADCP7-u(b) 3.74 .80 .40 2.05 2.64 3.32 4.19

23 ADCP7-u(s) 5.52 2.72 .48 3.17 6.31 6.71 8.02

24 ADCP7-v(b) 7.82 2.82 1.63 6.77 4.37 4.82 9.25

25 ADCP7-v(s) -1.23 -2.74 -3.23 -1.04 6.12 7.03 6.48
Simulation Results 25

Figure 14. Hourly values of lake elevation at Upper Klamath Lake, Oregon, measured at two gages, 2006. 

OR19_0141_UKL_fig_14

AREA COVERED

IN GRAPH ABOVE

15 15

25 29

1 1 15

15 20

1 15 1

JUNE JULY AUGUST

JULY

SEPTEMBER

1 14

OCTOBER

2006

1,261.6

1,261.4

1,261.8

1,262.0

1,262.2

1,262.4

1,262.6

1,262.8

1,263.0

ELEVATION OF WATER SURFACE IN METERS ABOVE NGVD29

Klamath Falls

Rocky Point

Klamath Falls–simulated

Rocky Point–simulated

1,261.8

1,262.0

1,262.2

1,262.4
26 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

OR19_0141_UKL_fig_15

0

20

40

10

30

50

0

20

40

10

30

50

15 1 15 1 15 1 15 1 15 30

MAY JUNE JULY AUGUST SEPTEMBER

2006

15 1 15 1 15 1 15 1 15 30

MAY JUNE JULY AUGUST SEPTEMBER

2005

SPEED, IN 

CENTIMETERS

PER SECOND

100

200

300

400

500

DIRECTION

100

200

300

400

500

DIRECTION

ADCP1

ADCP1

ADCP9

ADCP9

Observations

EXPLANATION

Simulations

RMSE RMSE

RMSE

SPEED, IN 

CENTIMETERS

PER SECOND

LOCATION

CHANGED

Figure 15. Observed and simulated current speed and direction at two Acoustic Doppler Current Profiler (ADCP) sites 

in Upper Klamath Lake, Oregon, 2006. The error bar placed on graphs of speed indicates the root mean square error 

(RMSE) between the observed and simulated currents over the length of the simulation. The current direction ranges 

from 100 to 460 degrees because 360 degrees have been added to values of current direction less than 100 degrees to 

avoid rapid numerical changes when the current direction moves from east to west of north, and vice versa. The current 

direction is referenced clockwise from north and is the direction to which the current is moving.

standard deviation) is typically reduced by averaging multiple 

pings to a single ensemble average that is then recorded. 

The (random) standard deviation is reduced by a factor of 

the square root of the number of pings. In the case of the 

ADCP measurements in UKL, the ADCP ensemble standard 

deviation or accuracy was about 0.73 cm/s for the horizontal 

velocity.

Overall, the calibration of the hydrodynamic model for 

2005 and the validation for 2006 was successful. In both years, 

the biggest discrepancies between the observed and simulated 

currents occurred near the surface of the water column and 

indicated difficulty with simulating the friction boundary 

layer, so the errors in the depth-averaged speed were larger at 

shallow sites where the boundary layer is a bigger proportion 

of the water column. The simulation of the highest velocities 

through the trench was accurate, however, and provides 

confidence that the overall circulation patterns produced are 

correct. The validation of the model in 2006 produced currents 

that were somewhat underestimated in the trench, however, 

which indicated that even though the currents were accurately 

simulated at the same location in 2005 by adjusting only 

two calibration parameters, the model was not necessarily 

optimized for multiple years. 
Simulation Results 27

Table 5. Goodness-of-fit velocity statistics for the UnTRIM model of Upper 

Klamath Lake, Oregon, 2006. 

[Errors are calculated half-hourly as observed value-simulated value. A spatially variable 

wind was used in the simulation. Full Season at ADCP1a is from May 24 to July 18 when 

the instrument was moved to a new position; full season at ADCP1b is from July 19 through 

September 25 at the new position; full season at ADCP9 is May 24 to September 25. Midsummer is from July 26 to August 31. Site modifier: -u indicates east/west component of 

velocity; -v indicates north/south component of velocity; (b) indicates 1 meter from bottom, 

(s) indicates 1 m from surface. Abbreviations: cm/s, centimeter per second; –, no data]

Row Site

Mean of 

observations

(cm/s)

Mean error

(cm/s)

Root mean 

square error 

(cm/s)

Full 

season

MidsummerFull 

season

MidsummerFull 

season

Midsummer

Depth-averaged speed

1 ADCP1a 10.13 – 1.32 – 3.92 –

2 ADCP1b 12.42 12.98 2.17 2.30 3.83 3.88

3 ADCP9 3.65 3.42 .60 0.44 2.49 2.50

Velocity components

4 ADCP1a-u(b) -6.01 – -0.46 – 3.44 –

5 ADCP1a-u(s) -5.16 – -.06 – 5.20 –

6 ADCP1a-v(b) 5.87 – .80 – 4.08 –

7 ADCP1a-v(s) 5.70 – 2.91 – 6.91 –

8 ADCP1b-u(b) -6.36 -6.63 -.13 -0.08 2.82 2.84

9 ADCP1b-u(s) -4.42 -4.84 1.61 1.41 4.92 4.72

10 ADCP1b-v(b) 8.51 8.89 1.81 1.92 3.85 3.92

11 ADCP1b-v(s) 4.65 5.71 1.17 1.96 4.95 4.86

12 ADCP9-u(b) -3.14 -3.05 -.22 -.43 2.53 2.50

13 ADCP9-u(s) -2.09 -2.08 -.28 .00 5.97 4.84

14 ADCP9-v(b) .34 .31 -.80 -.96 2.55 2.50

15 ADCP9-v(s) -.48 -.61 .62 .46 4.45 4.03

Heat Transport

The heat transport model was calibrated by adjusting one parameter, 

rSW, the amount of incoming solar radiation reflected at the water surface. 

Comparisons between the observed and simulated temperature at the 

shallow sites where there was one continuous monitor, placed 1 m from 

the bottom, are provided in fig. 16. The goodness-of-fit temperature 

statistics for the midsummer overlap period of July 26 to August 31 are 

provided in table 6 for all three simulations and are presented visually 

for a subset of the sites in fig. 17. The absolute value of the ME of the 

VW simulation is less than 1°C at all sites and less than the ME of the 

UW-MDL simulation at most sites. The RMSE of the VW simulation 

also is less than 1°C at most of the sites, has a maximum value of 1.12, 

and is less than the RMSE of the UW-MDL simulation at most sites. 

The ME statistics are both positive (indicating an underestimation of 

the temperature) and negative (indicating an overestimation of the 

temperature), but the negative errors are larger, so that there a slight 

overestimation of the overall temperature around the lake during August, 

which is the warmest month. Interestingly, the RMSE of the UW-MDN 

simulation is generally comparable to and sometimes less than the RMSE 

of the VW simulation, indicating that the correct simulation of the 

currents is not required for correct simulation 

of the temperature. This underscores the fact 

that the temperature at any point in the lake 

is primarily a function of the meteorological 

forcing at the surface rather than advection. 

After September 4, there is a noticeable 

increase in the diel range in temperature at 

some of the sites, which is due in part to 

the use of a land-based measurement of air 

temperature (from WMR-MET) to calculate 

incoming longwave radiation and evaporative 

heat loss for the boundary condition after 

the measurement of air temperature from the 

raft at site MDL became unavailable. The 

diel range in temperature over land is greater 

than that over water, primarily because the 

daily minimum is lower (fig. 18); the use 

of the measurement made over the water is 

preferable and generates more realistic results. 

Rapid drops in air temperature due to weather 

patterns such as those occurring around June 

17, September 22, and October 1 appear to 

produce a more extreme response in simulated 

water temperature than is realistic. During these 

times, cloud cover is likely underestimated in 

the model by the use of a uniform value, and the 

incoming longwave radiation is therefore likely 

underestimated as well. Thus the incorporation 

of a temporally variable cloud cover, or directly 

measured longwave radiation, would likely 

moderate the simulated temperature response of 

the lake to these weather patterns. 

Another important aspect of temperature 

simulation, in addition to correct simulation of 

the daily average temperature and temperature 

range, is the ability to predict the difference 

between the upper and the lower parts of the 

water column at those sites where the water 

is deep enough to stratify, including sites 

MDT, EPT, and MDN (Hoilman and others, 

2008). Comparisons between the observed and 

simulated temperature records at the deep sites 

where there was a near-surface and near-bottom 

monitor in 2005 are provided in fig. 19. The 

observed and simulated daily maximum and 

daily minimum range in temperature from top 

to bottom at sites MDT, EPT, and MDN are 

shown in fig. 20. In general the VW simulation 

seems to be able to correctly reproduce those 

periods following July 26 when there is some 

temperature stratification that persists through 

several days; these are most noticeable at 

site MDT, the deepest site. The UW-MDL 

simulation does not reproduce those periods 

of persistent temperature stratification prior to 
28 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

OR19_0141_TWood_fig_16

Observations

Simulations 10

5

20

30

10

5

20

30

10

5

20

30

10

5

20

30

10

5

20

30

10

5

20

30

10

5

20

30

10

5

20

30

10

5

20

30

10

5

20

30

1 15 1 15 1 15 1 15 1 14

JUNE JULY AUGUST SEPTEMBER OCTOBER 

RPT

WMR

FBS

HDB

NBI

BLB

MDL

EHB

SHB

MPT

UNIFORM WIND FORCING VARIABLE WIND FORCING

WATER TEMPERATURE, IN DEGREES CELSIUS 

2005

Figure 16. Observed and simulated temperatures at 10 shallow sites in Upper Klamath Lake, Oregon, 2005.
Simulation Results 29

Figure 17. Goodness-of-fit temperature statistics at 13 sites in Upper Klamath Lake, Oregon, 2005 and 2006.

OR19_0141_UKL_fig_17ab

-0.4

-0.2

-0.6

-0.8

0.0

0.8

0.6

0.4

0.2

1.0

1.2

1.6

1.8

1.4

2.0

MDT

-L

MDT

-U

EPT

-L

EPT

-U

MDN

-L

MDN

-U

BLB EHB FBS HDB NBI RPT WMR

MEAN ERROR (ME) IN DEGREES CELSIUS

2005 VW

2005 UW-MDL

2005 UW-MDN

2006 VW

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

ROOT MEAN SQUARED ERROR (RMSE),

IN DEGREES CELSIUS

MDT

-L

MDT

-U

EPT

-L

EPT

-U

MDN

-L

MDN

-U

BLB EHB FBS HDB NBI RPT WMR

Table 6. Goodness-of-fit temperature statistics for the UnTRIM model of Upper Klamath Lake, 

Oregon, 2005. 

[Errors are calculated hourly as observed value-simulated value. Site modifier: -L indicates monitor 1 meter 

from the bottom, -U indicates monitor 1 meter from the surface. Site names in boldface type indicate that 

temperature was measured at the site during both 2005 and 2006. Mean and Root Mean Square Error: VW 

indicates model run using spatially variable wind forcing; UW indicates model run using spatially uniform 

wind forcing as observed at the site indicated (MDN or MDL). Overall mean was calculated based on sites in 

boldface type only. Abbreviations: °C, degrees Celsius]

Row Site

Mean of 

observations

(°C)

Mean error

(°C)

Root mean 

square error 

(°C)

VW UWMDL

UWMDN VW UWMDLUWMDN

Sites with two water-quality monitors

1 MDT-L 21.53 -0.26 0.63 -0.21 0.49 0.75 0.47

2 MDT-U 22.86 .01 1.21 .20 .92 1.46 .90

3 EPT-L 21.42 -.50 .39 -.38 .71 .59 .60

4 EPT-U 22.02 -.35 .77 -.09 .64 .90 .49

5 MDN-L 21.94 -.23 1.07 .31 .48 1.14 .50

6 MDN-U 22.37 -.11 1.26 .46 .43 1.34 .63

Sites with one water-quality monitor at 1 meter from the bottom

7 BLB 22.57 -0.63 1.12 0.37 0.77 1.37 0.59

8 EHB 21.61 .21 .83 -.14 .67 .99 .68

9 FBS 22.78 -.42 1.63 .72 .77 1.78 .97

10 HDB 22.35 .73 1.41 .58 1.00 1.52 .86

11 MDL 21.26 -.87 .50 -.31 .99 .67 .62

12 MPT 21.50 -.94 1.05 -.03 1.12 1.27 .83

13 NBI 21.80 .10 .99 -.14 .87 1.17 .58

14 RPT 21.67 .03 1.18 .28 .57 1.27 .60

15 SHB 22.66 -.85 1.37 .50 1.08 1.52 .84

16 WMR 22.18 -.68 1.60 .78 1.04 1.78 1.07

17 Overall Mean -.12 1.08 .20 .72 1.22 .70
30 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

July 26. The model correctly reproduces the large diel 

variation in the near-surface layer while the near-bottom 

temperature remains more nearly constant during the day; 

this is true even at site MDN, the shallowest of the three 

sites, where two continuous monitors were deployed. The 

ME and RMSE goodness-of-fit statistics at the nearbottom monitors at the deep sites tend to be bigger than 

the goodness-of-fit statistics at the near-surface monitors 

at the same sites (table 6, rows 1–6). 

The model validation run used observed 2006 

meteorology between May 20 and October 15. The 

value of the calibration parameter, rSW, was the same as 

that established during calibration for the 2005 season. 

The observed and simulated temperature records at 

the four deep sites where there was a near-surface and 

near-bottom continuous temperature monitor in 2006 are 

shown in fig. 21. The observed and simulated temperature 

at the eight shallow sites where there was one continuous 

monitor, placed 1 m from the bottom, are shown in 

fig. 22. In addition, in 2006 five continuous monitors 

were placed in nearshore locations (within 100 m of the 

shoreline), and the observed and simulated temperatures 

at these sites are shown in fig. 23. Error statistics were 

calculated over the entire 149-day simulation, and over 

July 26 to August 31 in order to compare to the results 

from the 2005 calibration (table 7). The ME was bigger 

over the longer period at some sites and smaller at others, 

but the overall ME was smaller over the entire simulation. 

The overall RMSE grew slightly from 0.75 to 0.88 from 

the 36-day simulation to the 149-day simulation (table 7, 

row 22). As in 2005, the largest errors were not at the 

deepest sites, where a monitor was located 1 m from 

the bottom and 1 m from the surface (table 7, rows 1–8; 

fig. 22). The top-to-bottom range in temperature at the 

deeper sites was generally well-simulated over the entire 

149-day simulation (fig. 24). 

At sites where temperature was collected in both 

2005 and 2006, the RMSE was comparable over July 26 

to August 31 (rows 1–10, 13,14, and 16 in table 6, and 

rows 1–6, 9–12, 14, 15, and 17 in table 7; fig. 17). The 

average of the RMSE over July 26 to August 31 at sites 

where temperature was collected in both 2005 and 2006 

was 0.72 in 2005 and 0.75 in 2006. The ME for both 

the month of August and the entire simulation indicates 

that the model underestimated the temperatures, as there 

are more and larger positive errors than negative. Closer 

consideration shows that underestimation of temperature 

was more common at shallow and nearshore sites. The 

ME was greater than 1 at several sites (table 7)—all of 

these except site HDB were nearshore sites that were 

added in 2006 (fig. 23). The measured temperatures at 

the nearshore sites were generally higher than measured 

temperatures at the other sites (table 7). This may 

indicate that sediments play a role at these very shallow sites 

by absorbing incoming shortwave radiation that is transmitted 

through the entire water column and transferring some of 

the heat back to the water column. In the model, shortwave 

radiation that is not absorbed within the water column does 

not contribute heat to the water column. At some nearshore 

sites, most notably GBE and WDW, the simulated temperature 

shows greater diel swings than the observations (fig. 23); 

this may be an indication that the depth at these sites was not 

accurately portrayed in the model grid. At these very shallow 

nearshore sites, small errors in the bathymetry make a large 

proportional error in the depth of the water column. 

Table 7. Goodness-of-fit temperature statistics for the UnTRIM 

model of Upper Klamath Lake, Oregon, 2006. 

[Errors are calculated hourly as observed value-simulated value. Site modifier: 

-L indicates monitor 1 meter from the bottom; -U indicates monitor 1 meter from 

the surface. A spatially variable wind was used in the simulation. Site names in 

boldface type indicate that temperature was measured at the site during both 2005 

and 2006. Full Season is from May 20 to September 25. Mid-summer is from July 

26 to August 31. Overall mean is calculated based on sites in boldface type only. 

Abbreviations: °C, degrees Celsius]

Row Site

Mean of 

observations

(°C)

Mean error

(°C)

Root mean 

square error 

(°C)

Full 

season

MidsummerFull 

season

MidsummerFull 

season

Midsummer

Sites with two water-quality monitors

1 MDT-L 18.00 20.34 -0.16 0.14 0.59 0.41

2 MDT-U 19.15 21.53 -.11 .26 1.02 .90

3 EPT-L 17.81 20.18 -.29 -.06 .71 .45

4 EPT-U 18.45 21.05 -.08 .23 .67 .54

5 MDN-L 18.19 20.57 -.21 .04 .71 .53

6 MDN-U 18.54 20.95 -.14 .17 .67 .49

7 SET-L 18.16 20.35 .23 .52 .71 .68

8 SET-U 18.91 21.04 .37 .61 .78 .80

Sites with one water-quality monitor at 1 meter from the bottom

9 EBB 17.97 20.30 -0.65 -0.56 1.10 0.92

10 EHB 18.31 20.52 .13 .50 .85 .83

11 FBS 18.42 20.94 -.79 -.41 1.39 1.02

12 HDB 18.94 21.28 .61 1.16 1.21 1.35

13 MRM 18.06 20.49 -.59 -.18 1.01 .73

14 RPT 18.14 20.30 .00 .25 .76 .67

15 NBI 18.36 20.49 .13 .43 .91 .82

16 WMR 18.59 21.04 -.05 .08 1.07 .93

Nearshore sites with one water-quality monitor at mid water column

17 GBE 19.78 21.64 1.12 0.69 2.22 1.47

18 HPK 19.42 21.15 1.36 1.50 1.72 1.86

19 SHL 20.00 21.60 .86 1.09 1.23 1.31

20 WDW 19.51 21.51 .57 .13 1.66 1.21

21 SSR 20.49 22.30 .89 1.25 1.42 1.44

22 Overall Mean -0.08 0.23 0.88 0.75
Simulation Results 31

10

5

20

30

10

5

20

30

10

5

20

30

10

5

20

30

10

5

20

30

10

5

20

30

UNIFORM WIND FORCING VARIABLE WIND FORCING

WATER TEMPERATURE, IN DEGREES CELSIUS 

1 15 1 15 1 15 1 15 1 14

JUNE JULY AUGUST SEPTEMBER OCTOBER 

2005

MDT-observed

EPT-observed

MDN-observed

MDT-simulated

MDN-simulated

EPT-simulated

Observations

EXPLANATION

Simulations

Upper

Lower

Upper

Lower

Figure 19. Observed and simulated temperatures at three deep sites in Upper Klamath Lake, Oregon, 2005.

OR19_0141_UKL_fig_18

15 1 15 1 15 1 15 1 15 1 14

MAY JUNE JULY AUGUST SEPTEMBER OCTOBER 

2005

0

10

20

30

0

10

20

30

AIR TEMPERATURE, IN DEGREES CELSIUS

EXPLANATION

WMR-MET

MDL

2006

Figure 18. Air temperature at site MDL in 2005 and 2006, and air temperature at site WMR-MET in 2005, Upper Klamath 

Lake, Oregon.
32 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

OR19_0141_UKL_fig_20

0

4

2

8

6

0

4

2

8

6

0

4

2

8

6

0

4

2

8

6

0

4

2

8

6

0

4

2

8

6

WATER TEMPERATURE DIFFERENCE, IN DEGREES CELSIUS 

EXPLANATION

15 1 15 1 15 1 15 1 15 1 14

MAY JUNE JULY AUGUST SEPTEMBER OCTOBER 

2006

EPT-observed

EPT-simulated

MDT-observed

MDT-simulated

MDN-simulated

MDN-observed

Observed daily maximum vertical temperature difference

Observed daily minimum vertical temperature difference

Simulated daily maximum vertical temperature difference

Simulated daily minimum vertical temperature difference

Figure 20. Daily maximum and minimum difference in temperature between near-surface and near-bottom of the 

water column at three deep sites in Upper Klamath Lake, Oregon, 2005.
Simulation Results 33

Observations

EXPLANATION

Simulations

Near-surface

Near-bottom

Near-surface

Near-bottom

10

20

30

10

20

30

10

20

30

10

20

30

10

20

30

10

20

30

10

20

30

10

20

30

15 1 15 1 15 1 15 1 15 1 14

MAY JUNE JULY AUGUST SEPTEMBER OCTOBER 

MDT-observed

WATER TEMPERATURE, IN DEGREES CELSIUS 

2006

MDT-simulated

MDN-observed

MDN-simulated

EPT-observed

EPT-simulated

SET-observed

SET-simulated

Figure 21. Observed and simulated temperatures at four deep sites in Upper Klamath Lake, Oregon, 2006.
34 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

OR19_0141_UKL_fig_22

10

20

30

20

30

10

10

20

30

10

20

30

10

20

30

10

20

30

15 1 15 1 15 1 15 1 15 1 14

MAY JUNE JULY AUGUST SEPTEMBER OCTOBER 

WATER TEMPERATURE, IN DEGREES CELSIUS 

2006

Observations

Simulations

10

10

20

30

20

30

RPT

WMR

FBS

HDB

NBI

EBB

MRM

EHB

Figure 22. Observed and simulated temperatures at eight shallow sites in Upper Klamath Lake, Oregon, 2006.
Simulation Results 35

OR19_0141_UKL_fig_23

10

20

30

10

20

30

10

20

30

10

20

30

10

20

30

WATER TEMPERATURE, IN DEGREES CELSIUS 

Observations

Simulations

15 1 15 1 15 1 15 1 14

JUNE JULY AUGUST SEPTEMBER OCTOBER 

2006

GBE

HPK

SHL

SSR

WDW

5

15

25

5

15

25

5

15

25

5

15

25

5

15

25

Figure 23. Observed and simulated temperatures at five nearshore sites in Upper Klamath Lake, Oregon, 2006.

Overall, the calibration of the model with 2005 

temperature data and the validation with 2006 temperature 

data was satisfactory and adequate for the purpose of 

simulating the weekly to seasonal variations in temperature. 

The temperature is used in density calculations within 

the model and therefore determines the vertical mixing 

characteristics of the water column. The simulation of 

temperature captures both the daily thermal stability and those 

periods of time when the water column stratifies in deeper 

parts of the lake for several days at a time. Even though the 

lake is shallow, there has been discussion about the role that 

temporary stratification plays in water quality (Kann and 

Welch, 2005; Wood and others, 2006). Buoyant cyanobacteria 

like AFA are particularly suited to taking advantage of mild 

thermal stability to position colonies within the photic zone. In 

addition, as discussed further below, the tendency for the water 

column to stratify in the trench may provide a mechanism, 

particularly during a rapid bloom decline, for concentrating 

rising AFA colonies in the central part of the lake. 
36 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

OR19_0141_UKL_fig_24

0

4

2

8

6

0

4

2

8

6

0

4

2

8

6

0

4

2

8

6

0

4

2

8

6

0

4

2

8

6

WATER TEMPERATURE DIFFERENCE, IN DEGREES CELSIUS 

15 1 15 1 15 1 15 1 15 1 14

MAY JUNE JULY AUGUST SEPTEMBER OCTOBER 

2006

EPT-observed

EPT-simulated

MDT-observed

MDT-simulated

MDN-simulated

MDN-observed

EXPLANATION

Observed daily maximum vertical temperature difference

Observed daily minimum vertical temperature difference

Simulated daily maximum vertical temperature difference

Simulated daily minimum vertical temperature difference

Figure 24. Daily maximum and minimum difference in temperature between near-surface and near-bottom of the water 

column at three deep sites in Upper Klamath Lake, Oregon, 2006.
Implications for Water Quality 37

Implications for Water Quality

Because severe water-quality problems in UKL are 

detrimental to the survival of endangered suckers, the 

hydrodynamic model is particularly useful insofar as it can 

provide insight into the processes that control water quality. 

The simulations of hydrodynamics and heat transport in 

2005 and 2006 result in the following observations that have 

important implications for water quality: first, that water 

moving northward through the trench to the west of Bare 

Island is routed both to the north and to the south of the lake, 

and second, that the water originating from deeper in the 

trench is routed preferentially toward the north, while the 

water from nearer the surface is routed preferentially toward 

the south.

These observations were illustrated and quantified 

by using model simulations of conservative tracers that 

were designed to demonstrate how water from the trench is 

transported to the rest of the lake. First, the polygons in the 

numerical grid east of Eagle Ridge and south of Bare Island 

that were deeper than 4.5 m were identified (fig. 25). Within 

those polygons, and within layers deeper than 4.5 m, the 

concentration of tracer T1 was set to 1 (arbitrary units) at 

every time step, and within those layers shallower than 4.5 m 

the concentration of tracer T1 was set to 0 at every time step. 

Within the same polygons, and within layers deeper than 

4.5 m, the concentration of tracer T2 was set to 0 at every time 

step, and within layers shallower than 4.5 m the concentration 

of tracer T2 was set to 1 at every time step. Within polygons 

outside of those identified, the initial concentrations of T1 

and T2 were set to 0. Thus, the concentration of T1 and T2 

at any point in the grid is a measure of the fraction of water 

at that point that originated in the trench at depths deeper 

than 4.5 m and shallower than 4.5 m, respectively. Assigning 

the concentrations of two distinct tracers in this way allowed 

tracking of the surface water of the trench separately from 

the deeper water. The tracer experiments began on August 1, 

2005, and used the calibrated hydrodynamic and heat-transport 

model described above. The concentrations in the trench 

were assigned continuously, rather than as an initial condition 

at a single point in time, in order to observe the continuous 

transport of water through the trench and out into the lake for 

several days; as a result the total mass of each tracer increases 

throughout the simulation. For that reason, the simulations are 

most interesting over the first few days and were stopped after 

10 days. 

“Snapshots” of the concentration of tracers T1 and T2 

in the surface layer of the grid at the end of day 5 are shown 

in figure 26. The snapshots show that water that exits the 

trench at the west of Bare Island (fig. 1) flows both to the east 

around the island and to the west into the northern part of the 

lake. The bifurcation of the flow at Eagle Point results in the 

development of two different circulation patterns, the longer 

Figure 25. The location of polygons in the numerical grid 

where the depth was greater than 4.5 meters and where tracer 

concentrations were defined for numerical experiments.

of which carries water around the northern end of the lake and 

down the eastern shoreline, then turns around north of Buck 

Island. This circulation pattern is responsible for the direct 

influence of the water exiting the trench on water quality in 

the northern part of the lake. A second, shorter circulation 

pattern carries water around the north and east sides of Bare 

Island and turns around north of Howard Bay. Previous 

numerical experiments with this model have shown that the 

relative strengths of these two circulation patterns varies with 

the strength of the wind forcing, such that a stronger prevailing 

wind is more effective at pushing water into the northern part 

of the lake and enhances the larger pattern at the expense of 

the smaller (Wood and Cheng, 2006).

OR19_0141_UKL_fig25

0 1 2 3 4 5 MILES

0 1 2 3 4 5 KILOMETERS
38 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

Consideration of the sum of the concentration of the 

two tracers together allows some quantification of the first 

observation provided above. The sum of concentrations of 

tracers T1 and T2 shows that water leaving the trench at the 

beginning of the simulation is detected at site MDN, centrally 

located in the northern part of the lake, in about 3 days, and 

the fraction of the water at site MDN that can be traced back to 

the trench increases rapidly over 7 days to a total of about 60 

percent (fig. 27A). The concentrations were simulated in the 

middle of the water column, but there is little top-to-bottom 

difference in the simulated tracer concentration at either site 

MDN or MDL. At site MDL, located in the central part of 

the lake, the influence of the trench can be detected sooner, 

within about 1 day, and within 3 days the trench accounts for 

about 80 percent of the water at site MDL. Thus, it appears 

Figure 26. Concentration of tracers (A) T1 and (B) T2 in the surface layer of Upper Klamath Lake, Oregon, at day 5 of numerical 

experiments starting on August 1, 2005.

that the influence of the trench at site MDL through the shorter 

circulation pattern is both quicker and quantitatively larger 

than at site MDN. Note that these fractions of trench water 

and the travel time from the trench apply specifically to the 

wind conditions that were measured during the period of the 

simulation, from August 1 to 10, 2005, when the winds were 

weak to moderate and with no major reversals in direction; 

relative fractions of trench water and travel times would be 

expected to change with different wind conditions (Wood and 

Cheng, 2006). Although the relative strength of the shorter 

and longer circulation loops varies with wind conditions, the 

observation that water moving northward through the trench 

to the west of Bare Island is routed both to the north and to 

the south of the lake is robust under conditions of prevailing 

winds.

OR19_0141_UKL_fig26b

B.

fig26a

A.

1

0.8

0.6

0.4

0.2

0

Concentration of 

tracer T1 in the 

surface layer at day 5 

of simulation.

0 1 2 3 4 5 MILES

0 1 2 3 4 5 KILOMETERS

EXPLANATION
Implications for Water Quality 39

OR19_0141_UKL_fig_27

MDN

MDL

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0 2 4 6 8 10

0 2 4 6 8 10

DAY OF SIMULATION

0 2 4 6 8 10

TRACER CONCENTRATION, NO UNITS

A. T1+T2

B. T1

C. T2

Figure 27. Time series of the simulated concentration of tracers 

T1 and T2 during 10 days of numerical experiments starting 

on August 1, 2005, at sites MDN and MDL in Upper Klamath 

Lake, Oregon. The total fraction of trench water at each site is 

represented by the sum of (A) T1 and T2, the fraction of water 

originating from depths greater than 4.5 m in the trench is 

represented by (B) T1, and the fraction of water originating from 

depths less than 4.5 m in the trench is represented by (C) T2.

The second observation, that the water originating from 

deeper in the trench is routed preferentially toward the north, 

while the water from nearer the surface is routed preferentially 

toward the south, is quantified by considering the difference 

in the time series of tracers T1 and T2 at the two sites. When 

the flow through the trench west of Bare Island splits into two 

different pathways, more of the surface water flows to the east 

around the north side of Bare Island, and more of the bottom 

water flows to the west into the northern part of the lake. 

Thus, while the overall influence of the trench at site MDN is 

less than that at site MDL, the influence of water from below 

4.5 m is less at site MDL (fluctuating between about 0 and 

10 percent during days 5–10 of the simulation) than at site 

MDN (fluctuating between about 10 and 20 percent during 

days 5–10 of the simulation; fig. 27B). Similarly, the influence 

of water from above 4.5 m is less at site MDN (fluctuating 

between about 30 and 40 percent during days 5–10 of the 

simulation) than at site MDL (fluctuating between about 90 

and 100 percent during days 5–10 of the simulation; fig. 27C). 

The trench plays a significant role in the overall oxygen 

budget of UKL, a conclusion supported by at least two 

separate datasets. First, light and dark bottle experiments 

have shown that the trench (site MDT) is characterized by net 

oxygen consumption, whereas much of the rest of the lake is 

characterized by net oxygen production (fig. 28, Appendix A). 

The data shown in figure 28 were collected between May 

and September of 2006; the methods used are discussed in 

Hoilman and others (2008). Second, oxygen isotopes of the 

dissolved oxygen in the lake are consistent with respiratory 

processes dominating in the trench and with domination by 

photosynthetic production away from the trench (S.R. Silva, 

U.S. Geological Survey, unpub. data, 2007). In addition, 

dissolved oxygen measurements from continuous monitors 

show that the lowest dissolved oxygen events in the northern 

part of the lake are not restricted to that part of the lake, but 

consistently are concurrent with very low dissolved oxygen 

concentrations in the trench as well (Wood and others 2006, 

Hoilman and others, 2008). Of particular interest is the 

influence the trench has on dissolved oxygen concentrations in 

the prime adult sucker habitat in the northern part of the lake, 

to the north and west of Eagle Point (Wood and others, 2006).

Observations from continuous monitors in 2005 show that 

the low dissolved oxygen concentrations that accompanied 

the annual bloom decline at the end of July 2005 were more 

severe in the northern part of the lake, as measured at site 

MDN, and were closer to the concentrations (particularly 

the near-bottom concentrations) detected in the trench at site 

MDT, than the concentrations in the southern part of the lake, 

as measured at site MDL (fig. 29). This is consistent with 

the numerical tracer experiments that demonstrate a greater 

influence of the near-bottom water in the trench at site MDN 

than at site MDL and less influence of the surface water. 
40 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

OR19_0141_UKL_fig_28

-20

-10

0

10

0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500

CHANGE IN DISSOLVED OXYGEN CONCENTRATION,

IN MILLIGRAMS PER LITER

CHLOROPHYLL a, IN MICROGRAMS PER LITER 

A. MDN B. MDT C. RPT

Figure 28. The 24-hour change in dissolved oxygen concentration that would result from the oxygen production and 

consumption rates measured in light and dark bottle incubations at sites (A) MDN, (B) MDT, and (C) RPT, in Upper Klamath 

Lake, Oregon, between June and October, 2006, as a function of the chlorophyll a concentration measured at the same 

site on the same date. The 24-hour change in dissolved oxygen concentration was calculated at each site as described in 

Appendix A.

1 15 1 15 1 15 1 15 1 14

JUNE JULY AUGUST SEPTEMBER OCTOBER 

2005

MDN

MDL

MDT

0

5

10

15

20

25

0

5

10

15

20

25

DISSOLVED OXYGEN CONCENTRATION, IN MILLIGRAMS PER LITER

0

5

10

15

20

25

1 meter from bottom

1 meter from surface

Figure 29. Time series of dissolved oxygen concentration data collected at sites MDN, MDT, and MDL in Upper Klamath 

Lake, Oregon, 2005.
Implications for Water Quality 41

The implications are twofold. First, the near-bottom water 

has a lower dissolved oxygen concentration than the 

near-surface water. Dissolved oxygen is not a conservative 

quantity, however, and given that water at site MDT has 

several days of travel time to reach either site MDN or site 

MDL, the differing influence of water from the bottom of 

the trench on these sites does not fully explain why such 

different concentrations were measured by the continuous 

monitors at these sites. Another possibility is that as water 

passes through the trench, the surface and bottom of the 

water column develop different capacities for replacement 

of the oxygen consumed in the trench once the water exits, 

as discussed below. 

Vertical velocities measured by the current profiler 

at site ADCP1 in 2005 showed that at this site, but not at 

others, there was a clear distinction between the upper 

water column, where the AFA population was dominated 

by rising colonies, and the lower water column, where 

sinking colonies dominated the population (fig. 30; Gartner 

and others, 2007). The velocities measured by the ADCP 

are not water velocities, but rather the velocities of particles 

suspended in the water column. Vertical water velocities 

are expected to be one to two orders of magnitude lower 

than the vertical velocities measured by the ADCP. This 

suggests that under certain circumstances the effect of the 

two different circulation patterns is to concentrate rising 

colonies in the area between Bare Island and Howard 

Bay, and toward the western shoreline, and to send water 

relatively depleted of colonies (because they settled out 

in the trench), or dominated by sinking colonies, into the 

northern part of the lake. This may imply a greater capacity 

for photosynthetic production, on a per volume basis, in the 

central part of the lake between Bare Island and Howard 

Bay than in the northern part of the lake. Photosynthetic 

production of oxygen is a dominant term in the oxygen 

budget and dominates high respiratory and other 

consumptive demands over much of the lake. In the trench 

the opposite is true—consumption dominates production. If 

the water exiting the trench to the west around Eagle Point 

loses much of its capacity for photosynthetic production 

for at least a period of several days, while water column 

oxygen demands remain high, then the dissolved oxygen 

in the northern part of the lake, in the upper lobe of the 

larger clockwise circulation pattern, will be replenished 

0 2 4 6 8 10 12 14 16 18 20 22 24

HOUR OF DAY

-0.1

0.0

0.1

0.3

0.2

AVERAGE VERTICAL VELOCITY, IN CENTIMETERS PER SECOND

-0.5

-0.4

-0.3

-0.2

-0.1

A. NEAR SURFACE

B. NEAR BOTTOM

Figure 30. Average hourly vertical velocities at site 

ADCP1 (A) near surface and (B) near bottom, as measured 

by an Acoustic Doppler Current Profiler (ADCP), Upper 

Klamath Lake, Oregon. Velocities were averaged between 

June 21 and September 12, 2005.

more slowly than in the central part of the lake. This would be 

consistent with the bottle incubation experiments that showed 

that site RPT (like site MDL, located in the southern part of 

the lake) was consistently a site of net production of dissolved 

oxygen, whereas the results at site MDN were variable 

(fig. 28). Further development and testing of this hypothesis, 

both with numerical modeling and field work designed to 

better understand the physiology of the AFA blooms, is 

needed.
42 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

Acknowledgments

The authors thank Jerad Bales, Annett Sullivan, and 

Chauncey Anderson for their detailed and thoughtful reviews 

of the manuscript, which greatly improved the report. 

Discussions with Michael MacWilliams and Kutay Celebioglu 

also were very helpful.

Summary

Upper Klamath Lake (UKL) is a hypereutrophic lake that 

experiences annual blooms of the buoyant cyanobacterium 

Aphanizomenon flos aquae. The severe water-quality 

conditions associated with these annual blooms (high pH, 

supersaturated to undersaturated dissolved oxygen, high 

un-ionized ammonia) are detrimental to the survival of two 

species of endangered suckers in the lake. As part of an 

ongoing and multipronged effort to understand the ecosystem 

of the lake, the U.S. Geological Survey and the Bureau of 

Reclamation entered into a cooperative agreement to develop a 

three-dimensional numerical model of the hydrodynamics and 

heat transport in the lake. 

The model is built on the UnTRIM computational core. 

The important features of this model include a semi-implicit 

finite difference solution method for the governing equations, 

flexibility in the use of the hydrostatic assumption, the use of 

a mass-conserving scheme to solve the constituent transport 

equation that also ensures that the solution is bounded by the 

initial and boundary conditions, and the use of a flux limiter 

in the solution of the transport equation to preserve accuracy 

while relaxing the stability constraint on the size of the 

polygons in the numerical grid. UnTRIM solves the governing 

equations on an unstructured orthogonal grid. The advantage 

of this type of grid is that it allows the size of the polygons 

that make up the grid to vary over the domain according to the 

requirements of the bathymetry and geometry, and allows for a 

shoreline-fitting boundary. The solution schemes used to solve 

the governing equations within the model provide a robust, 

stable, and computationally efficient platform on which to 

develop the UKL model.

Source and sink terms to calculate surface heat fluxes 

were added to the UnTRIM computational core. These 

include incoming shortwave radiation, atmospheric longwave 

radiation, reflected longwave radiation, evaporative cooling, 

and conduction. Two additional enhancements were made to 

the computational core. The first was a submodel designed 

to estimate a spatially variable surface wind based on massconserving interpolation. The second was a two-equation 

turbulence closure model designed to calculate spatially 

and temporally varying vertical diffusivities as a function of 

turbulent kinetic energy and its rate of dissipation. 

The hydrodynamic model was calibrated for June through 

September 2005 by using two parameters to adjust surface and 

bottom shear stress. Two different spatially uniform surface 

winds were used to run the model from either June or July 

through August—one measured in the northern part of the lake 

and one measured in the central part of the lake. The model 

also was run using a spatially variable surface wind from 

August through September. The spatially variable surface wind 

was generated from the variable-wind submodel using wind 

data collected at two sites on the lake and four around the 

shoreline. Simulations resulting from all three surface winds 

could be compared during the midsummer overlap period. 

Datasets that were used to evaluate model simulations during 

the calibration process included surface elevation determined 

by 3 gages around the lake, currents measured by Acoustic 

Doppler Current Profilers (ADCPs) at 5 sites in the lake, and 

temperature measured at 14 sites by continuous monitors.

The model predicted wind-driven circulation patterns 

in the lake that reproduced the prevailing currents measured 

by ADCPs—a clockwise circulation consisting of broad, 

shallow flow with the wind (to the south-southeast) on the 

eastern side of the lake, and a narrow band of flow opposing 

the wind (to the north-northwest) through the deep trench 

along the western shoreline and passing to the west of Bare 

Island. The model simulations interpolate currents between the 

ADCP measurements sites; they reveal that under prevailing 

wind conditions there are two predominant modes to this 

circulation. The first is a smaller gyre that circulates water 

between Bare Island and Rattlesnake Point; the second is a 

larger gyre that circulates water between Rattlesnake Point and 

the northernmost part of UKL. 

In order to evaluate the performance of the model, the 

mean error (ME) and root mean squared error (RMSE) of 

the hourly measured and simulated currents and temperature 

were calculated. In the process of model calibration with 

data collected during June to September of 2005, emphasis 

was placed on correctly simulating the highest velocities 

through the trench. The ME and RMSE of the VW simulation 

(calculated over 37 days between July 26 and August 31) 

in the depth-averaged speed at site ADCP1 (located in the 

middle of the trench) were small (0.50 and 3.08 centimeters 

per second (cm/s), respectively) in comparison to the mean 

of the depth-averaged speed at that site over the same period 

(11.23 cm/s). The simulated depth-averaged speed at two 

sites that were located in the northern part of the lake where 

velocities were lower were biased low, as indicated by a 

larger, positive ME of 0.80 cm/s at site ADCP3 (mean depthaveraged speed 4.09 cm/s) and 2.58 cm/s at site ADCP5 

(mean depth-averaged speed 5.81 cm/s). Consideration of the 

velocity components near the surface and near the bottom 

showed that the goodness-of-fit errors were higher near the 

surface than near the bottom, indicating that the model may 

have difficulty in correctly simulating the surface boundary 

layer, particularly in those areas of the lake where the depthaveraged currents oppose the prevailing wind stress. The 

use of a spatially variable wind forcing that was interpolated 

between two meteorological sites on the lake and four on 

the shoreline generally resulted in improved error statistics. 

The improvement obtained by using the spatially variable 
Summary 43

surface wind was particularly evident upon closer inspection 

of the individual velocity components at a site in the northern 

part of the lake. At this site (ADCP3), the spatially variable 

surface wind correctly simulated the diel phase shift between 

the bottom and surface currents, whereas the uniform surface 

wind did not. In the process of comparing results among 

the spatially variable surface wind simulation and the two 

uniform surface wind simulations, it was demonstrated that, 

if a uniform wind must be used, the measurement should 

be made in the central part of the lake rather than in the 

northern part. In interpreting any of these results it is useful to 

remember that there is some inherent accuracy limitation to 

the ADCP measurements, and that accuracy may be a function 

of distance from the ADCP (in other words, near-bottom 

measurements may be somewhat more accurate than nearsurface measurements). 

The calibration of the heat transport was accomplished 

by adjusting one calibration parameter that specified the 

amount of incoming shortwave radiation reflected at the water 

surface. The error statistics indicated a small, high bias in the 

simulated temperatures over the lake, as the ME calculated 

between July 26 and August 31 at all the sites was more often 

negative than positive, although it was less than 1 degree 

Celsius (°C) at all sites, ranging from –0.94 to 0.73°C when a 

spatially variable surface wind was used. The RMSE was less 

than 1°C at all but one site and ranged from 0.40 to 1.12°C 

when a spatially variable surface wind was used. The model 

accurately reproduced thermal stability in the water column on 

a daily basis, and also correctly identified those periods when 

the stratification was maintained for several days.

The model validation simulations spanned 123 days 

between May 15 and October 15, 2006. Those simulations 

used a spatially variable surface wind and the same calibration 

parameters as established for 2005. Fewer ADCP velocity 

records were available to verify the model in 2006 (two sites), 

but one of these was site ADCP1, a site located in the deep 

trench where an ADCP also was deployed in 2005. This 

provided the opportunity to directly compare the performance 

of the model at this site for a comparable period, roughly 

the month of August, in both years. The ME and RMSE of 

the depth-averaged speed at this site, where the mean depthaveraged speed was 12.98 cm/s in 2006, were 2.30 and 

3.88 cm/s, compared to 0.50 and 3.08 in 2005, so the model 

did not perform as well in 2006. The larger positive ME 

indicates that the velocities in the trench were underestimated; 

model calibration may not be optimized for multiple years. 

A second ADCP site in 2006 was located in an area of very 

low velocities (mean depth-averaged speed 3.42 cm/s) in the 

northern part of the lake. This site was farther from the trench 

and the main circulation gyre than the two sites where ADCPs 

were placed in the northern part of the lake in 2005; the 

accuracy of water currents simulated at this site was poor. The 

error statistics showed that most of the error was in the surface 

currents, as might be expected because the measured velocities 

at this site largely opposed the prevailing wind stress. 

The results of the hydrodynamic and heat transport model 

for water quality in the lake indicate first, that water moving 

northward through the trench to the west of Bare Island is 

routed both to the north and to the south of the lake, and 

second, that the water originating from deeper in the trench 

is routed preferentially toward the north, while the water 

from nearer the surface is routed preferentially toward the 

south. Numerical tracer experiments were used to illustrate 

and quantify these observations. These tracer experiments 

are specific to the wind conditions, which were weak to 

moderate, during the dates of the simulation—from August 1 

to 10, 2005. The experiments showed that 10 percent or less 

of the water passing through site MDL, in the southern part 

of the lake, originated from depths greater than 4.5 meters 

in the trench within the previous 5 days, whereas as much 

as 20 percent of the water at site MDN, in the northern part 

of the lake, originated from there. Similarly, the amount of 

water passing through site MDL that originated from depths 

shallower than 4.5 meters in the trench was between 90 and 

100 percent, whereas only about 40 percent of the water 

passing through site MDN originated from there. These 

percentages were obtained for the specific conditions of early 

August 2005 and therefore may not be appropriate for other 

conditions. However, the more qualitative observation that 

surface water is circulated very effectively from the trench 

through the central part of the lake between Bare Island and 

Rattlesnake Point is robust and generally applicable, as is the 

observation that much less of the surface water but most of the 

deeper water in the trench moves instead toward the northern 

part of the lake. This means that during severe low dissolved 

oxygen events, most of the deeper water with a very low 

dissolved oxygen concentration moves into the prime adult 

sucker habitat in the northern part of the lake. 

Given that dissolved oxygen is not conservative and that 

travel times are several days, this does not entirely explain 

why dissolved oxygen concentrations tend to be much lower 

and closer to concentrations in the trench during these events 

in the northern part of the lake than in the central part of the 

lake. Some difference in the capacity to replenish dissolved 

oxygen concentration between surface and deep water also is 

implied. Vertical velocities collected at site ADCP1 suggest 

that the trench is deep enough to effectively separate rising 

and sinking Aphanizomenon flos aquae (AFA) colonies. In 

conjunction with the two circulation gyres, this may provide 

a mechanism to concentrate healthier, rising AFA colonies 

in the central part of the lake, while at the same time causing 

the northern part of the lake to become relatively depleted of 

AFA colonies capable of photosynthetic oxygen production. 

This hypothesis could help explain why the water-quality 

conditions associated with rapid AFA bloom declines are more 

severe in the northern part of the lake, but further testing of 

this idea is required. 
44 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

References Cited

Adams, B., Swigle, B., and Shively, R., 2003, Movements and 

behavior of radio-tagged adult Lost River and shortnose 

suckers with respect to water quality in Upper Klamath 

Lake: U.S. Geological Survey Draft Annual Report of 

research to the U.S. Bureau of Reclamation, Mid-Pacific 

Region, Klamath Area Office, Klamath Falls, Oregon.

Bowen, I.S., 1926, The ratio of heat losses by conduction and 

by evaporation from any water surface: Physics Review, 

v. 27, p. 779–787.

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 

(Second edition): U.S. Environmental Protection Agency, 

EPA/600/3–85/040, 455 p.

Bradbury, J.P., Colman, S.M., and Dean, W.E., 2004, 

Limnological and climatic environments at Upper Klamath 

Lake, Oregon during the past 45,000 years: Journal of 

Paleolimnology, v. 31, no. 2, p. 167–188.

Bradbury, J.P., Colman, S.M., and Reynolds, R.L., 2004, The 

history of recent limnologic changes and human impact on 

Upper Klamath Lake, Oregon: Journal of Paleolimnology, 

v. 31, p. 151–165.

Bureau of Reclamation, 2008, Pacific Northwest Cooperative 

Agricultural Weather Network, accessed April 17, 2008, at 

http://www.usbr.gov/pn/agrimet/wxdata.html

Casulli, Vincenzo, 1990, Semi-implicit finite-difference 

methods for the two-dimensional shallow water equations: 

Journal of Computational Physics, v. 86, p. 56–74.

Casulli, Vincenzo, 1999a, A semi-implicit finite difference 

method for non-hydrostatic, free-surface flows: 

International Journal for Numerical Methods in Fluids, 

v. 30, p. 425–440.

Casulli, Vincenzo, 1999b, A semi-implicit numerical method 

for non-hydrostatic free surface flows on unstructured grid, 

in Proceedings of the International Workshop on Numerical 

Modelling of Hydrodynamic Systems, Zaragoza, Spain, 

1999: Zaragoza, Spain, University of Zaragoza, p. 175–193.

Casulli, Vincenzo, and Cattani, E., 1994, Stability, accuracy 

and efficiency of a semi-implicit method for threedimensional shallow water flow: Computers & Mathematics 

with Applications, v. 27, p. 99–112. 

Casulli, Vincenzo, and Cheng, R.T., 1992, Semi-implicit 

finite difference methods for three-dimensional shallow 

water flow: International Journal for Numerical Methods in 

Fluids, v. 15, p. 629–648.

Casulli, Vincenzo, and Walters, R.A., 2000, An unstructured 

grid, three-dimensional model based on the shallow water 

equations: International Journal for Numerical Methods in 

Fluids, v. 32, p. 331–348.

Casulli, Vincenzo, and Zanolli, P., 1998, A three-dimensional 

semi-implicit algorithm for environmental flows on 

unstructured grids, in Baines, M.J., Proceedings of the 

Conference On Numerical Methods for Fluid Dynamics, 

Oxford, England, 1998: Oxford, England, University of 

Oxford, p. 57-70.

Casulli, Vincenzo, and Zanolli, P., 2002, Semi-implicit 

numerical modeling of nonhydrostatic free-surface flows 

for environmental problems: Mathematical and Computer 

Modelling, v. 36, p. 1131–1149.

Casulli, Vincenzo, and Zanolli, P., 2005, High resolution 

methods for multidimensional advection–diffusion 

problems in free-surface hydrodynamics: Ocean Modelling, 

v. 10, p. 137–151.

Celebioglu, T.K., and Piasecki, M., 2006a, Comparison of 

turbulence models with use of UnTRIM in the Delaware 

Bay, in Proceedings of the Ninth International Conference 

on Estuarine and Coastal Modeling, Charleston, S.C., 2005: 

American Society of Civil Engineers, p. 105–122.

Celebioglu, T.K., and Piasecki, M., 2006b, Simulation of 

salinity and suspended sediment dynamics in Delaware Bay 

by comparing turbulence closure models, in Proceedings 

of the 7th International Conference on Hydroscience and 

Engineering (ICHE–2006), Philadelphia, Pennsylvania, 

2006: Philadelphia, Drexel University, accessed April 30, 

2008, at http://idea.library.drexel.edu/handle/1860/1482

Cheng, R.T., Gartner, J.W., and Wood, T.M., 2005, Modeling 

and validation of wind-driven circulation in Upper Klamath 

Lake, Oregon, in Walton, Raymond, ed., Proceedings of the 

2005 World Water and Environmental Resources Congress: 

Anchorage, Alaska, Environmental and Water Resources 

Institute of American Society of Civil Engineers, accessed 

April 30, 2008, at http://scitation.aip.org/getabs/servlet/Geta

bsServlet?prog=normal&id=ASCECP00017304079200042

6000001&idtype=cvips&gifs=Yes

Colman, S.M., Bradbury, J.P., McGeehin, J.P., Holmes, 

C.W., Edginton, David, and Sarna-Wojcicki, A.M., 2004, 

Chronology of sediment deposition in Upper Klamath Lake, 

Oregon: Journal of Paleolimnology, v. 31, p. 139–149.

Colman, S.M., Bradbury, J.P., and Rosenbaum, J.G., 2004, 

Paleolimnology and paleoclimate studies in Upper Klamath 

Lake: Journal of Paleolimnology, v. 31, p. 129–138.

Dingman, S.L., 2002, Physical hydrology: Upper Saddle 

River, N.J., Prentice Hall, 646 p.
References Cited 45

Eilers, J.M., Kann, Jacob, Cornett, J., Moser, K., and St. 

Amand, A., 2004, Paleolimnological evidence of change 

in a shallow, hypereutrophic lake: Upper Klamath Lake, 

Oregon, USA: Hydrobiologia, v. 520, p. 7–18.

Eilers, J.M., Kann, Jacob, Cornett, J., Moser, K., St. Amand, 

A., and Gubala, C., 2001, Recent paleolimnology of Upper 

Klamath Lake, Oregon: Roseburg, Oregon, J.C. Headwaters, 

Inc., Report submitted to the Bureau of Reclamation, 44 p.

Fisher, L.H., and Wood, T.M., 2004, Effect of water-column 

pH on sediment-phosphorus release rates in Upper Klamath 

Lake, Oregon, 2001: U.S. Geological Survey WaterResources Investigations Report 03–4271, 25 p.

Garratt, J.R., 1977, Review of drag coefficients over oceans 

and continents: American Meteorological Society Monthly 

Weather Review, v. 105, p. 915–929.

Gartner, J.W., Wellman, R.E., Wood, T.M., and Cheng, R.T., 

2007, Water velocity and suspended solids measurements 

by in-situ instruments in Upper Klamath Lake, Oregon: U.S. 

Geological Survey Open-File Report 2007–1279, 136 p. 

plus data files.

Gill, A.E., 1982, Atmosphere–ocean dynamics, International 

Geophysics Series, v. 30: New York, Academic Press, 

662 p. 

Golub, G.H., and van Loan, C.F., 1996, Matrix computations 

(3rd ed.): Baltimore, Johns Hopkins University Press, 728 p. 

Herrett, T.A., Hess, G.W., Stewart, M.A., Ruppert, G.P., and 

Courts, M.L., 2006, Water resources data for Oregon, water 

year 2005: U.S. Geological Survey Water-Data Report 

OR–05–1, 942 p., accessed March 6, 2008, at 

http://pubs.usgs.gov/wdr/2005/wdr-or-05/

Hoilman, G.R., Lindenberg, M.K., and Wood, T.M., 2008, 

Water-quality conditions in Upper Klamath and Agency 

Lakes, Oregon, 2005: U.S. Geological Survey Scientific 

Investigations Report 2008–5026, 44 p., accessed April 30, 

2008, at http://pubs.usgs.gov/sir/2008/5026/

Hubbard, L.H., 1970, Water budget of Upper Klamath Lake 

southwestern Oregon: U.S. Geological Survey Hydrologic 

Investigations Atlas HA–351.

Janssen, K.D., 2006, Daily energy-budget and Penman 

evaporation from Upper Klamath Lake, Oregon: Portland, 

Oregon, Portland State University, M.S. thesis. 

Kann, Jacob, 1998, Ecology and water quality dynamics of a 

shallow hypereutrophic lake dominated by cyanobacteria 

(Aphanizomenon flos aquae): Chapel Hill, University of 

North Carolina, Ph.D. dissertation, 110 p.

Kann, Jacob, and Walker, W.W., 1999, Nutrient and 

hydrologic loading to Upper Klamath Lake, Oregon, 

1991–1998: Klamath Falls, Oregon, Report to Klamath 

Tribes Natural Resources Department [variously paged].

Kann, Jacob, and Welch, E.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.

Kantha, L.H., and Clayson, C.A., 1994, An improved mixed 

layer model for geophysical applications: Journal of 

Geophysical Research, v. 99, p. 25235–25266.

Lerman, Abraham, 1971, Time to chemical steady-states in 

lakes and oceans: Advances in Chemistry Series, v.106 , 

p. 30–76.

Lippert, C., and Sellerhoff, F., 2006, Efficient generation of 

orthogonal unstructured grids: Proceedings of the Seventh 

International Conference on Hydroscience and Engineering, 

Philadelphia, Pa., 2006: Philadelphia, Drexel University, 

accessed April 18, 2008, at http://idea.library.drexel.edu/

handle/1860/1472

Ludwig, F.L., Livingston, J.M., and Endlich, R.M., 1991, 

Use of mass conservation and critical dividing streamline 

concepts for efficient objective analysis of winds in 

complex terrain: Journal of Applied Meteorology, v. 30, 

p. 1490–1499.

Martin, J.L., and McCutcheon, S.C., 1999, Hydrodynamics 

and transport for water quality modeling: Boca Raton, 

Florida, CRC Press, 794 p.

McCutcheon, S.C., 1989, Water quality modeling, in French, 

R.H., ed., Volume I, Transport and surface exchange in 

rivers: Boca Raton, Florida, CRC Press, 334 p.

Miller, W.E., and Tash, J.C., 1967, Interim report Upper 

Klamath Lake studies: Corvallis, Oregon, Federal Water 

Pollution Control Administration, Eutrophication Research 

Branch, Pacific Northwest Water Laboratory [variously 

paged].

Oregon Climate Service, 2008, Zone 3—Climate data 

archives, accessed September 17, 2007, at http://www.ocs.

oregonstate.edu/page_links/climate_data_zones/monthly_

airport/nmez3.html

Perkins, D., Kann, Jacob, and Scoppettone, G.G., 2000, The 

role of poor water quality and fish kills in the decline of 

endangered Lost River and shortnose suckers in Upper 

Klamath Lake: U.S. Geological Survey, Biological 

Resources Division Report submitted to U.S. Bureau of 

Reclamation, Klamath Falls Project Office, Klamath Falls, 

Oregon, Contract 4–AA–29–12160, 39 p.
46 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

Phinney, H.K., and Peek, C.A., 1960, Klamath Lake, an 

instance of natural enrichment, in Transactions of the 

Seminar on Algae and Metropolitan Wastes: Cincinnati, 

Ohio, U.S. Public Health Service, 6 p.

Rodi, Wolfgang, 1993, Turbulence models and their 

application in hydraulics: Rotterdam, A.A. Balkema 

publishers, International Association for Hydraulic 

Research Monograph Series, 116 p.

Roe, P.L., 1986, Characteristic-based schemes for the Euler 

equations: Annual Review of Fluid Mechanics, v. 18, 

p. 337–365.

Scheffer, Marten, 1998, Ecology of shallow lakes: London, 

Chapman and Hall, 357 p.

Schlichting, Hermann, 1955, Boundary-layer theory: New 

York, McGraw-Hill, 535 p.

Smith, S.D., and Banke, E.G., 1975, Variation of the sea 

surface drag coefficient with wind speed: Quarterly Journal 

of the Royal Meteorological Society, v. 101, p. 665–673.

Swinbank, W.C., 1963, Longwave radiation from clear skies: 

Quarterly Journal of the Royal Meteorological Society, 

v. 89, p. 339–348.

Umlauf, Lars, and Burchard, H., 2003, A generic length-scale 

equation for geophysical turbulence models: Journal of 

Marine Research, v. 61, p. 235–265.

U.S. Fish and Wildlife Service, 2001, Biological/Conference 

opinion regarding the effects of operation of the Bureau 

of Reclamation’s Klamath Project on the endangered Lost 

River sucker, endangered shortnose sucker, threatened 

bald eagle, and proposed critical habitat for the Lost River/

shortnose suckers: Klamath Falls, Oregon, U.S. Fish and 

Wildlife Service, 188 p.

U.S Geological Survey, 2007, Water-resources data for the 

United States, water year 2006: U.S. Geological Survey 

Water-Data Report WDR–US–2006, accessed March 6, 

2008, at http://wdr.water.usgs.gov/wy2006/search.jsp

Walker, W.W., 2001, Development of a phosphorus TMDL 

for Upper Klamath Lake, Oregon: Concord, Massachusetts, 

Report prepared for Oregon Department of Environmental 

Quality, 80 p.

Warner, J.C., Sherwood, C.R., Arango, H.G., and Signell, 

R.P., 2005, Performance of four turbulence closure models 

implemented using a generic length scale method: Ocean 

Modeling, v. 8, p. 81–113.

Wood, T.M., and Cheng, R.T., 2006, Use of UnTRIM to 

investigate dissolved oxygen transport in Upper Klamath 

Lake, Oregon: Proceedings of the Seventh International 

Conference on Hydroscience and Engineering, Philadelphia, 

Pa., 2006: Philadelphia, Drexel University, accessed April 

18, 2008, at http://idea.library.drexel.edu/handle/1860/1424

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.

Wu, J., 1969, Wind stress and surface roughness at air–sea 

interface: Journal of Geophysical Research, v. 74, 

p. 444–455.

Wuest, Alfred, and Lorke, Andreas, 2003, Small-scale 

hydrodynamics in lakes: Annual Review of Fluid 

Mechanics, v. 35, p. 373–412.

Wunderlich, W.O., 1972, Heat and mass transfer between a 

water surface and the atmosphere, TN Report 14: Norris, 

Tennessee, Tennessee Valley Authority Water Resources 

Research Engineering Laboratory.
Appendix A 47

An estimate of the potential change in dissolved oxygen 

over a 24-hour period (ΔDO) was made from the measured 

oxygen production and consumption rates in the light and dark 

bottles at two levels in the water column and the measured 

extinction coefficient as follows. The gross production G 

in mg L-1 hr-1 was assumed to have the same exponential 

dependence with depth as light radiation. Then, in analogy 

with equation (21), the gross production at any depth in the 

water column Gz as a function of the gross production at the 

water surface G0 is given by

0 . e k z G G e z

− = (A1)

Solving for G0 in terms of the value obtained from 0.5 m depth 

yields

0.5

0 0.5 . e k G G e = (A2)

The gross production integrated over a water column of depth 

D is given by

( ) 0

0 3

0

1 

10 , e e

D

k z k D

e

G G G e dz e

k

− −

− = = − ∫ (A3)

where 10-3 is a proportionality constant. The units of G are 

mg m-2 hr-1.

Measured respiration rates (R) were very similar at both 

depths, and are therefore assumed to be constant throughout 

the water column for this calculation. Then the respiration 

integrated over the water column is given by

3 10 . RD R − = (A4)

where 10-3 is a proportionality constant and the units of R

are mg m-2 hr-1. The final simplification is to assume that the 

gross production rate is constant at G for 12 hours of daylight 

in every 24 hours and zero for the remaining 12 hours, and 

the respiration rate is constant for a full 24 hours. Then the 

estimated change in dissolved oxygen over a 24-hour period, 

expressed in terms of a concentration (mg L-1), is

 ( ) 3

0 10 12 (12 24 ) 1 24 . e k D

e

G DO G R e R

D k D

−

− ∆ = + = − + (A5)

Appendix A. Calculation of the Rate of Change in Dissolved Oxygen Over 

24 Hours from Light/Dark Bottle Incubations and the Extinction Coefficient
48 Modeling Hydronamics and Heat Transport in Upper Klamath Lake, Oregon, and Implications for Water Quality

This page intentionally left blank.
Manuscript approved for publication, May 5, 2008

Prepared by the USGS Publishing Network

Bill Gibbs

Ellen Hardy

Jackie Olson

Bobbie Jo Richey 

Linda Rogers

Sharon Wahlstrom

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://oregon.usgs.gov/
Wood and others —Modeling Hydrodynamics and Heat Transport in Upper Klamath Lake, Oregon—Scientific Investigations Report 2008–5076