QUAL2K:

A Modeling Framework for Simulating River and Stream Water Quality

(Version 2.12)

Documentation

The Mystic River at Medford, MA

Steve Chapra, Greg Pelletier and Hua Tao

May 29, 2012
Chapra, S.C., Pelletier, G.J. and Tao, H. 2012. QUAL2K: A Modeling Framework for Simulating
River and Stream Water Quality, Version 2.12: Documentation and Users Manual. Civil and
Environmental Engineering Dept., Tufts University, Medford, MA, [email protected]

QUAL2K

May 29, 2012

Disclaimer
The information in this document has been funded partly by the United States Environmental
Protection Agency. It is currently being subjected to the Agency's peer and administrative review
and has yet to be approved for publication as an EPA document. Mention of trade names or
commercial products does not constitute endorsement or recommendation for use by the U.S.
Environmental Protection Agency.
The QUAL2K model (Q2K) described in this manual must be used at the user's own risk.
Neither the U.S. Environmental Protection Agency, Tufts University, the Washington Dept. of
Ecology, nor the program authors can assume responsibility for model operation, output,
interpretation or usage.
The creators of this program have used their best efforts in preparing this code. It is not
absolutely guaranteed to be error free. The author/programmer makes no warrantees, expressed or
implied, including without limitation warrantees of merchantability or fitness for any particular
purpose. No liability is accepted in any event for any damages, including accidental or
consequential damages, lost of profits, costs of lost data or programming materials, or otherwise
in connection with or arising out of the use of this program.

QUAL2K

May 29, 2012

1 INTRODUCTION
QUAL2K (or Q2K) is a river and stream water quality model that is intended to represent a
modernized version of the QUAL2E (or Q2E) model (Brown and Barnwell 1987). Q2K is similar
to Q2E in the following respects:

One dimensional. The channel is well-mixed vertically and laterally.

Branching. The system can consist of a mainstem river with branched tributaries.
Steady state hydraulics. Non-uniform, steady flow is simulated.
Diel heat budget. The heat budget and temperature are simulated as a function of meteorology
on a diel time scale.
Diel water-quality kinetics. All water quality variables are simulated on a diel time scale.
Heat and mass inputs. Point and non-point loads and withdrawals are simulated.

The QUAL2K framework includes the following new elements:

Software Environment and Interface. Q2K is implemented within the Microsoft Windows
environment. Numerical computations are programmed in Fortran 90. Excel is used as the
graphical user interface. All interface operations are programmed in the Microsoft Office
macro language: Visual Basic for Applications (VBA).
Model segmentation. Q2E segments the system into river reaches comprised of equally
spaced elements. Q2K also divides the system into reaches and elements. However, in
contrast to Q2E, the element size for Q2K can vary from reach to reach. In addition, multiple
loadings and withdrawals can be input to any element.
Carbonaceous BOD speciation. Q2K uses two forms of carbonaceous BOD to represent
organic carbon. These forms are a slowly oxidizing form (slow CBOD) and a rapidly
oxidizing form (fast CBOD).
Anoxia. Q2K accommodates anoxia by reducing oxidation reactions to zero at low oxygen
levels. In addition, denitrification is modeled as a first-order reaction that becomes
pronounced at low oxygen concentrations.
Sediment-water interactions. Sediment-water fluxes of dissolved oxygen and nutrients can be
simulated internally rather than being prescribed. That is, oxygen (SOD) and nutrient fluxes
are simulated as a function of settling particulate organic matter, reactions within the
sediments, and the concentrations of soluble forms in the overlying waters.
Bottom algae. The model explicitly simulates attached bottom algae. These algae have
variable stoichiometry.
Light extinction. Light extinction is calculated as a function of algae, detritus and inorganic
solids.
pH. Both alkalinity and total inorganic carbon are simulated. The rivers pH is then computed
based on these two quantities.
Pathogens. A generic pathogen is simulated. Pathogen removal is determined as a function of
temperature, light, and settling.
Reach specific kinetic parameters. Q2K allows you to specify many of the kinetic parameters
on a reach-specific basis.
Weirs and waterfalls. The hydraulics of weirs, as well as the effect of weirs and waterfalls on
gas transfer, are explicitly included.

QUAL2K

May 29, 2012

2 GETTING STARTED
As presently configured, an Excel workbook serves as the interface for QUAL2K. That is, all
input and output as well as model execution are implemented from within Excel. All interface
functions are programmed in Excels macro language: Visual Basic for Applications (VBA). All
numerical calculations are implemented in Fortran 90 for speed of execution. The following
material provides a step-by-step description of how the model can be set up on your computer and
used to perform a simulation.
Step 1: Copy the file, Q2Kv2_12.zip, to a directory (e.g., C:\). When this file is unzipped, it will
set up a subdirectory, Q2Kv2_12 which includes an Excel file (Q2KMasterv2_12b1.xls), and an
executable file (Q2KFortran2_11.exe). The first is the Q2K interface that allows you to run Q2K
and display its results. The second is the Fortran executable that actually performs the model
computations. These two files must always be in the same directory for the model to run
properly. Note that after you run the model, some assisting files will be automatically created by
the Fortran executable file to exchange information with Excel.
NOTE: DO NOT DELETE THE .zip file. If for some reason, you modify Q2K in a way that
makes it unusable, you can always use the zip file to reinstall the model.
Step 2: Create a subdirectory off of C:\Q2Kv2_12 called DataFiles.
Step 3: Open Excel and make sure that your macro security level is set to medium (Figure 1).
This can be done using the menu commands: Tools Macro Security. Make certain that the
Medium radio button is selected.

Figure 1 The Excel Macro Security Level dialogue box. In order to run Q2K, the Medium
level of security should be selected.

Step 4: Open Q2KMasterFortranv2_12.xls. When you do this, the Macro Security Dialogue Box
will be displayed (Figure 2).

QUAL2K

May 29, 2012

Figure 2 The Excel Macro security dialogue box. In order to run Q2K, the Enable Macros
button must be selected.

Click on the Enable Macros button.

Step 5: On the QUAL2K Worksheet, go to cell B10 and enter the path to the DataFiles
directory: C:\QUAL2K\DataFiles as shown in Figure 3.

Figure 3 The QUAL2K Worksheet showing the entry of the file path into cell B10.

Step 6: Click on the Run Fortran button.

If the program does not work correctly

There are two primary reasons why the program would not work properly. First, you may be
using an old version of Microsoft Office. Although Excel is downwardly compatible for some
earlier versions, Q2K will not work with very old versions.

QUAL2K

May 29, 2012

Second, you may have made a mistake in implementing the preceding steps. A common
mistake is to have mistyped the file path that you entered in cell B10. For example, suppose that
you mistyped the path as C:\Q2KFortranv2_12\DataFles. If this is the case, you will receive an
error message (Figure 4).

Figure 4 An error message that will occur if you type the incorrect file path into cell B10 on
the QUAL2K Worksheet.

If this occurs, click OK. This will terminate the run and bring you back to the QUAL2K
Worksheet where you can correct the file path entry.

If the program works correctly

QUAL2K will begin to execute. A window will open showing the progress of the Fortran
computations (Figure 5).

Figure 5 This window is displayed showing the progress of the model computations as
executed in Fortran. It allows you to follow the progress of a model run.

The program is set up to simulate a fictitious river with a mainstem along with two
tributaries. If the program works properly, the following dialogue box will appear when the run is
completed:

QUAL2K

May 29, 2012

Press OK and the following dialogue box will be displayed:

This box allows you to choose the parts of the system that you want to plot. As shown, it defaults
to the rivers Mainstem. Press OK to see the travel time for the Mainstem. Note that all plots are
updated when you press OK.
To switch to see the plots for one of the tributaries, you would press the button on the upper
left of the screen

This causes the plot options dialogue box to be displayed. The pulldown can then be used to
select another tributary.
Step 7: On the QUAL2K Worksheet click on the Open Old File button. Browse to get to the
directory: C:\Q2KFortranv2_12\DataFiles. You should see that a new file has been created with
the name that was specified in cell B9 (in the case of the example in Figure 3,
BogusExample.q2k). Click on the Cancel button to return to Q2K.
Note that every time that Q2K is run, a data file will be created with the file name specified in
cell B9 on the QUAL2K Worksheet (Figure 3). The program automatically affixes the extension
.q2k to the file name. Since this will overwrite previous versions of the file, make certain to
change the file name when you perform a new application.
Now that you have successfully run Q2K on your computer, the following pages are devoted
to documenting the science that underlies the model.

3 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 6, the
reaches are numbered in ascending order starting from the headwater of the rivers main stem.

QUAL2K

May 29, 2012

Notice that both point and non-point sources and point and non-point withdrawals (abstractions)
can be positioned anywhere along the channels length.
Headwater boundary
1
Point source

Point withdrawal

Point withdrawal

3
4
5

Point source

Non-point
withdrawal

6
Non-point
source

7
Point source

Downstream boundary

Figure 6 QUAL2K segmentation scheme for a river with no tributaries.

For systems with tributaries (Figure 7), 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
numbering continues at that tributarys headwater. Observe that both the headwaters and the
tributaries are also numbered consecutively following a sequencing scheme similar to the reaches.
Note also that the major branches of the system (that is, the main stem and each of the tributaries)
are referred to as segments. This distinction has practical importance because the software
provides plots of model output on a segment basis. That is, the software generates individual plots
for the main stem as well as each of the tributaries.
HW#1
1

HW#2
6
8

HW#3

2
ib
11
Tr
10
9

Tr
ib

2
1

13

4
5

12
14

15 16
17
18

HW#4

22

Tr
ib

23

19
3

20

24

21
25 26

Main stem

27
28
29
(a) A river with tributaries

QUAL2K

(b) Q2K reach representation

May 29, 2012

Figure 7 QUAL2K segmentation scheme for (a) a river with tributaries. The Q2K reach
representation in (b) illustrates the reach, headwater and tributary numbering
schemes.

Finally, any model reach can be further divided into a series of equally-spaced elements. As
in Figure 8, this is done by merely specifying the number of elements that are desired.

n=4

Elements

Reach

Figure 8 If desired, any model reach can be further subdivided into a series of n equallength elements.

In summary, the nomenclature used to describe the way in which Q2K organizes river
topology is as follows:

3.1

Reach. A length of river with constant hydraulic characteristics.

Element. The models 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.

Flow Balance

As described in the last section, Q2Ks most fundamental unit is the element. A steady-state flow
balance is implemented for each model element as (Figure 9)
Qi Qi 1 Qin,i Qout ,i Qevap ,i

(1)

where Qi = outflow from element i into the downstream element i + 1 [m3/d], Qi1 = inflow from
the upstream element i 1 [m3/d], Qin,i is the total inflow into the element from point and
nonpoint sources [m3/d], Qout,i is the total outflow from the element due to point and nonpoint
withdrawals [m3/d], and Qevap,i is the outflow due to evaporation [m3/d]. Thus, the downstream
outflow is simply the difference between inflow and source gains minus withdrawal and
evaporative losses.

QUAL2K

10

May 29, 2012

Qin,i

Qout,i

Qevap,i

Qi1

i1

Qi
i

i+1

Figure 9 Element flow balance.

The total inflow from sources is computed as

Qin,i

psi

Q ps ,i , j

j 1

npsi

Qnps,i, j

(2)

j 1

where Qps,i,j is the jth point source inflow to element i [m3/d], psi = the total number of point
sources to element i, Qnps,i,j is the jth non-point source inflow to element i [m3/d], and npsi = the
total number of non-point source inflows to element i.
The total outflow from withdrawals is computed as
pai

Qout,i

j 1

npai

Q pa ,i , j

(3)

npa ,i , j

j 1

where Qpa,i,j is the jth point withdrawal outflow from element i [m3/d], pai = the total number of
point withdrawals from element i, Qnpa,i,j is the jth non-point withdrawal outflow from element i
[m3/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 10, the nonpoint 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.
Qnpt

25%

25%

50%

start

end

Figure 10 The manner in which non-point source flow is distributed to an element.

QUAL2K

11

May 29, 2012

3.1.1

Evaporative Losses

Although the flow balance includes evaporative water losses, it should be noted that these losses
are not computed, but must be specified as an optional input by the user. This is done by entering
evaporative rates for each reach in mm/d on the Reach Worksheet. The rates are then converted
to m/d and multiplied by each elements surface area to generate an evaporative flow rate in m3/d
for use in each elements water balance.

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

3.2.1

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

3.2.1.1 Sharp-crested weirs

Figure 11 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 11 are defined as: Hi =
the depth of the element upstream of the weir [m], Hi+1 = the depth of the element downstream of
the weir [m], elev2i = the elevation above sea level of the tail end of the upstream element [m],
elev1i+1 = the elevation above sea level of the head end of the downstream element [m], Hw = the
height of the weir above elev2i [m], Hd = the drop between the elevation above sea level of the
surface of element i and element i+1 [m], Hh = the head above the weir [m], Bw = the width of the
weir [m]. Note that the width of the weir can differ from the width of the element, Bi.
(a) Side

(b) Cross-section

Bw
Hi

Hh

Hd

Hi

Hw

Hw

Hi+1
elev2i

elev2i
elev1i+1

elev1i+1

Figure 11 A sharp-crested weir occurring at the boundary between two reaches.

For the case where Hh/Hw < 0.4, flow is related to head by (Finnemore and Franzini 2002)
Qi 1.83Bw H h3/2

(4)

where Qi is the outflow from the element upstream of the weir (m3/s), and Bw and Hh are in m.
Equation (1) can be solved for
QUAL2K

12

May 29, 2012

 Qi
Hh

1.83Bw

2/3

(5)

This result can then be used to compute the depth of element i,

Hi H w H h

(6)

and the drop over the weir

H d elev 2i H i elev1i 1 H i 1

(7)

Note that this drop is used to compute oxygen and carbon dioxide gas transfer due to the weir (see
pages 43 and 48).
The cross-sectional area, velocity, surface area and volume of element i can then be
computed as
Ac,i Bi H i
Ui

(8)

Qi
Ac ,i

(9)

As ,i Bi xi
Vi Bi H i xi

where Bi = the width of element i, xi = the length of element i. Note that for reaches with weirs,
the reach width must be entered. This value is entered in the column AA (labeled "Bottom
Width") of the Reach Worksheet.
3.2.1.2 Broad-crested weirs
Using the same nomenclature as for the sharp-quested weir (see Figure 11), the equations for the
broad-crested weir are (Munson et al. 2009)
1.5

2
Qi Cw Bw g
3

H h1.5

where Cw = weir coefficient (dimensionless), and g = the gravitational constant (m/s2). Cw is

determined using the weir height (Hw) as in
Hh
Hw
Cw 1.125
H
2 h
Hw
1

Substituting Eq. (8) into Eq. (7) and collecting terms gives

QUAL2K

13

May 29, 2012

1.5

Hw Hh
2
Qi 1.125
Bw g
2H w H h
3

H h1.5

Equation (9) cannot be solved explicitly for Hh. However, it can be solved with a fixed-point or
successive substitution approach (Chapra 2011). This is done by rearranging it so that the last Hh
is brought to the left-hand side,
H hj

Q
i

1.1252/3 Bw g

1.5

2 H w H hj 1

H w H hj 1

2/3

where H hj and H hj 1 = the head above the top of the weir for the present and the previous
iterations, respectively. By guessing an initial head ( H h0 ), this equation can be solved iteratively to
determine H hj for j = 1, 2, ..., n.
3.2.2

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,
U aQ b

H Q

(10)

(11)

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
U
A
B c
H
Ac

(12)
(13)

The surface area and volume of the element can then be computed as
As Bx
V BH x

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
b

U = aQ
H = Q

QUAL2K

Typical
value

Exponent

0.43
0.45

14

Range

0.40.6
0.30.5

May 29, 2012

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

Manning Equation

Each element in a particular reach can be idealized as a trapezoidal channel (Figure 12). Under
conditions of steady flow, the Manning equation can be used to express the relationship between
flow and depth as
Q

S01/2 Ac5/3
n P 2/3

(14)

where Q = flow [m3/s] 1, S0 = bottom slope [m/m], n = the Manning roughness coefficient, Ac =
the cross-sectional area [m2], and P = the wetted perimeter [m].

S0

B1
H

1
ss1

1
ss2

B0

Q, U

Figure 12 Trapezoidal channel.

The cross-sectional area of a trapezoidal channel is computed as

Ac B0 0.5( ss1 ss 2 ) H H

(15)

where B0 = bottom width [m], ss1 and ss2 = the two side slopes as shown in Figure 12 [m/m], and
H = element depth [m].
The wetted perimeter is computed as
P B0 H ss21 1 H ss22 1

(16)

After substituting Eqs. (15) and (16), Eq. (14) can be solved iteratively for depth (Chapra and
Canale 2006),

Hk

(Qn)3/5 B0 H k 1 ss21 1 H k 1 ss22 1

S 3/10 B0 0.5( ss1 ss 2 ) H k 1

2/5

(17)

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

15

May 29, 2012

where k = 1, 2, , n, where n = the number of iterations. An initial guess of H0 = 0 is employed.

The method is terminated when the estimated error falls below a specified value of 0.001%. The
estimated error is calculated as
H k 1 H k
100%
H k 1

(18)

The cross-sectional area is determined with Eq. (15) and the velocity can then be determined
from the continuity equation,
U

Q
Ac

(19)

The average element width, B [m], is computed as

B

Ac
H

(20)

The top width, B1 [m], is computed as

B1 B0 ( s s1 s s 2 ) H
The surface area and volume of the element can then be computed as
As B1 x

V BHx
Suggested values for the Manning coefficient are listed in Table 2. Mannings 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 Mannings 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. For example, in upland streams, rather than the type of bed
material, the roughness is heavily influenced by the pool-riffle structure and can be very large
(Beven et al. 1979).
Table 2 The Manning roughness coefficient for various open channel surfaces (from Chow
et al. 1988).
MATERIAL
Man-made channels
Concrete
Gravel bottom with sides:
Concrete
mortared stone
Riprap
Natural stream channels
Clean, straight

QUAL2K

n
0.012
0.020
0.023
0.033
0.025-0.04

16

May 29, 2012

Clean, winding and some weeds

Weeds and pools, winding
Mountain streams with boulders
Heavy brush, timber
Steep upland channels

3.2.4

0.03-0.05
0.05
0.04-0.10
0.05-0.20
0.075-?

Waterfalls

In Section 3.2.1, 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 13). Note that waterfalls can only occur at the end of a reach.

Hi

Hd

elev2i

Hi+1
elev1i+1

Figure 13 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, elev2i > elev1i+1.

3.3

Travel Time

The residence time of each element is computed as

k

Vk
Qk

(21)

where k = the residence time of the kth element [d]; Vk = the volume of the kth element [m3] =
Ac,kxk; Ac,k = the cross-sectional area of the kth element [m2]; and xk = the length of the kth
element [m]. These times are then accumulated to determine the travel time along each of the
rivers 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,
tt, j

(22)

k 1

where tt,j = the travel time [d].

QUAL2K

17

May 29, 2012

3.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
channels hydraulics (Fischer et al. 1979),
E p ,i 0.011

U i2 Bi2

(23)

H iU i*

where Ep,i = the longitudinal dispersion between elements i and i + 1 [m2/s], Ui = velocity [m/s],
Bi = width [m], Hi = mean depth [m], and Ui* = shear velocity [m/s], which is related to more
fundamental characteristics by
U i* gH i Si

(24)

where g = acceleration due to gravity [= 9.81 m/s2] and S = channel slope [dimensionless].
After computing or prescribing En,i, the numerical dispersion is computed as
En , i

U i xi
2

(25)

The model dispersion Ei (i.e., the value used in the model calculations) is then computed as
follows:

If En,i Ep,i, the model dispersion, Ei is set to Ep,i En,i.

If En,i > Ep,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 14, 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
Qout ,i
dTi Qi 1
Q
E'
E'
Ti 1 i Ti
Ti i 1 Ti 1 Ti i Ti 1 Ti

dt
Vi
Vi
Vi
Vi
Vi
J a ,i
m3

wC pwVi 106 cm3 wC pw H i
Wh,i

J s ,i
m

C
100
cm

w pw H i

(26)

100 cm

where Ti = temperature in element i [oC], t = time [d], Ei = the bulk dispersion coefficient
between elements i and i + 1 [m3/d], Wh,i = the net heat load from point and non-point sources into

QUAL2K

18

May 29, 2012

element i [cal/d], w = the density of water [g/cm3], Cpw = the specific heat of water [cal/(g oC)],
Ja,i = the air-water heat flux [cal/(cm2 d)], and Js,i = the sediment-water heat flux [cal/(cm2 d)].

heat load

atmospheric
transfer

inflow

heat withdrawal

outflow

dispersion

dispersion
sediment-water
transfer
sediment

Figure 14 Heat balance for an element.

The bulk dispersion coefficient is computed as

Ei'

Ei Ac ,i

(27)

xi xi 1 / 2

Note that two types of boundary condition are used at the rivers 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)
npsi
psi

Wh,i C p Q ps ,i , j Tps ,i , j
Qnps ,i , j Tnps ,i , j
j 1

j 1

(28)

where Tps,i,j is the temperature of the jth point source for element i [oC], and Tnps,i,j is the
temperature of the jth non-point source temperature for element i [oC].

4.1

Surface Heat Flux

As depicted in Figure 15, surface heat exchange is modeled as a combination of five processes:
J h I (0) J an J br J c J e

(29)

where I(0) = net solar shortwave radiation at the water surface, Jan = net atmospheric longwave
radiation, Jbr = longwave back radiation from the water, Jc = conduction, and Je = evaporation.
All fluxes are expressed as cal/cm2/d.

QUAL2K

19

May 29, 2012

non-radiation terms

radiation terms

air-water
interface
solar
shortwave
radiation

atmospheric
longwave
radiation

water
longwave
radiation

net absorbed radiation

conduction
and
convection

evaporation
and
condensation

water-dependent terms

Figure 15 The components of surface heat exchange.

4.1.1

Solar Radiation

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

I (0)

I0

extraterrestrial atmospheric
radiation

ac

(1 Rs )

cloud

reflection

at

(1 S f )
shading

(30)

attenuation attenuation

where I(0) = solar radiation at the water surface [cal/cm2/d], I0 = extraterrestrial radiation (i.e., at
the top of the earths atmosphere) [cal/cm2/d], at = atmospheric attenuation, ac = cloud
attenuation, Rs = albedo (fraction reflected), and Sf = effective shade (fraction blocked by
vegetation and topography).
Extraterrestrial radiation. The extraterrestrial radiation is computed as (TVA 1972)
I0

W0
r2

sin

(31)

where W0 = the solar constant [1367 W/m2 or 2823 cal/cm2/d], r = normalized radius of the
earths orbit (i.e., the ratio of actual earth-sun distance to mean earth-sun distance), and = the
suns altitude [radians], which can be computed as
sin sin sin Lat cos cos Lat cos

(32)

where = solar declination [radians], Lat = local latitude [radians], and = the local hour angle
of the sun [radians].
The local hour angle in radians is given by
trueSolarTime

180
4

180

(33)

where:
QUAL2K

20

May 29, 2012

trueSolarTime localTime eqtime 4 Llm 60 timezone

(34)

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), Llm 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 NOAAs 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 tss tsr

(35)

where tss = time of sunset [hours] and tsr = time of sunrise [hours].
Atmospheric attenuation. Various methods have been published to estimate the fraction of the
atmospheric attenuation from a clear sky (at). Two alternative methods are available in QUAL2K
to estimate at (Note that the solar radiation model is selected on the Light and Heat Worksheet
of QUAL2K):
1) Bras (default)

The Bras (1990) method computes at as:

at e

n fac a1m

(36)

where nfac 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 (a1) is calculated as
a1 0.128 0.054 log10 m

(37)

where m is the optical air mass, calculated as

m

(38)

sin 0.15( d 3.885) 1.253

where d is the suns altitude in degrees from the horizon = (180o/).

2) Ryan and Stolzenbach

QUAL2K

21

May 29, 2012

The Ryan and Stolzenbach (1972) model computes at from ground surface elevation and solar
altitude as:
at atc

288 0.0065elev
m

288

5.256

(39)

where atc 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,
NOAAs 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 RyanStolzenbach 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
ac 1 0.65CL2

(40)

where CL = fraction of the sky covered with clouds.

Reflectivity. Reflectivity is calculated as
Rs A d B

(41)

where A and B are coefficients related to cloud cover (Table 3).

Table 3 Coefficients used to calculate reflectivity based on cloud cover.
Cloudiness
CL
Coefficients

Clear
0
A
1.18

B
0.77

Scattered
0.1-0.5
A
B
2.20
0.97

Broken
0.5-0.9
A
0.95

B
0.75

Overcast
1.0
A
B
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.
4.1.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
J an

Tair 273

sky

1 RL

(42)

where = the Stefan-Boltzmann constant = 11.7108 cal/(cm2 d K4), Tair = air temperature [oC],
sky = effective emissivity of the atmosphere [dimensionless], and RL = longwave reflection
QUAL2K

22

May 29, 2012

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)

Brunts (1932) equation is an empirical model that has been commonly used in water-quality
models (Thomann and Mueller 1987),
clear Aa Ab eair

where Aa and Ab are empirical coefficients. Values of Aa have been reported to range from about
0.5 to 0.7 and values of Ab have been reported to range from about 0.031 to 0.076 mmHg0.5 for
a wide range of atmospheric conditions. QUAL2K uses a default mid-range value of Aa = 0.6
together with a value of Ab = 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).
1/7

1.333224eair
clear 1.24

Ta

where eair is the air vapor pressure [mm Hg], and Ta 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):
17.27Td

eair 4.596e 237.3Td

(43)

where Td = the dew-point temperature [oC].

3) Koberg

Koberg (1964) reported that the Aa in Brunts formula depends on both air temperature and the
ratio of the incident solar radiation to the clear-sky radiation (Rsc). As in Figure 16, he presented a
series of curves indicating that Aa increases with Tair and decreases with Rsc with Ab held constant
at 0.0263 millibars0.5 (about 0.031 mmHg0.5).
The following polynomial is used in Q2K to provide a continuous approximation of Kobergs
curves.

QUAL2K

23

May 29, 2012

2
Aa ak Tair
bk Tair ck

where
3
2
ak 0.00076437Rsc
0.00121134Rsc
0.00073087Rsc 0.0001106
3
2
bk 0.12796842 Rsc
0.2204455Rsc
0.13397992 Rsc 0.02586655
3
2
ck 3.25272249Rsc
5.65909609Rsc
3.43402413Rsc 1.43052757

The fit of this polynomial to points sampled from Kobergs curves are depicted in Figure 16.
Note that an upper limit of 0.735 is prescribed for Aa.

Kobergs Aa constant for

Brunts equation for longwave radiation

0.8

Ratio of incident to clear sky radiation Rsc

0.50
0.65
0.75
0.80
0.85

0.7

0.6

0.90
0.95

0.5
1.00
0.4

0.3
-4

12

16

Air temperature Tair

20

24

28

32

(oC)

Figure 16 The points are sampled from Kobergs family of curves for determining the
value of the Aa constant in Brunts 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 Aa
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 (CL) of the form (TVA,
1972):
sky clear (1 0.17CL2 )

(44)

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

Water Long-wave Radiation

QUAL2K

24

May 29, 2012

The back radiation from the water surface is represented by the Stefan-Boltzmann law,
J br T 273

(45)

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

4.1.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,
J c c1 f U w Ts Tair

(46)

where c1 = Bowen's coefficient (= 0.47 mmHg/oC). The term, f(Uw), defines the dependence of
the transfer on wind velocity over the water surface where Uw 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(Uw):

QUAL2K

25

May 29, 2012

1) Brady, Graves, and Geyer (default)

f (U w ) 19.0 0.95U w2

where Uw = 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.
0.05
2
f (U w ) 0.271 (22.4 v1/3 ) 2 (24.2 Aacres
,iU w, mph )

where Uw,mph is wind speed in mph and Aacres,i is surface area of element i in acres. The constant
0.271 converts the original units of BTU ft2 day mmHg to cal cm day mmHg.
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
f (U w ) 0.271 (22.4 v1/3 ) 2 (17U w, mph ) 2

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 (v in F) accounts for the buoyancy of the moist air above a heated water surface. The
virtual temperature difference is estimated from water temperature (Tw,f in F), air temperature
(Tair,f in F), vapor pressure of water and air (es and eair in mmHg), and the atmospheric pressure
(patm is estimated as standard atmospheric pressure of 760 mmHg in QUAL2K):
Tair , f 460
Tw, f 460

v
460
460
1 0.378es / patm
1 0.378eair / patm

(47)

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

QUAL2K

26

May 29, 2012

(Uw,z in m/s) taken at any height (zw 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):
z
U w U wz
zw

0.15

(48)

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

Evaporation and Condensation

The heat loss due to evaporation can be represented by Daltons law,

J e f (U w )(es eair )

(49)

where es = the saturation vapor pressure at the water surface [mmHg], and eair = the air vapor
pressure [mmHg]. The saturation vapor pressure is computed as
17.27T

eair 4.596e 237.3T

4.2

(50)

Sediment-Water Heat Transfer

A heat balance for bottom sediment underlying a water element i can be written as
dTs ,i
dt

J s ,i

(51)

s C ps H sed ,i

where Ts,i = the temperature of the bottom sediment below element i [oC], Js,i = the sedimentwater heat flux [cal/(cm2 d)], s = the density of the sediments [g/cm3], Cps = the specific heat of
the sediments [cal/(g oC)], and Hsed,i = the effective thickness of the sediment layer [cm].
The flux from the sediments to the water can be computed as
J s ,i s C ps

s
H sed ,i / 2

Tsi Ti

86, 400 s
d

(52)

where s = the sediment thermal diffusivity [cm2/s].

The thermal properties of some natural sediments along with its components are summarized
in Table 4. Note that soft, gelatinous sediments found in the deposition zones of lakes are very
porous and approach the values for water. Some very slow, impounded rivers may approach such
a state. However, rivers tend to have coarser sediments with significant fractions of sands, gravels
and stones. Upland streams can have bottoms that are dominated by boulders and rock substrates.
Table 4 Thermal properties for natural sediments and the materials that comprise natural
sediments.

QUAL2K

27

May 29, 2012

Table 4. Thermal properties of various materials

Type of material
Sediment samples
Mud Flat
Sand
Mud Sand
Mud
Wet Sand
Sand 23% saturation with water
Wet Peat
Rock
Loam 75% saturation with water
Lake, gelatinous sediments
Concrete canal
Average of sediment samples:
Miscellaneous measurements:
Lake, shoreline
Lake soft sediments
Lake, with sand
River, sand bed
Component materials:
Water
Clay
Soil, dry
Sand
Soil, wet
Granite
Average of composite materials:

thermal conductivity

thermal diffusivity

w/m/C

cal/s/cm/C

m /s

cm /s

1.82
2.50
1.80
1.70
1.67
1.82
0.36
1.76
1.78
0.46
1.55
1.57

0.0044
0.0060
0.0043
0.0041
0.0040
0.0044
0.0009
0.0042
0.0043
0.0011
0.0037
0.0037

4.80E-07
7.90E-07
5.10E-07
4.50E-07
7.00E-07
1.26E-06
1.20E-07
1.18E-06
6.00E-07
2.00E-07
8.00E-07
6.45E-07

0.0048
0.0079
0.0051
0.0045
0.0070
0.0126
0.0012
0.0118
0.0060
0.0020
0.0080
0.0064

0.59

0.0014
3.25E-07
4.00E-07
7.70E-07

0.0033
0.0040
0.0077

1.40E-07
9.80E-07
3.70E-07
4.70E-07
4.50E-07
1.27E-06
6.13E-07

0.0014
0.0098
0.0037
0.0047
0.0045
0.0127
0.0061

0.59
1.30
1.09
0.59
1.80
2.89
1.37

0.0014
0.0031
0.0026
0.0014
0.0043
0.0069
0.0033

Cp

Cp

g/cm3

cal/(g C)

cal/(cm^3 C)

2.200

0.210

0.906
0.757
0.844
0.903
0.570
0.345
0.717
0.357
0.709
0.550
0.460
0.647

reference

(1)
"
"
"
(2)
(3)
(2)
(4)
(3)
(5)
"

(5)
"
"
"

1.000
1.490
1.500
1.520
1.810
2.700
1.670

0.999
0.210
0.465
0.190
0.525
0.202
0.432

1.000
0.310
0.700
0.290
0.950
0.540
0.632

(6)
"
"
"
"
"

(1) Andrews and Rodvey (1980)

(2) Geiger (1965)
(3) Nakshabandi and Kohnke (1965)
(4) Chow (1964) and Carslaw and Jaeger (1959)
(5) Hutchinson 1957, Jobson 1977, and Likens and Johnson 1969
(6) Cengel, Grigull, Mills, Bejan, Kreith and Bohn

Inspection of the component properties of Table 4 suggests that the presence of solid material
in stream sediments leads to a higher coefficient of thermal diffusivity than that for water or
porous lake sediments. In Q2K, we suggest a default value of 0.005 cm2/s for this quantity.
In addition, specific heat tends to decrease with density. Thus, the product of these two
quantities tends to be more constant than the multiplicands. Nevertheless, it appears that the
presence of solid material in stream sediments leads to a lower product than that for water or
gelatinous lake sediments. In Q2K, we suggest default values of s = 1.6 g/cm3 and Cps = 0.4
cal/(g oC). This corresponds to a product of 0.64 cal/(cm3 oC) for this quantity. Finally, as derived
in Appendix C, the sediment thickness is set by default to 10 cm in order to capture the effect of
the sediments on the diel heat budget for the water.

QUAL2K

28

May 29, 2012

5 CONSTITUENT MODEL
5.1

Constituents and General Mass Balance

The model constituents are listed in Table 5.

Table 5 Model state variables
Variable
Symbol
Units*
Conductivity
mhos
s
Inorganic suspended solids
mgD/L
mi
Dissolved oxygen
mgO2/L
o
Slowly reacting CBOD
mgO2/L
cs
Fast reacting CBOD
mgO2/L
cf
Organic nitrogen
gN/L
no
Ammonia nitrogen
gN/L
na
Nitrate nitrogen
gN/L
nn
Organic phosphorus
gP/L
po
Inorganic phosphorus
gP/L
pi
Phytoplankton
gA/L
ap
Phytoplankton nitrogen
gN/L
INp
Phytoplankton phosphorus
gP/L
IPp
Detritus
mgD/L
mo
Pathogen
cfu/100 mL
X
Alkalinity
mgCaCO3/L
Alk
Total inorganic carbon
mole/L
cT
Bottom algae biomass
mgA/m2
ab
Bottom algae nitrogen
mgN/m2
INb
Bottom algae phosphorus
mgP/m2
IPb
Constituent i
Constituent ii
Constituent iii
* mg/L g/m3; In addition, the terms D, C, N, P, and A refer to dry weight, carbon, nitrogen, phosphorus, and
chlorophyll a, respectively. The term cfu stands for colony forming unit which is a measure of viable bacterial
numbers.

For all but the bottom algae variables, a general mass balance for a constituent in an element
is written as (Figure 17)
Qout ,i
dci Qi 1
Q
E'
E'
W
ci 1 i ci
ci i 1 ci 1 ci i ci 1 ci i Si

dt
Vi
Vi
Vi
Vi
Vi
Vi

(53)

where Wi = the external loading of the constituent to element i [g/d or mg/d], and Si = sources and
sinks of the constituent due to reactions and mass transfer mechanisms [g/m3/d or mg/m3/d].

QUAL2K

29

May 29, 2012

mass load

atmospheric
transfer

inflow

mass withdrawal

outflow

dispersion

bottom algae

dispersion

sediments

Figure 17 Mass balance.

The external load is computed as (recall Eq. 2),

psi

npsi

j 1

j 1

Wi Q ps ,i , j c psi , j Qnps ,i , j c npsi , j

(54)

where cps,i,j is the jth point source concentration for element i [mg/L or g/L], and cnps,i,j is the jth
non-point source concentration for element i [mg/L or g/L].
For bottom algae, the transport and loading terms are omitted,
dab,i
dt

Sb ,i

dIN b,i
dt
dIPb,i
dt

(55)

SbN ,i

(56)

SbP ,i

(57)

where Sb,i = sources and sinks of bottom algae biomass due to reactions [mgA/m2/d], SbN,i =
sources and sinks of bottom algae nitrogen due to reactions [mgN/m2/d], and SbP,i = sources and
sinks of bottom algae phosphorus due to reactions [mgP/m2/d].
The sources and sinks for the state variables are depicted in Figure 18 (note that the internal
levels of nitrogen and phosphorus in the bottom algae are not depicted). The mathematical
representations of these processes are presented in the following sections.

QUAL2K

30

May 29, 2012

re

rda
ds

mo

rod

cs

ox

cT

cf

cT

ox

se

rna

no

cT
h

na

s
nn

se

sod

cf

mi

dn

Alk

se

cT
d

u
IN

rpa

po

pi

se

IP

r
e

cT

Figure 18 Model kinetics and mass transfer processes. The state variables are defined in
Table 5. Kinetic processes are dissolution (ds), hydrolysis (h), oxidation (ox),
nitrification (n), denitrification (dn), photosynthesis (p), respiration (r), excretion
(e), death (d), respiration/excretion (rx). Mass transfer processes are reaeration
(re), settling (s), sediment oxygen demand (SOD), sediment exchange (se), and
sediment inorganic carbon flux (cf).

5.2

Reaction Fundamentals

5.2.1

Biochemical Reactions

The following chemical equations are used to represent the major biochemical reactions that take
place in the model (Stumm and Morgan 1996):
Plant Photosynthesis and Respiration:

Ammonium as substrate:
106CO 2 16NH 4 HPO 24 106H 2 O

C106 H 263O110 N16 P1 106O 2 14H

(58)

Nitrate as substrate:
106CO 2 16NO3 HPO 24 122H 2 O+18H +

C106 H 263O110 N16 P1 138O 2

(59)

Nitrification:
NH 4 2O 2 NO3 H 2 O 2H

(60)

Denitrification:

QUAL2K

31

May 29, 2012

5CH 2 O 4NO3 4H 5CO 2 2N 2 7H 2 O

(61)

Note that a number of additional reactions are used in the model such as those involved with
simulating pH and unionized ammonia. These will be outlined when these topics are discussed
later in this document.
5.2.2

Stoichiometry of Organic Matter

The model requires that the stoichiometry of organic matter (i.e., phytoplankton and detritus) be
specified by the user. The following representation is suggested as a first approximation (Redfield
et al. 1963, Chapra 1997),
100 gD : 40 gC : 7200 mgN : 1000 mgP : 1000 mgA

(62)

where gX = mass of element X [g] and mgY = mass of element Y [mg]. The terms D, C, N, P,
and A refer to dry weight, carbon, nitrogen, phosphorus, and chlorophyll a, respectively. It should
be noted that chlorophyll a is the most variable of these quantities with a range of approximately
500-2000 mgA (Laws and Chalup 1990, Chapra 1997).
These values are then combined to determine stoichiometric ratios as in
rxy

gX
gY

(63)

For example, the amount of detritus (in grams dry weight or gD) that is released due to the death
of a unit amount of phytoplankton (in milligrams of chlorophyll a or mgA) can be computed as
rda

100 gD
gD
0.1
1000 mgA
mgA

5.2.2.1 Oxygen Generation and Consumption

The model requires that the rates of oxygen generation and consumption be prescribed. If
ammonia is the substrate, the following ratio (based on Eq. 58) can be used to determine the
grams of oxygen generated for each gram of plant matter that is produced through photosynthesis.
roca

106 moleO 2 (32 gO 2 /moleO 2 )

gO
2.67 2
106 moleC(12 gC/moleC)
gC

(64)

If nitrate is the substrate, the following ratio (based on Eq. 59) applies
rocn

138 moleO 2 (32 gO 2 /moleO 2 )

gO
3.47 2
106 moleC(12 gC/moleC)
gC

(65)

Note that Eq. (58) is also used for the stoichiometry of the amount of oxygen consumed for plant
respiration.
For nitrification, the following ratio is based on Eq. (60)

QUAL2K

32

May 29, 2012

ron

2 moleO 2 (32 gO 2 /moleO 2 )

gO
4.57 2
1 moleN(14 gN/moleN)
gN

(66)

5.2.2.2 CBOD Utilization Due to Denitrification

As represented by Eq. (61), CBOD is utilized during denitrification,

rondn 2.67

5.2.3

gO 2 5 moleC 12 gC/moleC
gO 2
1 gN

0.00286
gC 4 moleN 14 gN/moleN 1000 mgN
mgN

(67)

Temperature Effects on Reactions

The temperature effect for all first-order reactions used in the model is represented by
k (T ) k (20) T 20

(68)

where k(T) = the reaction rate [/d] at temperature T [oC] and = the temperature coefficient for
the reaction.

5.3

Composite Variables

In addition to the model's state variables, Q2K also displays several composite variables that are
computed as follows:
Total Organic Carbon (mgC/L):
TOC

cs c f
roc

rca a p rcd mo

(69)

Total Nitrogen (gN/L):

TN no na nn IN p

(70)

Total Phosphorus (gP/L):

TP po pi IPp

(71)

Total Kjeldahl Nitrogen (gN/L):

TKN no na IN p

(72)

Total Suspended Solids (mgD/L):

TSS rda a p mo mi

(73)

Ultimate Carbonaceous BOD (mgO2/L):

QUAL2K

33

May 29, 2012

CBODu cs c f roc rca a p roc rcd mo

5.4

(74)

Relationship of Model Variables and Data

For all but slow and fast CBOD (cf and cs), there exists a relatively straightforward relationship
between the model state variables and standard water-quality measurements. These are outlined
next. Then we discuss issues related to the more difficult problem of measuring CBOD.
5.4.1

Non-CBOD Variables and Data

The following are measurements that are needed for comparison of non-BOD variables with
model output:
TEMP = temperature (oC)
TKN = total kjeldahl nitrogen (gN/L) or TN = total nitrogen (gN/L)
NH4 = ammonium nitrogen (gN/L)
NO2 = nitrite nitrogen (gN/L)
NO3 = nitrate nitrogen (gN/L)
CHLA = chlorophyll a (gA/L)
TP = total phosphorus (gP/L)
SRP = soluble reactive phosphorus (gP/L)
TSS = total suspended solids (mgD/L)
VSS = volatile suspended solids (mgD/L)
TOC = total organic carbon (mgC/L)
DOC = dissolved organic carbon (mgC/L)
DO = dissolved oxygen (mgO2/L)
PH = pH
ALK = alkalinity (mgCaCO3/L)
COND = specific conductance (mhos)
The model state variables can then be related to these measurements as follows:
s = COND
mi = TSS VSS or TSS rdc (TOC DOC)
o = DO
no = TKN NH4 rna CHLA
or
no = TN NO2 NO3 NH4 rna CHLA
na = NH4
nn = NO2 + NO3
po = TP SRP rpa CHLA
pi = SRP
ap = CHLA
mo = VSS rda CHLA or rdc (TOC DOC) rda CHLA
pH = PH
Alk = ALK
5.4.2

Carbonaceous BOD

The interpretation of BOD measurements in natural waters is complicated by three primary

factors:

QUAL2K

34

May 29, 2012

Filtered versus unfiltered. If the sample is unfiltered, the BOD will reflect oxidation of
both dissolved and particulate organic carbon. Since Q2K distinguishes between
dissolved (cs and cf) and particulate (mo and ap) organics, an unfiltered measurement
alone does not provide an adequate basis for distinguishing these individual forms. In
addition, one component of the particulate BOD, phytoplankton (ap) can further
complicate the test through photosynthetic oxygen generation.

Nitrogenous BOD. Along with the oxidation of organic carbon (CBOD), nitrification also
contributes to oxygen depletion (NBOD). Thus, if the sample (a) contains reduced
nitrogen and (b) nitrification is not inhibited, the measurement includes both types of
BOD.

Incubation time. Short-term, usually 5-day, BODs are typically performed. Because Q2K
uses ultimate CBOD, 5-day BODs must be converted to ultimate BODs in a sensible
fashion.

We suggest the following as practical ways to measure CBOD in a manner that accounts for
the above factors and results in measurements that are compatible with Q2K.
Filtration. The sample should be filtered prior to incubation in order to separate dissolved from
particulate organic carbon.
Nitrification inhibition. Nitrification can be suppressed by adding a chemical inhibiting agent
such as TCMP (2-chloro-6-(trichloro methyl) pyridine. The measurement then truly reflects
CBOD. In the event that inhibition is not possible, the measured value can be corrected for
nitrogen by subtracting the oxygen equivalents of the reduced nitrogen (= ron TKN) in the
sample. However, as with all such difference-based adjustments, this correction may exhibit
substantial error.
Incubation time. The model is based on ultimate CBOD, so two approaches are possible: (1) use
a sufficiently long period so that the ultimate value is measured, or (2) use a 5-day measurement
and extrapolate the result to the ultimate. The latter method is often computed with the formula
CBODFNU

CBODFN5

(75)

1 e k1 5

where CBODFNU 2 = the ultimate dissolved carbonaceous BOD [mgO2/L], CBODFN5 = the 5day dissolved carbonaceous BOD [mgO2/L], and k1 = the CBOD decomposition rate in the bottle
[/d].
It should be noted that, besides practical considerations of time and expense, there may be
other benefits from using the 5-day measurement with extrapolation, rather that performing the
longer-term CBOD. Although extrapolation does introduce some error, the 5-day value has the
advantage that it would tend to minimize possible nitrification effects which, even when
inhibited, can begin to be exerted on longer time frames.
If all the above provisions are implemented, the result should correspond to the model
variables by
2

The nomenclature FNU stands for Filtration, Nitrification inhibition and Ultimate

QUAL2K

35

May 29, 2012

cf + cs = CBODFNU
Slow versus Fast CBOD. The final question relates to discrimination between fast and slow
CBODU. Although we believe that there is currently no single, simple, economically-feasible
answer to this problem, we think that the following 2 strategies represent the best current
alternatives.

Option 1: Represent all the dissolved, oxidizable organic carbon with a single pool (fast CBOD).
The model includes parameters to bypass slow CBOD. If no slow CBOD inputs are entered, this
effectively drops it from the model. For this case,
cf = CBODFNU
cs = 0

Option 2: Use an ultimate CBOD measurement for the fast fraction and compute slow CBOD by
difference with a DOC measurement. For this case,
cf = CBODFNU
cs = rocDOC CBODFNU

Option 2 works very nicely for systems where two distinct types of CBOD are present. For
example, sewage effluent and autochthonous carbon from the aquatic food chain might be
considered as fast CBOD. In contrast, industrial wastewaters such as pulp and paper mill effluent
or allochthonous DOC from the watershed might be considered more recalcitrant and hence could
be lumped into the slow CBOD fraction. In such case, the hydrolysis rate converting the slow into
the fast fraction could be set to zero to make the two forms independent.
For both options, the CBODFNU can either be (a) measured directly using a long incubation
time or (b) computed by extrapolation with Eq. 75. In both situations, a time frame of several
weeks to a month (i.e., a 20- to 30-day CBOD) is probably a valid period in order to oxidize most
of the readily degradable organic carbon. We base this assumption on the fact that bottle rates for
sewage-derived organic carbon are on the order of 0.05 to 0.3/d (Chapra 1997). As in Figure 19,
such rates suggest that much of the readily oxidizable CBOD will be exerted in about 20 to 30
days.
CBOD
CBODu

CBODu

0.3/d
0.5

0.1/d
0.05/d

0
0

10

20
30
time (days)

40

50

Figure 19 Progression of CBOD test for various levels of the bottle decomposition rate.

QUAL2K

36

May 29, 2012

In addition, we believe that practitioners should consider conducting long-term CBOD tests at
30oC rather than at the commonly employed temperature of 20oC. The choice of 20oC originated
from the fact that the average daily temperature of most receiving waters and wastewater
treatment plants in the temperate zone in summer is approximately 20oC.
If the CBOD measurement is intended to be used for regulation or to assess treatment plant
performance, it makes sense to standardize the test at a particular temperature. And for such
purposes, 20oC is as reasonable a choice as any. However, if the intent is to measure an ultimate
CBOD, anything that speeds up the process while not jeopardizing the measurement's integrity
would seem beneficial.
The saprophytic bacteria that break down nonliving organic carbon in natural waters and
sewage thrive best at temperatures from 20C to 40C. Thus, a temperature of 30oC is not high
enough that the bacterial assemblage would shift to thermophilic organisms that are atypical of
natural waters and sewage. The benefit should be higher oxidation rates which would result in
shorter analysis times for CBOD measurements. Assuming that a Q10 2 is a valid approximation
for bacterial decomposition, a 20-day BOD at 30oC should be equivalent to a 30-day BOD at
20oC.

5.5

Constituent Reactions

The mathematical relationships that describe the individual reactions and concentrations of the
model state variables (Table 5) are presented in the following paragraphs.
5.5.1

Conservative Substance (s)

By definition, conservative substances are not subject to reactions:

Ss = 0
5.5.2

(76)

Phytoplankton (ap)

Phytoplankton increase due to photosynthesis. They are lost via respiration, death, and settling
Sap PhytoPhoto PhytoResp PhytoDeath PhytoSettl

(77)

5.5.2.1 Photosynthesis

Phytoplankton photosynthesis is computed as

PhytoPhoto p a p

(78)

where p = phytoplankton photosynthesis rate [/d] is a function of temperature, nutrients, and

light,
p k gp (T )NpLp

(79)

where kgp(T) = the maximum photosynthesis rate at temperature T [/d], Np = phytoplankton

nutrient attenuation factor [dimensionless number between 0 and 1], and Lp = the phytoplankton
light attenuation coefficient [dimensionless number between 0 and 1].
QUAL2K

37

May 29, 2012

Nutrient Limitation. The nutrient limitation due to inorganic carbon is represented by a

Michaelis-Menten. In contrast, for nitrogen and phosphorus, the photosynthesis rate depends on
intracellular nutrient levels using a formulation originally developed by Droop (1974). The
minimum value is then employed to compute the nutrient attenuation factor,

q0 Np

qNp

Np min 1

,1

q0 Pp

[H 2 CO*3 ] [HCO3 ]

[H 2 CO*3 ] [HCO3 ]

qPp k sCp

(80)

where qNp and qPp = the phytoplankton cell quotas of nitrogen [mgN mgA1] and phosphorus
[mgP mgA1], respectively, q0Np and q0Pp = the minimum phytoplankton cell quotas of nitrogen
[mgN mgA1] and phosphorus [mgP mgA1], respectively, ksCp = inorganic carbon half-saturation
constant for phytoplankton [mole/L], [H2CO3*] = dissolved carbon dioxide concentration
[mole/L], and [HCO3] = bicarbonate concentration [mole/L]. The minimum cell quotas are the
levels of intracellular nutrient at which growth ceases. Note that the nutrient limitation terms
cannot be negative. That is, if q < q0, the limitation term is set to 0.
The cell quotas represent the ratios of the intracellular nutrient to the phytoplankton biomass,
qNp
qPp

IN p

(81)

ap
IPp

(82)

ap

where INp = phytoplankton intracellular nitrogen concentration [gN/L] and IPp = phytoplankton
intracellular phosphorus concentration [gP/L].
Light Limitation. It is assumed that light attenuation through the water follows the Beer-Lambert
law,
PAR ( z ) PAR (0)e ke z

(83)

where PAR(z) = photosynthetically active radiation (PAR) at depth z below the water surface
[ly/d] 3, and ke = the light extinction coefficient [m1]. The PAR at the water surface is assumed to
be a fixed fraction of the solar radiation (Szeicz 1984, Baker and Frouin 1987):
PAR(0) = 0.47 I(0)

The extinction coefficient is related to model variables by

ke keb i mi o mo p a p pn a 2/3
p

(84)

where keb = the background coefficient accounting for extinction due to water and color [/m], i,
o, p, and pn, are constants accounting for the impacts of inorganic suspended solids
3

ly/d = langley per day. A langley is equal to a calorie per square centimeter. Note that a ly/d is related to
the E/m2/s by the following approximation: 1 E/m2/s 0.45 Langley/day (LIC-OR, Lincoln, NE).
QUAL2K

38

May 29, 2012

[L/mgD/m], particulate organic matter [L/mgD/m], and chlorophyll [L/gA/m and (L/gA)2/3/m],
respectively. Suggested values for these coefficients are listed in Table 6.
Table 6 Suggested values for light extinction coefficients
Symbol

Value

Reference

i
o
p
pn

0.052
0.174
0.0088
0.054

Di Toro (1978)
Di Toro (1978)
Riley (1956)
Riley (1956)

Three models are used to characterize the impact of light on phytoplankton photosynthesis
(Figure 20):
1

Smith

Steele
Half-saturation

0.5

0
0

100

200

300

400

500

Figure 20 The three models used for phytoplankton and bottom algae photosynthetic light
dependence. The plot shows growth attenuation versus PAR intensity [ly/d].

Half-Saturation (Michaelis-Menten) Light Model (Baly 1935):

FLp

PAR( z )
K Lp PAR( z )

(85)

where FLp = phytoplankton growth attenuation due to light and KLp = the phytoplankton light
parameter. In the case of the half-saturation model, the light parameter is a half-saturation
coefficient [ly/d]. This function can be combined with the Beer-Lambert law and integrated over
water depth, H [m], to yield the phytoplankton light attenuation coefficient
Lp

K Lp PAR (0)
1
ln
ke H K Lp PAR(0)e ke H

(86)

Smiths Function (Smith 1936):

FLp

PAR( z )
2
K Lp

(87)

PAR ( z ) 2

where KLp = the Smith parameter for phytoplankton [ly/d]; that is, the PAR at which growth is
70.7% of the maximum. This function can be combined with the Beer-Lambert law and
integrated over water depth to yield
QUAL2K

39

May 29, 2012

Lp

PAR (0) / K Lp 1 PAR (0) / K Lp 2

1

ln
ke H
ke H PAR (0) / K
1 PAR(0) / K Lp e ke H
Lp e

(88)

Steeles Equation (Steele 1962):

FLp

PAR ( z ) 1

e
K Lp

PAR ( z )
K Lp

(89)

where KLp = the PAR at which phytoplankton growth is optimal [ly/d]. This function can be
combined with the Beer-Lambert law and integrated over water depth to yield
Lp

2.718282

e
ke H

PAR (0) ke H
e
K Lp

PAR (0)
K Lp

(90)

5.5.2.2 Losses
Respiration. Phytoplankton respiration is represented as a first-order rate that is attenuated at low
oxygen concentration,
PhytoResp Foxp krp (T ) a p

(91)

where krp(T) = temperature-dependent phytoplankton respiration/excretion rate [/d] and Foxp =

attenuation due to low oxygen [dimensionless]. Oxygen attenuation is modeled by Eqs. (127) to
(129) with the oxygen dependency represented by the parameter Ksop.
Death. Phytoplankton death is represented as a first-order rate,
PhytoDeath kdp (T )a p

(92)

where kdp(T) = temperature-dependent phytoplankton death rate [/d].

Settling. Phytoplankton settling is represented as
PhytoSettl

va
ap
H

(93)

where va = phytoplankton settling velocity [m/d].

5.5.3

Phytoplankton Internal Nitrogen (INb)

The change in intracellular nitrogen in phytoplankton cells is calculated from

S pN PhytoUpN qNp PhytoDeath PhytoExN

QUAL2K

40

(94)

May 29, 2012

where PhytoUpN = the uptake rate of nitrogen by phytoplankton (gN/L/d), PhytoDeath =

phytoplankton death (gN/L/d), and PhytoExN = the phytoplankton excretion of nitrogen
(gN/L/d), which is computed as
PhytoExN qNp kep (T )a p

(95)

where kep(T) = the temperature-dependent phytoplankton excretion rate [/d].

The N uptake rate depends on both external and intracellular nutrients as in (Rhee 1973),
PhytoUpN mNp

K qNp
na nn
ap
k sNp na nn K qNp (qNp q0 Np )

(96)

where mNp = the maximum uptake rate for nitrogen [mgN/mgA/d], ksNp = half-saturation constant
for external nitrogen [gN/L] and KqNp = half-saturation constant for intracellular nitrogen [mgN
mgA1].
5.5.4

Phytoplankton Internal Phosphorus (IPb)

The change in intracellular phosphorus in phytoplankton cells is calculated from

S pP PhytoUpP qPp PhytoDeath PhytoExP

(97)

where PhytoUpP = the uptake rate of phosphorus by phytoplankton (gP/L/d), PhytoDeath =

phytoplankton death (gP/L/d), and PhytoExP = the phytoplankton excretion of phosphorus
(gP/L/d), which is computed as
PhytoExP qPp kep (T ) a p

(98)

where kep(T) = the temperature-dependent phytoplankton excretion rate [/d].

The P uptake rate depends on both external and intracellular nutrients as in (Rhee 1973),
PhytoUpP mPp

K qPp
pi
ap
k sPp pi K qPp (qPp q0 Pp )

(99)

where mPp = the maximum uptake rate for phosphorus [mgP/mgA/d], ksPp = half-saturation
constant for external phosphorus [gP/L] and KqPp = half-saturation constant for intracellular
phosphorus [mgP mgA1].
5.5.5

Bottom algae (ab)

Bottom algae increase due to photosynthesis. They are lost via respiration and death.
Sab BotAlgPhoto BotAlgResp BotAlgDeath

(100)

5.5.5.1 Photosynthesis

QUAL2K

41

May 29, 2012

Two representations can be used to model bottom algae photosynthesis. The first is based on a
temperature-corrected zero-order rate attenuated by nutrient and light limitation (McIntyre 1973,
Rutherford et al. 1999),
BotAlgPhoto C gb (T ) NbLb

(101)

where Cgb(T) = the zero-order temperature-dependent maximum photosynthesis rate [mgA/(m2

d)], Nb = bottom algae nutrient attenuation factor [dimensionless number between 0 and 1], and
Lb = the bottom algae light attenuation coefficient [dimensionless number between 0 and 1].
The second uses a first-order model,
BotAlgPhoto Cgb (T ) NbLbSb ab

(102)

where, for this case, Cgb(T) = the first-order temperature-dependent maximum photosynthesis rate
[d1], and Sb = bottom algae space limitation attenuation factor.
Temperature Effect. As for the first-order rates, an Arrhenius model is employed to quantify the
effect of temperature on bottom algae photosynthesis,
Cgb (T ) Cgb (20) T 20

(103)

Nutrient Limitation. The effect of nutrient limitation on bottom plant photosynthesis is modeled
in the same way as for the phytoplankton. That is, a Droop (1974) formulation is used for
nitrogen and phosphorus limitation whereas a Michaelis-Menten equation is employed for
inorganic carbon,

Nb min 1

q0 Nb
q
[H 2 CO*3 ] [HCO3 ]
,1 0 Pb ,

qNb
qPb k sCb [H 2 CO*3 ] [HCO3 ]

(104)

where qNb and qPb = the bottom algae cell quotas of nitrogen [mgN mgA1] and phosphorus [mgP
mgA1], respectively, q0Nb and q0Pb = the bottom algae minimum cell quotas of nitrogen [mgN
mgA1] and phosphorus [mgP mgA1], respectively, and ksCb = the bottom algae inorganic carbon
half-saturation constant [mole/L]. As was the case for phytoplankton, the nutrient limitation terms
cannot be negative.
The cell quotas represent the ratios of the intracellular nutrient to the bottom plants biomass,
qNb

INb
ab

(105)

qPb

IPb
ab

(106)

where INb = intracellular nitrogen concentration [mgN/m2] and IPb = intracellular phosphorus
concentration [mgP/m2].

QUAL2K

42

May 29, 2012

Light Limitation. In contrast to the phytoplankton, light limitation at any time is determined by
the amount of PAR reaching the bottom of the water column. This quantity is computed with the
Beer-Lambert law (recall Eq. 83) evaluated at the bottom of the river,
PAR ( H ) PAR(0)e ke H

(107)

As with the phytoplankton, three models (Eqs. 85, 87 and 89) are used to characterize the
impact of light on bottom algae photosynthesis. Substituting Eq. (107) into these models yields
the following formulas for the bottom algae light attenuation coefficient,
Half-Saturation Light Model (Baly 1935):
Lb

PAR (0)e ke H

(108)

K Lb PAR (0)e ke H

Smiths Function: (Smith 1936)

Lb

PAR (0)e ke H
2
K Lb

PAR (0)e

(109)

ke H 2

Steeles Equation (Steele 1962):

Lb

PAR(0)e ke H 1

e
K Lb

PAR (0) e ke H
K Lb

(110)

where KLb = the appropriate bottom algae light parameter for each light model.
Space Limitation. If a first-order growth model is used, a term must be included to impose a
space limitation on the bottom algae. A logistic model is used for this purpose as in
Sb 1

ab
ab,max

where ab,max = the carrying capacity [mgA/m2].

5.5.5.2 Losses
Respiration. Bottom algae respiration is represented as a first-order rate that is attenuated at low
oxygen concentration,
BotAlgResp Foxb krb (T ) ab

(111)

where krb(T) = temperature-dependent bottom algae respiration rate [/d] and Foxb = attenuation
due to low oxygen [dimensionless]. Oxygen attenuation is modeled by Eqs. (127) to (129) with
the oxygen dependency represented by the parameter Ksob.
Death. Bottom algae death is represented as a first-order rate,

QUAL2K

43

May 29, 2012

BotAlgDeath kdb (T )ab

(112)

where kdb(T) = the temperature-dependent bottom algae death rate [/d].

5.5.6

Bottom Algae Internal Nitrogen (INb)

The change in intracellular nitrogen in bottom algal cells is calculated from

SbN BotAlgUpN q Nb BotAlgDeath BotAlgExN

(113)

where BotAlgUpN = the uptake rate of nitrogen by bottom algae (mgN/m2/d), BotAlgDeath =
bottom algae death (mgA/m2/d), and BotAlgExN = the bottom algae excretion of nitrogen
(mgN/m2/d), which is computed as
BotAlgExN qNb keb (T )ab

(114)

where keb(T) = the temperature-dependent bottom algae excretion rate [/d].

The N uptake rate depends on both external and intracellular nutrients as in (Rhee 1973),
BotAlgUpN mNb

K qNb
na nn
ab
ksNb na nn K qNb (qNb q0 Nb )

(115)

where mNb = the maximum uptake rate for nitrogen [mgN/mgA/d], ksNb = half-saturation constant
for external nitrogen [gN/L] and KqNb = half-saturation constant for intracellular nitrogen
[mgN/mgA].
5.5.7

Bottom Algae Internal Phosphorus (IPb)

The change in intracellular phosphorus in bottom algal cells is calculated from

SbP BotAlgUpP qPb BotAlgDeath BotAlgExP

(116)

where BotAlgUpP = the uptake rate of phosphorus by bottom algae (mgP/m2/d), BotAlgDeath =
bottom algae death (mgA/m2/d), and BotAlgExP = the bottom algae excretion of phosphorus
(mgP/m2/d), which is computed as
BotAlgExP qPb keb (T )ab

(117)

where keb(T) = the temperature-dependent bottom algae excretion rate [/d].

The P uptake rate depends on both external and intracellular nutrients as in (Rhee 1973),
BotAlgUpP mPb

QUAL2K

K qPb
pi
ab
ksPb pi K qPb (qPb q0 Pb )

44

(118)

May 29, 2012

where mPb = the maximum uptake rate for phosphorus [mgP/mgA/d], ksPb = half-saturation
constant for external phosphorus [gP/L] and KqPb = half-saturation constant for intracellular
phosphorus [mgP mgA1].
5.5.8

Detritus (mo)

Detritus or particulate organic matter (POM) increases due to plant death. It is lost via dissolution
and settling
Smo rda PhytoDeath rda

BotAlgDeath
DetrDiss DetrSettl
H

(119)

where
DetrDiss kdt (T )mo

(120)

where kdt(T) = the temperature-dependent detritus dissolution rate [/d] and

DetrSettl

vdt
mo
H

(121)

where vdt = detritus settling velocity [m/d].

5.5.9

Slowly Reacting CBOD (cs)

Slowly reacting CBOD increases due to detritus dissolution. It is lost via hydrolysis and
oxidation,
Scs (1 F f )rod DetrDiss SlowCHydr SlowCOxid

(122)

where Ff = the fraction of detrital dissolution that goes to fast reacting CBOD [dimensionless],
and
SlowCHydr khc (T )cs

(123)

where khc(T) = the temperature-dependent slow CBOD hydrolysis rate [/d], and
SlowCOxid Foxc kdcs (T )cs

(124)

where kdcs(T) = the temperature-dependent slow CBOD oxidation rate [/d] and Foxc = attenuation
due to low oxygen [dimensionless].
5.5.10 Fast Reacting CBOD (cf)

Fast reacting CBOD is gained via the dissolution of detritus and the hydrolysis of slowly-reacting
CBOD. It is lost via oxidation and denitrification.
Scf F f rod DetrDiss SlowCHydr FastCOxid rondn Denitr

QUAL2K

45

(125)

May 29, 2012

where
FastCOxid Foxc kdc (T )c f

(126)

where kdc(T) = the temperature-dependent fast CBOD oxidation rate [/d] and Foxc = attenuation
due to low oxygen [dimensionless]. The parameter rondn is the ratio of oxygen equivalents lost per
nitrate nitrogen that is denitrified (Eq. 67). The term Denitr is the rate of denitrification
[gN/L/d]. It will be defined in Sec. 5.5.15 below.
Three formulations are used to represent the oxygen attenuation:
Half-Saturation:
Foxrp

o
K socf o

(127)

where Ksocf = half-saturation constant for the effect of oxygen on fast CBOD oxidation [mgO2/L].
Exponential:
Foxrp (1 e

K socf o

(128)

where Ksocf = exponential coefficient for the effect of oxygen on fast CBOD oxidation [L/mgO2].
Second-Order Half Saturation:
Foxrp

o2

(129)

K socf o 2

where Ksocf = half-saturation constant for second-order effect of oxygen on fast CBOD oxidation
[mgO22/L2].
5.5.11 Organic Nitrogen (no)

Organic nitrogen increases due to plant death. It is lost via hydrolysis and settling.
Sno f onp qNp PhytoDeath f onb q Nb

BotAlgDeath
ONHydr ONSettl
H

(130)

where fonp = the fraction of the phytoplankton internal nitrogen that is in organic form which is
calculated as
f onp
f onp 1

rna
qNp

if qNp rna

(131)

if qNp rna

The fraction of the bottom algae internal nitrogen that is in organic form, fonb, is calculated in a
similar fashion.

QUAL2K

46

May 29, 2012

The rate of organic nitrogen hydrolysis is computed as

ONHydr khn (T )no

(132)

where khn(T) = the temperature-dependent organic nitrogen hydrolysis rate [/d]. Organic nitrogen
settling is determined as
ONSettl

von
no
H

(133)

where von = organic nitrogen settling velocity [m/d].

5.5.12 Ammonia Nitrogen (na)

Ammonia nitrogen increases due to organic nitrogen hydrolysis and plant death and excretion. It
is lost via nitrification and plant photosynthesis:
BotAlgDeath
H
(134)
BotAlgExN
BotAlgUpN
PhytoExN
Nitrif Pap PhytoUpN Pab
NH3GasLoss
H
H

Sna ONHydr (1 f onp )q Np PhytoDeath (1 f onb ) qNb

The ammonia nitrification rate is computed as

Nitrif Foxna kn (T )na

(135)

where kn(T) = the temperature-dependent nitrification rate for ammonia nitrogen [/d] and Foxna =
attenuation due to low oxygen [dimensionless]. Oxygen attenuation is modeled by Eqs. (127) to
(129) with the oxygen dependency represented by the parameter Ksona.
The coefficients Pap and Pab are the preferences for ammonium as a nitrogen source for
phytoplankton and bottom algae, respectively,
Pap

Pab

na khnxp
na nn

na )(khnxp nn ) (na nn )(khnxp nn )

(136)

na nn
na khnxb

(khnxb na )(khnxb nn ) (na nn )(khnxb nn )

(137)

(khnxp

where khnxp = preference coefficient of phytoplankton for ammonium [mgN/m3] and khnxb =
preference coefficient of bottom algae for ammonium [mgN/m3].
5.5.13 Unionized Ammonia

The model simulates total ammonia. In water, the total ammonia consists of two forms:
ammonium ion, NH4+, and unionized ammonia, NH3. At normal pH (6 to 8), most of the total
ammonia will be in the ionic form. However at high pH, unionized ammonia predominates. The
amount of unionized ammonia can be computed as

QUAL2K

47

May 29, 2012

nau Fu na

(138)

where nau = the concentration of unionized ammonia [gN/L], and Fu = the fraction of the total
ammonia that is in unionized form,
Ka

Fu

10

pH

(139)

Ka

where Ka = the equilibrium coefficient for the ammonia dissociation reaction, which is related to
temperature by
pK a 0.09018

2729.92
Ta

(140)

where Ta is absolute temperature [K] and pKa = log10(Ka). Note that the fraction of the total
ammonia that is in ionized form, Fi, can be computed as 1 Fu or
Fi

10 pH
10

pH

(141)

Ka

5.5.14 Ammonia Gas Transfer

The loss of ammonia via gas transfer is computed as

NH3GasLoss

vnh3 (T )
naus (T ) nau
H

where vnh3(T) = the temperature-dependent ammonia gas-transfer coefficient [m/d], and naus(T) =
the saturation concentration of ammonia [gN/L] at temperature, T.
The transfer coefficient is calculated as
vnh3 Kl

He
K
H e RTa l
Kg

where vv = the mass-transfer coefficient (m/d), Kl and Kg = liquid and gas film exchange
coefficients [m/d], respectively, R = the universal gas constant (= 8.206105 atm m3/(K mole)),
Ta = absolute temperature [K], and He = Henrys constant (atm m3/mole).
The saturation concentration is calculated as
naus (T )

pnh3
CF
He

where pnh3 = the partial pressure of ammonia in the atmosphere (atm), and CF is a conversion
factor (gN/L per moleNH3/m3). The partial pressure of ammonia ranges from 1-10 ppb in rural

QUAL2K

48

May 29, 2012

and moderately polluted areas to 10-100 ppb in heavily polluted areas (Holland 1978, FinlaysonPitts and Pitts 1986). We will assume that a value of 2 ppb, which corresponds to 2109 atm,
represents a typical value. The conversion factor is
CF

m3
14gN
106 gN
gN/L

14 103
1000L moleNH3
gN
moleNH3 / m3

The liquid-film coefficient can be related to the oxygen reaeration rate by (Mills et al. 1982),
32
Kl ka H
17

0.25

1.171ka H

The gas-film coefficient is computed by

18
K g K g,H2O
17

0.25

Kg,H2O = a mass-transfer velocity for water vapor [m/d], which can be related to wind speed by
(Schwarzenbach et al. 1993)
K g ,H2O 0.2U w,10 0.3

m
86,400s

100cm
d

where Uw,10 = wind speed at a height of 10 m (m/s). Combining the above equations gives
K g 175.287U w,10 262.9305

The Henrys constant for ammonia at 20oC is 1.3678105 atm m3/mole (Kavanaugh and
Trussell 1980). The value at temperatures other than 20oC can be computed with
H e (T ) H e (20)1.052T 20

5.5.15 Nitrate Nitrogen (nn)

Nitrate nitrogen increases due to nitrification of ammonia. It is lost via denitrification and plant
uptake:
Sni Nitrif Denitr (1 Pab )

BotAlgUptakeN
H

(142)

The denitrification rate is computed as

Denitr (1 Foxdn )kdn (T )nn

(143)

where kdn(T) = the temperature-dependent denitrification rate of nitrate nitrogen [/d] and Foxdn =
effect of low oxygen on denitrification [dimensionless] as modeled by Eqs. (127) to (129) with
the oxygen dependency represented by the parameter Ksodn.

QUAL2K

49

May 29, 2012

5.5.16 Organic Phosphorus (po)

Organic phosphorus increases due to plant death and excretion. It is lost via hydrolysis and
settling.
PhytoDeath
BotAlgDeath
f opb qPb
OPHydr OPSettl
H
H

S po f opp qPp

(144)

where fopp = the fraction of the phytoplankton internal phosphorus that is in organic form which is
calculated as
f opp

rpa

if qPp rpa

qPp

f opp 1

(145)

if qPp rpa

The fraction of the bottom algae internal phosphorus that is in organic form, fopb, is calculated in a
similar fashion.
The rate of organic phosphorus hydrolysis is computed as
OPHydr khp (T ) po

(146)

where khp(T) = the temperature-dependent organic phosphorus hydrolysis rate [/d]. Organic
phosphorus settling is determined as
OPSettl

vop
H

po

(147)

where vop = organic phosphorus settling velocity [m/d].

5.5.17 Inorganic Phosphorus (pi)

Inorganic phosphorus increases due to organic phosphorus hydrolysis and plant excretion. It is
lost via plant uptake. In addition, a settling loss is included for cases in which inorganic
phosphorus is lost due to sorption onto settleable particulate matter such as iron oxyhydroxides:
BotAlgDeath
H
BotAlgExP
BotAlgUpP
PhytoExP
PhytoUpP
IPSettl
H
H

S pi OPHydr (1 f opp ) qPp PhytoDeath (1 f opb ) qPb

(148)

where
IPSettl

vip
H

pi

(149)

where vip = inorganic phosphorus mass transfer coefficient [m/d].

5.5.18 Inorganic Suspended Solids (mi)

QUAL2K

50

May 29, 2012

Inorganic suspended solids are lost via settling,

Smi = InorgSettl
where
InorgSettl

vi
mi
H

(150)

where vi = inorganic suspended solids settling velocity [m/d].

5.5.19 Dissolved Oxygen (o)

Dissolved oxygen increases due to plant photosynthesis. It is lost via fast CBOD oxidation,
nitrification and plant respiration. Depending on whether the water is undersaturated or
oversaturated it is gained or lost via reaeration,
So roa PhytoPhoto ro a

BotAlgPhoto
roc FastCOxid ron NH4Nitr
H
BotAlgResp
roa PhytoResp roa
OxReaer
H

(151)

where
OxReaer ka (T ) os (T , elev) o

(152)

where ka(T) = the temperature-dependent oxygen reaeration coefficient [/d], os(T, elev) = the
saturation concentration of oxygen [mgO2/L] at temperature, T, and elevation above sea level,
elev.
5.5.19.1 Oxygen Saturation

The following equation is used to represent the dependence of oxygen saturation on temperature
(APHA 1995)
ln os (T , 0) 139.34411

1.575701105 6.642308 107

Ta
Ta2

1.243800 1010
Ta3

8.621949 1011

(153)

Ta4

where os(T, 0) = the saturation concentration of dissolved oxygen in freshwater at 1 atm [mgO2/L]
and Ta = absolute temperature [K] where Ta = T +273.15.
The effect of elevation is accounted for by
os (T , elev) eln os (T , 0) (1 0.0001148elev)

QUAL2K

(154)

51

May 29, 2012

5.5.19.2 Reaeration Formulas

The reaeration coefficient (at 20 oC) can be prescribed on the Reach Worksheet. If reaeration is
not prescribed (that is, it is blank or zero for a particular reach), it is computed as a function of the
rivers hydraulics and (optionally) wind velocity,
ka (20) kah (20)

K Lw (20)
H

(155)

where kah(20) = the reaeration rate at 20oC computed based on the rivers hydraulic characteristics
[d1], KL(20) = the reaeration mass-transfer coefficient based on wind velocity [m/d], and H =
mean depth [m].
Hydraulic-based Formulas:

OConnor-Dobbins (OConnor and Dobbins 1958):

kah (20) 3.93

U 0.5

(156)

H 1.5

where U = mean water velocity [m/s] and H = mean water depth [m].
Churchill (Churchill et al. 1962):
kah (20) 5.026

(157)

H 1.67

Owens-Gibbs (Owens et al. 1964):

kah (20) 5.32

U 0.67

(158)

H 1.85

Tsivoglou-Neal (Tsivoglou and Neal 1976):

Low flow, Q = 0.0283 to 0.4247 cms (1 to 15 cfs):
kah (20) 31,183 US

(159)

High flow, Q = 0.4247 to 84.938 cms (15 to 3000 cfs):

kah (20) 15,308 US

(160)

where S = channel slope [m/m].

Thackston-Dawson (Thackston and Dawson 2001):
kah (20) 2.16(1 9 F 0.25 )

QUAL2K

U*
H

(161)

52

May 29, 2012

where U* = shear velocity [m/s], and F = the Froude number [dimensionless]. The shear velocity
and the Froude number are defined as
U * gRh S

(162)

and
U

(163)

gH d

where g = gravitational acceleration (= 9.81 m/s2), Rh = hydraulic radius [m], S = channel slope
[m/m], and Hd = the hydraulic depth [m]. The hydraulic depth is defined as
Hd

Ac
Bt

(164)

where Bt = the top width of the channel [m].

USGS (Pool-riffle) (Melching and Flores 1999):
Low flow, Q < 0.556 cms (< 19.64 cfs):
kah (20) 517(US )0.524 Q 0.242

(165)

High flow, Q > 0.556 cms (> 19.64 cfs):

kah (20) 596(US )0.528 Q 0.136

(166)

where Q = flow (cms).

USGS (Channel-control) (Melching and Flores 1999):
Low flow, Q < 0.556 cms (< 19.64 cfs):
kah (20) 88(US )0.313 H 0.353

(167)

High flow, Q > 0.556 cms (> 19.64 cfs):

kah (20) 142(US )0.333 H 0.66 Bt0.243

(168)

Internal (Covar 1976):

Reaeration can also be internally calculated based on the following scheme patterned after a plot
developed by Covar (1976) (Figure 21):

If H < 0.61 m, use the Owens-Gibbs formula

If H > 0.61 m and H > 3.45U2.5, use the OConnor-Dobbins formula
Otherwise, use the Churchill formula

QUAL2K

53

May 29, 2012

This is referred to as option Internal on the Rates Worksheet of Q2K. Note that if no option is
specified, the Internal option is the default.

10

OConnor
Dobbins
0.05
0.1

Churchill

Depth (m)

0.2
0.5
1

2
5
10
20
50
100

0.1
0.1

Owens
Gibbs

Velocity (mps)
Figure 21 Reaeration rate (/d) versus depth and velocity (Covar 1976).
Wind-based Formulas:

Three options are available to incorporate wind effects into the reaeration rate: (1) it can be
omitted, (2) the Banks-Herrera formula and (3) the Wanninkhof formula.
Banks-Herrera formula (Banks 1975, Banks and Herrera 1977):
Klw 0.728U w0.5,10 0.317U w,10 0.0372U w2 ,10

(169)

where Uw,10 = wind speed measured 10 meters above the water surface [m s1]
Wanninkhof formula (Wanninkhof 1991):

K lw 0.0986U w1.,64
10

(170)

Note that the model uses Eq. (48) to correct the wind velocity entered on the Meteorology
Worksheet (7 meters above the surface) so that it is scaled to the 10-m height.
5.5.19.3 Effect of Control Structures: Oxygen

QUAL2K

54

May 29, 2012

Oxygen transfer in streams is influenced by the presence of control structures such as weirs,
dams, locks, and waterfalls (Figure 22). Butts and Evans (1983) have reviewed efforts to
characterize this transfer and have suggested the following formula,
rd 1 0.38ad bd H d (1 0.11H d )(1 0.046T )

(171)

where rd = the ratio of the deficit above and below the dam, Hd = the difference in water elevation
[m] as calculated with Eq. (7), T = water temperature (C) and ad and bd are coefficients that
correct for water-quality and dam-type. Values of ad and bd are summarized in Table 7. If no
values are specified, Q2K uses the following default values for these coefficients: ad = 1.25 and
bd = 0.9.

Qi1 oi1
Qi1 oi1

Hd
i1

Figure 22 Water flowing over a river control structure.

Table 7 Coefficient values used to predict the effect of dams on stream reaeration.

(a) Water-quality coefficient

ad

Polluted state

0.65
1.0
1.6
1.8

Gross
Moderate
Slight
Clean
(b) Dam-type coefficient

bd

Dam type

Flat broad-crested regular step

Flat broad-crested irregular step
Flat broad-crested vertical face
Flat broad-crested straight-slope face
Flat broad-crested curved face
Round broad-crested curved face
Sharp-crested straight-slope face
Sharp-crested vertical face
Sluice gates

0.70
0.80
0.60
0.75
0.45
0.75
1.00
0.80
0.05

The oxygen mass balance for the element below the structure is written as
Qab,i
Wo,i
doi Qi 1
Q
E'
o 'i 1 i oi
oi i oi 1 oi

S o ,i
dt
Vi
Vi
Vi
Vi
Vi

(172)

where oi1 = the oxygen concentration entering the element [mgO2/L], where
QUAL2K

55

May 29, 2012

o 'i 1 os ,i 1

os ,i 1 oi 1

(173)

rd

5.5.20 Pathogen (X)

Pathogens are subject to death and settling,

S X PathDeath PathSettl

(174)

5.5.20.1 Death

Pathogen death is due to natural die-off and light (Chapra 1997). The death of pathogens in the
absence of light is modeled as a first-order temperature-dependent decay and the death rate due to
light is based on the Beer-Lambert law,
PathDeath kdX (T ) X path

I (0) / 24
1 e ke H X
ke H

(175)

where kdX(T) = temperature-dependent pathogen die-off rate [/d] and path = a light efficiency
factor [dimensionless].
5.5.20.2 Settling

Pathogen settling is represented as

PathSettl

vX
x
H

(176)

where vX = pathogen settling velocity [m/d].

5.5.21 pH

The following equilibrium, mass balance and electroneutrality equations define a freshwater
dominated by inorganic carbon (Stumm and Morgan 1996),
K1
K2

[HCO3 ][H ]

(177)

[H 2 CO*3 ]
[CO32 ][H ]

(178)

[HCO3 ]

K w [H ][OH ]

(179)

cT [H 2 CO*3 ] [HCO3 ] [CO32 ]

Alk [HCO3 ] 2[CO32 ] [OH ] [H ]

(180)
(181)

where K1, K2 and Kw are acidity constants, Alk = alkalinity [eq L1], H2CO3* = the sum of
dissolved carbon dioxide and carbonic acid, HCO3 = bicarbonate ion, CO32 = carbonate ion, H+
= hydronium ion, OH = hydroxyl ion, and cT = total inorganic carbon concentration [mole L1].
The brackets [ ] designate molar concentrations.

QUAL2K

56

May 29, 2012

Note that the alkalinity is expressed in units of eq/L for the internal calculations. For input
and output, it is expressed as mgCaCO3/L. The two units are related by
Alk (mgCaCO3 /L) 50, 000 Alk (eq/L)

(182)

The equilibrium constants are corrected for temperature by

Harned and Hamer (1933):
pK w =

4787.3
7.1321 log10 (Ta ) 0.010365Ta 22.80
Ta

(183)

Plummer and Busenberg (1982):

log K1 = 356.3094 0.06091964Ta 21834.37 / Ta
126.8339 log Ta 1, 684,915 / Ta2

(184)

Plummer and Busenberg (1982):

log K 2 = 107.8871 0.03252849Ta 5151.79 / Ta
38.92561log Ta 563, 713.9 / Ta2

(185)

The nonlinear system of five simultaneous equations (177 through 181) can be solved
numerically for the five unknowns: [H2CO3*], [HCO3], [CO32], [OH], and {H+}. An efficient
solution method can be derived by combining Eqs. (177), (178) and (180) to define the quantities
(Stumm and Morgan 1996)
[H ]2

(186)

[H ] K1[H ]+K1 K 2
K1[H ]

(187)

K1 K 2

(188)

[H ]2 K1[H ]+K1 K 2

[H ]2 K1[H ]+K1 K 2

where 0, 1, and 2 = the fraction of total inorganic carbon in carbon dioxide, bicarbonate, and
carbonate, respectively. Equations (179), (187), and (188) can be substituted into Eq. (181) to
yield,
Alk=(1 2 2 )cT

Kw

[H ]

[H ]

(189)

Thus, solving for pH reduces to determining the root, {H+}, of

f ([H ])=(1 2 2 )cT

QUAL2K

Kw

[H ]

[H ] Alk

(190)

57

May 29, 2012

where pH is then calculated with

pH log10 [H ]

(191)

The root of Eq. (190) is determined with a numerical method. The user can choose either
bisection, Newton-Raphson or Brents method (Chapra and Canale 2006, Chapra 2007) as
specified on the QUAL2K sheet. The Newton-Raphson is the fastest but can sometimes diverge.
In contrast, the bisection method is slower, but more reliable. Because it balances speed with
reliability, Brents method is the default.
5.5.22 Total Inorganic Carbon (cT)

Total inorganic carbon concentration increases due to fast carbon oxidation and plant respiration.
It is lost via plant photosynthesis. Depending on whether the water is undersaturated or
oversaturated with CO2, it is gained or lost via reaeration,
BotAlgResp
H
BotAlgPhoto
rcca PhytoPhoto rcca
CO2Reaer
H

ScT rcco FastCOxid rcca PhytoResp rcca

(192)

where
CO2Reaer kac (T ) [CO 2 ]s 0 cT

(193)

where kac(T) = the temperature-dependent carbon dioxide reaeration coefficient [/d], and [CO2]s =
the saturation concentration of carbon dioxide [mole/L].
The stoichiometric coefficients are computed as 4
gC moleC
m3
rcca rca

mgA 12 gC 1000 L
rcco

1
roc

gC

gO 2

(194)

moleC
m3

12 gC 1000 L

(195)

5.5.22.1 Carbon Dioxide Saturation

The CO2 saturation is computed with Henrys law,

[CO 2 ]s K H pCO2

(196)

where KH = Henry's constant [mole (L atm)1] and pCO2 = the partial pressure of carbon dioxide in
the atmosphere [atm]. Note that the partial pressure is input on the Rates Worksheet in units of
ppm. The program internally converts ppm to atm using the conversion: 106 atm/ppm.
4

The conversion, m3 = 1000 L is included because all mass balances express volume in m3, whereas total
inorganic carbon is expressed as mole/L.

QUAL2K

58

May 29, 2012

The value of KH can be computed as a function of temperature by (Edmond and Gieskes

1970)
pK H =

2385.73
0.0152642Ta 14.0184
Ta

(197)

The partial pressure of CO2 in the atmosphere has been increasing, largely due to the combustion
of fossil fuels (Figure 23). Values in 2011 are approximately 103.4073 atm (= 391.5 ppm).

concentration (ppm)

400
380
360
340
320
300
1950 1960 1970 1980 1990 2000 2010 2020

Figure 23 Concentration of carbon dioxide in the atmosphere as recorded at Mauna Loa

Observatory, Hawaii. 5

5.5.22.2 CO2 Gas Transfer

The CO2 reaeration coefficient can be computed from the oxygen reaeration rate by
32
kac (20)
44

0.25

0.923 ka (20)

(198)

5.5.22.3 Effect of Control Structures: CO2

As was the case for dissolved oxygen, carbon dioxide gas transfer in streams can be influenced by
the presence of control structures. Q2K assumes that carbon dioxide behaves similarly to
dissolved oxygen (recall Sec. 5.5.19.3). Thus, the inorganic carbon mass balance for the element
immediately downstream of the structure is written as
dcT ,i
dt

Qab,i
WcT ,i
Qi 1
Q
E'
c 'T ,i 1 i cT ,i
cT ,i i cT ,i 1 cT ,i
ScT ,i
Vi
Vi
Vi
Vi
Vi

(199)

where c'T,i1 = the concentration of inorganic carbon entering the element [mgO2/L], where
c 'T ,i 1 (1 2 )cT ,i 1 CO 2, s ,i 1

CO 2, s ,i 1 2 cT ,i 1
rd

(200)

http://www.esrl.noaa.gov/gmd/ccgg/trends/

QUAL2K

59

May 29, 2012

where rd is calculated with Eq. (171).

5.5.23 Alkalinity (Alk)

As summarized in the present model accounts for changes in alkalinity due to several
mechanisms:
Table 8 Processes that effect alkalinity.
Process
Nitrif
Denitr
OPHydr
ONHydr
PhytoPhoto

Utilize
NH4
NO3

SRP
NH4
NH4
NO3
SRP

PhytoResp
PhytoUpN
PhytoUpP
PhytoExcrN
PhytoExcrP
BotAlgUpN
BotAlgUpP
BotAlgExcrN
BotAlgExcrP

Create
NO3

NH4
SRP
NH4
NO3
SRP
NH4
SRP
NH4
NO3
SRP
NH4
SRP

Alkalinity change
Decrease
Increase
Decrease
Increase
Decrease
Increase
Increase
Increase
Decrease
Decrease
Increase
Increase
Increase
Decrease
Decrease
Increase
Increase
Increase
Decrease

5.5.23.1 Nitrification

According to Eq. (60), nitrification utilizes ammonium and creates nitrate. Hence, because a
positive ion is taken up and a negative ion is created, the alkalinity is decreased by two
equivalents. The change in alkalinity can be related to the nitrification rate by
Sa , nitr

2 eq
moleN
gN 50,000 mgCaCO3
gN
Nitrif

6
moleN 14.007 gN 10 gN
1 eq
Ld

(201)

5.5.23.2 Denitrification

According to Eq. (61), denitrification utilizes nitrate and creates nitrogen gas. Hence, because a
negative ion is taken up and a neutral compound is created, the alkalinity is increased by one
equivalent. The change in alkalinity can be related to the denitrification rate by
Sa , denitr

1 eq
moleN
gN 50,000 mgCaCO3
gN
Denitr

6
moleN 14.007 gN 10 gN
1 eq
Ld

(202)

5.5.23.3 Organic P Hydrolysis

Hydrolysis of organic P results in the creation of inorganic phosphate. Depending on the pH, the
phosphate will either have 1 (pH 2 to 7) or 2 (pH 7 to 12) negative charges. Hence, because

QUAL2K

60

May 29, 2012

negative ions are being created, the alkalinity is decreased by one or two equivalents,
respectively. The change in alkalinity can be related to the P hydrolysis rate by 6,
Sa ,OPh
( H 2 PO 4 2 HPO 4 3 PO 4 )

eq
moleP
gP 50,000 mgCaCO3
gP
OPHydr

moleP 30.974 gP 106 gP

1 eq
Ld

(203)

where
H 2 PO 4
HPO 4

PO 4

K p1[H ]2

(204)

[H ]3 K p1[H ]2 +K p1 K p 2 [H ]+K p1 K p 2 K p 3
K p1 K p 2 [H ]

(205)

[H ]3 K p1[H ]2 +K p1 K p 2 [H ]+K p1 K p 2 K p 3
K p1 K p 2 K p 3
3

(206)

[H ] K p1[H ] +K p1 K p 2 [H ]+K p1 K p 2 K p 3

where Kp1 = 102.15, Kp2 = 107.2, and Kp3 = 1012.35.

5.5.23.4 Organic N Hydrolysis

Hydrolysis of organic N results in the creation of ammonia. Depending on the pH, the ammonia
will either be in the form of ammonium ion with a single positive charge (pH < 9) or neutral
ammonia gas (pH > 9). Hence, when the positive ions are created, the alkalinity is increased by
one equivalent. The change in alkalinity can be related to the N hydrolysis rate by
Sa ,ONh Fi

1 eq
moleN
gN 50,000 mgCaCO3
gN
ONHydr

6
moleN 14.007 gN 10 gN
1 eq
Ld

(207)

5.5.23.5 Phytoplankton Photosynthesis

Phytoplankton photosynthesis takes up nitrogen as either ammonia or nitrate and phosphorus as

inorganic phosphate. If ammonia is the primary nitrogen source, this leads to a decrease in
alkalinity because the uptake of the positively charged ammonium ions is much greater than the
uptake of the negatively charged phosphate ions (recall the stoichiometry as described on p. 32).
If nitrate is the primary nitrogen source, this leads to an increase in alkalinity because both nitrate
and phosphate are negatively charged.
The following representation relates the change in alkalinity to phytoplankton photosynthesis
depending on the nutrient sources as well as their speciation as governed by the pH,

Note that although it will almost always be negligible, Eq. (203) includes PO43 for completeness.

QUAL2K

61

May 29, 2012

Sa , PhytP

50,000 mgCaCO3
1 eq

1 eq
moleN
gN
rna Pap Fi
6
moleN
14.007
gN
10 gN

1 eq
moleN
gN
rna 1 Pap
6
moleN 14.007 gN 10 gN

(208)

rpa ( H 2 PO 4 2 HPO 4 3 PO 4 )

1 eq moleP
gP

6
moleP 30.974 gP 10 gP
gA
PhytoPhoto

Ld

5.5.23.6 Phytoplankton Nutrient Uptake

Phytoplankton take up nitrogen as either ammonia or nitrate and phosphorus as inorganic

phosphate. The following representation relates the change in alkalinity to phytoplankton uptake
rates depending on the nutrient sources as well as their speciation as governed by the pH,
S a , PUp

50,000 mgCaCO3
1 eq

1 eq
moleN
gN
PhytoUpN(gN/L/d)
Pap Fi (1 Pap )

6
moleN
14.007
gN
H (m)
10
gN

( H 2 PO 4 2 HPO 4 3 PO 4 )

(209)

1 eq moleP
gP
PhytoUpP( gP/L/d)

moleP 30.974 gP 106 gP

H (m)

5.5.23.7 Phytoplankton Nutrient Excretion

Phytoplankton excrete ammonia and inorganic phosphate. The following representation relates
the change in alkalinity to phytoplankton excretion rates including the effect of pH on the
nutrients speciation,
S a , PEx

50,000 mgCaCO3
1 eq

1 eq
moleN
gN
PhytoExN(gN/L/d)
Fi

6
moleN
14.007
gN
H (m)
10
gN

( H 2 PO 4 2 HPO 4 3 PO 4 )

(210)

1 eq moleP
gP
PhytoExP(gP/L/d)

6
moleP 30.974 gP 10 gP
H (m)

5.5.23.8 Bottom Algae Nutrient Uptake

Bottom algae take up nitrogen as either ammonia or nitrate and phosphorus as inorganic
phosphate. The following representation relates the change in alkalinity to bottom algae uptake
rates depending on the nutrient sources as well as their speciation as governed by the pH,

QUAL2K

62

May 29, 2012

S a ,BAUp

50,000 mgCaCO3
1 eq

1 eq
moleN
gN
BotAlgUpN(mgN/m 2 /d)
Pab Fi (1 Pab )

H (m)
moleN 14.007 gN 10 gN

( H 2 PO 4 2 HPO 4 3 PO 4 )

(211)

1 eq moleP
gP
BotAlgUpP(mgP/m 2 /d)

H (m)
moleP 30.974 gP 10 gP

5.5.23.9 Bottom Algae Nutrient Excretion

Bottom algae excrete ammonia and inorganic phosphate. The following representation relates the
change in alkalinity to bottom algae excretion rates including the effect of pH on the nutrients
speciation,
S a , BAEx

50,000 mgCaCO3
1 eq

1 eq
moleN
gN
BotAlgExN(mgN/m 2 /d)
Fi
moleN 14.007 gN 106 gN
H (m)

( H 2 PO 4 2 HPO 4 3 PO 4 )

5.6

(212)

1 eq moleP
gP
BotAlgExP(mgP/m 2 /d)

moleP 30.974 gP 10 gP
H (m)

SOD/Nutrient Flux Model

Sediment nutrient fluxes and sediment oxygen demand (SOD) are based on a model developed by
Di Toro (Di Toro et al. 1991, Di Toro and Fitzpatrick. 1993, Di Toro 2001). The present version
also benefited from James Martins (Mississippi State University, personal communication)
efforts to incorporate the Di Toro approach into EPAs WASP modeling framework.
A schematic of the model is depicted in Figure 24. As can be seen, the approach allows
oxygen and nutrient sediment-water fluxes to be computed based on the downward flux of
particulate organic matter from the overlying water. The sediments are divided into 2 layers: a
thin ( 1 mm) surface aerobic layer underlain by a thicker (10 cm) lower anaerobic layer. Organic
carbon, nitrogen and phosphorus are delivered to the anaerobic sediments via the settling of
particulate organic matter (i.e., phytoplankton and detritus). There they are transformed by
mineralization reactions into dissolved methane, ammonium and inorganic phosphorus. These
constituents are then transported to the aerobic layer where some of the methane and ammonium
are oxidized. The flux of oxygen from the water required for these oxidations is the sediment
oxygen demand. The following sections provide details on how the model computes this SOD
along with the sediment-water fluxes of carbon, nitrogen and phosphorus that are also generated
in the process.

QUAL2K

63

May 29, 2012

cf o

AEROBIC

WATER

Jpom

CH4

CO2

nn

pi

NH4p

NH4d

NO3

N2

NH4p

NH4d

NO3

N2

PO4p

PO4d

PO4p

PO4d

CH4(gas)

POC
ANAEROBIC

na

PON

POP

DIAGENESIS

METHANE

AMMONIUM

NITRATE

PHOSPHATE

Figure 24 Schematic of SOD-nutrient flux model of the sediments.

5.6.1

Diagenesis

As summarized in Figure 25, the first step in the computation involves determining how much of
the downward flux of particulate organic matter (POM) is converted into soluble reactive forms
in the anaerobic sediments. This process is referred to as diagenesis. First the total downward flux
of carbon, nitrogen and prosphorus are computed as the sums of the fluxes of settling
phytoplankton and organic matter delivered to the sediments from the water column
J POC rca va a p rcd vdt mo
J PON va q Np a p von no

(213)

J POC va qPp a p vop po

where JPOC = the downward flux of POC [gC m2 d1], rca = the ratio of carbon to chlorophyll a
[gC/mgA], va = phytoplankton settling velocity [m/d], ap = phytoplankton concentration
[mgA/m3], rcd = the ratio of carbon to dry weight [gC/gD], vdt = detritus settling velocity [m/d],
mo = detritus concentration [gD/m3], JPON = the downward flux of PON [mgN m2 d1], qNp = the
phytoplankton nitrogen cell quota [mgN mgA1], von = organic nitrogen settling velocity [m/d], no
= organic nitrogen concentration [mgN m3], JPOP = the downward flux of POP [mgP m2 d1], qPp
= the phytoplankton phosphorus cell quota [mgP mgA1], vop = organic phosphorus settling
velocity [m/d], and po = organic phosphorus concentration [mgP m3] .

QUAL2K

64

May 29, 2012

vaap

vdtmo

JPOC, G1

POC2,G1

JPOC, G2

POC2,G2

JPOC, G3

POC2,G3

JPON, G1

PON2,G1

JPON, G2

PON2,G2

JPON, G3

PON2,G3

JPOP, G1

POP2,G1

JPOP, G2

POP2,G2

JPOP, G3

POP2,G3

fG1

rda

JPOC

fG2
fG3

roc

JPOM
rcd

fG1
rnd

JPON

fG2
fG3

rpd
fG1
JPOP

fG2
fG3

v dtmo

rcd
roc

vaap

JPOC

fC2
fC3

rca

vonno

fN1
JPON

fN2
fN3

vaqNpap

voppo

POC2,G1

JPOC, G2

POC2,G2

JPOC, G3

POC2,G3

JPON, G1

PON2,G1

JPON, G2

PON2,G2

JPON, G3

PON2,G3

JPOP, G1

POP2,G1

JPOP, G2

POP2,G2

JPOP, G3

POP2,G3

fP1
JPOP

v aqPpap

JPOC, G1
fC1

fP2
fP3

JC, G1
JC

JC, G2

JN, G1
JN

JN, G2

JP, G1
JP

JP, G2

JC, G1
JC, G2

JC

JN, G1
JN, G2

JN

JP, G1
JP, G2

JP

Figure 25 Representation of how settling particulate matter is transformed into fluxes of

dissolved carbon (JC) , nitrogen (JN) and phosphorus (JP) in the anaerobic
sediments.

Note that for convenience, we express the particulate organic carbon (POC) as oxygen
equivalents using the stoichiometric coefficient roc. Each of the nutrient fluxes is further broken
down into three reactive fractions: labile (G1), slowly reacting (G2) and non-reacting (G3).
These fluxes are then entered into mass balances to compute the concentration of each
fraction in the anaerobic layer. For example, for labile POC, a mass balance is written as
H2

dPOC2,G1
dt

T 20
J POC ,G1 k POC ,G1 POC
,G1 H 2 POC2,G1 w2 POC2,G1

(214)

where H2 = the thickness of the anaerobic layer [m], POC2,G1 = the concentration of the labile
fraction of POC in the anaerobic layer [gO2/m3], JPOC,G1 = the flux of labile POC delivered to the
anaerobic layer [gO2/m2/d], kPOC,G1 = the mineralization rate of labile POC [d1], POC,G1 =
temperature correction factor for labile POC mineralization [dimensionless], and w2 = the burial
velocity [m/d]. At steady state, Eq. (214) can be solved for

QUAL2K

65

May 29, 2012

POC2,G1

J POC ,G1

(215)

T 20
k POC ,G1 POC
,G1 H 2 vb

The flux of labile dissolved carbon, JC,G1 [gO2/m2/d], can then be computed as
T 20
J C ,G1 k POC ,G1 POC
,G1 H 2 POC2,G1

(216)

In a similar fashion, a mass balance can be written and solved for the slowly reacting
dissolved organic carbon. This result is then added to Eq. (216) to arrive at the total flux of
dissolved carbon generated in the anaerobic sediments.
J C J C ,G1 J C ,G 2

(217)

Similar equations are developed to compute the diagenesis fluxes of nitrogen, JN [gN/m2/d], and
phosphorus JP [gP/m2/d].
5.6.2

Ammonium

Based on the mechanisms depicted in Figure 24, mass balances can be written for total
ammonium in the aerobic layer and the anaerobic layers,
H1

dNH 4,1
dt

12 f pa 2 NH 4, 2 f pa1 NH 4,1 K L12 f da 2 NH 4, 2 f da1 NH 4,1 w2 NH 4,1

(218)

K NH 4
o
n
NH 4,1 T 20
s a f da1 NH 4,1
f da1 NH 4,1
NH 4
1000
s
K

NH
2
K

NH 4
4,1
NH 4 ,O 2 o
2

H2

dNH 4, 2
dt

J N 12 f pa1 NH 4,1 f pa 2 NH 4, 2 K L12 f da1 NH 4,1 f da 2 NH 4, 2

(219)

w2 NH 4,1 NH 4, 2

where H1 = the thickness of the aerobic layer [m], NH4,1 and NH4,2 = the concentration of total
ammonium in the aerobic layer and the anaerobic layers, respectively [gN/m3], na = the
ammonium concentration in the overlying water [mgN/m3], NH4,1 = the reaction velocity for
nitrification in the aerobic sediments [m/d], NH4 = temperature correction factor for nitrification
[dimensionless], KNH4 = ammonium half-saturation constant [gN/m3], o = the dissolved oxygen
concentration in the overlying water [gO2/m3], and KNH4,O2 = oxygen half-saturation constant
[mgO2/L], and JN = the diagenesis flux of ammonium [gN/m2/d].
The fraction of ammonium in dissolved (fdai) and particulate (fpai) form are computed as
1
1 mi ai
1 f dai

f dai

(220)

f pai

(221)

where mi = the solids concentration in layer i [gD/m3], and ai = the partition coefficient for
ammonium in layer i [m3/gD].

QUAL2K

66

May 29, 2012

The mass transfer coefficient for particle mixing due to bioturbation between the layers, 12
[m/d], is computed as
12

T 20
D p Dp
POC2,G1 / roc

H2

(222)

K M , Dp o

POCR

where Dp = diffusion coefficient for bioturbation [m2/d], Dp = temperature coefficient

[dimensionless], POCR = reference G1 concentration for bioturbation [gC/m3] and KM,Dp = oxygen
half-saturation constant for bioturbation [gO2/m3].
The mass transfer coefficient for pore water diffusion between the layers, KL12 [m/d], is
computed as,
K L12

T 20
Dd Dd
H2 / 2

(223)

where Dd = pore water diffusion coefficient [m2/d], and Dd = temperature coefficient

[dimensionless].
The mass transfer coefficient between the water and the aerobic sediments, s [m/d], is
computed as
s

SOD
o

(224)

where SOD = the sediment oxygen demand [gO2/m2/d].

At steady state, Eqs. (218) and (219) are two simultaneous nonlinear algebraic equations. The
equations can be linearized by assuming that the NH4,1 term in the Monod term for nitrification is
constant. The simultaneous linear equations can then be solved for NH4,1 and NH4,2. The flux of
ammonium to the overlying water can then be computed as
n

J NH 4 s f da1 NH 4,1 a
1000

5.6.3

(225)

Nitrate

Mass balances for nitrate can be written for the aerobic and anaerobic layers as
H1

dNO3,1
dt

K L12 NO3,2 NO3,1 w2 NO3,1 s n NO3,1

1000

H2

dNO3,2

QUAL2K

dt

2
NH
4,1

T 20
NH
4

NO3,1 T 20
K NH 4
o
NO3 NO3,1
f da1 NH 4,1
K NH 4 NH 4,1 2 K NH 4,O 2 o
s

T 20
J N K L12 NO3,1 NO3,2 w2 NO3,1 NO3,2 NO3,2 NO
3 NO3,2

67

(226)

(227)

May 29, 2012

where NO3,1 and NO3,2 = the concentration of nitrate in the aerobic layer and the anaerobic layers,
respectively [gN/m3], nn = the nitrate concentration in the overlying water [mgN/m3], NO3,1 and
NO3,2 = the reaction velocities for denitrification in the aerobic and anaerobic sediments,
respectively [m/d], and NO3 = temperature correction factor for denitrification [dimensionless].
In the same fashion as for Eqs. (218) and (219), Eqs. (226) and (227) can be linearized and
solved for NO3,1 and NO3,2. The flux of nitrate to the overlying water can then be computed as
n

J NO3 s NO3,1 n
1000

(228)

Denitrification requires a carbon source as represented by the following chemical equation,

5CH 2 O 4NO3 4H 5CO 2 2N 2 7H 2 O

(229)

The carbon requirement (expressed in oxygen equivalents per nitrogen) can therefore be
computed as
rondn 2.67

gO 2 5 moleC 12 gC/moleC
gO 2
1 gN

0.00286
gC 4 moleN 14 gN/moleN 1000 mgN
mgN

(230)

Therefore, the oxygen equivalents consumed during denitrification, JO2,dn [gO2/m2/d], can be
computed as
J O 2,dn 1000

5.6.4

2
NO

mgN
3,1 T 20
T 20
NO3 NO3,1 NO3,2 NO
rondn
3 NO3,2
s

gN

(231)

Methane

The dissolved carbon generated by diagenesis is converted to methane in the anaerobic sediments.
Because methane is relatively insoluble, its saturation can be exceeded and methane gas
produced. As a consequence, rather than write a mass balance for methane in the anaerobic layer,
an analytical model developed by Di Toro et al. (1990) is used to determine the steady-state flux
of dissolved methane corrected for gas loss delivered to the aerobic sediments.
First, the carbon diagenesis flux is corrected for the oxygen equivalents consumed during
denitrification,
J CH 4,T J C J O 2, dn

(232)

where JCH4,T = the carbon diagenesis flux corrected for denitrification [gO2/m2/d]. In other words,
this is the total anaerobic methane production flux expressed in oxygen equivalents.
If JCH4,T is sufficiently large ( 2KL12Cs), methane gas will form. In such cases, the flux can be
corrected for the gas loss,
J CH 4,d 2 K L12 Cs J CH 4,T

QUAL2K

(233)

68

May 29, 2012

where JCH4,d = the flux of dissolved methane (expressed in oxygen equivalents) that is generated
in the anaerobic sediments and delivered to the aerobic sediments [gO2/m2/d], Cs = the saturation
concentration of methane expressed in oxygen equivalents [mgO2/L]. If JCH4,T < 2KL12Cs, then no
gas forms and
J CH 4,d J CH 4,T

(234)

The methane saturation concentration is computed as

H
Cs 100 1 1.02420T
10

(235)

where H = water depth [m] and T = water temperature [oC].

A methane mass balance can then be written for the aerobic layer as
H1

dCH 4,1
dt

J CH 4,d s c f CH 4,1

2
CH
4,1

T 20
CH
4 CH 4,1

(236)

where CH4,1 = methane concentration in the aerobic layer [gO2/m3], cf = fast CBOD in the
overlying water [gO2/m3], CH4,1 = the reaction velocity for methane oxidation in the aerobic
sediments [m/d], and CH4 = temperature correction factor [dimensionless]. At steady, state, this
balance can be solved for

J CH 4,d sc f

CH 4,1

2
CH
4 ,1

(237)

T 20
CH
4

The flux of methane to the overlying water, JCH4 [gO2/m2/d], can then be computed as

J CH 4 s CH 4,1 c f

5.6.5

(238)

SOD

The SOD [gO2/m2/d] is equal to the sum of the oxygen consumed in methane oxidation and
nitrification,

SOD CSOD NSOD

(239)

where CSOD = the amount of oxygen demand generated by methane oxidation [gO2/m2/d] and
NSOD = the amount of oxygen demand generated by nitrification [gO2/m2/d]. These are
computed as
CSOD

QUAL2K

2
CH
4,1

T 20
CH
4 CH 4,1

(240)

69

May 29, 2012

NSOD ron

2
NH
4,1

T 20
NH
4

K NH 4
o
f da1 NH 4,1
K NH 4 NH 4,1 2 K NH 4,O 2 o

(241)

where ron = the ratio of oxygen to nitrogen consumed during nitrification [= 4.57 gO2/gN].
5.6.6

Inorganic Phosphorus

Mass balances can be written total inorganic phosphorus in the aerobic layer and the anaerobic
layers as
H1

dPO4,1
dt

12 f pp 2 PO4,2 f pp1 PO4,1 K L12 f dp 2 PO4,2 f dp1 PO4,1

(242)

w2 PO4,1 s i f da1 PO4,1

1000

H2

dPO4,2
dt

J P 12 f pp1 PO4,1 f pp 2 PO4,2 K L12 f dp1 PO4,1 f dp 2 PO4,2

(243)

w2 PO4,1 PO4,2

where PO4,1 and PO4,2 = the concentration of total inorganic phosphorus in the aerobic layer and
the anaerobic layers, respectively [gP/m3], pi = the inorganic phosphorus in the overlying water
[mgP/m3], and JP = the diagenesis flux of phosphorus [gP/m2/d].
The fraction of phosphorus in dissolved (fdpi) and particulate (fppi) form are computed as
f dpi

1
1 mi pi

(244)

f ppi 1 f dpi

(245)

where pi = the partition coefficient for inorganic phosphorus in layer i [m3/gD].

The partition coefficient in the anaerobic layer is set to an input value. For the aerobic layer,
if the oxygen concentration in the overlying water column exceeds a critical concentration, ocrit
[gO2/m3], then the partition coefficient is increased to represent the sorption of phosphorus onto
iron oxyhydroxides as in
p1 p 2 PO 4,1

(246)

where PO4,1 is a factor that increases the aerobic layer partition coefficient relative to the
anaerobic coefficient.
If the oxygen concentration falls below ocrit then the partition coefficient is decreased
smoothly until it reaches the anaerobic value at zero oxygen,

p1 p 2 PO 4,1 o / ocrit

(247)

Equations (242) and (243) can be solved for PO4,1 and PO4,2. The flux of phosphorus to the
overlying water can then be computed as
QUAL2K

70

May 29, 2012

J PO 4 s PO4,1 i
1000

5.6.7

(248)

Solution Scheme

Although the foregoing sequence of equations can be solved, a single computation will not yield a
correct result because of the interdependence of the equations. For example, the surface mass
transfer coefficient s depends on SOD. The SOD in turn depends on the ammonium and methane
concentrations which themselves are computed via mass balances that depend on s. Hence, an
iterative technique must be used. The procedure used in QUAL2K is
1. Determine the diagenesis fluxes: JC, JN and JP.
2. Start with an initial estimate of SOD,
'
SODinit J C ron
JN

(249)

where ron = the ratio of oxygen to nitrogen consumed for total conversion of ammonium
to nitrogen gas via nitrification/denitrification [= 1.714 gO2/gN]. This ratio accounts for
the carbon utilized for denitrification.
3. Compute s using
s

SODinit
o

(250)

4. Solve for ammonium, nitrate and methane, and compute the CSOD and NSOD.
5. Make a revised estimate of SOD using the following weighted average
SOD

SODinit CSOD NSOD

2

(251)

6. Check convergence by calculating an approximate relative error

a

SOD SODinit
100%
SOD

(252)

7. If a is greater than a prespecified stopping criterion s then set SODinit = SOD and return
to step 2.
8. If convergence is adequate (a s), then compute the inorganic phosphorus
concentrations.
9. Compute the ammonium, nitrate, methane and phosphate fluxes.
5.6.8

Supplementary Fluxes

QUAL2K

71

May 29, 2012

Because of the presence of organic matter deposited prior to the summer steady-state period (e.g.,
during spring runoff), it is possible that the downward flux of particulate organic matter is
insufficient to generate the observed SOD. In such cases, a supplementary SOD can be
prescribed,

SODt SOD SODs

(253)

where SODt = the total sediment oxygen demand [gO2/m2/d], and SODs = the supplementary SOD
[gO2/m2/d]. In addition, prescribed ammonia and methane fluxes can be used to supplement the
computed fluxes.

QUAL2K

72

May 29, 2012

REFERENCES
Adams, E.E., D. J. Cosler, and K.R. Helfrich. 1987. "Analysis of evaporation data from heated
ponds", cs-5171, research project 2385-1, Electric Power Research Institute, Palo Alto,
California 94304. April.
Andrews and Rodvey, 1980. Heat exchange between water and tidal flats. D.G.M 24(2). (in
German).
APHA, 1995. Standard methods for the examination of water and wastewater, 19th Edn.
American Public Health Association, American Water Works Association and Water
Environment Federation: Washington, D.C.
Asaeda, T., and T. Van Bon, Modelling the effects of macrophytes on algal blooming in eutrophic
shallow lakes. Ecol. Model., 104: 261-287, 1997.
Baker, K.S. and Frouin, R. 1987. Relation between photosynthetically available radiation and
total insolation at the ocean surface under clear skies. Limnol. Oceanogr. 1370-1377.
Baly, E.C.C. 1935. The Kinetics of Photosynthesis. Proc. Royal Soc. London Ser. B, 117:218239.
Banks, R. B. 1975. "Some Features of Wind Action on Shallow Lakes." J. Environ Engr. Div.
ASCE. 101(EE5): 813-827.
Banks, R. B. and Herrera, F. F. 1977. "Effect of Wind and Rain on Surface Reaeration." J.
Environ Engr. Div. ASCE. 103(EE3): 489-504.
Barnwell, T.O., Brown, L.C., and Mareck, W. 1989). Application of Expert Systems Technology
in Water Quality Modeling. Water Sci. Tech. 21(8-9):1045-1056.
Bejan, A. 1993. Heat Transfer. Wiley, New York, NY.
Beven, K.J., Gilman, K., and Newson, M. 1979. Flow and flow routing in upland channel
networks. Hydrol. Sci. Bull. 24:43-69.
Bowie, G.L., Mills, W.B., Porcella, D.B., Campbell, C.L., Pagenkopf, J.R., Rupp, G.L., Johnson,
K.M., Chan, P.W.H., Gherini, S.A. and Chamberlin, C.E. 1985. Rates, Constants, and Kinetic
Formulations in Surface Water Quality Modeling. U.S. Envir. Prot. Agency, ORD, Athens,
GA, ERL, EPA/600/3-85/040.
Brady, D.K., Graves, W.L., and Geyer, J.C. 1969. Surface Heat Exchange at Power Plant Cooling
Lakes, Cooling Water Discharge Project Report, No. 5, Edison Electric Inst. Pub. No. 69-901,
New York, NY.
Bras, R.L. 1990. Hydrology. Addison-Wesley, Reading, MA.
Brown, L.C., and Barnwell, T.O. 1987. The Enhanced Stream Water Quality Models QUAL2E
and QUAL2E-UNCAS, EPA/600/3-87-007, U.S. Environmental Protection Agency, Athens,
GA, 189 pp.
Brunt, D. 1932. Notes on radiation in the atmosphere: I. Quart J Royal Meteorol Soc 58:389-420.
Brutsaert, W. 1982. Evaporation into the atmosphere: theory, history, and applications. D.Reidel
Publishing Co., Hingham MA, 299 p.
Butts, T. A. and Evans, R. L. 1983. Effects of Channel Dams on Dissolved Oxygen
Concentrations in Northeastern Illinois Streams, Circular 132, State of Illinois, Dept. of Reg.
and Educ., Illinois Water Survey, Urbana, IL.
Carslaw, H.S. and Jaeger, J.C. 1959. Conduction of Heat in Solids, Oxford Press, Oxford, UK,
510 pp.
Cengel, Y.A. 1998 Heat Transfer: A Practical Approach. New York, McGraw-Hill.
Chapra and Canale 2006. Numerical Methods for Engineers, 5th Ed. New York, McGraw-Hill.
Chapra, S.C. 1997. Surface water quality modeling. New York, McGraw-Hill.
Chapra, S.C. 2007. Applied Numerical Methods with MATLAB for Engineering and Science,
2nd Ed., WCB/McGraw-Hill, New York, N.Y.

QUAL2K

73

May 29, 2012

Chow, V.T., Maidment, D.R., and Mays, L.W. 1988. Applied Hydrology. New York, McGrawHill, 592 pp.
Churchill, M.A., Elmore, H.L., and Buckingham, R.A. 1962. The prediction of stream reaeration
rates. J. Sanit. Engrg. Div. , ASCE, 88{4),1-46.
Coffaro, G. and A. Sfriso, Simulation model of Ulva rigida growth in shallow water of the
Lagoon of Venice. Ecol. Model., 102: 55-66, 1997.
Covar, A. P. 1976. "Selecting the Proper Reaeration Coefficient for Use in Water Quality
Models." Presented at the U.S. EPA Conference on Environmental Simulation and Modeling,
April 19-22, 1976, Cincinnati, OH.
Di Toro, D. M. and J. F. Fitzpatrick. 1993. Chesapeake Bay sediment flux model. Tech. Report
EL-93-2, U.S. Army Corps of Engineers, Waterways Experiment Station, Vicksburg,
Mississippi, 316 pp.
Di Toro, D.M, Paquin, P.R., Subburamu, K. and Gruber, D.A. 1991. Sediment Oxygen Demand
Model: Methane and Ammonia Oxidation. J. Environ. Eng., 116(5):945-986.
Di Toro, D.M. 1978. Optics of Turbid Estuarine Waters: Approximations and Applications.
Water Res. 12:1059-1068.
Di Toro, D.M. 2001. Sediment Flux Modeling. Wiley-Interscience, New York, NY.
Droop, M.R. 1974. The nutrient status of algal cells in continuous culture. J.Mar.Biol.Assoc. UK
54:825-855.
Ecology. 2003. Shade.xls - a tool for estimating shade from riparian vegetation. Washington State
Department of Ecology. http://www.ecy.wa.gov/programs/eap/models/
Edinger, J.E., Brady, D.K., and Geyer, J.C. 1974. Heat Exchange and Transport in the
Environment. Report No. 14, EPRI Pub. No. EA-74-049-00-3, Electric Power Research
Institute, Palo Alto, CA.
Finlayson-Pitts, B.J. and Pitts, J.N., Jr. 2000. Chemistry of the upper and lower atmosphere :
theory, experiments, and applications. Academic Press.
Finnemore, E.J. and Franzini, J.B. 2002. Fluid Mechanics with Engineering Applications, 10th Ed.
New York, McGraw, Hill.
Geiger, R. 1965. The climate near the ground. Harvard University Press. Cambridge, MA.
Gordon, N.D, T.A. McMahon, and B.L. Finlayson. 1992. Stream Hydrology, An Introduction for
Ecologists. Published by John Wiley and Sons.
Grigull, U. and Sandner, H. 1984. Heat Conduction. Springer-Verlag, New York, NY.
Hamilton, D.P., and S.G. Schladow, Prediction of water quality in lakes and reservoirs. 1. Model
description. Ecol. Model., 96: 91-110, 1997.
Harbeck, G. E., 1962, A practical field technique for measuring reservoir evaporation utilizing
mass-transfer theory. US Geological Survey Professional Paper 272-E, 101-5.
Harned, H.S., and Hamer, W.J. 1933. The Ionization Constant of Water. J. Am., Chem. Soc.
51:2194.
Helfrich, K.R., E.E. Adams, A.L. Godbey, and D.R.F. Harleman. 1982. Evaluation of models for
predicting evaporative water loss in cooling impoundments. Report CS-2325, Research
project 1260-17. Electric Power Research Institute, Pala Alto, CA. March 1982.
Hellweger, F/L., K.L. Farley, U. Lall, and D.M. DiToro. 2003. Greedy algae reduce arsenate.
Limnol.Oceanogr. 48(6), 2003, 2275-2288.
Holland, H.D. 1978. The Chemistry of the Atmosphere and Oceans. Wiley-Interscience, NY, NY.
Hutchinson, G.E. 1957. A Treatise on Limnology, Vol. 1, Physics and Chemistry. Wiley, New
York, NY.
Jobson, H.E. 1977. Bed Conduction Computation for Thermal Models. J. Hydraul. Div. ASCE.
103(10):1213-1217.
Koberg, G.E. 1964. Methods to compute long-wave radiation from the atmosphere and reflected
solar radiation from a water surface. US Geological Survey Professional Paper 272-F.

QUAL2K

74

May 29, 2012

Kreith, F. and Bohn, M.S. 1986. Principles of Heat Transfer, 4th Ed. Harper and Row, New
York, NY.
Laws, E. A. and Chalup, M. S. 1990. A Microalgal Growth Model. Limnol. Oceanogr. 35(3):597608.
LI-COR, 2003. Radiation Measurement Instruments, LI-COR, Lincoln, NE, 30 pp.
Likens, G. E., and Johnson, N. M. (1969). "Measurements and analysis of the annual heat budget
for sediments of two Wisconsin lakes." Limnol. Oceanogr., 14(1):115-135.
Mackay, D. and Yeun, A.T.K. 1983. Mass Transfer Coefficient Correlations for Volatilization of
Organic Solutes from Water. Environ. Sci. Technol. 17:211-233.
Marciano, J.K. and G.E. Harbeck. 1952. Mass transfer studies in water loss investigation: Lake
Hefner studies. Geological Circular 229. U.S. Geological Survey, Washington DC.
McIntyre, C. D. 1973. Periphyton dynamics in laboratory streams: A simulation model and its
implications. Ecol. Monogr. 43:399-420.
Meeus, J. 1999. Astronomical algorithms. Second edition. Willmann-Bell, Inc. Richmond, VA.
Melching, C.S. and Flores, H.E. 1999. Reaeration Equations from U.S. Geological Survey
Database. J. Environ. Engin., 125(5):407-414.
Mills, A.F. 1992. Heat Transfer. Irwin, Homewood, IL.
Moog, D.B., and Iirka, G.H. 1998. Analysis of reaeration equations using mean multiplicative
error. J. Environ. Engrg., ASCE, 124(2), 104-110.
Munson, B.R., Young, D.F., Okiishi, T.H., and W.D. Huebsch. 2009. Fundamentals of Fluid
Mechanics. 6th edition. John Wiley & Sons. Hoboken, NJ.
Nakshabandi, G.A. and H. Kohnke. 1965. Thermal conductivity and diffusivity of soils as related
to moisture tension and other physical properties. Agr. Met. Vol 2.
O'Connor, D.J., and Dobbins, W.E. 1958. Mechanism of reaeration in natural streams. Trans.
ASCE, 123, 641-684.
Owens, M., Edwards, R.W., and Gibbs, J.W. 1964. Some reaeration studies in streams. Int. J. Air
and Water Pollution, 8, 469-486.
Plummer, L.N. and Busenberg, E. 1982. The Solubilities of Calcite, Aragonite and Vaterite in
CO2-H2O Solutions Between 0 and 90 oC, and an Evaluation of the Aqueous Model for the
System CaCO3CO2H2O. Geochim. Cosmochim. 46:1011-1040.
Raudkivi, A. I. 1979. Hydrology. Pergamon, Oxford, England.
Redfield, A.C., Ketchum, B.H. and Richards, F.A. 1963. The Influence of Organisms on the
Composition of Seawater, in The Sea, M.N. Hill, ed. Vol. 2, pp. 27-46, Wiley-Interscience,
NY.
Riley, G.A. 1956. Oceanography of Long Island Sound 1952-1954. II. Physical Oceanography,
Bull. Bingham Oceanog. Collection 15, pp. 15-16.
Rosgen, D. 1996. Applied river morphology. Wildland Hydrology publishers. Pagosa Springs,
CO.
Rutherford, J.C., Scarsbrook, M.R. and Broekhuizen, N. 2000. Grazer Control of Stream Algae:
Modeling Temperature and Flood Effects. J. Environ. Eng. 126(4):331-339.
Ryan, P.J. and D.R.F. Harleman. 1971. Prediction of the annual cycle of temperature changes in a
stratified lake or reservoir. Mathematical model and users manual. Ralph M. Parsons
Laboratory Report No. 137. Massachusetts Institute of Technology. Cambridge, MA.
Ryan, P.J. and K.D. Stolzenbach.1972. Engineering aspects of heat disposal from power
generation, (D.R.F. Harleman, ed.). R.M. Parson Laboratory for Water Resources and
Hydrodynamics, Department of Civil Engineering, Massachusetts Institute of Technology,
Cambridge, MA
Schwarzenbach, R.P., Gschwend, P.M., and Imboden, D.M. 1993. Environmental Organic
Chemistry. Wiley-Interscience, 681 pp.

QUAL2K

75

May 29, 2012

Shanahan, P. 1984. Water temperature modeling: a practical guide. In: Proceedings of stormwater
and water quality model users group meeting, April 12-13, 1984. USEPA, EPA-600/9-85003. (users.rcn.com/shanahan.ma.ultranet/TempModeling.pdf).
Smith, E.L. 1936. Photosynthesis in Relation to Light and Carbon Dioxide. Proc. Natl. Acad. Sci.
22:504-511.
Steele, J.H. 1962. Environmental Control of Photosynthesis in the Sea. Limnol. Oceanogr. 7:137150.
Stumm, W. and Morgan, J.J. 1996. Aquatic Chemistry, 3rd Ed., New York, Wiley-Interscience,
1022 pp.
Szeicz, G. 1974. Solar radiation for plant growth. J. Appl. Ecol. 11:617-636.
Thackston, E. L., and Dawson, J.W. III. 2001. Recalibration of a Reaeration Equation. J. Environ.
Engrg., 127(4):317-320.
Thackston, E. L., and Krenkel, P. A. 1969. Reaeration-prediction in natural streams. J. Sanit.
Engrg. Div., ASCE, 95(1), 65-94.
Tsivoglou, E. C., and Neal, L.A. 1976. Tracer Measurement of Reaeration. III. Predicting the
Reaeration Capacity of Inland Streams. Journal of the Water Pollution Control Federation,
48(12):2669-2689.
Thomann, R.V. and Mueller, J.A. 1987. Principles of Surface Water Quality Modeling and
Control. New York, Harper-Collins.
TVA, 1972. Heat and mass transfer between a water surface and the atmosphere. Water
Resources Research, Laboratory Report No. 14. Engineering Laboratory, Division of Water
Control Planning, Tennessee Valley Authority, Norris TN.
Wanninkhof, R., Ledwell, I. R., and Crusius, I. 1991. "Gas Transfer Velocities on Lakes
Measured with Sulfur Hexafluoride." In Symposium Volume of the Second International
Conference on Gas Transfer at Water Surfaces, S.C. Wilhelms and I.S. Gulliver, eds.,
Minneapolis, MN .
Wood, K.G. 1974. Carbon Dioxide Diffusivity Across the Air-Water Interface. Archiv fr
Hydrobiol. 73(1):57-69.

QUAL2K

76

May 29, 2012

APPENDIX A: NOMENCLATURE
Symbol

Units
acres

Definition

Aacres ,i

surface area of element i

[CO2]s

saturation concentration of carbon dioxide

carbonate ion
hydronium ion
sum of dissolved carbon dioxide and carbonic acid
bicarbonate ion
hydroxyl ion
velocity rating curve coefficient
mean atmospheric transmission coefficient after
scattering and absorption
mean atmospheric transmission coefficient

[CO32]
[H+]
[H2CO3*]
[HCO3]
[OH]

a
a"
a
a1
Aa
Ab

ab
Ac
ad
Alk
ap
at
atc

B
b
B0
bd
c1
cf
Cgb(T)
CH4,1
CL
cnps,i,j
Cp
cps,i,j
cs
Cs
CSOD
cT
c'T,i1
d
Dd
Dp

QUAL2K

atmospheric molecular scattering coefficient for radiation

transmission
atmospheric long-wave radiation coefficient

atmospheric long-wave radiation coefficient

bottom algae
cross-sectional area
coefficient correcting dam reaeration for water quality
alkalinity
phytoplankton concentration
atmospheric attenuation
atmospheric transmission coefficient
average element width
velocity rating curve exponent
bottom width
coefficient correcting dam reaeration for dam type
Bowen's coefficient
fast reacting CBOD
temperature-dependent maximum photosynthesis rate
methane concentration in the aerobic sediment layer
fraction of sky covered with clouds
jth non-point source concentration for element i
specific heat of water
jth point source concentration for element i
slowly reacting CBOD
saturation concentration of methane
the amount of oxygen demand generated by methane
oxidation
total inorganic carbon
concentration of inorganic carbon entering element
below a dam
dust attenuation coefficient
pore water diffusion coefficient
diffusion coefficient for bioturbation

77

mole/L
mole/L
mole/L
mole/L
mole/L
mole/L
dimensionless
Dimensionless
Dimensionless
dimensionless
dimensionless
mmHg-0.5 or mb-0.5
mgA/m2
m2
dimensionless
eq L1 or
mgCaCO3/L
mgA/m3
dimensionless
dimensionless
m
dimensionless
m
dimensionless
0.47 mmHg/oC
mgO2/L
mgA/(m2 d)
gO2/m3
dimensionless
o
C
cal/(g oC)
o
C
mgO2/L
mgO2/L
gO2/m2/d
mole L1
mgO2/L
dimensionless
m2/d
m2/d

May 29, 2012

Ei
eair
elev
En
Ep,i

JCH4,T
Je
Ja
JN

bulk dispersion coefficient between elements i and i + 1

air vapor pressure
elevation above sea level
numerical dispersion
longitudinal dispersion between elements i and i + 1
equation of time: the difference between mean solar time
and true solar time when located on the reference
longitude of the time zone considered
saturation vapor pressure at water surface
photoperiod
fraction of ammonium in dissolved form in sediment
layer i
fraction of inorganic phosphorus in dissolved form in
sediment layer i
phytoplankton growth attenuation due to light
attenuation of CBOD oxidation due to low oxygen
enhancement of denitrification at low oxygen
concentration
attenuation due to low oxygen
attenuation due to low oxygen
attenuation due to low oxygen
fraction of ammonium in particulate form in sediment
layer i
fraction of inorganic phosphorus in particulate form in
sediment layer i
fraction unionized ammonia
acceleration due to gravity
mass of element X
water depth
drop in water elevation for a dam
thickness of sediment layer i
sediment thickness
solar radiation at water surface
extraterrestrial radiation
net atmospheric longwave radiation flux
longwave back radiation flux from water
conduction flux
flux of labile dissolved carbon
flux of methane from sediments to the overlying water
flux of dissolved methane generated in anaerobic
sediments corrected for methane gas formation
carbon diagenesis flux corrected for denitrification
evaporation flux
air-water heat flux
nitrogen diagenesis flux

JO2,dn

oxygen equivalents consumed during denitrification

JP
JPOC,G1

phosphorus diagenesis flux

the flux of labile POC delivered to the anaerobic

eqtime

es
f
fdai
fdpi
FLp
Foxc
Foxdn
Foxna
Foxrb
Foxrp
fpai
fppi
Fu

g
gX

H
Hd
Hi
Hsed
I(0)
I0
Jan
Jbr
Jc
JC,G1
JCH4
JCH4,d

QUAL2K

78

m3/d
mm Hg
m
m2/d

m2/s
minutes

mmHg
fraction of day
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless

= 9.81 m/s2
g
m
m
m
cm
cal/cm2/d
cal/cm2/d
cal/(cm2 d)
cal/(cm2 d)
cal/(cm2 d)
gO2/m2/d
gO2/m2/d
gO2/m2/d
gO2/m2/d
cal/(cm2 d)
cal/(cm2 d)
gN/m2/d
gO2/m2/d
gP/m2/d
gO2/m2/d

May 29, 2012

JPOM
Jsi
Jsn
k(T)
K1
K2
Ka

ka(T)
kac(T)
kdb(T)
kdc(T)
kdn(T)
kdp(T)
kdt(T)
kdX(T)
ke
keb
kgp(T)
KH
khc(T)
khn(T)
khnxb
khnxp
khp(T)
KL12

KLb
KLp
KM,Dp
kna(T)
KNH4
KNH4,O2
kPOC,G1
krb(T)
krp(T)
ksNb
ksNp
Ksocf
Ksodn
Ksona
ksPb
ksPp
Kw
QUAL2K

sediment layer
the downward flux of POM
sediment-water heat flux
net solar shortwave radiation flux
temperature dependent first-order reaction rate
acidity constant for dissociation of carbonic acid
acidity constant for dissociation of bicarbonate
equilibrium coefficient for ammonium dissociation
temperature-dependent oxygen reaeration coefficient
temperature-dependent carbon dioxide reaeration
coefficient
temperature-dependent bottom algae death rate
temperature-dependent fast CBOD oxidation rate
temperature-dependent denitrification rate
temperature-dependent phytoplankton death rate
temperature-dependent detritus dissolution rate
temperature-dependent pathogen die-off rate
light extinction coefficient
a background coefficient accounting for extinction due to
water and color
maximum photosynthesis rate at temperature T
Henry's constant
temperature-dependent slow CBOD hydrolysis rate
temperature-dependent organic nitrogen hydrolysis rate
preference coefficient of bottom algae for ammonium
preference coefficient of phytoplankton for ammonium
temperature-dependent organic phosphorus hydrolysis
rate
pore water diffusion mass transfer coefficient

bottom algae light parameter

phytoplankton light parameter
oxygen half-saturation constant for bioturbation
temperature-dependent nitrification rate for ammonia
nitrogen
ammonium half-saturation constant
oxygen half-saturation constant
the mineralization rate of labile POC
temperature-dependent bottom algae respiration rate
temperature-dependent phytoplankton
respiration/excretion rate
nitrogen half-saturation constant for bottom algae
nitrogen half-saturation constant for phytoplankton
parameter for oxygen dependency of fast CBOD
oxidation
parameter for oxygen dependency of denitrification
parameter for oxygen dependency of nitrification
phosphorus half-saturation constant for bottom algae
phosphorus half-saturation constant for phytoplankton
acidity constant for dissociation of water

79

gD m2 d1
cal/(cm2 d)
cal/(cm2 d)
/d

/d
/d
/d
/d
/d
/d
/d
/d
/m1
/m
/d
mole/(L atm)
/d
/d
mgN/m3
mgN/m3
/d
m/d
ly/d
gO2/m3
/d
gN/m3
mgO2/L
d1
/d
/d
gN/L
gN/L

gP/L
gP/L

May 29, 2012

Lat
Llm
localTime
Lsm
m
mgY

mi
mo
n
na
nau
nfac

NH4,i
nn
no
NO3,i
npai
npsi
NSOD
o
oi1
ocrit
os(T, elev)
P
Pab
pai
Pap
PAR(z)
patm

pCO2
pi
po
PO4,i

POC2,G1
POCR
psi
pwc
Q
Qout,i
Qi

QUAL2K

latitude
longitude of local meridian
local standard time
longitude of standard meridian
optical air mass
mass of element Y
inorganic suspended solids
detritus
Manning roughness coefficient
the ammonium concentration in the overlying water
unionized ammonia nitrogen
atmospheric turbidity factor
the concentration of total ammonium in sediment layer i
nitrate concentration in the overlying water
organic nitrogen
nitrate concentration in layer i
total number of non-point withdrawals outflows from
element i
total number of non-point sources inflows to element i
the amount of oxygen demand generated by nitrification
dissolved oxygen
oxygen concentration entering an element below a dam
critical oxygen concentration for sediment phosphorus
sorption
saturation concentration of oxygen at temperature, T, and
elevation above sea level, elev
wetted perimeter
preference for ammonium as a nitrogen source for
bottom algae
total number of point withdrawals from element i
preference for ammonium as a nitrogen source for
phytoplankton
photosynthetically active radiation (PAR) at depth z
below water surface
atmospheric pressure
atmospheric partial pressure of carbon dioxide
inorganic phosphorus
organic phosphorus
the concentration of total inorganic phosphorus in
sediment layer i
the concentration of the labile fraction of POC in the
anaerobic sediment layer
reference G1 concentration for bioturbation
total number of point sources to element i
mean daily atmospheric precipitable water content
flow
total outflow from element due to point and nonpoint
withdrawals
outflow from element i into element i + 1

80

radians
degrees
minutes
degrees
dimensionless
mg
mgD/L
mgD/L
mgN/m3
mgN/m3
dimensionless
gN/m3
mgN/m3
mgN/m3
gN/m3

dimensionless
dimensionless
gO2/m2/d
mgO2/L
mgO2/L
gO2/m3
mgO2/L
m
dimensionless

dimensionless
dimensionless
ly/d
mm Hg
atm
gP/L
gP/L
gP/m3
gO2/m3
gC/m3

dimensionless
m3/s or m3/d
m3/d
m3/d

May 29, 2012

Qin,i
Qnpa,i,j
Qnps,i,j
Qpa,i,j
Qps,i,j
r
rcndn
rd
rda
RL
roc
roca
rocn
ron
ron'
Rs
S
s
s

S0
Sb,i
Si
SOD
SODs
SODt
ss
t
T
Tw,f
Ta
Tair
Tair,f
Td

Ti
timezone

QUAL2K

total inflow into element from point and nonpoint

sources
jth non-point withdrawal outflow from element i
jth non-point source inflow to element i
jth point withdrawal outflow from element i
jth point source inflow to element i
normalized radius of earths orbit (i.e., ratio of actual
earth-sun distance to mean earth-sun distance
ratio of oxygen equivalents lost per nitrate nitrogen that
is denitrified
ratio of deficit above and below dam
the ratio of dry weight to chlorophyll a
longwave reflection coefficient
ratio of oxygen consumed per organic carbon oxidized to
carbon dioxide
ratio of oxygen generated per organic carbon produced
by photosynthesis when nitrate is taken up
ratio of oxygen generated per organic carbon produced
by photosynthesis when ammonia is taken up
the ratio of oxygen to nitrogen consumed during
nitrification
the ratio of oxygen to nitrogen consumed for total
conversion of ammonium to nitrogen gas via
nitrification/denitrification
albedo or reflectivity (fraction of solar radiation that is
reflected)
channel slope
chloride
mass transfer coefficient between the water and the
aerobic sediments
bottom slope
sources and sinks of constituent due to reactions for
bottom algae
sources and sinks of constituent due to reactions and
mass transfer mechanisms for water constituents
the sediment oxygen demand
the supplemental sediment oxygen demand
total sediment oxygen demand = SOD + SODs
channel side slope
time
water temperature
water temperature
absolute temperature
air temperature
air temperature
dew-point temperature
temperature in element i
time zone indicates local time zone in relation to
Greenwich Mean Time (GMT) (negative in western

81

m3/d
m3/d
m3/d
m3/d
m3/d
dimensionless
gO2/gN
dimensionless
gD/mgA
dimensionless
gO2/gC
gO2/gC
gO2/gC
= 4.57 gO2/gN
= 1.714 gO2/gN

dimensionless
dimensionless
mgCl/L
m/d
m/m
mgA/m2/d
g/m3/d or mg/m3/d
gO2/m2/d
gO2/m2/d
gO2/m2/d
m/m
d
o
C
o
F
K
o

C
F
o
C
o
C
hours
o

May 29, 2012

Tnps,i,j
Tps,i,j
trueSolarTime

Ts,i
tsr
tss
Tstd
tt,i
U
U*
Uw
Uw,mph
Uw,z
va
vdt
vi
Vi
vp
vpc
W0
w2
Wh,i
Wi
X
zw

hemisphere)
jth non-point source temperature for element i
jth point source temperature for element i
time determined from actual position of the sun in the
sky
temperature of bottom sediment
time of sunrise
time of sunset
standard time
travel time from headwater to end of element i
mean velocity
shear velocity
wind speed
wind speed
wind speed at height zw above water surface
phytoplankton settling velocity
detritus settling velocity
inorganic suspended solids settling velocity
volume of ith element
pathogen settling velocity
non-living particulate organic carbon settling velocity
solar constant
burial velocity
net heat load from point and non-point sources into
element i
external loading of constituent to element i
pathogen
height above water for wind speed measurements

C
C
minutes
o

C
Hr
Hr
Hr
D
m/s
m/s
m/s
mph
m/s
m/d
m/d
m/d
m3
m/d
m/d
2851 cal/cm2/d
m/d
cal/d
g/d or mg/d
cfu/100 mL
m

Greek:
Symbol

d
s

0
1
2
i
o
p
pn
path

PO4,1

QUAL2K

Definition
depth rating curve coefficient
suns altitude
suns altitude
sediment thermal diffusivity
fraction of total inorganic carbon in carbon dioxide
fraction of total inorganic carbon in bicarbonate
fraction of total inorganic carbon in carbonate
effect of inorganic suspended solids on light attenuation
effect of particulate organic matter on light attenuation
linear effect of chlorophyll on light attenuation
non-linear effect of chlorophyll on light attenuation
pathogen light efficiency factor
depth rating curve exponent
solar declination
factor that increases the aerobic sediment layer partition
coefficient relative to the anaerobic coefficient

82

Units
dimensionless
radians
degrees
cm2/s
dimensionless
dimensionless
dimensionless
L/mgD/m
L/mgD/m
L/gA/m
(L/gA)2/3/m
dimensionless
dimensionless
radians
dimensionless

May 29, 2012

v
ts
xi

clear
sky
a

Lb
Lp
Nb
Np
p

am
12
ai
pi
CH4
Dd
Dp
NH4
NO3
POC,G1
CH4,1
NH4,1
NO3,i

QUAL2K

virtual temperature difference between the water and air

difference between standard and local civic time
length of ith element
emissivity of water
emissivity of longwave radiation from the sky with no
clouds
emissivity of longwave radiation from the sky with
clouds
estimated error
bottom algae light attenuation
phytoplankton light attenuation
bottom algae nutrient attenuation factor
phytoplankton nutrient attenuation factor
phytoplankton photosynthesis rate
density of water
Stefan-Boltzmann constant
local hour angle of sun
residence time of ith element
temperature coefficient for zero and first-order reactions
elevation adjusted optical air mass
the bioturbation mass transfer coefficient between the
sediment layers
the partition coefficient for ammonium in sediment layer
i
the partition coefficient for inorganic phosphorus in
sediment layer i
temperature correction factor for sediment methane
oxidation
temperature coefficient for porewater diffusion
temperature coefficient for bioturbation diffusion
temperature correction factor for sediment nitrification
sediment denitrification temperature correction factor
temperature correction factor for labile POC
mineralization
the reaction velocity for methane oxidation in the
aerobic sediments
the reaction velocity for nitrification in the aerobic
sediments
denitrification reaction velocity sediment layer i

83

F
hr

m
dimensionless
0-1
0-1
%
0-1
0-1
0-1
0-1
/d
g/cm3
11.710-8 cal/(cm2 d K4)
radians
d
dimensionless
m/d
m3/gD
m3/gD
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
dimensionless
m/d
m/d
m/d

May 29, 2012

APPENDIX B: SOLAR POSITION, SUNRISE, AND SUNSET

CALCULATIONS
The sunrise/sunset and solar position functions are a VBA translation of NOAA's sunrise/sunset calculator
and NOAA's solar position calculator at the following web pages:
http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html
http://www.srrb.noaa.gov/highlights/sunrise/azel.html
The calculations in the NOAA Sunrise/Sunset and Solar Position Calculators are based on equations from
Astronomical Algorithms, by Jean Meeus. The sunrise and sunset results have been verified by NOAA to
be accurate to within a minute for locations between +/- 72 latitude, and within 10 minutes outside of
those latitudes.
Five main functions are included for use from Excel worksheets or VBA programs:
sunrise (lat, lon, year, month, day, timezone, dlstime) calculates the local time of sunrise for a location
and date
solarnoon (lat, lon, year, month, day, timezone, dlstime) calculates the local time of solar noon for a
location and date (the time when the sun crosses the meridian)
sunset (lat, lon, year, month, day, timezone, dlstime) calculates the local time of sunset for a location
and date
solarazimuth (lat, lon, year, month, day, hour, minute, second, timezone, dlstime) calculates the solar
azimuth for a location, date, and time (degrees clockwise from north to the point on the horizon
directly below the sun)
solarelevation (lat, lon, year, month, day, hour, minute, second, timezone, dlstime) calculates the solar
elevation for a location, date, and time (degrees vertically from horizon to the sun)
A subroutine is also provided that calculates solar azimuth (az), solar elevation (el):
solarposition (lat, lon, year, month, day, hour, minute, second, timezone, dlstime, az, el,
earthRadiusVector)
The sign convention for the main functions and subroutine is:
positive latitude decimal degrees for northern hemisphere
negative longitude degrees for western hemisphere
negative time zone hours for western hemisphere
The Excel/VBA functions and subroutines for solar position, sunrise, and sunset times are as follows:
Option Explicit
Function radToDeg(angleRad)
'// Convert radian angle to degrees
radToDeg = (180# * angleRad / Application.WorksheetFunction.Pi())
End Function
Function degToRad(angleDeg)
'// Convert degree angle to radians
degToRad = (Application.WorksheetFunction.Pi() * angleDeg / 180#)
End Function
Function calcJD(year, month, day)
'***********************************************************************/
'* Name:
calcJD
'* Type:
Function
'* Purpose: Julian day from calendar day
'* Arguments:
'*
year : 4 digit year
'*
month: January = 1
'*
day : 1 - 31
'* Return value:
'*
The Julian day corresponding to the date
'* Note:
'*
Number is returned for start of day. Fractional days should be
'*
added later.
'***********************************************************************/
Dim A As Double, B As Double, jd As Double

QUAL2K

84

May 29, 2012

If (month <= 2) Then

year = year - 1
month = month + 12
End If
A = Application.WorksheetFunction.Floor(year / 100, 1)
B = 2 - A + Application.WorksheetFunction.Floor(A / 4, 1)
jd = Application.WorksheetFunction.Floor(365.25 * (year + 4716), 1) + _
Application.WorksheetFunction.Floor(30.6001 * (month + 1), 1) + day + B - 1524.5
calcJD = jd
'gp put the year and month back where they belong
If month = 13 Then
month = 1
year = year + 1
End If
If month = 14 Then
month = 2
year = year + 1
End If
End Function
Function calcTimeJulianCent(jd)
'***********************************************************************/
'* Name:
calcTimeJulianCent
'* Type:
Function
'* Purpose: convert Julian Day to centuries since J2000.0.
'* Arguments:
'*
jd : the Julian Day to convert
'* Return value:
'*
the T value corresponding to the Julian Day
'***********************************************************************/
Dim t As Double
t = (jd - 2451545#) / 36525#
calcTimeJulianCent = t
End Function
Function calcJDFromJulianCent(t)
'***********************************************************************/
'* Name:
calcJDFromJulianCent
'* Type:
Function
'* Purpose: convert centuries since J2000.0 to Julian Day.
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
the Julian Day corresponding to the t value
'***********************************************************************/
Dim jd As Double
jd = t * 36525# + 2451545#
calcJDFromJulianCent = jd
End Function
Function calcGeomMeanLongSun(t)
'***********************************************************************/
'* Name:
calGeomMeanLongSun
'* Type:
Function
'* Purpose: calculate the Geometric Mean Longitude of the Sun
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
the Geometric Mean Longitude of the Sun in degrees
'***********************************************************************/
Dim l0 As Double
l0 = 280.46646 + t * (36000.76983 + 0.0003032 * t)
Do
If (l0 <= 360) And (l0 >= 0) Then Exit Do
If l0 > 360 Then l0 = l0 - 360
If l0 < 0 Then l0 = l0 + 360
Loop
calcGeomMeanLongSun = l0
End Function
Function calcGeomMeanAnomalySun(t)
'***********************************************************************/
'* Name:
calGeomAnomalySun
'* Type:
Function
'* Purpose: calculate the Geometric Mean Anomaly of the Sun
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
the Geometric Mean Anomaly of the Sun in degrees
'***********************************************************************/
Dim m As Double
m = 357.52911 + t * (35999.05029 - 0.0001537 * t)
calcGeomMeanAnomalySun = m
End Function
Function calcEccentricityEarthOrbit(t)
'***********************************************************************/
'* Name:
calcEccentricityEarthOrbit
'* Type:
Function
'* Purpose: calculate the eccentricity of earth's orbit
'* Arguments:

QUAL2K

85

May 29, 2012

'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
the unitless eccentricity
'***********************************************************************/
Dim e As Double
e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t)
calcEccentricityEarthOrbit = e
End Function
Function calcSunEqOfCenter(t)
'***********************************************************************/
'* Name:
calcSunEqOfCenter
'* Type:
Function
'* Purpose: calculate the equation of center for the sun
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
in degrees
'***********************************************************************/
Dim m As Double, mrad As Double, sinm As Double, sin2m As Double, sin3m As Double
Dim c As Double
m = calcGeomMeanAnomalySun(t)
mrad = degToRad(m)
sinm = Sin(mrad)
sin2m = Sin(mrad + mrad)
sin3m = Sin(mrad + mrad + mrad)
c = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) _
+ sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289
calcSunEqOfCenter = c
End Function
Function calcSunTrueLong(t)
'***********************************************************************/
'* Name:
calcSunTrueLong
'* Type:
Function
'* Purpose: calculate the true longitude of the sun
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
sun's true longitude in degrees
'***********************************************************************/
Dim l0 As Double, c As Double, O As Double
l0 = calcGeomMeanLongSun(t)
c = calcSunEqOfCenter(t)
O = l0 + c
calcSunTrueLong = O
End Function
Function calcSunTrueAnomaly(t)
'***********************************************************************/
'* Name:
calcSunTrueAnomaly (not used by sunrise, solarnoon, sunset)
'* Type:
Function
'* Purpose: calculate the true anamoly of the sun
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
sun's true anamoly in degrees
'***********************************************************************/
Dim m As Double, c As Double, v As Double
m = calcGeomMeanAnomalySun(t)
c = calcSunEqOfCenter(t)
v = m + c
calcSunTrueAnomaly = v
End Function
Function calcSunRadVector(t)
'***********************************************************************/
'* Name:
calcSunRadVector (not used by sunrise, solarnoon, sunset)
'* Type:
Function
'* Purpose: calculate the distance to the sun in AU
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
sun radius vector in AUs
'***********************************************************************/
Dim v As Double, e As Double, r As Double
v = calcSunTrueAnomaly(t)
e = calcEccentricityEarthOrbit(t)
r = (1.000001018 * (1 - e * e)) / (1 + e * Cos(degToRad(v)))
calcSunRadVector = r
End Function
Function calcSunApparentLong(t)
'***********************************************************************/
'* Name:
calcSunApparentLong (not used by sunrise, solarnoon, sunset)
'* Type:
Function
'* Purpose: calculate the apparent longitude of the sun
'* Arguments:

QUAL2K

86

May 29, 2012

'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
sun's apparent longitude in degrees
'***********************************************************************/
Dim O As Double, omega As Double, lambda As Double
O = calcSunTrueLong(t)
omega = 125.04 - 1934.136 * t
lambda = O - 0.00569 - 0.00478 * Sin(degToRad(omega))
calcSunApparentLong = lambda
End Function
Function calcMeanObliquityOfEcliptic(t)
'***********************************************************************/
'* Name:
calcMeanObliquityOfEcliptic
'* Type:
Function
'* Purpose: calculate the mean obliquity of the ecliptic
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
mean obliquity in degrees
'***********************************************************************/
Dim seconds As Double, e0 As Double
seconds = 21.448 - t * (46.815 + t * (0.00059 - t * (0.001813)))
e0 = 23# + (26# + (seconds / 60#)) / 60#
calcMeanObliquityOfEcliptic = e0
End Function
Function calcObliquityCorrection(t)
'***********************************************************************/
'* Name:
calcObliquityCorrection
'* Type:
Function
'* Purpose: calculate the corrected obliquity of the ecliptic
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
corrected obliquity in degrees
'***********************************************************************/
Dim e0 As Double, omega As Double, e As Double
e0 = calcMeanObliquityOfEcliptic(t)
omega = 125.04 - 1934.136 * t
e = e0 + 0.00256 * Cos(degToRad(omega))
calcObliquityCorrection = e
End Function
Function calcSunRtAscension(t)
'***********************************************************************/
'* Name:
calcSunRtAscension (not used by sunrise, solarnoon, sunset)
'* Type:
Function
'* Purpose: calculate the right ascension of the sun
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
sun's right ascension in degrees
'***********************************************************************/
Dim e As Double, lambda As Double, tananum As Double, tanadenom As Double
Dim alpha As Double
e = calcObliquityCorrection(t)
lambda = calcSunApparentLong(t)
tananum = (Cos(degToRad(e)) * Sin(degToRad(lambda)))
tanadenom = (Cos(degToRad(lambda)))
'original NOAA code using javascript Math.Atan2(y,x) convention:
'
var alpha = radToDeg(Math.atan2(tananum, tanadenom));
'
alpha = radToDeg(Application.WorksheetFunction.Atan2(tananum, tanadenom))
'translated using Excel VBA Application.WorksheetFunction.Atan2(x,y) convention:
alpha = radToDeg(Application.WorksheetFunction.Atan2(tanadenom, tananum))
calcSunRtAscension = alpha
End Function
Function calcSunDeclination(t)
'***********************************************************************/
'* Name:
calcSunDeclination
'* Type:
Function
'* Purpose: calculate the declination of the sun
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
sun's declination in degrees
'***********************************************************************/
Dim e As Double, lambda As Double, sint As Double, theta As Double
e = calcObliquityCorrection(t)
lambda = calcSunApparentLong(t)
sint = Sin(degToRad(e)) * Sin(degToRad(lambda))
theta = radToDeg(Application.WorksheetFunction.Asin(sint))
calcSunDeclination = theta
End Function
Function calcEquationOfTime(t)

QUAL2K

87

May 29, 2012

'***********************************************************************/
'* Name:
calcEquationOfTime
'* Type:
Function
'* Purpose: calculate the difference between true solar time and mean
'*
solar time
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'* Return value:
'*
equation of time in minutes of time
'***********************************************************************/
Dim epsilon As Double, l0 As Double, e As Double, m As Double
Dim y As Double, sin2l0 As Double, sinm As Double
Dim cos2l0 As Double, sin4l0 As Double, sin2m As Double, Etime As Double
epsilon = calcObliquityCorrection(t)
l0 = calcGeomMeanLongSun(t)
e = calcEccentricityEarthOrbit(t)
m = calcGeomMeanAnomalySun(t)
y = Tan(degToRad(epsilon) / 2#)
y = y ^ 2
sin2l0 = Sin(2# * degToRad(l0))
sinm = Sin(degToRad(m))
cos2l0 = Cos(2# * degToRad(l0))
sin4l0 = Sin(4# * degToRad(l0))
sin2m = Sin(2# * degToRad(m))
Etime = y * sin2l0 - 2# * e * sinm + 4# * e * y * sinm * cos2l0 _
- 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m
calcEquationOfTime = radToDeg(Etime) * 4#
End Function
Function calcHourAngleSunrise(lat, SolarDec)
'***********************************************************************/
'* Name:
calcHourAngleSunrise
'* Type:
Function
'* Purpose: calculate the hour angle of the sun at sunrise for the
'*
latitude
'* Arguments:
'*
lat : latitude of observer in degrees
'* solarDec : declination angle of sun in degrees
'* Return value:
'*
hour angle of sunrise in radians
'***********************************************************************/
Dim latrad As Double, sdRad As Double, HAarg As Double, HA As Double
latrad = degToRad(lat)
sdRad = degToRad(SolarDec)
HAarg = (Cos(degToRad(90.833)) / (Cos(latrad) * Cos(sdRad)) - Tan(latrad) * Tan(sdRad))
HA = (Application.WorksheetFunction.Acos(Cos(degToRad(90.833)) _
/ (Cos(latrad) * Cos(sdRad)) - Tan(latrad) * Tan(sdRad)))
calcHourAngleSunrise = HA
End Function
Function calcHourAngleSunset(lat, SolarDec)
'***********************************************************************/
'* Name:
calcHourAngleSunset
'* Type:
Function
'* Purpose: calculate the hour angle of the sun at sunset for the
'*
latitude
'* Arguments:
'*
lat : latitude of observer in degrees
'* solarDec : declination angle of sun in degrees
'* Return value:
'*
hour angle of sunset in radians
'***********************************************************************/
Dim latrad As Double, sdRad As Double, HAarg As Double, HA As Double
latrad = degToRad(lat)
sdRad = degToRad(SolarDec)
HAarg = (Cos(degToRad(90.833)) / (Cos(latrad) * Cos(sdRad)) - Tan(latrad) * Tan(sdRad))
HA = (Application.WorksheetFunction.Acos(Cos(degToRad(90.833)) _
/ (Cos(latrad) * Cos(sdRad)) - Tan(latrad) * Tan(sdRad)))
calcHourAngleSunset = -HA
End Function
Function calcSunriseUTC(jd, Latitude, longitude)
'***********************************************************************/
'* Name:
calcSunriseUTC
'* Type:
Function
'* Purpose: calculate the Universal Coordinated Time (UTC) of sunrise
'*
for the given day at the given location on earth
'* Arguments:
'*
JD : julian day
'*
latitude : latitude of observer in degrees
'*
longitude : longitude of observer in degrees
'* Return value:
'*
time in minutes from zero Z
'***********************************************************************/
Dim t As Double, eqtime As Double, SolarDec As Double, hourangle As Double
Dim delta As Double, timeDiff As Double, timeUTC As Double
Dim newt As Double
t = calcTimeJulianCent(jd)

QUAL2K

88

May 29, 2012

'

// *** First pass to approximate sunrise

eqtime = calcEquationOfTime(t)
SolarDec = calcSunDeclination(t)
hourangle = calcHourAngleSunrise(Latitude, SolarDec)

delta = longitude - radToDeg(hourangle)

timeDiff = 4 * delta
' in minutes of time
timeUTC = 720 + timeDiff - eqtime
' in minutes
' *** Second pass includes fractional jday in gamma calc
newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC / 1440#)
eqtime = calcEquationOfTime(newt)
SolarDec = calcSunDeclination(newt)
hourangle = calcHourAngleSunrise(Latitude, SolarDec)
delta = longitude - radToDeg(hourangle)
timeDiff = 4 * delta
timeUTC = 720 + timeDiff - eqtime
' in minutes
calcSunriseUTC = timeUTC
End Function
Function calcSolNoonUTC(t, longitude)
'***********************************************************************/
'* Name:
calcSolNoonUTC
'* Type:
Function
'* Purpose: calculate the Universal Coordinated Time (UTC) of solar
'*
noon for the given day at the given location on earth
'* Arguments:
'*
t : number of Julian centuries since J2000.0
'*
longitude : longitude of observer in degrees
'* Return value:
'*
time in minutes from zero Z
'***********************************************************************/
Dim newt As Double, eqtime As Double, solarNoonDec As Double, solNoonUTC As Double
newt = calcTimeJulianCent(calcJDFromJulianCent(t) + 0.5 + longitude / 360#)
eqtime = calcEquationOfTime(newt)
solarNoonDec = calcSunDeclination(newt)
solNoonUTC = 720 + (longitude * 4) - eqtime
calcSolNoonUTC = solNoonUTC
End Function
Function calcSunsetUTC(jd, Latitude, longitude)
'***********************************************************************/
'* Name:
calcSunsetUTC
'* Type:
Function
'* Purpose: calculate the Universal Coordinated Time (UTC) of sunset
'*
for the given day at the given location on earth
'* Arguments:
'*
JD : julian day
'*
latitude : latitude of observer in degrees
'*
longitude : longitude of observer in degrees
'* Return value:
'*
time in minutes from zero Z
'***********************************************************************/
Dim t As Double, eqtime As Double, SolarDec As Double, hourangle As Double
Dim delta As Double, timeDiff As Double, timeUTC As Double
Dim newt As Double
t = calcTimeJulianCent(jd)
'

// First calculates sunrise and approx length of day

eqtime = calcEquationOfTime(t)
SolarDec = calcSunDeclination(t)
hourangle = calcHourAngleSunset(Latitude, SolarDec)
delta = longitude - radToDeg(hourangle)
timeDiff = 4 * delta
timeUTC = 720 + timeDiff - eqtime

'

// first pass used to include fractional day in gamma calc

newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC / 1440#)
eqtime = calcEquationOfTime(newt)
SolarDec = calcSunDeclination(newt)
hourangle = calcHourAngleSunset(Latitude, SolarDec)

'

delta = longitude - radToDeg(hourangle)

timeDiff = 4 * delta
timeUTC = 720 + timeDiff - eqtime
// in minutes
calcSunsetUTC = timeUTC

End Function
Function sunrise(lat, lon, year, month, day, timezone, dlstime)
'***********************************************************************/
'* Name:
sunrise
'* Type:
Main Function called by spreadsheet
'* Purpose: calculate time of sunrise for the entered date
'*
and location.
'* For latitudes greater than 72 degrees N and S, calculations are
'* accurate to within 10 minutes. For latitudes less than +/- 72
'* accuracy is approximately one minute.
'* Arguments:
'
latitude = latitude (decimal degrees)
'
longitude = longitude (decimal degrees)
'
NOTE: longitude is negative for western hemisphere for input cells

QUAL2K

89

May 29, 2012

'
in the spreadsheet for calls to the functions named
'
sunrise, solarnoon, and sunset. Those functions convert the
'
longitude to positive for the western hemisphere for calls to
'
other functions using the original sign convention
'
from the NOAA javascript code.
'
year = year
'
month = month
'
day = day
'
timezone = time zone hours relative to GMT/UTC (hours)
'
dlstime = daylight savings time (0 = no, 1 = yes) (hours)
'* Return value:
'*
sunrise time in local time (days)
'***********************************************************************/
Dim longitude As Double, Latitude As Double, jd As Double
Dim riseTimeGMT As Double, riseTimeLST As Double
' change sign convention for longitude from negative to positive in western hemisphere
longitude = lon * -1
Latitude = lat
If (Latitude > 89.8) Then Latitude = 89.8
If (Latitude < -89.8) Then Latitude = -89.8
jd = calcJD(year, month, day)
'

// Calculate sunrise for this date

riseTimeGMT = calcSunriseUTC(jd, Latitude, longitude)

'

// adjust for time zone and daylight savings time in minutes

riseTimeLST = riseTimeGMT + (60 * timezone) + (dlstime * 60)

'

// convert to days
sunrise = riseTimeLST / 1440

End Function
Function solarnoon(lat, lon, year, month, day, timezone, dlstime)
'***********************************************************************/
'* Name:
solarnoon
'* Type:
Main Function called by spreadsheet
'* Purpose: calculate the Universal Coordinated Time (UTC) of solar
'*
noon for the given day at the given location on earth
'* Arguments:
'
year
'
month
'
day
'*
longitude : longitude of observer in degrees
'
NOTE: longitude is negative for western hemisphere for input cells
'
in the spreadsheet for calls to the functions named
'
sunrise, solarnoon, and sunset. Those functions convert the
'
longitude to positive for the western hemisphere for calls to
'
other functions using the original sign convention
'
from the NOAA javascript code.
'* Return value:
'*
time of solar noon in local time days
'***********************************************************************/
Dim longitude As Double, Latitude As Double, jd As Double
Dim t As Double, newt As Double, eqtime As Double
Dim solarNoonDec As Double, solNoonUTC As Double
' change sign convention for longitude from negative to positive in western hemisphere
longitude = lon * -1
Latitude = lat
If (Latitude > 89.8) Then Latitude = 89.8
If (Latitude < -89.8) Then Latitude = -89.8
jd = calcJD(year, month, day)
t = calcTimeJulianCent(jd)
newt = calcTimeJulianCent(calcJDFromJulianCent(t) + 0.5 + longitude / 360#)
eqtime = calcEquationOfTime(newt)
solarNoonDec = calcSunDeclination(newt)
solNoonUTC = 720 + (longitude * 4) - eqtime
'

// adjust for time zone and daylight savings time in minutes

solarnoon = solNoonUTC + (60 * timezone) + (dlstime * 60)

'

// convert to days
solarnoon = solarnoon / 1440

End Function
Function sunset(lat, lon, year, month, day, timezone, dlstime)
'***********************************************************************/
'* Name:
sunset
'* Type:
Main Function called by spreadsheet
'* Purpose: calculate time of sunrise and sunset for the entered date
'*
and location.
'* For latitudes greater than 72 degrees N and S, calculations are
'* accurate to within 10 minutes. For latitudes less than +/- 72
'* accuracy is approximately one minute.
'* Arguments:
'
latitude = latitude (decimal degrees)
'
longitude = longitude (decimal degrees)
'
NOTE: longitude is negative for western hemisphere for input cells
'
in the spreadsheet for calls to the functions named
'
sunrise, solarnoon, and sunset. Those functions convert the
'
longitude to positive for the western hemisphere for calls to
'
other functions using the original sign convention
'
from the NOAA javascript code.
'
year = year
'
month = month
'
day = day
'
timezone = time zone hours relative to GMT/UTC (hours)
'
dlstime = daylight savings time (0 = no, 1 = yes) (hours)
'* Return value:
'*
sunset time in local time (days)
'***********************************************************************/
Dim longitude As Double, Latitude As Double, jd As Double

QUAL2K

90

May 29, 2012

Dim setTimeGMT As Double, setTimeLST As Double

' change sign convention for longitude from negative to positive in western hemisphere
longitude = lon * -1
Latitude = lat
If (Latitude > 89.8) Then Latitude = 89.8
If (Latitude < -89.8) Then Latitude = -89.8
jd = calcJD(year, month, day)
'

// Calculate sunset for this date

setTimeGMT = calcSunsetUTC(jd, Latitude, longitude)

'

// adjust for time zone and daylight savings time in minutes

setTimeLST = setTimeGMT + (60 * timezone) + (dlstime * 60)

'

// convert to days
sunset = setTimeLST / 1440

End Function
Function solarazimuth(lat, lon, year, month, day, _
hours, minutes, seconds, timezone, dlstime)
'***********************************************************************/
'* Name:
solarazimuth
'* Type:
Main Function
'* Purpose: calculate solar azimuth (deg from north) for the entered
'*
date, time and location. Returns -999999 if darker than twilight
'*
'* Arguments:
'*
latitude, longitude, year, month, day, hour, minute, second,
'*
timezone, daylightsavingstime
'* Return value:
'*
solar azimuth in degrees from north
'*
'* Note: solarelevation and solarazimuth functions are identical
'*
and could be converted to a VBA subroutine that would return
'*
both values.
'*
'***********************************************************************/
Dim
Dim
Dim
Dim
Dim
Dim
Dim
Dim
Dim
Dim
Dim

longitude As Double, Latitude As Double

Zone As Double, daySavings As Double
hh As Double, mm As Double, ss As Double, timenow As Double
jd As Double, t As Double, r As Double
alpha As Double, theta As Double, Etime As Double, eqtime As Double
SolarDec As Double, earthRadVec As Double, solarTimeFix As Double
trueSolarTime As Double, hourangle As Double, harad As Double
csz As Double, zenith As Double, azDenom As Double, azRad As Double
azimuth As Double, exoatmElevation As Double
step1 As Double, step2 As Double, step3 As Double
refractionCorrection As Double, Te As Double, solarzen As Double
longitude = lon * -1
Latitude = lat
If (Latitude > 89.8) Then Latitude = 89.8
If (Latitude < -89.8) Then Latitude = -89.8
Zone = timezone * -1
daySavings = dlstime * 60
hh = hours - (daySavings / 60)
mm = minutes
ss = seconds

'//

timenow is GMT time for calculation in hours since 0Z

timenow = hh + mm / 60 + ss / 3600 + Zone
jd = calcJD(year, month, day)
t = calcTimeJulianCent(jd + timenow / 24#)
r = calcSunRadVector(t)
alpha = calcSunRtAscension(t)
theta = calcSunDeclination(t)
Etime = calcEquationOfTime(t)
eqtime = Etime
SolarDec = theta '//
earthRadVec = r

in degrees

solarTimeFix = eqtime - 4# * longitude + 60# * Zone

trueSolarTime = hh * 60# + mm + ss / 60# + solarTimeFix
'//
in minutes
Do While (trueSolarTime > 1440)
trueSolarTime = trueSolarTime - 1440
Loop
hourangle = trueSolarTime / 4# - 180#
'//
Thanks to Louis Schwarzmayr for the next line:
If (hourangle < -180) Then hourangle = hourangle + 360#
harad = degToRad(hourangle)
csz = Sin(degToRad(Latitude))
Sin(degToRad(SolarDec))
Cos(degToRad(Latitude))
Cos(degToRad(SolarDec))

*
+
*
*

_
_
_
Cos(harad)

If (csz > 1#) Then

csz = 1#
ElseIf (csz < -1#) Then
csz = -1#
End If
zenith = radToDeg(Application.WorksheetFunction.Acos(csz))
azDenom = (Cos(degToRad(Latitude)) * Sin(degToRad(zenith)))
If (Abs(azDenom) > 0.001) Then
azRad = ((Sin(degToRad(Latitude)) * _
Cos(degToRad(zenith))) - _
Sin(degToRad(SolarDec))) / azDenom
If (Abs(azRad) > 1#) Then
If (azRad < 0) Then

QUAL2K

91

May 29, 2012

azRad = -1#
Else
azRad = 1#
End If
End If
azimuth = 180# - radToDeg(Application.WorksheetFunction.Acos(azRad))
If (hourangle > 0#) Then
azimuth = -azimuth
End If
Else
If (Latitude > 0#) Then
azimuth = 180#
Else
azimuth = 0#
End If
End If
If (azimuth < 0#) Then
azimuth = azimuth + 360#
End If
exoatmElevation = 90# - zenith
If (exoatmElevation > 85#) Then
refractionCorrection = 0#
Else
Te = Tan(degToRad(exoatmElevation))
If (exoatmElevation > 5#) Then
refractionCorrection = 58.1 / Te - 0.07 / (Te * Te * Te) + _
0.000086 / (Te * Te * Te * Te * Te)
ElseIf (exoatmElevation > -0.575) Then
step1 = (-12.79 + exoatmElevation * 0.711)
step2 = (103.4 + exoatmElevation * (step1))
step3 = (-518.2 + exoatmElevation * (step2))
refractionCorrection = 1735# + exoatmElevation * (step3)
Else
refractionCorrection = -20.774 / Te
End If
refractionCorrection = refractionCorrection / 3600#
End If
solarzen = zenith - refractionCorrection
solarazimuth = azimuth
End Function
Function solarelevation(lat, lon, year, month, day, _
hours, minutes, seconds, timezone, dlstime)
'***********************************************************************/
'* Name:
solarazimuth
'* Type:
Main Function
'* Purpose: calculate solar azimuth (deg from north) for the entered
'*
date, time and location. Returns -999999 if darker than twilight
'*
'* Arguments:
'*
latitude, longitude, year, month, day, hour, minute, second,
'*
timezone, daylightsavingstime
'* Return value:
'*
solar azimuth in degrees from north
'*
'* Note: solarelevation and solarazimuth functions are identical
'*
and could converted to a VBA subroutine that would return
'*
both values.
'*
'***********************************************************************/
Dim
Dim
Dim
Dim
Dim
Dim
Dim
Dim
Dim
Dim
Dim

longitude As Double, Latitude As Double

Zone As Double, daySavings As Double
hh As Double, mm As Double, ss As Double, timenow As Double
jd As Double, t As Double, r As Double
alpha As Double, theta As Double, Etime As Double, eqtime As Double
SolarDec As Double, earthRadVec As Double, solarTimeFix As Double
trueSolarTime As Double, hourangle As Double, harad As Double
csz As Double, zenith As Double, azDenom As Double, azRad As Double
azimuth As Double, exoatmElevation As Double
step1 As Double, step2 As Double, step3 As Double
refractionCorrection As Double, Te As Double, solarzen As Double
longitude = lon * -1
Latitude = lat
If (Latitude > 89.8) Then Latitude = 89.8
If (Latitude < -89.8) Then Latitude = -89.8
Zone = timezone * -1
daySavings = dlstime * 60
hh = hours - (daySavings / 60)
mm = minutes
ss = seconds

'//

timenow is GMT time for calculation in hours since 0Z

timenow = hh + mm / 60 + ss / 3600 + Zone
jd = calcJD(year, month, day)
t = calcTimeJulianCent(jd + timenow / 24#)
r = calcSunRadVector(t)
alpha = calcSunRtAscension(t)
theta = calcSunDeclination(t)
Etime = calcEquationOfTime(t)
eqtime = Etime
SolarDec = theta '//
earthRadVec = r

in degrees

solarTimeFix = eqtime - 4# * longitude + 60# * Zone

trueSolarTime = hh * 60# + mm + ss / 60# + solarTimeFix
'//
in minutes
Do While (trueSolarTime > 1440)
trueSolarTime = trueSolarTime - 1440
Loop
hourangle = trueSolarTime / 4# - 180#

QUAL2K

92

May 29, 2012

'//
Thanks to Louis Schwarzmayr for the next line:
If (hourangle < -180) Then hourangle = hourangle + 360#
harad = degToRad(hourangle)
csz = Sin(degToRad(Latitude))
Sin(degToRad(SolarDec))
Cos(degToRad(Latitude))
Cos(degToRad(SolarDec))

*
+
*
*

_
_
_
Cos(harad)

If (csz > 1#) Then

csz = 1#
ElseIf (csz < -1#) Then
csz = -1#
End If
zenith = radToDeg(Application.WorksheetFunction.Acos(csz))
azDenom = (Cos(degToRad(Latitude)) * Sin(degToRad(zenith)))
If (Abs(azDenom) > 0.001) Then
azRad = ((Sin(degToRad(Latitude)) * _
Cos(degToRad(zenith))) - _
Sin(degToRad(SolarDec))) / azDenom
If (Abs(azRad) > 1#) Then
If (azRad < 0) Then
azRad = -1#
Else
azRad = 1#
End If
End If
azimuth = 180# - radToDeg(Application.WorksheetFunction.Acos(azRad))
If (hourangle > 0#) Then
azimuth = -azimuth
End If
Else
If (Latitude > 0#) Then
azimuth = 180#
Else
azimuth = 0#
End If
End If
If (azimuth < 0#) Then
azimuth = azimuth + 360#
End If
exoatmElevation = 90# - zenith
If (exoatmElevation > 85#) Then
refractionCorrection = 0#
Else
Te = Tan(degToRad(exoatmElevation))
If (exoatmElevation > 5#) Then
refractionCorrection = 58.1 / Te - 0.07 / (Te * Te * Te) + _
0.000086 / (Te * Te * Te * Te * Te)
ElseIf (exoatmElevation > -0.575) Then
step1 = (-12.79 + exoatmElevation * 0.711)
step2 = (103.4 + exoatmElevation * (step1))
step3 = (-518.2 + exoatmElevation * (step2))
refractionCorrection = 1735# + exoatmElevation * (step3)
Else
refractionCorrection = -20.774 / Te
End If
refractionCorrection = refractionCorrection / 3600#
End If
solarzen = zenith - refractionCorrection
solarelevation = 90# - solarzen
End Function
Sub solarposition(lat, lon, year, month, day, _
hours, minutes, seconds, timezone, dlstime, _
solarazimuth, solarelevation, earthRadVec)
'***********************************************************************/
'* Name:
solarposition
'* Type:
Subroutine
'* Purpose: calculate solar azimuth (deg from north)
'*
and elevation (deg from horizon) for the entered
'*
date, time and location.
'*
'* Arguments:
'*
latitude, longitude, year, month, day, hour, minute, second,
'*
timezone, daylightsavingstime
'* Return value:
'*
solar azimuth in degrees from north
'*
solar elevation in degrees from horizon
'*
earth radius vector (distance to the sun in AU)
'*
'* Note: solarelevation and solarazimuth functions are identical
'*
and could converted to a VBA subroutine that would return
'*
both values.
'*
'***********************************************************************/
Dim longitude As Double, Latitude As Double
Dim Zone As Double, daySavings As Double
Dim hh As Double, mm As Double, ss As Double, timenow As Double
Dim jd As Double, t As Double, r As Double
Dim alpha As Double, theta As Double, Etime As Double, eqtime As Double
'Dim SolarDec As Double, earthRadVec As Double, solarTimeFix As Double
Dim SolarDec As Double, solarTimeFix As Double
Dim trueSolarTime As Double, hourangle As Double, harad As Double
Dim csz As Double, zenith As Double, azDenom As Double, azRad As Double
Dim azimuth As Double, exoatmElevation As Double
Dim step1 As Double, step2 As Double, step3 As Double
Dim refractionCorrection As Double, Te As Double, solarzen As Double
longitude = lon * -1
Latitude = lat
If (Latitude > 89.8) Then Latitude = 89.8

QUAL2K

93

May 29, 2012

If (Latitude < -89.8) Then Latitude = -89.8

Zone = timezone * -1
daySavings = dlstime * 60
hh = hours - (daySavings / 60)
mm = minutes
ss = seconds
'//

timenow is GMT time for calculation in hours since 0Z

timenow = hh + mm / 60 + ss / 3600 + Zone
jd = calcJD(year, month, day)
t = calcTimeJulianCent(jd + timenow / 24#)
r = calcSunRadVector(t)
alpha = calcSunRtAscension(t)
theta = calcSunDeclination(t)
Etime = calcEquationOfTime(t)
eqtime = Etime
SolarDec = theta '//
earthRadVec = r

in degrees

solarTimeFix = eqtime - 4# * longitude + 60# * Zone

trueSolarTime = hh * 60# + mm + ss / 60# + solarTimeFix
'//
in minutes
Do While (trueSolarTime > 1440)
trueSolarTime = trueSolarTime - 1440
Loop
hourangle = trueSolarTime / 4# - 180#
'//
Thanks to Louis Schwarzmayr for the next line:
If (hourangle < -180) Then hourangle = hourangle + 360#
harad = degToRad(hourangle)
csz = Sin(degToRad(Latitude))
Sin(degToRad(SolarDec))
Cos(degToRad(Latitude))
Cos(degToRad(SolarDec))

*
+
*
*

_
_
_
Cos(harad)

If (csz > 1#) Then

csz = 1#
ElseIf (csz < -1#) Then
csz = -1#
End If
zenith = radToDeg(Application.WorksheetFunction.Acos(csz))
azDenom = (Cos(degToRad(Latitude)) * Sin(degToRad(zenith)))
If (Abs(azDenom) > 0.001) Then
azRad = ((Sin(degToRad(Latitude)) * _
Cos(degToRad(zenith))) - _
Sin(degToRad(SolarDec))) / azDenom
If (Abs(azRad) > 1#) Then
If (azRad < 0) Then
azRad = -1#
Else
azRad = 1#
End If
End If
azimuth = 180# - radToDeg(Application.WorksheetFunction.Acos(azRad))
If (hourangle > 0#) Then
azimuth = -azimuth
End If
Else
If (Latitude > 0#) Then
azimuth = 180#
Else
azimuth = 0#
End If
End If
If (azimuth < 0#) Then
azimuth = azimuth + 360#
End If
exoatmElevation = 90# - zenith
If (exoatmElevation > 85#) Then
refractionCorrection = 0#
Else
Te = Tan(degToRad(exoatmElevation))
If (exoatmElevation > 5#) Then
refractionCorrection = 58.1 / Te - 0.07 / (Te * Te * Te) + _
0.000086 / (Te * Te * Te * Te * Te)
ElseIf (exoatmElevation > -0.575) Then
step1 = (-12.79 + exoatmElevation * 0.711)
step2 = (103.4 + exoatmElevation * (step1))
step3 = (-518.2 + exoatmElevation * (step2))
refractionCorrection = 1735# + exoatmElevation * (step3)
Else
refractionCorrection = -20.774 / Te
End If
refractionCorrection = refractionCorrection / 3600#
End If
solarzen = zenith - refractionCorrection
solarazimuth = azimuth
solarelevation = 90# - solarzen
End Sub

QUAL2K

94

May 29, 2012

APPENDIX C: SEDIMENT-WATER HEAT EXCHANGE

Although the omission of sediment-water heat exchange is usually justified for deeper systems, it
can have a significant impact on the heat balance for shallower streams. Consequently, sedimentwater heat exchange is included in QUAL2K.
A major impediment to its inclusion is that incorporating sediment heat transfer often carries
a heavy computational burden. This is because the sediments are usually represented as a
vertically segmented distributed system. Thus, inclusion of the mechanism results in the addition
of numerous sediment segments for each overlying water element.
In the present appendix, I derive a computationally-efficient lumped approach that yields
comparable results to the distributed methods.
The conduction equation is typically used to simulate the vertical temperature distribution in
a distributed sediment (Figure 26a)
T
2T
2
t
x

(254)

This model can be subjected to the following boundary conditions:

T (0, t ) T Ta cos (t )
T (, t ) T

where T = sediment temperature [oC], t = time [s], = sediment thermal diffusivity [m2 s1], and z
= depth into the sediments [m], where z = 0 at the sediment-water interface and z increases
downward, T = mean temperature of overlying water [oC], Ta = amplitude of temperature of
overlying water [oC], = frequency [s1] = 2/Tp, Tp = period [s], and = phase lag [s]. The first
boundary condition specifies a sinusoidal Dirichlet boundary condition at the sediment-water
interface. The second specifies a constant temperature at infinite depth. Note that the mean of the
surface sinusoid and the lower fixed temperature are identical.
0

z
(a) distributed

(b) Lumped

Figure 26. Alternate representations of sediments: (a) distributed and (b) lumped.

Applying these boundary conditions, Eq. (254) can be solved for (Carslaw and Jaeger 1959)
QUAL2K

95

May 29, 2012

T ( z , t ) T Ta e ' z cos (t ) ' z

(255)

where [m1] is defined as

'

(256)

The heat flux at the sediment water interface can then be determined by substituting the
derivative of Eq. (255) into Fouriers law and evaluating the result at the sediment-water interface
(z = 0) to yield
J (0, t ) C p Ta cos (t ) / 4

(257)

where J(0, t) = flux [W/m2].

An alternative approach can be developed using a first-order lumped model (Figure 26b),
H s s C ps

dTs s s C ps
T Ta cos (t ) Ts

dt
Hs / 2

where Hsed = the thickness of the sediment layer [m], s = sediment density [kg/m3], and Cps =
sediment specific heat [joule (kg oC)]1]. Collecting terms gives,
dT
khT khT khTa cos (t )
dt

where kh = 2s/Hsed2. After initial transient have died out, this solution to this equation is
kh

T T

kh 2

Ta cos (t ) tan 1 ( / kh )

(258)

which can be used to determine the flux as

J

kh
2
C pTa cos (t )
cos (t ) tan 1 ( / kh )
H sed

kh 2 2

(259)

It can be shown that Eqs. (231) and (234) yield identical results if the depth of the single layer
is set at
H sed

(260)

'

Water quality models typically consider annual, weekly and diel variations. Using = 0.0035
cm2/s (Hutchinson 1957), the single-layer depth that would capture these frequencies can be
calculated as 2.2 m, 30 cm and 12 cm, respectively.

QUAL2K

96

May 29, 2012

Because QUAL2K resolves diel variations, a value on the order of 12 cm should be selected
for the sediment thickness. We have chosen of value of 10 cm as being an adequate first estimate
because of the uncertainties of the river sediment thermal properties (Table 4).

QUAL2K

97

May 29, 2012