Qual2K: a modeling Framework for Simulating River and Stream Water Quality (Version 11)

Travel Time

The residence time of each element is computed as

()
where ?_k = the residence time of the k^th element [d], V_k = the volume of the k^th element [m³] = A_c,k?x_k, A_c,k = the cross-sectional area of the k^th element [m²], and ?x_k = the length of the k^th element [m]. These times are then accumulated to determine the travel time along each of the river’s segments (that is, either the main stem or one of the tributaries). For example, the travel time from the headwater to the downstream end of the jth element in a segment is computed as,
()
where t_t,j = the travel time [d].

4. Longitudinal Dispersion

Two options are used to determine the longitudinal dispersion for a boundary between two elements. First, the user can simply enter estimated values on the Reach Worksheet. If the user does not enter values, a formula is employed to internally compute dispersion based on the channel’s hydraulics (Fischer et al. 1979),

()
where E_p,i = the longitudinal dispersion between elements i and i + 1 [m²/s], U_i = velocity [m/s], B_i = width [m], H_i = mean depth [m], and U_i^* = shear velocity [m/s], which is related to more fundamental characteristics by
()
where g = acceleration due to gravity [= 9.81 m/s²] and S = channel slope [dimensionless].
After computing or prescribing E_p,i, the numerical dispersion is computed as
()
The model dispersion E_i (i.e., the value used in the model calculations) is then computed as follows:

• 
  If E_n,i ? E_p,i, the model dispersion, E_i is set to E_p,i ? E_n,i.
  
• 
  If E_n,i > E_p,i, the model dispersion is set to zero.
  

For the latter case, the resulting dispersion will be greater than the physical dispersion. Thus, dispersive mixing will be higher than reality. It should be noted that for most steady-state rivers, the impact of this overestimation on concentration gradients will be negligible. If the discrepancy is significant, the only alternative is to make element lengths smaller so that the numerical dispersion becomes smaller than the physical dispersion.

4. TEMPERATURE MODEL

As in Figure , the heat balance takes into account heat transfers from adjacent elements, loads, withdrawals, the atmosphere, and the sediments. A heat balance can be written for element i as

()
where T_i = temperature in element i [^oC], t = time [d], E’_i = the bulk dispersion coefficient between elements i and i + 1 [m³/d], W_h,i = the net heat load from point and non-point sources into element i [cal/d], ?_w = the density of water [g/cm³], C_pw = the specific heat of water [cal/(g ^oC)], J_a,i = the air-water heat flux [cal/(cm² d)], and J_s,i = the sediment-water heat flux [cal/(cm² d)].

Figure Heat balance for an element.

The bulk dispersion coefficient is computed as

()
Note that two types of boundary condition are used at the river’s downstream terminus: (1) a zero dispersion condition (natural boundary condition) and (2) a prescribed downstream boundary condition (Dirichlet boundary condition). The choice between these options is made on the Downstream Worksheet.
The net heat load from sources is computed as (recall Eq. 2)
()
where T_ps,i,j is the temperature of the jth point source for element i [^oC], and T_nps,i,j is the temperature of the jth non-point source temperature for element i [^oC].

1. Surface Heat Flux

As depicted in Figure , surface heat exchange is modeled as a combination of five processes:

()
where I(0) = net solar shortwave radiation at the water surface, J_an = net atmospheric longwave radiation, J_br = longwave back radiation from the water, J_c = conduction, and J_e = evaporation. All fluxes are expressed as cal/cm²/d.

Figure The components of surface heat exchange.

1. Solar Radiation

The model computes the amount of solar radiation entering the water at a particular latitude (L_at) and longitude (L_lm) on the earth’s surface. This quantity is a function of the radiation at the top of the earth’s atmosphere which is attenuated by atmospheric transmission, cloud cover, reflection, and shade,

()
where I(0) = solar radiation at the water surface [cal/cm²/d], I₀ = extraterrestrial radiation (i.e., at the top of the earth’s atmosphere) [cal/cm²/d], a_t = atmospheric attenuation, a_c = cloud attenuation, R_s = albedo (fraction reflected), and S_f = effective shade (fraction blocked by vegetation and topography).
Extraterrestrial radiation. The extraterrestrial radiation is computed as (TVA 1972)
()
where W₀ = the solar constant [1367 W/m² or 2823 cal/cm²/d], r = normalized radius of the earth’s orbit (i.e., the ratio of actual earth-sun distance to mean earth-sun distance), and ? = the sun’s altitude [radians], which can be computed as
()
where ? = solar declination [radians], L_at = local latitude [radians], and ? = the local hour angle of the sun [radians].
The local hour angle in radians is given by
()
where:
()
where trueSolarTime is the solar time determined from the actual position of the sun in the sky [minutes], localTime is the local time in minutes (local standard time), L_lm is the local longitude (positive decimal degrees for the western hemisphere), and timezone is the local time zone in hours relative to Greenwich Mean Time (e.g. ?8 hours for Pacific Standard Time; the local time zone is selected on the QUAL2K Worksheet). The value of eqtime represents the difference between true solar time and mean solar time in minutes of time.
QUAL2K calculates the solar declination, hour angle, solar altitude, and normalized radius (distance between the earth and sun), as well as the times of sunrise and sunset using the Meeus (1999) algorithms as implemented by NOAA’s Surface Radiation Research Branch (www.srrb.noaa.gov/highlights/sunrise/azel.html). The NOAA method for solar position that is used in QUAL2K also includes a correction for the effect of atmospheric refraction. The complete calculation method that is used to determine the solar position, sunrise, and sunset is presented in Appendix B.
The photoperiod f [hours] is computed as
()
where t_ss = time of sunset [hours] and t_sr = time of sunrise [hours].
Atmospheric attenuation. Various methods have been published to estimate the fraction of the atmospheric attenuation from a clear sky (a_t). Two alternative methods are available in QUAL2K to estimate a_t (Note that the solar radiation model is selected on the Light and Heat Worksheet of QUAL2K):

1) Bras (default)

The Bras (1990) method computes a_t as:

()
where n_fac is an atmospheric turbidity factor that varies from approximately 2 for clear skies to 4 or 5 for smoggy urban areas. The molecular scattering coefficient (a₁) is calculated as
()
where m is the optical air mass, calculated as
()
where ?_d is the sun’s altitude in degrees from the horizon = ? ? (180^o/?).
2) Ryan and Stolzenbach
The Ryan and Stolzenbach (1972) model computes a_t from ground surface elevation and solar altitude as:
()
where a_tc is the atmospheric transmission coefficient (0.70-0.91, typically approximately 0.8), and elev is the ground surface elevation in meters.
Direct measurements of solar radiation are available at some locations. For example, NOAA’s Integrated Surface Irradiance Study (ISIS) has data from various stations across the United States (http://www.atdd.noaa.gov/isis.htm). The selection of either the Bras or Ryan-Stolzenbach solar radiation model and the appropriate atmospheric turbidity factor or atmospheric transmission coefficient for a particular application should ideally be guided by a comparison of predicted solar radiation with measured values at a reference location.
Cloud Attenuation. Attenuation of solar radiation due to cloud cover is computed with
()
where C_L = fraction of the sky covered with clouds.
Reflectivity. Reflectivity is calculated as
()
where A and B are coefficients related to cloud cover (Table ).
Table Coefficients used to calculate reflectivity based on cloud cover.

Cloudiness	Clear	Scattered		Broken		Overcast
CL		0		0.1-0.5		0.5-0.9		1.0
Coefficients		A	B	A	B	A	B	A	B
		1.18	?0.77	2.20	?0.97	0.95	?0.75	0.35	?0.45

Shade. Shade is an input variable for the QUAL2K model. Shade is defined as the fraction of potential solar radiation that is blocked by topography and vegetation. An Excel/VBA program named ‘Shade.xls’ is available from the Washington Department of Ecology to estimate shade from topography and riparian vegetation (Ecology 2003). Input values of integrated hourly estimates of shade for each reach are entered on the Shade Worksheet of QUAL2K.

2. Atmospheric Long-wave Radiation

The downward flux of longwave radiation from the atmosphere is one of the largest terms in the surface heat balance. This flux can be calculated using the Stefan-Boltzmann law

()
where ? = the Stefan-Boltzmann constant = 11.7x10^-8 cal/(cm² d K⁴), T_air = air temperature [^oC], ε_sky = effective emissivity of the atmosphere [dimensionless], and R_L = longwave reflection coefficient [dimensionless]. Emissivity is the ratio of the longwave radiation from an object compared with the radiation from a perfect emitter at the same temperature. The reflection coefficient is generally small and is assumed to equal 0.03.
The atmospheric longwave radiation model is selected on the Light and Heat Worksheet of QUAL2K. Three alternative methods are available for use in QUAL2K to represent the effective emissivity (?_sky):
1) Brunt (default)
Brunt’s (1932) equation is an empirical model that has been commonly used in water-quality models (Thomann and Mueller 1987),

where A_a and A_b are empirical coefficients. Values of A_a have been reported to range from about 0.5 to 0.7 and values of A_b have been reported to range from about 0.031 to 0.076 mmHg^?0.5 for a wide range of atmospheric conditions. QUAL2K uses a default mid-range value of A_a = 0.6 together with a value of A_b = 0.031 mmHg^-0.5 if the Brunt method is selected on the Light and Heat Worksheet.
2) Brutsaert
The Brutsaert equation is physically-based instead of empirically derived and has been shown to yield satisfactory results over a wide range of atmospheric conditions of air temperature and humidity at intermediate latitudes for conditions above freezing (Brutsaert, 1982).

where e_air is the air vapor pressure [mm Hg], and T_a is the air temperature in °K. The factor of 1.333224 converts the vapor pressure from mm Hg to millibars. The air vapor pressure [in mm Hg] is computed as (Raudkivi 1979):
()
where T_d = the dew-point temperature [^oC].
3) Koberg
Koberg (1964) reported that the A_a in Brunt’s formula depends on both air temperature and the ratio of the incident solar radiation to the clear-sky radiation (R_sc). As in Figure , he presented a series of curves indicating that A_a increases with T_air and decreases with R_sc with A_b held constant at 0.0263 millibars^?0.5 (about 0.031 mmHg^?0.5).
The following polynomial is used in Q2K to provide a continuous approximation of Koberg’s curves.

where





The fit of this polynomial to points sampled from Koberg’s curves are depicted in Figure . Note that an upper limit of 0.735 is prescribed for A_a.

Figure The points are sampled from Koberg’s family of curves for determining the value of the A_a constant in Brunt’s equation for atmospheric longwave radiation (Koberg, 1964). The lines are the functional representation used in Q2K.
For cloudy conditions the atmospheric emissivity may increase as a result of increased water vapor content. High cirrus clouds may have a negligible effect on atmospheric emissivity, but lower stratus and cumulus clouds may have a significant effect. The Koberg method accounts for the effect of clouds on the emissivity of longwave radiation in the determination of the A_a coefficient. The Brunt and Brutsaert methods determine emissivity of a clear sky and do not account for the effect of clouds. Therefore, if the Brunt or Brutsaert methods are selected, then the effective atmospheric emissivity for cloudy skies (ε_sky) is estimated from the clear sky emissivity by using a nonlinear function of the fractional cloud cover (C_L) of the form (TVA, 1972):
()
The selection of the longwave model for a particular application should ideally be guided by a comparison of predicted results with measured values at a reference location. However, direct measurements are rarely collected. The Brutsaert method is recommended to represent a wide range of atmospheric conditions.

3. Water Long-wave Radiation

The back radiation from the water surface is represented by the Stefan-Boltzmann law,

()
where ? = emissivity of water (= 0.97) and T = the water temperature [^oC].

4. Conduction and Convection

Conduction is the transfer of heat from molecule to molecule when matter of different temperatures are brought into contact. Convection is heat transfer that occurs due to mass movement of fluids. Both can occur at the air-water interface and can be described by,
()
where c₁ = Bowen's coefficient (= 0.47 mmHg/^oC). The term, f(U_w), defines the dependence of the transfer on wind velocity over the water surface where U_w is the wind speed measured a fixed distance above the water surface.
Many relationships exist to define the wind dependence. Bras (1990), Edinger et al. (1974), and Shanahan (1984) provide reviews of various methods. Some researchers have proposed that conduction/convection and evaporation are negligible in the absence of wind (e.g. Marciano and Harbeck, 1952), which is consistent with the assumption that only molecular processes contribute to the transfer of mass and heat without wind (Edinger et al. 1974). Others have shown that significant conduction/convection and evaporation can occur in the absence of wind (e.g. Brady Graves and Geyer 1969, Harbeck 1962, Ryan and Harleman 1971, Helfrich et al. 1982, and Adams et al. 1987). This latter opinion has gained favor (Edinger et al. 1974), especially for waterbodies that exhibit water temperatures that are greater than the air temperature.
Brady, Graves, and Geyer (1969) pointed out that if the water surface temperature is warmer than the air temperature, then “the air adjacent to the water surface would tend to become both warmer and more moist than that above it, thereby (due to both of these factors) becoming less dense. The resulting vertical convective air currents … might be expected to achieve much higher rates of heat and mass transfer from the water surface [even in the absence of wind] than would be possible by molecular diffusion alone” (Edinger et al. 1974). Water temperatures in natural waterbodies are frequently greater than the air temperature, especially at night.
Edinger et al. (1974) recommend that the relationship that was proposed by Brady, Graves and Geyer (1969) based on data from cooling ponds, could be representative of most environmental conditions. Shanahan (1984) recommends that the Lake Hefner equation (Marciano and Harbeck, 1952) is appropriate for natural waters in which the water temperature is less than the air temperature. Shanahan also recommends that the Ryan and Harleman (1971) equation as recalibrated by Helfrich et al. (1982) is best suited for waterbodies that experience water temperatures that are greater than the air temperature. Adams et al. (1987) revisited the Ryan and Harleman and Helfrich et al. models and proposed another recalibration using additional data for waterbodies that exhibit water temperatures that are greater than the air temperature.
Three options are available on the Light and Heat Worksheet in QUAL2K to calculate f(U_w):
1) Brady, Graves, and Geyer (default)

where U_w = wind speed at a height of 7 m [m/s].
2) Adams 1
Adams et al. (1987) updated the work of Ryan and Harleman (1971) and Helfrich et al. (1982) to derive an empirical model of the wind speed function for heated waters that accounts for the enhancement of convection currents when the virtual temperature difference between the water and air (??_v in degrees F) is greater than zero. Two wind functions reported by Adams et al., also known as the East Mesa method, are implemented in QUAL2K (wind speed in these equations is at a height of 2m).
This formulation uses an empirical function to estimate the effect of convection currents caused by virtual temperature differences between water and air, and the Harbeck (1962) equation is used to represent the contribution to conduction/convection and evaporation that is not due to convection currents caused by high virtual water temperature.

where U_w,mph is wind speed in mph and A_acres,i is surface area of element i in acres. The constant 0.271 converts the original units of BTU ft^?2 day^?1 mmHg^?1 to cal cm^?2 day^?1 mmHg^?1.
3) Adams 2
This formulation uses an empirical function of virtual temperature differences with the Marciano and Harbeck (1952) equation for the contribution to conduction/convection and evaporation that is not due to the high virtual water temperature

Virtual temperature is defined as the temperature of dry air that has the same density as air under the in situ conditions of humidity. The virtual temperature difference between the water and air ( in °F) accounts for the buoyancy of the moist air above a heated water surface. The virtual temperature difference is estimated from water temperature (T_w,f in °F), air temperature (T_air,f in °F), vapor pressure of water and air (e_s and e_air in mmHg), and the atmospheric pressure (p_atm is estimated as standard atmospheric pressure of 760 mmHg in QUAL2K):
()
The height of wind speed measurements is also an important consideration for estimating conduction/convection and evaporation. QUAL2K internally adjusts the wind speed to the correct height for the wind function that is selected on the Light and Heat Worksheet. The input values for wind speed on the Wind Speed Worksheet in QUAL2K are assumed to be representative of conditions at a height of 7 meters above the water surface. To convert wind speed measurements (U_w,z in m/s) taken at any height (z_w in meters) to the equivalent conditions at a height of z = 7 m for input to the Wind Speed Worksheet of QUAL2K, the exponential wind law equation may be used (TVA, 1972):
()
For example, if wind speed data were collected from a height of 2 m, then the wind speed at 7 m for input to the Wind Speed Worksheet of QUAL2K would be estimated by multiplying the measured wind speed by a factor of 1.2.

5. Evaporation and Condensation

The heat loss due to evaporation can be represented by Dalton’s law,

()
where e_s = the saturation vapor pressure at the water surface [mmHg], and e_air = the air vapor pressure [mmHg]. The saturation vapor pressure is computed as
()

2.