Modeling Thermal Stratification Effects in Lakes and Reservoirs

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.


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

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.

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 vector forces acting on control volume, m: mass within control volume, a ! : acceleration of fluid within 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.
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 • all velocities and pressure are turbulent time averages, i.e., u ¼ u þú, where udt andú 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 þẃ, and the pressure, p ¼ p þṕ The governing equations become after time averaging and simplifying: 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.

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 ∂η 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.

X-momentum equation
where: τ xx ¼ ρu 0 u 0 where τ xx is the turbulent shear stress acting in x direction on the x-face of control volume, τ xy ¼ ρu 0 v 0 where τ xy is the turbulent shear stress acting in x direction on the y-face of control volume, τ xz ¼ ρu 0 w 0 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: 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 where the term μ turbulent is the turbulent eddy viscosity analogous to molecular viscosity. The pressure is usually decomposed into the following terms: p ¼ where p a 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 where: τ yx ¼ ρv 0 u 0 where τ yx is the turbulent shear stress acting in y direction on the x-face of control volume, τ yy ¼ ρv 0 v 0 where τ yy is the turbulent shear stress acting in y direction on the y-face of control volume, τ yz ¼ ρv 0 w 0 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 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 where: τ zx ¼ ρw 0 u 0 where τ zx is the turbulent shear stress acting in z direction on the x-face of control volume, τ zy ¼ ρw 0 v 0 where τ zy is the turbulent shear stress acting in z direction on the y-face of control volume, τ zz ¼ ρw 0 w 0 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 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.

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 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, cdt and c 0 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: where the turbulent mass fluxes in x, y and z were assumed to be defined as a gradient, diffusion-type process, such asúć ∂z , E x is the turbulent mass diffusivity in x [L 2 /T], E y is the turbulent mass diffusivity in y[L 2 /T], E z is the turbulent mass diffusivity in z [L 2 /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, E x , E y , and E z >> 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, c p is the coefficient of specific heat at constant pressure and ρ is the density, such that the governing equation for temperature, T, becomes after simplification where D T is the molecular thermal conductivity for heat and E x , E y , and E z are the heat and mass turbulent eddy diffusivities assuming they are of the same order of magnitude.

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 Typical equations of state for fresh and saltwater have been published by Gill [8] and Ford and Johnson [9].

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 E x , E y and E z, (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.

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/m 3 /s and in the heat balance equation the units are [Energy L À3 T À1 ] with a typical unit of J/m 3 /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].

Typical source-sink term Description
Temperature S ¼ À ∂φ ∂z φ is the heat flux in units of W/m 2 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.
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 c BODd to N and δ CBODpN is the stoichiometric equivalent of c CBODp to N), and sediment oxygen demand release under anoxic conditions (where SOD N 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 k nitr [T À1 ], and c ammonia is the total ammonia concentration.
Dissolved oxygen S DO ¼ δ aO2 μ growth c algae À μ respiration c algae À k CBODd c CBODd À k CBODp c CBODp À SOD A V À δ NO2 k nitr c ammonia þ k reaeraton c s À c DO ð Þ The source/sink term includes algae production and respiration (where δ aO2 is the stoichiometric equivalent of dissolved oxygen to algae), CBOD particulate and dissolved watercolumn decay, sediment oxygen demand, nitrification demand (where δ NO2 is the stoichiometric equivalent of dissolved oxygen to N), and reaeration at the surface only (where k reaeration is the reaeration rate in [T À1 ] which is generally a function of wind speed in lake and reservoirs and c s 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 δ aNOx 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

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

State variable
Typical source-sink term Description under anoxic conditions only (where k denit is the denitrification rate under anoxic conditions). Other models also include terms for diffusion of nitrate into bottom muds. c NOX is the concentration of nitrite and nitrate.
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 c CBODd to P and δ CBODpP is the stoichiometric equivalent of c CBODp to P), and sediment oxygen demand release under anoxic conditions (where SOD P is the rate of P release in mass/ area/time and V is the volume of the computational cell, c algae is the algae concentration, and A is the area of the sediment). Other models include adsorption of P onto inorganic particles. 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 longterm changes in model boundary conditions.

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.
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 R 2 correlation between modeled and predicted temperature was 0.996.
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.

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