MORPHODYNAMICS OF SAND-BED RIVERS ENDING IN DELTAS

1
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
2
DELTA SHAPE
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.
2D progradation
1D progradation
Mangoky River, Malagasy
Colorado River at Lake Mead, USA
3
AN EXAMPLE OF A 1D DELTA
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 https//zulu.ssc.nasa.g
ov/mrsid/mrsid.pl
4
HISTORY OF SEDIMENTATION IN LAKE MEAD, 1936 - 1948
Image based on an original from Smith et al.
(1960)
5
WHERE DOES THE SAND GO? WHERE DOES THE MUD GO?
Image based on an original from Smith et al.
(1960)
6
STRUCTURE OF A DELTA TOPSET, FORESET AND
BOTTOMSET
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.
standing water
7
(Kostic and Parker, 2003a,b)
8
SIMPLIFICATION TOPSET AND FORESET ONLY
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).
standing water
9
DERIVATION OF THE 1D EXNER EQUATION OF SEDIMENT
CONSERVATION
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
or
control volume
or
10
BUT HOW ARE FLOODS ACCOUNTED FOR IN THE EXNER
EQUATION OF SEDIMENT CONTINUITY?
Fly River, PNG in flood
11
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
12
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 or thus
13
DOWNSTREAM BOUNDARY CONDITION FOR THE BACKWATER
FORMULATION
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 https//zulu.ssc.nasa.gov/
mrsid/mrsid.pl
14
THE MORPHODYNAMIC PROBLEM
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.
Exner equation
Bed material load equation
Backwater 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.
15
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
16
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 cant reach em very well either.
17
THE PROBLEM OF IMPULSIVELY RAISED WATER SURFACE
ELEVATION (BASE LEVEL) AT t 0
M1 backwater curve
qta
qta
Note the M1 backwater curve was introduced in
the lecture on hydraulics and sediment transport
18
RESPONSE TO IMPULSIVELY RAISED WATER SURFACE
ELEVATION A PROGRADING DELTA THAT FILLS THE
SPACE CREATED BY BACKWATER
See Hotchkiss and Parker (1991) for more details.
19
FLOW OVER ANTECEDENT BED
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
20
CHANGES IMPOSED AT t 0
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.
21
ULTIMATE EQUILIBRIUM
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
22
NUMERICAL MODEL INITIAL AND BOUNDARY CONDITIONS
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 gt 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.
23
NUMERICAL MODEL DISCRETIZATION AND BACKWATER
CURVE
The reach L is discretized into M intervals of
length ?x bounded by M1 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.
24
NUMERICAL MODEL SEDIMENT TRANSPORT AND EXNER
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
25
INTRODUCTION TO RTe-bookAgDegBWChezy.xls
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 gt the larger
of either Hna or Hnu to cause delta formation.
26
INTRODUCTION TO RTe-bookAgDegBWChezy.xls contd.
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.
27
INTRODUCTION TO RTe-bookAgDegBWChezy.xls contd.
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.
28
GRAPHICAL USER INTERFACE
The graphical user interface in worksheet
Calculator is shown below.
29
BACKWATER CALCULATION Sub Do_Fluvial_Backwater
Sub Do_Fluvial_Backwater() 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) Next i For
i 1 To M xi(i) eta(i) H(i)
Next i End Sub
30
LOAD CALCULATION Sub Find_Shields_Stress_and_Load
Sub Find_Shields_Stress_and_Load() 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 Next i End Sub
31
IMPLEMENTATION OF EXNER Sub Find_New_eta
Sub Find_New_eta() 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 Else
qtback qt(i - 1) End If
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 Next i
qtdif qt(M) - qt(M 1) eta(M
1) eta(M 1) dt / (1 - lamp) / dx qtdif
Inter time time dt End Sub
32
CASE A
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.
33
CASE B
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!
delta!
34
CASE C
Same conditions as Case B, but the time duration
is 30 years.
35
CASE D
Same conditions as Case B, but the time duration
is 60 years.
36
CASE E
Same conditions as Case B, but the time duration
is 120 years.
37
CASE F
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.
38
CASE G
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.
39
GENERALIZATIONS

• The basic formulation can be generalized to
  include
• a self-formed channel with varying width,
• channel sinuosity,
• a channel-floodplain complex in which mud as
  well as sand can deposit,
• a foreset of specified slope and
• a 2D geometry that yields a fan shape to the
  delta.
• These generalizations are implemented for the Wax
  Lake Delta shown below in the paper Large scale
  river morphodynamics application to the
  Mississippi Delta (Parker et al,. 2006) included
  on the CD for this course as file WaxLake.pdf.
  

40
REFERENCES
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.