DOC

# Q2KDocv2_11b8.doc

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

1 SEGMENTATION AND HYDRAULICS

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

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

1

Point source2Point withdrawal

Point withdrawal3

4Point source

5Non-point

withdrawal6

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

1818

1919Trib3Main stem

HW#4222020

23212124

252626

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

ElementsElementsReachReach

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.

QQin,iout,i

QQi1i

i1ii+ 1

Figure 4 Element flow balance.

The total inflow from sources is computed as

psinpsi

(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

painpai

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

Qnpt

25%25%50%

112

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.

Weirs

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

Bw

HhHdHiHiHHww

Hi+1

elev2elev2ii

elev1i+1elev1i+1

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

VBHx

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

S0B1

H11

sss2s1Q,UB0

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

BB(ss)H10s1s2

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

ABxs1

VBHx

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

MATERIAL n

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.

HiHd

elev2iHi+1

elev1i+1

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,

j

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

follows:

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

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

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

QUAL2K 7 December 16, 2008

2 TEMPERATURE MODEL

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

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

''QdTQQEEout,iii1ii1iTTTTTTT；？；？i1iii1ii1idtVVVVViiiii

(26)

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

inflowoutflowidispersiondispersion

sediment-water

transfer

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)

psinpsi??

??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

air-water

interface

solaratmosphericwaterconductionevaporation

shortwavelongwavelongwaveandand

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

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

topography).

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 local hour angle in radians is given by

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

where:

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 (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

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 (http://www.atdd.noaa.gov/isis.htm). The selection of either the Bras or Ryan-Stolzenbach solar radiation model and the appropriate atmospheric turbidity factor or atmospheric transmission coefficient for a particular application should ideally be guided by a comparison of predicted solar radiation with measured values at a reference location. Cloud Attenuation. Attenuation of solar radiation due to cloud cover is computed with

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
cust-service@docsford.com