Endo Progress in Earth and Planetary Science (2016) 3:28
DOI 10.1186/s40645-016-0102-9

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Simple stochastic cellular automaton model
for starved beds and implications about
formation of sand topographic features interms of sand flux
Noritaka Endo

Abstract
A simple stochastic cellular automaton model is proposed for simulating bedload transport, especially for caseswith a low transport rate and where available sediments are very sparse on substrates in a subaqueous system.
Numerical simulations show that the bed type changes from sheet flow through sand patches to ripples as theamount of sand increases; this is consistent with observations in flume experiments and in the field. Withoutchanges in external conditions, the sand flux calculated for a given amount of sand decreases over time as
bedforms develop from a flat bed. This appears to be inconsistent with the general understanding that sand flux
remains unchanged under the constant-fluid condition, but it is consistent with the previous experimental data. For
areas of low sand abundance, the sand flux versus sand amount (flux–density relation) in the simulation shows asingle peak with an abrupt decrease, followed by a long tail; this is very similar to the flux–density relation seen in
automobile traffic flow. This pattern (the relation between segments of the curve and the corresponding bed
states) suggests that sand sheets, sand patches, and sand ripples correspond respectively to the free-flow phase,
congested phase, and jam phase of traffic flows. This implies that sand topographic features on starved beds are
determined by the degree of interference between sand particles. Although the present study deals with simple
cases only, this can provide a simplified but effective modeling of the more complicated sediment transport
processes controlled by interference due to contact between grains, such as the pulsatory migration of grain-size
bimodal mixtures with repetition of clustering and scattering.
Keywords: Stochastic cellar automaton, Bedload transport, Sand patches, Sand sheets, Sand ripples

Correspondence: wisteria@staff.kanazawa-u.ac.jp
Department of Earth Sciences, Graduate School of Natural Science &Technology, Kanazawa University, Kakuma-machi, Kanazawa, Ishikawa
920-1192, Japan© 2016 The Author(s). Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0
International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to
the Creative Commons license, and indicate if changes were made.

Endo Progress in Earth and Planetary Science (2016) 3:28

Introduction
Bedforms, including both subaerial and subaqueous ones,
have been studied for more than half a century, mainly
starting with the pioneering field observations of Bagnold(1941) and the laboratory studies of Allen (1968) and continuing to recent comprehensive experiments (e.g., Perillo
et al. 2014). Theoretical studies also have been conducted,
employing various approaches, such as linear and nonlinear stability analyses (e.g., Melo et al. 2012; Andreotti et al.
2012); first-principle calculations for fluids and/or grains(e.g., Wippermann and Gross 1986; Schwämmle andHerrmann 2003; Duran et al. 2014); simulations with simplified continuous models, such as a defect-crest line
dynamics model (Werner and Kocurek 1999), a skeleton
model (Niiya et al. 2010), and a single-dune-based particle
model (Parteli and Herrmann 2003; Diniega et al. 2010);
and simulations with various cellular automaton (CA)models (e.g., Nishimori and Ouchi 1993; Werner 1995;
Momiji et al. 2000; Katsuki et al. 2005). Even with increased computing speeds, CA models of bedforms
remain useful for complicated settings that are time consuming when other methods are used (Barchyn and
Hugenholtz 2012; Barrio-Parra and Rodríguez-Santalla
2014; Rozier and Narteau 2014).
The transport of sand grains in water can be classified
into four modes: sliding, rolling, saltation, and suspension (Allen 1982). Sliding and rolling grains remain incontinuous contact with the bed, while suspended particles have little or no contact with the bed. Saltating
grains jump over the bed in a series of low trajectories
that are on the order of 10 to 20 particle diameters long
and a few diameters high; thus, this is an intermediate
mode between sliding and rolling and suspension. Reidand Frostick (1994) defined the bedload as the grainsthat move along the bed and that are essentially incontact with it at all times (i.e., including only sliding
and rolling particles) and independently defined the saltation load; however, others (e.g., Allen 1982; Bridge and
Demicco 2008) defined the bedload to include saltating
grains in addition to sliding and rolling grains. Regardless of the definition of bedload, all researchers havepointed out that in subaqueous environments, sliding
and rolling are most prevalent when the level of transport is low (low transport stage; Allen 1982; Reid andFrostick 1994; Bridge and Demicco 2008).
Allen (1982) pointed out that in air and water, transport processes differ due to the difference in viscosityand density of the medium. In the case of aeolian dunes,
saltation is the most important mode (Bagnold 1941),
whereas bedload transport is most important for subaqueous ripples (Allen 1968). Sand ripples can be generated by a shear stress that is just above the threshold for
bedload transport (Bridge and Demicco 2008). In subaqueous environments, sliding and rolling are mostPage 2 of 11

prevalent at low transport stages (Allen 1982; Reid andFrostick 1994; Bridge and Demicco 2008), and therefore,
instead of considering changes in saltation length, thesubaqueous bedform at low transport stages should be
explained by a model that considers grain movementsthat are in continuous contact with the bed. Althoughprevious CA models for aeolian dunes have attributed
the bedform generation to changes in the saltation
length, in the present study, which employs a stochastic
model, we consider changes of the migration speed ofsand particles that are in continuous contact with thebed (i.e., sand that moves no faster than one cell pertime step). Although there have been a few preliminary
attempts to employ the CA model and explicitly include
the continuous bed-contact load (e.g., Miyata and Iwasa
1997), the relation between sand flux and sand abundance
was not considered.
Historically, most sedimentological and geomorphological studies of subaqueous beds have considered beds
in which sand covered the entire substrate, i.e., sedimentbeds. It is known that on starved beds (nonerodible
substrates with sparse sand), there occur topographic features that are different from those that occur on sedimentbeds, for example, sand patches. For starved beds, therelation between the sand flux and the topographic features has not been sufficiently studied. Although Betat
et al. (2002) measured the bedload flux for sediment beds,
there are few such studies, because it is difficult to measure the sand flux without disturbing the system; this isespecially true for starved beds. Numerical simulations
can provide useful information about the generation andtransition of sand topographic features on starved beds,
such as sand patches.
In this paper, we consider sparse sediments under relatively slow, unidirectional water flow; we present a CA
model and discuss the relation between bed states andsand flux. We compare the results of our model with theresults of a jam simulated by a traffic flow model. The
jam-like transportation of sediments has been reported,
such as sediments passing a meander bend of a channel
with a pool (Kleinhans et al. 2002) and the pulsatory
migration of sediment mixtures with grain-size bimodality (Iseya and Ikeda 1987). Following perturbations inthe sediment supply, sparsely distributed sediments in a
bedrock channel cause both negative and positive feedback to the incision rate (Finnegan et al. 2007). Althoughapplications to such complicated natural phenomena are
beyond the scope of this study, we note that the CAmodel is a potential tool for a deeper understanding of
the dynamics of starved beds.
The purposes of the present study are as follows: (1)to propose a new CA model that includes the continuous bed-contact load in small-scale subaqueous beds
with sparsely available sediments; (2) to compare theEndo Progress in Earth and Planetary Science (2016) 3:28

results of our model to those of flume experiments toassess how well it models changes in topographic features with increasing levels of sand abundance, as wellas temporal changes in the sand flux under constantconditions (this has not previously attracted much attention); and (3) to gain a better understanding of the dynamics of bed states by investigating the dependence ofthe sand flux on the abundance sand and topographic
features; we do this by a comparison with other kinds ofphysical phenomena, such as traffic flow.

Methods/Experimental

Page 3 of 11

ignores the inner structure of the flow, which becomes important when the flow is relatively fast andthe bed roughness (bedforms) is large; for simplicity,
we ignore the effects of strong erosional vortices andsecondary flows.
Formulation of bedload transport

When simulating bedload transport in a CA model, sandslabs should be restricted to moving to the adjacentdownstream cell, unlike the process for saltation. The
topographic height is affected by bedload transport asfollows:

CA modelCA models include the coarse graining of sand by considering clusters of sand grains (sand slabs) instead ofindividual particles. In CA models, the topographic
height is taken to be the number of compiled slabs.
The topographic height h(i, j) at site (i, j) in a twodimensional field changes with time. Previous CAmodels for aeolian sand dunes used a phenomenological formulation of saltation; that is, it was not basedon fluid motion. Although there is some variation in
the formulations of saltation, the previous models determined situational-specific saltation distances, suchas defining the jumping length as a function of thestate of the sand bed.
Allen (1982) pointed out that the differences in air andwater transport stemmed from differences in the viscosity and density of the medium. On subaqueous beds
under relatively slow flows, the main transport mode ofsand grains is rolling and sliding, for which bed contact
is continuous (Allen 1982; Reid and Frostick 1994;
Bridge and Demicco 2008). Therefore, instead of using
the saltation length, low transport stages of subaqueous
bedforms can be better explained by incorporatingchanges in the migration speed of grains in continuous
contact with the bed. In the simulation used in thepresent study, the movement of particles was consideredseparately for two processes: transport due to fluid
motion, with continuous bed contact; and avalanching
due to gravity. For the remainder of this paper, we willuse “bedload” to mean the continuous bed-contact loadthat is only due to the flow (not due to gravity). Other
processes due to flow, such as saltation or suspension,
will not be considered so that we can focus on subaqueous environments with low transport rates that form
topographies such as sand patches. In the simulation,
which process (bedload transport or avalanching) occurs
is evaluated deterministically from the geomorphologicalstate of the position. However, the movement of sandslabs due to bedload transport is a stochastic variable
that expresses finite velocities for discrete positions;that is, slabs migrate according to a probability that
represents the velocity (see below). The simulation

hði; jÞ→hði; jÞ−qb and hði þ 1; jÞ→hði þ 1; jÞ þ qb ;ð1Þ

where site (i + 1, j) is downstream of site (i, j), and qb is
the number of slabs that migrate during a single calculation step. Here, we consider bedload transport in only
the direction of the flow; this is the positive x-directionof the two-dimensional (x, y) system.
If the bed is perfectly uniform, sand may be transported everywhere at a constant rate, but otherwise, therate of sand migration depends on the location. We
formulated a simple phenomenological rule for the
probability of sand movement (surrogate of velocity)that does not require calculating the fluid motion. In
each time step, a sand slab on the surface moves to thenext downstream cell or remains at the same location,
as determined stochastically; this probability corresponds to the velocity of the sand migration velocity(that is, a higher probability corresponds to a higher
velocity), and this depends on the state of each site.
We made the following assumptions. Sand slabs on
a stoss side move fastest at the local peak (state s1),and other sand slabs (for example, partway down theslope) slow down as the downstream gradient increases, due to the competing gravitational force.
Sand slabs on lee slopes move more slowly than thoseon stoss slopes, because to some extent, the bedformshields the downstream flow. For a slab on a leeslope, the probability of movement depends on theextent of the shielding effect of the bedform; we assume that this depends on the gradient of the leeslope. When the slope of a lee side is less than a certain value (state s3), slabs can move as a bedload;otherwise, if the slope is greater than or equal to a
threshold value (state s4), slabs are in a perfectly
shielded area (shadow zone), in which the flow is not
sufficient to convey sand as a bedload, and thus, thesand moves only by avalanching. From these assumptions, we obtain a simple formulation for the statedependent (s1 to s4) probability, Bi, j, that the sand atsite (i, j) will migrate by bedload transport to thenearest downstream cell:Endo Progress in Earth and Planetary Science (2016) 3:28

8
α
>
>
<
α= ðhði þ 1; jÞ−hði; jÞÞBi; j ¼
> γ fβ− ðhði−1; jÞ−hði; jÞÞg>:0Page 4 of 11

on a local peak; including a flat top ðs1Þ
midway along a stoss side ðs2Þ
on a lee slope; but not in a shadow zone ðs3Þin a shadow zone ðs4Þ;

where α, β, and γ are positive constants. The second expression on the right-hand side shows that the difficultyof climbing a stoss slope depends on the gradient. The
third expression on the right-hand side shows the extent
to which the bedform reduces the flow, as a function of
the gradient of the lee slope.
In Eq. (2), the condition of the local state (s1 to s4) isdetermined as follows. If hi − 1,j ≦ hi,j, the site (i, j) is on astoss slope or at a local peak. If it is also true that hi + 1,j− hi,j ≦ 0, then the site is a peak or a flat top (s1); otherwise, it is somewhere along a stoss slope (s2), and thevelocity is inversely proportional to the gradient. On theother hand, if hi − 1,j > hi,j, the site (i, j) is on a lee slope.
If it is also true that hi − 1,j − hi,j ≦ β, then site (i, j) is notin a shadow zone (s3); otherwise, it is in a shadow zone(s4). Note that β is the threshold value for the shadow
zone. The above determination of the site state for the
bedload (s1 to s4) can be summarized as follows:8(
if hiþ1;j −hi;j ≦ 0>>>> if hi−1;j ≦ hi;j and<if hiþ1;j −hi;j > 0(>
if hi−1;j − hi;j ≦ β>>> if hi−1;j >hi;j and:if hi−1;j − hi;j >β: ðs1Þ

fdownward slopeg > Ar with the condition ðs4Þ;

ð2Þ

ð5Þ

where Ar is a positive constant that represents the angleof repose, and the slope in the downstream (positive xdirection) is defined as
fdownward slopeg ¼ hði ; jÞ−hði þ 1; jÞ:ð6Þ

Unlike with bedload transport, because of gravity, avalanching can also occur in the lateral direction:(
hði; j þ 1Þ→hði; j þ 1Þ þ qa
hði; jÞ→hði; jÞ−qa andor hði; j−1Þ→hði; j−1Þ þ qa g:
ð7ÞNote that oblique movements between adjacent cells
are ignored in the present simulation. Lateral avalanching occurs under the same condition as in Eq. (5), buthere, the slopes of interest are those in the positive andnegative y-directions, which are defined as: ðs2Þ
:: ðs3Þ

fdownward slopeg ¼ hði ; jÞ−hði; j þ 1Þ; hði ; jÞ−hði; j−1Þ;: ðs4Þ

respectively. If the lateral slopes in both the positive andnegative y-directions are steeper than the angle of repose, only one of these two directions is selected, andthey have equal probability:

ð3Þ

Formulation of avalanching

Avalanching occurs as a result of gravity, and consequently, sand can be conveyed by avalanching in ashadow zone, even when bedload transport does notoccur. Avalanching occurs when the lee slope is steeper
than the angle of repose; the probability that avalanchingoccurs is unity, which implies deterministic and immediate movement. For avalanching in the streamwise direction at site (i, j), the topographic height at (i, j) develops
as follows:hði; jÞ→hði; jÞ−qa and hði þ 1; jÞ→hði þ 1; jÞ þ qa ;

ð4Þ

where qa is the number of slabs that move by avalanching during one time step. According to flume experiments, avalanching does not occur on a stoss side; thus,
we do not need to consider avalanching in the negative
x-direction. Avalanching always occurs if the following
condition is satisfied:ð8Þ

1
Alþ ¼ Al− ¼ ;2

ð9Þ

where Al + and Al − are the probabilities of avalanchingin the positive and negative y-direction, respectively.
In previous CA models for aeolian dunes, avalanching
sand moves in the steepest direction (e.g., Werner
1995). In the present model, for simplicity, we adopted
the assumption stated in Eq. (9). When averaged overtime, this rule is equivalent to the distribution due toavalanching of sand in a sandpile (e.g., Bak et al.
1987); that is, the sand is distributed equally in all
possible directions.
Simulation

The simulation in the present study was executed as follows. At each time step, the movement of sand slabs inthe x-direction, due to both bedload transport and avalanching, was determined for each cell and updated inEndo Progress in Earth and Planetary Science (2016) 3:28

parallel. Next, lateral avalanching (in the y-direction)
was determined for each cell and updated in parallel.
During a single step, avalanching occurred only once percell (i.e., it was not repeated until perfect relaxation).
This was done for simplicity and because we consider
bedload transport, and thus, it is not necessary to obtaina static slope (in a strict sense) after perfect relaxation.
In selecting values for the phenomenological parameters, we made the following assumptions. The topographic height was taken to be the number of sandslabs, and therefore, β and Ar were positive integers. Because the slope of a stoss side of a bedform is always
gentler than the corresponding lee slope, the inequality
1 < β should be satisfied. In addition, we impose the inequality α > γ(β − 1) because the migration of grains isfaster on a stoss side than it is on the corresponding lee
slope (see Eq. (2) for the gentlest lee slope). The value of
α is selected arbitrarily, subject to the condition that it
be sufficiently smaller than unity; this probability is as asurrogate for the velocity of sand particles. The simulation results are insensitive to α, apart from the rate ofdevelopment for a given time step. For a given value ofα, the values of γ should be sufficiently small that theabove inequalities are satisfied; otherwise, the model will
not adequately reproduce the shielding effect of lee
slopes. For simplicity and to take account of the above
inequalities, we used the following values for the parameters: α = 0.05, β = 2, γ = α/8(=6.25 × 10− 3), and Ar = 4.
We let qb = 1, which represents the assumption that the
bedload is a thin layer. We let qa = 2 so that there was afast relaxation of the avalanching; however, we note thatwhen qb = 1, the simulation results were very similar,
providing α was sufficiently smaller than unity. The initial state of the bed was almost flat, with some small
random variation, for the available sediment (sand slabs),
while the substrate was perfectly flat. The calculation
space was a 200 × 40 grid of cells, although we also used
a 400 × 80 grid of cells to make particular comparisons.

Page 5 of 11

Calculation runs were performed for more than 105 time
steps, using a periodic boundary condition. Because thesimulation was a minimal model and used the simplest
formulae and boundary condition, the main discussion
will be restricted to the case of sparse available sediments. We note that when sediments are plentiful, theinterference between opposing boundaries cannot be ignored, and thus, it is not appropriate to use a periodic
boundary condition.

Results and discussion
Development of sand topographic features

First, to evaluate the present model, we qualitatively
assessed the dependence of the topographic features on
the amount of sand, as expressed by the average sandthickness, which is defined as the total number of sandslabs divided by the number of cells in the system. In
our simulation, the following types of bed appeared (inorder of increasing amount of sand): sheet flows, sandpatches, isolated (barchan) ripples, and continuous ripples. When the average sand thickness was less than 0.6,
bedforms (ripples) were not generated, and sand migrated only as a veneer of fast-moving sand particles,
which is referred to as sheet flow (Bagnold 1941; Cookeet al. 1993; Venditti et al. 2005). When the average sandthickness was about 0.65 (0.6–0.7), sand patches (small
bunches of sand; Bagnold 1941; Cooke et al. 1993) appeared and coexisted with sheet flows, but no stable ripples occurred. Sand patches were repeatedly generated
and destroyed, returning entirely to sheet flows (Fig. 1).
Such a situation is commonly seen in flume experiments, and this contrasted to the case of ripples, which,once formed, were robust and maintained their shapeduring migration. Ripples were generated only when the
average sand thickness was greater than 0.7. When the
average sand thickness was less than 22, the ripples werebarchans, an isolated type of ripple that typically takes acrescent shape (Fig. 2a); this is consistent with theFig. 1 Sand patches repeatedly appearing and disappearing. An identical run (average sand thickness of 0.65) is shown at different times: a 35,000
steps, b 50,000 steps; the sand moves from the left to right. A sand patch (enclosed by the oval in a) vanished 15,000 steps later. Note that there is
only a little sand in the area downstream of the sand patch, which suggests that the sand patch gathers sand from upstream. Nevertheless, under thiscondition, the sand patch does not grow to an isolated ripple. The height of one slab is one-tenth the length of one cellEndo Progress in Earth and Planetary Science (2016) 3:28

Page 6 of 11

Fig. 2 Calculated topographies from the oblique view. a Isolated ripples (barchan ripples) and b continuous ripples; shown for the same time(105 steps), with average sand thicknesses of 1.0 and 30, respectively. The movement of particles is from the left to right, and in the graphics, thelight is cast from the upper left. The height of one slab is one-tenth the length of one cellfinding that barchans form when the amount of sand is
not sufficient to cover the entire area, as found by fieldstudies (Bagnold 1941; Allen 1968; Wasson and Hyde
1983; Cooke et al. 1993) and flume experiments (Mantz
1992; Niño and Barahona 1997; Kleinhans et al. 2002;
Endo et al. 2004, 2005; Dreano et al. 2010). As is wellknown, if there is abundant sand, a bed will have continuous ripples (e.g., Bagnold 1941; Allen 1968; Wasson
and Hyde 1983; Dreano et al. 2010). In our simulation,
when the average sand thickness was more than 22,
continuous ripples were generated (Fig. 2b). Thus, thepresent simple model qualitatively reproduces the development of sand topographies that are due to the size ofthe bedload. Below, we will consider the starved state(beds with sparse sediments).
Temporal changes in sand flux compared with flume
experiments

The sand flux is defined as the sand discharge per width(note that this differs from the general definition of flux
in terms of the unit cross section). The sand discharge
does not distinguish between that due to flow motions
and that due to avalanching, and therefore, the sand flux
includes both bedload transport and avalanching. If thesand discharge is quantified by bulk volume (including
both sand and intergrain voids) per unit time, the dimensions of sand flux are the square of length over time
(L2/T), because the volume per width is measured inunits of L2.
It is known that when a bedform develops from a flat
bed, the migration velocity decreases, but the size of the
bedform increases asymptotically (Baas 1994; Coleman
et al. 2003; Perillo et al. 2014). However, the actual timedependence of the sand flux has rarely been measured;
this is probably due to the difficulty of doing so. As one
example, Betat et al. (2002) used statistical methods tomeasure the sediment transport rate, the average height,
and the migration velocity of irregular bed states during
the development of bedforms. According to a graphpresented in Betat et al. (2002), the sediment transport
rate decreased monotonically with time, except at thevery start of the experiment. Although this showed thetemporal change in the sand flux, the definition of thesediment transport rate in that paper (Betat et al. 2002)
was different from the usual definition of sand flux. For
a comparison using the usual definition, the measured
values were transformed to the sand flux using the simple relation Q ∝ H × V, where Q is the sand flux, H is theheight of the bedform, and V is the migration speed. If
the sand discharge is measured by bulk volume insteadof weight, Q = H × V (Dreano et al. 2010). Although actual beds are irregular, especially during the development of bedforms (before reaching equilibrium), thisrelation can be applied to the data of Betat et al. (2002),
because these are spatial averages for each moment andare obtained through precise measurements and calculations; therefore, they are suitable for estimating the overall trend of the sand flux.
Figure 3 shows the temporal changes in the reconstructed sand flux and the bedform developments foundin experiments (Betat et al. 2002). The sand flux decreased monotonically as it transitioned from a flat bed toa stable state; however, note that such a decrease in thesand flux is not generally observed, because natural locations are generally regarded to be in a stable state relative
to the time scale of the observation. There has been little
discussion of the decrease in sand flux prior to a stablestate that is seen in mathematical models.
Figure 4 shows an example of the temporal changes
in sand flux, as seen in the present model; note thatthe flux is spatially averaged over momentary values,
and the standard deviation of the topographic heightis used as an index of bed development. The simulation showed that the flux decreased monotonically as
the bedform developed and enlarged. This change
qualitatively agrees with the results of the experiment(Fig. 3), which had not been explicitly reproduced inprevious models.

Endo Progress in Earth and Planetary Science (2016) 3:28

Fig. 3 Experiment: temporal development of velocity, height, andreconstructed sand flux of aqueous bedforms. The original data aretaken from Betat et al. (2002). The sand flux is calculated using
Bagnold’s relation with a constant of unity, velocity × height, toshow the relative change

Calculated sand flux in a stable state and its dependence
on the amount of sand

In Fig. 5, the long-term stable state of the sand flux, as
shown in Fig. 4, is compared to the average sand thickness (up to 3.5). For the range of average sand thicknesses shown in Fig. 5, the bed state changes from sandsheets through sand patches to isolated (barchan) ripples. Figure 5 shows that the sand flux increases andthen decreases as the average sand thickness increases.
The stage of increasing flux seems to directly reflect anincrease in the amount of sand in the system, which
likely corresponds to the supply-limited transport mode.
The stage of decreasing flux in Fig. 5, however, cannot
be explained by a supply limitation nor by a capacity
limitation, because the model in this study does not
explicitly calculate the fluid motion or the capacity ofPage 7 of 11

sediment transport due to flow. Additionally, the rangeof average sand thickness considered here (0–3.5) is a
supply limitation in the ordinary sense of geomorphology
and sedimentology, because the amount of sand is not
sufficient to cover the entire substrate.
Figure 5 shows that the sand flux in a stable stateincreases and then decreases as the average sandthickness increases and the bed state changes fromsand sheet through sand patches to isolated ripples.
Additionally, in the cases in which we quadrupled thegrid (400 × 80 cells), the total trend was identical, although the position and the value of the peak (ca.
5.5 × 10−1 and 8.6 × 10−3, respectively) were slightly
different (ca. 0.5 × 10−1 and 0.4 × 10−3 smaller, respectively). In the present model, Eq. (2), especially in thecase of s2, indicates that the movement of a sand slab
is suppressed by the topographic state of the adjacentdownstream cell. In a specific test case in which condition (s2) was not included (sand migrates with constant probability regardless of the slope), the bedforms
that formed did not have a plausible shape (stoss sides
were steeper than lee slopes). This indicates the importance of interference by nearby particles in the development of beds; note that this works negatively interms of sand flux.
Comparison with traffic flow and implications

As stated in the previous section, interaction between
sand in neighboring cells (Eq. (2)) is necessary to reproduce bedforms with plausible shapes. Dynamic interaction
between grains is a theme of physics of multiparticle systems or granular dynamics, and it is generally complicated. For specific topics, however, simplified models
based on granular dynamics can reproduce real phenomena; an example of this is traffic flow. Although theFig. 4 Temporal change of sand flux and standard deviation (std) of the topographic height. The standard deviation of height is an index of thedevelopment of the bed, which increases as ripples enlarge; it can be easily calculated for irregular beds and isolated bedforms. Average sandthickness is 1.0. The height of the sand is equal to the number of sand slabs in each cellEndo Progress in Earth and Planetary Science (2016) 3:28

Page 8 of 11

Fig. 5 Simulation: long-term sand flux against the average sand thickness. As shown in Fig. 4, the sand flux asymptotically decreases until thetime step is approximately 104, and then it becomes almost constant with only a small fluctuation. The bars indicate the standard deviation of
the fluctuation (the length of each bar is 2σ), which is obtained from data sampled every 500 steps from 1 × 105 to 5 × 105 steps for 10 runs
under the same conditions (except for small random differences in the initially almost flat bed). The sand thickness is equal to the number ofcompiled slabs. The sand flux is the number of transported sand slabs per unit width per unit time steppresent model for starved beds is built phenomenologically on the basis of the continuous bed-contact load, it has
similarities to a model for traffic flows. Traffic flow is ofinterest both in the fields of physics and in civil engineering, where cars are regarded as grains and their flow considered a type of granular flow (e.g., Chowdhury et al.
2000). In such studies, the relation between car-flow
phases (i.e., “grain-flow phases”) and the abundance ofcars (density of grains) has been investigated. The grainflow phase is a class or state of transportation in a traffic
flow, and the traffic viewpoint is different from the transport modes of geomorphology and sedimentology (suchas rolling, sliding, saltation, and suspension). Previous
studies of traffic flow have shown that car-flow phases areconnected with car transportation rates (Chowdhury et al.
2000). Car transportation rates correspond to sand flux
for sediment beds, which is a quantifiable factor. In geomorphology and sedimentology, while interconnectional
properties among different levels of interaction (for example, grain−fluid, bedform−fluid, and bedform−bedform)
have been identified (e.g., Kocurek et al. 2010), differencesin the grain-flow phase have not been clearly identified.
Here, we review a previous model of traffic flow and compare it with the present model in order to understand therelation between sand flux and the state of sand transportation on starved beds.
The stochastic optimal velocity (SOV) model (Kanai
et al. 2005) is a stochastic CA model of the traffic flow.
The simple SOV model considers the case in which aparticle (car) migrates at most one cell per step, as in

the present model for bedload transport. Before discussing the similarities between the present model and the
SOV model, we will consider the difference. The simplest SOV model assumes a single-lane road and prohibits overtaking, whereas the present model for bedloadtransport allows overtaking. The present model accepts
more than one particle (slab) per cell, which represents acompilation of sand, while the SOV model allows only
one particle per cell.
A flux–density diagram is used to show the dynamic
properties of car traffic flows. Because the average sandthickness in the present study is equivalent to the densityof particles (the number of particles per unit area), Fig. 5
can be compared to a flux–density diagram of traffic flow.
Despite the different research targets, the two models havecommon points. In both models, the probability of migration of particles corresponds to the velocity. The SOVmodel reproduces the relation between the flux and thedensity of particles, and it is asymmetric (fast increase andlong-tailed decrease). This agrees with observations of areal highway (Sugiyama et al. 2008) and resembles the results of the present model (Fig. 5), although car traffic
jams always propagate upstream. Interestingly, the curvein Fig. 5 shows an abrupt decrease around the peak before
a tail. A similar abrupt decrease is seen in observations ofreal traffic and is interpreted as an overshoot of the freephase against the jam phase, showing a metastability bifurcation (Kanai et al. 2005).
When comparing the model of the present study andthat of Kanai et al. (2005), we note that the stage ofEndo Progress in Earth and Planetary Science (2016) 3:28

increasing sand flux (average sand thickness of less than0.6; sheet flows occur; Fig. 5) corresponds to the freeflow phase of traffic flows. The concepts of supplylimited sediment transport and free-flow phase are inaccordance, in that they are free from interference fromother particles. The state of decreasing sand flux (average sand thickness over 0.7; isolated ripples are generated; Fig. 5) corresponds to the jam phase. The vicinityof the peak sand flux (average sand thickness 0.6–0.7;
sand patches repeatedly appear and disappear, whilesand sheets exist continuously; Fig. 5) is the transition
between the above two phases and likely corresponds tothe congested phase. Therefore, it is suggested that sandtopographic features on starved beds are dominated by
the interference between sand particles.
Limitations of the present model and future
developments for potential applications

It is important to note that the present minimal model
can be applied only to specific cases in which the bedloadtransport rate is very low. For wider application, formulae
should include other processes, such as erosional effects
occurring downstream of lee slopes, three-dimensional effects of flows, and a more realistic avalanching process.
These effects become important when bedform heights
are large; therefore, the use of the present minimal model
requires restricting the size of the bedload and bedforms.
It is interesting to consider a flume experiment thatwas conducted under conditions similar to those discussed above. Dreano et al. (2010) showed that in thestable state, the bedload flux due to bedform migration
increased from barchan ripples to transverse barchanoid
ripples if the sand supply was increased (however, wenote that the setup of these experiments was different
from the present simulation). Dreano et al. (2010) observed dynamic equilibrium states in which there was awell-controlled forced influx of sand into the system,
while the present simulation assumes a closed system. In
our simulation, we did not observe development of isolated ripples as the abundance of sand increased; instead,
we investigated the relation between sand flux and thetransition between sand sheets, sand patches, and isolated ripples, which had not been investigated sufficiently, compared with the properties of bedforms(ripples and dunes). In the personal experience of theauthor, during preliminary runs and at early stages of actual runs of flume experiments for different purposes, it
is often observed that sand patches occur with thin sandsheets under conditions of sparsely available sediments.
The sand patches migrate more slowly than do sandsheets, and this probably decreases the sand flux, although no measurements have been conducted (notethat precise measurement of this is difficult). During the
growth of sand patches, sand particles transported fromPage 9 of 11

upstream are incorporated into sand patches and decelerated (like cars on a congested road); this local state is
not in equilibrium, but in a stochastic (statistical) sense,
it does not violate the equilibrium of the entire system,
because sand patches repeatedly appear and disappear inmany places. Such behavior is compatible with the stochastic model. Although, the dynamics in the background of the stochastic process and the connection of
the numerical model to actual beds remain to be further
explored in future studies using new techniques formeasurement and observation, considering the change
in sand flux due to bed states through a grain-flow phase
is a viewpoint that leads to a deeper understanding of
bedforms and other sand topographies.
There have been jam-like phenomena reported, although those are more complicated than the present
model. When sediments pass a meander bend of a river
channel with a pool, the migration speed of sediments
decreases greatly as they enter the pool, and it thenincreases as they leave (Kleinhans et al. 2002); thisresembles the behavior of cars entering and leaving traffic jams. Another example is the pulsatory migration ofbimodal mixtures of grain-sized sediment. A mixture oflarge grains (gravel) and finer grains repeats the clustering and scattering behavior when they are transported
downstream as a bedload (Iseya and Ikeda 1987). It isconsidered that this behavior comes directly from theinterference caused by the presence of other grains (i.e.,jams), rather than from the water flow due to a moving
boundary condition that is induced by the other grains.
Such a situation is very complicated in the sense of thephysics of multiparticle systems, but it could be at least
partially explained by modeling it as a jam of transported
particles; however, it would be necessary to add some
new ideas to the present model.
Partially alluvial covers in bedrock river channels play
an important role in the incision rate. Those can either
promote incision by providing tools or inhibit incisionby acting as armor. Finnegan et al. (2007) performed
flume experiments with a straight channel and homogeneous sediments; they demonstrated that following sediment supply perturbations, the growth of roughness ofthe bedrock, which is caused by the deposition of alluvial
cover, resulted in both negative and positive feedback to
the average incision rate. Gravel in a channel, which can
act as armor, can help preserve potential topographies;
for example, inverted paleochannels on the floor of
Miyamoto crater (Mars) were able to resist later (modern) aeolian erosion, as compared to the area outside thechannel (Newsom et al. 2010).
When a bedrock channel is considerably rough, sparsely
distributed sediments move in complicated ways. There isa need for a numerical model that can predict, for any
shape of channel, the spatial and temporal variation of theEndo Progress in Earth and Planetary Science (2016) 3:28

incision rates induced by sediment supply perturbations;beyond that, there is a need for a model that can do thisfor bimodal sediments. Although the present study considered only flat, nonerodible substrates, it may be possible to adapt it so that it can incorporate changes inroughness of the substrate; this could be done by formulating phenomenological equations for erosion due topartial alluvial covers.

Conclusions
A simple CA model simulating continuous bed-contact
loads (rolling and sliding) is proposed; this model wasdeveloped phenomenologically and without the calculation of fluid motion. It is intended primarily for starved
beds in a subaqueous system. As the amount of sand increases, the model reproduces the following bed types inorder: sheet flows, sand patches, and ripples. While sandpatches coexist with sand sheets and are not sustainable,
ripples are robust once formed.
In the simulation, during bed development, the sand flux
showed a monotonic decrease. This is in qualitative agreement with data from a previous well-controlled flume
experiment that had not previously been reproduced by anumerical model.
In the stable state, the sand flux versus average sandthickness (flux–density relation) shows a steep increase
followed by a tailed decrease, which is similar to the pattern seen in car traffic. The similarities between sandbeds and traffic flow suggest that the changes in sandflux due to the bed state and the interference between
sediments are important factors in the generation of
bedforms, especially for the transition between topographies on starved beds. In the simulation, the change insand topographic features corresponds well to theflux–density diagram. Sand topographic features shift
at the fluxes of the boundaries between segments ofthe curve in the diagram (Fig. 5), which suggests thatsand sheets, sand patches, and sand ripples correspondto the traffic phases of free flow, congested, and jam,respectively, depending on the degree of interference
between sand particles.
Abbreviationsh(i, j): Topographic height at position (i, j); Bi, j: Probability of sand migration
by bedload transport, depending on the state at position (i, j); qb andqa: Numbers of slabs migrating by bedload transport and by avalanching,respectively, during a single calculation step; Sc: Critical slope of stoss side;
Ar: Positive constant representing the angle of repose; β: Critical value for the
shadow zone; α and γ: Positive constants that are parameters for Bi, j; Al+ and
Al–: Probabilities of avalanching in the positive and negative y-directions,
respectively, in the case where both directions are possible
Acknowledgements
Special thanks are extended to Prof. F. Masuda of Doshisha Univ. for many
hours of helpful discussion. Additionally, I would like to thank Dr. A. Katsuki
for advice on the early version of the manuscript.

Page 10 of 11

Funding
This study was financially supported by JSPS KAKENHI Grant Number 25350424.
Authors’ contributions
NE performed calculations and drafted the manuscript.
Competing interests
The author declares that he has no competing interests.
Received: 3 December 2015 Accepted: 18 August 2016

References
Allen JRL (1968) Current ripples: their relation to patterns of water andsedimentation motion. North-Holland, New York
Allen JRL (1982) Sediment structures: their character and physical basis, vol 1.
Elsevier, Amsterdam
Andreotti B, Claudin P, Devauchelle O, Durán O, Fourrière A (2012) Bedforms
in a turbulent stream: ripples, chevrons and antidunes. J Fluid Mech 690:
94–128Baas JH (1994) A flume study on the development and equilibrium morphology
of current ripples in very fine sand. Sedimentology 41:185–209
Bagnold RA (1941) The physics of blown sand and desert dunes. Methuen, London
Bak P, Tang C, Wiesenfeld K (1987) Self-organized criticality: an explanation of the
1/f noise. Phys Rev Lett 59:381–384
Barchyn TE, Hugenholtz CH (2012) A new tool for modeling dune field evolution
based on an accessible, GUI version of the Werner dune model. Geomorphology
138:415–419
Barrio-Parra F, Rodríguez-Santalla I (2014) A free cellular model of dune dynamics:
application to El Fangar spit dune system (Ebro Delta, Spain). Comput Geosci
62:187–197
Betat A, Kruelle C, Frette V, Rehberg I (2002) Long-time behavior of sand ripples
induced by water shear flow. Eur Phys J E: Soft Matter Biol Phys 8:465–476
Bridge JS, Demicco R (2008) Earth surface processes, landforms and sedimentdeposits. Cambridge University Press, Cambridge
Chowdhury D, Santen L, Schadschneider A (2000) Statistical physics of vehiculartraffic and some related systems. Phys Rep 329:199–329
Coleman SE, Fedele JJ, Garcıa MH (2003) Closed conduit bed-form initiation anddevelopment. J Hydraul Eng 129:956–965
Cooke RU, Warren A, Goudie AS (1993) Desert geomorphology. Taylor & Francis,
LondonDiniega S, Glasner K, Byrne S (2010) Long-time evolution of models of aeolian
sand dune fields: influence of dune formation and collision. Geomorphology
121:55–68
Dreano J, Valance A, Lague D, Cassar C (2010) Experimental study on transientand steady-state dynamics of bedforms in supply limited configuration. Earth
Surf Proc Land 35:1730–1743
Duran O, Claudin P, Andreotti B (2014) Direct numerical simulations of aeolian
sand ripples. Proc Natl Acad Sci U S A 111:15665–15668
Endo N, Taniguchi K, Katsuki A (2004) Observation of the whole process of interaction
between barchans by flume experiments. Geophys Res Lett 31:L12503
Endo N, Sunamura T, Takimoto H (2005) Barchan ripples under unidirectional
water flows in the laboratory: formation and planar morphology. Earth Surf
Proc Land 30:1675–1682
Finnegan NJ, Sklar LS, Fuller TK (2007) Interplay of sediment supply, river incision,and channel morphology revealed by the transient evolution of anexperimental bedrock channel. J Geophys Res 112:F03S11
Iseya F, Ikeda H (1987) Pulsations in bedload transport rates induced by alongitudinal sediment sorting: a flume study using sand and gravel
mixtures. Geogr Ann A 69:15–27
Kanai M, Nishinari K, Tokihiro T (2005) Stochastic optimal velocity model and its
long-lived metastability. Phys Rev E 72:035102
Katsuki A, Kikuchi M, Endo N (2005) Emergence of a barchan belt in a unidirectional
flow: experiment and numerical simulation. J Phys Soc Jpn 74:878–881
Kleinhans MG, Wilbers AWE, De Swaaf A, Van Den Berg JH (2002) Sedimentsupply-limited bedforms in sand–gravel bed rivers. J Sediment Res 72:
629–640
Kocurek G, Ewing RC, Mohrig D (2010) How do bedform patterns arise? New
views on the role of bedform interactions within a set of boundary conditions.
Earth Surf Proc Land 35:51–63

Endo Progress in Earth and Planetary Science (2016) 3:28

Page 11 of 11

Mantz PA (1992) Cohesionless, fine-sediment bed forms in shallow flows.
J Hydraul Eng 118:743–764
Melo HP, Parteli EJ, Andrade JS, Herrmann HJ (2012) Linear stability analysis oftransverse dunes. Physica A 391:4606–4614
Miyata Y, Iwasa S (1997) Ripple initiation by random grain-rolling: water-tank
experiment and numerical simulation. In: Abstracts of the 104th annual
meeting of the ecological society of Japan., p 319
Momiji H, Carretero-Gonzalez R, Bishop SR, Warren A (2000) Simulation of theeffect of wind speedup in the formation of transverse dune fields. Earth Surf
Proc Land 25:905–918
Newsom HE, Lanza NL, Ollila AM, Wiseman SM, Roush TL, Marzo AG, Tornabene
LL, Okubo CH, Osterloo MM, Hamilton VE, Crumpler LS (2010) Inverted
channel deposits on the floor of Miyamoto crater, Mars. Icarus 205:64–72
Niiya H, Awazu A, Nishimori H (2010) Three-dimensional dune skeleton model asa coupled dynamical system of two-dimensional cross sections. J Phys Soc
Jpn 79:063002
Niño YB, Barahona M (1997) Barchan-ripples: emergence, evolution and flow–
sediment interactions. Proc Int Assoc Hydraul Res 27B:1037–1042
Nishimori H, Ouchi N (1993) Formation of ripple patterns and dunes by windblown sand. Phys Rev Lett 71:197
Parteli EJR, Herrmann HJ (2003) A simple model for a transverse dune field.
Physica A 327:554–562
Perillo MM, Best JL, Yokokawa M, Sekiguchi T, Takagawa T, Garcia MH (2014) A
unified model for bedform development and equilibrium under unidirectional,
oscillatory and combined-flows. Sedimentology 61:2063–2085
Reid I, Frostick LE (1994) Fluvial sediment transport and deposition. In: Pye K (ed)Sediment transport and depositional processes. Blackwell Scientific, Oxford,
pp 89–155
Rozier O, Narteau C (2014) A real-space cellular automaton laboratory. Earth Surf
Proc Land 39:98–109
Schwämmle V, Herrmann HJ (2003) Geomorphology: solitary wave behaviour ofsand dunes: colliding dunes appear to traverse through one another andemerge unscathed. Nature 426:619–620
Sugiyama Y, Fukui M, Kikuchi M, Hasebe K, Nakayama A, Nishinari K, Tadaki S,
Yukawa S (2008) Traffic jams without bottlenecks: experimental evidence forthe physical mechanism of the formation of a jam. New J Phys 10:033001
Venditti JG, Church M, Bennett SJ (2005) Morphodynamics of small-scale
superimposed sand waves over migrating dune bed forms. Water Resour
Res 41:W10423
Wasson R, Hyde R (1983) Factors determining desert dune type. Nature 304:337–339
Werner B (1995) Eolian dunes: computer simulations and attractor interpretation.
Geology 23:1107–1110
Werner BT, Kocurek G (1999) Bedform spacing from defect dynamics.
Geology 27:727–730
Wippermann F, Gross G (1986) The wind-induced shaping and migration of an
isolated dune: a numerical experiment. Bound-Lay Meteorol 36:319–334

Submit your manuscript to ajournal and beneﬁt from:
7 Convenient online submission
7 Rigorous peer review
7 Immediate publication on acceptance
7 Open access: articles freely available online
7 High visibility within the ﬁeld
7 Retaining the copyright to your article

Submit your next manuscript at 7 springeropen.com