Open access peer-reviewed chapter

Modeling Thermal Stratification Effects in Lakes and Reservoirs

Written By

Scott A. Wells

Submitted: November 27th, 2019 Reviewed: February 14th, 2020 Published: April 1st, 2020

DOI: 10.5772/intechopen.91754

From the Edited Volume

Inland Waters

Edited by Adam Devlin, Jiayi Pan and Mohammad Manjur Shah

Chapter metrics overview

763 Chapter Downloads

View Full Metrics


A brief overview of characteristics of stratified water bodies is followed by an in-depth analysis of the governing equations for modeling hydrodynamics and water quality. Equations are presented for continuity or the fluid mass balance; x-momentum, y-momentum, and z-momentum equations; mass constituent balance equation; the heat balance equation for temperature; and the equation of state (relating density to temperature and concentration of dissolved and suspended solids). Additional equations and simplifications such as the water surface equation and changes to the pressure gradient term are shown. Many of the assumptions that are made in water quality models are discussed and shown. Typical water quality source-sink terms for temperature, dissolved oxygen, algae, and nutrients are listed. A summary of some typical water quality models for lakes and reservoirs is shown. Two case studies showing how models can predict temperature and dissolved oxygen dynamics in stratified reservoirs are shown. The brief summary looks at ways to improve water quality and hydrodynamic models of lakes and reservoirs.


  • water quality modeling
  • hydrodynamic modeling
  • temperature modeling
  • reservoir modeling
  • dissolved oxygen modeling
  • reservoir
  • lake
  • stratification

1. Characteristics of lakes and reservoirs

Lakes and reservoirs are bodies of water that often serve multiple beneficial uses, such as water supply for municipal and agricultural use, recreation use, fishery enhancement, flood control, and power generation. Their physical, biological and chemical characteristics determine to a large extent how those beneficial uses are met. Survey texts, such as Wetzel [1] and Hutchinson [2], describe the important limnological processes that affect lake and reservoir water quality. An overview of reservoir dynamics and water quality is well-summarized in Martin et al. [3].

Lakes are different from man-made reservoirs where outlet (and perhaps inlet) hydraulic structures regulate the flow rates and often internal hydrodynamics of the reservoir. Not only does this flow regulation affect the reservoir temperature stratification, but also in consequence affects its water quality. An important distinction between rivers and lakes/reservoirs is the cycle of stratification that can occur throughout the year since most rivers are well-mixed vertically.

In some river systems though, stratification can occur if there are natural pools. For example, in the Chehalis River basin in Washington, USA, the Chehalis River is usually well-mixed except in pools of slow-moving water. This is shown where a large area of the Chehalis river has little to no channel slope and exhibits lake-like characteristics in Figure 1.

Figure 1.

Elevation drop along the Chehalis River, WA, USA, showing a section that is lake-like where summer stratification occurs. Sampling sites (multi-colored dots) are also shown.

Stratification in turn is related to the density of water as a function of temperature and dissolved substances. The progression of stratification during a summer period is shown in Figure 2 in a mountain lake during a summer period where the upper well-mixed layer, the epilimnion, is separated from the lower layer, the hypolimnion, by the strong density (temperature) gradient. Figure 3 shows the typical inverse stratification in the wintertime. Oftentimes, ice formation on the surface can impede gas transfer and create winter-time oxygen deficits even though there is reduced biological activity as a result of the cold temperatures.

Figure 2.

Progression of stratification in summer of Bull Run Lake, OR, USA, during 1997.

Figure 3.

Bull Run Lake, OR, USA, temperature profile on January 19, 1993.

The progression of summer stratification can also influence the progression of dissolved oxygen depletion (see Figure 4 for Tenkiller Reservoir, OK, USA). This seasonal depletion in Figure 4 includes both the metalimnetic minimum (caused by hydrodynamic interflow of low-dissolved oxygen water at the base of the epilimnion) and the hypolimnetic depletion as a result of sediment oxygen demand.

Figure 4.

Tenkiller reservoir dissolved oxygen profiles in 2006 showing progression of summer oxygen depletion.

Also, as a result of internal seiching, wind dynamics, surface cooling, and solar radiation input, the vertical profiles for water quality parameters can vary during the day. For example, Hemlock Lake temperature and dissolved oxygen vertical profiles are shown in Figures 5 and 6, respectively, for the morning (9 am) and early afternoon (1 pm). Variation of 1–2°C and 4–5 mg/l dissolved oxygen concentrations were noted over the 4-hour time difference between profiles.

Figure 5.

Hemlock Lake, NY, USA temperature profile July 13, 2013 at 9 am and 1 pm.

Figure 6.

Hemlock Lake, NY, USA dissolved oxygen profile July 13, 2013 at 9 am and 1 pm.

Showing the effect of diurnal wind on seiching dynamics, Figure 7 shows a temperature buoy at a depth of 15 m in Chester Morse Lake, WA, USA, where variations of 2–3°C can be common diurnally as wind-induced seiching occurs.

Figure 7.

Internal seiching as evident in temperature dynamics at a depth of 15 m in Chester Morse Lake, WA, USA. Variations of 2°C occur at a diurnal time scale are evident during the later spring and summer as a result of wind seiching and closeness to vertical temperature gradient.

In order to describe these changes in water quality in a lake or reservoir, the next section describes the mathematical framework for modeling lakes and reservoirs.


2. Governing equations for lake and reservoir water quality modeling

The basic governing equations for hydrodynamics and water quality were discussed by Wells et al. [4] and summarized and simplified here. The hydrodynamic governing equations include conservation of water mass and momentum. The water quality governing equations include conservation of constituent mass and heat including processes such as advection, turbulent diffusion, molecular diffusion (and dispersion if there is spatial averaging). An equation of state is used to relate the water density to salinity, temperature, and suspended solids that can affect fluid momentum.

2.1 Governing equations for mass, momentum, constituent mass and heat conservation

The equations for fluid motion are based on mass and momentum conservation. The development of the governing equations is based on a control volume of homogeneous properties. The conservation of fluid mass is the change in fluid mass within the control volume equaling the sum of mass inflows to the control volume and the sum of mass outflows from the control volume. The conservation of momentum is based on evaluating the sum of forces acting on a control volume in x, y, or z (for a Cartesian system) and equating these to the acceleration of a control volume as shown in Figure 8. Mathematically, conservation of momentum is described as F = m a , where F : vector forces acting on control volume, m: mass within control volume, a : acceleration of fluid within control volume.

Figure 8.

Example of a force acting on a control volume resulting in the acceleration of the fluid within the control volume.

The general coordinate system used in the development of the governing equations is shown in Figure 9. The rotation of the coordinate system can result in significant horizontal accelerations of fluids. This is usually restricted to large water bodies such as large lakes (such as the Great Lakes in the USA) and oceanic systems. The body force that causes horizontal accelerations because of the spinning coordinate system is termed the Coriolis force.

Figure 9.

Definition sketch of coordinate system for governing equations where x is oriented east, y is oriented north, and z is oriented upward opposite gravity, Ω is the angular velocity of the earth spinning on its axis and ϕ is the latitude.

The continuity (or conservation of fluid mass) and the conservation of momentum equations for a rotating coordinate system [5, 6, 7] are the governing equations used to determine the velocity field and water level.

The final form of the governing equations is obtained by making the following assumptions:

  • the fluid is incompressible, where Δρ ρ < < 1 where ρ is the fluid density and Δρ is the change in density,

  • the centripetal acceleration is a correction to gravitational acceleration,

  • the Boussinesq approximation (which is related to the incompressibility assumption) is applied to all terms in the momentum equation except those dealing with density gradient induced accelerations, i.e. 1 ρ = 1 ρ o + ρ 1 ρ o where ρ = ρ o + ρ , ρ o is a base value,

  • all velocities and pressure are turbulent time averages, i.e., u = u ¯ + u ´ , where u ¯ = 1 T t t + T udt and u ´ is the temporal fluctuation of u about the mean, and similarly for the velocity in the y-direction, v = v ¯ + v ´ , the velocity in the z direction w = w ¯ + w ´ , and the pressure, p = p ¯ + p ´

The governing equations become after time averaging and simplifying:

2.1.1 Continuity

u ¯ x + v ¯ y + w ¯ z = 0 E1

where u ¯ : temporal mean velocity in the x-direction, v ¯ : temporal mean velocity in the y-direction, w ¯ : temporal mean velocity in the z-direction. The continuity equation is usually also integrated vertically to provide the water surface equation, such that η ¯ t = x η h u ¯ dz + y η h v ¯ dz η h qdz where q is removal from or inflow to a model cell in units of flow rate per unit length, z = h is the location of the bottom referenced to a datum, and z=η is the water surface level referenced to a datum. This equation is used to solve for the water surface elevation.

2.1.2 X-momentum equation

u ¯ t unsteady acceleration + u ¯ u ¯ x + v ¯ u ¯ y + w ¯ u ¯ z convective acceleration 2 Ω z v ¯ Coriolis acceleration = 1 ρ p ¯ x pressure gradient + μ ρ 2 u ¯ x 2 + 2 u ¯ y 2 + 2 u ¯ z 2 viscous stresses + 1 ρ τ xx x + τ xy y + τ xz z turbulent stresses E2

where: τ xx = ρ u u ¯ where τxx is the turbulent shear stress acting in x direction on the x-face of control volume, τ xy = ρ u v ¯ where τxy is the turbulent shear stress acting in x direction on the y-face of control volume, τ xz = ρ u w ¯ where τxz is the turbulent shear stress acting in x direction on the z-face of control volume, μ = dynamic viscosity, Ω = component of Coriolis acceleration where: Ω z = Ω E sin φ , Ω y = Ω E cos φ , ϕ = latitude , Ω E = earth s rotation rate , and assuming 2 Ω y w ¯ is negligible. In general, the molecular viscous stresses are negligible except at boundaries. Analogous to laminar shear stress, the turbulent shear stresses are often parameterized as τ xx = μ turbulent xx u ¯ x = ρ u u ¯ , τ xx = μ turbulent xy u ¯ y = ρ u v ¯ , τ xz = μ turbulent xz u ¯ z = ρ u w ¯ where the term μ turbulent is the turbulent eddy viscosity analogous to molecular viscosity. The pressure is usually decomposed into the following terms: p ¯ = p ¯ a + g η z ρ ¯ dz where pa is the atmospheric pressure on the water surface and g is the acceleration due to gravity. The pressure gradient in the x-momentum then becomes after simplification 1 ρ p ¯ x = 1 ρ p ¯ a x + g η ¯ x g ρ η z ρ ¯ x dz .

2.1.3 Y-momentum equation

v ¯ t + u ¯ v ¯ x + v ¯ v ¯ y + w ¯ v ¯ z + 2 Ω z u ¯ = 1 ρ p ¯ y + μ ρ 2 v ¯ x 2 + 2 v ¯ y 2 + 2 v ¯ z 2 + 1 ρ τ yx x + τ yy y + τ yz z E3

where: τ yx = ρ v u ¯ where τyx is the turbulent shear stress acting in y direction on the x-face of control volume, τ yy = ρ v v ¯ where τyy is the turbulent shear stress acting in y direction on the y-face of control volume, τ yz = ρ v w ¯ where τyz is the turbulent shear stress acting in y direction on the z-face of control volume, and assuming 2 Ω x w ¯ is negligible. Analogous to laminar shear stress, the turbulent shear stresses are often parameterized as τ yx = μ turbulent yx v ¯ x = ρ v u ¯ , τ yy = μ turbulent yy v ¯ y = ρ v v ¯ , τ yz = μ turbulent xz v ¯ z = ρ v w ¯ . The pressure is usually decomposed into the following terms: p ¯ = p ¯ a + g η z ρ ¯ dz , and the pressure gradient in the y-momentum then becomes after simplification 1 ρ p ¯ y = 1 ρ p ¯ a y + g η ¯ y g ρ η z ρ ¯ y dz .

2.1.4 Z-momentum equation

w ¯ t + u ¯ w ¯ x + v ¯ w ¯ y + w ¯ w ¯ z = g 1 ρ p ¯ z + μ ρ 2 w ¯ x 2 + 2 w ¯ y 2 + 2 w ¯ z 2 + 1 ρ τ zx x + τ zy y + τ zz z E4

where: τ zx = ρ w u ¯ where τzx is the turbulent shear stress acting in z direction on the x-face of control volume, τ zy = ρ w v ¯ where τzy is the turbulent shear stress acting in z direction on the y-face of control volume, τ zz = ρ w w ¯ where τzz is the turbulent shear stress acting in z direction on the z-face of control volume, and neglecting the Coriolis terms 2 Ω y u ¯ + 2 Ω x v ¯ . Analogous to laminar shear stress, the turbulent shear stresses are often parameterized as τ zx = μ turbulent zx w ¯ x = ρ w u ¯ , τ zy = μ turbulent zy w ¯ y = ρ w v ¯ , τ zz = μ turbulent zz w ¯ z = ρ w w ¯ . In cases where vertical accelerations are much less than horizontal accelerations, this equation can be reduced to the hydrostatic equation, i.e., 1 ρ p ¯ z = g .

2.2 Conservation of constituent mass and heat: the ADVECTIVE diffusion equation

The conservation of constituent mass in a control volume is a sum of all the fluxes (advective and diffusive) into and out from the control volume plus sources and sinks (chemistry, biology, physics, withdrawals, inputs) within the control volume. Summing up the fluxes in each direction, assuming that the fluid is incompressible and that the molecular diffusivity, D, is homogeneous and isotropic, the advective diffusion equation becomes

c t unsteady change in concentration + u c x + v c y + w c z advective mass transpor t = D 2 c x 2 + 2 c y 2 + 2 c z 2 diffusive mass transpor t + S sources / sinks E5

where c is the concentration [M/L−3], S is the sources and sinks of reactions occurring in the control volume, or the reaction rate [ML−3 T−1].

This equation is a 3-D, unsteady equation that applies to all flow conditions: laminar and turbulent. Since we cannot determine the instantaneous velocity field, the x-y-and z momentum equations were time averaged and hence were only able to practically predict the temporal mean velocity. Similarly, we time average the conservation of mass/heat equation using time averages of the velocity field.

The instantaneous velocity and concentration are decomposed into a mean and an unsteady component. Similar to the velocity field shown earlier, for concentration, c, this becomes c = c ¯ + c where c ¯ = 1 T t t + T cdt and c is the fluctuation about the mean.

Substituting the time average and fluctuating components of concentration and velocities into the 3D governing equation and time averaging we obtain:

c ¯ t unsteady change in concentration + u ¯ c ¯ x + v ¯ c ¯ y + w ¯ c ¯ z mean advective mass transpor t = D 2 c ¯ x 2 + 2 c ¯ y 2 + 2 c ¯ z 2 molecular diffusive mass transpor t + x E x c ¯ x + y E y c ¯ y + z E z c ¯ z turbulent diffusive mass transpor t + S ¯ sources / sinks E6

where the turbulent mass fluxes in x, y and z were assumed to be defined as a gradient, diffusion-type process, such as u ´ c ´ ¯ = E x c ¯ x , v ´ c ´ ¯ = E y c ¯ y , w ´ c ´ ¯ = E z c ¯ z , Ex is the turbulent mass diffusivity in x [L2/T], Ey is the turbulent mass diffusivity in y[L2/T], Ez is the turbulent mass diffusivity in z [L2/T]. The new terms in the governing equation represent mass transport by turbulent eddies. As the intensity of turbulence increases, turbulent mass transport increases.

In turbulent fluids, Ex, Ey, and Ez >> D, and D can be neglected (except at boundaries or density interfaces where turbulent intensity may approach zero). The turbulent diffusion coefficients can be thought of as the product of the velocity scale of turbulence and the length scale of that turbulence. These coefficients are related to the turbulent eddy viscosity. In general, these turbulent diffusion coefficients are non-isotropic and non-homogeneous.

Spatial averaging of this equation leads to the introduction of “dispersion” coefficients which account for the transport of mass as a result of spatial irregularities in the velocity field.

These equations are also valid for heat transport and temperature modeling by substituting the concentration of heat, ρ c p T , where T is temperature, cp is the coefficient of specific heat at constant pressure and ρ is the density, such that the governing equation for temperature, T, becomes after simplification

T ¯ t unsteady change in temperature + u ¯ T ¯ x + v ¯ T ¯ y + w ¯ T ¯ z mean advective heat transpor t = D T 2 T ¯ x 2 + 2 T ¯ y 2 + 2 T ¯ z 2 molecular diffusive heat transpor t + x E x T ¯ x + y E y T ¯ y + z E z T ¯ z turbulent diffusive heat transpor t + S ¯ ρ c p heat flux E7

where DT is the molecular thermal conductivity for heat and Ex, Ey, and Ez are the heat and mass turbulent eddy diffusivities assuming they are of the same order of magnitude.

2.3 Equation of state

Since density is an important variable for the momentum equation to account for density-driven flows, the computation of density is accomplished through an equation of state where density is computed from dissolved and suspended solids concentrations ( c dissolved solids , c suspended solids ) and temperature, T, such as

ρ ¯ = f T ¯ c dissolved solids ¯ c suspended solids ¯ E8

Typical equations of state for fresh and saltwater have been published by Gill [8] and Ford and Johnson [9].

2.4 Solution of governing equations

There are six equations (continuity or conservation of fluid mass, conservation of momentum in x, y and z, and conservation of constituent mass or heat, equation of state) that we are solving for six unknowns: turbulent time average concentration (or temperature), velocities in x, y, and z, density and turbulent time average pressure (or water surface), i.e. c ¯ or T ¯ , u ¯ , v ¯ , w ¯ , ρ ¯ , and η ¯ or   p ¯ . The mathematical solution is dependent on specifying the following: (1) turbulent shear stresses or Reynolds stresses by specification of the turbulent eddy viscosities, (2) turbulent mass (heat) fluxes by specification of Ex, Ey and Ez, (3) initial and boundary conditions, (4) dynamic molecular viscosity and molecular diffusivity for computations at interfaces or boundaries (otherwise, they are usually neglected since all natural water bodies are highly turbulent), and (5) the Coriolis acceleration (if 2D horizontal or 3D for large water bodies).

Determination of the turbulent eddy viscosities and eddy diffusivities is often based on what are termed closure models that are based on the turbulent Schmidt number ( Sc = ratio of turbulent viscosity to turbulent diffusivity of mass ) and the turbulent Prandtl number ( Pr = ratio of turbulent viscosity to turbulent conductivity of heat ). Most experimental evidence suggests that the turbulent Sc and Pr numbers are close to unity for turbulent flows and that turbulent Sc or Pr numbers vary only little between flows. Even though many models use a constant value of these ratios such that mass and heat transfer turbulent coefficients are approximately equal, buoyancy affects that value [10, 11, 12].

Determination of turbulent eddy viscosities have been based on multiple approaches: (1) eddy viscosity models as a function of water stability [13, 14, 15, 16], (2) Mixing length models [17, 18], (3) One equation models for turbulent kinetic energy [19], (4) Two-equation k-ε models for turbulent kinetic energy and dissipation [11] and (5) Reynolds stress and algebraic stress models [11]. In many models, once the turbulent eddy viscosity is known, then the turbulent diffusion coefficients are computed from E μ turbulent ρ where the approximation is based on typical Sc or Pr numbers. Many water quality and temperature models for lakes and reservoirs use some form of a k-ε turbulence model [20].

Vertical boundary conditions for the hydrodynamic model usually involve a surface shear stress condition for the wind and a bottom shear stress condition for frictional resistance based on a specified friction coefficient (for example, Chezy or Manning’s). Vertical boundary conditions for temperature and water quality constituents are assumed to be known fluxes at the surface and bottom.

Horizontal boundary conditions for mass or heat include mass or heat fluxes as a result of advection and for hydrodynamics include water level (or head) or flow conditions. The flow conditions in outlets to stratified reservoirs can be complicated because of local vertical accelerations in the vicinity of the outlet. In many models, the vertical acceleration of a fluid parcel is assumed to be much less than the horizontal accelerations and hence the vertical momentum equation simplifies to the hydrostatic equation. In order to model the complicated outlet hydraulics in a reservoir, special selective withdrawal algorithms are often used [21, 22]. These allow the computation of flow from multiple vertical layers without having to solve the full-vertical momentum equation.

Typical assumptions of the flow field and water quality model are related to the dimensionality of the system (one, two or three-dimensions), whether the flow field is dynamic or steady-state, and the turbulence closure approximation. Based on the model assumptions, the model grid is developed where the governing equations are satisfied at points (differential equation representation) or over control volumes (integral representation). The resulting equations are then solved using numerical methods.

2.5 Sources-sinks for water quality and temperature

The source-sink term in the mass and heat conservation equation can be either positive or negative and is determined by each water quality state variable. The units of S ¯ in the mass conservation equation are [ML−3 T−1] with a typical unit of g/m3/s and in the heat balance equation the units are [Energy L−3 T−1] with a typical unit of J/m3/s. Table 1 shows some of the typical source sink terms for several water quality state variables. Details of these can be found in Wells [20] and Chapra [23].

State variable Typical source-sink term Description
Temperature S ¯ = φ z φ is the heat flux in units of W/m2 transmitted through the water body. This is the short-wave solar radiation transmitted through the water and is a function of light extinction. The variable z is assumed to be positive downward.
Salinity or conservative substance S ¯ = 0 No sources and sinks
Suspended solids S ¯ SS = w SS c SS z wss is the settling velocity of particles as a positive velocity, cSS is the concentration of suspended solids of a given size fraction. Often multiple size fractions are modeled independently using Stokes’ law for settling velocity, wss. The variable z is assumed to be positive downward.
CBOD S particulate ¯ = k CBODp c CBODp
w CBOD c CBODp z
S dissolved ¯ = k CBODd c CBODd
Source/sink terms are shown for dissolved CBOD (cCBODd) and particulate CBOD (cCBODp), kCBOD is a BOD decay rate for dissolved and particulate CBOD, and wCBOD is the settling velocity for particulate BOD. Models of CBOD usually use CBODultimate. Many models also track the P and N associated with this organic matter. Many models track multiple CBOD groups.
Algae S algae ¯ = μ growth c algae μ respiration c algae μ excretion c algae μ mortality c algae w algae c algae z Source sink terms include the algae growth rate μgrowth [T−1] (this is a complicated function of light, limiting nutrient and temperature), μrespiration [T−1] the “dark” respiration rate, μexcretion [T−1] the rate of excretion or biomass loss, μmortality [T−1] the mortality rate (which often can include zooplankton grazing as a separate loss rate based on zooplankton populations and zooplankton food preferences), and walgae the algae settling rate (this also can have complicated expressions especially for cyanobacteria and other species which migrate up and down in the water column). Often models include multiple algae groups. Calgae is the concentration of algae.
Ammonia-N S ammonia ¯ = δ aN μ growth c algae + μ respiration c algae + δ CBODdN k CBODd c CBODd + δ CBODpN k CBODp c CBODp + SOD N A V k nitr c ammonia The source/sink terms shown include algae uptake and release (where δaN is the stoichiometric equivalent of algae to ammonia-N, but the N source can be nitrate), organic matter release as particulate and dissolved CBOD decay (where δCBODdN is the stoichiometric equivalent of cBODd to N and δCBODpN is the stoichiometric equivalent of cCBODp to N), and sediment oxygen demand release under anoxic conditions (where SODN is the rate of N release in mass/area/time and V is the volume of the computational cell and A is the area of the sediment), nitrification decay rate knitr [T−1], and cammonia is the total ammonia concentration.
Dissolved oxygen S DO ¯ = δ aO 2 μ growth c algae μ respiration c algae k CBODd c CBODd k CBODp c CBODp SOD A V δ NO 2 k nitr c ammonia + k reaeraton c s c DO The source/sink term includes algae production and respiration (where δ aO 2 is the stoichiometric equivalent of dissolved oxygen to algae), CBOD particulate and dissolved water-column decay, sediment oxygen demand, nitrification demand (where δ NO 2 is the stoichiometric equivalent of dissolved oxygen to N), and reaeration at the surface only (where kreaeration is the reaeration rate in [T−1] which is generally a function of wind speed in lake and reservoirs and cs is the saturation value of dissolved oxygen). Other models also include terms for metal oxidation, methane oxidation, and oxidation of hydrogen sulfide.
Nitrate-Nitrite-N S NOx ¯ = k nitr c ammonia δ aNOx μ growth c algae k denit c NOx The source/sink terms include algae uptake (where δ aNO x is the stoichiometric equivalent of algae to nitrate-N since each algal group can have a preference for ammonia or nitrate as a N source), nitrification source, and a denitrification rate under anoxic conditions only (where kdenit is the denitrification rate under anoxic conditions). Other models also include terms for diffusion of nitrate into bottom muds. cNOX is the concentration of nitrite and nitrate.
PO4-P S PO 4 ¯ = δ aP μ growth c algae + μ respiration c algae + δ CBODdP k CBODd c CBODd + δ CBODpP k CBODp c CBODp + SOD P A V The source/sink terms shown include algae uptake and release (where δaP is the stoichiometric equivalent of algae to P), organic matter release as particulate and dissolved CBOD decay (where δCBODdP is the stoichiometric equivalent of cCBODd to P and δCBODpP is the stoichiometric equivalent of cCBODp to P), and sediment oxygen demand release under anoxic conditions (where SODP is the rate of P release in mass/area/time and V is the volume of the computational cell, calgae is the algae concentration, and A is the area of the sediment). Other models include adsorption of P onto inorganic particles.

Table 1.

Typical source-sink terms for temperature and some eutrophication water quality state variables.


3. Lake and reservoir water quality models

There are many models used to simulate reservoir and lake water quality. A summary of modeling approaches for lakes is shown in Mooij et al. [24] and Janssen et al. [25]. Table 2 shows a listing of some common lake and reservoir models.

Model name Description Reference
DYRESM and CAEDYM 1D model based on mixed layer dynamics, separate temperature and water quality models Tanentzap et al. [26]
CE-QUAL-W2 2D longitudinal-vertical, open source, eutrophication model, hydrodynamics and water quality solved together Wells [20]
CE-QUAL-R1 1D vertical Environmental Laboratory [27]
W3 3D, hydrodynamics and water quality solved together Al-Zubaidi and Wells [28]
EFDC and WASP 3D, hydrodynamics and water quality solved separately, both sigma stretch and z coordinate models Hamrick [29], Tetra Tech [30]
GLM 1D Hipsey et al. [31]
ELCOM and CAEDYM 3D-mixed layer dynamic model, hydrodynamics and water quality solved separately Hipsey et al. [32], Hodges and Dallimore [33]

Table 2.

List of common lake and reservoir water quality models.

The choice of a correct framework is dependent on several considerations: (1) dimensionality of the lake/reservoir system (even though all water bodies are in essence 3D, 2D and 1D models can often represent the important processes of water quality and temperature gradients), (2) documentation (up-to-date user manual with example problems), (3) ease of use and expertise required (all models require a degree of file manipulation and many include GUI interfaces that often facilitate running the model for new users), (4) established record of successful projects (as documented in papers and conference proceedings and technical reports) and (5) model processes represent important lake/reservoir processes (for example, if macrophyte growth is an important ecological consideration, does the model represent macrophytes).

In many cases, 3D models do not often do better than other model frameworks. One reason may be that the data and parameter uncertainty increase in higher dimensional models [34]. In a comparison of 2D and 3D models, many examples have shown [28, 35, 36] that 2D models often better represent temperature profiles than some 3D models. There may be many reasons for this, but the important message is that more complicated models do not necessarily mean better model predictions. Another issue with 3D models is the excessive computational time compared to lower dimensional models. In one comparison between a 2D and 3D model, the 3D model took 30× longer than the 2D model. This will vary depending on model configuration and model. This is becoming more of an issue as models are being used for multiple-decade simulations evaluating climate change and long-term changes in model boundary conditions.


4. Typical results of lake and reservoir modeling

Using the CE-QUAL-W2 model [20] as an example, consider an application to Folsom Reservoir, CA, USA, as presented in Martinez et al. [37].

Folsom lake, located near Sacramento California USA, is a deep-storage reservoir that provides municipal water, power generation and cold water for primarily salmonid fish in the lower American River (see Figure 10). The reservoir has multiple outlets that allow the operator to choose different water levels for downstream temperature control.

Figure 10.

Folsom reservoir bathymetry showing the north fork and south fork of the American River channels. Axes are labeled in m.

The model was set-up and calibrated to a 10-year period between January 1, 2001 and December 31, 2011. Boundary conditions for flow, meteorological data, and outflow during this period were developed. A very detailed approach for filling in data gaps was undertaken to provide a good set of boundary conditions. Typical model predictions compared to field data are shown for temperature in Figure 11 in 2002 and 2007 at multiple longitudinal stations in the reservoir. Error statistics for temperature profiles over the 10-year period using about 27,000 data comparisons were an average mean error of 0.004°C, an average absolute mean error (AME, average absolute value of the error) of 0.56°C, and a root mean square (RMS) average error of 0.71°C. The R2 correlation between modeled and predicted temperature was 0.996.

Figure 11.

Folsom reservoir model temperature predictions compared to field data on August 20, 2002 (left) and October 31, 2007 (right) at 6 different stations in Folsom reservoir.

In other examples of predicting the thermal regime, Cole [38] has shown that typical errors (AME, RMS) for temperature should often be well less than 1°C with a mean error of close to zero with minimal calibration if the boundary condition data are well-specified.

Oftentimes, the success of modeling other water quality state variables is first dependent on obtaining good temperature calibration results. For example, in a higher elevation pristine lake, Chester Morse Lake, WA, USA, Ceravich and Wells [39] have shown dissolved oxygen profiles mimicking the unusual behavior of the dissolved oxygen profile in a lake with little algae growth as shown in Figure 12. Error statistics for dissolved oxygen, which integrates all the water quality processes, were a ME of 0.15 mg/l, a AME of 0.42 mg/l, and a RMS error of 0.49 mg/l for 551 data-model comparisons.

Figure 12.

Predictions (solid lines) and field data (dots) of dissolved oxygen at one sampling site for Chester Morse Lake in 2015. Dates shown are Julian days since January 1, 2015.


5. Conclusions in hydrodynamic and water quality modeling

The complexity of existing models has often exceeded our capacity in the field to verify model coefficients usually because of cost and time. Deterministic water quality models require an incredible amount of information that is rarely measured. In the CE-QUAL-W2 model, for each algal group the model user must specify approximately 25 values describing rate coefficients for growth, respiration, excretion, mortality, stoichiometry, temperature preferences, N preferences, light saturation limits, and settling velocities. Even though this model has no limit to the number of algal groups one can represent mathematically, in a practical sense modeling living populations and their impact on nutrients, organic matter, pH, temperature, and oxygen is very complex. In the end, the model user tries to balance the known field data with literature values of the coefficients with the goal that if the boundary conditions are well-specified, the model requires little calibration.

If one cannot understand and interpret field data, then it will be challenging for a model to match field measurements. Hence, knowing and understanding the field data as one is setting up the model is important for making sure the model is agreeing with field data trends.

In other cases though, the model is able to discern complex interactions between water quality state variables that may be difficult for the model user to piece together a priori. For example, the unusual dissolved oxygen profiles in the field data and model shown in Figure 12 is one example where it was unclear the reasons for the unusual vertical profile until the combination of a sharp thermocline, algae growth within the metalimnion, and slow sediment oxygen demand caused the model to match the field data vertical trend.

Water quality models are adding more and more complex algorithms to reproduce admittedly complex phenomena. But this increasing complexity does not necessarily mean a better model or one that better reproduces field data. One example is the use of a complex model of bacterial populations on the Snake River in ID/OR, USA, from Harrison [40]. The bacterial populations were modeled based on Reichert et al. [41] as shown in Figure 13 and compared to a model with only a first order decay rate for organic matter decay (basically neglecting all the complex bacterial dynamics). In predicting the impact of organic matter on dissolved oxygen, the simpler model neglecting bacterial dynamics performed better. This does not mean that complex models may not be useful for research purposes, but more complicated does not mean a better model.

Figure 13.

Bacterial dynamics model compartments in the Snake River from Harrison [40].

Hence, to improve water quality models, one of the most fruitful areas is working on obtaining better boundary condition data by “smart” filling in of data gaps in time series of field data. This is still a critical component of modeling lakes and reservoirs. In addition, measuring field data on-site for lakes and reservoirs helps tremendously in understanding better the impact of hydrodynamics on water quality.


  1. 1. Wetzel R. Limnology Lake and River Ecosystems. NY: Academic Press; 2001
  2. 2. Hutchinson GE. A Treatise on Limnology. New York: John Wiley & Sons; 1967
  3. 3. Martin J, Higgins J, Edinger J, Gordon J. Energy Production and Reservoir Water Quality. Reston, VA: ASCE; 2007
  4. 4. Wells SA, Manson JR, Martin JL. Numerical hydrodynamic and transport models for reservoirs. In: Martin J, Higgins J, Edinger J, Gordon J, editors. Energy Production and Reservoir Water Quality. Reston, VA: ASCE; 2007. Chapter 4
  5. 5. Batchelor GK. An Introduction to Fluid Dynamics. NY: Cambridge University Press; 1967
  6. 6. Sabersky R, Acosta A, Haupmann E. Fluid Flow A First Course in Fluid Mechanics. NY: Macmillan Publishing Co.; 1989
  7. 7. Cushman-Roisin B. Introduction to Geophysical Fluid Dynamics. Englewood Cliffs, NJ: Prentice-Hall; 1994
  8. 8. Gill AE. “Appendix 3, Properties of Seawater”, Atmosphere Ocean Dynamics. Vol. 599. New York, NY: Academic Press; 1982. p. 600
  9. 9. Ford DE, Johnson MC. An Assessment of Reservoir Density Currents and Inflow Processes. Technical Rpt. E 83-7. Vicksburg, MS: US Army, Engineer Waterways Experiment Station; 1983
  10. 10. Turner JS. Buoyancy Effects in Fluids. Cambridge: Cambridge University Press; 1973
  11. 11. Rodi W. Turbulence models and their application in hydraulics. In: IAHR. 3rd ed. Rotterdam: A.A. Balkema; 1993
  12. 12. Venayagamoorthy SK, Koseff JR, Ferziger JH, Shih LH. Testing of RANS turbulence models for stratified flows based on DNS data. In: Annual Research Report Briefs. Center for Turbulence Research; 2003. pp. 127-138
  13. 13. Munk WH, Anderson ER. Notes on the theory of the thermocline. Journal of Marine Research. 1948;1:276
  14. 14. Mamayev OI. The influence of stratification on vertical turbulent mixing in the sea. Izv. Geophys. Ser. 1958:870-875
  15. 15. Fischer HB, List J, Koh R, Imberger J, Brooks N. Mixing in Inland and Coastal Waters. New York: Academic Press; 1979
  16. 16. Holland PR, Kay A, Botte V. A numerical study of the dynamics of the riverine thermal bar in a deep lake. Environmental Fluid Mechanics. 2001;1:311-332
  17. 17. Nezu I, Nakagawa H. Turbulence in Open-Channel Flows. Rotterdam: A.A. Balkema; 1993
  18. 18. Bloss S, Patterson JC. Modeling turbulent transport in a stratified estuary. ASCE Journal of Hydraulic Engineering. 1988;114(9):1115-1133
  19. 19. Martin JL, McCutcheon SC. Hydrodynamics and Transport for Water Quality Modeling. NY: Lewis Publishers; 1999
  20. 20. Wells SA. CE-QUAL-W2: A Two-Dimensional, Laterally Averaged, Hydrodynamic and Water Quality Model, Version 4.2, User Manual Part 2, Hydrodynamic and Water Quality Model Theory. Portland, OR: Department of Civil and Environmental Engineering, Portland State University; 2019
  21. 21. Brooks NH, Koh RC. Selective withdrawal from density-stratified reservoirs. Journal of Hydraulics Division, ASCE. 1969;95(HY4):1369-1397
  22. 22. Davis JE, Holland JP, Schneider ML, Wilhelms SC. SELECT: A Numerical, One-Dimensional Model for Selective Withdrawal, Instruction Report E-87-2. Waterways Experiment Station, Hydraulics Laboratory, Vicksburg, MS: Corps of Engineers; 1987
  23. 23. Chapra S. Surface Water Quality Modeling. NY: McGraw-Hill; 1997
  24. 24. Mooij WM, Trolle D, Jeppesen E, Arhonditsis G, Belolipetsky PV, Chitamwebwa DBR, et al. Challenges and opportunities for integrating lake ecosystem modelling approaches. Aquatic Ecology. 2010;44(3):633-667
  25. 25. Janssen ABG, Arhonditsis GB, Beusen A, Bolding K, Bruce L, Bruggeman J, et al. Exploring, exploiting and evolving diversity of aquatic ecosystem models: A community perspective. Aquatic Ecology. 2015;49(4):513-548. DOI: 10.1007/s10452-015-9544-1
  26. 26. Tanentzap AJ, Hamilton DP, Yan ND. Calibrating the dynamic reservoir simulation model (DYRESM) and filling required data gaps for one-dimensional thermal profile predictions in a boreal lake. Limnology and Oceanography: Methods. 2007;5:484-494
  27. 27. Environmental Laboratory. CE-QUAL-R1: A Numerical One-Dimensional Model of Reservoir Water Quality: Users Manual, Technical Report E-82-1. Vicksburg, MS: Corps of Engineers, Waterways Experiment Station, Environmental Laboratory; 1982
  28. 28. Al-Zubaidi HAM, Wells SA. Analytical and field verification of a 3D hydrodynamic and water quality numerical scheme based on the 2D formulation in CE-QUAL-W2. Journal of Hydraulic Research. 2020;58(1):152-171. DOI: 10.1080/00221686.2018.1499051
  29. 29. Hamrick JM. A three-dimensional environmental fluid dynamics computer code: Theoretical and computational aspects. In: Special Report. Vol. 317. The College of William and Mary, Virginia Institute of Marine Science; 1992. p. 63
  30. 30. Tetra Tech. The Environmental Fluid Dynamics Code User Manual US EPA Version 1.01. Fairfax, VA; 2007
  31. 31. Hipsey MR, Bruce LC, Hamilton DP. GLM—General lake model: Model overview and user information. In: AED Report #26. The University of Western Australia Technical Manual. 2014. p. 22
  32. 32. Hipsey MR, Romero J, Antenucci J, Hamilton D. Computational Aquatic Ecosystem Dynamics Model CAEDYM v2.0 Science Manual. Centre for Water Research, University of Western Australia; 2003
  33. 33. Hodges B, Dallimore C. Estuary and Lake Computer Model Science Manual Code Version 1.5.0. University of Western Australia: Centre for Water Research; 2001
  34. 34. EPA. Guidance on the development, evaluation, and application of environmental models. In: EPA/100/K-09/003. Washington, DC: Office of the Science Advisor; 2009
  35. 35. DeGasperi C. Modeling Lake Sammamish, WA Using a 2-D and 3-D Model, Presentation. Portland State University: Department of Civil Engineering; 2005
  36. 36. Geologic Survey Israel and Tahal. Red Sea Dead Sea Conveyance Study Program. GSI Report Number: GSI/10/2011. Israel, Jerusalem; 2011
  37. 37. Martinez VI, Wells SA, Addley RC. Meeting temperature requirements for fisheries downstream of Folsom Reservoir, California. In: EWRI, ASCE, Portland, OR: Proceedings World Environmental and Water Resources Congress. 2014. pp. 1081-1092
  38. 38. Cole T. Examples of model applications. In: Wells S, editor. CE-QUAL-W2: A Two-Dimensional, Laterally Averaged, Hydrodynamic and Water Quality Model, Version 4.2, User Manual Part 4, Model Examples. Portland, OR: Department of Civil and Environmental Engineering, Portland State University; 2019
  39. 39. Ceravich A, Wells SA. Water quality and hydrodynamic modeling of Chester Morse Lake. In: Technical Report Prepared for Seattle Public Utilities. WA, Seattle: Department of Civil and Environmental Engineering, Portland State University; 2020
  40. 40. Harrison J. Partitioning snake river organic matter and modeling aerobic oxidation. Brownlee Reservoir [PhD dissertation]. Idaho: Department of Civil Engineering, University of Idaho; 2005
  41. 41. Reichert P, Borchardt D, Henze M, Rauch W, Shanahan P, Somlyody L, et al. River water quality model no. 1. In: IWA Task Group on River Water Quality Modeling, editor. Scientific and Technical Report No. 12. London, UK: IWA Publishing; 2001. p. 131

Written By

Scott A. Wells

Submitted: November 27th, 2019 Reviewed: February 14th, 2020 Published: April 1st, 2020