of 14
JOURNAL OF GEOPHYSICAL RESEARCH
Supporting Information for “Modeling deltaic lobe-building
cycles and channel avulsions for the Yellow River delta,
China”
Andrew J. Moodie
1
, Jeffrey A. Nittrouer
1
, Honbgo Ma
1
, Brandee N. Carlson
1
,
Austin J. Chadwick
2
, Michael P. Lamb
2
, and Gary Parker
3,4
Contents of this file
1. Text S1 to S7
2. Figures S1 to S6
3. Table S1
1
Department of Earth, Environmental and
Planetary Sciences, Rice University, Houston,
Texas, USA
2
Division of Geological and Planetary
Sciences, California Institute of Technology,
Pasadena, California, USA
3
Department of Civil and Environmental
Engineering, University of Illinois at
Urbana-Champaign, Urbana, Illinois, USA
4
Department of Geology, University of
Illinois at Urbana-Champaign, Champaign,
Illinois, USA
1
X - 2
MOODIE ET AL.: MODELING DELTAIC LOBE-BUILDING CYCLES
1. Yellow River delta growth measurements
The mapped record of Yellow River delta growth through multiple avulsion cycles is pre-
sented in the main text in Figure 2b (Pang & Si, 1979; van Gelder et al., 1994; Ganti et al.,
2014). Ganti et al. (2014) made measurements of the timing and location of avulsion from this
same record, as well as a map of Yellow River deltaic avulsion from (Fan et al., 2006) (their
supplementary material, Figure S1 herein). Ganti et al. (2014) find the avulsion time scale for
the natural Yellow River delta system to be (
T
A
,
YR
=
7
±
2 yr); this result is used in the main text.
However, the values for “avulsion length” reported by Ganti et al. (2014) in their Supporting
Information are actually measurements of the streamwise distance from the datum of Lijin to
the avulsion location (pers. communication with V. Ganti), reflecting values tabulated for their
avulsion node forward-stepping rate calculation.
Avulsion length is measured for the Yellow River delta record in this study using both Figure 3
from the main text, as well as Figure S1. The measurements were made by mapping the figures
to a Cartesian coordinate system scaled by the figure’s scale bar, and distance was calculated
streamwise from avulsion location to contemporaneous shoreline (Table S1). Importantly, the
values measured differ from the measurements made by Ganti et al. (2014), though it is unclear
what leads to these differences. The mean and standard deviation calculated from measurements
of Figure 2b is used for the discussion of the Yellow River delta avulsion length in the main text
(
L
A
,
YR
=
52
.
5
±
12
.
3 km) because Figure S1 has additional distributary channels and poorly
constrained avulsion timing, thus increasing the overall uncertainty associated with these data.
The measurements from Figure 2b are characterized in the boxplot in Figure 3 in the main text.
2. Redistributing delta lobe sediment
The sediment mass distribution of the delta lobe is estimated numerically following a formu-
lation that incorporates the three-dimensionality of the quasi-2D delta by considering a delta that
is symmetric about the channel. In this way the formulation is cast in a streamwise coordinate
x
and vertical coordinate
η
.
The volume of the delta at time
t
is found by the intersection of three geometric volumes:
The “cylinder” (Figure S2) volume is given by a cylinder sector centered at
x
=
0
,
η
=
h
l
with opening angle
Γ
, radius
r
t
, and height
h
l
.
The “total cone” (Figure S2) volume is given by a cone sector concentric to the cylinder
with opening angle
Γ
, but with height
h
l
+
h
u
, where
h
u
is the elevation of the basin surface at
x
=
0.
The “upper cone” (Figure S2) volume is given by a cone sector centered at
x
=
0
,
η
=
0
with opening angle
Γ
, radius
x
0
, and height
h
u
, where
x
0
is the length from
x
=
0 to
x
(
η
=
0
)
.
Let
Γ
rad
the opening angle in radians. The volume of the delta at time
t
is then formally given:
V
(
t
1
) =
cylinder
(
total cone
upper cone
)
(1)
V
(
t
1
) =
Γ
rad
π
r
2
1
h
l
1
(
Γ
rad
π
r
2
1
h
l
1
3
Γ
rad
π
x
2
0
h
u
3
)
(2)
The change in volume
V
is given by the different in delta volumes at times
t
2
and
t
1
:
V
=
V
(
t
2
)
V
(
t
1
)
(3)
MOODIE ET AL.: MODELING DELTAIC LOBE-BUILDING CYCLES
X - 3
V
=
[
Γ
rad
π
r
2
2
h
l
2
(
Γ
rad
π
r
2
2
h
l
2
3
Γ
rad
π
x
2
0
h
u
3
)
]
[
Γ
rad
π
r
2
1
h
l
1
(
Γ
rad
π
r
2
1
h
l
1
3
Γ
rad
π
x
2
0
h
u
3
)
]
(4)
Rearranging and simplifying:
Γ
rad
π
(
r
2
2
h
l
2
r
2
2
h
l
2
3
)
=
V
+
Γ
rad
π
(
r
2
1
h
l
1
r
2
1
h
l
1
3
)
(5)
r
2
2
h
l
2
=
3
2
V
Γ
rad
π
+
r
2
1
h
l
1
(6)
which leaves one equation and two unknowns. The formulation for the new delta radius
r
2
given
a change in volume
V
, is solved numerically to a tolerance of 1 m.
3. Deltaic lobe width estimate
It is necessary to parameterize the width of deltaic lobes in this model because there is no
existing theoretical framework predicting the width of lobes. The width of a deltaic lobe likely
changes over time and space, but the model formulation requires a constant width parameteri-
zation and makes an implicit assumption that lobes are symmetric across the channel axis. This
parameterization is determined by measuring the downstream variation in lobe width of 3 lobes
from the Yellow River delta just before the lobes were abandoned by avulsion. Additionally,
the modern configuration of the active Yellow River delta lobe is measured. The half width of
each delta lobe is measured at 1 km intervals moving downstream from the lobe apex along the
channel centerline when it can be identified, otherwise the lobe centerline is used. Satellite im-
ages used in measurements and results are presented in Figure S3. Measurements are tabulated
in Table S2.
The lobe half-width measurements show variation at upstream width, but consistently decease
in width moving downstream to the end of the lobe. There is a mode of lobe half-widths in the
range of 4–5 km, and so the model is parameterized to a lobe width of 9 km (2
×
half-width
=
2
×
4
.
5).
4. Channel bed slope estimate
The initial channel-bed slope for the model runs is set as a constant (
S
0
=
6
.
4
×
10
5
). Slopes
reported for the Yellow River span from 5
×
10
5
to 1
.
0
×
10
4
(Chunhong et al., 2005; Ganti et
al., 2014). The initial channel-bed slope is derived from measurement of the water-surface slope
in the normal flow reach, that is, upstream of the backwater. A record of shipboard navigation
20 km long was collected between the cities of Kenli and Lijin (80–100 km upstream of
modern mouth). This navigation record is processed by projecting all navigation data onto the
channel centerline. This allows the record to be presented as elevation of the water surface
moving downstream from the city of Lijin (Figure S4).
Substantial river engineering along the navigated river course results in changes in flow ve-
locity that affect the draft of the vessel during transit. This produces significant excursions in the
water surface slope around bends in the channel course. To minimize the impact of bends on the
extracted water surface slope, only the elevation data from the longest straight-reach river sec-
tion are fit with a linear regression (black section, Figure S4). The resulting regression value is
X - 4
MOODIE ET AL.: MODELING DELTAIC LOBE-BUILDING CYCLES
6
.
35
×
10
5
, which is on the lower side of all reported slope values for the Yellow River, though
there is considerable noise in this section of the river as well (r
2
= 0.09). A sensitivity analysis
varying the upstream and downstream extent of the river used for the regression determines that
this a robust fit.
The fit slope used is rounded up to 6
.
4
×
10
5
as the initial channel-bed slope in all model
runs. This slope is within the range of the reported slopes for the Yellow River delta. Fur-
thermore, variations in slope are not expected to have a significant effect on the overall trends
predicted by the model, however, the precise values of predicted avulsion time and location will
vary as a result of initial-slope changes.
5. Bankfull and formative discharges
An important input to the model is the parameterization of the bankfull and formative dis-
charges of the lower Yellow River and delta. Bankfull discharge may be easily measured in the
modern river by simply observing the discharge for which the flow just begins to go overbank.
The bankfull discharge of the Yellow River at the city of Zhengzhou is reported by (Wang &
Liang, 2000) as
Q
w
,
b f
=
3000 m
3
/s. However, further downstream near the city of Lijin and the
delta, the bankfull discharge is higher and flow is contained within the channel for a discharge
close to 3900 m
3
/s (Wang & Liang, 2000).
The lower bankfull discharge (i.e., measured near Zhengzhou) is used in the model because
this channel configuration is interpreted to reflect the natural incision of the channel, though both
measurements reflect a highly engineered channel. The model incision depth is particularly
relevant because the bankfull normal flow depth (combined with the basin slope) is used to
determine the channel depth everywhere at the onset of the model run.
Maximum geomorphic work is often used to provide an explanation for the self-forming equi-
librium geometry of channels, whereby the bankfull discharge “is the flood whose combination
of frequency times magnitude moves the most sediment,” (Wolman & Miller, 1960; Jerolmack
& Brzinski, 2010). Small floods are frequent but transport little sediment, in contrast to large
floods which transport a substantial volume of sediment, yet are relatively rare. For the Yellow
River, (Ma et al., n.d.) demonstrate that significant sediment transport occurs even at low dis-
charges, inviting the possibility that the flood doing the maximum geomorphic work is not the
bankfull flood.
The geomorphic work (
G
) is defined by the product of the probability that a flow exceeds a
certain flood (
P
) and the magnitude of sediment transport occurring during that flow (
Q
s
). The
exceedance probability for a flow is determined by ranking the daily water discharge data from
the Lijin Hydrological Station (e.g., Figure 6a) and producing a complementary cumulative
distribution function (Figure S5a) following (Jerolmack & Brzinski, 2010). An exponential
function is transformed so that a linear model is fit to the distribution:
P
=
e
Q
w
log
(
P
) =
mQ
w
(7)
where
Q
w
is the ranked water discharge measurements and
m
is the regressed slope.
The sediment transport (
Q
s
) for each unique discharge measurement is then calculated using
the transport equation from Ma et al. (n.d.) (i.e., Equation 4), dimensionalized and adjusted to
MOODIE ET AL.: MODELING DELTAIC LOBE-BUILDING CYCLES
X - 5
reflect transport integrated across the flow width:
Q
s
=
B
RgD
50
3
α
C
f
τ
n
(8)
where
τ
=
HS
/
RD
50
is the Shields’ number,
H
is the flow depth calculated by a generalized
Darcy-Weisbach equation (Ganti et al., 2014), median grain diameter
D
50
is set to 90
μ
m,
α
=
0
.
895 and
n
=
1
.
678 are adjusted coefficients for sediment transport at Lijin (Ma et al., n.d.),
B
=
400 m is the flow width,
C
f
=
0
.
001 is the friction coefficient, slope is set to
S
=
6
.
4
×
10
5
,
submerged specific gravity is
R
=
1
.
65, and the gravitational acceleration constant is
g
=
9
.
81.
The sediment transport predictions are then used to calculated geomorphic work (
G
) as:
G
=
P
×
Q
s
(9)
and plotted versus the corresponding water discharge in Figure S5b. The linear model obtained
in Equation 7 is transformed to geomorphic work in a similar way, using the sediment transport
relation (Equation 8) to properly scale water discharges to corresponding geomorphic work
units (Figure S5b).
The regression model exhibits a well defined peak corresponding to a water discharge of
1300 m
3
/s. However, the regression model over-predicts the data in the water discharge range
from 1000 to 1500 m
3
/s, because flows in this range have roughly the same geomorphic work
potential. Theoretically, a poorly defined peak in the data suggests that all discharges in this
range should have a similar effect on the river system. However, a sensitivity test reveals that any
formative discharge within this range does not affect the outcome of the model considerably. for
example, a higher formative discharge leads to a mouth that progrades earlier in the model run,
yet the rate of lobe progradation is unaffected because the lobe growth rate is still dominantly a
function of sediment input. Therefore, the formative discharge is identified as the well defined
peak of the regression model:
Q
w
=
1300 m
3
/s.
6. Coastline processing code explanation
The method for processing the coastline used in this analysis is an adaptation of meth-
ods described by many authors for a variety of applications. The method involves thresh-
olding a reflectance (intensity) image to convert to a binary raster and a series of pro-
cessing steps to extract the coastline vector.
The code is written for Matlab 2015a
on Ubuntu 16.04; the code is available at https://github.com/amoodie/paper
resources and
https://github.com/amoodie/shoreline
extract.
The goal of this paper is not to present a novel and improved coastline extraction method, but
rather to describe the method employed for the purpose of complete information to the reader.
The general approach follows the methods of countless other coastline extraction routines (e.g.,
Lee and Jurkevich (1990); Alesheikh, Ghorbanali, and Nouri (2007)), but the precise sequence
of operations and parameterization were tweaked until a visual inspection indicated satisfactory
results. In particular, this approach worked well for the spatial extent of the Yellow River delta
system used in this study, but modifications are likely necessary for other applications.
The first, and arguably most important, step in the coastline extraction process is to convert
the reflectance image to a binary raster. A threshold reflectance value was calculated for each
Landsat Band 7 image, and was obtained by a series of steps (Figure S6):
X - 6
MOODIE ET AL.: MODELING DELTAIC LOBE-BUILDING CYCLES
1. build a histogram (
C
=
count) of binned pixel reflectance values (
R
) in Band 7 image
2. calculate the gradient of histogram bin counts (change in count over a change in reflectance
value,
Gr
=
dC
/
dR
)
3. identify the reflectance value with maximum count (
R
max
), and the smallest reflectance value
larger than
R
max
with an increasing gradient (
R
ig
)
4. find the mean gradient (
Gr
mean
) of all reflectance values between
R
max
and
R
ig
5. determine the x-intercept for a line from the point [
R
max
,
C
(
R
max
)
] with a slope of
Gr
mean
The value of the x-intercept is the threshold (
T
i
) used to binarize the
i
th image.
An image was binarized by declaring pixels with reflectance values
T
i
as water, and pixels
with reflectance values
>
T
i
as land. Then, a series of image processing operations were applied
to extract the coastline:
1. flood-fill operation to fill isolated water within land which cannot be reached by flooding
from edges
2. remove small isolated water-on-land objects (area
<
30000 px)
3. remove small isolated land-in-water objects (area
<
500 px)
4. morphological opening operation with disk-shaped filter (radius
=
5 px)
5. morphological closing operation with disk-shaped filter (radius
=
50 px)
6. remove small isolated land-in-water objects (area
<
10000 px)
7. flood-fill operation to fill water pixels, with padded edges (final clean-up operation)
8. retain only the largest contiguous land portion (i.e., the delta)
9. Sobel edge detection operation to identify coastline
Uncertainty estimates on the extracted shoreline positions are made by assuming a tidal am-
plitude of 1.5 m near the delta edge (B. Carlson, personal communication on unpublished data),
and a topset slope in the vicinity of the delta edge of 3
×
10
4
. This produces an error estimate
of
∼±
3 km.
7. Derivation of gradually varying flow equation
Given the steady one-dimensional mass- and momentum- conservation equations for an
arbitrarily-shaped channel (e.g., Chow (1959)):
UH
x
=
0
(10)
U
2
H
x
=
1
2
g
H
2
x
gH
∂ η
x
C
f
U
2
,
(11)
where
U
is depth-averaged streamwise flow velocity,
H
is flow depth,
x
is the streamwise spatial
coordinate,
g
is the gravitational acceleration constant,
η
is the channel bed elevation, and
C
f
is
the friction coefficient of the flow. Applying the chain rule and equation 10 to equation 11, the
left hand side becomes:
U
2
H
x
U
UH
x
+
UH
U
x
UH
U
x
,
(12)
and given the substitutive definition
∂ η
/
x
=
S
0
, the right hand side becomes:
1
2
g
H
2
x
→−
gH
H
x
+
gHS
0
C
f
U
2
,
(13)
MOODIE ET AL.: MODELING DELTAIC LOBE-BUILDING CYCLES
X - 7
which combines to give the expression:
UH
U
x
=
gH
H
x
+
gHS
0
C
f
U
2
.
(14)
Assuming a rectangular channel that is uniform across the width gives the following substitutive
definitions:
Q
=
UHB
(15)
Fr
2
=
U
2
(
gA
/
B
)
=
Q
2
gB
2
H
3
=
U
2
gH
,
(16)
where
Q
is the flow discharge rate,
B
is the width of channel, and Fr is the Froude number. Thus,
equation 14 can be further simplified by rearranging:
U
g
U
x
=
H
x
+
S
0
C
f
Fr
2
.
(17)
The spatial derivative of velocity (
U
/
x
) can be recast in terms of the discharge (
Q
), given
equation 15:
U
x
=
x
(
Q
HB
)
(18)
=
Q
H
2
B
H
x
Q
HB
2
B
x
(19)
=
U
H
H
x
U
B
B
x
,
(20)
which allows equation 17 to be simplified after substituting equation 20:
U
g
(
U
H
H
x
U
B
B
x
)
=
H
x
+
S
0
C
f
Fr
2
(21)
U
2
gH
H
x
U
2
gB
B
x
=
H
x
+
S
0
C
f
Fr
2
(22)
U
2
gH
H
x
+
H
x
=
S
0
C
f
Fr
2
+
U
2
gB
B
x
(23)
H
x
(
U
2
gH
+
1
)
=
S
0
C
f
Fr
2
+
U
2
gB
B
x
.
(24)
Applying equation 16 to equation 24 above yields:
H
x
(
Fr
2
+
1
)
=
S
0
C
f
Fr
2
+
Fr
2
H
B
B
x
,
(25)
which is rearranged to give a formulation for gradually varying flow depth through a rectangular
open-channel flow as function of the streamwise coordinate:
H
x
=
S
0
C
f
Fr
2
(
1
Fr
2
)
+
Fr
2
(
1
Fr
2
)
H
B
B
x
,
(26)
X - 8
MOODIE ET AL.: MODELING DELTAIC LOBE-BUILDING CYCLES
which is the same as equation 3 in the main manuscript.
If the assumption of a wide and shallow channel is imposed (i.e.,
B

H
), equation 26 can
be further simplified because
H
/
B
0, yielding the more common formulation for gradually
varying flow in delta hydrodynamics:
H
x
=
S
0
C
f
Fr
2
(
1
Fr
2
)
.
(27)
Alternatively, if the width is assumed to not change in the downstream direction (
B
/
x
0,
a similar simplification can be made, again yielding equation 27. Please take care to use the
definition for the Froude number that includes the width term when applying equation 26 (i.e.,
Fr
2
=
Q
2
/
B
2
H
3
).
References
Alesheikh, A. A., Ghorbanali, A., & Nouri, N. (2007, December). Coastline change detection
using remote sensing.
International Journal of Environmental Science & Technology
,
4
(1),
61–66. doi: 10.1007/BF03325962
Chow, V. T. (1959).
Open Channel Hydaulics
.
Chunhong, H., Yangui, W., Yanjing, Z., Yuling, T., Hongling, S., Cheng, L., . . . Deyi, W.
(2005, November).
CASE STUDY ON THE YELLOW RIVER SEDIMENTATION
(case
study). Beijing, China: INTERNATIONAL RESEARCH AND TRAINING CENTER
ON EROSION AND SEDIMENTATION.
Fan, H., Huang, H., Zeng, T. Q., & Wang, K. (2006, March). River mouth bar formation,
riverbed aggradation and channel migration in the modern Huanghe (Yellow) River delta,
China.
Geomorphology
,
74
(1-4), 124–136. doi: 10.1016/j.geomorph.2005.08.015
Ganti, V., Chu, Z., Lamb, M. P., Nittrouer, J. A., & Parker, G. (2014, November). Testing
morphodynamic controls on the location and frequency of river avulsions on fans versus
deltas: Huanghe (Yellow River), China: Avulsion drivers on fans versus deltas.
Geophysi-
cal Research Letters
,
41
(22), 7882–7890. doi: 10.1002/2014GL061918
Jerolmack, D. J., & Brzinski, T. A. (2010, August). Equivalence of abrupt grain-size transitions
in alluvial rivers and eolian sand seas: A hypothesis.
Geology
,
38
(8), 719–722. doi:
10.1130/G30922.1
Lee, J.-s., & Jurkevich, I. (1990, July). Coastline Detection And Tracing In SAr Images.
IEEE Transactions on Geoscience and Remote Sensing
,
28
(4), 662–668. doi: 10.1109/
TGRS.1990.572976
Ma, H., Nittrouer, J. A., Naito, K., Fu, X., Zhang, Y., Moodie, A. J., . . . Parker, G. (n.d.). The
exceptional sediment load of fine-grain dispersal systems: Example of the yellow river,
china. ,
3
, 7. doi: 10.1126/sciadv.1603114
Pang, J. Z., & Si, S. H. (1979). Evolution of the Yellow River mouth: I. Historical shifts.
Oceonologia et Limnologia Sinica
,
10
(2), 136–141.
van Gelder, A., van den Berg, J. H., Cheng, G., & Xue, C. (1994, May). Overbank and
channelfill deposits of the modern Yellow River delta.
Sedimentary Geology
,
90
(3–4),
293–305. doi: 10.1016/0037-0738(94)90044-2
Wang, Z.-Y., & Liang, Z.-Y. (2000, July). Dynamic characteristics of the Yellow River
mouth.
Earth Surface Processes and Landforms
,
25
(7), 765–782. doi: 10.1002/1096
-9837(200007)25:7$
$765::AID-ESP98$
$3.0.CO;2-K
MOODIE ET AL.: MODELING DELTAIC LOBE-BUILDING CYCLES
X - 9
Wolman, M. G., & Miller, J. P. (1960, January). Magnitude and Frequency of Forces in
Geomorphic Processes.
The Journal of Geology
,
68
(1), 54–74. doi: 10.1086/626637
do not specify file extension
X - 10
MOODIE ET AL.: MODELING DELTAIC LOBE-BUILDING CYCLES
Figure S1.
Map reproduced from Ganti et al. (2014) Supporting Information, figure after Fan et al.
(2006). Measurement of avulsion length and timing are made from this record and tabulated in Table
S1.
basin
channel bed
topset
x=0
x=r
x=x
0
η
=0
η
0
h
u
h
l
r
=
r
x
0
Figure S2.
Schematic of delta in dip section. The change in delta radius is calculated by the intersec-
tion of three-dimensional geometric volumes, for a given delta lobe volume
V
.