By Christine Reed,2014-02-14 17:03
13 views 0


    The model represents a river as a series of reaches. These represent stretches of river that have constant hydraulic characteristics (e.g., slope, bottom width, etc.). As depicted in Figure 1, the reaches are numbered in ascending order starting from the headwater of the river‟s main stem.

    Notice that both point and non-point sources and point and non-point withdrawals (abstractions) can be positioned anywhere along the channel‟s length.

    Headwater boundary


    Point source2Point withdrawal

    Point withdrawal3

    4Point source



    Non-point7source Figure 1 QUAL2K segmentation scheme for a river with no tributaries. Point source8For systems with tributaries (Figure 2), the reaches are numbered in ascending order starting at reach 1 at the headwater of the main stem. When a junction with a tributary is reached, the Downstream boundaryHW#1numbering continues at that tributary‟s headwater. Observe that both the headwaters and the

    tributaries are also numbered consecutively following a sequencing scheme similar to the reaches. 11HW#2Note also that the major branches of the system (that is, the main stem and each of the tributaries) 6622are referred to as segments. This distinction has practical importance because the software

    77provides plots of model output on a segment basis. That is, the software generates individual plots 3388for the main stem as well as each of the tributaries. 12124413131111551414Trib110101515991616HW#3Trib21717


    1919Trib3Main stem




    2727 2828Figure 2 QUAL2K segmentation scheme for (a) a river with tributaries. The Q2K reach

    representation in (b) illustrates the reach, headwater and tributary numbering 2929schemes.

    Finally, any model reach can be further divided into a series of equally-spaced elements. As in Figure 3, this is done by merely specifying the number of elements that are desired. (a) A river with tributaries(b) Q2K reach representation

QUAL2K 1 December 16, 2008

    nn= 4= 4


    Figure 3 If desired, any model reach can be further subdivided into a series of n equal-

    length elements.

    In summary, the nomenclature used to describe the way in which Q2K organizes river

    topology is as follows:

    ? Reach. A length of river with constant hydraulic characteristics.

    ? Element. The model‟s fundamental computational unit which consists of an equal length

    subdivision of a reach.

    ? Segment. A collection of reaches representing a branch of the system. These consist of

    the main stem as well as each tributary.

    ? Headwater. The upper boundary of a model segment. 1.1 Flow Balance

    As described in the last section, Q2K‟s most fundamental unit is the element. A steady-state flow

    balance is implemented for each model element as (Figure 4)

    QQQQ (1) ii1in,iout,i3where Q = outflow from element i into the downstream element i + 1 [m/d], Q= inflow from ii1 3the upstream element i 1 [m/d], Q is the total inflow into the element from point and in,i3nonpoint sources [m/d], and Q is the total outflow from the element due to point and nonpoint out,i3withdrawals [m/d]. Thus, the downstream outflow is simply the difference between inflow and

    source gains minus withdrawal losses.



    i1ii+ 1

    Figure 4 Element flow balance.

    The total inflow from sources is computed as


     (2) QQQ,,,,,??inipsijnpsij11jj3where Q is the jth point source inflow to element i [m/d], psi = the total number of point ps,i,j3sources to element i, Q is the jth non-point source inflow to element i [m/d], and npsi = the nps,i,j

    total number of non-point source inflows to element i.

    The total outflow from withdrawals is computed as


     (3) QQQut,,,,,??oipaijnpaij11jj

    QUAL2K 2 December 16, 2008

    3where Q is the jth point withdrawal outflow from element i [m/d], pai = the total number of pa,i,j

    point withdrawals from element i, Q is the jth non-point withdrawal outflow from element i npa,i,j3[m/d], and npai = the total number of non-point withdrawal flows from element i.

    The non-point sources and withdrawals are modeled as line sources. As in Figure 5, the non-

    point source or withdrawal is demarcated by its starting and ending kilometer points. Its flow is

    then distributed to or from each element in a length-weighted fashion.




    startend Figure 5 The manner in which non-point source flow is distributed to an element. 1.2 Hydraulic Characteristics

    Once the outflow for each element is computed, the depth and velocity are calculated in one of

    three ways: weirs, rating curves, and Manning equations. The program decides among these options in the following manner:

    ? If weir height and width are entered, the weir option is implemented.

    ? If the weir height and width are zero and rating curve coefficients are entered (a and ),

    the rating curve option is implemented.

    ? If neither of the previous conditions is met, Q2K uses the Manning equation.


    Figure 6 shows how weirs are represented in Q2K. Note that a weir can only occur at the end of a reach consisting of a single element. The symbols shown in Figure 6 are defined as: H = i

    the depth of the element upstream of the weir [m], H = the depth of the element downstream of i+1

    the weir [m], elev2 = the elevation above sea level of the tail end of the upstream element [m], i

    elev1 = the elevation above sea level of the head end of the downstream element [m], H = the i+1w

    height of the weir above elev2 [m], H = the drop between the elevation above sea level of the id

    surface of element i and element i+1 [m], H = the head above the weir [m], B = the width of the hw

    weir [m]. Note that the width of the weir can differ from the width of the element, B. i

    (a) Side(b) Cross-section






    Figure 6 A sharp-crested weir occurring at the boundary between two reaches. QUAL2K 3 December 16, 2008

    For a sharp-crested weir where H/H < 0.4, flow is related to head by (Finnemore and hw

    Franzini 2002)

    3/2Q1.83BH (4) iwh3where Q is the outflow from the element upstream of the weir in m/s, and B and H are in m. iwh

    Equation (4) can be solved for

    2/3??Qi?? (5) Hh??1.83Bw??

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

     (6) HHHiwh

    and the drop over the weir

    Helev2Helev1H (7) diii1i1

    Note that this drop is used to compute oxygen and carbon dioxide gas transfer due to the weir (see

    pages 43 and 49).

    The cross-sectional area, velocity, surface area and volume of element i can then be

    computed as

    ABH (8) c,iii

    Qi (9) UiAc,i

    ABx s,iii

    VBHx iiii

    where B = the width of element i, x = the length of element i. Note that for reaches with weirs, ii

    the reach width must be entered. This value is entered in the column AA (labeled "Bottom

    Width") of the Reach Worksheet.

    1.2.1 Rating Curves

    Power equations (sometimes called Leopold-Maddox relationships) can be used to relate mean

    velocity and depth to flow for the elements in a reach,

    bUaQ (10)

    ~HQ (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 (12) AcU

    Ac (13) BH

    The surface area and volume of the element can then be computed as

    ABx s


    The exponents b and ~ typically take on values listed in Table 1. Note that the sum of b and ~

    must be less than or equal to 1. If this is not the case, the width will decrease with increasing flow.

    If their sum equals 1, the channel is rectangular.

    Table 1 Typical values for the exponents of rating curves used to determine velocity and

    depth from flow (Barnwell et al. 1989).

    Equation Exponent Typical value Range

    bb 0.43 0.40.6 UaQ

    QUAL2K 4 December 16, 2008

    ~0.45 0.30.5 ~;HQ

    In some applications, you might want to specify constant values of depth and velocity that do not vary with flow. This can be done by setting the exponents b and ~ to zero and setting a equal

    to the desired velocity and equal to the desired depth.

    1.2.2 Manning Equation

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

    1/25/3SA0cQ (14) 2/3nP31where Q = flow [m/s], S = bottom slope [m/m], n = the Manning roughness coefficient, A = 0c2the cross-sectional area [m], and P = the wetted perimeter [m].




    Figure 7 Trapezoidal channel.

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

    ,:AB0.5(ss)HH (15) c0s1s2

    where B = bottom width [m], s and s = the two side slopes as shown in Figure 7 [m/m], and H 0s1s2

    = element depth [m].

    The wetted perimeter is computed as

    22 (16) PBHs1Hs1ss012

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

    2/53/522??(Qn)BHs1Hs1??0k1s1k1s2??H (17) k3/10,:SB0.5(ss)H0s1s2k1

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

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

    HHk1k??100% (18) aHk1

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

    QU (19) Ac

     1 Notice that time is measured in seconds in this and other formulas used to characterize hydraulics. This is how the computations are implemented within Q2K. However, once the hydraulic characteristics are determined they are converted to day units to be compatible with other computations. QUAL2K 5 December 16, 2008

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

    Ac (20) BH

    The top width, B [m], is computed as 1


    The surface area and volume of the element can then be computed as



    Suggested values for the Manning coefficient are listed in Table 2. Manning‟s n typically

    varies with flow and depth (Gordon et al. 1992). As the depth decreases at low flow, the relative roughness usually increases. Typical published values of Manning‟s n, which range from about 0.015 for smooth channels to about 0.15 for rough natural channels, are representative of conditions when the flow is at the bankfull capacity (Rosgen, 1996). Critical conditions of depth for evaluating water quality are generally much less than bankfull depth, and the relative roughness may be much higher.

    Table 2 The Manning roughness coefficient for various open channel surfaces (from Chow

    et al. 1988).


    Man-made channels

    Concrete 0.012

    Gravel bottom with sides:

    Concrete 0.020

    mortared stone 0.023

    Riprap 0.033

    Natural stream channels

    Clean, straight 0.025-0.04

    Clean, winding and some weeds 0.03-0.05

    Weeds and pools, winding 0.05

    Mountain streams with boulders 0.04-0.10

    Heavy brush, timber 0.05-0.20

    1.2.3 Waterfalls

    In Section 0, the drop of water over a weir was computed. This value is needed in order to compute the enhanced reaeration that occurs in such cases. In addition to weirs, such drops can also occur at waterfalls (Figure 8). Note that waterfalls can only occur at the end of a reach.




    Figure 8 A waterfall occurring at the boundary between two reaches.

    QUAL2K computes such drops for cases where the elevation above sea level drops abruptly at the boundary between two reaches. Equation (7) is used to compute the drop. It should be noted that the drop is only calculated when the elevation above sea level at the downstream end of a reach is greater than at the beginning of the next downstream reach; that is, elev2 > elev1. ii+1

QUAL2K 6 December 16, 2008

1.3 Travel Time

    The residence time of each element is computed as

    Vk (21) kQkthth3where = the residence time of the k element [d], V = the volume of the k element [m] = kkth2thAx, A = the cross-sectional area of the k element [m], and x = the length of the k c,kkc,kk

    element [m]. These times are then accumulated to determine the travel time along each of the river‟s segments (that is, either the main stem or one of the tributaries). For example, the travel time from the headwater to the downstream end of the jth element in a segment is computed as,


    t (22) tjk,?k1

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

    1.4 Longitudinal Dispersion

    Two options are used to determine the longitudinal dispersion for a boundary between two elements. First, the user can simply enter estimated values on the Reach Worksheet. If the user

    does not enter values, a formula is employed to internally compute dispersion based on the channel‟s hydraulics (Fischer et al. 1979),

    22UBiiE0.011 (23) p,i*HUii

    2where E = the longitudinal dispersion between elements i and i + 1 [m/s], U = velocity [m/s], p,ii*B = width [m], H = mean depth [m], and U = shear velocity [m/s], which is related to more iii

    fundamental characteristics by

    * (24) UgHSiii2where g = acceleration due to gravity [= 9.81 m/s] and S = channel slope [dimensionless].

    After computing or prescribing E, the numerical dispersion is computed as p,i

    Uxii (25) En,i2

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


    ? If E ? E, the model dispersion, E is set to E E. n,ip,iip,in,i

    ? If E > E, the model dispersion is set to zero. n,ip,i

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

    QUAL2K 7 December 16, 2008


    As in Figure 9, the heat balance takes into account heat transfers from adjacent elements, loads,

    withdrawals, the atmosphere, and the sediments. A heat balance can be written for element i as



    3WJJ??mmm????h,ia,is,i?? ????63??CVCH100 cmCH100 cm10 cm??????wpwiwpwiwpwiowhere T = temperature in element i [C], t = time [d], E = the bulk dispersion coefficient ii3between elements i and i + 1 [m/d], W = the net heat load from point and non-point sources into h,i3oelement i [cal/d], = the density of water [g/cm], C = the specific heat of water [cal/(g C)], wpw22J = the air-water heat flux [cal/(cm d)], and J = the sediment-water heat flux [cal/(cm d)]. a,is,iatmospheric

    transferheat loadheat withdrawal




    sediment Figure 9 Heat balance for an element.

    The bulk dispersion coefficient is computed as

    EAic,i'E (27) i;?xx/2ii1

    Note that two types of boundary condition are used at the river‟s downstream terminus: (1) a zero dispersion condition (natural boundary condition) and (2) a prescribed downstream boundary

    condition (Dirichlet boundary condition). The choice between these options is made on the

    Downstream Worksheet.

    The net heat load from sources is computed as (recall Eq. 2)


    ??WCQTQT (28) h,ipps,i,jpsi,jnps,i,jnpsi,j????j1j1??owhere T is the temperature of the jth point source for element i [C], and T is the ps,i,jnps,i,jotemperature of the jth non-point source temperature for element i [C].

    2.1 Surface Heat Flux

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

    J I(0)JJJJ (29) hanbrce

    where I(0) = net solar shortwave radiation at the water surface, J = net atmospheric longwave an

    radiation, J = longwave back radiation from the water, J = conduction, and J = evaporation. brce2All fluxes are expressed as cal/cm/d.

QUAL2K 8 December 16, 2008

    non-radiation termsradiation terms






    net absorbed radiationwater-dependent terms Figure 10 The components of surface heat exchange.

    2.1.1 Solar Radiation

    The model computes the amount of solar radiation entering the water at a particular latitude (L) at

    and longitude (L) on the earth‟s surface. This quantity is a function of the radiation at the top of lm

    the earth‟s atmosphere which is attenuated by atmospheric transmission, cloud cover, reflection, and shade,

    (0) (1) (1)IIaaRS0tcsf

     (30) extraterrestrial atmospheric cloud reflection shading

     radiation attenuation attenuation2where I(0) = solar radiation at the water surface [cal/cm/d], I = extraterrestrial radiation (i.e., at 02the top of the earth‟s atmosphere) [cal/cm/d], a = atmospheric attenuation, a = cloud attenuation, tc

    R = albedo (fraction reflected), and S = effective shade (fraction blocked by vegetation and sf


    Extraterrestrial radiation. The extraterrestrial radiation is computed as (TVA 1972)

    W0Isin (31) 02r22where W = the solar constant [1367 W/m or 2823 cal/cm/d], r = normalized radius of the 0

    earth‟s orbit (i.e., the ratio of actual earth-sun distance to mean earth-sun distance), and = the sun‟s altitude [radians], which can be computed as

    ;?sinsin?sinLcos?cosLcos (32) atat

    where ? = solar declination [radians], L = local latitude [radians], and = the local hour angle of at

    the sun [radians].

    The local hour angle in radians is given by

    trueSolarTime)??180 (33) ??4180??


    trueSolarTimelocalTimeeqtime4?L60?timezone (34) lm

    where trueSolarTime is the solar time determined from the actual position of the sun in the sky

    [minutes], localTime is the local time in minutes (local standard time), L is the local longitude lm

    (positive decimal degrees for the western hemisphere), and timezone is the local time zone in hours relative to Greenwich Mean Time (e.g. 8 hours for Pacific Standard Time; the local time

    zone is selected on the QUAL2K Worksheet). The value of eqtime represents the difference between true solar time and mean solar time in minutes of time.

QUAL2K 9 December 16, 2008

    QUAL2K calculates the solar declination, hour angle, solar altitude, and normalized radius (distance between the earth and sun), as well as the times of sunrise and sunset using the Meeus (1999) algorithms as implemented by NOAA‟s Surface Radiation Research Branch ( 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

    ftt (35) sssr

    where t = time of sunset [hours] and t = time of sunrise [hours]. sssr

    Atmospheric attenuation. Various methods have been published to estimate the fraction of the atmospheric attenuation from a clear sky (a). Two alternative methods are available in QUAL2K t

    to estimate a (Note that the solar radiation model is selected on the Light and Heat Worksheet t

    of QUAL2K):

    1) Bras (default)

    The Bras (1990) method computes a as: t

    namfac1 (36) aet

    where n is an atmospheric turbidity factor that varies from approximately 2 for clear skies to 4 fac

    or 5 for smoggy urban areas. The molecular scattering coefficient (a) is calculated as 1

    a0.1280.054logm (37) 110

    where m is the optical air mass, calculated as

    1m (38) 1.253sin0.15(3.885)dowhere is the sun‟s altitude in degrees from the horizon = ? (180/)). d

    2) Ryan and Stolzenbach

    The Ryan and Stolzenbach (1972) model computes a from ground surface elevation and solar t

    altitude as: 5.2562880.0065elev??m??aa288 (39) ??ttc

    where a is the atmospheric transmission coefficient (0.70-0.91, typically approximately 0.8), and tc

    elev is the ground surface elevation in meters.

    Direct measurements of solar radiation are available at some locations. For example, NOAA‟s Integrated Surface Irradiance Study (ISIS) has data from various stations across the United States ( The selection of either the Bras or Ryan-Stolzenbach solar radiation model and the appropriate atmospheric turbidity factor or atmospheric transmission coefficient for a particular application should ideally be guided by a comparison of predicted solar radiation with measured values at a reference location. Cloud Attenuation. Attenuation of solar radiation due to cloud cover is computed with

    2a10.65C (40) cL

    where C = fraction of the sky covered with clouds. L

    Reflectivity. Reflectivity is calculated as

    B (41) RAsd

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

    Table 3 Coefficients used to calculate reflectivity based on cloud cover. Cloudiness Clear Scattered Broken Overcast

    C 0 0.1-0.5 0.5-0.9 1.0 L

    A B A B A B A B Coefficients

     1.18 0.77 2.20 0.97 0.95 0.75 0.35 0.45

    QUAL2K 10 December 16, 2008

Report this document

For any questions or suggestions please email