Typical source-sink terms for temperature and some eutrophication water quality state variables.
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
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  and Hutchinson , 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. .
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.
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.  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 , where : vector forces acting on control volume, m: mass within control volume, : 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 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. where , is a base value,
all velocities and pressure are turbulent time averages, i.e., , where and is the temporal fluctuation of u about the mean, and similarly for the velocity in the y-direction, , the velocity in the z direction , and the pressure,
The governing equations become after time averaging and simplifying:
where : temporal mean velocity in the x-direction, : temporal mean velocity in the y-direction, : 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.
2.1.2 X-momentum equation
where: where τxx is the turbulent shear stress acting in x direction on the x-face of control volume, where τxy is the turbulent shear stress acting in x direction on the y-face of control volume, 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: , , , , and assuming 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 is the turbulent eddy viscosity analogous to molecular viscosity. The pressure is usually decomposed into the following terms: 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 .
2.1.3 Y-momentum equation
where: where τyx is the turbulent shear stress acting in y direction on the x-face of control volume, where τyy is the turbulent shear stress acting in y direction on the y-face of control volume, where τyz is the turbulent shear stress acting in y direction on the z-face of control volume, and assuming is negligible. Analogous to laminar shear stress, the turbulent shear stresses are often parameterized as , , . The pressure is usually decomposed into the following terms: , and the pressure gradient in the y-momentum then becomes after simplification .
2.1.4 Z-momentum equation
where: where τzx is the turbulent shear stress acting in z direction on the x-face of control volume, where τzy is the turbulent shear stress acting in z direction on the y-face of control volume, where τzz is the turbulent shear stress acting in z direction on the z-face of control volume, and neglecting the Coriolis terms . 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., .
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
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 where and 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 , , , 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, , 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
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 (and temperature, T, such as
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. . 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 () and the turbulent Prandtl number (). 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 , (4) Two-equation k-ε models for turbulent kinetic energy and dissipation  and (5) Reynolds stress and algebraic stress models . In many models, once the turbulent eddy viscosity is known, then the turbulent diffusion coefficients are computed from 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 .
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 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  and Chapra .
|State variable||Typical source-sink term||Description|
|Temperature||φ 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||No sources and sinks|
|Suspended solids||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||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||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||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||The source/sink term includes algae production and respiration (where is the stoichiometric equivalent of dissolved oxygen to algae), CBOD particulate and dissolved water-column decay, sediment oxygen demand, nitrification demand (where 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||The source/sink terms include algae uptake (where 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||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.|
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.  and Janssen et al. . Table 2 shows a listing of some common lake and reservoir models.
|DYRESM and CAEDYM||1D model based on mixed layer dynamics, separate temperature and water quality models||Tanentzap et al. |
|CE-QUAL-W2||2D longitudinal-vertical, open source, eutrophication model, hydrodynamics and water quality solved together||Wells |
|CE-QUAL-R1||1D vertical||Environmental Laboratory |
|W3||3D, hydrodynamics and water quality solved together||Al-Zubaidi and Wells |
|EFDC and WASP||3D, hydrodynamics and water quality solved separately, both sigma stretch and z coordinate models||Hamrick , Tetra Tech |
|GLM||1D||Hipsey et al. |
|ELCOM and CAEDYM||3D-mixed layer dynamic model, hydrodynamics and water quality solved separately||Hipsey et al. , Hodges and Dallimore |
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 . 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
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 R2 correlation between modeled and predicted temperature was 0.996.
In other examples of predicting the thermal regime, Cole  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  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.
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 . The bacterial populations were modeled based on Reichert et al.  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.