1 SEGMENTATION AND HYDRAULICS
The model represents a river as a series of reaches. These represent stretches of river that have constant hydraulic characteristics (e.g., slope, bottom width, etc.). As depicted in Figure 1, the reaches are numbered in ascending order starting from the headwater of the river‟s main stem.
Notice that both point and non-point sources and point and non-point withdrawals (abstractions) can be positioned anywhere along the channel‟s length.
Point source2Point withdrawal
Non-point7source Figure 1 QUAL2K segmentation scheme for a river with no tributaries. Point source8For systems with tributaries (Figure 2), the reaches are numbered in ascending order starting at reach 1 at the headwater of the main stem. When a junction with a tributary is reached, the Downstream boundaryHW#1numbering continues at that tributary‟s headwater. Observe that both the headwaters and the
tributaries are also numbered consecutively following a sequencing scheme similar to the reaches. 11HW#2Note also that the major branches of the system (that is, the main stem and each of the tributaries) 6622are referred to as segments. This distinction has practical importance because the software
77provides plots of model output on a segment basis. That is, the software generates individual plots 3388for the main stem as well as each of the tributaries. 12124413131111551414Trib110101515991616HW#3Trib21717
2727 2828Figure 2 QUAL2K segmentation scheme for (a) a river with tributaries. The Q2K reach
representation in (b) illustrates the reach, headwater and tributary numbering 2929schemes.
Finally, any model reach can be further divided into a series of equally-spaced elements. As in Figure 3, this is done by merely specifying the number of elements that are desired. (a) A river with tributaries(b) Q2K reach representation
QUAL2K 1 December 16, 2008
nn= 4= 4
Figure 3 If desired, any model reach can be further subdivided into a series of n equal-
In summary, the nomenclature used to describe the way in which Q2K organizes river
topology is as follows:
? Reach. A length of river with constant hydraulic characteristics.
? Element. The model‟s fundamental computational unit which consists of an equal length
subdivision of a reach.
? Segment. A collection of reaches representing a branch of the system. These consist of
the main stem as well as each tributary.
? Headwater. The upper boundary of a model segment. 1.1 Flow Balance
As described in the last section, Q2K‟s most fundamental unit is the element. A steady-state flow
balance is implemented for each model element as (Figure 4)
Q！Q？Q，Q (1) ii，1in,iout,i3where Q = outflow from element i into the downstream element i + 1 [m/d], Q= inflow from ii–1 3the upstream element i – 1 [m/d], Q is the total inflow into the element from point and in,i3nonpoint sources [m/d], and Q is the total outflow from the element due to point and nonpoint out,i3withdrawals [m/d]. Thus, the downstream outflow is simply the difference between inflow and
source gains minus withdrawal losses.
Figure 4 Element flow balance.
The total inflow from sources is computed as
(2) Q！Q？Q,,,,,??inipsijnpsij！1！1jj3where Q is the jth point source inflow to element i [m/d], psi = the total number of point ps,i,j3sources to element i, Q is the jth non-point source inflow to element i [m/d], and npsi = the nps,i,j
total number of non-point source inflows to element i.
The total outflow from withdrawals is computed as
QUAL2K 2 December 16, 2008
3where Q is the jth point withdrawal outflow from element i [m/d], pai = the total number of pa,i,j
point withdrawals from element i, Q is the jth non-point withdrawal outflow from element i npa,i,j3[m/d], and npai = the total number of non-point withdrawal flows from element i.
The non-point sources and withdrawals are modeled as line sources. As in Figure 5, the non-
point source or withdrawal is demarcated by its starting and ending kilometer points. Its flow is
then distributed to or from each element in a length-weighted fashion.
startend Figure 5 The manner in which non-point source flow is distributed to an element. 1.2 Hydraulic Characteristics
Once the outflow for each element is computed, the depth and velocity are calculated in one of
three ways: weirs, rating curves, and Manning equations. The program decides among these options in the following manner:
? If weir height and width are entered, the weir option is implemented.
? If the weir height and width are zero and rating curve coefficients are entered (a and ；),
the rating curve option is implemented.
? If neither of the previous conditions is met, Q2K uses the Manning equation.
Figure 6 shows how weirs are represented in Q2K. Note that a weir can only occur at the end of a reach consisting of a single element. The symbols shown in Figure 6 are defined as: H = i
the depth of the element upstream of the weir [m], H = the depth of the element downstream of i+1
the weir [m], elev2 = the elevation above sea level of the tail end of the upstream element [m], i
elev1 = the elevation above sea level of the head end of the downstream element [m], H = the i+1w
height of the weir above elev2 [m], H = the drop between the elevation above sea level of the id
surface of element i and element i+1 [m], H = the head above the weir [m], B = the width of the hw
weir [m]. Note that the width of the weir can differ from the width of the element, B. i
(a) Side(b) Cross-section
Figure 6 A sharp-crested weir occurring at the boundary between two reaches. QUAL2K 3 December 16, 2008
For a sharp-crested weir where H/H < 0.4, flow is related to head by (Finnemore and hw
3/2Q！1.83BH (4) iwh3where Q is the outflow from the element upstream of the weir in m/s, and B and H are in m. iwh
Equation (4) can be solved for
2/3??Qi?? (5) H！h??1.83Bw??
This result can then be used to compute the depth of element i,
and the drop over the weir
H！elev2？H，elev1，H (7) diii？1i？1
Note that this drop is used to compute oxygen and carbon dioxide gas transfer due to the weir (see
pages 43 and 49).
The cross-sectional area, velocity, surface area and volume of element i can then be
A！BH (8) c,iii
Qi！ (9) UiAc,i
where B = the width of element i, ；x = the length of element i. Note that for reaches with weirs, ii
the reach width must be entered. This value is entered in the column AA (labeled "Bottom
Width") of the Reach Worksheet.
1.2.1 Rating Curves
Power equations (sometimes called Leopold-Maddox relationships) can be used to relate mean
velocity and depth to flow for the elements in a reach,
where a, b, ； and ~ are empirical coefficients that are determined from velocity-discharge and stage-discharge rating curves, respectively. The values of velocity and depth can then be
employed to determine the cross-sectional area and width by
Q (12) A！cU
Ac (13) ！BH
The surface area and volume of the element can then be computed as
The exponents b and ~ typically take on values listed in Table 1. Note that the sum of b and ~
must be less than or equal to 1. If this is not the case, the width will decrease with increasing flow.
If their sum equals 1, the channel is rectangular.
Table 1 Typical values for the exponents of rating curves used to determine velocity and
depth from flow (Barnwell et al. 1989).
Equation Exponent Typical value Range
bb 0.43 0.4，0.6 U！aQ
QUAL2K 4 December 16, 2008
~0.45 0.3，0.5 ~;H！；Q
In some applications, you might want to specify constant values of depth and velocity that do not vary with flow. This can be done by setting the exponents b and ~ to zero and setting a equal
to the desired velocity and ； equal to the desired depth.
1.2.2 Manning Equation
Each element in a particular reach can be idealized as a trapezoidal channel (Figure 7). Under conditions of steady flow, the Manning equation can be used to express the relationship between flow and depth as
1/25/3SA0c！Q (14) 2/3nP31where Q = flow [m/s], S = bottom slope [m/m], n = the Manning roughness coefficient, A = 0c2the cross-sectional area [m], and P = the wetted perimeter [m].
Figure 7 Trapezoidal channel.
The cross-sectional area of a trapezoidal channel is computed as
，：A！B？0.5(s？s)HH (15) c0s1s2
where B = bottom width [m], s and s = the two side slopes as shown in Figure 7 [m/m], and H 0s1s2
= element depth [m].
The wetted perimeter is computed as
22 (16) P！B？Hs？1？Hs？1ss012
After substituting Eqs. (15) and (16), Eq. (14) can be solved iteratively for depth (Chapra and Canale 2006),
2/53/522??(Qn)B？Hs？1？Hs？1??0k，1s1k，1s2??H！ (17) k3/10，：SB？0.5(s？s)H0s1s2k，1
where k = 1, 2, …, n, where n = the number of iterations. An initial guess of H = 0 is employed. 0
The method is terminated when the estimated error falls below a specified value of 0.001%. The estimated error is calculated as
H，Hk？1k?！?100% (18) aHk？1
The cross-sectional area is determined with Eq. (15) and the velocity can then be determined from the continuity equation,
QU！ (19) Ac
1 Notice that time is measured in seconds in this and other formulas used to characterize hydraulics. This is how the computations are implemented within Q2K. However, once the hydraulic characteristics are determined they are converted to day units to be compatible with other computations. QUAL2K 5 December 16, 2008
The average element width, B [m], is computed as
Ac (20) ！BH
The top width, B [m], is computed as 1
The surface area and volume of the element can then be computed as
Suggested values for the Manning coefficient are listed in Table 2. Manning‟s n typically
varies with flow and depth (Gordon et al. 1992). As the depth decreases at low flow, the relative roughness usually increases. Typical published values of Manning‟s n, which range from about 0.015 for smooth channels to about 0.15 for rough natural channels, are representative of conditions when the flow is at the bankfull capacity (Rosgen, 1996). Critical conditions of depth for evaluating water quality are generally much less than bankfull depth, and the relative roughness may be much higher.
Table 2 The Manning roughness coefficient for various open channel surfaces (from Chow
et al. 1988).
Gravel bottom with sides:
mortared stone 0.023
Natural stream channels
Clean, straight 0.025-0.04
Clean, winding and some weeds 0.03-0.05
Weeds and pools, winding 0.05
Mountain streams with boulders 0.04-0.10
Heavy brush, timber 0.05-0.20
In Section 0, the drop of water over a weir was computed. This value is needed in order to compute the enhanced reaeration that occurs in such cases. In addition to weirs, such drops can also occur at waterfalls (Figure 8). Note that waterfalls can only occur at the end of a reach.
Figure 8 A waterfall occurring at the boundary between two reaches.
QUAL2K computes such drops for cases where the elevation above sea level drops abruptly at the boundary between two reaches. Equation (7) is used to compute the drop. It should be noted that the drop is only calculated when the elevation above sea level at the downstream end of a reach is greater than at the beginning of the next downstream reach; that is, elev2 > elev1. ii+1
QUAL2K 6 December 16, 2008
1.3 Travel Time
The residence time of each element is computed as
Vk： (21) ！kQkthth3where ： = the residence time of the k element [d], V = the volume of the k element [m] = kkth2thA；x, A = the cross-sectional area of the k element [m], and ；x = the length of the k c,kkc,kk
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,
t！： (22) tjk,?k！1
where t = the travel time [d]. t,j
1.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),
22UBiiE！0.011 (23) p,i*HUii
2where E = the longitudinal dispersion between elements i and i + 1 [m/s], U = velocity [m/s], p,ii*B = width [m], H = mean depth [m], and U = shear velocity [m/s], which is related to more iii
fundamental characteristics by
* (24) U！gHSiii2where g = acceleration due to gravity [= 9.81 m/s] and S = channel slope [dimensionless].
After computing or prescribing E, the numerical dispersion is computed as p,i
U；xii (25) E！n,i2
The model dispersion E (i.e., the value used in the model calculations) is then computed as i
? If E ? E, the model dispersion, E is set to E ， E. n,ip,iip,in,i
? If E > E, the model dispersion is set to zero. n,ip,i
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.
QUAL2K 7 December 16, 2008
2 TEMPERATURE MODEL
As in Figure 9, 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
3WJJ??mmm????h,ia,is,i?? ？ ？？????63??，CV，CH100 cm，CH100 cm10 cm??????wpwiwpwiwpwiowhere T = temperature in element i [C], t = time [d], E‟ = the bulk dispersion coefficient ii3between elements i and i + 1 [m/d], W = the net heat load from point and non-point sources into h,i3oelement i [cal/d], ， = the density of water [g/cm], C = the specific heat of water [cal/(g C)], wpw22J = the air-water heat flux [cal/(cm d)], and J = the sediment-water heat flux [cal/(cm d)]. a,is,iatmospheric
transferheat loadheat withdrawal
sediment Figure 9 Heat balance for an element.
The bulk dispersion coefficient is computed as
EAic,i'E！ (27) i；？；x？；x/2ii？1
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
The net heat load from sources is computed as (recall Eq. 2)
??W！，CQT？QT (28) h,ipps,i,jpsi,jnps,i,jnpsi,j????j！1j！1??owhere T is the temperature of the jth point source for element i [C], and T is the ps,i,jnps,i,jotemperature of the jth non-point source temperature for element i [C].
2.1 Surface Heat Flux
As depicted in Figure 10, surface heat exchange is modeled as a combination of five processes:
J！ I(0)？J，J，J，J (29) hanbrce
where I(0) = net solar shortwave radiation at the water surface, J = net atmospheric longwave an
radiation, J = longwave back radiation from the water, J = conduction, and J = evaporation. brce2All fluxes are expressed as cal/cm/d.
QUAL2K 8 December 16, 2008
non-radiation termsradiation terms
net absorbed radiationwater-dependent terms Figure 10 The components of surface heat exchange.
2.1.1 Solar Radiation
The model computes the amount of solar radiation entering the water at a particular latitude (L) at
and longitude (L) on the earth‟s surface. This quantity is a function of the radiation at the top of lm
the earth‟s atmosphere which is attenuated by atmospheric transmission, cloud cover, reflection, and shade,
(0) (1) (1)I！Iaa，R，S0tcsf
(30) extraterrestrial atmospheric cloud reflection shading
radiation attenuation attenuation2where I(0) = solar radiation at the water surface [cal/cm/d], I = extraterrestrial radiation (i.e., at 02the top of the earth‟s atmosphere) [cal/cm/d], a = atmospheric attenuation, a = cloud attenuation, tc
R = albedo (fraction reflected), and S = effective shade (fraction blocked by vegetation and sf
Extraterrestrial radiation. The extraterrestrial radiation is computed as (TVA 1972)
W0I！sin； (31) 02r22where W = the solar constant [1367 W/m or 2823 cal/cm/d], r = normalized radius of the 0
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
；？sin；！sin?sinL？cos?cosLcos： (32) atat
where ? = solar declination [radians], L = local latitude [radians], and ： = the local hour angle of at
the sun [radians].
The local hour angle in radians is given by
trueSolarTime)??180！， (33) ：??4180??
trueSolarTime！localTime？eqtime，4?L，60?timezone (34) lm
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 is the local longitude lm
(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 9 December 16, 2008
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
f！t，t (35) sssr
where t = time of sunset [hours] and t = time of sunrise [hours]. sssr
Atmospheric attenuation. Various methods have been published to estimate the fraction of the atmospheric attenuation from a clear sky (a). Two alternative methods are available in QUAL2K t
to estimate a (Note that the solar radiation model is selected on the Light and Heat Worksheet t
1) Bras (default)
The Bras (1990) method computes a as: t
，namfac1 (36) a！et
where n is an atmospheric turbidity factor that varies from approximately 2 for clear skies to 4 fac
or 5 for smoggy urban areas. The molecular scattering coefficient (a) is calculated as 1
a！0.128，0.054logm (37) 110
where m is the optical air mass, calculated as
1m！ (38) ，1.253sin？0.15(？3.885)；；dowhere ； is the sun‟s altitude in degrees from the horizon = ； ? (180/)). d
2) Ryan and Stolzenbach
The Ryan and Stolzenbach (1972) model computes a from ground surface elevation and solar t
altitude as: 5.256288，0.0065elev??m??a！a288 (39) ??ttc
where a is the atmospheric transmission coefficient (0.70-0.91, typically approximately 0.8), and tc
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
2a！1，0.65C (40) cL
where C = fraction of the sky covered with clouds. L
Reflectivity. Reflectivity is calculated as
B (41) R！A；sd
where A and B are coefficients related to cloud cover (Table 3).
Table 3 Coefficients used to calculate reflectivity based on cloud cover. Cloudiness Clear Scattered Broken Overcast
C 0 0.1-0.5 0.5-0.9 1.0 L
A B A B A B A B Coefficients
1.18 ，0.77 2.20 ，0.97 0.95 ，0.75 0.35 ，0.45
QUAL2K 10 December 16, 2008