0.49For values of the discharge up to the bankfhll discharge Qbp Rh and the mean boundary shearstress t both increase with increasing Q. At Qbf the mean bank and bed shear stresses aregreater than for any lower Q values. This is important from bank stability considerationsbecause if the sediment which forms the channel banks is not being eroded at Qbf it will bestable at all lower discharges.A more complex case exits for flows in excess of Qbf The flow is spread over a wide area andthere are complex interactions between the overbank flow on the floodplain, and the channelflow. This floodplain-channel interaction was first demonstrated experimentally by Sellin (1964;although see Cruff (1965) for earlier unpublished reports) who identified vortices with verticalaxes at the interface between the channel and floodplain flows. These vortices result inmomentum exchange between the relatively fast channel flow, and the slower floodplain flowand produce additional flow resistance effects.The simplest case for analysis is the channel with an infinitely wide floodplain. The excess flowabove bankfi.ill will be spread across an infinitely wide area, and therefore the flow depth on thefloodplain is essentially zero, while the depth within the channel does not increase beyond thebankfull depth. For discharges in excess of the Qbp value of f3 remains equal to zero, and thefloodplain-channel interaction can be assumed to be negligible. The bank and bed shear stresseswhich developed at Qbf therefore represent maximum values which are not exceeded at higherdischarges. A channel which has stable banks at Qbf will be stable at all other discharges, andthe sediment transport rate also reaches the maximum value atFor channels with floodplains of finite width, when Q exceeds values of f3 > 0 can develop.Under this condition complex three dimensional flow patterns can develop due to interactionsbetween the floodplain and main channel flow. These complex flow patterns have been thefocus of a considerable research effort. Several investigators have demonstrated the influence of50the floodplain flow on the main channel hydraulics by constructing physical models, measuringthe hydraulic parameters of the main channel with floodplain flow, and then isolating thefloodplain from the main channel by smooth impervious barriers such as glass sheets. Anexample is shown in Fig 3.6 from Sellin (1964). The isovels for the channel with floodplainflow indicate the mean velocity and discharge within the main channel are reduced verysignificantly when compared to the same channel with the flood plain isolated from the mainchannel. From his experiments Sellin (1964) determined that the velocity and discharge in themain channel are reduced by approximately 30 %.It can be observed from the velocity gradients in Fig 3.6 that the bed and bank shear stressesare also reduced significantly due to the floodplain-channel interaction. Barishnikov (1967)performed experiments similar to those of Sellin (1964) to determine the effect of thefloodplain-channel interaction on the bedload transport capacity of the channel and determinedthat the bedload capacity of the main channel could be reduced by 75 - 80 % due to overbankflow interactions. Similarly Meyers and Elsawy (1975) investigated the effect of the floodplain-channel interaction on the boundary shear stress and found that the mean and maximum valuesof the bed and bank shear stresses could be reduced by 20 - 30% of the bankfiill values forsmall positive values of3.In general, it has been determined that the effects of the floodplain-channel interaction aregreatest when the difference between the flow velocities on the floodplain and channel isgreatest. The interaction effects increase with increasing floodplain roughness (Barishnikov,1967; Knight and Hamed, 1984), for larger ratios of floodplain width to channel width (Knightand Demetriou, 1983; Holden and James, 1989), and for small non-zero values of (3 (Knightand Demitriou, 1983). Several investigators have determined that the floodplain-channelinteraction effects on the main channel flow become negligible for values of (3 greater thanabout 0.3 - 0.5 (Posey, 1967; Meyers and Elsawy, 1975; Bhowmik and Demisse, 1982; Pascheand Rouve, 1985).51In conclusion it is evident that for values of 13 between 0 and about 0.5, the bed and bank shearstresses within the main channel are in general somewhat less than the values at Qbp and that achannel which is stable at the banlcfull stage will remain stable for larger discharges.Furthermore the sediment transporting capacity for discharges above bankfl.ill may be nogreater, or even less than at Qbf This conclusion does not necessarily apply to very large floodevents where 13 is greater than about 0.5 and the floodplain-channel interactions becomenegligible.In the present model formulation the infinite floodplain model will be assumed where (3 remainsequal to zero for discharges equal to and exceeding Qbf The values of the bed and bank shearstresses will be assumed to be equal to the bankfull values for discharges in excess of bankfiill.The stability of the channel will be assessed at Qbp and if the channel banks are stable at Qbfthey will be assumed to be stable with respect to fluvial erosion of the bank sediment for allother discharges. Furthermore the value of the bed shear stress calculated at Qbf will be used tomodel the sediment transporting capacity of the channel for flows which exceed Qbf3.3.1 Recurrence Interval of Bankfull FlowSeveral investigators in the past have concluded that Qbf corresponds to a characteristic returnperiod. Wolman and Leopold (1957) and Leopold et al. (1964) that an average recurrenceinterval of 1.5 years based on the annual maximum series is a good estimate of Qbf SimilarlyDury (1973) suggests a value of 1.58 years. Nixon (1959) based his analysis on the partialduration series and found a mean recurrence interval of about 0.5 years. According toHenderson (1966, p. 465) the result of Wolman and Leopold (1957) is equivalent to arecurrence interval of 0.9 years based on the partial duration series.Williams (1978) completed an analysis of 36 stations to compute the return period of Qbf forthe active floodplain and found a mean recurrence interval of 0.9 and 1.5 years based upon the52partial duration and annual maximum series, respectively. A plot of the frequency distributionfor the analysis ofWilliams (1978) is shown in Fig 3.7. Disregard the valley flat curve as Qbfhasbeen defined herein as the height of the active floodplain. While the mean recurrence values inFig 3.7 agree closely with the 1.5 year average of Wolman and Leopold (1957) and Leopold etaL (1964), there is a wide spread of values between 1.01 to 32 years based on the annualmaximum series (0.25 to 32 years for the partial duration series), and the standard errorassociated with this mean value is 0.277 log units. Williams (1978) concluded that because ofthe wide range in recurrence intervals an average recurrence interval has little meaning and is apoor estimate of QbfThe suggestion of a characteristic recurrence interval has resulted in the widespread assumptionthat Qbf is a primary feature of the watershed, and can be treated as an independent variablewhereby the value of Qbf is imposed on the channel. This has been a prime assumption inseveral analytical studies of river channel development (Parker, 1978; Chang, 1979, 1980;White et a!., 1982; Millar and Quick, 1 993b) and empirical regime studies of rivers (Charlton eta!., 1978; Andrews, 1984; and Hey and Thorne, 1986). However if one agrees with theconclusion of Williams (1978) that there is in fact no characteristic recurrence interval forthe assumption of Qbf as an independent variable becomes somewhat tenuous. Furthermore theactual value of Qbf is defined by the width, depth, slope, and roughness of the channel, all ofwhich are dependent variables. The optimisation model will be used in Chapter 7 to show thatthe value of Qbf can be viewed as a dependent variable.3.4 SUBDiVISION OFfINTO BED AND BANK COMPONENTSThe preceding discussion requires the assumption that the values off for the bed and banksections of the channel are equal. Alternatively the wide channel approximation may beassumed whereby the channel roughness is equal to the bed roughness only, and thecontribution of the bank sections to! is negligible. The distribution of the shear force (per unit53channel length) which is defined as the product of shear stress and wetted perimeter can bedivided into bed and bank components:bed.d bnk’3a (3.18)where P is the wetted perimeter, and the subscripts bed and bank indicate the bed and bankcomponents of the total (Fig 3.5b). Eqn (3.18) forms the basis for the empirical methoddeveloped by Knight (1981) and Knight eta!. (1984) for calculating the mean values of the bedand bank shear stresses, and is used in Chapters 4 and 5 to calculate the bedload transportcapacity and to assess the bank stability.By using Eqn (3.5) to express t in terms ofJ dividing throughout by P and cancelling liketerms, Eqn (3.18) can be expressed in terms offP Pff bedf bankJ jbed ibankwherefbed is the total bed friction factor, and fbaflk is the total bank friction factor.For the wide channels typical of natural rivers bed P, and bank <1.59G[c1]=exp[14.2(c1_ 1) — 9.28(— 1)2] 1 1.59 (4.32)cF’42I-Equations (4.30) to (4.32) will now be generalised to apply to a bed comprised of gradedsediment. A key assumption is that the mobility of the armour layer is determined by themobility of the geometric mean grain diameter D. The function G is calculated as a function ofD. Further modifications are necessary to account for the hiding effects in a graded armour86layer, and the relation is generalised to a bed of arbitrary composition by an order-one strainingfunction. The value of W is given by:= 0.00218G[cI g0(o) c] (4.33)where the function G is unchanged from (4.32), 1 is calculated from Eqn (4.31) with DDsg divided by a dimensionless reference shear stress = 0.0386. The functiong0(4) is areduced hiding function, and a is a straining function. These two modifying functions will bediscussed below.The reduced hiding function accounts for the hiding of smaller grains and the over exposure oflarger grains relative to Dsg. The function 4. is a modification ofEqn. (4.25) and is defined as:D.Sf =_—‘ (4.34)The value ofDsg is given by:Dsg = expi1. lnDjJ (4.35)The function g0(4) was developed from Oak Creek data and was assumed to have the valueg0(1)1:go(S.) = (4.36)For grains smaller than Dsg the value ofg0(4) is less than one and the effective shear stressacting on the grain D. is reduced. The converse is true for grains larger than Dsg.87The order-one straining function a was introduced based upon an analysis of Oak Creek datato account for the grain size distribution of the armour layer sediment. From a numericalanalysis using Oak Creek data the function awas found to vary from an asymptotic value closeto 1 for low values ofc13g0 (and low values ofg) for which a coarse, well sorted equilibriumsurface prevails, to 0.453 at high values of where a finer and more poorly sorted surfacelayer is predicted. From this result it was postulated that a generalised straining function comight be a function of and some measure of the armour grain size dispersion. The measureof the grain size dispersion utilised is(4.37)In the case of uniform sediment = 0, and co should equal 1. For Oak Creek co = Co0. Thesimplest hypothesis assumes a linear variation in co which is consistent with uniform sedimentand Oak Creek:(4.38)where is the value of o obtained numerically from the Oak Creek data Fig 4.10.4.5.1 Calculation of Armour Grainsize DistributionThe usual application of the Parker surface-based relation is to input the values ofD3g, D1, aridand then together with TbedtO calculate the bedload transport rate per unit channel bed widthwhich is denoted g, when expressed in dry mass units, and qb in volumetric units. However ifDsg and o, together with the values of are unknown initially, it will be necessary to calculatethe value of these variables for imposed values of Gb. It is important to note here that while thetotal load Gb is an imposed value, the transport rates gb and q, are in effect dependent variables88as their values will be different for different channel geometries. By definition gb = Gb / Pbed, andsince the value of Pbed is considered in this thesis to be a dependent variable, the value of gb isalso a dependent variable. Furthermore changes in the channel geometry will also affect thevalues of md which in turn will affect the values ofgb and qb.By equating Eqns (4.28) and (4.33) and rearranging the following relation for is obtained:(s—i) gqb1,. = 15 (4.39)0.00218 G[sgo g0(s1)a] (Tb, I p)The fraction volume content of the transported bedload sediment in the jth range is denotedwhich is defined as:qb• (4.40)Combining Eqns (4.40) and (4.39) and rearranging yields:F—_____________(s—1)gq441— G[sgo go(s1)] 0.00218(rb / p)’5Since q andf are known for an imposed Gb and a known value ofPbed, the values ofF; and thecomposition of the equilibrium armour layer can be calculated for a value of r,ed from Eqn(4.41).However the optimization model to be developed in Chapter 7 evaluates a range of trial channelgeometries for which the value of q, is not known in advance. The calculation of an equilibriumarmour layer for this more general case will now be considered.89For a specific discharge the second term in Eqn (4.41) is constant, although the value of qb maynot be known initially. Eqn (4.41) can therefore be simplified to:(4.42)G[CI)so g0(31)a]Now by definition the following must hold:F;.=L0 (4.43)Therefore Eqn (4.42) can be presented as an equality by normalising over the summation of allfj/G[(I)sgo g0(o) a)](4.44)fj/G[sgo g0(o) wJEqn (4.44) permits the calculation of the values of F. without knowing the value of q,,beforehand.Since the values ofF, and therefore D and are unknown initially, an iterative procedure isrequired in order to solve Eqns (4.31), (4.36), (4.37), (4.38), and (4.44) in order to determinethe composition of the equilibrium armour layer, as well as q,. This will be discussed in Section7.2.5.By inspection between Eqns (4.44) and (4.41) it is evident that the value ofq is given by:90O.00218(r ip)15qb = (4.45)(s—i) gfj/G[go g0(o,) a]4.5.2 Modification of the Parker Surface-Based Relation for Natural Rivers withVariable DischargeThe discussion and theory presented in this section is an extension of the work by Parker(1990).The “equal mobility” hypothesis was derived originally from flume experiments which wereperformed using single, steady discharges (Parker et a!., 1982b). In its original form it impliesthat equal mobility occurs for all discharges above the critical for armour mobility, and thereforeassumes that the grain size distribution of the transported bedload is equal to the subarmour forall flows in excess of critical. This is known to be incorrect as is indicated by data from OakCreek (Milhous, 1972) which is shown in Fig 4.11. Note that below the critical Q of about 1m3/sec the value of d50 is approximately 2 mm and represents sandy throughput load movingabove an immobile armour. For flows in excess of 1 m3/s the value of d50 increases withincreasing Q.In its original form the equal mobility hypothesis supposes that d50 of the transported loadwould be a constant for discharges in excess of the critical, and for Oak Creek this constantvalue would be equal to 20 mm. As Fig 4.11 shows this is clearly incorrect. This wasacknowledged by Parker et aL (1982a) who state that equal mobility is at best a crude, butuseful, first-order approximation.However this hypothesis can be correctly applied to natural channels with varying discharge byrecognising that for a channel in equilibrium, the constraint of equal mobility must be fulfilledover a significant duration. Parker (1990) acknowledges that the size distribution of the annual91yield of gravel is similar to the subarmour. Similarly Church et al. (1991), although referring tofine sediment, state that “equal mobility is at best a statistical phenomenon which holds over a(significant duration)”.The experiments of Parker et a!. (1982b) and others have demonstrated that the development ofan equilibrium armour layer even under controlled conditions requires many hours. The flows inmany natural rivers, in particular ‘flashy’ gravel-bed rivers, are probably too variable for anarmour layer to develop an equilibrium state with a single discharge. However over a significantduration, which is usually at least one complete water year, a channel which is in equilibriummust have transported the bedload imposed on it without appreciable net deposition or scour,and therefore must have transported all size fractions at just the rate at which they were suppliedto the channel. That is all sediment sizes are rendered equally mobile. Therefore the armourlayer which develops in an equilibrium channel represents the grain size distribution necessary todevelop equal mobility for the complete range of flows experienced by the channel over asignificant duration. However any single discharge is likely to be in disequilibrium with thearmour because the sediment size distribution of the sediment being transported by a given flowat any given time is probably different than the subarmour.For each value Q, a value of 1•bed• is calculated from the equations presented in Section 4.3.2,from which can be determined. The value is the volumetric sediment transport rate perunit channel width for flow Q,. Eqns (4.28) and (4.33) must be modified as follows:(s—1)gq,.= (4.46)‘ (r,/p) ‘J =O.OO218G[cF0g0(81)a] (4.47)92where qb13 is the unit volumetric bedload transport rate for size fractionj for flow Q,, and thevalues Vbedj, Io1, and Co., and G are not constant, but now vary with Q.By equating Eqns (4.46) and (4.47) and rearranging the following relation for qb11 is obtained:O.00218F.qb= (s—l)gp’5 G[3goj g0(s) a]i, (4.48)By definition:= Pi (4.49)where p is the probability of Q,.Therefore from Eqn (4.48) and (4.49):O.00218F.(s— 1)g p15 1v G[goj g0(s) a] z, (4.50)Rearranging (4.50) together with (4.40) gives the following relationship for 13 which isanalogous to (4.41):F — fJ(s—1)gp5q (451)— p1G[0. g(o) ,] 0.002 18Since the second term of Eqn (4.51) is constant (although the value is not necessarily known inadvance), with Eqn (4.43) the following equation analogous to (4.44) is obtained:93f/P G[sgoj g0(o) c]r.F. = ‘‘ (4.52)G[so. g0(o) w}rBy inspection Eqns (4.51) and (4.52), it is evident that the total volumetric unit bedloadtransport rate per unit bed width q, is given by:0.00218qb = n m(4.53)(s- 1)g p’5 G[sgo. g0(o)€]If the SI system of units is used, qb will have units ofm2/sec.4.6 TOTAL BEDLOAD CONSTRAINTThe total bedload constraint for the channel subject to the complete range of flows defined inthe flow duration curve (Fig 4.7) is represented by:(4.54)where g,, is the total bedload transport rate per unit channel width in mass units, and Gb is anindependent variable and represents the total load imposed onto the channel from upstream. Thevalue ofgb is given by:g=Ispq (4.55)where r is a constant which is related to the time base of Gb. The value of qb in the totalbedload constraint is calculated using Eqn (4.53). The units ofg and Gb must reflect the loadimposed over a significant duration. For instance under steady flow conditions such as in alaboratory experiment Gb can have units kg / s as the rate is constant and will be the same for alltime scales. In natural rivers with highly variable and seasonal flows the logical time base would94be one water year, and therefore Gb would have units of kg / y. In general the modeling to beundertaken in Chapter 7 will assume a time base or significant duration equal to 1 year.Therefore, the value Gb represents the mean annual bedload supplied to the channel, and F takesa value equal to 365.25 24 3600 which is the number of seconds in one year.4.7 SUMMARYThe bedload constraint has been developed and consists of two components: boundary shearstress and sediment transport. In this thesis, the mean bed and bank shear stress (bank shearstress is necessary to assess the stability of the banks, see Chapter 5) are estimated using theempirical relations of Flintham and Caning (1988). These formulae were derived fromexperimental data in straight rectangular and trapezoidal flumes, and they permit the estimationof the mean values of rb€d and bank for trapezoidal channels with different bed and bankroughness. Flintham and Carling caution that these results may not be confidently applieddirectly to large scale natural channels or for bank angles less than 45°. This warning isacknowledged, however in the absence of anything better, the relations of Flintham and Caningwill be assumed to be valid for large natural channels with low angle banks.The effects on the shear stress distribution of weak, straight-channel secondary currents areimplicitly accounted for in the empirical shear stress relations. The strong secondary currentswhich are associated with curved channels are not accounted for in this thesis. The secondarycurrents are considered to exert only a second-order influence on the channel geometry and willresult principally in variations of the local geometry from the reach-averaged values. There issome support for this assumption in the regime equations ofHey and Thorne who found that thechannel width at pool and riffles sequences varied by less than 5% from the reach-averagedvalue.The Parker (1990) surface-based bedload transport relation will be used in this thesis to modelthe sediment transporting capacity of the channel in Chapter 7. This bedload transport relation95was selected because it explicitly accounts for the influence of the armour layer on the sedimenttransport, and the relation can be inverted to calculate the equilibrium armour layer grain sizedistribution that is required to transport an imposed bedload for given flow conditions. Thisallows the grain size distribution of the armour layer, together with the value of D50 to betreated as dependent variables in the optimization model. Therefore 135 can adjust together withW, Y, and S to define an optimum channel geometry.The bedload transporting capacity of the modelled channel is calculated using the completerange of flows represented by the flow-duration curve, rather than at a single discharge such asQbf The concept of “equal mobility” is restated whereby the armour layer adjusts such that allsize fractions are equally mobile over a given duration such as one year. The Parker (1990)relation is modified to apply to the full range of flows.96(a)(b) V/\\\‘s:==E— — —— Velocity ContoursWeak Secondary CurrentsFigure 4.1. Simplified velocity distribution in a straight open channel. (a) Momentum diffusiononly. (b) Momentum diffusion and secondary currents.I I I II I II II — I gf I g‘ /‘ — /I———“-‘97(a)(b)VdPFigure 4.2. Area methods. (a) Normals to isovels (Leighly, 1932). (b) Normals to bed(Lundgren and Jonsson, 1964).98OwyS’0.75OwySQ.97OwyS1.0 —0.90,V‘.‘ 0.78C0•00. 02__________DFigure 4.3. Boundaiy shear stress distribution from Lane (1955b). For Fig 4.3 only, b isequivalent to Pbed, y is equivalent to Y, and z is a measure of the bank angle which is equal tocotTropezois, z 2z I 50,0.8Q.0.7CCa.;V0.”‘iffl‘Troezods,z2ondI.5—/ecton9esI1 —“‘Tropezoids Z—,J.___‘——...l4— [IvLAH::iz:- - TflTTII —2 3 4 67b/y8 9 0 _0 2 .5b/y7 8 9 1099(a) 10.6a)C.0 0.4U,Ca)Ec 0.2C0z0(b) 10.80.6Cucr.iCC(0)C0)E0.2C0z0Figure 4.4. Boundary shear stress variation from Eqns (4.8) and (4.9) for variable bank angle6 and uniform boundary roughness in a straight trapezoidal channel. (a) Non-dimensional bedshear stress = Vbed / yJ’S. (b) Non-dimensional bank shear stress = Tbank / rYS.wIY100(a)10.80.6a)a)o 0.4Cl)Ca)E130.2e0z0(b) 1h 0.8Cl)Ca)a)C0Cl)Ca)UC0z0Figure 4.5. Boundary shear stress variation from Eqns (4.8) and (4.9) for constant 8(8 = 90°)and non-uniform boundary roughness. (a) Non-dimensional bed shear stress = id / yYS. (b)Non-dimensional batik shear stress = z / yYS.5 10 15 20 25w/Y101140zCl)C,)a)120100806040200Total Shear Stress, r (N/rn2)Figure 4.6. Grain shear stress versus total shear stress. The value of r is computed from theobserved channel geometry, and r’ is calculated from Eqn (4.15).0 20 40 60 80 100 120 140102C)a)EFigure 4.7. Discretized Flow-Duration Curve. The probability of exceedence refers to aspecified time base and period of record. The value Q is the characteristic discharge valueassigned to the interval, and Pi is the probability of a discharge being within a particularinterval.0Probability of Exceedence103(a) 0.120.1 -D‘ 0.08-0.06 -00.04- 0DDO0.02- DCB0 I0 1 2 3 4Q (me/s)(b) -0.01 -0.001 -0.0001 -1E-05-1E-06-f1E-07 -°1E-080 2 3 4Q (m3/s)Figure 4.8. Bedload transport rates versus Q for Oak Creek, Oregon, USA. Data from Milhous(1972). (a) Linear Y-axis. (b) Logarithmic Y-axis.104(a) 500400300G20010000.5(b)1001010.10.010.0010.00011E-050.5Figure 4.9. Variation of the function G from Eqn (4.32). (a) Linear Y-axis. (b) Logarithmic Yaxis.1 1.5 2•sgO1 1.5 2øsgO1052COo1.5—1 2 3 41’sgOFigure 4.10. Variation of the functions a and with q5 (from Parker, 1990). See text fordiscussion.10630002500d subarmour2000E 0000005 000E0,0:00023 4Q (m3/s)Figure 4.11. Variation of d501 d with Q for Oak Creek, Oregon, USA. Data from Milhous(1972). The value ofd50 from the subarmour sediment is approximately 20 mm.107CHAPTER 5BANK STABILITY CONSTRAINT5.1 INTRODUCTIONIn this chapter the bank stability constraints which apply to gravel-bed rivers will be developed.It is an obvious requirement that a stable channel must have stable banks. The bank stabilityconstraint will be formulated for both cohesive and noncohesive bank sediments.There are two fundamental processes of bank erosion which must be considered: mass failureand the fiuvial erosion of discrete grains or aggregates (Wolman, 1959; Thorne, 1982;Gressinger, 1982).5.2 COHESIVE BANK SEDIMENTCohesive channel banks are those composed of clay, silt and fine sand and have typicallyresulted from over-bank deposition of fine sediment from suspension. Cohesive sedimentscontain significant amounts of clay minerals which strongly influence the physical properties ofthe soil. According to Raudkivi (1990, p. 300) when the clay content of a soil is 10% orgreater, the physical properties of the soil are dominated by the clay fraction.5.2.1 Mass FailureMass failure of a river bank occurs when the driving force FD, which is due to the weight ofthe failed soil mass w; exceeds the resisting force FR, which is a result of the cohesion c andinternal friction angle çS of the soil acting along the failure surface. The unit shear strength of asoil is usually given by the Mohr-Coulomb failure criterion:108SR=c+NRtanØ (5.1)where SR is the shearing resistance for a unit area of soil, and is the normal upward forceacting on the base of the soil unit.Potential failure will occur along the surface where the factor of safety, F,, which is defined bythe ratio FR /FD, is a minimum. Failure will occur when F, is less than one. The critical heightof a bank corresponds to F, = 1, and can be expressed by the following equation fromTerzaghi and Peck (1948):Hmt=N (5.2)where y is the bulk unit weight of the soil and includes the weight of the pore fluid, and N isthe stability factor. Note that N as presented in Terzaghi and Peck (1948) and in Eqn (5.2) isthe inverse of Taylor’s (1948) stability number. The value of N, is a fi.mction of and thebank angle 6.The simplest model of slope failure is the Culmann analysis (Taylor, 1948; Spangler andHandy, 1982) which assumes that the failure surface is planar and passes through the toe of thebank (Fig 5.1). This method yields the following analytical solution for N:4 sinO cosçSN= 1—cos(6—.) (5.3)The slope of the failed surface, a, is obtained by differentiation:6+ça=2 (5.4)(See Fig 5.1).109Application of the Culmann solution is limited. It can only be applied with reasonableconfidence to completely drained slopes close to vertical. For lower angle slopes the failuresurfaces are strongly curved and the Culmann analysis under predicts N5, and thereforeThis will be demonstrated subsequently.No simple analytical solution exists for non-planar failure surfaces. The friction circle method(Taylor 1948), or the various methods of slices (Fellenius, 1936; Bishop, 1955; Morgensternand Price, 1965; Spencer, 1967) all require trial failure surfaces to be assessed. The frictioncircle method has in general been superseded by the method of slices and is seldom used today.For the method of slices the soil mass above each trial failure surface is divided into a numberof vertical slices of convenient width (Fig 5.2(a)), and the driving and resisting forces areresolved for each soil slice (Fig 5.2(b)). The factor of safety is established for that particularsurface by summation of the forces acting on all of the slices. A number of trial failure surfacesare examined until the surface with the lowest value of F5 is determined. If the value ofF3 isless than 1 the slope will fail.The two most commonly used methods of slices are the ordinary or Swedish method(Fellenius, 1936) and the Bishop simplified method (Bishop, 1955). For these and othermethods of slices, the number of unknowns exceeds the number of equations available forsolution and therefore some simplifying assumptions are necessary. In the ordinary method ofslices the vertical shearing forces between the slices are neglected, and the resultant force on aslice is assumed to act parallel to the base of the slice. This can lead to an under estimation ofthe factor of safety by as much as 60% (see Nash, 1987). In Bishop’s approach the verticalshearing forces between adjacent slices are considered in the analysis and are assumed to beequal and opposite. A key feature of the Bishop simplified method is an increase in the normalforce at the base of each slice which tends to compensate for the side forces on the slices.110The solution for the factor of safety using the ordinary method of slices and assuming acircular failure surface in terms of effective stresses is given by:FR (c’l+(wcosa—crl)tanç5’)SF wsinawhere o- is the mean pore pressure for the slice, and 1 is the basal length of each slice. The soilstrength parameters c’ and çS’ are given in terms of the effective stresses. The solution for theBishop simplified method is more complex and requires an iterative solution for F3.The ordinary and Bishop methods can also be applied to non-circular failure surfaces. Furtherdiscussion of slope stability methods can be found in soil mechanics texts such as Chowdhury(1978), Spangler and Handy (1982), or in a review by Nash (1987).The various methods of slices, as well as other more recent limit analysis methods whichrequire numerical techniques such as finite element analysis for solution (eg. Chen 1975; Chenand Lui, 1990), can be used to construct stability curves which show the variation ofN5 with &Examples of stability curves are shown in Fig 5.3(a). Those solutions indicated by the solidlines were derived from the graphs and tabulated data published in Taylor (1948, p.457) whichwere obtained using the friction circle method. The dashed lines are values of N from theCulmann solution, Eqn (5.3). For slopes close to vertical the Culmann solution agrees quiteclosely with the results of the friction circle method, however the two solutions diverge forlower values of 6The friction circle method agrees closely with results obtained from the ordinary method ofslices for simple failure geometries with homogeneous soils and no seepage (Taylor, 1948).However the friction circle method cannot be applied readily to more complex geometries,heterogeneous soils, or where seepage is significant, and therefore the method is not widely111used today. However the stability curves developed by Taylor (1948) can be applied to somespecial cases and these will be discussed in a following section.With the stability curves the value of can be calculated for any value of 9 and particularsoil values y c’ and ‘. Alternatively stability curves of versus 9 can be constructed forknown values of y,, and c’ as in Fig 5.3(b).To develop stability curves for an applied river bank analysis, field investigations are necessaryto determine the bank stratigraphy, the properties of the bank sediments, and the likely shapeof the failure surface. Also other variables such as the presence of surface tension cracks andthe influence of bank vegetation on the stability of the bank can be considered. Groundwatermonitoring might be necessary to determine the effect of pore pressures and seepage forces.Construction of the stability curves may require considerable effort, however once available thecurves can be approximated by empirical equations or piece-wise linearized segments, and thenthis information can be accessed quickly by the optimization model. This is necessary as theoptimization model must assess the stability of hundreds or thousands of H-9 combinationsbefore arriving at the optimum solution.5.2.1.1 Bank Stability and SubmergenceThe effect of submergence and emergence of the channel banks due to flow variability plays akey role in the stability of the channel banks. Bank submergence results in increased soil weightdue to saturation, and a decrease in the friction resistance due to increased pore waterpressures. These two effects may combine to reduce the stability of the banks. This decrease inbank stability may however be countered by increased hydrostatic loading along the banksurface.Three limiting cases will now be considered that can be analysed using the stability curvesdeveloped by Taylor (1948). These cases are shown in Fig 5.4 and will be used to demonstrate112the effect of submergence, saturation, and emergence on the stability of a river bank. Thevalues of c’ and ‘ to be used in the following discussion are assumed to equal the developedvalues, that is the factor of safety applied to the values c’ and qS’ equals 1, and the values c’ andare assumed to be known exactly.Case I corresponds to a condition of zero or low flow in the channel, together with drainedbanks. The value of N3 can be determined directly from the stability curves for a particularvalue of 0, and Hcrjt calculated for the known values of y and c’.Case II corresponds to that of a fhlly submerged bank. The soil mass is fully saturated and isbeing supported by buoyancy forces about its perimeter as indicated in Fig 5.3. For this casewhen the soil is fully saturated, and there is no seepage, the buoyancy forces are hydrostatic.The value ofN3 can be determined directly from the stability curves, and then calculatedusing the buoyant unit weight, Yb’ which is equal to saturated Yt minus y.Case III corresponds to a condition following complete and instantaneous drawdown of theflow in the channel. In this case the soil in the bank remains saturated, the pore pressuresunchanged from the fully submerged case. However the hydrostatic forces which werepreviously acting to stabilise the bank are now absent. Taylor (1948) determined that the valueofN for case III could be approximated by modifying the value of ‘to account for the porepressures acting along the failure surface which reduce the frictional resistance. It wasdetermined by Taylor (1948, p. 467) that the modified friction angle 9S,,, can be approximatedby:(5.6)Where y1 is the saturated value. This relation, which applies only to complete and instantaneousdrawdown, has been confirmed by Morgenstern (1963) who developed stability curves for113complete and partial instantaneous drawdown using the Bishop simplified method of slices. ForCase III the value of N., is determined from the stability curve corresponding to qSm, andcalculated using y1.The values ofN., and 11cr 1 will now be determined using reasonable soil values to demonstratethe relative stability of the three cases. A value of q’ = 25° together with the drained Yt = 20.80kN/m3, saturated y = 22.76 kNIm3,and Yb = 12.95 kN/m3 for 0= 700. From Eqn (5.6) thesevalues give a value of Øm = 14.2°. The values for ]V are read from the stability curve fromTaylor (1948) in Fig 5.3a. The results are presented below in Table 5.1.Table 5.1. Parameters for 3 cases of bank stability.Parameters Case I Case II Case IIIYr OT Yb (kN/m3) 20.80 12.95 22.76qY (°) 25.0 25.0 14.2c’ (kN/m2) 10.0 10.0 10.0M, 9.8 9.8 7.2(m) 4.71 7.56 3.16The values of in Table 5.1 indicate the role of submergence and emergence on the stabilityof the channel banks. Note from Case II that the effect of increasing the flow depth is tostabilise the channel banks and increase Therefore the limiting condition for bank stabilitymust be satisfied at low or zero flows. Cases I and ifi defined upper and lower limits for thebank stability. The critical period for bank stability is during flow recession. Case I correspondsto a river with very slow runoff response and/or very permeable bank sediment. As the flowrecedes the banks must be dewatering so that the bank above the water surface in the channelis fhlly drained, and the bank below the water surface fi.illy saturated. This continues until zero114flow condition and the banks are totally drained (Case I). The other limiting conditioncorresponds to a river with flashy runoff and/or highly impermeable channel banks. The flowrecedes at a rate at which no dewatering of the bank sediment occurs and the banks remainfl.illy saturated at the zero flow condition (Case III). The stability of natural rivers would liebetween the two limiting cases (I and III).5.2.1 2 Bank Height ConstraintThe bank height constraint requires that the bank height H, be less than or equal to the criticalbank heightHH (5.7)where H is equal to the bankfull flow depth 1’, the value for is calculated by Eqn (5.2) withN obtained from the stability curve. The value of should be evaluated for both Case I andCase III conditions. Examples of stability curves derived from field studies can be found inMarch et a!. (1993) and are shown in Fig 6.9. The stability curves and analysis of March et a!.are discussed in Section 6.4.5.5.2.2 Fluvial ErosionFluvial erosion is the process whereby individual grains or aggregates on the surface of thebank are removed and entrained by the flow. The driving force is the shearing action of thefluid on the bank sediment. The mechanics of fluvial erosion of cohesive sediment remainpoorly understood despite considerable research. Reviews of the research conducted on theerodibility of cohesive sediment are found in ASCE (1967), Graf (1971), Gressinger (1982),and Raudkivi (1990). The resultant interparticle force for cohesive sediment is the net result ofseveral forces of attraction and repulsion. These forces result from very complex electrochemical processes which include the clay mineralogy and content, and the temperature andchemistry of the pore and eroding fluids (Arulanandan et a!., 1980; Gressinger, 1982;115Raudkivi, 1990). Unlike noncohesive sediment, the weight component of the cohesive grains isquite insignificant when compared to the electro-chemical forces. The electro-chemical forcesmay be orders of magnitude larger than the weight forces of the individual grains (Raudkivi,1990, P. 310).Numerous experiments to determine the erodibiity of cohesive sediment have been performedusing a variety of techniques such as straight and circular flumes, rotating cylinders, submergedjets, impellers, and even portable flumes that can be used to test undisturbed soil in situ.Attempts to quantify the fluvial erosion of cohesive sediment often use the concept of criticalshear stress, critical velocity, or critical stream power which must be exceeded before erosioncommences.Attempts have been made to use the results from these erosion studies to develop correlationsbetween the erodibility of the cohesive sediment and other primary physical properties of thesoil which are more readily obtained. These physical properties include clay content, claychemistry such as the Ca/Na ratio and pH, grainsize, vane shear strength, plasticity index,organic matter content, as well as the temperature and chemistry of the eroding fluid. Ingeneral this approach has not been overly successful and has not yielded results which could bereadily used for predictive purposes. However the approach of Arulanandan et al. (1980) hasproduced charts which relate soil erodibility to the electrical conductivity and the magnitude ofthe dielectric dispersion of the soil, properties which reflect the type and amount of clay. Itaddition it was shown that the chemistry of the eroding fluid, expressed as the sodiumadsorption ratio, also significantly influences the soil erodibility. According to Arulanandan etal. (1980), while this study “did not yield a quantitative method to predict critical shear stress”,the charts “should give a reasonable estimate of r for a natural undisturbed soil”. Notethat the same soil can have different values of r for different eroding fluids, and therefore thevalue of corresponds to a particular soil-water system.116The critical shear stress concept will be adopted in this thesis. The critical shear stress of aparticular soil ;,,, can be evaluated directly using one of the erosion devices mentionedpreviously, preferably on undisturbed in situ samples, or it can be estimated indirectly bymethods such as Arulanandan eta!. (1980).5.2.2.1 Bank Shear ConstraintThe bank shear constraint requires that the bank shear stress be less than the critical valuerequired for fiuvial erosion of the bank sediment:;aflk (5.8)where the quantity Vbank is calculated from Eqn (4.8). This constraint assumes that the gravitycomponent of the eroded grains or aggregates is negligible compared to the cohesive forces,and therefore the value of z,rjt is independent of 6 and can be considered to be an independentproperty of the soil-water system.The assumption that a soil water system has a single definable value of belies thecomplexity of the erosion process. For example freeze-thaw action, ice wedging, anddesiccation and tension cracks can all result in aggregates of cohesive soil which may have asignificant gravity component and a much lower value of r than a smooth homogeneous soilmass. Nonetheless for modeling expedience it will be assumed that a soil-water system can becharacterised by a single value of z which is treated as an independent variable.5.3 NONCOHESIVE BANK SEDIMENTAs with cohesive sediments both mass failure and fiuvial erosion will be considered fornoncohesive sediment.1175.3.1 Mass FailureFor a bank composed of noncohesive sediment the resisting force that balances the drivingforce is due only to the friction term in Eqn (5.1). It is an elementary exercise to show that thefactor of safety for an infinite, drained slope of noncohesive sediment is given by:F = tanq5 (5.9)tan8from which the limiting stability occurs when the bank angle 8 is equal to the angle of internalfriction.Therefore for a bank to be stable the following constraint must be satisfied:8qS (5.10)If 8 exceeds 0, bank failure will result.When a bank is fully submerged as in Case II from the previous section, the stability isunchanged from the fully drained case as the buoyancy forces reduce the driving force and theresisting force by the same amount. For a condition of rapid drawdown seepage forces candevelop which reduces the stability of a noncohesive bank. Since only frictional resisting forcesare involved, the decrease in bank stability due to seepage forces can be simulated numericallyby modifying (ie. reducing) the value of ç. For the case where seepage is parallel to the slope,the modified value is Om is given by: (eg. see Spangler and Handy, 1982; p. 492):tan0m=tan0 (5.11)‘FtIn this case the value of is analogous that obtained by Taylor (1948) for cohesive sediment(Eqn 5.6).118However in many banks composed of noncohesive sediment, particularly those formed ofcoarse gravel, significant seepage forces will not develop due to their high permeability. Theriver with noncohesive banks will tend to behave as in Case I for cohesive banks, and Eqn(5 10) is sufficient to assess the bank stability with respect to mass failure.5.3.2 Fluvial ErosionIn flowing water the fluid forces exerted on the grains as well as the down slope gravitycomponent of the grains must be resisted by the frictional forces for a grain to remain stable.The well-known bank stability analysis described below was originally developed by the UnitedStates Bureau of Reclamation (USBR) to determine the stability of noncohesive banksediment. Its development is summarised in Lane (1955b) and it is presented in the followingsimplified form from Chow (1959, p.171) and Henderson (1966, p.419) for the limiting bankstability:Tbaflkc=—2 6 (5.12)Tbedc J sin2Øwhere rbaflkc is the critical shear stress for a grain located on a sloping bank, and Vbed is thecritical shear stress for the same grain located on the sub-horizontal channel bed.The original equation can be rewritten in a dimensionless form as follows:Vbank I S1fl6( _in 2 (5.13)r\s sinwhere s is the specific gravity of the sediment, D5Obaflk is the median bank grain diameter whichis assumed to be representative of the bank sediment, bed is the critical dimensionless(Shields) shear stress for a grain equivalent to D50bank on the channel bed.119Unlike the mobile bed sediment, the gravel comprising the channel banks may be consolidated,cemented by fine silt and clay, and stabilised by roots which have penetrated the gravel banks.Eqn (5.13) was developed for unconsolidated noncohesive sediment and must be modified toaccount for the properties of the bank sediment.USBR data (in Lane, 1955b) indicate that qS for coarse gravel approaches a maximum ofapproximately 400. However consolidation and imbrication of the bank sediments, togetherwith the effect of cohesive silt and clay cementing the grains and other stabilising influencescan increase the in situ friction angle above 0. Therefore çS in Eqn (5.13) can be replaced by thein situ friction angle, Ø, which can be allowed to take a maximum value up to 900.The value of t*bedc must also be adjusted because the critical shear stress for the stabilisedbank material would be higher. Henderson (1966, p. 414) has shown from the work of White(1940) that the dimensionless shear stress is related to 0:r=ktan (5.14)where k is an empirical constant. The value of k will be determined from field data in Chapter6.To include the influence of bank vegetation, consolidation, and cementation of the banksediments, 0 is replaced with Substituting Eqn (5.14) into (5.13), and using , theresulting constraint for noncohesive banks is:bank ktanqS1ll— Sifl 6 (5.15)sin120Note that Eqn (5.15) is presented as an inequality constraint with the mean bank shear stress rbank’ rather than ;ank in the numerator on the left hand side. The value 6,, can take a maximumvalue up to but not including 900, as tan 90° is infinity.By analogy Eqn (5.10) should be written in terms of the in situ angle of repose:(5.16)Eqn (5.15) automatically satisfies Eqn (5.16) for all values of çb< 90° because the square-rootterm becomes undefined for values of 9>. Eqn (5.15) therefore simultaneously constrainsthe bank stability with respect to both mass failure and fluvial erosion, unlike the noncohesivebanks where two separate constraints are required.The influence of the bank stability on the channel geometry will be demonstrated in Chapter 6.5.4 SUMMARYThe bank stability constraints were developed for both cohesive and noncohesive banksediment. For both sediment types there are two mechanisms of bank erosion that must beconsidered, namely mass failure and fluvial erosion of individual particles.For cohesive sediment the stability of the bank with respect to mass failure can be assessedusing slope stability techniques. To incorporate the slope stability data into the optimizationmodel stability curves are constructed. These curves can be accessed rapidly by theoptimization model. The banks are most susceptible to failure at low or zero flows when thehydrostatic supporting force from the flow is absent. Two limiting cases are recognised, fullydrained, and fully saturated bank conditions. For bank stability the cohesive bank sedimentmust also be stable with respect to fluvial erosion. The concept of critical shear stress forcohesive sediment is used. The bank shear constraint requires that the value of Tbank be less121than or equal to the critical shear stress for the bank sediment. The value ;ank represents themean bank shear stress. In practice it is the maximum bank shear stress that would likely occurnear the toe of the bank, rather than ank that would determine the bank stability with respectto fluvial erosion. Nonetheless ;ank will be used.To assess the bank stability of noncohesive bank sediment the USBR bank stability presentedby Lane (1955b) will be used in a modified form. The original equation is modified to reflectthe increase in the stability of the bank sediment due to packing, consolidation, imbrication,cementing by fines, and binding by root masses. This increase in stability is reflected in the insitu friction angle,that can take a value up to 900. The stability of noncohesive banksediment with respect to both mass failure and fluvial erosion is assessed by the single USBRrelation except in the case of large seepage forces following rapid drawdown. The influence oflarge seepage pressures in noncohesive sediment will not be considered as the high hydraulicconductivity of these sediments, particularly for coarse gravel sediment, would facilitate rapiddraining of the bank.Gravel-bed rivers commonly have channel banks that are composite in nature and typicallyhave a lower noncohesive unit which is overlain by an upper cohesive unit. In the case ofcomposite banks it must be determined which of the layers is controlling the bank stability.Often the lower noncohesive unit is eroded by fiuvial activity, and the overlying cohesive layerthen undergoes mass failure as it becomes undermined (Thorne, 1982). In this example itappears that the stability of lower noncohesive unit determines the stability of the bank.The influence of bank stability on the channel geometry will be examined in Chapter 6.122Figure 5.1. Bank stability analysis for planar failure surface.HFD = wsina= cL + NRtan123(a)(b)Figure 5.2. Bank stability analysis for the method of slices. (a) Definition sketch. (b) Forcesresolved on a single slice. The lateral forces not indicated.HL124(a) 252Oa,.0Ez 15.01o5(b)4Figure 5.3. Stability curves used to determine Hrit. (a) Stability number as a function of (b)Critical height as a function of Ofor prescribed values of c’ and rt.8(i)1086220 40 60 80e ()125Case IPore Water PressuresHydrostaticPressuresFigure 5.4. Three special cases for bank stability analysis. (a) Case I: Drained banks. (b) CaseII: Fully submerged and saturated. (c) Case ifi: Instantaneous and complete drawdown, bankfi.illy saturated. The hydrostatic and pore water pressures are indicated.Case IIACase III126CHAPTER 6EFFECT OF BANK STABILITYON CHANNEL GEOMETRY6.1 INTRODUCTIONIn Chapter 6 simplified versions of the optimization model will be presented whereby thesediment transporting capacity of the channel will be calculated at the bankfhll discharge only,and rather than the modified Parker (1990) relation developed in Chapter 4, a simple bedloadtransport relation will be used. This simplified model is inadequate to assess the response of anatural channel to variations in the discharge or sediment supply, but is sufficient to investigatethe influence of the bank stability on the channel geometry, and to introduce the optimizationmethodology to be fully developed in Chapter 7.6.2 BANKFULL : fiXED-CIIANNEL-SLOPE OPTIMIZATION MODELIn the bankfiill model the natural variation in the flows is reduced to a single dominant orchannel-forming discharge which is assumed to be equal to the bankfull discharge, Qbf Thesediment transporting capacity of the channel is indexed by the sediment transport ratecorresponding to the bankfull discharge, Gbf This simplification has been used by Chang127p(1980), White et a!. (1982), Millar (1991), and Millar and Quick (1993 a, b) in theiroptimization analyses of gravel rivers, and Hey and Thorne (1986) in their empirical regimeanalysis. The limitations of this approach with respect to calculating the sediment transportingcapacity of the channel were discussed in Chapter 4. The bankfiull model is however adequateto investigate the effect of the bank stability on the channel geometry.The banlcfull model will be developed for both noncohesive and cohesive channel banks. In thefirst version a further simplification will be made in that the channel slope S will be treated as anindependent variable. This is done to assess the effect of the bank stability on the values of Wand Y without the interference of varying channel slope. In Section 6.5 the model will begeneralised to include variable S.The bankfull:fixed-channel-slope model is analogous to an experimental setup where the slopeis fixed, and the channel width, depth, and sediment transport rate adjust to the imposeddischarge (eg Wolman and Brush, 1961). The model is formulated below.6.2.1. Independent VariablesThe variables that are assumed to be independent with known values are Qb S, d50,D50, k3bj,the bank stability variables which are D5Oba and ç for noncohesive banks, and c’, ‘, y,and for cohesive banks, and the following which are generally constant: g, p, s, v(kinematic viscosity), and y Since the value of S is specified, the value of Gb.,- is not required inthis formulation.1286.2.2. Dependent VariablesThe primary dependent variables to solve for are Pbed, Pbank, and 6 From these primaryvariables the secondary dependent variables such as j U Rh, W, Y, gb, and others can bedetermined.6.2.3. Objective functionThe objective function in this formulation is given by:maxf(,1,6)=77 (6.la)where:(61b)pQbfSAs the values of p, Qb, and S are prescribed, maximization of i is equivalent to themaximization of the product Pbed times gb. The value ofgb will be calculated using the Einstein-Brown sediment transport relation as presented in Vanoni (1975, p. 170). This relation requiresonly the value ofd50 to calculate gb:* 12.15 exp(_0.391 / <0.093gb (6.2a)L4o r 0.09350 50where gb*= dimensionless bedload transport rate per unit width given by:129* gbgb=1 ps,.J(s—1)gd50The dimensionless bed shear stress for the median bedload-grain diameter is given by:*— Tbed 62‘Cdr(sl)dCThe value of ‘Cbed is calculated from Eqn (4.9).The parameter F1 is related to the fall velocity of the sediment d,0 and is given by:12 36v2 I 36v21j3gd503(s—1) 1Jgd503(s—1) (6.2d)where v = kinematic viscosity of water. For gravel-sized sediment, F1 is approximately equal to0.82.6.2.4 ConstraintsThe only two constraints that are required in this formulation for uniform flow in a prismaticchannel are continuity and bank stability. As the value of S is prescribed the bedload constraintis not required.1306.2.4.1 ContinuityThe continuity constraint determines the dimensions of a channel which can just convey theimposed QbfThe continuity constraint for uniform floe in a prismatic channel is defined as:UA=Qbf (3.24)where U is the mean velocity, and A is the cross-sectional area of the flow. The value of U isobtained from the Darcy-Weisbach equation:I8gRSU=f (3.25)The value of the friction factorfis calculated using Eqn (3.19).6.2.4.2 Bank Stability ConstraintThe single bank stability constraint for noncohesive bank sediment is given by:Tbank I SiflOr(s—l)D501,ktanq1_Sfl2 (5.15)131where the value of Thank is given by Eqn (4.8). The value of the constant k is to be evaluatedfrom field data in Section 6.3.1.For cohesive bank sediments the bank stability constraints consists of the bank-height and thebank-shear components:HHmt (5.7)Tbank r, (5.8)Both constraints must be satisfied for cohesive banks to be stable.6.2.5. Optimization Scheme: Fixed-Channel-SlopeThe optimization flowchart is shown in Fig 6.1 and the source code for the computer programis given in Appendix A. The program is a step-wise iterative procedure that varies only onedependent variable at each stage. The program is initialised by selecting trial values of the threeprimary dependent variables, Pbed, Pbank, 0. The continuity constraint is satisfied by varying thevalue OfPbank for trial values OfPbed and When the continuity constraint has been satisfied forQbj the dimensions of the trial channel have been established.The bank stability is then assessed. The maximum value of 6 for which the banks are stable,6,,, is determined. Other values of 0 less than 0,,, may also satisf’ the bank stability132constraint, however these lower values do not represent the optimum bank angle. Fig 6.2shows the typical variation of Tj and r,ank for a range of values 6 for which the continuityconstraint has been satisfied. Note that the value of rd decreases monotonically withdecreasing 6 Therefore while values of 6 less than 9,, may satisfy the bank stability constraint,these values are associated with lower values of r, and is therefore the channel is less efficientwith respect to transporting sediment than a channel with a bank angle equal to O,,,. The finaloptimal solution must have a value of 8 that corresponds to a value of O,,,. Trial values of 8are assessed until 6m is determined. For each trial value of 6, a new value ofPbank must alsobe determined which satisfies the continuity constraint.Once the value of 8,,, has been determined, the sediment transport efficiency i is assessed. Thevariation of i over a range Of Pled values for a prescribed value of S is shown in Fig 6.3. Eachpoint on the solution curves in Fig 6.3 satisfies the continuity and bank stability constraints. Thevalue of i is reduced to zero either when Pbed equals zero, or when Pbed is very large and Ybecomes very small, and the value of r,,ed, and hence gb, both approach zero. The optimal valueof i is located between these two extremes, and is found by varying Pbed.The bisectrix method was found to be the most suitable for satisfying the constraints andlocating the optimum. This method requires that the initial upper and lower bounds of thesearch be initially specified for the variable in question. The midpoint between these twobounds is evaluated and, depending upon the outcome, becomes the upper or lower bound forthe next stage of the search, thus reducing the search area by half For example to satisfy the133continuity constraint upper and lower values of Pbank are established which are sure to containthe final value of Pbank. The midpoint between these two values is evaluated. If the calculateddischarge capacity of the channel at the midpoint is greater than Qj the midpoint Pbank value istoo large, and this value then becomes the upper bound for the next stage. Conversely if thedischarge capacity is less than Qb1 the midpoint value of Pbank then becomes the lower boundfor the next stage. Convergence is obtained when the calculated discharge capacity falls withina selected tolerance of Qb1 say ± 0.1%.A similar procedure is used to evaluate the optimal value ofPbed. The first derivative dr7 / dPbedis evaluated numerically by cental differencing at the midpoint between the upper and lowerbounds to Pbed. For values of di7 / dPd> 0 the midpoint becomes the lower bound, and forvalues ofdi7 / dPd <0 the midpoint becomes the upper bound for the next stage. Convergenceis attained when the separation between the upper and lower bounds of the search is reducedbelow a preset limit, say 0.01 m.The bisectrix method was found to be the most robust technique for obtaining convergence,although it is computationally intensive. Other convergence schemes such as the NewtonRapson method were generally more computationally efficient, however they were prone toinstability and occasionally caused the program to crash by returning negative values for thedependent variables.This iterative search scheme assumes that the functions are well behaved with no local optima.134This is a necessary requirement for this type of iterative search procedure, as well as othernonlinear optimization techniques such as the reduced gradient method, which can get caughtup in local optima. At a local optimum all of the convergence criteria can be satisfied and yet itnot possible to determine whether a local or global optimum has been found. Different optimalsolutions can be obtained from different starting points for the search procedure. The variationof the value of the objective function i with Pbed for the fixed slope bankfull model is shown inFig 6.3. The objective function is smooth with only a single global optimum. Evidence of localoptimum has not been observed in any of the modelling undertaken.In addition to , the variation of 6 with Pbed is also shown in Fig 6.3. The values of 6 are smallfor low values OfPbed because the high shear stresses acting on the banks. As Pbed increases, thevalue of Tbank decreases and banks with higher values of 6 are stable. In the example used in Fig6.3 the value q = 400 was used. As Pbed becomes very wide the value of 6 approaches .The optimization model will now be tested using observed gravel river data from channels withnoncohesive and cohesive bank sediment.6.3 NONCOHESWE BANK SEDIMENTThe analysis to be presented in this section is similar to that described in detail in Millar andQuick (1 993b), however there are significant differences which result in differences between thenumerical values obtained herein, and those from Millar and Quick. The general conclusions arehowever unchanged.135The model will be tested on the published gravel river data collected by Andrews (1984) andHey and Thorne (1986). The published data will be used as input to the model, and the outputgeometry will be compared to the observed geometry. The relevant data and modeffing resultsare tabulated in Appendix D.These rivers are described as stable single-thread channels with mobile beds that activelytransport sediment at higher stages. The banks are either comprised of noncohesive gravelsimilar to the material being transported, or have composite banks with a lower noncohesiveunit, overlain by a cohesive silty unit. For the composite banks it will be assumed in this analysisthat the stability of the banks is determined by the lower noncohesive unit. The banks arecharacterised in terms of their vegetation density. Andrews (1983) subdivided his data set intothose channels with thin and thick vegetation. Similarly, Hey and Thorne (1986) subdividedtheir data into four bank Vegetation Types ranging from Vegetation Type I, grass with no treesor bushes, to Vegetation Type IV, with> 50% trees and bushes.The empirical regime analyses of Andrews (1984) and Hey and Thome (1986) determined thatthe effect of the bank vegetation was to increase the stability of the channel banks whichresulted in narrower channels for the same value of Qbf In the first stage of the present analysisonly the channels which have low densities of bank vegetation will be assessed. It will beassumed that the bank stability of these channels is determined by the bank sediment propertiesalone, and is unaffected by the bank vegetation. This preliminary analysis will permit the136estimation of the bank stability parameters. In subsequent analyses (Section 6.3.2 and 6.3.3) thechannels with higher densities of bank vegetation will be examined and the effects of increasedbank stability on the channel geometry due to the bank vegetation will be demonstrated.The values of several key parameters were not explicitly given in the data sets. The value ofd50was not given by Hey and Thorne, and as an approximation D50 I 3 will be used. Thisapproximation will affect the absolute value of g calculated for the channel, but will notinfluence the location of the optimum significantly.Neither Andrews nor Hey and Thorne specified the values of D50k. In this analysis the valueOfD5Oba, will be assumed to equal D50, the median grain diameter of the armour layer sediment.Banks composed of sediment similar in size to the transported bedload sediment would beexpected to have a value of D5Oba,,k similar to d50, the median bedload or subarmour graindiameter. However previous work (Millar, 1991; Millar and Quick, 1993b) has shown that verypoor results were obtained using D5ob = d50, and much better results were obtained usingD5Qba = D50. This is consistent with the development of a coarse static armour layer on thechannel banks as the finer sediment is preferentially removed during bank erosion.Alternatively, the banks may be stabiised through the accumulation of coarse gravel at the toeof the bank.The assumption ofD5o0,g = d50 will be shown to be valid in Section 6.3.1.137The relative roughness values of the bed and banks are not known. Therefore since the bed andbanks are composed of similar sized gravel, the values of ksd and ks,flk will be assumed to bethe same, and both equal to the value of k which was calculated from the observed channelgeometry by inverting Eqn (3.2).6.3.1 Low Density Bank VegetationThe channels which will be analysed in this section are the channels described as having thinbank vegetation from Andrews (1984), and the Vegetation Type I from Hey and Thorne(1986). There are a total of 27 channels.The bank stability constraint, Eqn (5.15) contains two stability parameters: k which is related tothe critical dimensionless shear stress V*be4 by Eqn (5.14), and the insitu friction angle of thebank sediment.A value of çS, = 400 will be used as input which assumes that the banks arecomprised of loose, noncohesive gravel. A wide range of values for VCbe are cited in theliterature, with most faffing between 0.03 (Neil, 1968) and 0.06 (Egiazaroff 1965). For a valueof ç = 400 this translates into a range of values for k between 0.036 and 0.072. The value of kwithin this range that gives the best agreement between the observed and modelled channelwidths for the 27 channels will be determined, and then assumed to be a constant for allchannels.The model was run for a range of values of k to determine the best fit between the observedand modelled channel widths. The optimum value was determined to be k = 0.048 which is138equivalent to a value of Vbe4 = 0.04, for the value Ø = 400.The value of k obtained in the preceding analysis supports the assumption that D5Oba D50,rather than D5obQ,,k d50. Typically in gravel-bed rivers the value of 1),o of the armour layer isapproximately three times d50 from the subarmour sediment. If the assumption that Dsob d50is valid, then the value of k required to force an agreement between Wobs and would beapproximately three times as large than the value that is obtained when using D,0= D50. Thisvalue of k would correspond to a value of equal to about 0.12, which is well outside therange of values usually reported in the literature.Comparisons between the modelled and observed channel widths and mean depths are shown inFig 6.4(a, b). The agreement is good, the mean value of W0b I Wmod is 1.00 together with acoefficient of variation of 28.4%, and the mean value of Y0b/ Ymod is 1.05, and the coefficientof variation is equal to 14.5%.6.3.2 Effect of bank vegetationThe model was run with = 40° for the remaining 58 channels; those described as havingthick bank vegetation by Andrews, and the Vegetation Type II- IV channels from Hey andThorne. Reach 9074800 from Andrews’ data set was excluded from the analysis because astable channel width could not be obtained with çz5 = 40°.The results including the channels with the low density of bank vegetation are shown139graphically in Fig 6.5(a, b). Considerable scatter is evident away from the line of perfectagreement. Note the strong asymmetry in the scatter. The modelled channels are consistentlywider and shallower than the observed channels. The tendency for the modelled channels to bewider and shallower than the observed channels increases with the density of the bankvegetation. The channels with thin bank vegetation and those classified as Vegetation Type Iare scattered symmetrically about the line of perfect agreement as they were forced to do by theselection of the value for k in the previous section. The channels classified as having thick bankvegetation, and the Vegetation Type IV channels tend to lie furthest from the line of perfectagreement.The values of the ratios W0b/ Wmod and Yobs / Ymod for each Vegetation Type are summarisedin Table 6.1. Despite the wide scatter within each group it is evident that these two ratios showa systematic variation with the density of the bank vegetation with W0b/ W,,10d decreasing, andY0b3/ Ymod increasing with increasing bank vegetation density.Table 6.1. Ratios between the observed and modelled values of W and Y for fixed-slopeoptimization model with q = 400.Vegetation W,b / W I Y / Y, , No. ofType Mm Mean Max Mm Mean Max RiversI 0.64 0.97 1.77 0.70 1.05 1.29 13II 0.48 0.83 1.34 0.83 1.13 1.50 16III 0.28 0.72 1.16 0.90 1.28 2.05 13IV 0.25 0.59 1.00 0.98 1.48 2.16 20Thin 0.74 1.03 1.44 0.74 0.99 1.16 14Thick 0.17 0.58 0.82 1.11 1.31 1.44 9140The results obtained from the optimization modelling will now be compared to the those resultsobtained by Andrews (1984) and Hey and Thorne (1986) using empirical regime analyses.These researchers used regression analyses to derive empirical regime equations for thehydraulic geometry. Separate equations were obtained for each bank vegetation density. Thewidth equations were of the form:W=aQb/ (6.3)where Qbf is the bankfull discharge in Hey and Thorne, and is the dimensionless bankfulldischarge in Andrews. The exponent b shows little variation and typically takes a value ofapproximately 0.5. The coefficient a was found to vary with the density of the bank vegetationbecoming smaller as the bank vegetation density increased.The ratios formed by dividing the coefficient a from each of the regime equations, by thecoefficient a from the equation which represents the channels with the lowest density of bankvegetation is equal to W0b3 / Wunveg. where the subscript unveg refers to the unvegetated channelwidth. The ratio W0b / Wunveg is an index of the influence of the bank vegetation on the channelwidth. For example the value of a for the Vegetation Type I channels of Hey and Thorne is4.33, and the value of a for the heavily vegetated Type IV channels is 2.34. Therefore theVegetation Type IV channels are on average 2.34 / 4.33 = 0.54 times as wide as the weaklyvegetated Type I channels.141The ratio W0b / Wunveg is directly analogous to the preceding ratio W0b3 / Wmod where Wmod isthe calculated width assuming qS = 400. The values of W0b / Wunveg for the Andrews’ data setand the Hey and Thorne data set are summarised in Table 6.2 together with the values for W0bI Wmoj obtained from the modelling in this thesis with = 40°. Clearly there is reasonably goodagreement between the results of the optimization modelling with = 40°, and the empiricalregime analyses of Andrews and Hey and Thorne.Table 6.2. Summary of ratios of observed channel width divided by the unvegetated channelwidth. These ratios are an index of the effect of bank vegetation on the channel width. Thevalues from Thorne et a!. (1988) were calculated using their four regression equations and anassumed value of Wmod equal to 25m. See text for discussion.Vegetation W0b/ Wmod W0b/ Wunveg Wob / Wunveg Wobs / WmodType This Thesis Andrews (1984) Hey and Thorne Thome et a!.(1986) (1988)I 0.97 - 1.00 1.45II 0.83 - 0.77 1.20III 0.72- 0.63 0.96IV 0.59- 0.54 0.84Thin 1.03 1.00 - -Thick 0.58 0.79 - -Also in Table 6.2 are results from Thorne et a!. (1988) who tested the optimization model ofChang (1980) on the data set of Hey and Thorne (1986). The essential difference between theoptimization model developed by Chang, and the fixed-slope model from this chapter, is thatChang assumes a constant value for the bank angle Ofor a given channel, and does not considerbank stability.142A direct comparison between the results obtained from Chang’s optimization model and theresults obtained in this thesis is not possible as the values of Wmod were not given by Thorne etal. However an indirect comparison is possible. To account for the influence of the bankvegetation Thorne ci al. developed four regression equations, one corresponding to eachVegetation Type. This allowed Wmod obtained from Chang’s model to be corrected for theinfluence of the bank vegetation. Although the values of Wmo€i from Chang’s model are notknown, the four regression equations can be used to back calculate a value of Wmod for eachVegetation Type for a selected value of W0b3. The values of W0b / Wmod for a selected value ofWmod = 25 m are listed in Table 6.2.The results from the Thorne et al. analysis using Chang’s optimization model are significantlydifferent from the results of Hey and Thorne and in this thesis. This demonstrates that theinclusion of the bank stability constraint makes a very significant improvement in theperformance of the optimization model.6.3.3 Influence of ç6 on Channel GeometryThe preceding section indicates that the bank vegetation has a large influence on the bankstability, which in turn has a correspondingly large influence on the channel geometry. Theeffect of the bank vegetation on the value of will now be examined.One of the effects of the bank vegetation is to stabiise the bank sediment by binding of thegrains by the root masses. A simple method of accounting for this effect is by modifying the143value of which would take larger values for channels more strongly affected by the bankvegetation. This approach ignores the cohesive effects that would be introduced by the rootmasses, and is simply a device to account for the influence of the roots that can be incorporatedinto the noncohesive bank-stability constraint (Eqn 5.15).To demonstrate this effect the optimization model was programmed to run for a series of trialvalues of ç& for each of the channels from Andrews (1984) and Hey and Thorne (1986) untilthe value of Wmod was equal to W0b within a tolerance of ±1%. For example, Reach 13 fromHey and Thome has a value of Wobs = 18.4 m; when the optimization model was run with q5 =400 a value of Wmod = 74.4 m was obtained, which is over four times the observed value. Aftersuccessive trials it was determined that a value of = 73.10 forced an agreement betweenWmod and Wobs.The values of Ø obtained from this analysis depend on the value ofD50 assumed. Becausethe bank stability constraint (Eqn 5.15) is a function of two parameters, qS,. and D50 (thevalues of k, s, and yare assumed constant), the value of can only be determined if D5obO isknown, or in this case is assumed to equal D50. Using different values ofD50 will result indifferent estimates of Ø,.The maximum value of çS., obtained was 90.0°. From the bank stability constraint, Eqn. (5.15)the value of tan 90° is equal to infinity. However rounding-off errors in the computer program144result in a large real value being returned for tan 900.The values of ç5, obtained from the analysis are summarised in Table 6.3. Despite the largescatter within each vegetation grouping there is an observed increase in the mean values of q5with increasing bank vegetation density. This supports the general approach of accounting forthe effect of the bank vegetation by adjusting.The optimization model was rerun for all of the channels using the mean values of g fromTable 6.3 for each vegetation subdivision. The results are shown in Fig 6.6(a, b). There is areduction in the scatter when compared to the results shown in Fig 6.5(a, b) for = 40°. Usingthe mean values of for each Vegetation Type, for the complete set of channels the ratio W0b/ Wmod has a mean value of 1.00 and a coefficient of variation of 28.4%, and Yobs / Ymod has amean value of 1.05 and a coefficient of variation of 14.5%.Table 6.3. Summary of çt values obtained analytically for the data sets of Andrews (1984) andHey and Thorne (1986).Vegetation No. of çS valuesType Rivers Mm (°) Mean (°) Max (°I 13 21.1 42.0 52.8II 16 31.5 47.0 61.0III 13 34.2 53.0 72.5IV 20 39.7 60.1 90.0Thin 14 30.4 40.2 48.4Thick 9 45.7 55.7 67.6The residual scatter in Fig 6.6(a, b) can be attributed to several reasons. These include:1451. The bank vegetation categories are subjective, and a continuous gradation of the vegetationdensity exists within and between the vegetation types, and therefore a corresponding gradationin q5,. would be expected.2. The effect of the bank vegetation is not simply a function of the vegetation density, but isundoubtably dependent upon vegetation type, age, and rooting depth, as well as the thicknessof the overlying cohesive unit.3. Imbrication, packing, and cementing of the gravel by fine sediment is independent of theVegetation Type.4. The assumption that D5Obaflk = D50 is an additional source of uncertainty.6.4 COHESIVE BANK SEDIMENTThe model will now be formulated for cohesive bank sediments and tested on hypothetical andreal river data. The only change from the formulation for cohesive sediments is the bankstability constraint.As was discussed in Section 5.2 the bank stability constraint for cohesive sediment is composedof two individual constraints, namely the bank-height and the bank-shear constraints whichcorrespond to the erosion mechanisms of mass failure and fluvial erosion respectively.1466.4.1 Bank Stability RoutineThe bank stability routine satisfies the bank-height and bank-shear constraints consecutively.First the bank-height constraint is satisfied. The values ofH and Hrrt are initially assessed for 8= 900. For 6= 90°, and for subsequent values of 8, the stability number N3, is obtained fromstability curves such as those in Fig 5.3(a), and the value of Hcr,t is calculated from Eqn (5.2).The value H is by definition equal to the flow depth Y which is obtained by satisfying thecontinuity constraint for Qbf If the value ofH Hcrjt then the bank-height constraint is satisfied,and the routine then moves on to the bank-shear constraint. If the bank-height constraint is notsatisfied, the value of 8 is reduced until H = Hcrjt which gives the maximum bank angle 8 forwhich the bank-height constraint is just satisfied. For values of 9 < 8,,, the bank-heightconstraint is satisfied, but the channel will be less efficient with respect to transporting sedimentas was discussed in Section 6.2.5.Unlike the noncohesive case, the stability of cohesive banks with respect to fluvial entrainmentdoes not necessarily increase with a reduction of 8. Fig 6.2 shows the typical variation of Tbankwith 8 where continuity has been satisfied. As the value of 8 is reduced from the maximumvalue, Tbank increases to a maximum, and then decreases for low values of 9. Therefore exceptfor small values of 8, a reduction in 8 is accompanied by an increase in Tb07, and a reduction inthe stability of cohesive banks with respect to fluvial erosion. The bank stability routineassesses the bank-shear constraint for the value which satisfies the bank-height constraint,and if Tb0J r then the bank-stability constraint has been satisfied. If rit for 8,,,, thevalue of 8 is reduced until Thank = T. In practice it has been found that if the bank-147shear constraint is not satisfied for then the value of 8 which satisfies the constraint isabout 200 or less.The full data requirements for model testing were not available, however the model will betested on the published data from Chariton et a!. (1978) which lists many of the required inputvalues.6.4.2 Data From Chariton et aL (1978)Fourteen rivers from the Chariton et a!. (1978) data set which were described as having bankscomposed of fine sediment will be used. The rivers, together with the hydraulic geometry andother relevant information are listed in Appendix E. The value of d50 for the subarmoursediment which is used to calculate gb was not given in the data set, and as with thenoncohesive channels in Section 6.3, the value d50 = D50 / 3 will be used as an approximation.The value of ksbed was set equal to k5 calculated from the observed channel geometry. Unlikegravel-bed rivers with noncohesive gravel banks, those with cohesive banks may havesignificantly different values of and ksbank.The value of is unknown, so a value equalto 0.1 m was assumed for all channels. This assumption will not significantly affect the result asthe total channel roughness is determined largely by ksbed•The value of the unconfined compressive strength q was given for each channel. This value isan average of several samples from each river. The friction angle for the bank sediment was not148given and a value g’ = 25° is assumed for all channels which is a reasonable value for themoderately plastic silty soils described by Chariton et a!. (1978). With the values of qu and 0’,the value of c’ can be estimated from the following equation (eg. Spangler and Handy, 1982):c’=-tan(45__J (6.4)The values for the bulk unit soil weight were not given. A value of 7d 20.0 kNIm3 will beassumed for all channels for drained bank conditions, and rt 22.5 kN/m3 for saturated bankconditions. For the saturated bank conditions the value of ‘ will be modified by Eqn (5.6) toreflect the reduction in bank strength. This gives a modified value of the bank friction angle=(22.5-9.8)122.5 *250= 14.1°.The value ofN will be calculated from the stability curves in Fig 5.3(a). These stability curveswere approximated by piece-wise linearized segments over the range 0’ 6 90°. Althoughthese curves were developed for homogeneous, drained soils which are unlikely to exist inactual field situations, these curves will be used here for illustrative purposes. More realisticstability curves can be developed following field investigations. The stability curve for ‘ = 15°will be used for modelling the saturated conditions as it is sufficiently close to 14.10.The values for r were not given by Chariton et aL (1978), and it is not possible to directlyestimate the values for z. However in a manner similar to the estimation technique for for149the noncohesive channels in Section 6.3.3, the optimization model can be used to indirectlyestimate values forBecause the value of Trif is not known, initially the model will be run with a very large value forTcrjt (crit = 1000 N/rn2) to ensure that the bank-shear constraint is not influencing the solutionand the bank-height constraint alone is analysed. From this analysis the channels that arepotentially bank-height constrained can be determined. For the remaining bank-shearconstrained channels the values of rrjt can be estimated by forcing an agreement between Wmodand W0b5 as was done for Ø in Section 6.3.3 for the channels with noncohesive banks.The values obtained from the model assuming rrgt = 1000 N/rn2 are presented in Table 6.4 andFig 6.7(a, b). The model was run for both drained and saturated bank conditions. From thevalues of H / Hcrjt it is evident that channels 5, 7, 10, 12, and 13 (and 8 for saturated bankconditions) are constrained by the bank height. The remaining channels are bank-heightdegenerate. (The term degenerate is used in optimization to denote constraints that are notactively constraining the solution.) The bank-shear constraint is degenerate for all of the 14channels due to the large imposed value of Thrit. The bank-height constrained channels will nowbe considered. Note that different solutions are obtained for the drained and saturated bankconditions. The saturated banks are inherently less stable than the drained banks due to theincreased unit weight of the soil, and the increased pore pressures. This effect is well knownand is discussed in Section 5.2.1.1. The second feature evident from this analysis is that thebank-height constrained channels tend to be the larger channels with values of Qbf greater than150Table6.4.Resultsfromtheanalysisof datafromCharitonetaL(1978).Themodellingwith=1000N/m2forcesthebank-shearconstraint tobedegenerate, andonlythebank-height constraint isanalysed.TheestimatedvaluesofwereobtainedbyvaryingtrialvaluesofVcrjtuntilanagreementwasobtainedbetweenW,,,andW0b3.Reachnumbers10and13appeartobebank-height constrainedandthereforenoestimateof rrircouldbeobtained.ReachObservedModelledwith=1000N/rn2Estim-BankNo.DrainedSaturatedatedVegeWW‘H/Vbak/WjH/VbankrFtation(m)(m)(m)(m)HcrgtTent(m)(m)/rnjt(N/rn2)Type517.61.7911.92.361000.0412.22.241.000.0532.8T731.01.7712.53.081.000.0215.82.591.000.0210.0G828.71.6315.12.630.950.0215.02.601.000.0211.9T1019.02.4717.12.601.000.0220.72.261.000.02-T1239.32.6420.53.891.000.0225.83.291.000.0210.8G1359.44.1948.04.831.000.0263.34.001.000.02-G117.41.7811.72.380.650.0411.72.380.890.0429.7T214.00.736.11.180.500.036.11.180.680.0313.9G39.80.736.10.980.460.056.10.980.620.0532.4T413.71.348.91.820.380.038.91.820.520.0322.3T619.01.369.32.190.330.099.32.190.460.0950.0T916.70.696.81.260.530.026.81.260.720.029.4G115.20.653.30.930.280.013.30.930.380.017.0G1419.51.6712.72.230.610.0212.72.230.830.0211.1about 65 m3/s. This is to be expected as the smaller or lower discharge channels will developvalues ofH that are lower than Hcrjt.The remaining channels which are bank-height degenerate have the same optimal solution fordrained and saturated bank conditions. It is assumed that the value of Tcrjg is unaffected by banksaturation. For these channels the modelled values are all narrower and deeper than theirobserved counterparts. This suggests that the values of Thank from the modelled channels aregreater than can be sustained by the observed channel. By reducing the value of rrgt, andtherefore reducing the resistance of the banks to withstand applied shear, a wider and shallowerchannel will result. In this way the modelled channels can be brought into agreement with theobserved geometries, and estimates of zrjt obtained.Also note that channels 5, 7, 8, and 12, which appear to be bank-height constrained when thevalue r = 1000 N/rn2 is used, are also much narrower and deeper than their observedcounterparts. For smaller values of z,rit, these channels will become bank-shear constrained andthe value of Wmod can be brought into agreement with W0b. Therefore only channels 10 and 13are probably truly bank-height constrained as the values of W0b and Tobs lie within the rangedefined by the limiting cases of fully saturated and fully drained bank conditions.The value of TCrIt was varied until Wmod = W0 within a tolerance of ±1% for all channels withthe exception of numbers 10 and 13. These estimated values of are also shown in Table 6.4.The values range between 7.0 - 50.0 N/rn2,with a mean value of 20.1 N/rn2.1526.4.3 Effect of Bank Vegetation.Chariton et a!. (1978) have categorised the rivers in their data set on the basis of bankvegetation into channels with grassed (G) or treed (T) banks. The values of TCrjt obtained fromthe previous analysis show a strong influence by the bank vegetation. The mean value of Thrit forthe channels described as having grassed bank vegetation is 10.4 N/rn2, and 29.8 N/rn2 for thetreed banks.The values of Tcrjt are plotted in Fig 6.8 and this indicates a strong division between the twoVegetation Types. With the exception of one treed channel all of the treed channels plot above20 N/rn2, and all of the grassed channels below 14 N/rn2. This result suggests that the effect ofthe bank vegetation is to increase the value of Thrit either by binding the sediment by the rootmasses, or conversely by affording protection of the bank and effectively reducing the value ofThank acting on the bank sediment. Furthermore the bank vegetation may bind the bank sedimentas to increase the stability of the banks with respect to mass failure. In this way the roots act asinternal reinforcement, and have the effect of increasing the effective values of c’ and qS’ abovethe values obtained from the analysis of small samples of the bank material.Clearly there is a need for additional field work to identify the role of vegetation in stabiisingthe cohesive banks, a conclusion which also applies to noncohesive banks.6.4.4 Discussion of March etaL (1993)A recent paper published by March et a!. (1993) dealing with bank stability will now be153discussed as it highlights several important aspects of the above analysis. In this study stabilitycurves analogous to Fig 5.3(b) were constructed for the Long Creek drainage basin in northernMississippi using averaged values of c, Yt, and qS which were measured during previous fieldsurveys (Thorne er a!., 1981). Stability curves were developed for drained bank conditionsusing the measured averaged values of c, , and and for “worst case” soil conditions whichreflects saturated bank conditions. For values of 9 greater than 600 the Osman - Thorne slabfailure analysis (Osman and Thorne, 1988; Thorne and Osman, 1988) was used to determinethe stability curves that correspond to F = 1.0, from which Hcrjt can then be determined for anyvalue of The Osman - Thorne slab failure analysis is similar to the Culmann analysispresented in Fig 5.1 and Eqns (5.3) and (5.4), but is modified to include the effect of tensioncracks. For values of 0 less than 60° the stability curves were developed using the Bishop(1955) simplified method of slices. The stability curves from March et a!. (1993) arereproduced in Fig 6.9.With the exception of one data point, the points all plot below the “worst case” stability curvewhich represents saturated bank conditions at failure and corresponds to the CASE Ill examplein Chapter 5 (Fig 5.4). Three of the points lie on, or very close to the saturated curve. Thisdistribution of the observed H - 9 combinations is consistent with the bank-height constraint,Eqn (5.7), in that the bank heights must be less than or equal to the critical bank height.However note that 13 out of the 16 data points which plot in the stable field actually failedduring the previous winter period. Since these channels are all stable with respect to the “worst154case” saturated bank conditions, this may at first seem puzzling. However the bank-shearconstraint has not been directly addressed in the March et a!. study. A bank can only beconsidered to be stable when it satisfies both the bank-height and the bank-shear constraints.The role of the bank-shear constraint will now be addressed and it will be shown that this canexplain the failure of the supposedly stable banks given in March et a!.The two erosion mechanisms, mass failure and fluvial erosion, do not operate independently.The interrelation of the two mechanisms is illustrated in Figs 6. 10(a-d). These photographs areof the banks of a small creek flowing across a beach at low tide. The bank material is wellsorted fine to medium sand which is normally non-cohesive, but has developed apparentcohesion due to the moisture content and surface tension phenomena.Fig 6.10(a) was taken at an upstream location where the vertical, stable channel bank isapproximately 480 mm high. The flow at this location was not actively undercutting the bank,ie. Tbank There was no evidence of mass instability of the bank even when an additionalloading of 75 kg (the author’s body weight) was applied, and therefore H < Both bankstability constraints are satisfied.Figs 6. 10(b-d) are a sequence of photographs that were taken at a site approximately 20 metresdownstream from the location in Fig 6.10(a). This downstream location was experiencing rapidlateral channel shifting, and the bank is unstable. The time between each photograph is of theorder of 10 seconds. At this location the channel bank is approximately 100 mm high, much less155than at the stable upstream location.In Fig 6.10(b) the shallow, supercritical flow is beginning to undercut the vertical bank. In Fig6.10(c) the bank has been destabilised by undercutting and mass failure results. In Fig 6.10(d)the failed block has been washed away and a new cycle of undercutting is about to commence.If a bank stability analysis was conducted on the bank geometry in Fig 6.10(d), which may bepreserved at low flows, the bank-height constraint would appear satisfied as the bank is muchlower than the stable, vertical bank upstream. However the bank at this location is clearlyunstable as the bank-shear constraint is not satisfied, that is Zbank > z, and bank is beingundercut which leads to eventual mass failure. Therefore the primary erosion mechanism is infact fluvial erosion of sediment from the toe of the bank, and the observed mass failure is only asecondary effect.March er a!. do not explicitly recognise the requirement of the bank-shear constraint, althoughthey do state that the observed H - 8 combination may have been different at the time of failureand that channel migration may have over steepened the banks until failure occurred.Alternatively bed scouring may have increased H until Hcrit was exceeded.In conclusion it has been demonstrated that a bank composed of cohesive sediment can only beconsidered to be stable when both the bank-height and bank-shear constraints are satisfied.1566.5 BANKFULL MODEL: VARIABLE-CHANNEL-SLOPEThe model will now be modified to allow for variable channel slope. In this formulation thechannel slope S will be treated as a dependent variable. The sediment discharge capacity of thechannel at the bankfull discharge Gbf will be used as an independent variable. The continuity andbank stability constraints remain unchanged from the previous formulation in Section 6.2.There is an additional requirement for a bedload constraint:ed b = Gbf (6.5)where gb is calculated from the Eqn (6.2), and the value of Gbf is prescribed as an independentvariable.The objective function is modified to:maxf(Jd,],O,S)=1l (6.6a)where:Gb(6.6b)PQbfsince Gb1 p, and Qbf are all prescribed in this formulation, the maximization of i is equivalent to157a minimization of S.The flowchart for the modified optimization model is shown in Fig 6.11, and the source codefor the computer program in Appendix B. The only difference between Fig 6.11 and theflowchart for the previous formulation in Fig 6.1 is addition of the bedload constraint whichmust be satisfied before assessing the objective function. The bedload constraint is satisfied byvarying S for trial values OfPbed.An example of the variation of i with Pbed is shown Fig 6.12. Each point on the solution curvesin Fig 6.12 satisfies the continuity, bank stability, and bedload constraints. As with the fixedslope model the objective function curve is smooth with only a single global optimum. Alsoshown in Fig 6.12 is the variation of 5, which indicates the optimal solution corresponds to aminimum slope condition.6.5.1 Noncohesive Bank SedimentThe variable slope optimization model will now be run using the data from Andrews (1984) andHey and Thorne (1986). The value of Gbf was calculated from the observed geometry usingEqns (6.2) and (6.5). All other input values and assumptions are unchanged from Section 6.2.The program was run for all the data for Ø = 400 to determine the effect of the bank vegetationon the channel geometry, including the variation of S. As in Section 6.3.2 modelled valuesobtained assuming çS = 40° are assumed to represent the unvegetated channel dimension. Theresults are summarised in Table 6.5.158The results summarised in Table 6.5 indicate that as the density of the bank vegetation increasesthe observed channels become progressively narrower, deeper, and less steep. The observedand modelled channels summarised in Table 6.5 have the same sediment transporting capacities,therefore if more resistant banks permit narrower and deeper channels to be stable, the bedloadconstraint is satisfied by reducing the channel slope.Table 6.5. Ratios between the observed and modelled values of W, Y, and S for variable slopemodel with çb 400.Vegetation W0b5I Wmod Y0b5I Ymoci Sobs / Sm€jjType Mm Mean Max Mm Mean Max Mm Mean MaxI 0.59 0.98 1.84 0.64 1.05 1.43 0.84 0.99 1.26II 0.43 0.80 1.35 0.81 1.19 1.71 0.81 0.93 1.08III 0.23 0.70 1.17 0.88 1.39 2.61 0.66 0.90 1.03IV 0.20 0.55 0.97 1.01 1.67 2.60 0.68 0.85 0.96Thin 0.70 1.03 1.51 0.75 1.00 1.24 0.88 0.98 1.10Thick 0.16 0.54 0.78 1.16 1.37 1.59 0.86 0.95 1.38This effect of the bank stability will now be demonstrated using Reach 13 from Hey andThorne. The optimization model was run for Reach 13 using the values of ç& = 40° and c =73.1°, the latter value was found in Section 6.3.3 to force an agreement between the modelledand observed channel widths. The output together with the observed channel dimensions areshown in Table 6.6. Note the large influence that exerts on the channel geometry, includingthe channel slope.The results of this analysis indicate that W is most sensitive to variations in the bank stability,159followed by Y, with S the least sensitive. For the combined Vegetation Type IV and thickvegetation channels the channel widths, depths and slopes are in the order of 0.5, 1.6, and 0.9times their respective unvegetated channel dimension.Table 6.6. Effect of SL on Reach 13 from Hey and Thorne (1986). The effect of increasing thevalue of is for the channel to become narrower, deeper, and less steep.Channel Surface Width Mean Depth Channel Slope(m) (m)Observed 18.4 1.14 0.0133Modelled qS=40° 91.2 0.44 0.0183Modelledq5,=73.l° 18.4 1.18 0.0133By forcing an agreement between Wmôd and W0b, the values of Ymoci and Vobs must also showclose agreement as once Wmod is fixed, the value of Ym,ci is constrained by continuity. Similarlythe value of S,,, is determined largely by the bedload constraint, when the channel width isforced to agree with the observed value the modelled channel slope must also show closeagreement with the observed value. In other words when one of the values of any of thedependent variables is fixed, there is only one combination of the remaining dependent variablesthat can satisfy the constraints and fulfil the objective function.Recall from Section 6.3.3 that when the value ç5,. = 40° was used to model reach 13 with thefixed-slope-model, the computed surface width was equal to 73.5 m, in contrast to 90.0 m withthe variable-slope-model. The larger value of Wm obtained with the variable-slope-model is160due to the increase in the channel slope to S = 0,0183 from S = 0.0133 in order to satisf’j thebedload constraint. This illustrates the interrelationship between the dependent variables, andthat the adjustment of one variable such as W, cannot be viewed in isolation from theadjustments of the other dependent variables.6.5.2. Cohesive Bank SedimentThe model was run using hypothetical data; firstly Qbf was varied, and the value of Gbf heldconstant; then the value of Qbf was held constant and Gbf varied. For both analyses thefollowing values for the independent variables were used: D50 = 0.075 m, d50 = 0.025 m, ksbed =0,1 m, ksb nk = 0.1 m, c’ = 10 kN/m3, = 20 kNIm3, ç’ = 25°, and rjt = 25 N/rn2. When Qbfwas varied the value Gbf = 5 kg/s was held constant, and when Gbf was varied the value Qbf =100 rn3/s was held constant. The results are presented in Figs 6.13 and 6.14. The key resultfrom this analysis is the change in the active bank stability constraint.For small values of Qbf approximately less than 100 m3/s the optimum geometry is bank-shearconstrained and bank-height degenerate, and those in excess of 250 rn3/s or so the channelbecomes bank-height constrained and bank-shear degenerate (Fig 6.13(c)). Between about 100to 250 m3/s the channel is actively constrained by both the bank-shear and the bank-heightconstraints.When Qbf was held constant and Gbf varied, a similar result was obtained whereby the channelswith values of Gbf less than about 3.0 kg/s are bank-height constrained and bank-shear161degenerate, while those in excess of 5.0 kg/s are bank-shear constrained and bank-heightdegenerate (Fig 6.14(c)). Within the range between 3.0 to 5.0 kg/s the modelled channels bothbank-height and bank-shear constrained.To explain this change in the active constraint, first consider the case where Qbf is held constantso the results do not become confused by the different channel sizes which are associated withdifferent Qbf values. For low values of Gbf the values of S, ; and Thank are at their minimum, andY is a maximum (Fig 6.14). The channels corresponding to the low values of Gbf are bank-height constrained due to the relatively large values of Y, and the small values of Thank.As the slope increases with increasing Gbf, the channels become shallower due to the increasein both the mean velocity, U and W. The depth decreases at a slower rate than the value of Sincreases, and therefore the values of r and Thank both increase with increasing Gbf Thiscontinues until the value of Thank becomes equal to Tcrjt, at which point the channels becomebank-shear constrained. The value of S continues to increase more quickly than the depthdecreases, and therefore the channel remains bank-shear constrained.When Qbf is varied and Gbf is held constant, the values of S decrease with increasing Qbf (Fig6.13 (c)), which is a feature of natural rivers that is well known from field observations. Forsmall values of Qbf the channels are bank-shear constrained due to the combination of the high Tand Thank values that are associated with steep channel slopes, in addition to the small values ofH which are well below Hcr:t.162As Qbf is increased the value of Y increases and S decreases, however the values of v and rnkremain initially constant because the channel is bank-shear constrained. Eventually the value ofY (and by definition H as H = Y for Qbf) increases to the point where the channel becomes bank-height constrained. For increasing Qbf beyond this point, Y tends to increase at a lower rate thanthe decrease in S and therefore the shear stress values decrease, and this ensures that thechannel remains bank-height constrained.The variation of 6 has not been considered up to this point. In Fig 6.13 the value of Y is equalto 2.95 m at Qbf= 100 m3/s. At this point H Hrit. Yet for larger Qbf, the values of Y increaseand approaches 12.7 m for Qbf 2,500 m3/s where H is still equal to Hcrjt. This is possible dueto the change in 8 and the influence that this has on the value of Hcrjt. At Qbf = 100 m3/s thevalue of 6 = 90° and Hcrjt = 2.95 m. For Qbf = 2,500 m3/s the value of 6 has decreased to 44°,and has increased to 12.7 m. The influence of 9on the value ofHcrit is evident from Fig 5.3and Eqn (5.2).The general conclusion from the modelling in this section is that channels with large values of Stend to be bank-shear constrained. The large S may be associated with large sediment loads, orsmall values of Qbf, or a combination of both. Bank-height constrained channels tend to beassociated with low S values.1636.6 SUMMARYBank stability has been shown to exert a strong control on the optimal geometry of alluvialgravel-bed rivers. The bank stability constraint was formulated for noncohesive channel banksand the theory tested on the published data of Andrews (1984) and Hey and Thome (1986).The results indicate that the bank stability procedure significantly improves the modelperformance. The results are in good agreement with the results that Andrews and Hey andThorne obtained from empirical regime analyses. The influence of the bank vegetation appearsto stabilise the bank sediment, allowing the banks to withstand higher shear stresses. Invegetated channels this results in channels that are narrower, deeper and less steep than theirunvegetated counterparts. The change in W is the largest, followed by Y, and the smallest is inS.It was proposed that the effect of bank vegetation on channels with noncohesive banks can berepresented by,which was found to increase consistently with the bank vegetation density.Model formulations were also developed for channels with cohesive banks. The model wastested on data from Charlton et al. (1978). From this analysis it was shown that channels can beeither bank-shear or bank-height constrained. The values of TCr,t was shown to be higher inchannels with treed banks, than for there grassed counterparts.Modelling results using hypothetical data indicate that bank-shear constrained channels areassociated with large values of S which can result from large imposed sediment loads, and/or164small values of Qbf Conversely bank-height constrained channels tend to be associated withlower values of S which can result from large values of Qbf, low sediment loads, or acombination of both.In general only one bank stability constraint is active in channels with cohesive banks, andrecognition of this is essential when assessing channel stability.165Figure 6.1. Flow chart for the bankfull:fixed-slope optimization model.16614e(°)Tbank bedFigure 6.2. Variation of Thank and Tbed with 6 Each point on the curves satisfies the continuityconstraint for a constant value ofPhd.2018oJE161220 40 60 801670.2 401Figure 6.3. Variation of i and 6 with Fbed. Each point on the solution curves satisfies thecontinuity and bank-stability constraints.C.)C.)tCCl,0.180.160.140.120.10.080.060.0430102010010020 40 60 80bed (m)168(a) 1005020Va,:i0V52.(b)Ea,UVa,a)V0PerfectType I Thin Agreement Io •Figure 6.4. Comparison of Modelled and observed values of W and Y for rivers withnoncohesive banks and low densities of bank vegetation. The data points denoted “thin” arefrom Andrews (1984), and those denoted “Type r’ are from Hey and Thorne (1986). Themodelled values were calculated using the bankfull: fixed-slope optimization model.Observed Width (m)0.3 0.5 1Observed Depth (m)169(a) 100ro20.1o0522(b)2-Ea)DE0.5a)oO.30.20.101 0.2 0.3 0.5 1 2 3Observed Depth (m)Ripe I Type II Type III Type IV Thin Thick Perfect IAgreemento * CFigure 6.5. Comparison of Modelled and observed values of W and 7 for rivers withnoncohesive banks and variable densities of bank vegetation. The data points denoted “thick”are from Andrews (1984), and those denoted “Type r’ to “Type IV” are from Hey and Thome(1986). The modelled values were calculated using the bank1ill: fixed-slope optimizationmodel.Observed Width (m)0.,. 00D*•4.00•0ci0 00ci• IIII———I170(a)(b)-504-20.V032E4-a)0.5a)Vo 0.30.20.101Type I Type II Type III Type IV Thin Thick PerfectAgreement0 *Figure 6.6. Comparison of modelled and observed values of W and 7 for rivers withnoncohesive banks and for all categories of bank vegetation density using the mean valuesfrom Table 6.3. The data points denoted “thin” are from Andrews (1984), and those denoted“Type r’ to “Type IV” are from Hey and Thorne (1986). The modelled values were calculatedusing the bankfull: fixed-slope optimization model.1100Observed Width (m)0.2 0.3 0.5 1 2 3Observed Depth (m)171(a)(b)7060-c0.20CCu0,0)V04030100 0 1054210070Figure 6.7. Comparison of Modelled and observed values of Wand Y for rivers with cohesivebanks from Chariton et a!. (1978). The modelled values were calculated using the bankfull:fixed-slope optimization model for values of rrjt = 1000 N/rn2.20 30 40 50 60Observed Width (m)31 2 3 4Observed Mean Depth (m)Bank-Shear Constrained Bank-Height ConstrainedDrained Banks4Bank-Height Constrained PerfectSaturated Banks Agreement1726050 TC.jE40Z T30----T-T20 -Mean value for grassed banks10 T----G eT G .0— I I I I I I I I I0 2 4 6 8 10 12 14Reach NumberFigure 6.8. The values of z calculated from the data of Chariton ci a!. (1978) for treed (T)and grassed (G) channel banks. The mean values for the treed and grassed banks are indicated.173LUwII — C STASLEDURINGI992Figure 6.9. Bank stability curves from March et a!. (1993). The lower “worst case” curvecorresponds to the Case III example from Fig 5.4, and the upper curve for average soilconditions corresponds to Case I from Fig 5.4.BANK ANcLE (DEcEES)174(a) : Z4 ir .. -‘•I:..A.__ ‘—.(b)Figure 6.10 (a) and (b). Photographs of stream bank erosion from a small creek flowing acrossa beach at low tide. The lens cap in the photographs is for scale and is approximately 45 mm indiameter. See text for discussion.175a‘a(c)(d)Figure 6.10 (c) and (d). Photographs of stream bank erosion from a small creek flowing acrossa beach at low tide. The lens cap in the photographs is for scale and is approximately 45 mm indiameter. See text for discussion.r --:a.-a. -•‘176Figure 6.11. Flow chart for the banlcfull: variable-slope optimization model.N1770.055 0.00477Figure 6.12. Variation of 77 and S with Pbed. Each point on the solution curves satisfies thecontinuity, bank stability, and bedload constraints.C)C)C)0cjHC)C)clD0.050.0450.040.0350.030.0250.00350.0030.00250.0020.00151000 20 40 60 80bed (m)178(a) 100 141280100CD60 8.3: 40 6420 20 0(b) 0.0035 450.003400.0025a)o. 0.002o0.0015 30 30.001250.00050 20Bank-Shear Constrained Duai Constrained : Bank-Height Constrained(c)_______________________________1 25-i0.g -.-----------0.8 H/H - 200.7bank’ crlt0.6-153tbank0.50.450 100 200 500 1,000 2,000Q (me/sec)Figure 6.13. Variation of the optimal values of selected dependent variables as a function ofQbf (a) Wand Y. (b) S and i (c) HIH0rjt, Thank / Zcrit, and z. The value of Qbf is constant andequal to 100 m3/s. The active bank-stability constraint is indicated in Fig 6.13(c).50 100 200 500 1,000 2,00050 100 200 500 1,000 2,000179(a) .101 2 51020501(b) 0.004 500.0035 - Slope0.0005 250.5 1 2 5 10 20 50DuConsqedBank-HeghtConstrairied Ban honsJnedG bf (kglsec)Figure 6.14. Variation of the optimal values of selected dependent variables as a function ofGbf (a) W and Y. (b) S and z (c) H I Thani I Tcrit, and Thank. The value of Qbf is constant andequal to 100 m3/s. The active bank-stability constraint is indicated in Fig 6.13(c).180CHAPTER 7FULL MODEL FORMULATION7.1 INTRODUCTIONIn this chapter the final formulation of the optimization model is presented. The equilibriumgeometry is modelled using the full range of flows, and not just the bankfull discharge.Furthermore the composition of the bed surface is permitted to adjust by using the modifiedParker (1990) surface-based bedload transport relation.7.2 MODEL FORMULATIONThe full model formulation is similar to the bankfull:variable-slope model presented in Chapter6 except that the sediment transporting capacity of the channel will be estimated using theParker (1990) surface-based bedload transport relation, rather than the Einstein-Brownrelation. The Parker surface-based relation was modified in Chapter 4 to apply to a naturalchannel with variable flows. This relation can be inverted to calculate the composition of thebed surface, or armour layer. The mathematical formulation of the full model is presentedbelow.1817.2.1 Independent VariablesIn this analysis the independent variables whose values are presumed to be known are the fillrange of flows and their durations (which can be represented by a flow-duration curve), thegrain size mixture of the imposed bedload sediment, Gb, kjb& ksbk, the bank stability variablesD5Obank, $, for noncohesive banks, and c 6’, Yt, and for cohesive banks, and the constantsg, p, s, v, yThe flow-duration curve and the grain size distribution are subdivided into a finite number ofintervals which are assigned representative values. The flow-duration curve is divided into mintervals and the grain size distribution curve into n intervals (See Figs. 7.1 and 7.2). Becauseboth the flow-duration curves and the sediment size distributions are typically approximatelylog-normal (Shaw, 1988 p. 276; Parker, 1990) the representative value for each interval isgiven by the geometric mean:Qj=JQ,1.Q (7.1)D1=JD’.D (7.2)where Q’ and D are the geometric mean values of intervals i and j respectively, and thesuperscripts 1 and u indicate respectively the lower and upper bounds of the interval. The valuesof i and j range from 1 to m and 1 to n respectively. For each interval the proportion of thetotal is determined. For the flow-duration curve the proportion of the total flows which fallwithin a specified discharge interval is denoted by i. Similarly for the grain size distribution182curve, the proportion of the sediment by volume, or the fraction content within a specified grainsize interval is denoted F. Summations of i and F over all i andj respectively are both equalto 1.0.The imposed bedload, Gb, represents the total mass of sediment in excess of 2 mm in diameter,which is supplied to the channel reach over some significant duration. In this chapter thesignificant duration is assumed to be one year, and therefore the Gb represents the mean annualbedload supply. The units of Gb will be kg/y, and the conversion factor, T will have a valueequal to 365.25 * 24 * 3600 = 31.56 X 106 s/y (see Section 4.6).Typically the value of the bankfull discharge Qbf is considered to be an independent variable.This was questioned in Chapter 3 where it was suggested that Qbf may be considered adependent variable. This will be investigated in Section 7.6.7.2.2 Dependent VariablesThe primary dependent variables which are to be solved are Pbed, Pbank, O S, and the grain sizedistribution of the bed surface which will be represented by D50. Other secondary dependentvariables which include Rh, (J W, and Y are readily obtained once the primary dependentvariables have been determined.1837.2.3 Objective FunctionThe objective function in the full-model formulation is give by:maxf(Jd,Pb,6,S,DSO) = (7.3a)where the coefficient of efficiency i is given by:Gb Gb77=m = — (7.3b)Tpp,S ‘°where the mean annual flow rate with units of m3/s. The total stream power which isrepresented by pQS for steady, uniform flow is modified in the denominator of Eqn (7.3b) forthe full range of flows as given by the flow-duration curve (Fig 7.1). The total stream powerexpended over one year per unit channel length is given by I’ p S.Since the values of Gb and are imposed, a maximization of i is equal to a minimization of S.7.2.4 ConstraintsThe constraints for the full model formulation are given below.7.2.4.1 ContinuityThe continuity constraint is unchanged from previous formulations and requires that thedischarge capacity of the channel is equal to the value of the bankfull discharge, Qbf:184UA=Qbf (3.24)where A is the channel cross-sectional area at bankflill, and U is the mean channel velocity atbankfull which is given by:(3.35)the value offis given by Eqn (3.19).7.2.4.2 BedloadThe bedload constraint requires that the bedload transporting capacity of the channel be equalto the imposed sediment load Gb:Fg=G (4.54)where:g=rspq (4.55)where qb is given by Eqn (4.53). The value qb is the average volumetric bedload transport rateper metre channel width, and has units of m2/s. The units of Gb and gb are kg/y and kg/y/m,respectively.The bedload constraint implicitly contains the “equal mobility constraint” which requires that allgrain sizes be transported equally over the one year duration (see Chapter 4).1857.2.4.3 Bank StabilityThe bank stability constraint can be formulated for both noncohesive and cohesive banks.The noncohesive bank stability constraint is given by:y(s—i)D5kfl1_s: (5.15)The constant k was evaluated from field data in Chapter 4 as 0.048.The cohesive bank stability constraint is comprised of two components namely the bank-heightconstraint, and the bank shear constraint, which are, respectively:HH1 (5.7)rbflk 1crit (5.8)7.2.5 Optimization SchemeThe optimization flowchart is shown in Figure 7.3 and is unchanged from Fig. 6.11. Howeverthe bedload constraint is substantially more complex than previous formulations and is shown inFigure 7.4.186In order to calculate qb the values of the geometric mean grain diameter of the bed surface Dsg(Eqn 4.35) and the standard deviation on the sedimentological phi scale o (Eqn 4.37) must beknown. However these are dependent variables whose values are not known in advance.Therefore an iterative scheme is necessary which is initiated with trial values. The values of Fare calculated which satisfy “equal mobility” using Eqn (4.52) from which updated values ofDsg and o are obtained. This is repeated until convergence is attained on the values ofD ando. From the values of13. at convergence, the value ofD50 can be readily determined.7.3 VARIABLE FLOWSThe Parker surface-based bedload transport relation was modified in Chapter 4 toaccommodate the variable flows which are inherent in natural river systems. The armour layerdevelops such that “equal mobility” applies to the total load transported over a significantduration (which will be taken as 1 year) as a result of the total range of flows. The sizedistribution of the annual transported sediment load is constrained to be equal to the subarmoursediment.To illustrate the process of modelling a range of flows a hypothetical channel with noncohesivebanks will be examined. The flow-duration curve presented in Fig. 7.1 will be used as input.The mean annual flow, , is equal to 17.5 m3/s. The grain size distribution of the bedload andsubarmour sediment is shown in Fig. 7.2, and both have values ofd50= 0.025 m and o = 0.43.The value ag is the geometric standard deviation which is given by:187ir (Dl2usg=a[logLJj F (7.4)The values of and ksbk will both be set equal to 0.40 m, D5Oba,, set equal to 0.07 m whichis about the anticipated value ofD50, and ç = 400. The total annual imposed sediment load is2.5 X 106 kg/year. The value of the bankfiill discharge, Qb., is set equal to 85 m3/s which isequalled or exceeded 1.7% of the time. This value of Qbf is not arbitrary but it will be shown inSection 7.6 that it corresponds to an optimum. The flows that exceed Qbf are set equal to Qbffollowing the assumption of an infinitely wide flood plain. This assumption is discussed inSection 3.3.The optimum values of W, 7, S, and D50 obtained from modelling are 37.6 m, 1.16 m, 0.00426,and 0.073 m respectively, and the bankfull value of is equal to 0.043. All of thesemodelled values are typical for natural gravel rivers with little bank vegetation (eg. see regimeequations and field data from Hey and Thorne, 1986).The sediment transport rate Gb1’ for each value of Q’ is shown in Fig. 7.5(a). Here the prime (‘)is used to denote the rate in kg/s, and Gb, without the prime denotes the total transported in oneyear, or the rate in kg/year. Note that the maximum value of Gb,’ at bankfiill is only 1.89 kg/s.In Fig. 7.5(b) the total load transported by each flow interval is shown. The lowest flow interval(which occurs 75% of the year) transports a negiigible volume of sediment (about 4 kg), while188the flows that equal or exceed Qbf, transport over 40% of the annual sediment load over 1.7%of the year, which is equal to about 6 days.The variation of the median grain diameter of the transported sediment diOload, with discharge isshown in Fig. 7.6 and illustrates the definition of “equal mobility” as applied herein. The valueof dSOJoad varies with discharge; the lower discharges are associated with values of d5Oload lessthan d50, while higher discharges result in values of d5Oload greater than d50. The net result is thatover the entire year the value of dSOIoad for the total annual bedload is equal to d50. Note thatbelow a discharge of about 35 m3/s the median grain diameter is constant. For these lowdischarges the bed surface is essentially immobile (Fig. 7.5(a)), and in most natural rivers anysediment transport would be restricted to “sandy throughput load”. In Oak Creek the sedimenttrapped for discharges less than the threshold value of about 1 m3/s had a mean grain diameterof about 2 mm (see Fig. 4.12). Once this threshold value was exceeded the value of d5Oloadshowed a consistent increase with discharge.7.4 ADJUSTMENT OF THE BED SURFACE COMPOSITIONIn previous model formulations (Chapter 6; Millar and Quick, 1993; Miflar, 1991; as well asWhite et a!., 1982; Chang, 1979, 1980) the armour layer grain size distribution as representedby D0 was considered to be a fixed independent variable. In the full-model formulation thegrain size distribution of the armour layer is able to adjust by using the Parker (1990) surface-based transport relation. The value ofD50 will therefore now be treated as a dependent variable.189The degree of armour layer development will be represented by the value of 650 which is theratio of the median armour, to the median subarmour grain diameters:(4.25)The value of 8, observed in natural gravel rivers ranges from 1.0 with no armour layerdevelopment, to greater than 6. Most rivers have values of 650 around 3 when measured at lowflows. It is assumed in this thesis that the composition of the armour layer which is observed atlow flows is maintained at high discharges (Sect 4.4.1). In addition to D50, a characterisation ofthe armour layer should also include a measure of the dispersion of grain sizes about the medianvalue, however this is of secondary importance and will not be considered herein.7.4.1 Sso - T*D50 Solution CurvesParker (1990) has shown that for an imposed bedload grain size distribution a solution curvecan be calculated to predict the variation ofD with q5sgo. In this section 5o - T.D50 solutioncurves, which are analogous to Parker’s, will be developed to demonstrate the adjustment ofthe armour layer. The 850-TD50 solution curve for the Oak Creek sediment is shown in Fig.7.7. For values of less than about 0.03, which is normally considered the “threshold”value for gravel rivers, a static, immobile armour layer is developed. For values of r’D50 greaterthan 0.03 a mobile armour is developed (Parker, 1990). The value of 8o is seen to decrease asT’D50 increases, and eventually disappears for very large values of V*D50 in excess of 0.3.190The value of ö observed at Oak Creek is about 2.2 which corresponds to a value of v”D50equal to about 0.06. The bed surface at Oak Creek, as with any natural river, is subject to awide range of shear stresses. Therefore the value r*D50 = 0.06 for Oak Creek can be consideredto be the “dominant” or “channel forming” dimensionless shear stress which is analogous to thedominant discharge concept. This dominant dimensionless shear stress is likely to correspond tobankfull, or near bankfhll values as most of the sediment transport occurs at these discharges.Once a channel has reached equilibrium it is presumed that the composition of the armour layerdoes not change appreciably despite the daily and seasonal variation in the flows and shearstresses.The grain size distribution curves for the Oak Creek subarmour and armour sediment areshown in Fig 7.8. These values are from Parker (1990) and exclude the sediment less than 2.0mm in diameter. The curve for the subarmour sediment has been slightly modified at the upperend because Parker reports that 0% of the subarmour sediment sampled was in the range 102-203 mm, whereas 8% of the armour layer falls within this range. The upper end of thesubarmour curve was smoothed to approximate a log-normal distribution with the maximumgrain size equal to 203 mm. This modified curve now has 3% of the subarmour sediment in therange 102-203 mm.The composition of the armour layer for an imposed bed shear stress of 45 N/rn2was calculatedusing Eqn (4.44). The selected value of the bed shear stress gives a value = 0.06 andtherefore observed and calculated values ofD50 match. The results are plotted along with theobserved grain size distribution in Fig. 7.8 and there is an excellent match. However the Parker191surface-based bedload relation that is used to calculate the composition of the armour layer wasdeveloped from bedload transport data collected on Oak Creek, and therefore the good matchis not surprising.As an independent test of the Parker theory the experimental data of Dietrich et aL (1989) willbe used. Dietrich et at. performed three flume experiments to investigate the effect of sedimentsupply on the development of the armour layer. They initially started with very high sedimentfeed rates which at equilibrium resulted in no armour layer development (8o = 1). The sedimentsupply rate was progressively reduced and the bed surface was observed to coarsen to values of850 = 1.16, and 1.32. The 850 V*D5 solution curve was calculated using the sediment feedgrain size distribution, and is shown together with the three observed armour-layer values inFig. 7.9. There is close agreement with the theory for two of the points, however the data pointcorresponding to the highest sediment feed rate plots well below the theoretical curve.The Parker surface-based bedload transport relation was developed over a fairly limited rangeof data. In particular the reduced hiding function (Eqn. 4.36) was developed for r’j50 0.05.At this value the armour layer is only moderately mobile. However for higher sustained shearstresses where the bed surface is much more mobile and the value of öo approaches 1, it isunlikely that Eqn (4.36) would still be valid. The results from the analysis of the data fromDietrich et a!., which is admittedly based only on three data points, suggest that the Parkertheory is applicable only up to moderate values of z*D50, say up to a maximum of 0.07. Beyondthis value the mobility of the bed sediment and the hiding relations appear to be considerablydifferent from those observed at Oak Creek. Fortunately natural gravel rivers generally have192bankfull values of T*D50 less than 0.07, with most falling in the 0.03 - 0.06 range where theParker theory can be applied.7.4.2 Effect of Sediment Gradation onTheoretical ö-r’D50 solution curves were calculated for a range of sediment mixtures fromuniform to poorly sorted which are shown in Fig. 7.10. The sediments all have the same mediangrain size (d50 = 25 mm). The geometric standard deviation, 0g, was used as a measure of thedispersion about d50. The resulting 650 - V*D50 curves are shown in Fig. 7.11. The sedimentgradation is shown to have a large influence on the shape of the curves and the values of 850.The sediment size parameters from Figs 7.10 and 7.11 are shown in Table 7.1. The value ofD50 is approximately equal to d85 in this example.Table 7.1. Summary of sediment size parameters from Figs 7.10 and 7.11. The values of Ugand d90 are from the subarmour sediment, and 850 and D50 are from the computed armour layer.The value ofD50 is approximately equal to d85.(Jcg d85 D50(mm) (mm)0 25 1 250.28 45 1.6 400.43 75 3.1 780.50 100 4.4 1107.5 EFFECT OF SEDIMENT LOADThe effect of increasing sediment load on the channel will be modelled in this section and thetheoretical results compared to experimental and field observations, and qualitative relationsobtained from these observations.193The effect of a 10-fold increase in the imposed sediment load from 2.5 X 106 to 2.5 X i07kg/year on the optimal hydraulic geometry for the set of independent variables described inSection 7.3 is shown in Fig. 7.12(a) and (b). The values of W and S increase in value withincreasing sediment load by 32% and 56% respectively, while the values of 7, and ]3 show adecrease of 24% and 15%.These modelled adjustments above are in general agreement with observed adjustments. Recallthe qualitative proportionality of Lane (1955a) which can be rearranged as follows:QdSG (1.1)where Qd is the characteristic or dominant discharge which is here assumed to be equal to Qb.,and D is a characteristic sediment grain diameter, which is presumed here to be D50.Since the value of Qd is assumed to be constant (cf Section 7.6.1), Formula (1.1) indicates thatan increase in Gb will produce an increase in the value of S, and a decrease in D. The decreasein the value ofTho with increasing sediment load was discussed in Section 7.4 and is supportedby numerous experiments and field observations (eg Dietrich et at., 1989; Lisle and Madej,1989; Kuhnle, 1989).A similar qualitative relation proposed by Schumm (1969):194W2ScxG,, (1.3)where 2= the meander wavelength, and = sinuosity, which is equal to the valley slope dividedby the channel slope, S / S.Formula (1.3) was derived from field observations, and indicates that an increase in Gb willresult in an increase in W, S, and 2, and a decrease in Y and . By definition an increase in S isequivalent to a decrease in as the value of S is assumed to remain constant. The adjustmentof 2 is not considered in the optimization modelling.For an increase in the imposed sediment load the results obtained from the optimizationmodelling show good qualitative agreement with the proportionalities developed by Lane(195 5a) and Schumm (1969). These relations are widely used as guidelines for interpreting andpredicting river adjustments, and the agreement between the modeffing results, and the abovequalitative formulas lend general support to the optimization approach.7.5.1 Valley Slope ConstraintThe adjustments to increasing sediment loading modelled above do not consider the valleyslope constraint which represents a physical bound that defines the maximum channel slope thatcan be attained over engineering time scales. The valley slope is considered to have developedover geologic time, and therefore over engineering time scales can be considered to be anindependent variable.195The initial hypothetical channel from Section 7.5 with a slope of 0.00426 and a sedimenttransporting capacity equal to 2.5 X 106 kg/year, would require a value of 1.57 if it were toadjust to an imposed load of 2.5 X kg/year which requires an increase in the channel slopeof 57% to 0.00668 (Fig. 7.12(b)). The 62 rivers listed in Hey and Thome (1986) have anaverage sinuosity of only 1.34, and therefore the majority of these could not increase theirchannel slope by as much as 57%. When the required channel slope exceeds the valley slope,the valley slope constraint is violated and the solution is not feasible.The valley slope defines the maximum channel slope that can be attained, and thereforedetermines the maximum sediment load which can be accommodated by a river system and stillmaintain a stable equilibrium channel. For example if the valley slope for the river system beingmodelled in Section 7.5 were 0.0056, then from Fig. 7.12(b) it can be seen that the maximumimposed load that can be maintained is about 1.0 X i0 kg/year. For a sediment load in excessof this value the channel cannot develop the required slope and therefore sediment continuitycannot be maintained. An aggrading, unstable, and possibly braided channel would be expected.7.6 BANKFULL DISCHARGE AS A DEPENDENT VARIABLEUp to this point in this thesis, and generally throughout the literature, the value of Qbf isconsidered to be an independent variable which is a function of the catchment hydrology. Thiswas questioned in Section 3.3. The optimization model will now be used to show that there isan optimum value of Qbf which suggests that, in the case of a river with an active floodplain, the196bankilill discharge should more correctly be considered a dependent variable. This analysis doesnot apply to incised channels.The equilibrium geometry for the hypothetical channel described in Section 7.3 was calculatedfor a range of Qbf values from 55 to 155 m3/s. The equilibrium values of S for the various Qbfvalues are shown in Fig 7.13(a). It is evident that an optimal value of Qbf which corresponds toa minimum value of S (and hence maximum ii), occurs at about Qbf= 85 m3/s.The reason for the existence of this optimum value of Qbf can be realised from examination ofthe transport rates and total annual loads transported over each flow interval. Contrast thetransport rates for values of Qbf equal to 65 and 135 m3/s which are given in Table 7.2. Forflow interval Q which equals 59.8 m3/s the sediment transport rates are 0.84 and 0.26 kg/s foreach channel respectively. Over one year Q transports 540,000 kg or 21.6% of the total for avalue of Qbf = 65 m3/s, while for Qbf = 135 m3/s the total transported over one year by Q4 isonly 170, 000 kg, or 6.8% of the total. The same flow in the smaller channel has a higher rateof transport because the depth of flow and therefore the value of Thed are greater. For Qbf = 65m3/s the depth of flow and value of Thed for flow Q are respectively 1.18 m and 46.1 N/rn2,while for a value of Qbf = 135 m3/s the respective values are 0.85 m and 35.6 N/rn2.As the value of Qbf increases the sediment transport rate at bankfull also increases. This is aresult of the greater discharge as well as the wider channel bed across which bedload transportcan occur. Compare the bankfl.ill sediment transport rates in both channels. For Qbf = 65 rn3/sthe bankflill transport rate is 1.13 kg/s, while for Qbf 135 m3/s it is 6.26 kg/s. However as Qbf197Table7.1Comparisonofflowparametersfor Qbf=65m3/sandQbf135m3/s.00I 2 3 4 58.130.448.159.865.0Qbf=65m3/s=135m3/s1PiYGb’GbIPiYGb’Gbms/smkg/skg/ym3/smkg/skg/y0.7530.1410.0500.0200.0370.040.81 1.041.181.2300.0260.3950.8391.1317.01.2X6.2X5.4Xi051.3X106I 2 3 4 5 6 7 8 9 10 11 128.130.448.159.869.879.889.999.9109.9119.9129.91350.7530.1410.0500.0200.0120.0080.0040.0030.0030.0020.0010.0040.270.590.760.850.93 1.001.071.131.191.251.311.3400.0010.0470.2640.5570.8801.331.982.833.965.396.260.60490074,0001.7X2.1X2.1X1.8X105 2.1X2.2X2.OX2.2X105 8.1Xincreases the duration which the bankflill discharge is equalled or exceeded also decreases.From Qbf= 65 to 135 m3/s the duration of the bankfull flow decreases from 0.0365 to 0.0041.Therefore the optimal value of Qbf results from the opposing tendencies of high sedimenttransport rates for low flows in the smaller channels, and higher effective discharges and greaterbed widths which combine to produce larger transport rates in the larger channels.7.6.1 Effect of Sediment Load on Bankfull DischargeThe effect of sediment load on the channel geometry will be used to demonstrate evidence thatQbf should be considered a dependent variable. It was shown in Fig 7.13(a) that for an imposedbedload of Gb = 2.5 X 106 kg/year the model returns an optimal value of Qbf = 85 m3/s. Theresults for an imposed load of 2.5 X i07 kg/year are shown in Fig. 7.13(b) and indicate anoptimal value of about Qi,= 135 m3/s. Therefore the value of Qbf is seen to adjust together withthe other dependent variables to define an optimal solution.There is some support for this result. The qualitative formula of Lane (1955a; see above Eqn(1.1)) indicates that an increase in Gb will cause an increase in Qd (or Qbf) which is in agreementwith the model predictions. Williams (1978) suggests that the recurrence interval (and hencethe magnitude) of the bankfull discharge is possibly affected by the sediment load. It is aninteresting result that requires follow up work to determine its significance.This result can have considerable influence on the solution. For instance Fig. 7.14 shows theeffect on the channel width of a 10 fold increase in sediment load from Gb = 2.5 X 106 kg/year199using the independent variables from Section 7.3. The variation in Wassuming a constant valueof Qbf = 85 m3/s, as well as variable Qbf are shown. A much wider channel is predictedassuming a variable QbrIt can be seen from Figs 7.13(a) and (b) that the optimum value of S, and hence i, is relativelyinsensitive to changes in Qbf For instance it can be seen from Fig. 7.13(a) that a 55% increasein Qbf from 55 to 85 m3/s produces a decrease in S (or an increase in i) of only 3.4%. A widerange of values of Qbf are therefore “near optimal”. It would seem reasonable that the “drivingforce” for a channel to attain the optimum value of Qbf would be proportional to ÔS / ôQbf.Since the value of this partial derivative appears to be small, the “driving force” for the channelto attain the optimum value would be correspondingly small. Therefore while a theoreticaloptimal value of Qbf can be determined, in reality there is a wide range of values of Qbf whichcome close to satisfying the objective function, and therefore in natural systems it could beexpected that the observed bankfull flows would be randomly scattered over a wide rangeabout the theoretical optimum.7.7 APPLICATION OF THE MODELThe formulation of the model allows it to be applied to any situation where there is alteration inthe volume and size distribution of the imposed sediment load, the volume and timing of theflows, or the properties of the bank sediment. For example the construction of a dam will affectthe sediment supply and flows. The sediment supply will be reduced dramatically, ofteneffectively to zero directly downstream from the dam. The area of interest may be somedistance downstream and the reduction of the sediment load from the dam may represent only a200fraction of the total transported at this point. The flows are usually affected to a large extent,typically by truncating the higher flows, and increasing the proportion of low flows. The totalrunoff volumes may or may not be affected. The bank stability parameters may not be affected.Other potential applications relate to land-use changes. For instance removal of forest cover bylogging or for agricultural development may increase the sediment production from thecatchment, and yet may or may not affect the runoff to any great extent. Often the riparianvegetation is affected by these developments, and as has been demonstrated in Chapter 6, thiscan have a profound influence on the bank stability.Regardless of the type of development or catchment disturbance, the input required for themodel is a flow-duration curve, an estimate of the volume of sediment load and grain sizedistribution, and estimates of the bank stability parameters. Hydrologic modelling and sedimentbudget studies may be necessary together with field observations to determine the appropriatevalues to use as input to the model. These values may be difficult to measure in an intact systemlet alone to predict the values for a disturbed catchment. Estimates of the current sediment loadof a stable river system can be obtained using the observed hydraulic geometry, which includesthe grain size distribution of the bed surface, together with the measured or estimated flowduration curve.The model can be used to perform sensitivity analyses to determine the effect of a potentialdevelopment on the river channel where the post-development inputs cannot be accuratelydetermined. For instance some studies have indicated that typical logging practises in coastal201B.C. and the Pacific north-west can result in an 8 to 10 fold increase in the yield of coarse bed-material sediment as a result of increased mass-wasting processes (eg Madej, 1978). The modelcan be calibrated using the observed hydraulic geometry, and then the sensitivity of the channelto increased sediment load can be determined.When using this optimization model it must be realised that the optimum value is only atheoretical value which the river may show a tendency to adjust towards. The optimizationmodel is only an adjunct to other techniques such as air photo and field monitoring of observedchannel adjustments of the river of interest and nearby channels that may have been subjectedto similar disturbances. Any modelling results must be tempered with sound engineeringjudgement, and must recognise the location of geologic controls that may limit the computedadjustment.7.8 HOW AND WHY DO ALLUVIAL CHANNELS OPTIMIZE?Up to this point the issue of how and why do alluvial channels reach the optimal solution hasnot been addressed. The value of the modelling approach has relied on its empirical success inpredicting channel geometry. It seems that the general reluctance of the engineering communityto more fully accept this approach to river adjustments is largely due to the “lack of a physicalbasis” for this approach.Yang and Song (1979, 1986) have attempted to explain their minimum energy dissipation rateand minimum unit stream power hypotheses by arguing from fundamental fluid mechanicsprecepts. Yang and Song suggest that many aspects of the behaviour of fluids, and not just202river channel adjustments, can be explained by their minimum energy dissipation ratehypothesis. However their explanations have not met with widespread acceptance.The argument to be presented in this section will attempt to demonstrate that river channelsdisplay the tendency to move towards the optimal configuration as a result of randomperturbations which exist in abundance in natural river systems. It is the asymmetry in riverchannel response to the random perturbations (which will be explained below) that drives thesystem towards the optimum.Figure 7.15 shows a single solution curve that satisfies the continuity, bedload, and bankstability constraints for a particular set of independent variables. The optimal geometry islocated at Point A. Any channels that are located in the area above the solution curve are overcapacity with respect to transporting sediment, and have the capacity to transport sediment inexcess of the imposed load. Similarly channels located in the area below the curve do not havethe capacity to transport the imposed sediment load. The further a point is above or below thesolution curve, the greater the sediment transport capacity of the channel diverges from theimposed load.Consider a channel that is initially located at Point B. This channel is narrower and steeper thanthe optimum. Ideally this channel could remain indefinitely at Point B in apparent stability as itsatisfies all of the constraints. However in nature there are abundant opportunities for thechannel geometry, at least locally, to be perturbed from B. Consider a local widening whichcould result from an event such as tree falling into the river which acts as a locus for deposition,203from which a gravel bar is developed. This could deflect the current locally against the bankresulting in local bank erosion and widening of the channel to Point C. The channel is now, atleast locally, over capacity, and will attempt to fulfil its sediment transporting capacity throughbed and bank erosion.It is a primary assumption throughout this thesis that bed erosion and degradation is negligiblein comparison to bank erosion, and that bank erosion is concentrated along the outer banks ofmeanders. When the sediment transporting capacity of the channel becomes greater than isrequired to transport the volume of sediment imposed from upstream, net erosion will occuralong the reach as the channel seeks to fulfil its sediment transporting capacity. Since theerosion is presumed to be concentrated along the outer banks of the meanders, this erosion willresult in further development of the meanders, therefore causing a reduction in the channelslope. There may also be some limited channel widening. The net result is that the channeleventually returns to the solution curve at Point D. The result of the initial perturbation istherefore to drive the channel geometry closer to the optimum.Perturbations which displace the channel geometry below the solution curve will produce alocal channel geometry that has insufficient capacity to transport all of the imposed sedimentload. The channel will not have the potential to modif,’ the channel boundary through erosion,and although the channel can modify its dimensions through deposition, it is difficult toenvisage how significant changes would result.204It is therefore possible that an asymmetry exists whereby perturbations which displace thechannel geometry above the solution curve result in a modification of the channel geometrywhich moves it towards the optimum. Perturbations which displace the channel geometry belowthe solution curve into the zone of under capacity will be damped out and have lesser effect onthe geometry.Once the channel reaches the optimum, any changes in the channel width will reduce thesediment transport capacity and these perturbations will be damped out. Any steepening of thechannel that produces an excess of sediment transporting capacity (such as a meander cut-off)will result in erosion principally along the outer bends of the meanders, which will result in areduction in the channel slope returning the channel to the optimal configuration.The above explanation of the optimization process is admittedly incomplete. The role of thesecondary currents is not considered. This explanation does indicate however, that the tendencyfor alluvial rivers to attain an optimum configuration is not a systematic process, but resultsfrom random heterogeneities in the channel.7.9 COMPARISON WiTH OTHER NUMERICAL MODELSThe modelling strategy developed in this thesis is fundamentally different from most othernumerical models which have been developed to predict river adjustments. The most widelyused modelling approaches are the dynamic models which were mentioned in Section 1.4.4,2.and summarised in Hey (1988). The HEC-6 model HEC, 1974; Thomas and Prashun, 1977) isthe most well known. These models assume that the channel width remains constant and205changes in the bed elevation are determined for each channel interval at each time step.Changes in the channel slope are assumed to occur as a result of bed aggradation ordegradation. Osman and Thorne (1988) and Thorne and Osman (1988) modified this approachto include a bank erosion algorithm for cohesive channels. Chang (1982) proposes that cross-sectional changes, which include a reduction or increase in the channel width, can bedetermined for each time step by using the minimum stream power hypothesis. These modelscan be run until steady state conditions have been attained.In contrast the optimization model proposed herein does not consider the transient, dynamicchannel adjustments and gives only the final steady-state channel geometry. Changes in channelwidth are modelled, and the slope adjustments over engineering time scales are assumed tooccur due to changes in the channel sinuosity, and not through aggradation or degradation.The results obtained by the different modelling approaches are illustrated in Fig. 7.16. Anydifferences that would arise between the models due to the different flow resistance andbedload transport equations used in each model are ignored. Two solution curves similar to Fig7.15 are shown. The two curves in this case can be taken to correspond to two values ofimposed bedload, the upper curve representing the higher value of Gb. Consider a stablechannel located at Point A on the upper curve which is subject to a reduction in sediment loadto the value represented by the lower curve. The optimization model predicts that the newstable equilibrium geometry will be given by Point B on the lower curve which represents areduction in the channel slope and a decrease in the width. There is no indication of thepathway or the time required for the channel to adjust from Point A to B.206A fixed width model such as HEC-6 will return a progressive decrease in the channel slope dueto bed degradation as the sediment transporting capacity of the channel will initially be greaterthan the sediment supply following the reduction. Assuming that the degradation is notterminated by the development of a static armour layer (which is almost certain to occur in agravel river), the HEC-6 type models will eventually arrive at a steady state condition close toPoint C on the lower curve. Since the degradation has not involved any changes in channelwidth the final channel must be deeper than the original and it is probable that the bank stabilityconstraint will not be satisfied at steady state. However with the exception of the OsmanThorne model to be discussed below, these dynamic models do not consider the stability of thebanks.The Osman-Thorne model recognises that the bank stability and erosion must be addressed. Intheir model they calculate the net bank erosion in each time step and add this volume of erodedbank sediment into the sediment continuity calculation for each time step. The eroded banksediment partially fulfils the sediment transporting capacity of the particular channel interval,and therefore reduces the net bed degradation. The Osman-Thorne model will reach a steadystate solution at some location between Points C and D on the lower curve. The exact locationwill depend upon the relative erodibility of the bed and bank sediment, for highly resistantbanks the final solution will be close to Point C (slope adjustment only), and for easily erodiblebanks relative to the bed the final solution will be closer to Point D (width adjustment only).207The numerical model of Chang (1982) uses sediment routing to compute the change in thecross-sectional area of the channel, AA, which can be positive or negative depending uponwhether the channel is aggrading or degrading. The minimum stream power (Section 2.2.3) isused to assign i4 to the bed, or the banks, or both. The channel width is first varied until S,and therefore i.QS is minimized. Once the banks are adjusted, the remaining i4 is applied to thebed. One feature of the model is that uniform flow is not assumed, and therefore the slope ofthe bed and the energy gradient are not necessarily the same. The dynamic optimization modelof Chang predicts that the channel in Fig 7.15 will move along a direct line from Point A toPoint B.Despite the capacity for width adjustment, the numerical model of Chang does not containalgorithms for calculating the distribution of the boundary shear stress, or for bank erosion ordeposition. For example the model formulation does not recognise that the rate of bank erosionis dependent upon bank shear stress and the erodibility of the bank sediment. Furthermore it isassumed that the channel adjustment is a continuous process of dynamic equilibrium, and thatthe channel maintains equilibrium at each time step. However in the static optimization modelproposed in this thesis, while the channel is assumed to ultimately reach an optimum, thetransient adjustments are not presumed to adhere to any principle of equilibrium.Clearly the models appear to give very different results, particularly the optimization model andthe Osman-Thorne model which both consider the bank stability and therefore one wouldassume that they would be in closer agreement than the fixed-width, HEC-6 type models.Importantly the Osman-Thorne model does not allow for bank deposition and therefore channel208narrowing in this case is not possible. The optimization model developed in this thesis does notexplicitly include bank deposition in its formulation, however bank deposition is implied whenthe value of falls below the critical value for bank stability since the banks of an optimalchannel develop at the limiting stability.One other important difference between the optimization model and the dynamic models is thelatter’s dependence upon the initial conditions. It is evident from Fig. 7.16 that the finalsolution for the dynamic models is dependent upon the initial conditions, different initialconditions will yield different final solutions. However the optimization model is independent ofthe initial conditions. It was argued in Section 7.8 that random perturbations in the course of ariver channel’s adjustment were necessary for it to attain the optimal geometry. These randomevents are not accounted for in the dynamic models and are important as they have the effect ofresetting to some degree the initial conditions. While final steady-state solution from thedynamic optimization model of Chang (1982) is not dependent on the initial conditions, thepath and the time taken during adjustment are.An interesting possibility would be to include a random component into the non optimizationdynamic models to see if there is a tendency for the channel to adjust to the optimalconfiguration.2097.10 AN EXAMPLE OF RiVER BEI{AVIOIJR WITHIN AN OPTIMIZATIONFRAMEWORKThe optimization framework can be used to interpret the response of an alluvial river followingriver training activities. A solution curve calculated using the optimization model is shownschematically in Fig 7.17. The valley slope constraint is also indicated. The feasible part of thesolution curve, that is the section of the curve that satisfies all of the constraints including thevalley slope constraint, is located below the valley slope constraint where S S.Assume that the channel is originally at the optimum (Point A). In order to reduce flooding,improve navigation, and to “stabilise” shifting channels, a meandering reach is commonlystraightened. If the straightened channel is oriented parallel to the valley axis, it will have aslope equal to S. The channellized reach can be designed such that it possesses approximatelythe same sediment and discharge capacity as the natural meandering channel, and also hasstable banks. A straight channel which is oriented parallel to the valley axis and also is designedso that is satisfies the discharge, bedload, bank stability and valley slope constraints is indicatedat Point B. Note that this represents the narrower of two possible options; a wider channel thatalso satisfies all of the constraints is located at Point C.However, although all of the constraints are satisfied at Points B and C, by straightening thereach the channel slope is increased, and therefore the channel is no longer at the minimumslope that corresponds to the optimal geometry. These straightened channelled reaches mustusually be maintained as any slight channel curvature or irregularity is often amplified by theflow, and this can result in the development of meanders which cause a reduction in the channel210slope. The tendency of a straightened channellized reach to develop meanders is interpreted asthe channel attempting to reestablish the optimum geometry.7.11 SUMMARYThe complete model formulation was presented in this chapter. The model was formulated touse the fill range of flows, which are represented by a flow-duration curve, as inputs to themodel. Furthermore the composition of the bed surface is able to adjust to the optimalconfiguration which is based on the maximization of the coefficient of sediment transportefficiency, i. This represents a significant advance over previous optimization formulationssuch as those of Chang (1980), White et al. (1982), and Millar and Quick (1993b).Rigorous verification of the model was not attempted, however some well known features ofgravel rivers were demonstrated. The effect of sediment load on the equilibrium channel wasshown to agree well with field observation. The model predicts that an increase in the sedimentload will result in the channel becoming wider, shallower, steeper, with a decrease in the valueofD50 of the bed surface.The common assumption of the bankfull discharge as a fixed independent variable wasexamined. It was demonstrated that there is a value of Qbf for rivers with active floodplains thatcorresponds to an optimum. However the value of i is relatively insensitive to different Qbfvalues, and it is probable that natural rivers show a range of Qbf values that are scatteredrandomly about this optimal value.211It is argued that natural channels tend towards an optimal geometry as a result of randomperturbations which exist in abundance in natural rivers. These random events are notconsidered in the dynamic, numerical models of river channel adjustment. Furthermore thedynamic models all attribute changes in the channel slope to either aggradation or degradation,and do not consider changes in the channel sinuosity.2122001501100500Fraction of Time Discharge is Equalled or ExceededFigure 7.1. Flow-duration curve showing numerical approximation. The flow-duration curve inthis example is subdivided into 18 discharge intervals. A maximum discharge value must beselected that is equalled or exceeded all of the time which is the maximum discharge on record,or an approximation.0.001 0.01 0.1 1213a)CLL100806040200Grain Size (mm)200Figure 7.2. Sediment gradation curve showing numerical approximation. The sedimentgradation curve is subdivided in this example into 7 sediment intervals.2 5 10 20 50 100214Figure 7.3. Flow-chart for fhll optimization model formulation.N215Initialize TrialD and a ValuessgCalculate TransportRates For AllGrainsizesCalculate all F That Use CalculatedSatisify Equal Mobility D and a ValuesEqn(4.52) sg •Calculate NewD and a ValuessgHavesg and a, ValuesConverged?CalculateGb D50End BedloadSub FunctionFigure 7.4. Flow-chart for bedload subfunction. An iterative scheme is required to calculate thevalues ofD and u.216(a)2a,t0Cl,I4-a,a,Co0100(b)1 E + 071E+061E+05- 1E+041E+031E+021E+01I-1 E + 00Figure 7.5. Calculated sediment transport rates as a function of discharge. (a) The values of Gb’which have units of kg/s. (b) The values of Gb which have units of kg/y.0 20 40 60 80Discharge (m3Is)Discharge (m3/s)2173025ENU)a). 20a)CVa,Ca)Va,15100 100Figure 7.6. Variation Of d5Oload with discharge. The value of the subarmour d50 is 25 mm.20 40 60 80Discharge (m3Is)2182.62.42.220‘01.81.61.41.210.01Figure 7.7. The calculated 850-solution curve for Oak Creek, Oregon. Data fromMilhous (1972). The solution curve was calculated using the Parker (1990) surface-basedsediment transport relation.0.02 0.03 0.05 0.1 0.2 0.3T*D50219100ICU806040200Grain Size (mm)Figure 7.8. Armour and subarmour grain size distributions for Oak Creek, Oregon. Data fromMilhous (1972). The calculated values were determined using the Parker (1990) surface-basedsediment transport relation.2 5 10 20 50 100 2002201.51.41.30U,‘01.21.110.01Figure 7.9. Calculated öo-TD50 solution curve for data from Dietrich et a!. (1989). Thesolution curve was calculated using the Parker (1990) surface-based sediment transportrelation.0.02 0.03 0.05 0.1 0.2 0.3T*D502211008060a)CU402002Figure 7.10. Sediment gradation curves for a range of crg values.5 10 20 50 100Grain Size (mm)200222540210.01T*D50Figure 7.11. Calculated c5o - TD50 curve for a range of ag values. The curves were calculatedusing the Parker (1990) surface-based sediment transport relation.0.02 0.03 0.05 0.1 0.22230.0070.00650.0060)0.00550.0050.00450.004 3.OE+06 5.OE+06 1.OE+07Sediment Load (kg/y)1.20.90.80.080.0750.070.065 ,,0.060.0550.056055__50EV 45403530WidthDepth3.OE+06 5.OE+06 1.OE+07 2.OE+07Sediment Load (kgly)2.OE+07Figure 7.12. Effect of sediment load on channel geometry. (a) W, Y. (b) S, D50. The valueswere calculated using the optimization model.2240.00480.0046 -a)0.0•0.0042 —0.00460 80 100 120 140Bankfull Discharge (m3/s)0.00720.007 -::::0.006460 80 100 120 140 160Bankfull Discharge (m3/s)Figure 7.13. Calculated variation of S with Qbf (a) Gb = 2.5 X 106 kg/y. (b) Gb = 2.5 Xkg/y. The values were calculated using the optimization model.22570-. Variabre bfConstant bf602.•130 I I3.OE+06 5.OE+06 1.OE+07 2.OE+07Sediment Load (kgly)Figure 7.14. The calculated variation of W with Gb for constant and variable Qbf The valueswere calculated using the optimization model.226Over Capacity0C’)WidthFigure 7.15. Solution curve used to demonstrate channel optimization. The solution curve isshown schematically. Point A is at the optimum, Points B and D satisfy the constraints, but arenon-optimal, and Point C is overcapacity with respect to transporting sediment.CSolutionCurveUnder Capacity2270Cl)WidthFigure 7.16. Solution curves used for a comparison of numerical models. The solution curvesare schematic. The solution curves represent two different imposed sediment loads, with theupper curve corresponding to the higher imposed load. Points A and B are at the optima oftheir respective solution curves.CB228Non-Feasible\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\a /B\ /C.2 s=S \ Feasible /CD VAWidthFigure 7.17. A solution curve used to interpret the behaviour of a meandering river that hasbeen straightened. The valley constraint is also indicated with shading on the non-feasible side.The initial equilibrium geometry is given by Point A. This solution is opimal. The straightenedchannel will be located at Points B or C, both ofwhich are non-optimal.229CHAPTER 8CONCLUSIONS AND RECOMMENDATIONS8.1 SUMMARYIn this thesis an optimization model has been developed to predict the stable, equilibriumgeometry of alluvial gravel-bed rivers for a given set of independent variables. The independentvariables are the discharges, both the magnitude and duration which are represented by a flow-duration curve; the mean annual bedload, both total mass and size distribution, which isimposed on to the channel reach from upstream; and the geotechnical properties of the banksediment.The unknown dependent or decision variables to be solved for include the channel width, depth,bank angle, roughness, and grain size distribution of the bed surface. The dependent variablesadjust subject to the constraints of discharge, bedload and bank stability to determine a channelgeometry which is optimal as defined by , which is the coefficient of sediment transportefficiency (Eqn 2.2).The work in this thesis is an extension of work by Chang (1979, 1980) and White et a!. (1982)whose models have predicted the geometry of sand and gravel rivers with reasonable success,230however the degree of scatter associated with these models limited their application toquantitative engineering applications. Others such as Yang (1971, 1976) and Yang et al. (1981)have proposed similar models which describe various features of alluvial rivers. The advances inthis thesis over the earlier optimization models are presented below:8.1.1 Inclusion Bank Stability Analysis.The bank stability analysis was introduced in Millar (1991) for channels with noncohesive banksediment that was unaffected by bank vegetation. It has subsequently been extended to includethe influence of bank vegetation and also for channels with cohesive banks.For noncohesive bank sediment the model (in a simplified form) was tested on the publisheddata of Andrews (1984) and Hey and Thome (1986). It was shown that the inclusion of thebank stability procedure significantly improves the model performance. The results of theinfluence of the bank vegetation on the channel stability are in good agreement with the resultsthat Andrews and Hey and Thorne obtained from empirical regime analyses. The influence ofthe bank vegetation is to stabilise the bank sediment, allowing the banks to withstand highershear stresses. This results in channels that are narrower, deeper, and less steep than theirunvegetated counterparts. The change in W is the largest, followed by Y, and the smallestchange is in S.The analysis of channels with cohesive banks indicates that channels tend to be either bankheight or bank shear constrained, and less commonly both bank stability constraints are active.In general steeper channels with high sediment loads and low bankfull discharges tend to be231bank shear constrained, while conversely channels with flatter gradients, low sediment loads,and large bankfhll discharges tend to be bank height constrained. The terms “low” and “large”in the previous sentence are relative, and depend largely on the properties of the bank sediment.In a typical channel with only one active constraint, the other constraint is degenerate in that isdoes not influence the channel geometry. It is therefore necessary to determine which of thetwo possible constraints is active in order to assess the stability of the channel. Failure to do socan lead to an assessment of the channel stability based on the degenerate constraint which willlead to erroneous results with regard to the stability of the banks.An analysis of published data from Charlton et al. (1978) for channels with cohesive banksindicates that the bank vegetation has a strong effect on the value of Tcrir (Fig 6.8). Channelswith treed banks have higher values of rcrit than channels with grassed banks. This resultsindicates that the effect of bank vegetation on channels with cohesive banks is similar to thenoncohesive case and that rivers with heavily vegetated banks would be narrower, deeper, andless steep than their unvegetated counterparts.These results regarding the bank vegetation have important consequences for streammanagement as the removal of the bank vegetation can reduce çS or zrjt and destabilise thechannel. Published observations have shown that streamside logging can result in increasedwidth and destabilisation of alluvial channels (Roberts and Church, 1986; Hartman andScrivener, 1990). The optimization modelling could be used to determine the sensitivity of achannel to streamside logging.2328.1.2. Modelling Using The Full Flow-Duration DataPrevious optimization models such as Chang (1979, 1980) and White et aL (1982), as well asthe simplified optimization models presented in Chapter 6 calculated the sediment transportingcapacity of the channel based only on the bankfull rate. It was argued in Chapter 4 that this isinadequate as flows less than bankfull can transport significant volumes of sediment.Furthermore not only is the bankfull sediment transport rate important, but the fraction of thetime where Qbf is equalled or exceeded determines the total volume of sediment that istransported by the channel over one year. The final version of the optimization model presentedin Chapter 7 uses the complete range of flows which are represented by a flow-duration curve.This curve is discretised and used as input to the model. Over-bank flows are also consideredand as a first approximation an infinitely wide flood plain is assumed. The depth of flow withinthe channel, and more importantly, the values of vE,d and Thank reach their maximum value atQbj; and are presumed to remain constant at this value for all larger discharges.8.1.3 Adjustment of the Bed Surface.The grain size distribution of the bed surface or armour layer, which is represented most simplyby D50, is considered to be a dependent variable. The size distribution of the bed surface and thebedload transporting capacity of the channel are calculated using the Parker (1990) surface-based bedload transport relation which is modified in Chapter 4 to apply to variable flows.In this thesis the concept of “equal mobility” is presumed to apply over some significantduration which is generally assumed to be one water year. Therefore the “equal mobility”233concept as applied herein is that the bed surface or armour layer forms so as to render all sizefractions equally mobile over a period of one year. Therefore the grain size distribution of themean annual load transported by a channel in equilibrium is exactly the same as the distributionof the imposed load. It is assumed that the grain size distribution of the mean annual imposedbedload is identical to the that of the subarmour sediment (Parker, 1990).The definition of equal mobility developed in this thesis is significantly different from theoriginal definition by Parker et a!. (1982a) who assumed that equal mobility was achieved foreach flow. This is known to be incorrect from field observations (eg Fig. 4.11), and the originaldefinition of equal mobility was acknowledged by Parker et a!. as a first-order approximationonly.8.1.4 Bankfull Discharge as a Dependent VariableIt was shown in Chapter 7 that there is an optimal value of Qbf for a given set of independentvariables. This result indicates that the value of Qbf can be considered a dependent variable,rather than an imposed value which is commonly presumed.However the limited modelling undertaken in Chapter 7 indicates that the value of i isrelatively insensitive to Qbf, and from this it would be expected that the observed values fromnatural rivers would be randomly scattered across a wide range about the optimal value. This isan interesting result which deserves fhrther attention to determine its significance.2348.1.5 Verification of the Optimization ModelModelling results of Chang (1979, 1980) and White et at. (1982) showed reasonably goodagreement between modelled and observed values of W and 7, although the degree of scatterwas too great for quantitative engineering projects. This scatter has been reduced in Millar(1991) and Millar and Quick (1993a, b) and in Chapter 6 through the inclusion of the bankstability constraint. In addition Yang et a!. (1981) have shown that their optimization modelyields the exponents ofwell-known empirical regime equations.Furthermore good quantitative agreement was obtained between the modelling in Chapter 6and Millar and Quick (1993a, b) for the effect of bank vegetation on channels with noncohesivechannel banks, when compared to the empirical regime equations of Andrews (1984) and Heyand Thorne (1986).In addition general qualitative agreement between the optimization modelling and the formulasof Lane (1955a) and Schumm (1969) has been demonstrated by Chang (1980), Millar (1991),and in Chapter 7.These results represent only partial verification of the optimization theory that has beenpresented in this thesis. This theory was advanced beyond that which could be fully verifiedwithin the scope of this thesis. A program for full verification of the optimization model isoutlined in Section 8.2.6.2358.2 RECOMMENDATIONS FOR FUTURE WORKThe optimization model presented herein is formulated using the most suitable equations thatwere available from the literature. Some modifications were made, however no experimental orfield work was carried out during this study to develop new, or improve existing relations.Several aspects of the optimization model that can benefit from additional research arepresented below.8.2.1 Channel Form or Bar RoughnessIn Section 3.2 it was argued that the grain roughness could be calculated using Eqn (3.8),however no existing methods were available to calculate the form or bar roughness component.The similar methods of Einstein and Barbarossa (1952) and Parker and Peterson (1980) wereshown to be a result of spurious correlation and therefore do not yield meaningful estimates ofthe form roughness. The approach adopted by Prestegaard (1983) is promising although nopredictive equations have yet been developed.The nature of the form roughness can be studied using an approach similar to that adopted byPrestegaard (1983). The total channel roughness,J can be determined from channel surveys ata location where good quality flow measurements have been obtained, preferably at ahydrometric station. The grain roughness,f’, can be calculated using Eqn (3.8) and the value off”is equal to!- f’ from Eqn (3.7). Data obtained from detailed geomorphic mapping of thechannel, such as bar amplitude and spacing, vegetation etc, can be then used in an effort todevelop empirical relations forf”.2368.2.2 Surface-Based Transport Relation Based on t’bdThe Parker (1990) surface-based transport relation uses the total bed shear stress Tbed,, and notthe grain shear stress r ‘L,ed,, to calculate the transport rate. The use of the total bed shear stressand not the grain component is a result of Parker and Peterson (1980) who concluded that theform component of roughness and shear stress in gravel rivers is negligible, particularly for highin-bank flows. This was shown to be incorrect in Section 3.2, and furthermore it was shown inSection 4.2.4 that as little as 11% of the total shear stress can be assigned to the graincomponent (See Fig 4.6).Because Thed appears in Parker’s dimensionless bedload transport parameter W (Eqn 4.28),conversion of the transport relation to one based on r is not simply a matter of changing acoefficient, but rather all of the relations which include the hiding and straining functions needto be recalculated.8.2.3 Bank StabilityThe results in Chapter 6 indicated that for both noncohesive and cohesive banks the bankvegetation exerts a strong influence on the bank stability, and hence the channel geometry. Theeffect of the bank vegetation on the bank stability is reduced to a single parameter which is çfor noncohesive bank sediment, or rent for cohesive bank sediment. In both cases the effect ofthe bank vegetation was interpreted as an increase in and znjt which resulted from thebinding and reinforcement of the sediment by the root masses. An alternative explanation is thephysical shielding of the bank sediment by the vegetation which reduces the magnitude ofacting on the bank sediment. Field investigation is necessary to determine the actual role of the237bank vegetation and how it relates to bank stability, and whether the effect is an increase in gor rrjt, or a decrease in Tbank. Furthermore it is necessary to develop field-based methods fordetermining the insitu values of g and rit8.2.4 Secondary CurrentsThe “strong” secondary currents that occur as a result of channel curvature have not beenconsidered in this thesis. The “weak” secondary currents which are present in straight channelsare implicitly accounted for in the empirical boundary shear stress equations of Flintham andCaning (1988) which were developed from experimental data.One important effect of secondary currents is the redistribution of the boundary shear stress. Inthis thesis the bank and bed shear stresses are assumed to be uniform across their respectivechannel perimeters. However redistribution of the shear stress can have a large effect forexample on the sediment transporting capacity of the channel.Consider a channel that at bankfull has a mean bed shear stress value that is just below thecritical value for bed mobilisation. This channel will, based upon the mean value, transport anegligible volume of sediment. However if secondary currents redistribute the shear stress suchthat one half of the channel bed has a shear stress value about 1.5 times the critical value, andthe other half is only about 0.5 times the critical value, significant sediment transport will occuracross one half of the channel.238From the optimization model it is assumed that a reduction in the channel slope implies adecrease in the sediment transporting capacity of the channel. However this reduction inchannel slope is presumed to be accompanied by increased sinuosity, resulting in strongersecondary currents which could conceivably increase the sediment transporting capacity of thechannel.The optimization model calculates the mean value of Vbank, and this value must be less than theshear stress required for fluvial erosion of the bank sediment in order to satisify the bank-stability constraint. Locally however, the value of Tbank may exceed the value required for bankerosion, particularly in areas such as along the outer bank of meanders. Along the convex innerbank the shear stress is lower and deposition of sediment can occur forming a point bar deposit.Therefore even though bank erosion and lateral channel migration is occurring, the channel isstill considered to be in equilibrium, and “statistically stable”, if there is no net change in themean channel geometry over a representative channel reach.8.2.5 Formulation for Sand-Bed RiversThe equations used in this thesis have been selected and developed for application to gravel-bedrivers, however the model can easily be reformulated for application to sand-bed rivers. Forexample the Einstein transport relation could be used in place of the modified Parker (1990)relation.2398.2.6 Model VerificationBefore the optimization model can be used for quantitative engineering studies to predict theimpact of engineering structures or land-use practices on the channel geometry a fhll programof model verification needs to be undertaken to determine the reliability in the model.Initially data should be collected from stable rivers which appear to be in approximateequilibrium. Existing data sets such as Kellerhals et a!. (1972) were found to lack certain keydata. However these data sets could be updated with some field work to determine the valuesof the bank stability parameters, and review of the original records such as flow data, sedimentanalyses etc.The values of the independent variables obtained from field surveys can be input into the modeland a comparison made between the modelled and observed value of the dependent variables.This should include an analysis of the bankfull discharge to determine if the optimal value of Qbffound from modelling has any physical significance.If the model was found to predict the geometry of existing rivers, this would then improve theconfidence of using these models for predicting the response of a river to catchmentdisturbance.Despite the limited quantitative verification that has been done, the model in its present form isvaluable as it gives the user some insight into, and understanding of the nature of river channel240adjustments which can be used to improve the judgement of a scientist or engineer interested inthe response of alluvial rivers.241REFERENCESAckers, P., and White, W.R., 1973: “Sediment transport: new approach and analysis.” I Hydr.Div., ASCE, 99 (11), 2041 -2060.Alam, A.M.Z., and Kennedy, J.F., 1969: “Friction factors for flow in sand bed channels.” IHydr. Div., ASCE, 95 (6), 1973- 1992.Andrews, E.D., 1984: “Bed material entrainment and the hydraulic geometry of gravel-bedrivers in Colorado.” Geol. Soc. Am. Bull., 95, 371-378.Andrews, E.D., and Erman D.C., 1986: “Persistence in the size distribution of surflcial bedmaterial during an extreme snowmelt flood.” Wat. Res. Res., 22 (4), 191-197.Andrews, E.D., and Parker, G., 1987: “Formation of a coarse surface layer as the response togravel mobility.” In Sediment Transport in Gravel Rivers, Thorne, C.R., Bathurst, J.C., and,Hey, RD. (Eds.), John Wiley and Sons., Chichester, England, 269-300.Arulanandan, K., Gillogley, E., and Tully, R., 1980: “Development of a quantitative method topredict critical shear stress and rate of erosion of natural undisturbed cohesive soils.”, Tech.Rept. GL-80-5, US Army Eng. Wat. Exp. Stn., Vicksburg, Miss., 99 p.ASCE Task Commitee on Erosion of Cohesive Materials (Masch, F.D., Chair), 1967: “Erosionof cohesive sediments.” I Hydr. Div., ASCE, 94 (4), 1017-1049.242Bagnold, R.A., 1956: “Flow of cohesionless grains in fluids.” Roy. Soc. Lonti Phil. Trans.,249A (12), 23 5-279.Bagnold, R.A., 1966: “Approach to the sediment transport problem from general physics.”Prof Pap.,U.S. Geol. Surv., 422-I, 37 p.Barishnikov, N.B., 1967: “Sediment transport in channels with flood plains.” Pub., lASH, 75,404-412.Bazin, H.E., 1865: “Recherches hydrauliques, memoires presentes par divers savant.s science,mathematiques etphysiques.” Series 2, 19, Paris, France.Beschta, R.L., 1978: “Long term patterns of sediment production following road constructionand logging in the Oregon coast range.” Wat. Res. Res., 14(9), 1011-1016.Bettess, R., and White, W.R., 1987: “Extremal Hypotheses applied to river regime” InSediment Transport in Gravel Rivers, Thorne, C.R., Bathurst, J.C., and, Hey, R.D. (Eds.),John Wiley and Sons., Chichester, England, 15-45.Bhowmik, N.G., and Demisse, M., 1982: “Carrying capacity of flood plains.” .1 Hydr. Div.,ASCE, 108 (3), 443 - 452.Bishop, A.W., 1955: “The use of the slip circle in the stability analysis of slopes.” Geotech., 5(1), 7-17.Blench, T., 1957: Regime behaviour ofrivers. Butterworths, London, England, 138 p.243Blench T., 1969: Mobile bed fluviology, 2nd Ed., University of Alberta Press, Edmonton,Alberta, 164 p..Bray, D.I., 1979: “Estimating average velocity in gravel-bed rivers.” J Hydr. Div., ASCE, 105(9), 1103-1122.Bray, D.I., 1982a: “Flow resistance in gravel-bed rivers.” In Gravel-Bed Rivers, Hey, R.D.,Bathurst, J.C., and Thorne, C.R., (Eds.). John Wiley and Sons., Chichester, England, 109-133.Bray, D.I., 1982b: “Regime equations for gravel-bed rivers.” In Gravel-Bed River Hey, R.D.,Bathurst, J.C., and Thorne, C.R., (Eds.). John Wiley and Sons., Chichester, England, 517-552.Chang, H.H., 1979: “Minimum stream power and river channel patterns.” .1 Hydrol.(Amsterdam), 41, 303-327.Chang, H.H., 1980: “Geometry of gravel streams”. J Hydr. Div., ASCE, 106 (9), 1443-1456.Chang, H.H., 1982: “Mathematical model for erodible channels.” J Hydr. Div., ASCE, 108(5), 678-689.Chariton, F.G., Brown, P.M., and Benson, R.W., 1978: “The hydraulic geometry of somegravel rivers in Britian.” Rept. IT 180, Hydr. Res. Stn., Wallingford, England, 48 p.Chen, W.F., 1975: LimitAnalysis and Soil Plasticny. Elsevier, New York, 638 p.Chen, W.F., and Lui, X.L., 1990: Limit Analysis in Soil Mechanics. Elsevier, New York, 477p.244Chorley, R.J., 1978: “Basis for theory in geomorphology” In Geomorphology, PresentProblems and Future Prospects, Brunsden, D., and Jones, D.K.C. (Eds), Oxford Univ. Press,Oxford, England, 1-13.Chorley, R.J., and Kennedy, B.A., 1971: Physical Geography: A Systems Approach. Prentice-Hall International Inc., London, England, 370 p.Chow, V.T., 1959: Open Channel Hydraulics. McGraw-Hill, New York, 586 p.Chowdhury, R.N., 1978: Slope Analysis. Elsevier, New York, 423 p.Church, M., Wolcott, J.F., and Fletcher, W.K., 1991: “A test of equal mobility in fluvialsediment transport: Behaviour of the sand fraction.” Wat. Res. Res., 27 (11), 2941-295 1.Cruff R.W., 1965: “Cross-channel transfer of linear momentum in smooth rectangularchannels.” Water-supplypaper, 1592-B, U.S. Geol. Surv., B1-B26.Davies, T.R.H., 1980: “Bedform spacing and flow resistance.” .J Hydr. Div., ASCE, 106 (3),423-433.Davies, T.R.H., 1987: “Channel boundary shape-evolution and equlibrium.” In River channels,K. Richards (Ed.), Basil Blackwell, New York, 228-248.Davies, T.R.H., and Sutherland, A.J., 1983: “Extremal hypotheses for river behavior.” WaterRes. Res., 19 (1), 141-148.245Dietrich, W.E., Kirchner, J.W., Ikeda, H., and Iseya, F., 1989: “Sediment supply and thedevelopment of the coarse surface layer in gravel-bedded rivers.” Nature (London), 340,(6230), 215-217.Dury, G.H., 1973: “Magnitude - frequency analysis and channel morphometry.” In Fluvialgeomorphology, Morisawa. M., (Ed.), State Univiversity of New York, Binghamton, NewYork., 91-121.Egia.zaroff I.V., 1965: “Calculation of nonuniforn sediment concentrations.” J Hydr. Div.,ASCE, 91(4), 225-247.Einstein, H.A., 1942: “Formulas for the transportation of bed load.” Trans., ASCE, 107, 561-597.Einstein, H.A., 1950: “The bedload fhnction for sediment transport in open channel flows.”Tech. Bull., 1026, U.S. Dept. Agr., Washington D.C., 71 p.Einstein, H.A., and Barbarossa, N., 1952: “River channel roughness.” Trans. ASCE, 117, Pap.No. 2528, 1121-1146.Einstein, H.A., and Li, H., 1958: “Secondary currents in straight channels.” Trans., Am. Geo.Union, 39 (6), 1085-1088.Engulund, F., 1966: “Hydraulic resistance of alluvial streams.” I Hydr. Div., ASCE, 92 (2),315-327.Engulund, F., 1967: closure to “Hydraulic resistance of alluvial streams.” I Hydr. Div., ASCE,93 (4), 287-296.246Felenius, W., 1936: “Calculation of stability of earth dams” Trans. 2nd In!. Cong. Large Dams,4, 445.Flintham, T.P., and Caning, P.A., 1988: “The prediction of mean bed and wall boundary shearin uniform and compositely roughened channels.” In International Conference on RiverRegime, W.P. White (Ed.). John Wiley and Sons, Chichester, England, 267-287.Gilbert, G.K., 1914: “The transportation of debris by running water.” Prof Pap., U.S. Geol.Surv., 86.Gosh, S.N., and Roy, H., 1970: “Boundary shear distribution in open channel flow.” J Hydr.Div., ASCE, 96 (4), 967-994.Gosh, S.N., and Roy, H., 1972: “Boundary shear distribution in channels with varying wallroughness.” Proc. Inst. Civil Eng., London, England, 53 (2), 529-544.Graf W.H., 1971: Hydraulics ofSediment Transport. McGraw-Hill, New York, 512 p.Gressinger, E.H., 1982: “Bank erosion of cohesive materials.” In Gravel-Bed Rivers; Hey,R.D., Bathurst, J.C., and Thorne, C.R., (Eds.). John Wiley and Sons., Chichester, England,273-287.Griffith, M.W., 1927: “A theory of silt and scour.” Proc., Inst. Civil Eng., London, England,223, 243-3 14.Gniffiths, G.A., 1984: “Extremal Hypotheses for river regime: an illusion of progress.” Wa!.Res. Res. 20(1), 113-118.247Hammer, T.R., 1972: “Stream channel enlargement due to urbanization.” Wat. Res. Res. 8,1530-1540.Hartman, G.F., and Scrivener, J.C., 1990: “Impacts of forest practices on a coastal streamecosystem, Carnation Creek, British Columbia.” Can. Bull. Fish. Aquat. Sd., No. 223, 148 p.HEC, 1974: “Scour and deposition in rivers.” HEC-6 users manual Hydrol. Eng. Cen., DavisCalifHenderson, F. M., 1966: Open channelflow. Macmillan Pub. Co., New York, 522 p.Hey, R,D, 1978: “Determinate hydraulic geometry of river channels”. .1 Hydr. Div., ASCE,104 (6), 869-885.Hey, R.D., 1979: “Flow resistance in gravel-bed rivers.” I Hydr. Div., ASCE, 105 (4), 356-379.Hey, R.D., 1982: “Gravel-bed rivers: form and process.” In Gravel-Bed Rivers, Hey, R.D.,Bathurst, J.C., and Thorne, C.R., (Eds.). John Wiley and Sons., Chichester, England, 5-13.Hey, R.D., 1988: “Mathematical mdels of channel morphology.” In ModellingGeomorphological Systems, M.G. Anderson (Ed.), John Wiley and Sons, Chichester, England,99-125.Hey, R.D., and Thorne, C.R., 1986: “Stable channels with mobile gravel beds.” I Hydr. Div.,ASCE, 112, (8), 671-689.248Hickin, E.J., 1983: “River channel changes: retropsect and prospect”. Spec. Pub., Tnt. Assoc.Sed. 6, 61-83.Holden, A.P., and James, C.S., 1989: “Boundary shear distribution on flood plains.” .J Hydr.Res., 27(1), 75-89.Johnson, J.W., 1942: “The importance of considering side-wall friction in bedloadinvestigations.” Civil Eng. 12 (6), 329-331.Kamphius, J.W., 1974: “Determination of sand roughness for fixed beds.” I Hydr. Res., IAHR,12(2), 193-203.Kartha, V.C., and Leutheusser, H.J., 1970: “Distribution of tractive force in open channels.” IHydr. Div., ASCE, 96 (7), 1469-1483.Kellerhals, R., 1967: “Stable channels with paved gravel beds.” I Water. Harb. Div., ASCE,93 (1), 63 -83.Kellerhals, R., and Bray, D.I., 1971: “Sampling procedures for coarse fluvial sediments.” IHydr. Div., ASCE, 97(8), 1165-1180.Kellerhals, R., Neill, C.R., and Bray, D.I., 1972: “Hydraulic and geomorphic characteristics ofrivers in Alberta.” River Eng. and Surf Hydrol. Rept., 72-I, Res. Co. Alberta, Edmonton,Alberta., 54 p.Keulegan, G.H., 1938: “Laws of turbulent flow in open channels.” INat. Bur. Stand, 21, Res.Pap., 1151, 707-741.249Kirkby, M.J., 1977: “Maximum sediment efficiency as a criterion for alluvial channels”. InRiver channel changes, K.J. Gregory (Ed.), John Wiley and Sons, New York, 429-442.Knight, D.W., 1981: “Boundary shear in smooth and rough channels.” .1 Hydi. Div., ASCE,107, (7), 839-851.Knight, D.W., and Demetriou, J.D., 1983: “Flood plain and main channel flow interaction.” .1Hydr. Div., ASCE, 109 (8), 1073-1092.Knight, D.W., and Hamed, M.E., 1984: “Boundary shear in symmetrical compound channels.”J Hydr. Div., ASCE, 110 (10), 1412-1430.Knight, D.W., Demetrious, and J.D., Hamed, M.E., 1984: “Boundary shear in smoothrectangular channels.” I Hydi. Div., ASCE, 110 (4), 405-422.Knight, D.W., and Macdonald, J.A. 1979a: “Hydraulic resistance of artificial strip roughness.”I Hydr. Div., ASCE, 105 (6), 675-690.Knight, D.W., and Macdonald, J.A. 1979b: “Open channel flow with varying bed roughness.”I Hyck. Div., ASCE, 105 (9), 1167-1183.Konopinski, E.J., 1969: Classical descriptions of motion. W.H. Freeman and Co., SanFrancisco, USA, 504 p.Kuhnle, R.A., 1989: “Bed surface size changes in a gravel-bed channel.” I Hydr. Div., ASCE,115 (6), 73 1-743.250Kyrala, A., 1967: Theoretical Physics: applicattions of vectors matrice tensors andquaternions. W.B. Saunders Co., Philideiphia, USA, 359 p.Lacey, G., 1930: “Stable channels in alluvium.” Proc. ICE., 29, 259-292.Lane E.W., 1955a: “The importance of fluvial morphology in hydraulic engineering.” Proc.ASCE, 81, 1-17.Lane, E.W., 1955b: “The design of stable channels.” Trans. ASCE, 120, 1234-1279.Langbein, W.B., and Leopold, L.B., 1968: “River channel bars and dunes-theory of kinematicwaves.” Prof Pap. US Geol Sun’., 422L.Leighly, J.B., 1932: “Toward a theory of the morphologic significance of turbulence in the flowof water in streams.” Univ. Cal. Pub. Geog., 6, 1-22.Leopold, L.B., and Maddock, T. Jr., 1953: “The hydraulic geometry of stream channels andsome physiographic implications.” Prof Pap., U.S. Geol. Sun’., 252, 57 p.Leopold, L.B., Wolman, M.G., and Miller, J.P., 1964: Fluvial processes in geomorphology.W.H. Freeman, San Francisco, California, 522 p.Lighthill, M.J., and Whitham, G.B., 1955: “On kinematic waves. II. A theory of traffic flow onlong crowded roads.” Proc. Roy. Soc., 229 A, 317-345.Linsley, R.K., Franzini, J.B., Freyberg, D.L., and Tchobanoglous, G., 1992: Water-resourcesengineering, Fourth Edition. McGraw-Hill Inc., New York, USA, 841.251Lisle, T.E., and Madej, M.A., 1989: “Weak bed armouring in a channel with high bedloadsupply.” Trans. Am. Geoph. Union., 70 (15), 320.Lundgren, H., and Jonsson, I.G., 1964: “Shear and velocity distribution in shallow channels.” IHydr. Div., ASCE, 90 (1), 1-21.Mackin, J.H., 1948: “Concept of a graded river.” Bull. Geol. Soc. Am., 59 (5), 463-5 12.March, D.E., Abt, S.R., and Thorne, C.R., 1993: “Bank stability analyses versus fieldobservations.” In Hydr. Eng. ‘93, Proc. 1993 ASCE Hydr. Conf, Shen H.W., Su, S.T, andWen, F. (Eds), 88 1-886.Madej, M. A., 1978: “Sediment transport and channel changes in an aggrading stream in thePuget lowland, Washington.” Gen. Tech. Rept. US.D.A., PNW-141, 97-108.Meyer - Peter, E., and Muller, R., 1948: “Formulas for bed-load transport.” Rept. 2ndMeeting, IAHR, Stokholm, Sweden, 39-64.Meyers, R.C., 1978: “Momentum transfer in a compound channel.” I Hydr. Res., IAHR, 16(2), 139-150.Meyers, R.C., and Elsawy, E.M., 1975: “Boundary shear in channel with floodplain.” I Hydr.Div., ASCE, 101 (7), 933-946.Milhous, R.T., 1972: Sediment transport in a gravel-bottomed stream. Unpub. Ph.D. thesis,Oregon State University, Corvallis, Oregon, 232 p.252Mithous, R.T., and Klingeman, P.C., 1973: “Influence of the bed surface layer in gravel-bedriver dynamics.” In Hydraulic engineering and the environment. Proc. 21 Annual Hydr. Div.Spec. Conf, ASCE, Bozema, Montana, 293-303.Millar, R.G., 1991: Development of an analytical model of river response. Unpub. M.A.Scthesis, University ofBritish Columbia, Vancouver, Canada, 135 p.Millar, R.G., and Quick, M.C., 1993a: “An optimization approach to river adjustments.” InHydr. Eng. ‘93, Proc. 1993 ASCE Hydr. Conf, ShenH.W., Su, S.T, and Wen, F. (Eds), 1635-1640.Millar, R.G., and Quick, M.C., 1993b: “Effect of bank stability on geometry of gravel rivers.”I Hydr. Div., ASCE, 119 (12), 1343-1363.Morgenstern, N.R., 1963: “Stability charts for earth slopes during rapid drawdown.” Geotech.,13 (2), 121—13 1.Morgenstem, N.R., and Price, V.E., 1965: “The analysis of the stability of general slipsurfaces.” Geotech., 15(1), 79-93.Nanson, G.C., and Erskine, W.D., 1988: “Episodic changes of channels and floodplains oncoastal rivers in New South Wales.” In Fluvial Geomorphology of Australia, R.F. Warner(Ed), Academic Press, Sydney, 201-221.Nash, D., 1987: “A comparative review of limit equilibrium methods of stability analysis.” InSlope Stability, Anderson, M.G., and Richards, K.S. (Eds), John Wiley and Sons, Chichester,England, 11-75.253Neill, C.R., 1968: “A re-examination of the beginning of movement of coarse granular bedmaterials.” Hydr. Res. Stn. Rept., No. TNT 68, 37p.Newson, M.D., 1986: “River basin engineering - fluvial geomorphology.” J Inst. Wa. Eng.Sd., 40 (4), 307-324.Nikuradse, 3., 1933: “Stromungsgesetze in rauhen rohren.” (Laws of flow in rough pipes),Forschungshefl verein deutscher ingenieure, No. 361.Nixon, M., 1959: “A study of bankfi.ill discharges of the rivers in England and Wales.” Proc.,Inst. Civil Eng., London, England, 12, 157-174.Noutsopoulos, G.C., and Hadjipanos, P.A., 1982: Discussion of “Boundary shear in smoothand rough channels,” by D.W. Knight, I Hydr. Div., ASCE, 108 (6), 809-8 12.Osman, A.K., and Thorne, C.R., 1988: “Riverbank stability analysis. I: Theory.” I Hydr. Div.,ASCE, 114 (2), 134-150.Park, CC., 1977: “Man-induced changes in stream channel capacity.” In River channelchanges, Gregory, K.J. (Ed), John Wiley and Sons., Chichester, England, 12 1-144.Parker, G., 1978: “Self-formed rivers with equilibrium banks and mobile bed. Part 2. The gravelriver.” I FluidMech., 89 (1), 127-146.Parker, G., 1990: “Surface-based bedload transport relation for gravel rivers.” I Hydr. Res.,20(4), 417-436.254Parker, G., and Peterson, A.W., 1980: “Bar resistance of gravel bed rivers.” .] Hydr. Div.,ASCE, 106 (10), 1559-1575.Parker, G., and Klingeman, P.c., 1982: “On why gravel-bed streams are paved.” Wat. Res.Res., 18 (5), 1409-1423.Parker, G., Klingeman, P.C., and McLean, D.G., 1982a: “Bedload and size distribution inpaved gravel-bed streams.” J Hydr. Div., ASCE, 108 (4), 544-571.Parker, P., Dhamotharan, S., and Stefan, H., 1982b: “Model experiments on mobile, pavedgravel-bed streams.” Wat. Res. Res., 18 (5), 1395-1408.Pasche, E., and Rouve, G., 1985: “Overbank flow with vegetatively roughened flood plains.” IHydr. Div., ASCE, 111 (9), 1262-1278.Phinney, R., 1975: Stage-discharge relations for gravel-bed rivers. Unpub. M.Eng. Report,University ofNew Brunswick, Fredericton, NB, Canada, 165 p..Posey, C.J., 1967: “Computation of discharge including over-bank flow.” Civil Eng., 37(4), 62-63.Prestegaard, K.L., 1983: “Bar resistance in gravel bed streams at bankfull discharge.” Wat. Res.Res., 19 (2), 472-476.Rajaratnam, N., and Muralidhar, D., 1969: Boundary shear stress distribution in rectangularopen channels.” Houille Blanche, 24 (6), 603-609.Raudkivi, A.J., 1990: Loose Boundary Hydraulics, 3rd Ed. Permagon Press, New York, 538 p.255Raynov, S., Pechinov, D., and Kopaliany, Z., 1986: River response to hydraulic structures.UNESCO, Paris, France, 115 p.Roberts, R.G., and Church, M., 1986: “The sediment budget in severly disturbed watersheds,Queen Charlotte Ranges, B.C.” Can. I For. Res., 16(5), 1092-1106.Schumm, S. A., 1969: “River metamorphosis.” I Hydr. Div., ASCE, 95 (1), 255-273.Schumm, S.A., 1977: The Fluvial System. John Wiley and Sons, New York, 338 p.Schumm, S.A., and Lichty, R.W., 1965: “Time, space and causality in geomorphology.” Am. ISci., 263 (2), 110-119.Sellin, R.H.J., 1964: “A laboratory investigation into the interaction between flow in thechannel of a river and that over its flood plain.” Houille Blanche, 19(7), 793-802.Shaw, E.M., 1988: Hydrology in practice, 2nd Ed. Van Nostrand Reinhold Int, London, 539 p.Shen, H.W., 1962: “Development of bed roughness in alluvial channels.” I Hydr. Div., ASCE,88 (3), 45-58.Simons, S.B., and Albertson, M.L., 1963: “Uniform water conveyance channels in alluvialmaterial.” Trans., ASCE, 128 (1), No. 3399, 65-167.Singh, B., 1961: “Bed load transport in channels.” Irrig. and Power, J. Cent. Board of Irrig.Power, 18(5), 41 1-430.256Song, C.C.S. and Yang, C.T., 1979: “Velocity profiles and minimum stream power.” .1 Hydr.Div., ASCE, 105 (8), 981-998.Spangler, M.G., and Handy, R.L., 1982: Soil Engineering, 4th Ed. Harper and Row, NewYork, 819p.Spencer, E., 1967: “A method of analysis for stability of embankments using parallel inter-sliceforces., Geotech., 17(1), 11-26.Taylor, D.W., 1948: Fundamentals of Soil Mechanics. John Wiley and Sons, New York, 700p.Terzaghi, K., and Peck, R.B., 1948: Soil Mechanics in Engineering Practice. John Wiley andSons, New York, 729 p.Thomas, W.A., and Prasuhn, A.L., 1977: “Mathematical modeling of scour and deposition.” JHydr. Div., ASCE, 103 (8), 85 1-863.Thorne C.R., 1982: “Processes and mechanisms of river bank erosion” In Gravel-Bed Rivers,Hey, R.D., Bathurst, J.C., and Thorne, C.R., (Eds.). John Wiley and Sons., Chichester,England, 227-271.Thorne C.R., Murphey, J.B., and Little, W.C., 1981: “Bank stability and bank materialproperties in the Bluffline streams of Northwest Mississippi.” Rept. US Army Eng., USDASed. Lab., Oxford, Miss.257Thorne, C.R., Hey, RD., and Chang, H.H., 1988: “Prediction of hydraulic geometry of gravel-bed streams using the minimum streampower concept.” In International Conference on RiverRegime, W.P. White (Ed.). John Wiley and Sons., Chichester, England, 29-40.Thorne, C.R., and Osman, A.M., 1988: “Riverbank stability analysis. II: Applications.” I Hydr.Div., ASCE, 114 (2), 15 1-172.Vanoni, V.A., 1975: Sedimentation engineering. A.S.C.E. Manuals and Reports onEngineering Practice, 54, 745p.Vanoni, V.A., and Hwang, L. S., 1967: “Relation between bed forms and ffiction in streams.”I Hydr. Div., ASCE, 93 (3), 121-144White, C.M., 1940: “Equilibrium of grains on the bed of a stream.” Proc. Roy. Soc. Lon. (A),174, 322-338.White, W.R., Bettess, R., Paris, E., 1982: “An analytical approach to river regime.” .1 Hydr.Div., ASCE, 108 (10), 1179-1193.Williams, G.P., 1978: “Bank-fi.ill discharge of rivers.” Wat. Res. Res., 14(6), 1141-1154.Wolman, M.G., 1955: “The natural channel of Brandywine creek.” Prof Pap., U.S. Geol.Sun’., 271 p.Wolman, M.G., 1959: “Factors influencing the erosion of cohesive river banks.” Am .1 Sci.,257 (3), 204-2 16.258Wolman, M.G., 1967: “A cycle of sedimentation and erosion in urban river channels.”GeografislcaAnnaler, 49A, 385-395.Wolman, M.G., and Brush, L.M., 1961: “Factors controlling the size and shape of streamchannels in coarse non-cohesive sands.” Prof Pap., U.S. Geol. Sun’., 282-G,Wolman, M.G., and Leopold, L.B., 1957: “River flood plains: Some observations on theirformation.” Prof Pap., U.S. Geol. Sun’., 282-C.Wolman, M.G., and Williams, G.P., 1984: Downstream effects of dams on alluvial rivers. USGov. Print. Office, Washington, DC.Womack W.R., and Schumm, S.A., 1977: “Terraces on Douglas Creek, north-westernColorado: an example of episodic erosion”. Geology, 5, 72-76.Yang C.T., 1971: “On river meanders.” J Hydrol., 13 (3), 231-235.Yang, C.T., 1976: “Minimum unit stream power and fluvial hydraulics.” I Hydr. Div., ASCE,102 (7), 919-934.Yang, C.T., and Song, C.C.S., 1979: “Theory of minimum rate of energy dissipation.” I Hyd.r.Div., ASCE, 105 (7), 769-784.Yang, C.T., Song, C.C.S., and Woldenberg, M.J., 1981: “Hydraulic geometry and minimumrate of energy dissipation.” Water Res. Res., 17(4), 1014-1018.259Yang, C.T., and Song, C.C.S., 1986: “Theory of minimum energy and energy dissipation rate.”In Encyclopedia of fluid mechanics, N.P. Cheremisinoff (Ed.), Vol. 1, Chapt. 11, GulfPublishing Co., Houston, Texas, 3 53-399.260BANKFULL: FIXED SLOPE MODELThe Bankfhll:Fixed Slope model is the simplest version of the optimization model presented inthis thesis. In this version of the model the channel slope is treated as an independent variable.The source code for the computer program (fixslope.bas) used in the thesis is presented below.The programming was encoded and run using Microsoft Quick Basic (Version 3.0 or later).This program can be used for channels with cohesive or noncohesive bank sediment.The program is designed to use an input data file named fixslope.dat, and will output theoptimal geometry in a data file named fixslope.out. The data file can contain the data for morethan one channel. The data for each channel must be input as follows:Qbf S ksbd ksbk d50 D50 D5Obkfor channels with noncohesive banks, and:Qbf S ksbd ksbk d50 D50for channels with cohesive banks.The data relating to each channel must start on a new line. The data values must be separatedby at least one space, although the number of spaces is not important. An example of an inputdat file is presented below:50 0.003100 0.0030.1 0.1 0.025 0.075 0.075 400.1 0.1 0.025 0.075 0.05 40261APPENDIX A‘isisrt C’ TcritThis file contains the input data for two channels with noncohesive banks. The user is requiredto input the bank sediment type (n for noncohesive or c for cohesive) as prompted. The outputdata will written to a data file named fixslope.out in the following order:W H 17* S 0The output file resulting from the example input file is:W H * S theta18.0 1.40 1.21 0.0030 30.037.7 1.28 1.17 0.0030 21.7The Bankfiull:Fixed Slope model was developed principally for illustrative purposes. It can beused in a field investigation to estimate or ‘calibrate’ the bank stability parameters that may bedifficult to measure. For example with noncohesive bank sediment the value of Ø., may bedifficult to estimate in the field. The measured independent variables can be input into themodel together with trial values of until the modelled width is equal to the observed channelwidth. This procedure was used in Chapter 6 to demonstrate the relationship between andbank vegetation. The value of which gives the best agreement between the modelled andobserved geometry can then be used as input into the fully developed model (Appendix C).The source code for the Bankfull: Fixed Slope Model is presented below.262‘BANKFULL:FIXED SLOPE MODEL FOR BOTH COHESIVE AND NONCOHESIVE BANK‘SEDIMENT. (FIXSLOPE.BAS)DECLARE SUB stabcurve 0DECLARE SUB bankstabilitycohesive 0DECLARE SUB optimum 0DECLARE SUB bankstabilitynoncohesive 0DECLARE SUB ebrown 0DECLARE FUNCTION bankshear! 0DECLARE FUNCTION bedshear! 0DECLARE FUNCTION shearforce! 0DECLARE SUB continuity 0DECLARE FUNCTION velocity! 0DECLARE FUNCTION hydrad! 0DECLARE FUNCTION meandepth! 0DECLARE FUNCTION area? 0DECLARE FUNCTION surfwidth! 0DECLARE FUNCTION depth? 0I*********************************************************************************program to calculate curves to demonstrate the optimumchannel geometiy for a given set of independent variables‘bankfull: fixed slope model for channels with‘both cohesive and noncohesive bank sediment‘Global Variables‘gravity!= gravitational constant gamma? unit weight of water‘ss! specific gravity of sediment viscosity?=kinematic viscosity of water‘density!=density of water slope?=channel slope‘pbed!=bed perimeter pbank?—bank perimeter‘theta!=bank angle (radians) dd5O?=median armour grain diameter‘d50!median subarmour grain diameter d5Obank?=median grain diameter of bank sediment‘discharge!=bankfull discharge ‘gbcalc!=valculated bedload transport rate ‘neta!=coefficient ofefficiency‘phipnme?=friction angle of bank sediment (radians)‘phidegrees!=friction angle of bank sediment (degrees)‘ksbed?=measure of bed roughness ksbank!=measure of bank roughness‘banktype$=type of bank sediment Wointerger counter‘gammat!=unit weight of cohesive bank sediment‘stabnum’=stability number for cohesive sediment‘cohesion!=soil cohesion taucrit?=critical shear stressCOMMON SHARED gravity?, gamma!, ss?, viscosity?, density!COMMON SHARED slope!, pbed?, pbank?, theta!, ddSO’, d5Obank?COMMON SHARED discharge!, gbcalc!, phiprime!COMMON SHARED d50!, bankcond$, pi?COMMON SHARED ksbed?, ksbank!, neta!, bankte$, f%COMMON SHARED ganimat?, stabnum?, cohesion!, taucrit?, phidegrees?OPEN “fixslope.dat” FOR INPUT AS #1OPEN “flxslope.out” FOR OUTPUT AS #2263‘SET VALUES OF INDEPENDENT CONSTANTSpi! = 3.14159gravity!= 9.81gamma! = 9810viscosity! = .000001density! = 1000ss! = 2.65CLSPRINTPRINTPRINT “BANKFULL:FIXED SLOPE OPTIMIZATION MODEL”PRINT : INPUT “Hit to Continue “,dummy$=0PRINT #2,” W H * S theta”CLS10PRINTPRINTPRINTPRiNT : INPUT “INPUT BANK SEDIMENT TYPE (n/c) “; banktype$iF banktype$ = “c” OR banktype$ = “C” THENbanktype$ = “cohesive”ELSEIF bankte$ = “N” OR banktype$ = “n” THENbanktype$ = “noncohesive”ELSEBEEPPRINT “Bank Type Unknown”GOTO 10END IFCLSDO WHILE NOT EOF(1)PRINTPRINT=1% + 1 ‘counterPRiNT “Channel Number“; f’%iNPUT #1, discharge!, slope!, ksbed!, ksbank!, d50!, dd5O!‘Input the independent variables from data fileIF banktype$ “noncohesive” THEN ‘input bank stability parametersINPUT #1, d50bank!, phidegrees!ELSE ‘banktype$=”cohesive”INPUT #1, gamxnat!, phidegrees!, cohesion!, taucrit!END IF264phiprime! = phidegrees! * 2 * pi! / 360 ‘convert friction angle from degrees to radiansCALL optimum‘calculates the optimal geometryPRINTPRINT” W H * S theta”PRINT USiNG “###.# ##.## ##.## #.#### ##.#“; surfwidth!; depth!; meandepth!; slope!; theta! * 360/2 /pi!PRINT #2, USING “###.# ##.## ##.## #.#### ##.#“; surfwidth!; depth!; meandepth!; slope!; theta! * 360 / 2/ pi!LOOP‘loops until end of data fileCLOSE #1CLOSE #2ENDFUNCTION area‘function to calculate the cross-sectional areaarea! = .5 * (pbed! + surfwidth!) * depth!END FUNCTIONFUNCTION bankshear‘calculates the bank shear stressbankshear! = gamma! * depth! * slope! * shearforce! * ((surfwidth! + pbed!) * SIN(theta!) / (4 * depth!))END FUNCTIONSUB bankstabilitycohesive‘satisfies bank stability constraint for cohesive banksheightcond$ = “unknown”shearcond$ = “unknown”bankcond$ = “unknown”thetamax! =90*2 *pi! /360thetanün! =0‘initialise search and convergence criteria‘First test ifvertical bank is stable wrt bank heighttheta! = thetamax!CALL continuity‘satisfy continuity265CALL stabcurve‘catculate stability numbercriticalheight! = stabnuml * cohesionl I gaxnmat’calculate the critical heightstability 1! = depth / criticaiheight!‘calculates stability criteria wit bank height‘if stability 1 I <= criticaiheight then bank is stable‘with respect to bank heightIF stabilityll <= 1.001 THENheightcond$ = “stable”ELSE ‘vertical bank is not stable therefore reduce the bank angle.‘Determine the maximum bank angle that is just stable.This is theta max from Chapter 6.DO UNTIL heightcond$ = ‘just stable”theta1 = (thetamax? + thetamin’) /2‘calculate midpoint of range for‘bisectrix convergence schemeCALL continuity‘satisfy continuity constraintCALL stabcurve‘calculate stability numbercriticaiheight! = stabnum! * cohesion! / gaimnat!stabilityll = depth / criticaiheight’‘calculates critical height and stabilty number‘for the bank height constraintIF stabilityl! >= 1.001 THENheightcond$ = “unstable”thetamax! theta?ELSEIF stabilityll <= .999 THENheightcond$ = “understable”thetamin! = theta!ELSEheightcond$ = ‘just stable”‘Bank Height = Critical HeightEND IFIF (thetamax! - thetamin’) / theta! <.001 THENheightcond$ = ‘Just stable”stabilityl! = 1‘Secondary Convergence Criterion in case of‘convergence problems due to numerical schemeEND IFLOOP266END IFbank height constraint now satisified‘now assess the bank shear constraintthetaniin?=0thetamax! = theta!‘thetamax is reset to maximum bankangle which satisfied the‘bank height constraint abovestability2! = bankshear! / taucritsatbility criterion for bank shear constraintIF stability2! < 1.001 THEN‘bank shear constraint is satisifledshearcond$ = “stable”bankcond$ = “stable”ELSE ‘must reduce theta!DO UNTIL shearcond$ = “just stable”theta! = (thetamax! + thetaniin!) /2‘calculate mid pointIFtheta! < 10*2*3.14159/36OTHEN‘bank angles less than 10 degrees‘nominally assumed to be unstableshearcond$ = “unstable”bankcond$ = “unstable”EXIT DOEND IFCALL continuity‘satisfy continuitystability2! = bankshear! / taucrit!‘calculate stability criterion for bank shear constraintIF stability2! >= 1.001 THEN‘bankshear > taucritshearcond$ = “unstable”thetamax! = theta!ELSEIF stability2! <= .999 THENshearcond$ = “understable”thetamin! theta!ELSEshearcond$ = “just stable”bankcond$ = “stable”267END IFIF (thetaxnax! - thetanun!) / theta! <.001 THENshearcond$ ‘just stable”bankcond$ = “stable”stability2? = 1END IFLOOPEND IFEND SUBSUB bankstabilitynoncohesive‘calculates theta where banks just stablebankcond$ = “unknown”thetamax! = phiprime!thetamin! =0‘initialise search bounds and convergence criteriaDO UNTIL bankcond$ = ‘just stable”theta! (thetaniax! + thetaniin!) / 2‘determine midpoint of search rangeIF theta! <5 * 2 * 3.14159/360 THENbankcond$ = “unstable”EXiT DO‘if theta is less than 5 degrees the channel is assumed‘to be unstableEND IFCALL continuity‘satisfies continuity constraintstability! = (bankshear! / (gamma! * (ss! - 1) * d5obank!)) I (.048 * TAN(phiprime!) * (1 - SIN(theta!) A 2 /SIN(phiprime!) A 2) A .5)‘calculates stability criterionIF stability! >= 1.001 THENbankcond$ = “unstable”thetamax! = theta!ELSEIF stability! <= .999 THENbankcond$ = “stable”thetamin! = theta!ELSEbankcond$ = ‘just stable”‘primaiy convergence criterionEND IF268IF (thetamax! - thetamin!) I theta! <.001 THENbankcond$ = ‘just stable”stability! = 1‘secondaiy convergence criterionEND IFLOOPEND SUBFUNCTION bedshear‘calculates the bed shear stressbedshear! gamma! * depth! * slope! * (1 - shearforcel) * ((surfwidth! / (2 * pbed!) + .5))END FUNCTIONSUB continuity‘varies pbank! for trial values of Pbed!, slope!, and theta! to‘satisfy the contiuity constraintpbankcond$ = “unknown”errorcalc! = 1000minpbank! =0maxpbank! =20 * discharge! A •35‘initialise search and convergence criteriaDO UNTIL pbankcond$ = “OK”pbank! = (minpbank! + maxpbank!) /2‘calculate midpointerrorcalcl = (area * velocity / discharge!)‘calculate the normalised errorIF errorcaic! > 1.001 THENpbankcond$ “too large”maxpbank! = pbank!ELSEIF errorcalc! <.999 THENpbankcond$ = “too small”minpbank! = pbank!ELSEpbankcond$ = “OK”END IFIF (maxpbank! - minpbank!) / pbank! <.0001 AND pbankcond$ = “too small” THEN‘resets maxpbank! if too smallmaxpbank! =2 * maxpbank!END IF269LOOPEND SUBFUNCTION depth‘function to calculate flow depth of a trapezoidal channeldepth! = .5 * SIN(theta?) * pbank!END FUNCTIONSUB ebrown‘calculate the sed trans capacity using em-brown formuladimlessbedshear! = bedshear! / (gamma! * (ss! - 1) * d50!)IF dinflessbedshear! <=0 THENdimlessbedload! =0ELSEIF dimlessbedshear! <.093 THENdimlessbedload! = 2.15 * EXP(-.391 / dimlessbedshear!)ELSEdimlessbedload! =40 * dimlessbedshear! “3END IFfi!= (2 / 3 + 36 * viscosity? “2/ (gravity! * d50! “3 * (ss? - 1))) “.5 -(36 * viscosity? “2/ (gravity! * d.50? “3*(ss!_1)))gbcalc! = pbed! * fi! * ss? * density! * ((ss! - 1) * gravity! * d50! A 3) A 5 * dimlessbedload!END SUBFUNCTION hydrad!‘function to calculate the hydraulic radiushydrad! = area! / (pbed? + pbank!)END FUNCTIONFUNCTION meandepth‘function to calculate the mean depthmeandepth! = area / surfwidthEND FUNCTIONSUB optimum‘determines the optimal geometzyoptcond$ = “unknown”‘initialises optimality test condition270lowerpbed! =0upperpbed? = 10 * discharge! “.5‘set bounds for Pbed based on Regime Eqns‘typical optimal value of Pbed is 2 to 5 * discharge” .5minpbed! = lowerpbed!maxpbed! = upperpbed!PRINTDO WHILE optcond$ <> “optimum”pbed! = (minpbed! + maxpbed!) /2‘calculate mid point of search rangeIF pbed! > .95 * upperpbed! THENupperpbed!2*upperpbed!maxpbed! = upperpbed!‘reset upper bound of searchEND IFPRINT “Assessing Trial Bed Perimeter “;PRiNT USiNG “####.##“; pbed’;PRiNT “ m “;pbed! = pbed! * .975‘calculate backwards difference value of PbedIF banktype$ = “noncohesive” THENCALL bankstabilitynoncohesive‘solve bank stability constraint for noncohesive banksELSE ‘banktype$ = “cohesive”CALL bankstabilitycohesive‘solve bank stability constraint for cohesive banksEND IFIF bankcond$ <> “unstable” THEN‘test conditions indicate that stable channel geometry‘has been determined and the bedload constraint‘has been satisifiedCALL ebrown‘calculate sdiment tarnsporting capacity of the channelnetal! = gbcalc! / (density! * discharge! * slope!)‘evaluate neta for backward difference pointpbed! pbed! 1.975 * 1.025‘calculate pbed for forward duffIF banktype$ = “noncohesive” THENCALL bankstabilitynoncohesive‘solve bank stability constraint for noncohesive banksELSE ‘banlctype$ = “cohesive”CALL bankstabilitycohesive‘solve bank stability constraint for cohesive banksEND IFCALL ebrown271‘calculate sdiment tarnsporting capacity of the channelneta2! = gbcalcl / (density’ * discharge! * slope!)‘calculates neta for forward difference valuepbed! =pbed! / 1.025‘reset pbed to midpoint valuednetabydpbed! = (neta2! - netal!) I (.05 * pbed!)‘calculate first derivative by finite difference‘to assess optimum conditionIF dnetabydpbed? <0 THEN ‘Trial Pbed is too demaxpbed! = pbed’ ‘Reduce upper bound of Pbedoptcond$ = “too wide”ELSEIF dnetabydpbed! >0 THEN ‘Trial Pbed is too narrowminpbed! = pbed! ‘Increase lower bound of Pbedoptcond$ = “too narrow”ELSEoptcond$ = “optimum” ‘Optimum AchievedEND IFIF (maxpbed! - minpbed!) I pbed! <.001 THEN optcond$ “optimum”‘Convergence attained. Second optimality criterionELSE ‘trial Pbed too small for stable geometiy: Increase minimum Pbedpbed! =pbed!/.975‘reset from backward difference to correct valueminpbed! = pbed!‘set lower limit at current trial value‘(as the optimum value must be greater)optcond$ = “too narrow”END IFPRiNT optcond$LOOP‘Evaluate geometiy at exact optimal Pbed‘(not at forward or backward difference value)IF banktype$ = “noncohesive” THENCALL bankstabilitynoncohesive‘solve bank stability constraint for noncohesive banksELSE banktype$ = “cohesive”CALL bankstabilitycohesive‘solve bank stability constraint for cohesive banksEND IFCALL ebrown‘calculate sdiment tarnsporting capacity of the channelneta! = (pbed! * gbcalc!) / (density! * discharge! * slope!)272‘evaluate objective functionEND SUB‘*********************************************************************************FUNCTION shearforce‘calculates the proportion of shear force on the banksshearforce! = 1.766 * (pbed! I pbank! + 1.5) A -1.4026END FUNCTION‘*********************************************************************************SUB stabcurve‘stability curves are from Figure 5.3thetadegrees! =theta! /2/3.14159 * 360IF phidegrees! <20 THEN‘use phiprime = 15 degree stability curve from Figure 5.3‘from Taylor (1948)SELECT CASE thetadegrees!CASE 48.65 TO 90stabnum! = -.141 * thetadegrees! + 17.46CASE 29.36 TO 48.6499stabnum! = -.533 * thetadegrees! + 36.53CASE 20.12 TO 29.3599stabnum! = -3.07 * thetadegrees! + 111CASE 15.001 TO 20.11999stabnum! = -189.6 * thetadegrees! + 3844CASEIS<= 15stabnum! = 1000000END SELECTELSE‘use phiprime =25 degrees curve from Figure 5.3‘After Taylor, (1948)SELECT CASE thetadegrees!CASE 58.62 TO 90stabnum! = -.212 * thetadegreesi + 24.86CASE 39.11 TO 58.61999stabnum! = -.8975 * thetadegrees! + 65.04CASE 25 TO 39. 10999stabnuml = -6.97 * thetadegrees! + 302.5273CASE IS <=25stabnuml = 1000000END SELECTEND IFEND SUB‘*********************************************************************************FUNCTION surfwidth‘function to calculate the surface widthsurfwidth! = pbed! + COS(theta?) * pbank?END FUNCTIONFUNCTION velocity!‘calculates the mean velocityrhbank! = bankshear! I (gamma! * slope!)rhbed! = bedshear! I (gamma! * slope!)Thank! = (2.03 * LOG(12.2 * rhbank! I ksbank!) I LOG(10)) A -2Ibed! = (2.03 * LOG(12.2 * rhbed! / ksbed!) I LOG(10)) A -2ifactor! = (ibed! * pbed! I (pbed! + pbank!) + Ibank! * pbankl / (pbank! + pbed!)) Avelocity! = (8 * gravity! * hydrad! * slope!) A 5 * ifactor!END FUNCTION274APPENDIX BBANKFULL: VARIABLE SLOPE MODELThe Bankfull:Variable Slope differs from the Fixed Slope model in that the channel slope isnow treated as a dependent variable.The source code for the computer program (varslope.bas) used in the thesis is presented below.The programming was encoded and run using Microsoft Quick Basic (Version 3.0 or later).This program can be used for channels with cohesive or noncohesive bank sediment.The program is designed to use an input data file named varslope.dat, and will write the outputto a data file named varslope.out. The data file can contain data for more tahn one channel. Thedata for each channel must be input as follows:Qbf Gbf ksd ksbk d50 D50 D5Obk sfor channels with noncohesive banks, and:Qbf Gbf ksbd k,bk d50 D50 rt c’ rcritfor channels with cohesive banks.The data relating to each channel must start on a new line. The data values must be separatedby at least one space, although the number of spaces is not important. An example of an inputdat file is presented below:50 5 0.1 0.1 0.025 0.075 0.075 40100 20 0.1 0.1 0.025 0.075 0.05 40275This file contains the input data for two channels with noncohesive banks. The user is requiredto input the bank sediment type (n for noncohesive or c for cohesive) as prompted. The outputdata will written to a data file named varslope.out in the following order:W H 17* S 6The output file resulting from the example input file is:W H * S theta16.5 1.66 1.39 0.0023 32.035.3 1.39 1.26 0.0027 22.2The Bankfiill:Variable Slope model was developed principally for illustrative purposes. Thesource code for the Bankfull:Variable Slope Model is presented below.276BANKFULL:VAR1ABLE SLOPE MODEL FOR BOTH COHESIVE AND NONCOHESIVE BANKSEDIMENT‘(VARSLOPE.BAS)DECLARE SUB stabcurve 0DECLARE SUB bankstabilitycohesive 0DECLARE SUB bedload2 0DECLARE SUB optimum 0DECLARE SUB bankstabilitynoncohesive 0DECLARE SUB bedload 0DECLARE SUB ebrown 0DECLARE FUNCTION bankshear! 0DECLARE FUNCTION bedshear? 0DECLARE FUNCTION shearforce? 0DECLARE SUB continuity 0DECLARE FUNCTION velocity! 0DECLARE FUNCTION hydrad! 0DECLARE FUNCTION meandepth! 0DECLARE FUNCTION area? 0DECLARE FUNCTION surfwidth! 0DECLARE FUNCTION depth! ()program to calculate curves to demonstrate the optimumchannel geometiy for a given set of independent variables‘bankfull: variable slope model for channels with‘both cohesive and noncohesive bank sediment‘Global Variables‘gravity!= gravitational constant gamma?= unit weight of water‘ss1= specific gravity of sediment viscosity!=kinematic viscosity of water‘density!=density of water slope!=channel slope‘pbedl=bed perimeter pbank!=bank perimeter‘theta!=bank angle (radians) dd5O!=median armour grain diameter‘d50!inedian subarinour grain diameter d5Obank!=median grain diameter of bank sediment‘discharge!=bankfull discharge gbload!=imposed bedload transport rate‘gbload!=calculated bedload transport rate ‘neta!coefficient of efficiency‘phiprime!=friction angle of bank sediment (radians)‘phidegrees!=friction angle of bank sediment (degrees)‘ksbcd!=measure of bed roughness ksbank!measure of bank roughness‘banktype$=type of bank sediment Wo=interger counter‘gammat!=unit weight of cohesive bank sediment‘stabnuxn!=stability number for cohesive sediment‘cohesion!=soil cohesion taucrit?=cntical shear stressCOMMON SHARED gravity!, gamma!, ss!, viscosity?, density?COMMON SHARED slope!, pbed?, pbank!, theta!, dd5O!, d5Obank!COMMON SHARED discharge!, gbcalc!, gbload!, phiprime?COMMON SHARED d50!, bankcond$, pi!COMMON SHARED ksbed!, ksbank!, neta!, slopecond$, banktype$, 1%COMN’ION SHARED ganimat!, stabnum!, cohesion!, taucrit?, phidegrees!277OPEN “varslope.dat” FOR iNPUT AS #1OPEN “varslope.out” FOR OUTPUT AS #2‘SET VALUES OF INDEPENDENT CONSTANTSpi! = 3.14159gravity! = 9.81gamma! = 9810viscosity! = .000001density! = 1000ss! =2.65CLSPRINTPRINTPRINT “BANKFULL:VARIABLE SLOPE OPTIMIZATION MODEL”PRINT : INPUT “Hit to Continue “, dummy$=0PRINT #2,” W H * S theta”CLS10PRiNTPRINTPRINTPRINT : INPUT “INPUT BANK SEDIMENT TYPE (a/c) “;banktype$IF banktype$ = “c” OR banktype$ = “C” THENbanktype$ = “cohesive”ELSEIF banktype$ = “N” OR banktype$ = “n” THENbanktype$ = “noncohesive”ELSEBEEPPRINT “Bank Type Unkno”GOTO 10END IFCLSDO WHILE NOT EOF(l)PRINTPRINTI% =1% + 1 ‘counterPRiNT “Channel Number “; WeINPUT #1, discharge!, gbload!, ksbed!, ksbank!, d50!, dd5O!‘Input the independent variables from data fileIF banktype$ = “noncohesive” THEN ‘input bank stability parametersINPUT #1, d50bank!, phidegrees!ELSE ‘banktypeS=”cohesive”iNPUT #1, ganunat!, phidegrees!, cohesion!, taucrit!END IF278phiprimel = phidegrees! * 2 * pit / 360 ‘convert friction angle from degrees to radiansCALL optimum‘calculates the optimal geometiyPRiNTPRINT” W H * S theta”PRiNT USiNG “##*.# #1.## #.## #.## #.#“; surfwidth!; depth!; meandepth!; slope!; theta! * 360 / 2 /pi!PRINT #2, USING “###.# ##.## #.#### ##.#“; surfwidth!; depth!; meandepth!; slope!; theta! * 360 /2/ pi!‘prints out the values of the dependent variablesLOOP‘loops until end of data fileCLOSE #1CLOSE #2ENDFUNCTION area‘function to calculate the cross-sectional areaarea! = .5 * (pbed! + surfwidth!) * depth!END FUNCTIONFUNCTION bankshear‘calculates the bank shear stressbankshear! = gamma! * depth! * slope! * shear.force! * ((surfwidth! + pbed!) * SIN(theta!) / (4 * depth!))END FUNCTIONSUB bankstabilitycohesive‘satisfies bank stability constraint for cohesive banksheightcond$ = “unknown”shearcondS = “unknown”bankcond$ = “unknown”thetaniax! =90*2*pi!/360thetamin! =0‘imtialise search and convergence criteria‘First test ifvertical bank is stable wrt bank heighttheta! = thetamax!CALL continuity‘satisfy continuityCALL stabcurve‘catculate stability number279criticalheight! = stabnum? * cohesion! I ganunat!calculate the critical heightstability!! = depth I criticalheight!‘calculates stability criteria wrt bank height‘if stabilityl! <= criticaiheight then bank is stable‘with respect to bank heightIF stabilityl! < 1.001 THENheightcond$ = “stable”ELSE ‘vertical bank is not stable therefore reduce the bank angle.‘Determine the maximum bank angle that is just stable.‘This is theta max from Chapter 6.DO UNTIL heightcond$ = ‘just stable”theta! = (thetamax! + thetanun!) /2‘calculate midpoint of range for‘bisectrix convergence schemeCALL continuity‘satisfy continuity constraintCALL stabcurve‘calculate stability numbercriticaiheight! = stabnum! * cohesion! / ganunat!stability 1! = depth / criticaiheight!‘calculates critical height and stabilty numberTor the bank height constraintIF stabilityl! >= 1.001 THENheightcond$ = “unstable”thetamax! =theta!ELSEIF stabilityl! < .999 THENheightcond$ = “understable”thetamin! = theta!ELSEheightcond$ = “just stable”‘Bank Height = Critical HeightEND IFIF (thetamax! - thetamin!) / theta! <.001 THENheightcond$ = “just stable”stability!! = 1‘Secondary Convergence Criterion in case of‘convergence problems due to numerical schemeEND IFLOOPEND IF280‘bank height constraint now satisified‘now assess the bank shear constraintthetamin! =0thetaniax! = theta!‘thetamax is reset to maximum bankangle which satisfied the‘bank height constraint abovestability2! = bankshear! / taucrit!‘satbility criterion for bank shear constraintIF stability2! <1.001 THEN‘bank shear constraint is satisifiedshearcond$ = “stable”bankcond$ = “stable”ELSE ‘must reduce theta!DO UNTIL shearcond$ = “just stable”theta? = (thetamax! + thetamin!) /2‘calculate mid pointIFtheta! <10*2*3.14159/360 THEN‘bank angles less than 10 degrees‘nominally assumed to be unstableshearcond$ = “unstable”bankcond$ = “unstable”EXIT DOEND IFCALL continuity‘satisfy continuitystability2! = bankshear! / taucrit!‘calculate stability criterion for bank shear constraintIF stability2! >= 1.001 THEN‘bankshear > taucritshearcondS = “unstable”thetamax! =theta!ELSEIF stability2! <= .999 THENshearcond$ = “understable”thetainin! =theta!ELSEshearcond$ = “just stable”bankcond$ = “stable”END IFIF (thetamax! - thetamin!) / theta! <.001 THEN281shearcond$ = “just stable”bankcond$ = “stable”stability2! = 1END IFLOOPEND IFEND SUBSUB bankstabilitynoncohesive‘calculates theta where banks just stablebankcond$ = “unknown”thetamax! = phiprime’thetamin! =0‘initialise search bounds and convergence criteriaDO UNTIL bankcond$ = ‘just stable”theta! = (thetamax’ ÷ thetamin!) /2‘determine midpoint of search rangelFtheta! <5 * 2 * 3.14159/360 THENbankcond$ = “unstable”EXIT DO‘if theta is less than 5 degrees the channel is assumed‘to be unstableEND IFCALL continuity‘satisfies continuity constraintstability! = (bankshear! / (gamma! * (ss! - 1) * d50bank!)) / (.048 * TAN(phiprime!) * (1- SIN(theta!) A 2/S1N(phiprime!) A 2) A .5)‘calculates stability criterionIF stability! >= 1.00 1 THENbankcond$ = “unstable”thetamax! =theta!ELSEIF stability! < .999 THENbankcond$ = “stable”thetamin! = theta!ELSEbankcond$ = “just stable”‘primary convergence criterionENDIFIF (thetamax! - thetainin!) I theta! <.001 THENbankcond$ = ‘just stable”stability! = 1282‘secondary convergence criterionEND IFLOOPEND SUBSUB bedload‘satisfies the bedload constraint using the Einstein Brown sedimenttransport relationgbcalc!=0minslope’ =0maxslope! = .05‘set bound of slope searchslopecond$ “unknown”‘initialise slope conditionDO WHILE slopecond$ <> “satisfied”slope? = (minslope! + maxslope!) I 2‘calculate mid point of range for bisectrix convergence schemeIF banktype$ = “noncohesive” THENCALL bankstabilitynoncohesive‘solve bank stability constraint for noncohesive banksELSE ‘banktype$ = “cohesive”CALL bankstabilitycohesive‘solve bank stability constraint for cohesive banksEND IFIF bankcond$ <> “unstable” THEN‘bank stability is satisfiedCALL ebrown‘calculate sediment transporting capacity of trial channelIF gbcalc! <.999 * gbload! THEN‘compares calculated to imposed sediment loadslopecond$ = “too flat”minslope! = slope!ELSEIF gbcalc!> 1.001 * gbload? THENslopecondS “too steep”maxslope! = slope!ELSEslopecond$ = “satisfied”‘primary convergence criterion‘sediment transporting capacity of the channel‘equals the imposed load bedload constraint satisfiedEND IFIF (maxslope! - minslope!) / slope! <.0001 THENIF gbcalc! > .99 * gbload! THENslopecondS = “satisfied”‘secondary convergence criterion283‘when bounds of slope very small: convergence attained‘occasionally a problem with convergence due to numerical‘approximationsELSEslopecond$ = “not satisfied”EXIT DO‘indicates that trial Pbed value is too narrowEND IFEND IFELSEmaxslope! = slope!‘trial slope too steep: reduce upper boundslopecond$ = “unstable”END IFLOOPEND SUBSUB bedload2‘identical to SUB bedload except search range for slope!‘reduced to save on computations. Optimum slope! for forward‘difference will be very close to that determined forbackward differenceminslope! = 9 * slope!maxslope! = 1.1 * slope!‘set bound of slope searchslopecond$ = “unknown”‘initialise slope conditionDO WHILE slopecond$ <> “satisfied”slope! = (minslope! + maxslope!) /2‘calculate mid point of range for bisectrix convergence schemeIF banktype$ = “noncohesive” THENCALL bankstabilitynoncohesive‘solve bank stability constraint for noncohesive banksELSE banktype$ = “cohesive”CALL bankstabilitycohesive‘solve bank stability constraint for cohesive banksEND IFIF bankcond$ <> “unstable” THEN‘ifbank stability is satisfiedCALL ebrown‘calculate sediment transporting capacity of trial channelIF gbcalc! <.999 * gbload! THEN‘compares calculated to imposed sediment loadslopecond$ = “too flat”minslope! = slope!ELSEIF gbcalc!> 1.001 * gbload! THEN284slopecondS = “too steep”maxslope? = slope!ELSEslopecond$ = “satisfied”‘sediment transporting capacity of the channel‘equals the imposed load : bedload constraint satisfiedEND IFIF (maxslope! - minslope!) / slope! <.0001 THEN slopecond$ = “satisfied”‘secondary convergence criterion‘when bounds of slope veiy small: convergence attained‘occasionally a problem with convergence due to numerical‘approximationsELSE‘if bank stability constraint not satisfiedmaxslope! = slope!‘trial slope too steep: reduce upper boundslopecond$ = “unstable”END IFLOOPEND SUBFUNCTION bedshear‘calculates the bed shear stressbedshear! = gamma! * depth! * slope! * (1 - shearforce!) * ((surfwidth! / (2 * pbed!) + .5))END FUNCTIONSUB continuity‘varies pbank! for trial values of Pbed!, slope!, and theta! to‘satisfy the contiuity constraintpbankcond$ = “unknown”errorcalci = 1000minpbank! 0maxpbank! =20 * discharge! “.35‘initialise search and convergence criteriaDO UNTIL pbankcond$ = “OK”pbank’ = (minpbank! + maxpbank!) /2‘calculate midpointerrorcalc! = (area * velocity / discharge!)‘calculate the normalised errorIF errorcalc! > 1.001 THENpbankcond$ = “too large”285maxpbank! = pbank!ELSEIF errorcaic! <.999 THENpbankcond$ = “too small”minpbank! = pbank!ELSEpbankcond$ = “OK”END IFIF (maxpbankl - minpbank’) / pbank? <.0001 AND pbankcond$ “too small” THEN‘resets maxpbank’ if too smallmaxpbank? =2* maxpbank!END IFLOOPEND SUBFUNCTION depth‘function to calculate flow depth of a trapezoidal channeldepth! = .5 * STN(theta!) * pbank’END FUNCTIONSUB ebrown‘calculate the sed trans capacity using em-brown formuladimlessbedshear! = bedshear! / (gamma! * (ss! - 1) * d50t)IF dimlessbedshear! <=0 THENdimlessbedload? = 0ELSEIF dimlessbedshear! <.093 THENdimlessbedload! = 2.15 * EXP(-.391 / dimlessbedshear’)ELSEdimlessbedload! = 40 * dimlessbedshear! “ 3END IFfi! = (2/3 + 36 * viscosity! A 2 / (gravity! * d50’ A 3 * (ss! - 1))) A 5 -(36 * viscosity! A 2 / (gravity! * dSO? A 3*(ss!_1)))gbcalc! =pbed! *fl! *! *deusity! *((ssi - 1)*gravity! *d5o! A3)A.5 *dje55dicadfEND SUBFUNCTION hydrad!‘function to calculate the hydraulic radiushydrad! area! / (pbed! + pbank!)END FUNCTIONFUNCTION meandepth‘function to calculate the mean depth286meandepth! = area / surfwidthEND FUNCTIONSUB optimum‘determines the optimal geometryoptcond$ = “unknown”‘initialises optimality test conditionlowerpbedl =0upperpbedl = 10 * discharge? A‘set bounds for Pbed based on Regime Eqns‘typical optimal value of Pbed is 2 to 5 * discharge Aminpbed! = lowerpbed?maxpbed! = upperpbedlPRINTPRINTDO WHILE optcond$ <> “optimum”pbed? = (minpbedl + niaxpbed!) /2‘calculate mid point of search rangeIFpbed’ > .95 * upperpbed! TI-lENupperpbedl =2 * upperpbed?maxpbed! = upperpbed?‘reset upper bound if necessaryEND IFPRiNT “Assessing Trial Bed Perimeter “;PRINT USING “#####“; pbed?;PRINT “ mpbed’ =pbed? * .975‘calculate backwards difference value of PbedCALL bedload‘satisfies continuity, bank stability and bedload‘constraints for the trail Pbed valueIF slopecond$ = “satisfied” AND bankcond$ <> “unstable” THEN‘test conditions indicate that stable channel geometry‘has been determined and the bedload constraint‘has been satisifiednetal! = gbload! / (density! * discharge! * slope!)‘evaluate neta for backward difference pointpbed! =pbed! / .975 * 1.025‘calculate pbed for forward duff287CALL bedload2‘satisfies continuity, bank stability and bedload‘constraints for the trail Pbed value‘indentical to Sub Bedload except the bounds of‘the slope search has been reduced to reduce computationsneta2 = gbload! / (density! * discharge! * slope!)‘calculates neta for forward difference valuepbed’ pbed! / 1.025‘reset pbed to midpoint valuednetabydpbed! = (neta2! - netal!) / (.05 * pbed!)‘calculate first derivative by finite difference‘to assess optimum conditionIF dnetabydpbed! <0 THEN ‘Trial Pbed is too widemaxpbed! = pbed! ‘Reduce upper bound of Pbedoptcond$ = “too wide”ELSETF dnetabydpbed! >0 THEN ‘Thai Pbed is too narrowminpbedt= pbed! ‘Increase lower bound of Pbedoptcond$ = “too narrow”ELSEoptcond$ = “optimum” ‘Optimum AchievedEND IFIF (maxpbed! - minpbed!) / pbed! <.001 THEN optcond$ = “optimum”‘Convergence attained. Second optimality criterionELSE ‘tnal Pbed too small for stable geomeuy: Increase minimum Pbedpbed! = pbed! I .975‘reset from backward difference to correct valueminpbed! =pbed!‘set lower limit at current trial value‘(as the optimum value must be greater)optcond$ = “too narrow”END IFPRINT optcond$LOOP‘Evaluate geometiy at exact optimal Pbed‘(not at forward or backward difference value)IF banktype$ = “noncohesive” THENCALL bankstabilitynoncohesive‘solve bank stability constraint for noncohesive banksELSE ‘banktype$ = “cohesive”CALL bankstabilitycohesive‘solve bank stability constraint for cohesive banksEND IF288CALL bedload2neta = gbload! / (density! * discharge! * slope!)‘evaluate objective functionEND SUBFUNCTION shearforce‘calculates the proportion of shear force on the banksshearforce! = 1.766 * (pbed! / pbank! + 1.5) A -1.4026END FUNCTIONSUB stabcurve‘stability curves are from Figure 5.3thetadegrees! = theta! / 2 / 3.14159 * 360iF phidegrees! <20 THEN‘use phiprime = 15 degree stability curve from Figure 5.3‘from Taylor (1948)SELECT CASE thetadegrees!CASE 48.65 TO 90stabnum! = -.141 * thetadegrees! + 17.46CASE 29.36 TO 48.6499stabnum! = -.533 * thetadegrees! + 36.53CASE 20.12 TO 29.3599stabnum! = -3.07 * thetadegrees! + 111CASE 15.001 TO 20.11999stabnum! = -189.6 * thetadegrees! + 3844CASE IS <= 15stabnum! = 1000000END SELECTELSE‘use phiprime =25 degrees curve from Figure 5.3After Taylor, (1948)SELECT CASE thetadegrees!CASE 58.62 TO 90stabnum! = -.212 * thetadegrees! + 24.86CASE 39.11 TO 58.61999stabnum? = -.8975 * thetadegrees! + 65.04CASE 25 TO 39. 10999289stabnum = -6.97 * thetadegrees! + 302.5CASE IS <=25stabnum? = 1000000END SELECTEND IFEND SUBFUNCTION surfwidthTunction to calculate the surface widthsurfwidth? = pbed? + COS(theta’) * pbankEND FUNCTIONFUNCTION velocity!‘calculates the mean velocityrhbank! = bankshear! / (gamma! * slope!)rhbed! = bedshear! / (gamma! * slope!)Ibank! = (2.03 * LOG(12.2 * rhbank! / ksbank!) / LOG(10)) ‘ -2ibed! = (2.03 * LOG(12.2 * rhbed! / ksbed!) / LOG(10)) A -2ifactor! = (ibed! * pbed! / (pbed! + pbank!) + fipJç! * pbank! / (pbank! + pbed!)) Avelocity! = (8 * gravity! * hydrad! * slope!) A •5 * ifactor!END FUNCTION290APPENDIX CFULL MODEL FORMULATIONThe Full Model Formulation differs from the Bankfull:Variable Slope Model in that the fillflow duration data is used as input, the sediment transporting capacity of the channel iscalculated over one year using the modified Parker (1990) surface-based sediment transportrelation, and the size distribution of the armour layer is treated as a dependent variable.The source code for the computer program (rivmod6.bas) used in the thesis is presented below.The programming was encoded and run using Microsoft Quick Basic (Version 3.0 or later).This program can be used for channels with cohesive or noncohesive bank sediment.The program is designed to use an input data ifie named rivmod.dat, and will write the outputdata to a data file named rivmod.out. The input file contains data for one channel only and mustbe set-up as follows (for noncohesive bank sediment):291m n Qbf Gb k3d k$bk banktype1.0Q Probability of Q1U Being Equalled or ExceededQ Probability of Q Being Equalled or ExceededQ Probability of Q Being Equalled or ExceededQ. Probability of Q’ Being Equalled or Exceeded0D 0D1 Proportion Finer than DD Proportion Finer than DD Proportion Finer than D’D Proportion Finer than DD 1.0D5Ob,For noncohesive channel banks the banktype is equal to “n”, and for cohesive channels it isequal to “c”.For cohesive banks the last line of the data file becomes:Ti C’ critIn rivmod.dat m is the number of flow intervals in the numerical approximation of a flowduration curve (See Figure 7.1); n is the number of intervals in the numerical approximation ofa sediment gradation curve (See Figure 7.2); Q is the discharge at the lower bound of flowinterval 1 which is the lowest flow on record (or an approximate value), Q, is the discharge at292the upper bound of flow interval m and therefore the maximum flow on record (or anapproximation); D is equal to D0, and D is equal to D1. In both cases the superscript 1 or urefer to the lower or upper bound of the interval, and the subscript indicates the intervalnumber. The upper bound of one interval is equal to the lower bound of the next interval. Thedata for the input file is designed to be read directly from flow duration and sediment gradationcurves. The geometric mean values for each interval are calculated in the program.The other symbols used above have been defined numerous times in the text and are given inthe List of Symbols.The data values must be separated by at least one space although the number of spaces is notimportant. An example of an input file is given below:11 7 85 2500000 0.4 0.4 ii3 122 0.247342 0.106455 0.056965 0.036575 0.024785 0,017290 0.01505205 0.0095115 0.007125 0.0054205 00.002 00.0064 0.10.01 0.20.018 0.40.032 0.60.06 0.80.091 0.90.15 10.07 40The output data will be output to rivmod. out in the following order:293W H 17* S D50 9 T*D50The output data from the sample input data file given above is:W H * S D50 theta T*D5037.7 1.29 1.16 0.00427 0.073 19.1 0.043The source code for the Full Optimisation Model is presented below.‘OPTIMISATION MODEL BASED ON FULLY DEVELOPED MODELDESCRIBED IN PhD THESIS BY RG MILLAR‘UNIVERSITY OF BRiTISH COLUMBIA‘SEPTEMBER 1994‘CALCULATES THE OPTIMUM (OR EQUILIBRIUM) CHANNEL GEOMETRYDECLARE SUB stabcurve 0DECLARE SUB bedload2 0DECLARE SUB bankstabilitynoncohesive 0DECLARE SUB bankstabiitycohesive 0DECLARE SUB parkerl99O 0DECLARE SUB capG 0DECLARE SUB vaiytheta 0DECLARE SUB optimum 0DECLARE SUB bedload 0DECLARE FUNCTION bankshear! 0DECLARE FUNCTION bedshear! 0DECLARE FUNCTION shearforce! 0DECLARE SUB continuity (discharge!)DECLARE FUNCTION velocity! 0DECLARE FUNCTION hydrad! 0DECLARE FUNCTION meandepth! 0DECLARE FUNCTION area! 0DECLARE FUNCTION surfwidth! 0DECLARE FUNCTION depth! 0program to calculate curves to demonstrate the optimumCOMMON SHARED gravity!, gamma!, ss!, viscosity!, density!, meandischargelCOMMON SHARED slope!, pbed!, pbank!, theta!, ddSO!, d50bank!COMMON SHARED qbfl, gbcalc!, gbload!, phiprime!, stability!, phidegreesCOMMON SHARED thetacond$, xnincond$, alpha!, dSO!, bankcond$, pi!COMMON SHARED ksbed!, ksbank!, phisgo!, sigmaphi!, sigxnaphitiial!, dsgttial!COMMON SHARED capgamma!, capitalg!, dsg!, i%, j%, m%, n%, iefi%COMMON SHARED banktype$, slopecond$‘Global Variables‘gravlty!= gravitational constant gamma! unit weight of water‘ss!= specific gravity of sediment viscosity!=kinematic viscosity of water‘density!=density of water slope!channel slope‘pbed!=bed perimeter pbank!bank perimeter294‘theta!=bank angle (radians) ddSO!=inedian armour grain diameter‘dSO!—median subarmour grain diameter d5Obank!median grain diameter of bank sediment‘discharge!=bankfufl discharge gbload’=imposed bedload transport rate‘gbload!=calculated bedload transport rate ‘neta’=coefficient of efficiency‘phiprime!=friction angle of bank sediment (radians)‘phidegrees!=friction angle ofbank sediment (degrees)‘ksbed!measure ofbed roughness ksbank!=measure ofbank roughness‘banktype$=type of bank sediment Wo=interger counter‘gammat!=unit weight of cohesive bank sediment‘stabnum!=stability number for cohesive sediment‘cohesion!=soil cohesion taucrit!=critical shear stress‘meandischarge!=mean annual discharge‘phisgol=dimensionless bed shear stress (Eqn Eqn 4.31)‘sigmaphi!=grainsize dispersion of armour (Eqn 4.37)‘sigmaphitrial!=trial value of sigmaphil‘dsg!=geometric mean grain diameter of pavement layer (Eqn 4.35)‘capgamtna! = constant converts kg/sec to kg/year‘capitalgl=dimensionless bedload function Eqn (4.32)‘j%ith discharge j%jth sediment size‘m%=numder of flow intervals n%=number of sediment size intervals‘SET VALUES OF INDEPENDENT CONSTANTSpi! = 3.14159gravity! = 9.81gamma! = 9810viscosity! = .000001density! = 1000ss! =2.65capganuna! = 365.25 * 24 * 3600CLSPRINTPRINTPRINTPRINT” RIVERMOD AN OPTIMISATION MODEL FOR”PRINT” PREDICTING THE ADJUSTMENT OF ALLUVIAL RIVERS”PRINTPRINTPRINTPRJNT” WR1TEN BY”PRINT” ROBERT MILLAR”PRINTPRINT” UNIVERSITY OF BRITISH COLUMBIA”PRINTPRINTPRINTPRINTPRINT” HiT TO CONTINUE”SLEEPINPUT dummyOPEN “nvmod.dat” FOR INPUT AS #1295OPEN “rivmod.out” FOR OUTPUT AS #2iNPUT #1, m%, n%, qbfl, gbload!, ksbed!, ksbanld, banktype$‘input number of flow and sediment size intervals respectively‘bankfufl discharge, mean annual sediment loadbed and bank roughness heights‘and bank typeIF bankt,e$ = “c” OR banktype$ = “C” THENbanktype$ = “cohesive”ELSEIF banktype$ = “N” OR banktype$ = “n” THENbanktype$ = “noncohesive”ELSEBEEPCLSPRINTPRINTPRiNTPRiNT “ERROR iN DATA FILE”ENDEND IF‘*4c************************************4*******************************************DIM SHARED sediment(n%), f(n%), ff(n%), discharge(m%), qprob(m%), upperd(n%)DIM SHARED propfiner(n%)DIM SHARED qb(m%, n%), qbi(m%), upperq(m%), fload(n%), bedsheari(m%)DIM SHARED probexceed(m%), depthi(m%)‘declare shared array variablesi% =0DOWHILEi% TO CONTINUE”SLEEPINPUT dununy$CLSPRINT “DATA SUMMARY”PRINTPRINTPRINT” Sediment Data”PRINT”j Di fl”‘calculate geometric mean grain diameter for each grain interval‘and volumetric proportion within each intervalj% = 1f(O) = 0DO WHTLEj% <= n%sediment(j%) = (upperd(j%) * upperd(j% - 1)) Af(j%) = propfiner(j%) - propfiner(j% - 1)PRINT USING “#41 #.###41 #.###“; j%; sediment(j%); f(j%)j%=j%+ 1LOOPPRINT‘calculate mean annual dischargei% = 1totalq!=ODOWHILEi%<=m%totalq! = totalq! + discharge’(i%) * ciprob(i%)i% = i% + 1LOOPmeandischarge! = totalq!PRINT USING ‘Mean Annual Discharge ##*#.# m”3/sec”; meandischargePRINT USING “Mean Annual Sediment Load II 111111111111 .41 kg/year”; gbload’‘input bank stability parametersIF banktype$ = “cohesive” OR banktype$ = “COHESIVE” THENbanktype$ = “cohesive”ELSEIF banktype$ = “NONCOHESIVE” OR banktype$ = “noncohesive” THENbanktype$ = “noncohesive”ELSEPRINT “Banktype unknoii: Error in Data File”BEEPCLOSE #1297CLOSE #2ENDEND IFPRINTPRINT “Bank Sediment Type:IF banktype$ = “noncohesive” THEN ‘input bank stability parametersINPUT #1, d5Obank!, phidegrees!PRINT USING “D50 bankPRINT USING “Phi#.### rn”; d50bank!##.# degrees”; phidegrees!ELSE ‘banktype$=”cohesive”INPUT #1, gzmmt!, phidegrees!, cohesion!, taucrit!PRiNT USING “ganunatPRINT USING “PhiPRINT USING “cohesionPRiNT USING “Tau critEND IF##.## kN/m”2”; gammat! / 1000#W.# degrees”; phidegrees!#### J/jA”; cohesion! / 1000##.# N/m’2”; taucrit!phiprime! = phidegrees! * 2 * pi! I 360PRINTPRINT : INPUT “IS INPUT DATA OK (y/n) “; booleanSIF boolean$ = “n” OR boolean$ = “N” THENCLOSE #1CLOSE #2ENDEND IFCLSPRINTPRINTdsgtrial! = (upperd(0) + upperd(n%)) /2‘set trial dsg value at midpoint of sediment gradationsigmaphitrial! = 1‘set trial sigamphi value‘COMMENCE OPTIMIZATION ROUTINECALL optimumPRINTPRINT” W H * S D50 theta T*D50”PRINT USING “###.# #4.## ##.#4 #.IIihf liii #.### ##.# #.###“; surfwidth!; depth!; meandepth!; slope!;dd5O!; theta! * 360 / 2 / pi!; bedshear! / (gamma! * dd5O! * (! - 1))PRINT #2,” W H Y* S D50 theta T*D50”298PRiNT #2, USING “##4t.# ##.## ##.## #.ihIIIiII( #.### ##.# #.###“; surfwidth?; depth!; meandepth!; slope!;ddSO!; theta! * 360 / 2 /pi!; bedshear / (gamma! * dd5O! * (ss! - 1))‘prints out the values of the dependent variablesCLOSE #1CLOSE #2END,**********************************************************************************FUNCTION area‘function to calculate the cross-sectional areaarea!=.5*(pbed!+surfwidth!)*depth!END FUNCTIONFUNCTION bankshear‘calculates the bank shear stressbankshear! = gamma! * depth! * slope! * shearforce! * ((surfwidth! + pbed!) * S1N(theta!) / (4 * depth!))END FUNCTIONSUB bankstabilitycohesive‘satisfies bank stability constraint for cohesive banksheightcond$ = “unknown”shearcond$ = “unknown”bankcond$ = “unknown”thetamax!=90*2*pi!/360thetamin! =0‘initialise search and convergence criteriaTirst test ifvertical bank is stable wrt bank heighttheta! = thetamax!CALL continuity(qbf!)‘satisfy continuityCALL stabcuive‘catculate stability numbercriticalheight! = stabnum! * cohesion! / gmmatI‘calculate the critical heightstabilityl! = depth / criticaiheight!‘calculates stability criteria wit bank height‘if stabilityl! <= criticaiheight then bank is stable‘with respect to bank heightIF stabilityl! <= 1.001 THENheightcondS = “stable”299ELSE ‘vertical bank is not stable therefore reduce the bank angle.Determine the maximum bank angle that is just stable.‘This is theta max from Chapter 6.DO UNTIL heightcond$ = ‘just stable”theta! = (thetamax! + thetamin!) /2‘calculate midpoint of range for‘bisectrix convergence schemeCALL continuity(qbfl)‘satisfy continuity constraintCALL stabcurve‘calculate stability numbercriticaiheight! = stabnum! * cohesion! / ganimat!stability 1! = depth / criticalheight!‘calculates critical height and stabilty number‘for the bank height constraintIF stabilityl! > 1.001 THENheightcond$ = “unstable”thetamax! = theta!ELSEIF stabiityl! <= .999 THENheightcond$ = “understable”thetamin! = theta!ELSEheightcond$ = ‘just stable”‘Bank Height = Critical HeightEN])IFIF (thetainax! - thetainin!) / theta! <.001 THENheightcond$ = ‘just stable”stabilityl! = 1‘Secondary Convergence Criterion in case of‘convergence problems due to numerical schemeEND IFLOOPEND IF‘bank height constraint now satisifled‘now assess the bank shear constraintthetamin!=0thetaniax! = theta!‘thetaniax is reset to maximum bankangle which satisfied the‘bank height constraint abovestabilily2! = bankshear! / taucrit!‘satbility criterion for bank shear constraint300IF stability2! < 1.001 THEN‘bank shear constraint is satisifiedshearcond$ = “stable”bankcond$ = “stable”ELSE ‘must reduce theta!DO UNTIL shearcond$ = “just stable”theta! =(thetaniax! +thetaniiu!)/2‘calculate mid pointIF theta! <10 * 2 * 3.14159/360 THENbank angles less than 10 degrees‘nominally assumed to be unstableshearcond$ = “unstable”bankcond$ = “unstable”EXIT DOEND IFCALL continuity(qbfi)‘satisfy continuitystability2! = bankshear! / taucrit!‘calculate stability criterion for bank shear constraintIF stability2! >= 1.001 THEN‘bankshear > taucntshearcond$ = “unstable”thetainax! =theta!ELSEIF stabiity2! <= .999 THENshearcond$ “understable”thetaniin! = theta!ELSEshearcondS = “just stable”bankcond$ = “stable”END IFIF (thetamax! - thetaniin!) / theta! <.001 THENshearcond$ = “just stable”bankcond$ = “stable”stability2! =1END iFLOOPEND IFEND SUB301‘**********************************************************************************SUB bankstabiitynoncohesive‘calculates theta where banks just stablebankcond$ = “unknown”thetamax! = phiprime’thetaniin! =0‘initialise search bounds and convergence criteriaDO UNTIL bankcond$ = ‘just stable”theta =(thetamax! +thetaniin!)/2‘determine midpoint of search rangeIF theta! <5*2*3.14159/360 THENbankcond$ = “unstable”EXIT DO‘if theta is less than 5 degrees the channel is assumed‘to be unstableEND IFCALL continuity(qbfl)‘satisfies continuity constraint for qbfstability! = (bankshear! / (gamma! * (ss! - 1) * d50bank!)) / (.048 * TAN(phiprime’) * (1 - S1N(theta!) “2/SIN(phiprime!) A 2) .5)‘calculates stability criterionIF stability! >= 1.001 THENbankcond$ = “unstable”thetamax! =theta!ELSEEF stability! <= .999 THENbankcond$ = “stable”thetamin! = theta!ELSEbankcond$ = “just stable”‘pnmaiy convergence criterionEND IFIF (thetamax! - thetamin!) / theta! <.001 THENbankcond$ = “just stable”stability! = 1‘secondaiy convergence criterionENDIFLOOPEND SUBSUB bedload‘satisfies the bedload constraint using the Parker 1990302‘surface-based sediment transport relationgbcalc! =0minslope! =0maxslope! = .05‘set bound of slope searchslopecond$ = “unknown”‘initialise slope conditionDO WHILE slopecond$ “satisfied”slope! = (minslope! + maxslope!) /2‘calculate mid point of range for bisectrix convergence schemeIF banktype$ = “noncohesive” THENCALL bankstabilitynoncohesive‘solve bank stability constraint for noncohesive banksELSE ‘banktype$ = “cohesive”CALL bankstabiitycohesive‘solve bank stability constraint for cohesive banksEND IFIF bankcond$ “unstable” THEN‘bank stability is satisfiedCALL parker 1990‘calculate sediment transporting capacity of trial channelIF gbcalc! <.999 * gbload! THEN‘compares calculated to imposed sediment loadslopecond$ = “too flat”minslope! = slope!ELSEIF gbcalc! > 1.001 * gbload! THENslopecond$ = “too steep”maxslope! = slope!ELSEslopecond$ = “satisfied”‘primary convergence criterion‘sediment transporting capacity of the channel‘equals the imposed load : bedload constraint satisfiedEN]) IFIF (maxslope! - minslope!) / slope! <.00001 THENlFgbcalc! >.99*gbload! THENslopecond$ = “satisfied”‘secondary convergence criterion‘when bounds of slope very small: convergence attained‘occasionally a problem with convergence due to numerical‘approximationsELSEslopecond$ = “not satisfied”EXIT DO‘indicates that trial Pbed value is too narrowEND IF303END IFELSEmaxslope! = slope!‘trial slope too steep: reduce upper boundslopecond$ = “unstable”END IFLOOPEND SUBSUB bedload2‘identical to SUB bedload except search range for slope!‘reduced to save on computations. Optimum slope! for forward‘difference will be very close to that determined for‘backward differenceminslope! = •9 * slope!maxslope! = 1.1 * slope!‘set bound of slope searchslopecond$ = “unknown”‘initialise slope conditionDO WHILE slopecond$ <> “satisfied”slope! = (minslope! + maxslope!) /2‘calculate mid point of range for bisectrix convergence schemeIF banktype$ = “noncohesive” THENCALL bankstabiitynoncohesive‘solve bank stability constraint for noncohesive banksELSE banktype$ = “cohesive”CALL bankstabiitycohesive‘solve bank stability constraint for cohesive banksEND IFIF bankcond$ “unstable” THEN‘ifbank stability is satisfiedCALL parker 1990‘calculate sediment transporting capacity of trial channelIFgbcalc! <.999*gbload! THEN‘compares calculated to imposed sediment loadslopecond$ = Htoo flat”minslope! = slope!ELSEIF gbcalc! > 1.001 * gbload! THENslopecond$ = “too steep”maxslopel = slope!ELSEslopecondS = “satisfied”‘sediment transporting capacity of the channel‘equals the imposed load : bedload constraint satisfiedEND IF304IF (maxslope! - minslope!) / slope! <.00001 THEN slopecond$ = “satisfied”‘secondary convergence criterion‘when bounds of slope very small: convergence attained‘occasionally a problem with convergence due to numerical‘approximationsELSE‘ifbank stability constraint not satisfiedrnaxslope! = slope!‘trial slope too steep: reduce upper boundslopecond$ = “unstable”END iFLOOPEND SUB‘**********************************************************************************FUNCTION bedshear‘calculates the bed shear stressbedshear! = ganuna! * depth? * slope! * (1 - shearforce!) * ((surfwidth! / (2 * pbed!) + .5))END FUNCTIONSUB capG‘calculate straining functionSELECT CASE phisgo!‘calculate omegao! based on piecewise linearisation of‘function given in Parker (1990) (See figure 4.10 in thesis)CASE IS <= 1.033omegao! = 1.011CASE IS <= 5.44omegao! = 1.027 * phisgo! A -.483CASE IS > 5.44omegao! = .453END SELECTSELECT CASE phisgo!‘calculate sigmaphio based on piecewise linearisation of‘function given in Parker (1990) (See figure 4.10 in thesis)CASE IS < .985sigmaphio = .816CASE IS <= 2.19sigmaphio! = .395 * phisgo! + .426CASE IS <=7sigiuaphio’ = .044 * phisgo! + 1.194CASE IS>7signiaphio! = 1.501END SELECT305straining! = 1 + sigmaphil I sigmaphio! * (omegao! - 1)‘Eqn4.38‘calculate hiding functionhiding! = (sediment(j%) / dsg!) A -.095 1‘Eqn 4,34 in thesismodifiedphi! = phisgo! * hiding! * straining!‘modifiedphi! is the product inside the square brackets‘in Eqn 4.33SELECT CASE modifiedphi!‘calculate Capital G: dimensionless bedad function‘Eqn 4.32 in thesisCASE IS < 1capitaig! =modifiedphi! A 14.2CASE IS <= 1.59capitalg! = EXP(14.2 * (modifiedphi! - 1) - 9.28 * (modifiedphi! - 1) A 2)CASE IS> 1.59capitaig! =5474*(1853/mjfipbj!)A45END SELECTEND SUBSUB continuity (discharge!)‘varies pbank! for trial values of Pbed!, slope!, and theta! to‘satisfy the contiuity constraintpbankcond$ = “unknown”errorcalcl = 1000minpbank! =0maxpbank! = 20 * discharge! A 35‘initialise search and convergence criteriaDO UNTIL pbankcond$ = “OK”pbank! = (minpbank! + maxpbank!) /2‘calculate midpointerrorcalc! = (area * velocity / discharge!)‘calculate the normalised errorIF errorcaic! > 1.001 THENpbankcond$ = “too large”niaxpbank! = pbank!ELSEIF errorcalc! <.999 THENpbankcond$ = “too small”minpbank! = pbank!306ELSEpbankcond$ = “OK”END IFIF (maxpbank! - minpbank!) / pbank! <.0001 AND pbankcond$ = “too small” THEN‘resets maxpbánk! if too smallmaxpbank! =2* maxpbank!END IFLOOPEND SUBFUNCTION depth‘ftmction to calculate flow depth of a trapezoidal channeldepth! = .5 * S1N(theta!) * phank!END FUNCTIONFUNCTION hydrad!‘function to calculate the hydraulic radiushydrad! =area! /(pbed! +pbank!)END FUNCTIONFUNCTION meandepth‘function to calculate the mean depthmeandepth! = area / surfwidthEND FUNCTIONSUB optimum‘determines the optimal geometryoptcond$ = “unknown”‘initialises optimality test conditionlowerpbed! =0upperpbed! = 10 * qbf! “.5‘set bounds for Pbed based on Regime Eqns‘typical optimal value of Pbed is 2 to 5 * discharge” .5minpbed! = lowerpbed!maxpbed! = upperpbed!PRINTPRINTDO WHILE optcond$ <> “optimum”307pbed! = (minpbed! + maxpbed!) /2‘calculate mid point of search rangePRINT “Assessing Trial Bed Perimeter “;PRINT USING “#### ##“; pbedPRINT “m “;pbed! = pbed! * .975‘calculate backwards difference value of PbedIFp > 95 * upperpbed! THENupperpbed! =2 * upperpbedlmaxpbed’ = upperpbed!‘reset upper bound if necessalyEND IFCALL bedload‘satisfies continuity, bank stability and bedload‘constraints for the trail Pbed valueIF slopecond$ = “satisfied” AND bankcond$ <> “unstable” THEN‘test conditions indicate that stable channel geometty‘has been determined and the bedload constraint‘has been satisifiednetal! = gbload! I (capgamma! * density! * meandischarge! * slope!)‘evaluate nets for backward difference pointpbed! =pbed! 1.975 * 1.025‘calculate pbed for forward duffCALL bedload2‘satisfies continuity, bank stability and bedload‘constraints for the trail Pbed value‘indentical to Sub Bedload except the bounds of‘the slope search has been reduced to reduce computationsneta2! = gbload! / (capgainma! * density! * meandischarge! * slope!)‘calculates nets for forward difference valuepbed! = pbed! / 1.025‘reset pbed to midpoint valuednetabydpbed! = (neta2!- netal!) / (.05 * phed!)‘calculate first derivative by finite difference‘to assess optimum conditionIF dnetabydpbed! <0 THEN ‘Trial Pbed is too widemaxpbed! =pbedl ‘ReduceupperboundofPbedoptcond$ “too wide”ELSEIF dnetabydpbed! >0 THEN ‘Trial Pbed is too narrowminpbed! = pbed! ‘Increase lower bound of Pbedoptcond$ = “too narrow”ELSE308optcond$ = “optimum” ‘Optimum AchievedEND IFIF (maxpbed - nunpbedl) / pbed’ <.001 THEN optcond$ = “optimum”‘Convergence attained. Second optimality criterionELSEtrial Pbed too small for stable geometry: Increase minimum Pbedpbed! =pbed’ 1.975‘reset from backward difference to correct valueminpbed! = pbed!‘set lower limit at current trial value‘(as the optimum value must be greater)optcond$ = “too narrow”END IFPRINT optcond$LOOP‘Evaluate geometry at exact optimal Pbed‘(not at forward or backward difference value)IF banktype$ = “noncohesive” THENCALL bankstabilitynoncohesive‘solve bank stability constraint for noncohesive banksELSE ‘banktype$ = “cohesive”CALL bankstabilitycohesive‘solve bank stability constraint for cohesive banksEND IFCALL bedload2‘satisfies continuity, bank stability and bedload‘constraints for the midpoint Pbed value‘indentical to Sub Bedload except the bounds of‘the slope search has been reduced to reduce computationsEND SUB‘*********************************************************************************,SUB parkerl99ODIM OMEGAIJ(m% n%), OMEGAJ(n%)‘OMEGAJ() and OmegalJ 0 are used in the snmmations in‘Eqn 4.53 onwards. See Below for explanation‘these are key functions in the modified Parker (1990)‘bedload equation= 1DOWH[LEi%<=m%discharge! = discharge(i%)IF discharge! > qbfl THEN discharge! = qbfl309‘sets all flows in excess of qbf equal to qbfl‘in order to calculate sediment transport rate‘for overbank flowsCALL continuity(dischargel)‘satisify continuity for flow i%bedsheari(i%) = bedshear!depthi(i%) = depth!‘calculate depth of flow and bedshear for each flow‘this saves recalculation below‘records depths and bedshear values in an array= i% + 1LOOPdsgcond$ = “unknown”sigmaphicondS = “unknown”‘initialise convergence criteria for dsg! and sigmaphi!DO WHILE dsgcond$ “converged” OR sigmaphicondS “converged”dsg!=dsgtrial!sigmaphi! = sigmaphitrial!OMEGA! =0j% = 1DOWIIILEj%<=n%= 1OMEGAJj%) =0DOWHILEi%<=m%‘calculate for sediment(j%) over range of flowsphisgo! = bedsheari(i%) / (gamma! * (ss! - 1) * dsg!) / .0386CALL capGOMEGAIJ(i%, j%) = qprob(i%) * capitaig! * bedsheari(i%) A 1.5‘OMEGAIJ(m%) is used in the summation over i% in the numerator‘of Eqn 4.52OMEGAJ(j%) = OMEGAJ(j%) + OMEGAIJ(i%, j%)‘summation for all flows for sediment(j%) Eqn 4.52‘OMEGAJ(n%) is the numerator in Eqn 4.52= i% + 1LOOPOMEGA! = OMEGA! + f(j%) I OMEGAJ(j%)‘OMEGA is the denomenator of Eqn 4.52j%=j%+ 1LOOP‘calculate Fj and Dsg and sigmaphi!j% = 1lndsg! =0310sigmaphi2! =0ff1 =0DOWHILEj%<=n%ff(j%) = f(j%) / OMEGAJ(j%) / OMEGA!Indsg! = lndsg! + ff(j%) * LOG(sediment(j%))sigmaphi2! = sigiuaphi2! + (LOG(sediment(j%) / dsgl) /LOG(2)) “2 * ff(j%)j%=j%+ 1LOOPdsg! = EXP(lndsg!)sigmaphi! = sigmaphi2! “.5‘test for convergence for dsg!IF (dsg! / dsgtrial!) < .99 THENdsgcond$ = “too small”ELSEIF (dsg! I dsgtrial!)> 1.01 THENdsgcond$ = “too large”ELSEdsgcond$ = “converged”END IF‘test for convergence for siginaphi!IF (sigmaphi! / sigmaphitrial!) < .99 THENsigmaphicond$ = “too small”ELSEIF (signiaphi! / sigmaphitrial!)> 1.01 THENsigmaphicondS = “too large”ELSEsigmaphicond$ = “converged”END IFsigmaphitrial! = sigmaphidsgtnal! dsg!LOOP‘loop until convergence for dsg! and signaphi!‘calculate the median grain diameterff(0) = 0sediment(0) =0j% =0totif! =0DO WHILE totifi <.5j%=j%+ 1totifi = totfP + ff(j%)LOOPdeldd50! = (totifi - .5) * (upperd(j%) - upperd(j% - 1)) / ff(j%)ddSO! = upperd(j%) - deldd50!‘calculate D50 by interpolationqbcalc’ = .00218 I ((ss! - 1) * gravity! * densityV 1.5 * OMEGA!)‘calculate the mean unit volumetric sediment transport rate‘in units of m”2/sec311gbcalc! = qbcalc? * pbed! * density! * ss! * capganuna!‘cacluate the total load in mass units‘transported over one year.END SUBFUNCTION shearforce‘calculates the proportion of shear force on the banksshearforce! = 1.766 * (pbed! /pbank! + 1.5) A -1.4026END FUNCTION,**********************************************************************************SUB stabcurve‘stability curves are from Figure 5.3‘NB ‘!!!? THE STABILITY CURVES USED IN ThIS SUB ARE‘FROM TAYLOR (1948) AND APPLY TO HOMOGENEOUS FULLY DRAINED‘SOIL:- THESE ARE FOR ILLUSTRATIVE PURPOSES ONLY!!!thetadegrees! = theta! / 2/3.14159 * 360IF phidegrees! <20 THEN‘use phiprime = 15 degree stability curve from Figure 5.3‘from Taylor (1948)SELECT CASE thetadegrees!CASE 48.65 TO 90stabnum! = -.141 * thetadegrees? + 17.46CASE 29.36 TO 48.6499stabnum! = -.533 * thetadegrees? + 36.53CASE 20.12 TO 29.3599stabnum! = -3.07 * thetadegrees! + 111CASE 15.001 TO 20.11999stabnum! = -189.6 * thetadegrees! + 3844CASE IS <=15stabnum! = 1000000‘bank of infinite height is stableEND SELECTELSE‘use phiprinie =25 degrees curve from Figure 5.3‘After Taylor, (1948)SELECT CASE thetadegrees!CASE 58.62 TO 90stabnum! = -.2 12 * thetadegrees! + 24.86312CASE 39.11 TO 58.61999stabnum! = -.8975 * thetadegrees! + 65.04CASE 25 TO 39.10999stabnum! = -6.97 * thetadegrees! + 302.5CASE IS <=25stabnuml = 1000000bank of infinite height is stableEND SELECTEND IFEND SUBFUNCTION surfwidth‘function to calculate the surface widthsurfwidth! = pbed’ + COS(theta!) * pbank!END FUNCTIONFUNCTION velocity!‘calculates the mean velocityrhbank! = bankshear! / (gamma! * slope?)rhbed? = bedshear! / (gamma! * slope!)Ibank! = (2.03 * LOG(12.2 * rhbank! / ksbank!) / LOG(10)) “-2Ibed! = (2.03 * LOG(12.2 * rhbed! / ksbed!) / LOG(10)) A -2ffactor2! = (Ibed! * pbed! / (pbed! + pbank!) + thank! * pbank! / (pbank! + pbed!)) Avelocity! = (8 * gravity! * hydrad! * slope!) A 5 * fl’actor2!END FUNCTION313APPENDIX DDATA FROM HEY AND THORNE (1986)AND ANDREWS (1984)Table D-1 contains the data from Hey and Thorne (1986). The “Observed” columns contain theactual published data: Qbfi W, Y, D50 and Vegetation Type. The values of k3 and Gbf arecalculated from the observed hydraulic geometry using Eqn (3.2) and the Einstein-BrownEquation (Eqn 6.2) respectively. The value Of is set equal to D50, and d50 is set equal toD50/3.Table D-2 contains the data from Andrews (1984). The “Observed” columns contain the actualpublished data: QbJ; W, Y, d50 and Vegetation Type. The values of k and Gbf arecalculated from the observed hydraulic geometry using Eqn (3.2) and the Einstein-BrownEquation (Eqn 6.2) respectively. The value of is set equal to 1)5o. Channel Number9074800 from Andrews (1984) was excluded from the analysis. (See Section 6.3.3).314TableD-1. DatafromHeyandThorne(1986).CalculatedChannelObservedFixedSVariableSNumber,=40°Est.Meanççb,—40°Qbfk3GbfWYSD50d50D5obOVegWIWIWIS(mi/s)(m)(kg/s)(m)(m)(m)(m)(m)Type(m)(m)(°)(m)(m)(m)(m)11040.16101.5525.21.520.00350.0570.0190.057III41.11.1555.527.11.4845.11.040.004121540.1194.4342.71.400.00310.0610.0200.061II45.41.3541.9361.5646.21.320.003231240.2640.9632.71.700.00250.0820.0270.082III33.91.6941.325.12.0434.21.660.00264170.604.9712.80.910.00450.1090.0360.109II11.60.9935.9101.0811.410.00455630.20293.8620.01.270.00420.0240.0080.024III71.70.6272.544.60.8189.80.480.00646680.0713.8321.61.360.00220.0610.0200.061IV21.41.3839.715.91.722.51.320.00237701.34450.7527.51.150.01270.1760.0590.176I37.60.9948.729.21.1339.40.950.013781480.0424.2929.41.940.00140.0410.0140.041IV381.6849.523.52.3139.21.60.0015976.10.24430.2415.81.470.00640.0460.0150.046IV59.70.6975.828.21.0772.60.560.0092101720.0722.7528.72.190.00160,0700.0230.070IV312.1246.323.52.5831.82,040.0017113040.1135.1248.52.490.00130.0660.0220,066III47.22.5837.536.13.0847.52.570.001312450.5427.5918.21.410.00300.0670.0220.067III25.61.185518.41.4426.41.140.003213500.851341.018.41.140.01330.0600.0200.060IV74.40,527333.20.7991.30.440.018414460.0719.9620.11.050.00270.0430.0140.043II240.964619.51.0924.60.920.002915530.1930.1725.21.020.00390.0750.0250.075II23.81.073819.81.223.91.060.003916530,63100.3224.41.160.00510.0710.0240,0711134.20.9750.126.81.12360.920.005617690.1061.0829.71.020.00360.0450.0150.045II35.90,9245.228.11.0637.40.870.003918480.12170.8713.01.070.00820.0650.0220.065IV32.20.646716.20.96370.550.010519240.2121.2417.70.800.00400.0480.0160.048II20.30.7544.116.20.8520.90.720.004220240.6641.8921.10.920.00340.0340.0110.034131.60.7551.229.10.7833.40.710.003821220.5330.6510.21.420.00220.0190.0060.019IV370.79018.21.0540.60.640.002622140.1450.4610.60.750.00390.0140.0050.014IV37.40.3772.517.60.5746.30.290.0058237.10.4613.0213.70.480.00610.0480.0160.048I13.70.4939.712.80.5113.70.480.006224100.2190.595.50.590.02140.1090.0360.109IV13.40.3673.97.20.5214.40.330.024925280.3433.0017.00.990.00360.0430.0140.043II25.40852,8200.922670.750.0040TableD- Icontinued.CalculatedChannelObservedFixedSVariableSNumber=40°Est.Meanç6,Qbfk,GbfWYSD50d.,0D50VegW1’WYWIS(mi/s)(m)(kg/s)(m)(m)(m)(m)(m)Type(m)(m)(°)(m)(m)(m)(m)26270.3240.8915.00.950.00480.0550.0180.055II23.10.7554.518.40.8624.30.710.005427340.328.3814.41.610.00130.0270.0090.027III26.31.1570.318.41.4426.71.120.0014281000.2161.9725.91.520.00350.0810.0270.081IV32.21.3547.319.71.8533.21.290.0038291000.1731.4023.01.700.00270.0840.0280.084IV27.11.5749.519.51.9628.91.490.0029307.50.456.9215.60.480.00520.0550.0180.055112.10.5630.911.60.5711.70.580.0048317.50.300.0712.60.560.00340.0860.0290.08617.10.821.16.90.8270.860.0027321700.17116.4231.81.820.00330.0790.0260.079IV43.11.5350.125.12.1645.41.440.0037337.10.211.919.60.750.00150.0170.0060.017II13.70.625511.30.714.10.60.0016349.70.3239.8312.30.480.01080.0660.0220.066115.40.4346.314.30.45160.410.011535350.44671.8212.50.930.01520.0680.0230.068IV46.70.4673.621.50.756.60.380.02123611.21.5122.7114.10.810.00550.0740.0250.0741180.7349.817.20.7518.10.720.005737530.38302.5815.21.180.00920.0770.0260.077IV39.30.767.719.11.0645.30.610.0117382371.1279.1357.82.600.00140.0560.0190.056I65.52.4544.161.62.5466.12.410.0015394241.47316.9176.53.100.00160.0500.0170.050I117.92.4651.71092.57126.82.290.001840610.1023.2922.91.200.00270.0530.0180.053IV261.1344.116.11.5426.41.10.002941610.3690.4827.71.290.00260.0280.0090.028II52.10.9157.740.21.0558.70.80.003242450.29107.%18.61.120.00460.0420.0140.042II38.80.746129.90.86440.650.005743750.5087.7320.41.580.00420.0810.0270.081IV33.51.2158.319.31.735.41.140.0047441200.59348.1041.71.250.00660.0910.0300.091I53.51.146.849.51.1456.41.040.0072453480.38438.0746.62.510.00300.0660.0220.066IV87.21.7658.343.82.6697.61.560.00374660.40.133.703221.310.00110.0360.0120.036I23.51.6130.422.31.6622.81.690.00104736.51.28126.7228.00.990.00570.0640.0210.064I380.8648.435.10.8939.70.820.00614874.70.87666.5819.31.530.00730.0580.0190.058IV69.40.7873.631.81.1882.60.660.009849950.4229.1332.21.400.00380.1300.0430.130II24.51.6631.522.31.7723.71.740.003550500.3498.1317.51.200.00550.0700.0230.070III30.40.8858.320.11.1334.20790.0064TableD-1continued.CalculatedChannelObservedFixedSVariableSNumberçó.=40°Est.Mean=40°Qbfk3Gk,WYSD50d50D50VegWY,..WYW‘(ms/s)(m)(kg/s)(m)(m)(m)(m)(m)Type(m)(m)(9(m)(m)(m)(m)5190.50.5732.6424.61.730.00330.1230.0410.123IV24.81.7639.719.82.0525.11.720.003452126.81.25227.4031.32.030.00360.0710.0240.071III59.11.4659.137.41.8863.91.340.004253358.30.1370.3177.12.010.00160.0600.0200.060I612.3233.457.42.4159.72.390.00155491.30.1026.7731.21.410.00190.0410.0140.041II34.91.3343.528.11.5235.61.290.002055153.30.099.5533.41.980.00140.0690.0230.069III28.62.2334.222.32.6428.42.260.001456192.30.6786.2145.22.080.00250.0920.0310.092111472.0640.833.62.5347.22.030.00265717.60.1218.189.10.730.00840.1140.0380.114IV10.70.6852.38.10.8211.20.650.0089581960.37334.7541.01.700.00470.0840.0280.084III60.51.3751.238.51.7965.11.270.005459190.2424.4111.61.050.00240.0180.0060.018111310.6170.319.40.7935.50.530.003060200.0523.2217.50.620.00350.0200.0070.020I27.40.4852.825.50.530.10.430.0041613.90.150.676.50.590.00190.0230.0080.023117.40.5646.36.40.627.60.540.002062280.052.2812.31.140.00180.0490.0160.049III131.1343010.31.3313.11.110.0019TableD-2.DatafromAndrews(1984).00CalculatedChannelObservedFixedSVariableSNumberçS.=40°Est.Mean—=400Qbfk3GbfWYSD.sod,0D50VegWY.WYWYS(mIs)(m)(kg!s)(m)(in)(m)(m)(m)Type(m)(m)(0)(m)(m)(in)(m)66119006.710.032.9010.40.480.00230.0230.0050.023Thin9.70.537.59.60.59.60.50.002466148000.70.177.852.30.200.02600.0490.0190.049Thick5.10.1367.630.185.40.130.0293662000085.20,0213.0147.20.970.00140.0240.0080.024Thin36.11.1633.1361.1634.71.220.001370830007.080.83129.829.10.520.01700.0790.0230.079Thick18.70.3760.510.50.4920.20.350.019070910009.770.45222.0111.50.460.01900.0850.0190.085Thick18.60.3652.810.50.49200.340.0213901350012.20.4012.1211.90.730.00460.0610.0190.061Thin140.6846.313.90.68140.670.004890180002.210.6210.617.30.340.01100.0520.0190.052Thin9.50.347.99.50.39.70.290.011590220002.690.2430.527.00.290.01400.0340.0110.034Thick130.2156.67.20.2914.40.190.016590348001.870.132.165.20.310.00610.0240.0110.024Thick7.30.2652.84.80.347.80.250.006790359008.360.18123.259.20.430.01500.0730.0130.073Thick13.40.3550.680.4714.40.330.0168903600022.60.1619.6117.50.730.00440.0430.0170.043Thick21.10.6645.7130.8922.10.630.004890781003.170.3410.106.30.390.01000.0610.0190.061Thick8.30.3450.65.50.448.40.330.010690782002.520.5010.885.60.410.00920.0370.0150.037Thick11.40.2963.86.70.38120.270.01029081600490.46169.6934.10.840.00580.0980.0150.098Thin25.21.0132.025.31.0124.91.030.0055911250037.50.8366.2226.00.910.00670.0910.0360.091Thin27.60.8941.927.50.928.10.880.006991155007.080.297.0611.60.520.00460.0460.0140.046Thin11.40.5339.711.40.5311.50.520.00479124500420.34151.6624.90.880.00580.0790.0140.079Thin25.40.8840.825.30.8825.70.860.006092425001010.64114.8336.61.450.00370.0640.0240.064Thin48.51.2547.947.91.26521.180.004092444101670.2146.7753.31.630.00180.0580.0170.058Thin43.41.8533.9431.8742.91.890.0017924950046.71.3916.8224.41.620.00200.0450.0220.045Thin32.61.4148.432.51.4135.41.30.002392510002550.1251.5383.81.850.00090.0340.0070.034Thin58.32.330.457.92.356.32.410.0008925300072.20.692184.7530.51.130.00710.1220.0070.122Thin32.41.141.932.31.1133.41.070.007692570001140.4342.6636.61.650.00240.0700.0250.070Thin37.51.6540.837.51.6538.21.610.0025APPENDIX EDATA FROM CHARLTON ETAL. (1978)Table E-1 contains the data from Chariton eta!. (1978) for channels with cohesive banks whichwas analysed in Chapter 6. The actual published data are: Q W, 7, S, D.50, q, and theVegetation Type. The values for k are calculated from the observed hydraulic geometry usingEqn (3.2), and c’ is calculated from qu using Eqn (6.4) with the assumption of ç4’ = 25° for allchannels. The value ofd50 is set equal toD50/3, and the drained and saturated values of the unitsoil weight are assumed to be 20.0 and 22.45 kN/m3respectively for all channels.The modelled output data is presented in Table 6.4.319t%)CTableE-1.DatafromCharitonela!. (1978).ReachRiverWYSD,d,oVegNumber(m’/s)(m)(m)(m)(in)(m)kN/m2kN/m2(°)kN/m3kN/m3Type1AforLwyd6417.41.780.00441.090.0750.02539.512.62522.4520T2Aiwen1-510140.730.00641.420.1060.03525.68.22522.4520G3Alwen5-910.79.80.730.01301.250.1130.03823.27.42522.4520T4Arrow29.513.71.340.00451.160.0410.01451.416.42522.4520T5Ceirog1-36617.61.790.00481.210.0710.02423.37.42522.4520T6Ceirog4-766191.360.01051.070.0820.02770.922.62522.4520T7Exe67311.770.00181.350.0430.01426.18.32522.4520G8IrfonC8128.71.630.00140.190.0550.01830.19.62522.4520T9Otter14.216.70.690.00320.290.0570.01925.88.22522.4520G10Teign66192.470.00141.410.0510.01715.75.02522.4520T11Trent2.75.20.650.00230.560.0330.01136.711.72522.4520G12Usk15739.32.640.00090.600.0720.02421.36.82522.4520G13Wye55059.44.190.00070.280.0280.00916.05.12522.4520G14Wrye3819.51.670.00201.540.0400.01339.912.72522.4520G