MORPHODYNAMICS OF SAND-BED RIVERS ENDING IN DELTAS. Wax Lake Delta, Louisiana. Delta from Iron Ore Mine into Lake Wabush, Labrador. Delta of Eau Claire River at Lake Altoona, Wisconsin. DELTA SHAPE.
Wax Lake Delta, Louisiana
Delta from Iron Ore Mine into Lake Wabush, Labrador
Delta of Eau Claire River at Lake Altoona, Wisconsin
Most deltas show a 2D spread in the lateral direction (fan-delta). The channel(s) migrate and avulse to deposit over the entire surface of the fan-delta as they prograde into standing water.
But some deltas forming in canyons can be approximated as 1D. The delta of the Colorado River at Lake Mead was confined within a canyon until recently.
Mangoky River, Malagasy
Colorado River at Lake Mead, USA
Hoover Dam was closed in 1936. Backwater from the dam created Lake Mead. Initially backwater extended well into the Grand Canyon. For much of the history of Lake Mead, the delta at the upstream end has been so confined by the canyon that is has propagated downstream as a 1D delta. As is seen in the image, the delta is now spreading laterally into Lake Mead, forming a 2D fan-delta.
View of the Colorado River at the upstream end of Lake Mead.
Image from NASA
Image based on an original from Smith et al. (1960)
Image based on an original from Smith et al. (1960)
A typical delta deposit can be divided into a topset, foreset and bottomset. The topset is coarse-grained (sand or sand and gravel), and is emplaced by fluvial deposition. The foreset is also coarse-grained, and is emplaced by avalanching. The bottomset is fine-grained (mud, e.g. silt and clay) and is emplaced by either plunging turbidity currents are rain from sediment-laden surface plumes.
Here the problem is simplified by considering a topset and foreset only. That is, the effect of the mud is ignored. It is not difficult to include mud: see Kostic and Parker (2003a,b).
The channel has a constant width. Let x = streamwise distance, t = time, qt = the volume sediment transport rate per unit width and p = bed porosity (fraction of bed volume that is pores rather than sediment). The mass sediment transport rate per unit width is then sqt, where s is the material density of sediment. Mass conservation within a control volume with length x and a unit width (Exner, 1920, 1925) requires that:
/t (sediment mass in bed of control volume) = mass sediment inflow rate – mass sediment outflow rate
BUT HOW ARE FLOODS ACCOUNTED FOR IN THE EXNER EQUATION OF SEDIMENT CONTINUITY?
fraction of time flow is in flood
SIMPLE ADAPTATION TO ACCOUNT FOR FLOODING
Rivers move the great majority of their sediment, and are morphodynamically active during floods. Paola et al. (1992) represent this in terms of a flood intermittency If. They characterize floods in terms of bankfull flow, which carries a volume bed material transport rate per unit width qt for whatever fraction of time is necessary is necessary to carry the mean annual load of the river. That is, where Gma is the mean annual bed material load of the river, Bbf is bankfull width and Ta is the time of one year, If is adjusted so that
where = water density. The river is assumed to be morphodynamically inactive at other times. Wright and Parker (2005a,b) offer a specific methodology to estimate If.
The Exner equation is thus modified to
KEY FEATURES OF THE MORPHODYNAMICS OF THE SELF-EVOLUTION OF 1D SAND-BED DELTAS UNDER THE INFLUENCE WITH BACKWATER
Water discharge per unit width qw is conserved, and is given by the relation
and shear stress is related to flow velocity using a Chezy (recall Cf-1/2 = C) or Manning-Strickler formulation;
In low-slope sand-bed streams boundary shear stress cannot be computed from from the depth-slope product, but instead must be obtained from the full backwater equation;
Base level is specified in terms of a prescribed downstream water surface elevation d(t) = (L, t) +H(L, t) rather than downstream bed elevation (L, t).
The base level of the Athabasca River, Canada is controlled by the water surface elevation of Lake Athabasca.
Delta of the Athabasca River at Lake Athabasca, Canada.
Image from NASA
The formulation below assumes a single characteristic bed grain size D, a constant Chezy resistance coefficient Cz = Cf-1/2 and the Engelund-Hansen (1967) total bed material load formulation as an example.
Bed material load equation
Specified initial bed profile
Specified upstream bed material feed rate
Downstream boundary condition, or base level set in terms of specified water surface elevation.
NUMERICAL SOLUTION TO THE BACKWATER FORMULATION OF MORPHODYNAMICS
The case of subcritical flow is considered here.
At any given time t the bed profile (x, t) is known. Solve the backwater equation
upstream from x = L over this bed subject to the boundary condition
Evaluate the Shields number and the bed material transport rate from the relations
Find the new bed at time t + t
Repeat using the bed at t + t
IN MORPHODYNAMICS THE FLOW AND THE BED TALK TO AND INTERACT WITH EACH OTHER
Sure, sweetie, but could you cut my toenails for me afterward? I can’t reach ‘em very well either.
THE PROBLEM OF IMPULSIVELY RAISED WATER SURFACE ELEVATION (BASE LEVEL) AT t = 0
M1 backwater curve
Note: the M1 backwater curve was introduced in the lecture on hydraulics and sediment transport
A PROGRADING DELTA THAT FILLS THE SPACE CREATED BY BACKWATER
See Hotchkiss and Parker (1991) for more details.
Before the water surface is raised, the is assume to be at normal, mobile-bed equilibrium with antecedent bed slope Sa and water dischage per unit width during floods qw. It is useful to compute the characteristics of the normal flow that would prevail with the specified flow over the specified bed. Let Hna and qtna denote the flow depth and total volume bed material load per unit width that prevail at antecedent normal mobile-bed equilibrium with flood water discharge per unit width qw and bed slope Sa. From Slides 19 and 21 of the lecture on hydraulics and sediment transport, these are given as
The reach has length L and the bed has antecedent bed slope Sa. It is assumed for simplicity that the antecedent bed elevation at the downsteam end of the reach (x = L) is da = 0, so that antecedent water surface elevation da is given as
At time t = 0 the water surface elevation (base level) at the downstream end d is impulsively raised to a value higher than da and maintained indefinitely at the new level. In addition, the sediment feed rate from the upstream end is changed from the antecedent value qta to a new feed value qtf.
Eventually the backwater zone (e.g. reservoir zone behind a dam) fills with sediment and the river established a new normal, mobile-bed equilibrium in consonance with the flood, water discharge per unit width qw, the sediment supply rate qtf (which becomes equal to the transport rate qtu at ultimate equilibrium) and the water surface elevation d. The bed slope Su and flow depth Hnu at this ultimate normal equilibrium can be determined by solving the two equations below for them.
A morphodynamic numerical model can then be used to describe the evolution
from the antecedent equilibrium to the ultimate equilibrium
The channel is assumed to have uniform grain size D and some constant ambient slope Sa (before changing conditions at t = 0) which is in equilibrium with an ambient transport rate qta. The reach of interest has length L. The antecedent bed profile (which serves as the initial condition for the calculation) is then
where here da can be set equal to zero. The boundary condition at the upstream end is the changed feed rate qtf for t > 0, i.e.
where qtf(t) is a specified function (but here taken as a constant). The downstream boundary condition, however, differs from that used in the normal flow calculation, and takes the form
where d(t) is in general a specified function, but is here taken to be a constant.
Note that downstream bed elevation (L,t) is not specified, and is free to vary during morphodynamic evolution.
The reach L is discretized into M intervals of length x bounded by M+1 nodes. In addition, sediment is fed in at a “ghost” node where bed elevation is not tracked.
The backwater calculation over a given bed proceeds as in the lecture on hydraulics and sediment transport:
where i proceeds downward from M to 1.
Note that the difference scheme used to compute qt,i/x is a central difference scheme only for au= 0.5. With a backwater formulation takes an advectional-diffusional form and a value of au greater than 0.5 (upwinding) is necessary for numerical stability. The difference qt,1 is computed using the sediment feed rate at the ghost node:
The worksheet RTe-bookAgDegBWChezy.xls provides both a grapical user interface and code in Visual Basic for Applications (VBA). A tutorial on VBA is provided in the workbook Rte-bookIntroVBA; it introduces the concept of modules. The code for the morphodynamic model is contained in Module 1 of Rte-bookAgDegBWChezy.xls. It can be seen by clicking “Tools”, “Macros”, “Visual Basic Editor” from Excel, and then double-clicking “Module1” in the VBA Project Window at the upper left of the Screen. The Security Level (“Tools”, “Macro”, “Security”) must be set to no higher than “medium” in order to run the code.
Most of the input is specified in worksheet “Calculator”. The first set of input includes: water discharge per unit width qw at flood, flood intermittency If, grain size D, reach length L, Chezy resistance coefficient Cz, antecedent bed slope Sa and volume total bed material feed rate per unit width during floods qtf. The specified numbers are then used to compute the normal flow depth Hna at antecedent conditions, the final ultimate equilibrium bed slope Su and the final ultimate equilibrium normal flow depth Hnu.
The user then specifies a downstream water surface elevation d. This
value should be > the larger of either Hna or Hnu to cause delta formation.
The following input parameters are then specified on worksheet “Calculator” by the user: reach length L, time step t, the number of time steps until data is generated for output (by printing it onto another worksheet in the workbook) Ntoprint, the number of times data is generated Nprint, number of spatial intervals M and upwinding parameter au. The total duration of the calculation is thus equal to t x Ntoprint x Nprint, and the spatial step length x = equal to L/M.
The parameter R is specified in worksheet “AuxiliaryParameter”.
Once all the input parameters are specified, the code is executed by clicking the button “Do the Calculation” in worksheet “Calculator”.
The numerical output is printed onto worksheet “ResultsofCalc”. The output consists of the position x, bed elevation , water surface elevation and flow depth H at every node for time t = 0 and Nprint subsequent times. The bed elevations and final water surface elevations are plotted on worksheet “PlottheData”.
In worksheet “Calculator” the flow discharge qw (m2/s) and bed material feed rate at flood flow qtf (m2/s) are specified per unit channel width.
In worksheet “MeanAnnualFeedRate” the user can specify a channel width Bf at flood flow (e.g. bankfull width). The flood discharge Qf = qf Bf in m3/s and the mean annual bed material feed rate Gma are then computed directly on the worksheet.
The input for all the cases (Cases A ~ G) illustrated subsequently in this presentation is given in worksheet “WorkedCases”.
As noted in Slide 25, the code itself can be viewed by clicking “Tools”, “Macros” and “Visual Basic Editor”, and then double-clicking Module 1 in the VBA Project Window in the upper left of the screen. Each unit of the code is termed a “Sub” or a “Function” in VBA. Three of these units are illustrated in Slides 29, 30 and 31.
The graphical user interface in worksheet “Calculator” is shown below.
Dim Hpred As Double: Dim fr2p As Double: Dim fr2 As Double:
Dim fnp As Double: Dim fn As Double: Dim Cf As Double
Dim i As Integer
H(M + 1) = xio - eta(M + 1)
For i = 1 To M
fr2p = qw ^ 2 / g / H(M + 2 - i) ^ 3
Cf = (1 / alr ^ 2) * (H(M + 2 - i) / kc) ^ (-1 / 3)
fnp = (eta(M + 1 - i) - eta(M + 2 - i) - Cf * fr2p * dx) / (1 - fr2p)
Hpred = H(M + 2 - i) - fnp
fr2 = qw ^ 2 / g / Hpred ^ 3
fn = (eta(M + 1 - i) - eta(M + 2 - i) - Cf * fr2 * dx) / (1 - fr2)
H(M + 1 - i) = H(M + 2 - i) - 0.5 * (fnp + fn)
For i = 1 To M
xi(i) = eta(i) + H(i)
Dim i As Integer
Dim taux As Double: Dim qstarx As Double: Dim Cfx As Double
For i = 1 To M + 1
taux = Cfx * (qw / H(i)) ^ 2 / (Rr * g * D)
qstarx = alt/cFS *taux ^ 2.5
qt(i) = ((Rr * g * D) ^ 0.5) * D * qstarx
Dim i As Integer
Dim qtback As Double: Dim qtit As Double: Dim qtfrnt As Double: Dim qtdif As Double
For i = 1 To M
If i = 1 Then
qtback = qtf
qtback = qt(i - 1)
qtit = qt(i)
qtfrnt = qt(i + 1)
qtdif = au * (qtback - qtit) + (1 - au) * (qtit - qtfrnt)
eta(i) = eta(i) + dt / (1 - lamp) / dx * qtdif * Inter
qtdif = qt(M) - qt(M + 1)
eta(M + 1) = eta(M + 1) + dt / (1 - lamp) / dx * qtdif * Inter
time = time + dt
This is a case for which a) base level d is unaltered from the antecedent value, and b) the bed material feed rate qtf is equal to the antecedent normal equilibrium value qtna. Condition a) is ensured by setting d equal to the antecedent normal equilibrium depth Hna, and condition b) is ensured by setting the bed material feed rate qtf so that the ultimate equilibrium bed slope Su is equal to the antecedent bed slope Sa. The result is the intended one: nothing happens over the 600-year duration of the calculation.
This case has the same input as Case A, except that the downstream water surface elevation is raised from 8.89 m to 20 m. The duration of the calculation is 6 years. A delta starts to form at the upstream end!
Same conditions as Case B, but the time duration is 30 years.
Same conditions as Case B, but the time duration is 60 years.
Same conditions as Case B, but the time duration is 120 years.
Same conditions as Case B, but the time duration is 600 years.
The “dam” is filled and overtopped, and ultimate normal mobile-bed equilibrium is achieved.
Same conditions as Case D, but the upstream feed rate qtf is tripled compared to the antecedent normal value qtna, so ensuring that Su is greater than Sa
This slope is steeper than this slope
Note that the delta front has prograded out much farther than Case D because the sediment feed rate is three times higher.
Exner, F. M., 1920, Zur Physik der Dunen, Sitzber. Akad. Wiss Wien, Part IIa, Bd. 129 (in German).
Exner, F. M., 1925, Uber die Wechselwirkung zwischen Wasser und Geschiebe in Flussen, Sitzber. Akad. Wiss Wien, Part IIa, Bd. 134 (in German).
Hotchkiss, R. H. and Parker, G., 1991, Shock fitting of aggradational profiles due to backwater, Journal of Hydraulic Engineering, 117(9): 1129‑1144.
Kostic, S. and Parker, G., 2003a, Progradational sand-mud deltas in lakes and reservoirs. Part 1. Theory and numerical modeling, Journal of Hydraulic Research, 41(2), pp. 127-140.
Kostic, S. and Parker, G.. 2003b, Progradational sand-mud deltas in lakes and reservoirs. Part 2. Experiment and numerical simulation, Journal of Hydraulic Research, 41(2), pp. 141-152.
Paola, C., Heller, P. L. & Angevine, C. L., 1992, The large-scale dynamics of grain-size variation in alluvial basins. I: Theory, Basin Research, 4, 73-90.
Parker, G., Sequeiros, O. and River Morphodynamics Class of Spring 2006, 2006, Large scale river morphodynamics: application to the Mississippi Delta, Proceedings, River Flow 2006 Conference, Lisbon, Portugal, September 6 – 8, 2006, Balkema.
Smith, W. O., Vetter, C.P. and Cummings, G. B., 1960, Comprehensive survey of Lake Mead, 1948-1949: Professional Paper 295, U.S. Geological Survey, 254 p.
Wright, S. and Parker, G., 2005a, Modeling downstream fining in sand-bed rivers. I: Formulation, Journal of Hydraulic Research, 43(6), 612-619.
Wright, S. and Parker, G., 2005b, Modeling downstream fining in sand-bed rivers. II: Application, Journal of Hydraulic Research, 43(6), 620-630.