**Research article**
17 Jul 2018

**Research article** | 17 Jul 2018

# A lattice grain model of hillslope evolution

Gregory E. Tucker Scott W. McCoy and Daniel E. J. Hobley

^{1},

^{2},

^{3}

**Gregory E. Tucker et al.**Gregory E. Tucker Scott W. McCoy and Daniel E. J. Hobley

^{1},

^{2},

^{3}

^{1}Cooperative Institute for Research in Environmental Science (CIRES) and Department of Geological Sciences, University of Colorado, Boulder, CO 80305, USA^{2}Department of Geological Sciences and Engineering, University of Nevada, Reno, NV 89557, USA^{3}School of Earth and Ocean Sciences, Cardiff University, Cardiff, CF10 3AT, Wales

^{1}Cooperative Institute for Research in Environmental Science (CIRES) and Department of Geological Sciences, University of Colorado, Boulder, CO 80305, USA^{2}Department of Geological Sciences and Engineering, University of Nevada, Reno, NV 89557, USA^{3}School of Earth and Ocean Sciences, Cardiff University, Cardiff, CF10 3AT, Wales

**Correspondence**: Greg Tucker (gtucker@colorado.edu)

**Correspondence**: Greg Tucker (gtucker@colorado.edu)

Received: 20 Jan 2018 – Discussion started: 07 Feb 2018 – Revised: 21 May 2018 – Accepted: 22 Jun 2018 – Published: 17 Jul 2018

This paper describes and explores a new continuous-time stochastic cellular automaton model of hillslope evolution. The Grain Hill model provides a computational framework with which to study slope forms that arise from stochastic disturbance and rock weathering events. The model operates on a hexagonal lattice, with cell states representing fluid, rock, and grain aggregates that are either stationary or in a state of motion in one of the six cardinal lattice directions. Cells representing near-surface soil material undergo stochastic disturbance events, in which initially stationary material is put into motion. Net downslope transport emerges from the greater likelihood for disturbed material to move downhill than to move uphill. Cells representing rock undergo stochastic weathering events in which the rock is converted into regolith. The model can reproduce a range of common slope forms, from fully soil mantled to rocky or partially mantled, and from convex-upward to planar shapes. An optional additional state represents large blocks that cannot be displaced upward by disturbance events. With the addition of this state, the model captures the morphology of hogbacks, scarps, and similar features. In its simplest form, the model has only three process parameters, which represent disturbance frequency, characteristic disturbance depth, and base-level lowering rate, respectively. Incorporating physical weathering of rock adds one additional parameter, representing the characteristic rock weathering rate. These parameters are not arbitrary but rather have a direct link with corresponding parameters in continuum theory. Comparison between observed and modeled slope forms demonstrates that the model can reproduce both the shape and scale of real hillslope profiles. Model experiments highlight the importance of regolith cover fraction in governing both the downslope mass transport rate and the rate of physical weathering. Equilibrium rocky hillslope profiles are possible even when the rate of base-level lowering exceeds the nominal bare-rock weathering rate, because increases in both slope gradient and roughness can allow for rock weathering rates that are greater than the flat-surface maximum. Examples of transient relaxation of steep, rocky slopes predict the formation of a regolith-mantled pediment that migrates headward through time while maintaining a sharp slope break.

- Article
(21871 KB) -
Supplement
(102131 KB) - BibTeX
- EndNote

Hillslopes take on a rich variety of forms. Their profile shapes may be convex-upward, concave-upward, planar, or some combination of these. Some slopes are completely mantled with soil, whereas others are bare rock, and still others draped in a discontinuous layer of mobile regolith. The processes understood to be responsible for shaping them are equally varied, ranging from disturbance-driven creep to dissolution to large-scale mass movement events.

Considerable research has been devoted to understanding the evolution of soil-mantled slopes that are primarily governed by disturbance-driven creep, such as downslope soil transport by biotic and abiotic soil-mixing processes. As a result, the geomorphology community has mathematical models that account well for observed slope forms and patterns of regolith thickness (Roering, 2008). Furthermore, stochastic-transport theory provides a mechanistic link between the statistics of particle motion, the resultant average rates of downslope transport, and the emergence of convex-upward, soil-mantled slope forms (Culling, 1963; Foufoula-Georgiou et al., 2010; Furbish and Haff, 2010; Furbish et al., 2009; Roering, 2004; Tucker and Bradley, 2010).

One gap that remains, however, lies in understanding steep, rocky slopes (Fig. 1). “Rocky” implies slopes that lack a continuous soil cover (Howard and Selby, 1994); here, a transport law that assumes the existence of such a cover no longer applies. “Steep” implies angles approaching or exceeding the effective angle of repose for loose, granular material, so that ravel may be an important transport mode (Gabet, 2003; Gabet and Mendoza, 2012; Lamb et al., 2011; Roering and Gerber, 2005) and particles have the potential to fall as soon as they are released from bedrock. This type of relatively fast, long-distance transport does not fit comfortably in the framework of standard diffusion-based models of hillslope soil transport, which derive from an underlying assumption that the characteristic length scale of motion is short relative to the length of the slope.

Rocky slopes are rarely completely barren. More commonly, they have a patchy cover of loose material, which may either retard rock weathering by shielding the rock surface from moisture or temperature fluctuations, or enhance it by trapping water and allowing limited plant growth. A discontinuous cover does not fit easily within the popular exponential-decay regolith-production models (Heimsath et al., 2012; Lamb et al., 2013), which assume an essentially continuous soil mantle.

An additional issue, which pertains to both rocky and soil-mantled slopes, is the connection between sediment movement at the scale of individual “motion events,” and the resulting longer-term average sediment flux, which forms the basis for continuum models of hillslope evolution. Recent theoretical and experiment work has begun to forge a mechanistic connection between these scales (Culling, 1963, 1965; Furbish and Haff, 2010; Furbish et al., 2009; Gabet and Mendoza, 2012; Lamb et al., 2013; Tucker and Bradley, 2010). However, the community's resources for computational analysis of particle-level dynamics remain limited, lagging behind developments in understanding sediment transport in coastal environments (Drake and Calantoni, 2001) and rivers (Furbish and Schmeeckle, 2013; MacVicar et al., 2006; McEwan and Heald, 2001; Schmeeckle, 2014).

To further our understanding of how grain-level weathering and transport processes translate into hillslope evolution, both for hillslopes in general and rocky slopes in particular, it would be useful to have a computational framework with which to conduct experiments. Ideally, such a framework should be sophisticated enough to capture the essence of weathering and granular mechanics, while remaining simple enough to involve only a small number of parameters and provide reasonable computational efficiency.

Our aim in this paper is to describe one such computational framework, test whether it is capable of reproducing commonly observed hillslope-profile forms, and examine how its parameters relate to the bulk-behavior parameters used in conventional continuum models of soil creep and regolith production. The model uses a pairwise, continuous-time stochastic (CTS) approach to combine a lattice grain model with rules for stochastic bedrock-to-regolith conversion (“weathering”) and disturbance of surface regolith particles. One goal of this event-based approach is to study how bulk behavior, such as the diffusion-like net downslope transport of soil, can emerge from a large ensemble of stochastic events. In this paper, we present the “Grain Hill” model and examine its ability to reproduce three common types of slope profile: (1) convex-upward, soil-mantled slopes (Fig. 2a, b), (2) quasi-planar rocky slopes (Fig. 2c, d), and (3) cliff-rampart morphology in layered strata (Fig. 2e, f).

We begin with a description of the modeling technique. We then present results that illustrate the macroscopic behavior of the model under a variety of boundary conditions, and define the relationship between the cellular model's parameters and the parameters of conventional continuum mechanics models for hillslope evolution.

The model combines a cellular automaton representation of granular mechanics with rules for weathering of rock to regolith and for episodic disturbance of regolith. Cellular automata are widely used in the granular mechanics community, because they can represent the essential physics of granular materials at a reasonably low computational cost. Because the principles are often similar to those of lattice-gas automata in fluid dynamics (Chen and Doolen, 1998), cellular automata for granular mechanics are sometimes referred to as lattice grain models (LGrMs) (Alonso and Herrmann, 1996; Cottenceau and Désérable, 2010; Désérable, 2002; Désérable et al., 2011; Gutt and Haff, 1990; Károlyi and Kertész, 1998, 1999; Károlyi et al., 1998; Martinez and Masson, 1998; Peng and Herrmann, 1994).

## 2.1 CTS lattice grain model

Our approach starts with a two-dimensional (2-D) CTS lattice grain cellular automaton. A cellular automaton can be
broadly defined as a computational model that consists of a lattice of cells,
with each cell taking on one of *N* discrete states (represented by an
integer value). These states evolve over time according to a set of rules
that describe transitions from one state to another as a function of a
particular cell's immediate neighborhood. A continuous-time stochastic model
is one in which the timing of transitions is probabilistic rather than
deterministic. Whereas transitions in traditional cellular automata occur in
discrete time steps, in a CTS model, they are both stochastic and
asynchronous. A CTS model can be viewed as a type of Boolean delay equation
(Ghil et al., 2008), though the number of possible states is not
necessarily limited to just two.

The method we present here, which we will refer to as the Grain Hill model, is implemented in the Landlab modeling framework (Hobley et al., 2017). The lattice grain component, on which Grain Hill builds, is described in detail by Tucker et al. (2016). Here, we present only a brief overview of the lattice grain model's rules and behavior. The framework is based on the pairwise (“doublet”) method developed by Narteau and colleagues (Rozier and Narteau, 2014), which has been applied to problems as diverse as eolian dune dynamics (Narteau et al., 2009; Zhang et al., 2010, 2012) and the core-mantle interface (Narteau et al., 2001).

In the basic CTS lattice grain model, the domain consists of a lattice of hexagonal cells. Each cell is assigned one of eight states (Table 1, states 0–7). These states represent the nature and motion status of the material: state 0 represents fluid (an “empty” cell into which a solid particle can move), states 1–6 represent a grain moving in one of the six lattice directions, and state 7 indicates a stationary grain (or aggregate of grains, as discussed below). For purposes of modeling hillslope evolution, we add an additional state (8) to represent rock, which is immobile until converted to granular material, representing regolith. An optional additional state (9) is used to model large blocks, as described below. Figure 3 shows several of these states in the form of a time sequence of transition events. Note that the timing of transition events is purely stochastic; there are no time steps in the usual sense.

Like other lattice grain models, the CTS lattice grain model is designed to represent, in a simple way, the motion and interaction of an ensemble of grains in a gravitational field. The physics of the material are represented by a set of transition rules, in which a given adjacent pair of states is assigned a certain probability per unit of time of undergoing a transition to a different pair. For example, consider a vertically aligned pair of cells in which the top cell has state 4 (moving downward) and the bottom cell state 0 (empty/fluid) (Fig. 3f). Downward motion (falling) is represented by a transition in which the two states switch places (Fig. 3g).

The stochastic pairwise transitions in the CTS lattice grain model are
treated as Poisson processes. The probability density function for the
waiting time, *t*, to the next transition event at a particular pair is given
by an exponential function with a rate parameter *r*, which has dimensions of
inverse time:

Each transition type is associated with a rate parameter that represents the speed of whichever process the transition is designed to represent. To implement these transitions, the CTS lattice grain model steps from one transition to the next, rather than iterating through time steps of fixed duration. Whenever the state of one or both cells in a particular pair changes, if the new pair is subject to a transition, the time at which the transition is scheduled to occur is added to a queue of pending events. The soonest among all pending events is chosen for processing, and the process repeats until either the desired run time has completed or there are no further events in the queue. Further details on the implementation and algorithms are provided in Tucker et al. (2016).

Grain motion through fluid is represented by a transition involving a moving grain and an adjacent fluid cell in the direction of the grain's motion: the two cells exchange states, representing the motion of the grain into the fluid-filled cell, and the replacement of the grain's former location with fluid (Fig. 3c, d and f, g). During this transition, the grain's motion direction remains unchanged (Fig. 4, top left). Note that the lattice itself never moves; rather, material motion is represented simply by an exchange of grain and fluid states between an adjacent pair of cells.

Gravity is represented by transitions in which a rising grain decelerates to
become stationary, a stationary grain accelerates downward to become a
falling particle, and a grain moving upward at an angle accelerates downward
to move downward at an angle (Fig. 5). An additional rule
allows for acceleration of a particle resting on a slope: a stationary
particle adjacent to a fluid cell below it and to one side may transition to
a moving particle (Fig. 5, bottom row). Importantly for our
purposes, this latter rule effectively imposes an angle of repose at
30^{∘}.

For gravitational transitions, the rate parameter, *r*_{g}, is determined by
considering the time it would take for an initially stationary object to fall
a distance of one cell width under gravitational acceleration without fluid
drag. This works out to be

where *δ* is cell width and g is gravitational acceleration. This rate
parameter is used for all of the gravitational transitions illustrated in
Fig. 5.

Because of the stochastic treatment of all transitions – including gravitational ones – it is possible for grains in the model to hover in mid-air for a brief period of time before plunging downward (Coyote, 1949). For purposes of modeling hillslope evolution, this is fine; what matters most is that there is a distinct timescale gap between “fast” (large rate constant) processes associated with grain motion and “slow” (small rate constant) processes associated with weathering and disturbance, which are described below. First, however, we must consider frictional interactions among moving particles.

We assume that biophysical disturbance events such as the growth of roots and
burrowing by animals, and the settling motions that follow, tend to impart
low kinetic energy, with “low” defined as ballistic displacement lengths
that are short relative to hillslope length and comparable to or less than
the characteristic disturbance-zone thickness. We consider such motions to be
dominated by frictional dissipation rather than by transfer of kinetic energy
by elastic impacts. This view is similar to the reasoning of
Furbish et al. (2009) that the mean free path of mobile grains will
typically be short relative to hillslope length, scaling with the grain
radius and particle concentration. For this reason, unlike the original
lattice grain model of Tucker et al. (2016), the present formulation
includes only inelastic collisions (Fig. 4). These inelastic
(frictional) collisions are represented by a set of rules in which one or
both colliding particles become stationary, representing loss of momentum and
kinetic energy as a result of the collision. The particular choices for
frictional interaction are motivated simply by the geometry of the problem.
They are non-unique in the sense that one could imagine reasonable
alternatives to the rules illustrated in Fig. 4; however,
the details of frictional interactions have little influence on the outcomes
of the Grain Hill model. In the general lattice grain CTS model, the rate
parameter for frictional transitions is set equal to the product of the
gravitational parameter and a dimensionless friction factor, $f\in [\mathrm{0},\mathrm{1}]$
(there is also a corresponding elastic factor equal to 1−*f*). In the Grain
Hill implementation explored in this paper, *f*=1, such that particle
collisions are purely frictional.

One limitation of the CTS lattice grain model is that falling grains do not accelerate through time; instead, they have a fixed transition probability that implies a statistically uniform downward fall velocity. This treatment is obviously unrealistic for particles falling in a vacuum, though it is consistent with a terminal settling velocity for grains immersed in fluid. Consistent with the above reasoning, the relatively short ballistic displacement lengths asserted for the modeled hillslopes also reduce the importance of this assumption, as a particle would typically have little time to accelerate before impacting another particle.

Tests of the CTS lattice grain model show that it reproduces several basic aspects of granular behavior (Tucker et al., 2016). For example, when gravity and friction are deactivated, the model conserves kinetic energy. When gravity and friction are active, the model reproduces some of the common behaviors observed with granular materials. For example, Fig. 6 illustrates a simulation of the emptying of a silo to form an angle-of-repose grain pile. For our purposes, what matters most is simply that the model captures, in a reasonable way, the response of particles on a slope to episodic disturbance events.

## 2.2 Weathering and soil creep

Weathering of rock to form mobile regolith is modeled with a transition rule:
when a rock cell lies adjacent to a fluid cell (which here is assumed to be
air), there is a specified probability per unit of time, *w* (1∕*T*), of
transition to a grain–air pair (Figs. 3a–b, and 7, top). In other words, *w* is the Poisson rate constant
for the weathering transition process. This treatment means that the
effective maximum expected weathering rate, in terms of the propagation of a
weathering front, is cell diameter, *δ*, multiplied by *w*. An indirect
consequence of this approach is that the weathering rate declines with
increasing regolith thickness. As average regolith thickness increases, the
fraction of the surface where rock is in contact with air diminishes, and
consequently so does the average transition rate. A limitation of the
approach is that when the rock is completely mantled, no further weathering
can take place. We explore the consequences of this rule below, and compare
it with the behavior of continuum regolith-production models.

Soil creep is modeled by a transition rule that mimics the process of
episodic disturbance of the mobile regolith (which we use here as a generic
term that includes various forms of unconsolidated granular material, such as
soil, colluvium, and scree). For each resting grain that is adjacent to an
air cell, there is a specified probability per unit of time, *d* (1∕*T*), that the
regolith and air will exchange places, representing movement
(Fig. 3b, c). The regolith cell is also converted from a
stationary state to a state of motion (Fig. 7b). An
advantage of this approach is that it mimics, in a general way, the
effectively stochastic disturbance processes that are understood to drive
soil creep.

Our definition of *d* is closely related to the activation rate, *N*_{a}, in
the probabilistic theory for soil creep developed by
Furbish et al. (2009). When combined with the lattice grain
gravitational rules, the resulting cellular model captures both the
scattering (disturbance) and settling (gravitational) behavior articulated by
Furbish et al. (2009). In the Grain Hill cellular model, as in their
theory, downslope regolith flux arises because, on average, scattering occurs
perpendicular to the local surface while setting is vertical. The Grain Hill model
includes an additional element not present in the
Furbish et al. (2009) theory: an increase in (downward) scattering
distance for particles on slopes steeper than 30^{∘}. This behavior, as
illustrated below, promotes a non-linear relationship between gradient and
flux, and leads to the possibility of threshold slopes.

Note that the weathering and disturbance rate constants (*w* and *d*,
respectively) are understood to be considerably smaller than the
gravitational rate constant, *r*_{g}. As noted above, a key concept here is
that there are two distinct timescales: a short timescale associated with
grain motion, and a much longer scale associated with weathering and
disturbance frequency.

## 2.3 Cells as grain aggregates

Natural regolith disturbance events usually impact many grains at once.
Raindrop impacts on bare sediment typically dislodge several grains at once
(Furbish et al., 2007). Excavation of an animal burrow disturbs a volume of
grains equal to the volume of the burrow, and the fall of a tree mobilizes a
volume of regolith similar to the volume of the tree's root mound.
Observations of such processes suggest that there may be a characteristic
volume of disturbance that in some cases may be much larger than the volume
of a single grain. For this reason, we envision regolith cells as being grain
aggregates, with a length scale (width of a cell) *δ* and a volume scale
*δ*^{3}.

## 2.4 Initial and boundary conditions

The 2-D model domain represents the cross section of a hypothetical hillslope,
on which particles move within the cross-sectional plane. Any regolith cells
that reach the model's side or top boundaries disappear. This treatment is
meant to represent the presence of a stream channel at the base of each side
of the model hillslope; particles reaching these channels are assumed to be
eroded. Progressive lowering of the base level at the two model boundaries is
treated by moving the interior cells upward away from the lower boundary, and
adding a new row of rock or regolith cells along the bottom row. A new row of
cells is added at time intervals of *τ*.

Cells around the lattice perimeter retain their initial states. If, for example, a transition occurs in which a grain “moves” into a fluid cell on the lattice perimeter, its former location will correctly transition to a fluid cell, but the perimeter cell itself will retain its status as a fluid cell. Effectively, this treatment means that grains or blocks reaching either of the two vertical boundaries are instantly eroded.

The initial condition for most runs presented here has the bottom two rows filled with regolith grains. The lower left and lower right cells are assigned to be rock, which represents the base-level (and incidentally helps keep a consistent color scheme among different model configurations, because the rock state is always present). The rest of the domain is initialized as air cells.

## 2.5 Scaling and non-dimensionalization

The basic model has four parameters: the disturbance rate, *d* (cells ∕ *T*),
weathering rate, *w* (cells ∕ *T*), base-level lowering interval, *τ* (*T*), and
width of domain, *λ* (cells). The base-level lowering timescale *τ*
represents the time interval between episodes of relative uplift in which the
interior domain is lifted by one cell relative to its side boundaries. The
domain width might properly be considered a boundary condition rather than a
parameter, but we include it here with an eye toward examining how slope
width impacts hillslope properties such as mean height. Once we define the
width of a cell, *δ* (*L*), we can define versions of these four
parameters that explicitly incorporate this length scale:

Consider the case of dynamic equilibrium, in which the rate of base-level
lowering is balanced by the hillslope's rate of erosion. The mean height of
this steady state hillslope, *H*, is a function of the above four parameters
plus the characteristic length scale *δ*, such that we end up with a
total of six variables:

Buckingham's Pi theorem dictates that these six variables, which collectively include dimensions of length and time, may be grouped into four dimensionless quantities:

The ratio ${d}^{\prime}=d\mathit{\tau}=D/U$ is a dimensionless disturbance rate. Similarly, ${w}^{\prime}=w\mathit{\tau}=W/U$ is a dimensionless weathering rate. Noting the definitions above, Eq. (8) is equivalent to

where $h=H/\mathit{\delta}$ is dimensionless hillslope height. Hence, we have a
dimensionless property of the hillslope, *h*, that depends uniquely on three
other non-dimensional variables, representing disturbance rate, weathering
rate, and length.

One can similarly define a dimensionless regolith thickness, $r=R/\mathit{\delta}$,
where *R* is the dimensional equivalent; it too should depend on the three
dimensionless parameters that represent disturbance rate *d*^{′}, weathering
rate, *w*^{′}, and hillslope length, *λ*, respectively. For a hillslope
composed entirely of regolith, *r* and *h* depend solely on *d*^{′} and
*λ*. Finally, we define a fractional regolith cover *F*_{r}. In the Grain
Hill model, *F*_{r} is calculated as the number of air–regolith cell pairs
divided by the total number of cell pairs that juxtapose air with either
regolith or rock.

## 2.6 Blocks

The foregoing model is designed to represent regolith as grain aggregates composed of gravel-sized and finer grains: material fine enough that it is susceptible to being moved by processes such as animal burrowing, frost heave, tree throw, and so on. Some hillslopes, however, are adorned with grains that are simply too large to be displaced significantly by such processes. For example, Glade et al. (2017) presented a case study and model of slopes formed beneath a resistant rock unit that periodically sheds meter-scale or larger blocks. On at least some of these types of slope, the distance between surface blocks and their source unit is considerably greater than the distance that they could roll during an initial release event (Duszyński and Migoń, 2015). This observation implies that the blocks are transported downslope by a process of repeated undermining. Glade et al. (2017) hypothesized that erosion of soil beneath and immediately downhill can cause a block to topple and hence move a distance comparable to its own diameter in each such event.

We wish to capture this form of “too big to disturb” behavior in the Grain
Hill model. The CTS approach, at least as it is defined here, does not lend
itself to variations in grain size or geometry. Instead, we introduce an
additional type of particle that represents the behavior of blocks rather
than treating their difference in size explicitly. In a sense, the approach
can be viewed as treating blocks as having greater density, rather than
greater size, than other grains. A block particle differs from a normal
regolith cell in that it cannot be scattered upward by disturbance. Motion of
a block particle can only occur under two circumstances: when it lies
directly above an air cell (in which case it falls vertically, trading places
with the air cell) and when it lies above and to the side of an air cell (in
which case it falls downslope at a 30^{∘} angle, with probability per
time *d*). These rules mimic the undermining process discussed by
Glade et al. (2017).

As in the Glade et al. (2017) model, block particles can also undergo
weathering. Here, weathering is again treated in a probabilistic fashion:
blocks form from weathering of bedrock, at probability per time *w*. Once
created, a block can undergo a conversion to normal regolith with probability
*w* when it sits adjacent to an air cell. This treatment of blocks captures,
in a simple way, the weathering of blocks as they move downslope. For
purposes of this paper, the block component is included simply to test
whether a cellular automaton treatment produces results that are
qualitatively consistent with observations, and also consistent with the
hybrid continuum–discrete model of Glade et al. (2017) and
Glade and Anderson (2017).

## 3.1 Fully soil-mantled hillslope

We start by considering the case of fully soil-mantled hillslopes, in which
the supply of mobile regolith is effectively unlimited
(Fig. 2a, b). Under this condition, the Grain Hill model
represents a testable mechanistic hypothesis: that a transport limited,
soil-mantled hillslope behaves essentially as a granular medium subject to
periodic, quasi-random disturbance events. This concept was also the essence
of the acoustic-disturbance experiments by Roering et al. (2001). To
test the hypothesis, we run the Grain Hill model with a constant rate of
material uplift relative to base-level until the system reaches quasi-steady
state, to determine whether its steady form is smoothly convex upward (when
the gradient is below the failure threshold) to planar (when the gradient
lies at or near the failure threshold). Model runs were performed using a
251-row by 580-column lattice. Disturbance rates were varied from 0.001 to 0.1 yr^{−1} and intervals between relative-uplift events from
100 to 10 000 years.

Results show that the Grain Hill model produces parabolic to planar hillslope
forms, depending on the ratio of disturbance to uplift rates, which is
encapsulated in the dimensionless parameter *d*^{′} (Fig. 8). At high
*d*^{′} (frequent disturbance and/or slow base-level fall), hillslope relief is
low and the form is smoothly convex upward (Fig. 8, lower right
panels). At somewhat lower *d*^{′}, the lower part of the slope approaches a
threshold angle while the upper part remains smoothly convex
(Fig. 8, middle diagonal panels). At low *d*^{′}, the form becomes
predominantly planar and achieves a threshold relief that is insensitive to
further increases in *d*^{′} (Fig. 8, upper left panels).

Scaling of mean height as a function of *d*^{′} is shown in
Fig. 9. The figure shows results for 125 model runs
spanning 2 orders of magnitude in each parameter (*d*, *τ*, and
*λ*) in half-decade intervals. The 125 runs represent a $\mathrm{5}\times \mathrm{5}\times \mathrm{5}$ grid of experiments, in which each grid point represents a
particular combination of the three parameters (*d*, *τ*, and *λ*).

For any given hillslope length, there are three regimes of behavior. Low *d*^{′}
(upper left of graph) leads to threshold hillslopes, in which relief depends
only on hillslope length. Under moderate *d*^{′}, mean height scales inversely
with *d*^{′}, as expected from linear diffusion theory. At high *d*^{′}, we have a
finite-size regime in which dimensionless hillslope mean height is comparable
to the disturbance scale, *δ* (cell size in the model); in other words,
the hill is only one or a few cells high.

The behavior of the Grain Hill model in its simple, transport-limited
configuration can be compared to diffusion theory, which relates volumetric
sediment flux per unit contour length, **q**_{s}, to topographic gradient:

where *η* is land-surface height, *x* is horizontal distance, and *D*_{s} is
an effective transport coefficient. The Furbish et al. (2009)
probabilistic theory for transport due to particle scattering and settling
formulates *D*_{s} as

where *k* is a dimensionless coefficient, *r*_{p} is particle radius, *R*_{a} is
active regolith thickness, *N*_{a} is the activation rate, *θ* is slope
angle, *c* is particle concentration, and *c*_{m} is a maximum concentration.
The over-bar denotes an average over the active regolith thickness. For the
Grain Hill model, *R*_{a} scales with the characteristic disturbance depth,
*δ*. Further, because we treat grain aggregates, we may also assume
*r*_{p}∼*δ*. Therefore, we have the prediction that

where *a* is a dimensionless proportionality constant.

The mean expected activation rate, *N*_{a}, is closely related to the Grain
Hill model's disturbance frequency parameter, *d*. To relate the two
quantitatively, one needs to make a trivial lattice-geometry correction. A
straight-as-possible cut through the hex lattice exposes on average two faces
per cell, both of which are susceptible to a disturbance event. Because *d*
is the expected disturbance frequency per cell face, and because independent
Poisson events are additive, the resultant disturbance frequency for each
cell exposed along a quasi-horizontal surface is *N*_{a}=2*d*.

A more important difference is that whereas *N*_{a} is defined as activation
rate per unit horizontal area, *d* represents the rate per unit surface area
regardless of orientation. For a given *d*, *N*_{a} will increase with surface
roughness (because there is more exposed area of regolith–air contact), and
with gradient (because the slope length increases).

An additional effect arises from the model's effective 30^{∘} angle of
repose. On slopes steeper than this, the expected disturbance rate increases
substantially because gravitational dislodgement becomes activated
(Fig. 5, bottom row). Thus, the Grain Hill model incorporates
an additional non-linear relationship between flux and gradient inasmuch as
*N*_{a} depends on gradient.

We can derive an effective diffusivity, *D*_{e}, from the modeled topography by
applying the expected relationship between mean elevation and diffusivity.
Here, *D*_{e} is defined as that value which, if it were spatially uniform,
would yield the same mean steady-state elevation as that produced by the
particle model. Framing it this way allows us to interrogate how the
effective transport coefficient varies as a function of mean slope gradient.
At steady state, mass balance implies that

where *E* is the rate of erosion – equal to the rate of material uplift
relative to base level – and *x* is horizontal distance from the ridge top.
Substituting Eq. (10) and rearranging gives

Integrating and then averaging over *x*, we can solve for the average
elevation, $\stackrel{\mathrm{\u203e}}{\mathit{\eta}}$:

where ${L}_{h}=L/\mathrm{2}$ is the length of the slope from ridge top to base (in other
words, half the total length of the domain). We can then rearrange this to
find *D*_{e}:

To examine how *D*_{e} scales, we can define a dimensionless form, normalizing
by the disturbance frequency, *d*, and the square of active regolith
thickness (equal to particle diameter), *δ*^{2}:

Noting that $E=\mathit{\delta}/\mathit{\tau}$, *L*=2*L*_{h}, and $L/\mathit{\delta}=\mathit{\lambda}$, this is equivalent to

where $\stackrel{\mathrm{\u203e}}{h}$ is the mean hillslope height in particle diameters.

As expected, ${D}_{\mathrm{e}}^{\prime}$ increases with hillslope gradient
(Fig. 10). The effective diffusivity approaches an asymptote at
30^{∘} (mean gradient ≈ 0.6), representing an angle of repose. The
pattern resembles the family of non-linear flux–gradient curves introduced by
Andrews and Bucknam (1987) and explored further by
Howard (1994) and Roering et al. (1999). At low
gradients, ${D}_{\mathrm{e}}^{\prime}$ approaches a value of about 60. (This method of estimating
${D}_{\mathrm{e}}^{\prime}$ is similar to fitting the standard theoretical parabolic curve to the
experimental profiles, except that here we use the integral of the profiles.)

The link between *D*_{e} and *d* provides a way to scale the Grain Hill model
to field-derived estimates of *D*_{s} and *R*_{a}. Here, we equate the theoretical
effective diffusivity, *D*_{e}, with the definition of the transport
coefficient *D*_{s} of Furbish et al. (2009). Noting that, at low
gradients, cos^{2}*θ* in Eq. (12) approaches unity, and using
the prior relation *N*_{a}=2*d*, we may write *D*_{s} for low slope angle as

In the Grain Hill model, the fact that low-angle ${D}_{\mathrm{e}}^{\prime}\approx \mathrm{60}$ implies
that the dimensional equivalent ${D}_{\mathrm{e}}(\mathit{\theta}\to \mathrm{0})\approx \mathrm{60}\phantom{\rule{0.125em}{0ex}}{\mathit{\delta}}^{\mathrm{2}}d$. Equating *D*_{s} (the transport coefficient derived by
Furbish et al., 2009) and *D*_{e} (the effective transport
coefficient derived from the Grain Hill model),

This relation can be used to scale the parameters in the Grain Hill model
with field data. For example, if one were to assume an active regolith
thickness of 0.4 m and a low-gradient transport coefficient of *D*_{s}=0.01 m^{2} yr^{−1}, and set *δ* to the active regolith thickness, then

Here, *d* represents the frequency with which a given exposed patch of
regolith of width and depth *δ* is disturbed upward. With the above
values, the simulated hills in Fig. 8 would be 232 m long
(valley to valley), with height ranging from 1.6 to 57.6 m and base-level
lowering rate from 0.04 to 4 mm yr^{−1}.

## 3.2 Hillslope with regolith production from rock

Having established that the Grain Hill model reproduces classic soil-mantled
hillslope forms and has parameters that can be related to the parameters in
commonly used continuum hillslope transport theories, we turn now to the case
in which regolith is generated from bedrock with a production rate that may
(or may not) limit the rate of erosion. We explore the role of regolith
production with a series of model runs in which *w*^{′} varies from 0.4 to 40.
The upper end of this range represents a condition in which the potential
maximum rate of regolith production greatly exceeds the rate of base-level
lowering. The lower end, 0.4, is less than the rate of base-level fall, and
would seem to be insufficient to allow for equilibrium to occur, and yet
nonetheless it does. Examples of equilibrium hillslope forms found in this
parameter space are shown in Fig. 11.

Relationships among mean gradient, fractional regolith cover, dimensionless
disturbance rate *d*^{′}, and dimensionless weathering rate *w*^{′} are illustrated
in Fig. 12. For ${w}^{\prime}>\mathrm{1}$, the gradient–*d*^{′} relation
(Fig. 12a) has the same shape as in the purely regolith models:
a threshold regime at lower *d*^{′} transitioning to an inverse gradient–*d*^{′}
relation at higher *d*^{′}. This indicates that when the maximum weathering rate
(for a flat surface) is substantially greater than the rate of base-level
fall, we recapture transport-limited conditions. With ${w}^{\prime}<\mathrm{1}$, however, the
hillslope achieves an equilibrium gradient that is greater than that for the
transport-limited case, and at lower *d*^{′}, is greater than the threshold
angle (Fig. 12a, b).

We can also examine the fractional regolith cover, which is defined here as
the number of rock–air cell pairs divided by the total number of cell pairs
at which air meets either regolith or rock (Fig. 12c, d). The
fractional regolith cover shows relatively little sensitivity to *d*^{′}
(Fig. 12c). The cover hovers around unity for high *w*^{′} and
*d*^{′} but systematically declines with *w*^{′} when *w*^{′} is below about 10.
(Note that the data points representing ${d}^{\prime}=\mathrm{1000}$ and ${w}^{\prime}>\mathrm{1}$ have
hillslope heights of only a few particles and are therefore sensitive to
finite-size effects.)

The models with ${w}^{\prime}<\mathrm{1}$ present a seeming paradox: how is it possible to
achieve an equilibrium form when the maximum weathering rate appears to be
lower than the rate of uplift relative to base level? The solution to the
paradox lies in surface area. The surface area of rock that is exposed to
weathering is not fixed but rather depends on the overall slope length, the
terrain roughness, and the fractional regolith cover. To appreciate the first
effect, consider a planar slope at angle *θ* with no regolith cover. If
*w**δ* represents the maximum slope-normal bedrock weathering rate, then
the vertical rate is simply *w**δ*∕cos*θ*. All else equal,
increasing gradient will increase vertical weathering rate, thereby providing
a feedback between gradient and rock lowering rate. A second feedback relates
to topographic roughness: all else equal, a rougher surface will experience a
greater weathering rate because it provides more surface area. The third
feedback, which is embedded in the depth-dependent regolith production
hypothesis (Ahnert, 1967) lies in regolith cover: the greater the
exposure of rock (or the thinner the cover), the faster the average rate of
rock-to-regolith conversion. In the Grain Hill model, this third feedback is
represented by fractional bedrock exposure (since weathering only occurs when
rock cells are juxtaposed with air cells).

To test whether these are indeed the feedbacks responsible for equilibrium
topography in the Grain Hill model, we can compare the rate of material
influx (uplift relative to base-level) with the expected rate of
rock-to-regolith conversion. In the Grain Hill model, the expected rate of
regolith production, *P*, in cross-sectional area per time, is the product of
weathering rate per cell face, *w*, the cross-sectional area of a cell, *A*,
and the number of rock–air cell faces, *n*_{ra},

The rate of material addition due to uplift relative to base-level, *U*, again
in cross-sectional area per time, is the area of a cell, *A*, multiplied by the
horizontal width of the domain in cells, *n*_{H}, divided by the interval
between uplift events, *τ*:

Equality between rock uplift and weathering can be expressed as

The ratio on the right side represents the surface-area effect, in the form
of surface area exposed to weathering per unit horizontal area. The balance
is illustrated in Fig. 13, which compares the left-hand and
right-hand terms for each of the 125 model runs with weathering. Each data
point represents a single snapshot in time, and so scatter is to be expected.
To help diagnose the scatter around the 1:1 line, the data are divided into
quintiles by fractional regolith cover, *F*_{r} (note that some of the points
in the lower quintiles are obscured by being over-plotted along the 1:1
line). Many of the points that fall off the 1:1 line, especially at the high
end (higher 1∕*τ*), come from runs with *F*_{r}>80 *%*; with very few exposed
rock–air pairs, a small fluctuation in the *n*_{ra} can produce a relatively
large change in predicted weathering rate. At the low end, many of the points
above the 1:1 line come from runs with a maximum height of only a few cells,
which are subject to finite-size effects.

The main message of Fig. 13 is that the Grain Hill model demonstrates an equilibrium adjustment between rock uplift and rock weathering. The weathering rate does not have a fixed upper “speed limit,” but rather is set by the exposed surface area, which in turn is a function of gradient, roughness, and regolith cover. Solutions with a discontinuous regolith cover are indicative of this adjustment. Slopes can grow arbitrarily steep, with weathering and erosion increasingly attacking from the sides as the gradient rises.

## 3.3 Comparison between weathering rule and inverse-exponential model

The most popular function to describe regolith production from bedrock is the decaying exponential formula proposed by Ahnert (1967), which has proved consistent with estimates of production rate obtained using cosmogenic radionuclides (Heimsath et al., 1997; Small et al., 1999). The production rate is given by

where *P*_{0} is the maximum (bare bedrock) production rate, *R* is regolith
thickness, and *R*_{*} is a depth-decay scale on the order of decimeters. On a
flat surface, assuming no erosion or deposition, the expected rate of change
of *R* over time is

where *ρ*_{r} and *ρ*_{s} are the bulk densities of
parent material and regolith, respectively, and *ω* is the fraction of
parent material removed in solution upon weathering. Starting from a bare
surface, and assuming isovolumetric weathering (in which case
${\mathit{\rho}}_{\mathrm{s}}=(\mathrm{1}-\mathit{\omega}){\mathit{\rho}}_{\mathrm{r}}$), the expected regolith
thickness as a function of time can be found by integrating Eq. (26):

We can compare this with the behavior of the cellular weathering rule by
running the case of a flat, initially bare-rock surface from which weathered
material may neither enter nor leave (Fig. 14, case $d\phantom{\rule{0.125em}{0ex}}{w}^{-\mathrm{1}}=\mathrm{0}$).
When the disturbance rate is zero, the cellular weathering model
asymptotically approaches a steady regolith thickness of exactly one cell
(thickness equal to *δ*). This is so because the model allows weathering to
occur only when rock cells are exposed to air cells, and there is no
disturbance process that would juxtapose rock and air once the initial
weathered layer has formed. When disturbance rate is non-zero, however,
regolith continues to form even after the mean thickness *r* exceeds unity
(representing one characteristic disturbance depth). Continuation of regolith
production occurs because the disturbance process intermittently exposes
rock, at which point it becomes subject to weathering. The greater the
disturbance rate, the more frequent the exposure and hence the more rapid
weathering (Fig. 14). For any ratio *d* *w*^{−1}, the model's
weathering behavior clearly differs from the logarithmic growth in thickness
predicted by exponential theory. This represents both a strength and a
weakness in the Grain Hill model. On the one hand, the model under its
present configuration cannot account for rock-to-regolith conversion
resulting from processes that penetrate more than one characteristic
disturbance depth *δ* into the subsurface. For example, the model
neglects the possibility that some plant roots may penetrate deeply and
contribute to disaggregation, or that an unusually deep freezing front in a
cold winter might cause rock fracture and displacement of the resulting
fragments (Anderson et al., 2012). On the other hand, the model
honors the likelihood that soil disturbance and regolith production are
closely linked processes, rather than independent: all else equal, a
greater disturbance rate will tend to produce faster rates of both regolith
production and downslope soil movement.

## 3.4 Rock collapse and vertical cliffs

Some rock slopes display a cliff-and-rampart morphology in which a vertical or near-vertical rock face stands above an inclined, often sediment-mantled buttress (Figs. 1 and 15). Although common in sedimentary rocks where a resistant unit forms the cliff and a weaker unit the buttress, the same morphology is sometimes found in apparently homogeneous lithology (Fig. 15b). The cliff portion of such slopes suggests a process of undermining and collapse, with the cliff-forming material being cohesive enough to maintain a vertical face but too weak to support overhangs.

To explore the origins of ramp-and-cliff morphology, we consider a version of the Grain Hill model that adds an extra rule to represent collapse: any rock particle that directly overlies air has the possibility to transition to a falling regolith particle, with the same rate as gravitational transition from resting to falling; in other words, as soon as a rock particle has been undermined, it behaves like cohesionless material.

Under dynamic equilibrium, this rule produces a morphology with slopes that
are roughly planar, with alternating vertical and sloping sections and patchy
regolith cover (Fig. 16). With ${w}^{\prime}\le \mathrm{1}$, gradient and
regolith cover depend strongly on *w*^{′} and show little or no sensitivity to
*d*^{′}. When ${w}^{\prime}\ll \mathrm{1}$, the hillslope forms resemble pinnacles. These examples
demonstrate two combined feedbacks between weathering and base-level fall: the
surface area susceptible to weathering, and the frequency and magnitude of
material collapse through undermining.

The case of transient evolution under a stable base level leads to the formation of a regolith-mantled, angle-of-repose ramp (Fig. 17). The slope break remains relatively sharp as it retreats headward. The ramp forms as a transport slope. The angle of repose is an attractor state: if the angle were steeper, weathered material would be rapidly removed as a result of gravitational instability; if it were substantially lower, material would accumulate, because transport would be limited to the (much lower) rate of disturbance-driven creep motions. Hence, the Grain Hill model predicts that formation of a sediment-mantled ramp beneath a steeper, actively weathering rock slope is an expected outcome for a steep rock slope under stable base level.

## 3.5 Blocks

Weathering and erosion in landscapes underlain by relatively massive, fracture- or joint-bounded rock can sometimes produce large “blocks” of rock, defined here as clasts that are too large to be displaced upward by normal hillslope processes. The release of blocks from dipping sedimentary or volcanic strata can alter both the shape and relief of hillslopes (Glade et al., 2017). When blocks are delivered to streams, they can influence the channel's roughness, gradient, erosion rate, and longitudinal profile shape (Shobe et al., 2016).

As discussed in Sect. 2.6, the Grain Hill model can be modified to honor blocks by defining an additional cell type that represents blocks. The weathering process is modified such that a rock cell now weathers into a block, and the block in turn may weather to form regolith. When a block is undermined directly from below, it will fall just as a normal regolith particle would. When a block particle lies adjacent to and above an air cell, a disturbance event may occur that causes the block to shift downward on the slope. By these means, blocks in the model may move downward, or downward and laterally, but never upward. An implicit assumption in this treatment is that blocks do not roll long distances (further than their own diameter) upon release.

We examine model runs in which a resistant rock layer is embedded in a weak sedimentary material that is soft enough to be treated as regolith (Fig. 18). The modeled hillslopes are qualitatively consistent both with field observations and with the mixed continuum–discrete model of Glade et al. (2017) and Glade and Anderson (2017) in that block-mantled slopes are generally concave upward, reflecting a downslope decrease in the flux of blocks as weathering progressively transmutes them into regolith.

We perform a basic validation of the Grain Hill model by comparing its output
to real field sites, testing whether the model is capable of reproducing
realistic hillslope forms at the correct spatial scale under known boundary
conditions. Field sites were chosen such that model boundary conditions could
be derived from independent field estimates of rate parameters such as *D*_{s}
and the rate of base-level fall. To perform this test, we consider two
examples: a convex-upward, soil-mantled hillslope in Gabilan Mesa,
California, USA (Fig. 2a, b), and a steep, quasi-planar,
discontinuously mantled hillslope in the Yucaipa Ridge, California, USA
(Fig. 2c, d). For each of these two case studies, the
hillslopes appear to be approximately at steady state, and independent
estimates exist for the rate of base-level fall, *U*
(Binnie et al., 2007; Perron et al., 2009, 2012). We estimated
the effective transport coefficient, *D*_{s}, for the profiles shown in
Fig. 2a, c by measuring the second derivative of the
one-dimensional hillslope elevation profiles, $\frac{{\partial}^{\mathrm{2}}\mathit{\eta}}{\partial {x}^{\mathrm{2}}}$, and solving for *D*_{s} using

For the Gabilan Mesa profile, we estimated the profile-averaged effective
transport coefficient as 0.0345 m^{2}yr^{−1}. The effective rate of
base-level lowering has been estimated at $U\approx \mathrm{1.47}\times {\mathrm{10}}^{-\mathrm{4}}$ m yr^{−1} (Perron et al., 2012). To construct a Grain Hill model for the
Gabilan profile, we begin by assuming a characteristic disturbance depth of
*δ*=1 m. This value was chosen to be consistent with measured soil
depths that typically range between 0.2 and 1.2 m (Johnstone et al., 2017).
We treat the system as transport-limited, consisting of mobile material, so
that weathering is not explicitly modeled. The disturbance parameter, *d*, is
then calculated from the independently estimated value of *D*_{s} using
Eq. (20). The interval between uplift events is $\mathit{\tau}=\mathit{\delta}/U\approx \mathrm{6800}$ years. The resulting modeled equilibrium profile provides a
reasonably good match to the observed Gabilan profile, with a convex-upward
shape and a hilltop height of about 45 m above the slope base
(Fig. 19a).

For Yucaipa Ridge, we estimated the transport coefficient at
*D*_{s}∼0.028 m^{2} yr^{−1} on the basis of hilltop curvature and an
estimated effective rate of base-level lowering of
≈ 0.0027 m yr^{−1} (Binnie et al., 2007). Using Eq. (20), this equates to a
disturbance-rate parameter *d*=0.00468 yr^{−1} and an uplift interval of
370 years in the Grain Hill model. Bedrock outcrops are common on the Yucaipa
Ridge hillslopes, implying a thin, discontinuous regolith cover. We therefore
treat the system as consisting of bedrock that must be weathered before it
can become mobile. Because we do not have independent information on the
effective maximum rock weathering rate, the Yucaipa case is a somewhat weaker
test: we can only ask whether there exists a geologically reasonable value of
*w* such that the model reproduces the observed relief and shape of the slope
profile. Through trial and error, we find that with a weathering rate
parameter *w*=0.002 yr^{−1} (which corresponds to a maximum regolith
production rate of 2 mm yr^{−1}), the model does a credible job of capturing
the shape and size of the Yucaipa profile (compare Fig. 2c
with 19b). Although this particular value was obtained
through a simple calibration process, it is at least both geologically
reasonable and, as one might expect, somewhat lower than the rate of
base-level lowering.

To test the sensitivity of the Yucaipa example to the assumed characteristic
disturbance depth, we ran a second experiment in which *δ* was reduced
to 0.75 m, and the weathering, disturbance rate, and uplift parameters were
rescaled accordingly. The relief and mean gradients of the two cases are
nearly identical, with planar slopes and a relief in both cases of
∼100 m.

These two examples demonstrate that the Grain Hill model parameters are not arbitrary but instead can be linked through straightforward reasoning to field estimates of transport efficiency and base-level lowering. When one does so, the model successfully reproduces both the shape and scale of observed slopes.

With just three parameters – disturbance frequency (*d*), characteristic
disturbance depth (*δ*), and base-level fall frequency (*u*) – the Grain
Hill algorithm can reproduce the convex-upward to quasi-planar forms
associated with soil-mantled hillslopes (Fig. 8). With the
addition of a parameter that represents rock-to-regolith conversion rate, the
algorithm accommodates partly mantled, rocky hillslopes
(Figs. 11, 16, 17). By adding a rule
for detachment of blocks from resistant rock, the model reproduces hillslope
forms associated with hogbacks and ledge-forming escarpments
(Fig. 18).

A common criticism of cellular automaton models is that they involve
arbitrary rules and/or parameters that can neither be measured nor verified
in the real world. That is not the case for the Grain Hill model, for which
the parameters are tied to measurable physical quantities. For example, the
disturbance frequency *d* is directly related to the frequency parameter
*N*_{a} in statistical theory of soil transport developed by
Furbish et al. (2009), and through that theory to the diffusion-like
transport coefficient *D*_{s} that is commonly estimated in field studies. This
connection between model parameters and field measurements is illustrated by
the model's ability to reproduce the correct shape and scale of observed
hillslope forms when estimates of *D*_{s} and *U* are available
(Figs. 2, 19). In the transport-limited
case, there are no tunable parameters: given independent estimates of *D*_{s}
and *U*, the correct morphology is recovered (Figs. 2a,
19a). In the case where rock weathering appears to play a
role, and an independent estimate of *P*_{0} is not available, the model
requires an estimation of maximum weathering rate *w*. Nonetheless, a
plausible value of *w* (0.002 m yr^{−1}), somewhat smaller than the rate of
base-level fall (0.0027 m yr^{−1}), reproduces the observed shape and relief
in the Yucaipa Ridge case study.

The transport dynamics predicted by the Grain Hill model are consistent with continuum soil-transport theory, which treats soil as a fluid with a downslope flow rate that depends on slope gradient. Like the popular Andrews–Bucknam non-linear transport law (Andrews and Bucknam, 1987; Howard, 1994; Roering et al., 1999), the transport-limited form of the Grain Hill model predicts diffusion-like behavior in which the effective diffusivity increases with slope gradient, with an asymptote at a threshold angle (Fig. 10). In one sense, the Grain Hill model is actually closer to the process level than fluid-like continuum models, because net downslope mass flux arises from a sequence of stochastic disturbance events rather than being dictated by a macroscopic transport law.

One limitation of the Grain Hill model is that its threshold-like behavior
arises from the lattice geometry: regolith cells perched at a 30^{∘}
angle above and to one side of an air cell are treated as unconditionally
unstable. Whereas the timing of motion is treated as a stochastic process,
the occurrence of motion is inevitable (unless some other event occurs
first). This treatment neglects the possibility of frictional locking among
noncohesive grains at angles somewhat above 30^{∘}, as well as the
possibility of cohesion. This limitation could be overcome by introducing a
probabilistic treatment of grain stability: a grain aggregate will be stable
with a given probability *p*, and unstable with probability 1−*p*. Such a
treatment would introduce an additional parameter, but this parameter could
in principle be estimated from physical experiments. The addition of a
“sticking rule” like this might also make it possible for models with
alternative lattice geometries to manifest the same dynamics, thereby
decoupling the basic model framework from the geometry of the lattice on
which it is implemented.

The inclusion of rock-to-regolith conversion enables the Grain Hill model to predict a continuum of slope forms from fully soil mantled to intermittently covered to bare. However, there are several limitations in the treatment of regolith production that could be improved on. The weathering rule assumes that regolith production can only occur when rock is exposed to air, which obviously neglects the role of shallow subsurface processes such as root or frost wedging. The effective weathering depth scale is the same as the disturbance scale, and equal to the cell size. This assumption is probably reasonable if the processes responsible for weathering and disturbance were one and the same, but not if they are distinct processes with different length scales. The Grain Hill model also does not account explicitly for chemical weathering, which in some cases can extend well below the surface. Finally, the model's effective regolith-production behavior does not follow the log-growth curve predicted by inverse-exponential theory for a stable surface (Fig. 14). With these caveats in mind, one advantage of the stochastic model of regolith production is that it effectively treats the disturbance and regolith-production processes as being closely linked: all else equal, the production rate is higher when disturbance is more frequent.

The popular inverse-exponential model for regolith production implies the existence of a speed limit to landscape evolution: in the absence of rock landsliding, erosion rate cannot exceed the maximum rate of rock-to-regolith conversion. Moreover, the model implies the existence of a bare landscape once the rate of erosion exceeds the maximum rate of regolith production. Heimsath et al. (2012) found evidence, however, that in fact there are additional stabilizing mechanisms, and that these manifest in landscapes with thin, patchy soils. The Grain Hill model is consistent with these observations in that it predicts the natural emergence of a discontinuous regolith cover, with the fractional cover exerting an influence on the average rate of weathering and erosion. Furthermore, the model behavior highlights the importance of slope length and roughness in modulating the regolith production rate: all else equal, steeper or rougher slopes allow higher production rates, leading to an additional feedback between relief and erosion rate for rocky hillslopes. The possibility of rock collapse upon undermining by weathering provides another feedback mechanism that may allow rates of erosion to exceed the flat-surface maximum regolith production rate (Fig. 16).

The Grain Hill model also provides insight into transient evolution of rocky slopes. Experiments on the relaxation of rocky slopes that are steeper than the threshold angle predict the formation of a regolith-mantled pediment at the angle of repose, which extends upslope as the steep upper slope gradually recedes (Fig. 17). This scarp–pediment morphology emerges without any variation in material strength, requiring only a period of base-level stability.

As a computational framework for exploring hillslope forms, the Grain Hill model has the advantage that it provides a mechanistic link between events (disturbance and weathering) and long-term morphologic evolution, without the need to specify a flux law. The model has the further advantage of being fully two dimensional, allowing disturbance and weathering events to initiate from the side as well as vertically. A further key element is that the model can mix timescales: a short timescale associated with grain motion, an intermediate timescale associated with disturbance events, and a much longer timescale for slope evolution. Mixing these disparate timescales in a single computer model is made possible by the fact that most of the time grains are stationary: the algorithm operates on small (stochastic) time steps during those moments when grains are moving, and on much longer steps when no grains are in motion (for further information on the discrete-event algorithm behind the model, see Tucker et al., 2016).

The Grain Hill framework has several important limitations. It is not
practical to simulate motion of individual grains unless the spatial scale is
quite limited (e.g., Fig. 6) or the grains are unusually large
(Fig. 18). If one wished to model individual grains (of order, say,
10^{−3} m) at the scale of a hillslope (of order 10^{2} m), a much more
efficient solution algorithm would be needed. Furthermore, the nature of a
cellular automaton is such that physical interactions are limited to adjacent
cells only; long-distance effects such as stress transmission cannot easily
be represented. In one sense, the restriction to short-range influence could
be seen as an advantage, in that it forces one to think about how it is that
mass or energy is actually transmitted in a granular medium. But the
restriction means that well-known principles such as solid-state stress
cannot easily be represented. On the other hand, the model does capture
non-local transport, in which particles set in motion can travel a distance
comparable to the slope length
(Foufoula-Georgiou et al., 2010; Furbish and Roering, 2013; Tucker and Bradley, 2010). Non-local
transport emerges in the Grain Hill model when the slope angle is near or
above 30^{∘}, such that there is a high probability that a disturbed
particle will land in an unstable location and continue moving without the
need for a second disturbance event.

A further limitation concerns the fixed cell size. Because the model is restricted to a fixed cell size, the Grain Hill framework does not lend itself to treatment of multiple grain sizes (apart from the simple “aggregates and blocks” approach illustrated in Fig. 18). Despite these limitations, the Grain Hill model provides a useful framework for exploring hillslope process and form in the context of stochastic events.

A continuous-time stochastic cellular automaton model known as the Grain Hill model allows for computational simulation of two-dimensional slope forms that arise from stochastic disturbance and (possibly) weathering events. The model operates on a hexagonal lattice, with cell states representing fluid, rock, and grain aggregates that are either stationary or in a state of motion in one of the six cardinal lattice directions. An optional additional state represents unusually large grains (“blocks”) that cannot be displaced upward by disturbance events.

The Grain Hill model is able to reproduce a range of common slope forms, from fully soil mantled to rocky and partially mantled. The bestiary of forms that the model can produce includes convex-upward soil mantled slopes, planar slopes (bare, soil mantled, or in between), and cliffs with basal ramparts. When the model is configured to include a resistant rock layer that decomposes into blocks, the model reproduces observed hogback-like slope forms and qualitatively matches the behavior predicted by a recent continuum–discrete model (Glade and Anderson, 2017; Glade et al., 2017).

In its simplest guise, the model has only three process parameters, which represent disturbance frequency, characteristic disturbance depth, and base-level lowering rate, respectively. Incorporating physical weathering of rock adds one additional parameter, representing the characteristic rock weathering rate. These parameters are not arbitrary but rather have a direct link with corresponding parameters in continuum theory. Comparison between observed and modeled slope forms demonstrates that the model can reproduce both the shape and scale of real hillslope profiles.

Experiments with the Grain Hill model highlight the importance of regolith cover fraction in governing both the downslope mass transport rate and the rate of physical weathering. Equilibrium rocky hillslope profiles are possible even when the rate of base-level lowering exceeds the nominal bare-rock weathering rate, because increases in both slope gradient and roughness can allow for rock weathering rates that are greater than the flat-surface maximum. Finally, experiments in transient relaxation of steep, rocky slopes predict the formation of a regolith-mantled pediment that migrates headward through time while maintaining a sharp slope break.

The Grain Hill model is freely available in an open-source repository (Tucker, 2018b), as is the Landlab toolkit (Hobley et al., 2017; Hutton et al., 2016). Videos showing simulations of mesa and hogback evolution are also available as an online resource (Tucker, 2018a).

The supplement related to this article is available online at: https://doi.org/10.5194/esurf-6-563-2018-supplement.

The idea to develop a 2-D cellular rock-slope model arose from conversations among all three authors. The model code was written in Landlab by GT. Both GT and DEJH contributed to the underlying grid data structures and Python code. SWM extracted the hillslope profiles and estimated the parameters for the two field sites. GT performed the computational experiments and wrote the paper, with input and editing from SWM and DEJH.

The authors declare that they have no conflict of interest.

This research was supported by the US National Science Foundation
(EAR-1349390 and ACI-1450409 to Gregory E. Tucker and Daniel E. J. Hobley, and EAR-1349229 to Scott W. McCoy). Daniel E. J. Hobley's
participation was also supported in part by the National Center for Earth
Surface Dynamics (EAR-1246761). Support for high-performance computing and
software development was provided by the Community Surface Dynamics Modeling
System (CSDMS) (EAR-1226297). High-resolution topographic data were downloaded
from Open Topography (http://www.opentopography.org/, last access: 17 October 2017).

Edited by: Greg Hancock

Reviewed by: Efi Foufoula-Georgiou and Joshua Roering

Ahnert, F.: The role of the equilibrium concept in the interpretation of landforms of fluvial erosion and deposition, Proc. symp. l'évolution des versants (Liége), 50, 23–51, 1967. a, b

Alonso, J. and Herrmann, H.: Shape of the tail of a two-dimensional sandpile, Phys. Rev. Lett., 76, 4911, https://doi.org/10.1103/PhysRevLett.76.4911, 1996. a

Anderson, R. S., Anderson, S. P., and Tucker, G. E.: Rock damage and regolith transport by frost: An example of climate modulation of the geomorphology of the critical zone, Earth Surf. Proc. Land., 38, 299–316, 2012. a

Andrews, D. and Bucknam, R. C.: Fitting degradation of shoreline scarps by a nonlinear diffusion model, J. Geophys. Res., 92, 12857–12867, 1987. a, b

Binnie, S. A., Phillips, W. M., Summerfield, M. A., and Fifield, L. K.: Tectonic uplift, threshold hillslopes, and denudation rates in a developing mountain range, Geology, 35, 743–746, https://doi.org/10.1130/G23641A.1, 2007. a, b

Chen, S. and Doolen, G. D.: Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech., 30, 329–364, 1998. a

Cottenceau, G. and Désérable, D.: Open Environment for 2d Lattice-Grain CA, in: Cellular Automata, ACRI, edited by: Bandini S., Manzoni S., Umeo H., and Vizzari G., Lecture Notes in Computer Science, 6350, Springer, Berlin, Heidelberg, 2010. a

Coyote, W. E.: Fast and furry-ous: exploring the links between gravitational forces and situational awareness, PhD thesis, Acme Technical College, Tombstone, Arizona, USA, 1949. a

Culling, W.: Soil creep and the development of hillside slopes, J. Geol., 71, 127–161, 1963. a, b

Culling, W.: Theory of erosion on soil-covered slopes, J. Geol., 73, 230–254, 1965. a

Désérable, D.: A versatile two-dimensional cellular automata network for granular flow, SIAM J. Appl. Math., 62, 1414–1436, 2002. a

Désérable, D., Dupont, P., Hellou, M., and Kamali-Bernard, S.: Cellular automata in complex matter, Aip. Conf. Proc., 20, 67, 2011. a

Drake, T. G. and Calantoni, J.: Discrete particle model for sheet flow sediment transport in the nearshore, J. Geophys. Res.-Oceans, 106, 19859–19868, 2001. a

Duszyński, F. and Migoń, P.: Boulder aprons indicate long-term gradual and non-catastrophic evolution of cliffed escarpments, Stołowe Mts, Poland, Geomorphology, 250, 63–77, 2015. a

Foufoula-Georgiou, E., Ganti, V., and Dietrich, W.: A nonlocal theory of sediment transport on hillslopes, J. Geophys. Res., 115, F00A16, https://doi.org/10.1029/2009JF001280, 2010. a, b

Furbish, D. and Haff, P.: From divots to swales: Hillslope sediment transport across divers length scales, J. Geophys. Res., 115, F03001, https://doi.org/10.1029/2009JF001576, 2010. a, b

Furbish, D. J. and Roering, J. J.: Sediment disentrainment and the concept of local versus nonlocal transport on hillslopes, J. Geophys. Res.-Earth, 118, 937–952, 2013. a

Furbish, D. J. and Schmeeckle, M. W.: A probabilistic derivation of the exponential-like distribution of bed load particle velocities, Water Resour. Res., 49, 1537–1551, 2013. a

Furbish, D., Hamner, K., Schmeeckle, M., Borosund, M., and Mudd, S.: Rain splash of dry sand revealed by high-speed imaging and sticky paper splash targets, J. Geophys. Res., 112, F01001, https://doi.org/10.1029/2006JF000498, 2007. a

Furbish, D., Haff, P., Dietrich, W., and Heimsath, A.: Statistical description of slope-dependent soil transport and the diffusion-like coefficient, J. Geophys. Res., 114, F00A05, https://doi.org/10.1029/2009JF001267, 2009. a, b, c, d, e, f, g, h, i, j

Gabet, E.: Sediment transport by dry ravel, J. Geophys. Res., 108, 2049, https://doi.org/10.1029/2001JB001686, 2003. a

Gabet, E. J. and Mendoza, M. K.: Particle transport over rough hillslope surfaces by dry ravel: Experiments and simulations with implications for nonlocal sediment flux, J. Geophys. Res., 117, F01019, https://doi.org/10.1029/2011JF002229, 2012. a, b

Ghil, M., Zaliapin, I., and Coluzzi, B.: Boolean delay equations: A simple way of looking at complex systems, Physica D., 237, 2967–2986, 2008. a

Glade, R. and Anderson, R.: Quasi-steady evolution of hillslopes in layered landscapes: An analytic approach, J. Geophys. Res., 123, 26–45, https://doi.org/10.1002/2017JF004466, 2017. a, b, c

Glade, R. C., Anderson, R. S., and Tucker, G. E.: Block-controlled hillslope form and persistence of topography in rocky landscapes, Geology, 45, 311–314, 2017. a, b, c, d, e, f, g, h

Gutt, G. and Haff, P.: An automata model of granular materials, in: Proceedings of the fifth distributed memory computing conference, Charleston, SC, USA, 1990. a

Heimsath, A., Dietrich, W., Nishiizumi, K., and Finkel, R.: The soil production function and landscape equilibrium, Nature, 388, 358–361, 1997. a

Heimsath, A., DiBiase, R., and Whipple, K.: Soil production limits and the transition to bedrock-dominated landscapes, Nat. Geosci., 5, 210–214, 2012. a, b

Hobley, D. E. J., Adams, J. M., Nudurupati, S. S., Hutton, E. W. H., Gasparini, N. M., Istanbulluoglu, E., and Tucker, G. E.: Creative computing with Landlab: an open-source toolkit for building, coupling, and exploring two-dimensional numerical models of Earth-surface dynamics, Earth Surf. Dynam., 5, 21–46, https://doi.org/10.5194/esurf-5-21-2017, 2017. a

Howard, A. D.: A detachment-limited model of drainage basin evolution, Water Resour. Res., 30, 2261–2285, 1994. a, b

Howard, A. D. and Selby, M. J.: Rock Slopes, in: Geomorphology of Desert Environments, 123–172, Springer, Dordrecht, 1994. a

Hutton, E., Hobley, D. E. J., Tucker, G. E., Nudurupati, S. S., Adams, J. M., Gasparini, N. M., Knuth, J. S., Strauch, R., Shobe, C. M., Barnhart, K. R., Rengers, F. K., and Istanbulluoglu, E.: Landlab version 1.0., https://doi.org/10.5281/zenodo.154179, 2016.

Johnstone, S. A., Chadwick, K. D., Frias, M., Tagliaro, G., and Hilley, G. E.: Soil Development over Mud-Rich Rocks Produces Landscape-Scale Erosional Instabilities in the Northern Gabilan Mesa, Geol. Soc. Am. Bull., 129, 1266–79, 2017. a

Károlyi, A. and Kertész, J.: Lattice-gas model of avalanches in a granular pile, Phys. Rev. A., 57, 852, https://doi.org/10.1103/PhysRevE.57.852, 1998. a

Károlyi, A. and Kertész, J.: Granular medium lattice gas model: the algorithm, Comput. Phys. Commun., 121, 290–293, 1999. a

Károlyi, A., Kertész, J., Havlin, S., Makse, H. A., and Stanley, H. E.: Filling a silo with a mixture of grains: friction-induced segregation, Europhys. Lett., 44, 386, 1998. a

Lamb, M. P., Scheingross, J. S., Amidon, W. H., Swanson, E., and Limaye, A.: A model for fire-induced sediment yield by dry ravel in steep landscapes, J. Geophys. Res.-Earth, 116, F03006, https://doi.org/10.1029/2010JF001878, 2011. a

Lamb, M. P., Levina, M., DiBiase, R. A., and Fuller, B. M.: Sediment storage by vegetation in steep bedrock landscapes: Theory, experiments, and implications for postfire sediment yield, J. Geophys. Res.-Earth, 118, 1147–1160, 2013. a, b

MacVicar, B., Parrott, L., and Roy, A.: A two-dimensional discrete particle model of gravel bed river systems, J. Geophys. Res.-Earth, 111, F3, https://doi.org/10.1029/2005JF000316, 2006. a

Martinez, J. and Masson, S.: Lattice grain models, in: Silos, edited by: Brown, C. and Nielsen, J., London, CRC Press, 1998. a

McEwan, I. and Heald, J.: Discrete particle modeling of entrainment from flat uniformly sized sediment beds, J. Hydraul. Eng., 127, 588–597, 2001. a

Narteau, C., Le Mouël, J., Poirier, J., Sepúlveda, E., and Shnirman, M.: On a small-scale roughness of the core–mantle boundary, Earth Planet. Sc. Lett., 191, 49–60, 2001. a

Narteau, C., Zhang, D., Rozier, O., and Claudin, P.: Setting the length and time scales of a cellular automaton dune model from the analysis of superimposed bed forms, J. Geophys. Res.-Earth, 114, F03006, https://doi.org/10.1029/2008JF001127, 2009. a

Peng, G. and Herrmann, H. J.: Density waves of granular flow in a pipe using lattice-gas automata, Phys. Rev. A., 49, R1796, https://doi.org/10.1103/PhysRevE.49.R1796, 1994. a

Perron, J. T., Kirchner, J. W., and Dietrich, W. E.: Formation of evenly spaced ridges and valleys, Nature, 460, 502–505, 2009. a

Perron, J. T., Richardson, P. W., Ferrier, K. L., and Lapôtre, M.: The Root of BBranching River Networks, Nature, 492, 100–103, 2012. a, b

Roering, J.: Soil creep and convex-upward velocity profiles: Theoretical and experimental investigation of disturbance-driven sediment transport on hillslopes, Earth Surf. Proc. Land., 29, 1597–1612, 2004. a

Roering, J.: How well can hillslope evolution models “explain” topography? Simulating soil transport and production with high-resolution topographic data, Geol. Soc. Am. Bull., 120, 1248–1262, 2008. a

Roering, J., Kirchner, J., and Dietrich, W.: Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology, Water Resour. Res., 35, 853–870, 1999. a, b

Roering, J., Kirchner, J., Sklar, L., and Dietrich, W.: Hillslope evolution by nonlinear creep and landsliding: An experimental study, Geology, 29, 143–146, 2001. a

Roering, J. J. and Gerber, M.: Fire and the evolution of steep, soil-mantled landscapes, Geology, 33, 349–352, 2005. a

Rozier, O. and Narteau, C.: A real-space cellular automaton laboratory, Earth Surf. Proc. Land, 39, 98–109, 2014. a

Schmeeckle, M. W.: Numerical simulation of turbulence and sediment transport of medium sand, J. Geophys. Res.-Earth, 119, 1240–1262, 2014. a

Shobe, C. M., Tucker, G. E., and Anderson, R. S.: Hillslope-derived blocks retard river incision, Geophys. Res. Lett., 43, 5070–5078, 2016. a

Small, E., Anderson, R., and Hancock, G.: Estimates of the rate of regolith production using 10Be and 26Al from an alpine hillslope, Geomorphology, 27, 131–150, 1999. a

Tucker, G. and Bradley, D.: Trouble with diffusion: Reassessing hillslope erosion laws with a particle-based model, J. Geophys. Res., 115, F1, https://doi.org/10.1029/2009JF001264, 2010. a, b, c

Tucker, G. E., Hobley, D. E. J., Hutton, E., Gasparini, N. M., Istanbulluoglu, E., Adams, J. M., and Nudurupati, S. S.: CellLab-CTS 2015: continuous-time stochastic cellular automaton modeling using Landlab, Geosci. Model Dev., 9, 823–839, https://doi.org/10.5194/gmd-9-823-2016, 2016. a, b, c, d, e, f, g

Tucker, G. E.: GrainHill cellular hillslope model: GIF animations of hillslope evolution, https://doi.org/10.6084/m9.figshare.6720476, 2018a.

Tucker, G. E.: GrainHill version 1.0, https://doi.org/10.5281/zenodo.1306961, 2018b.

Zhang, D., Narteau, C., and Rozier, O.: Morphodynamics of barchan and transverse dunes using a cellular automaton model, J. Geophys. Res., 115, F3, https://doi.org/10.1029/2009JF001620, 2010. a

Zhang, D., Narteau, C., Rozier, O., and du Pont, S. C.: Morphology and dynamics of star dunes from numerical modelling, Nat. Geosci., 5, 463–467, 2012. a