Open access peer-reviewed chapter

Modeling Thermal Stratification Effects in Lakes and Reservoirs

By Scott A. Wells

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

DOI: 10.5772/intechopen.91754

Downloaded: 360


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=ma, 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 Δρρ<<1where ρ 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ρowhere ρ=ρo+ρ, ρois a base value,

  • all velocities and pressure are turbulent time averages, i.e., u=u¯+u´, where u¯=1Ttt+Tudtand 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


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ηhu¯dz+yηhv¯dzηhqdzwhere 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


where: τxx=ρuu¯where τxx is the turbulent shear stress acting in x direction on the x-face of control volume, τxy=ρuv¯where τxy is the turbulent shear stress acting in x direction on the y-face of control volume, τxz=ρuw¯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=ΩEsinφ, Ωy=ΩEcosφ, ϕ=latitude, ΩE=earthsrotation rate, and assuming 2Ωyw¯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=μturbulentxxu¯x=ρuu¯, τxx=μturbulentxyu¯y=ρuv¯, τxz=μturbulentxzu¯z=ρuw¯where the term μturbulentis the turbulent eddy viscosity analogous to molecular viscosity. The pressure is usually decomposed into the following terms: p¯=p¯a+gηzρ¯dzwhere 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¯ax+gη¯xgρηzρ¯xdz.

2.1.3 Y-momentum equation


where: τyx=ρvu¯where τyx is the turbulent shear stress acting in y direction on the x-face of control volume, τyy=ρvv¯where τyy is the turbulent shear stress acting in y direction on the y-face of control volume, τyz=ρvw¯where τyz is the turbulent shear stress acting in y direction on the z-face of control volume, and assuming 2Ωxw¯is negligible. Analogous to laminar shear stress, the turbulent shear stresses are often parameterized as τyx=μturbulentyxv¯x=ρvu¯, τyy=μturbulentyyv¯y=ρvv¯, τyz=μturbulentxzv¯z=ρvw¯. 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¯ay+gη¯ygρηzρ¯ydz.

2.1.4 Z-momentum equation


where: τzx=ρwu¯where τzx is the turbulent shear stress acting in z direction on the x-face of control volume, τzy=ρwv¯where τzy is the turbulent shear stress acting in z direction on the y-face of control volume, τzz=ρww¯where τzz is the turbulent shear stress acting in z direction on the z-face of control volume, and neglecting the Coriolis terms 2Ωyu¯+2Ωxv¯. Analogous to laminar shear stress, the turbulent shear stresses are often parameterized as τzx=μturbulentzxw¯x=ρwu¯, τzy=μturbulentzyw¯y=ρwv¯, τzz=μturbulentzzw¯z=ρww¯. 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

ctunsteady changein concentration+ucx+vcy+wczadvective mass transport=D2cx2+2cy2+2cz2diffusive mass transport+Ssources/sinksE5

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¯+cwhere c¯=1Ttt+Tcdtand cis 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¯tunsteady changein concentration+u¯c¯x+v¯c¯y+w¯c¯zmean advective mass transport=D2c¯x2+2c¯y2+2c¯z2molecular diffusive mass transport+xExc¯x+yEyc¯y+zEzc¯zturbulent diffusive mass transport+S¯sources/sinksE6

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´¯=Exc¯x, v´c´¯=Eyc¯y, w´c´¯=Ezc¯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, ρcpT, 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¯tunsteady changein temperature+u¯T¯x+v¯T¯y+w¯T¯zmean advective heat transport=DT2T¯x2+2T¯y2+2T¯z2molecular diffusive heat transport+xExT¯x+yEyT¯y+zEzT¯zturbulent diffusive heat transport+S¯ρcpheat fluxE7

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 (cdissolved solids,csuspended solids)and temperature, T, such as

ρ¯=fT¯cdissolved solids¯csuspended 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¯orT¯,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 turbulentconductivity 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 variableTypical source-sink termDescription
TemperatureS¯=φ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 substanceS¯=0No sources and sinks
Suspended solidsS¯SS=wSScSSzwss 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.
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.
AlgaeSalgae¯=μgrowthcalgaeμrespirationcalgaeμexcretioncalgaeμmortalitycalgaewalgaecalgaezSource 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-NSammonia¯=δaNμgrowthcalgae+μrespirationcalgae+δCBODdNkCBODdcCBODd+δCBODpNkCBODpcCBODp+SODNAVknitrcammoniaThe 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 oxygenSDO¯=δaO2μgrowthcalgaeμrespirationcalgaekCBODdcCBODdkCBODpcCBODpSODAVδNO2knitrcammonia+kreaeratoncscDOThe source/sink term includes algae production and respiration (where δaO2is the stoichiometric equivalent of dissolved oxygen to algae), CBOD particulate and dissolved water-column decay, sediment oxygen demand, nitrification demand (where δNO2is 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-NSNOx¯=knitrcammoniaδaNOxμgrowthcalgaekdenitcNOxThe source/sink terms include algae uptake (where δaNOxis 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-PSPO4¯=δaPμgrowthcalgae+μrespirationcalgae+δCBODdPkCBODdcCBODd+δCBODpPkCBODpcCBODp+SODPAVThe 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 nameDescriptionReference
DYRESM and CAEDYM1D model based on mixed layer dynamics, separate temperature and water quality modelsTanentzap et al. [26]
CE-QUAL-W22D longitudinal-vertical, open source, eutrophication model, hydrodynamics and water quality solved togetherWells [20]
CE-QUAL-R11D verticalEnvironmental Laboratory [27]
W33D, hydrodynamics and water quality solved togetherAl-Zubaidi and Wells [28]
EFDC and WASP3D, hydrodynamics and water quality solved separately, both sigma stretch and z coordinate modelsHamrick [29], Tetra Tech [30]
GLM1DHipsey et al. [31]
ELCOM and CAEDYM3D-mixed layer dynamic model, hydrodynamics and water quality solved separatelyHipsey 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.

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Scott A. Wells (April 1st 2020). Modeling Thermal Stratification Effects in Lakes and Reservoirs, Inland Waters - Dynamics and Ecology, Adam Devlin, Jiayi Pan and Mohammad Manjur Shah, IntechOpen, DOI: 10.5772/intechopen.91754. Available from:

chapter statistics

360total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Assessment of the CHIRPS-Based Satellite Precipitation Estimates

By Franklin Paredes-Trejo, Humberto Alves Barbosa, Tumuluru Venkata Lakshmi Kumar, Manoj Kumar Thakur and Catarina de Oliveira Buriti

Related Book

First chapter

Climatic Change in a Large Shallow Tropical Lake Chapala, Mexico

By Filonov Anatoliy, Iryna Tereshchenko, Cesar Monzon, David Avalos-Cueva and Diego Pantoja-González

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us