Water Temperature Simulation Module

Water temperature is one of the most important physical characteristics of aquatic environments. In addition to its own direct effects, temperature influences many biological and chemical reactions in water quality models; almost all kinetic rates are temperature dependent. In the Water Temperature Simulation Module (WTSM), water temperature is calculated using an energy balance approach based on conservation of energy. This energy balance accounts for heat exchange at the water surface and at the sediment–water interface. The schematic below illustrates the sources and sinks of heat for a multi-layer water column. These include short-wave solar radiation, long-wave atmospheric radiation, long-wave radiation emitted by the water, latent (evaporative) and sensible heat fluxes, and heat exchange with the sediment layer. In the WTSM, the temperature of the sediment layer may be explicitly modeled, or set to a constant value.

Schematic of sources and sinks of heat exchange in the multi-layer water column

As illustrated above, the net heat flux ( $\begin{array}{l}q_{net}\end{array}$ ) for the surface layer is calculated as follows:

$\begin{array}{l}\displaystyle q_{net} = q_{sw} - (1 - \beta) q_{sw} e^{-\lambda h} + q_{atm} - q_b + q_h - q_l + q_{sed}\end{array}$

where

$\begin{array}{l}q_{sw}\end{array}$ is the short-wave solar radiation flux (W/m²),
$\begin{array}{l}q_{atm}\end{array}$ is the atmospheric (downwelling) long-wave radiation flux (W/m²),
$\begin{array}{l}q_b\end{array}$ is the back (upwelling) long-wave radiation flux (W/m²),
$\begin{array}{l}q_h\end{array}$ is the sensible heat flux (W/m²),
$\begin{array}{l}q_l\end{array}$ is the latent heat flux (W/m²),
$\begin{array}{l}q_{sed}\end{array}$ is the sediment–water heat flux (W/m²),
$\begin{array}{l}\beta\end{array}$ is the fraction of shortwave radiation absorbed at the water surface (unitless),
$\begin{array}{l}\lambda\end{array}$ is the light extinction coefficient (1/m), and
$\begin{array}{l}h\end{array}$ is the water depth of the surface layer (m).

A detailed discussion of each of terms and mechanisms can be found in Water Resources Engineers Inc. (1967), Brown and Barnwell (1987), and Deas and Lowney (2000). The radiative heat fluxes ( $\begin{array}{l}q_{sw}\end{array}$ , $\begin{array}{l}q_{atm}\end{array}$ , and $\begin{array}{l}q_b\end{array}$ ) will always be calculated as positive quantities. The remaining fluxes may be positive or negative, depending on the air vapor pressure deficit ( $\begin{array}{l}q_l\end{array}$ ), the difference between air and water temperatures ( $\begin{array}{l}q_h\end{array}$ ), or the difference between sediment and water temperatures ( $\begin{array}{l}q_{sed}\end{array}$ ).

The change in water temperature due to a change of the net heat flux is expressed as:

$\begin{array}{l}\displaystyle \frac{dT_w}{dt} = \frac{A_s}{V \rho_w C_{pw}} q_{net}\end{array}$

where

$\begin{array}{l}T_w\end{array}$ is the water temperature (°C),
$\begin{array}{l}\rho_w\end{array}$ is the density of water (kg/m³),
$\begin{array}{l}C_{pw}\end{array}$ is the specific heat capacity of water (J/kg/°C),
$\begin{array}{l}V\end{array}$ is the volume of a water quality cell or layer (m³), and
$\begin{array}{l}A_s\end{array}$ is the surface area of the water quality cell or layer (m²).

Shortwave Radiation Fluxes

Shortwave radiation refers to incoming solar radiation spanning roughly from 0.4 to 4.0 µm in wavelength. The magnitude of solar radiation reaching the water’s surface depends on the position of the sun in the Earth’s sky, the season, the time of the day, the site location, atmospheric attenuation, and cloud cover. Solar radiation is always positive during the day and zero during nighttime hours. The solar radiation flux that reaches the surface of the Earth can be directly measured. Once solar radiation reaches the water surface, a fraction is reflected back into the atmosphere, while the remainder penetrates the surface and is absorbed, changing the heat content of the water column.

Within the water column, short wave solar radiation decreases exponentially with depth according to Beer's Law

$\begin{array}{l}\displaystyle q_{sw}(d) = (1 - \beta) q_{sw} e^{-\lambda z}\end{array}$

where

$\begin{array}{l}q_{sw}\end{array}$ is the incident short-wave solar radiation (W/m²), and
$\begin{array}{l}q_{sw}(h)\end{array}$ is the short wave radiation at a depth $\begin{array}{l}h\end{array}$ below the water surface (W/m²).

In the WTSM, the light extinction coefficient is broken into three components

$\begin{array}{l}\displaystyle \lambda = \lambda_o + \lambda_s SSC + \lambda_p A_p\end{array}$

where

$\begin{array}{l}\lambda_o\end{array}$ is the (typically small) background light attenuation of pure water (≈ 0.04 1/m),
$\begin{array}{l}\lambda_s\end{array}$ is the light attenuation due to suspended sediment (L/mg/m),
$\begin{array}{l}SSC\end{array}$ is the suspended sediment concentration (mg/L),
$\begin{array}{l}\lambda_p\end{array}$ is the light attenuation due to phytoplankton (L/mg/m), and
$\begin{array}{l}A_p\end{array}$ is the phytoplankton concentration (mg/L).

In HEC-ResSim, the phytoplankton light attenuation component is ignored unless phytoplankton are explicitly simulated. Suspended sediment concentrations may be entered as a constant, defined with a seasonal table for values, or explicitly simulated using a state variable.

Longwave Radiation Fluxes

Longwave radiation includes wavelengths greater than 4.0 µm and varies according to temperature and humidity at low latitudes (John et al. 2006). Longwave radiation is categorized into two types: downwelling radiation ( $\begin{array}{l}q_{atm}\end{array}$ ), which is emitted by the atmosphere, and upwelling radiation ( $\begin{array}{l}q_b\end{array}$ ), which is emitted by the water surface. The amount of atmospheric longwave radiation is influenced by clouds and particles in the atmosphere. $\begin{array}{l}q_{atm}\end{array}$ is typically calculated using an empirical equation developed by Wunderlich (1972):

$\begin{array}{l}\displaystyle q_{atm} = 0.937 \times 10^{-5} (1 + 0.17 C_l^2) \, \sigma \, T^6_{ak}\end{array}$

where

$\begin{array}{l}\sigma\end{array}$ is the Stefan-Boltzman constant (W/m²/K⁴),
$\begin{array}{l}C_l\end{array}$ is the cloud cover (unitless, range 0-1), and
$\begin{array}{l}T_{ak}\end{array}$ is the air temperature in Kelvin.

Upwelling radiation emitted by the water’s surface represents a loss of heat from the water. The back (upwelling) long-wave radiation flux ( $\begin{array}{l}q_b\end{array}$ ) is typically calculated using the Stefan-Boltzman radiation law:

$\begin{array}{l}\displaystyle q_b = 0.97 \, \sigma \, T_{wk}^4\end{array}$

where

$\begin{array}{l}T_{wk}\end{array}$ is the water temperature in Kelvin.

Latent Heat Flux

Energy associated with a phase change is referred to as latent heat. A gain or loss of energy occurs during phase changes such as condensation or evaporation. The latent heat flux ( $\begin{array}{l}q_l\end{array}$ ) is primarily a function of the difference between the water vapor pressure in the air at saturation and the actual water vapor pressure in the air. It is modulated by a wind speed function. Latent heat flux is calculated using an equation provided by Brady et al. (1969):

$\begin{array}{l}\displaystyle q_l = 0.622 \frac{L_v \rho_w}{P} (e_s - e_a) f(u_w)\end{array}$

where

$\begin{array}{l}P\end{array}$ is the atmospheric pressure (mb),
$\begin{array}{l}L_v\end{array}$ is the latent heat of vaporization (J/kg), which is a function of water temperature,
$\begin{array}{l}e_s\end{array}$ is the saturation vapor pressure (mb),
$\begin{array}{l}e_a\end{array}$ is the vapor pressure of overlying air (mb),
$\begin{array}{l}u_w\end{array}$ is the wind speed measured at a fixed height above the water surface (m/s), and
$\begin{array}{l}f(u_w)\end{array}$ is a wind function (unitless).

The saturation vapor pressure is the maximum pressure of water vapor that can exist in equilibrium above a flat, free water surface at a given temperature. It is a function of water temperature and can be calculated using an empirical equation provided by Brutsaert (1982):

$\begin{array}{l}\displaystyle e_s = a_o + T_{wk} (a_1 + T_{wk} (a_2 + T_{wk} (a_3 + T_{wk} (a_4 + T_{wk} (a_5 + T_{wk} a_6)))))\end{array}$

The vapor pressure of air expresses its moisture content. Actual vapor pressure can be measured directly or calculated from the wet bulb temperature ( $\begin{array}{l}T_{wb}\end{array}$ ).

$\begin{array}{l}e_a = e_s(T_{wb}) - 0.00066(1 + 0.00115T_{wb} ) (T_a - T_{wb} ) P_{air} \\ e_s(T_{wb}) = 6.1078 e^{(17.269 T_{wb}) / (237.3 + T_{wb})}\end{array}$

If the dew point temperature ( $\begin{array}{l}T_{dp}\end{array}$ ) is known, the vapor pressure can also be calculated from:

$\begin{array}{l}\displaystyle e_a = 6.1078 e^{(17.269 T_{dp})/(237.3 + T_{dp})}\end{array}$

Relative humidity ( $\begin{array}{l}RH\end{array}$ ) is related to vapor pressure and saturation vapor pressure. If relative humidity is known, the vapor pressure is calculated from

$\begin{array}{l}\displaystyle e_a = RH \, e_s\end{array}$

The relative humidity increases the rate of heat transfer between the air and water while decreasing the rate of evaporation. The difference $\begin{array}{l}e_s - e_a\end{array}$ can be shown to be proportional to the difference $\begin{array}{l}T_w - T_{dp}\end{array}$ . When $\begin{array}{l}T_w > T_{dp}\end{array}$ , water evaporates and $\begin{array}{l}q_l\end{array}$ is positive (indicating a loss of heat from the water). Conversely, if $\begin{array}{l}T_w < T_{dp}\end{array}$ , water condenses on the surface and $\begin{array}{l}q_l\end{array}$ is negative (indicating a heat gain).

Sensible Heat Flux

Sensible heat, or heat conduction, describes the flux of heat through molecular or turbulent transfer between the air and the water surface. The sensible heat flux is primarily influenced by temperature and pressure differences, which are modulated by wind speed. The amount of heat gained or lost through sensible heat depends on the temperature gradient in the vertical direction. The sensible heat flux ( $\begin{array}{l}q_h\end{array}$ ) is typically computed as

$\begin{array}{l}\displaystyle q_h = \left( \frac{K_h}{K_w} \right) C_{pa} \rho_w (T_a - T_w) f(u_w)\end{array}$

where

$\begin{array}{l}C_{pa}\end{array}$ is the specific heat capacity of air at constant pressure (J/kg/°C)
$\begin{array}{l}T_a\end{array}$ is the air temperature (°C), and
$\begin{array}{l}\frac{K_h}{K_w}\end{array}$ is the diffusivity ratio (unitless).

The diffusivity ratio is a parameter that allows users to partition flux between latent and sensible heat. This ratio is generally set to unity, but the model allows it to range between 0.5 and 1.5, with a recommended range of 0.9 to 1.1.

The wind speed function aims to characterize the turbulent exchange characteristic between the water surface and the overlying air mass. The different formulations arise from the empirical determination of $\begin{array}{l}f(u_w)\end{array}$ for water bodies of varying size and shape, using data from different locations. While the wind function varies slightly for sensible and latent heat, the same function can typically be used for most water temperature applications. It is usually expressed as an empirical equation that can be adjusted using coefficients $\begin{array}{l}a\end{array}$ , $\begin{array}{l}b\end{array}$ , and $\begin{array}{l}c\end{array}$ .

$\begin{array}{l}\displaystyle f(u_w) = f(R_i) (a + b \, u^c_w)\end{array}$

where

$\begin{array}{l}a, b, c\end{array}$ are user-defined coefficients (unitless),
$\begin{array}{l}R_i\end{array}$ is the Richardson number, and
$\begin{array}{l}f(R_i)\end{array}$ is a function of the Richardson number.

The coefficient $\begin{array}{l}a\end{array}$ represents vertical convection that occurs even when wind speed is zero; it is typically small, and generally becomes significant only for artificially heated waters. In contrast, the coefficient $\begin{array}{l}b\end{array}$ increases with increasing turbulence and decreases in a stable atmosphere, varying by more than 50% (Fischer et al. 1979).

In the semi-empirical parameterizations of heat flux, wind speed is conventionally measured at a height of 2 meters. The following equation converts wind speed from any measurement height to 2 meters:

$\begin{array}{l}\displaystyle u_{w2} = u_w \frac{\ln \left( 2.0 / z_o \right)}{\ln \left( z / z_o \right)}\end{array}$

where

$\begin{array}{l}u_{w2}\end{array}$ is the wind speed measured at 2 m height (m/s)
$\begin{array}{l}u_w\end{array}$ is the wind speed measured at height z (m/s)
$\begin{array}{l}z\end{array}$ is the station height for wind velocity (m)
$\begin{array}{l}z_o\end{array}$ is the wind roughness height (m).

The Richardson number function ( $\begin{array}{l}f(R_i)\end{array}$ ) depends on air temperature, water temperature, and wind speed, varying from .03 under very stable conditions to 12.3 under unstable conditions. Without the Richardson number ( $\begin{array}{l}R_i\end{array}$ ) included in the wind function, the model tends to underestimate mixing processes under unstable atmospheric conditions, leading to an underprediction of heat flux, and overpredict surface fluxes under stable conditions.

The $\begin{array}{l}R_i\end{array}$ is a measure of atmospheric stability and can be computed as:

$\begin{array}{l}\displaystyle R_i = - \frac{2g(\rho_{air}-\rho_{sat})}{\rho_{air}u_w^2}\end{array}$

where

$\begin{array}{l}g\end{array}$ is the acceleration of gravity (9.806 m/s²),
$\begin{array}{l}\rho_{air}\end{array}$ is the density of the air (at the air temperature) (kg/m³), and
$\begin{array}{l}\rho_{sat}\end{array}$ is the density of saturated air (at the water temperature) (kg/m³).

The $\begin{array}{l}R_i\end{array}$ is positive for stable atmospheric conditions, negative for unstable conditions, and near zero for neutral conditions.

For an unstable atmosphere $\begin{array}{l}\rho_{air} > \rho_{sat}\end{array}$ and

$\begin{array}{l}f(R_i) = \begin{cases} 12.3 \quad &R_i \leq -1\\ (1 - 22 R_i)^{0.8} \quad & -1 < R_i \le -0.01 \end{cases}\end{array}$

For a neutral atmosphere $\begin{array}{l}\rho_{air} \approx \rho_{sat}\end{array}$ and

$\begin{array}{l}\displaystyle f(R_i) = 1 \quad -0.01 < R_i < 0.01\end{array}$

For a stable atmosphere $\begin{array}{l}\rho_{sat} > \rho_{air}\end{array}$ and

$\begin{array}{l}f(R_i) = \begin{cases} (1-34 R_i)^{-0.8} \quad & 0.01 \le R_i < 2 \\ 0.03 \quad &R_i \ge 2 \end{cases}\end{array}$

At least one complete meteorological data set must be provided when applying the full energy balance temperature simulation module. As a minimum, a time series of the following information at the meteorological station is required: short-wave solar radiation, atmospheric pressure, air temperature, humidity, wind speed, and cloud cover. Due to significant fluctuations in air temperature and solar radiation, hourly meteorological data are typically necessary.

Sediment Heat Exchange

Sediment heat exchange with the water is generally small compared to surface heat exchange and has often been neglected in modeling water temperature. However, heat exchange between benthic sediments and the water column is significant in shallow waters, and it is included in the full energy temperature simulation module. Sediment temperature is calculated from the sediment-water interface down to a user-defined depth. The only source or sink of the sediment layer temperature is the sediment-water heat flux. The change in temperature for the upper sediment layer is expressed as:

$\begin{array}{l}\displaystyle \frac{dT_{sed}}{dt} = - \frac{q_{sed}}{h_s \rho_s C_{ps}}\end{array}$

The sediment–water heat flux ( $\begin{array}{l}q_{sed}\end{array}$ ) is determined by

$\begin{array}{l}\displaystyle q_{sed} = \rho_s C_{ps} \frac{\alpha_s}{0.5 \, h_s}\left(T_{sed} - T_w \right)\end{array}$

where

$\begin{array}{l}T_{sed}\end{array}$ is the sediment temperature (°C),
$\begin{array}{l}h_s\end{array}$ is the active sediment layer thickness (m),
$\begin{array}{l}\alpha_s\end{array}$ is the sediment thermal diffusivity (m²/s),
$\begin{array}{l}\rho_s\end{array}$ is the density of the sediment layer (kg/m³), and
$\begin{array}{l}C_{ps}\end{array}$ is the specific heat capacity of the sediment layer (J/kg/°C).

Sediment–water heat flux ( $\begin{array}{l}q_{sed}\end{array}$ ) is a function of water temperature, sediment temperature, the heat storage capacity of sediment material, and the thermal diffusivity of sediment. The active sediment layer thickness ( $\begin{array}{l}h_s\end{array}$ ) is an user-defined input parameter, with typical values from 5 mm to 15 cm. Thermal conductivity ( $\begin{array}{l}\alpha_s\end{array}$ ) measures the ability of the bed sediment to conduct heat. It is not a constant for a particular bed but varies with both depth and time, depending upon sediment porosity and moisture content. Chapra et al. (2008) provided values for $\begin{array}{l}\alpha_s\end{array}$ ranging from 0.002 to 0.012 cm²/s, with a recommended value of 0.005 cm²/s.

In HEC-ResSim, water quality cells for river reaches are longitudinal segments of the river, and are assumed to have no water quality variation in the vertical. For these cells, sediment heat exchange is calculated using an area equal to the cell area used for surface heat exchange. Reservoirs are represented using vertical stacks of water quality cells (layers). For these cells, sediment heat exchange is calculated using an area that is the difference between the layer surface area and the surface area of the layer immediately below it. This can be seen in the illustration at the beginning of this section as small brown rectangles along the sides of the reservoir layers.

Impacts of Vegetated Stream Corridors

Stream corridors are often vegetated. Depending on the density of the streamside vegetation, there may be significant impacts on:

shortwave radiation, by vegetation shading the water from the sun,
evaporative and sensible heat fluxes, by vegetation sheltering the water surface from the wind, and
longwave radiation, by vegetation having an insulation effect on longwave fluxes.

Shading and sheltering factors have been frequently employed in stream temperature modeling (see, for example, Sinokrot and Stefan, 1993). The longwave radiation insulation factor is perhaps newer, and has been described by Rutherford, et al. (1997). In short, surface waters in riparian zones that have steep banks or are highly vegetated are exposed to longwave emissions from such surface features rather than the open sky. Consequently, their longwave emissions are shaded or reflected by those surface features, and longwave heat exchange - both atmospheric and upwelling - is dampened.

The impacts of vegetation on the various surface heat flux components have been integrated into the WTSM as user-defined coefficients that may be constant for a simulation, or may vary seasonally or as a function of a simulation state variable.

Shortwave shading is taken into account with a shading coefficient, $\begin{array}{l}f_{shade}\end{array}$

$\begin{array}{l}\displaystyle \hat{q}_{sw} = q_{sw} (1 - f_{shade})\end{array}$

where the $\begin{array}{l}\hat{}\end{array}$ represents the effect of vegetation on the component. Similarly, wind sheltering is taken into account with a sheltering coefficient, $\begin{array}{l}f_{shelt}\end{array}$

$\begin{array}{l}\displaystyle \hat{u}_w = u_w (1 - f_{shelt})\end{array}$

and longwave insulation is taken into account with an insulation coefficient, $\begin{array}{l}f_{insul}\end{array}$

$\begin{array}{l}\displaystyle \hat{q}_{atm} = q_{atm} (1 - f_{insul}), \quad \hat{q}_b = q_b (1 - f_{insul})\end{array}$

In all of these parameterizations, the vegetative factor $\begin{array}{l}f\end{array}$ is unitless and in the range 0-1, where a value of 0 represents no impacts from vegetation.