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.

Table 1States in the Grain Hill model.

Download Print Version | Download XLSX

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.

Figure 4Rules for motion and frictional (inelastic) collisions, illustrated here for one of the six lattice directions.

Download

Figure 5Illustration of gravitational rules. The bottom row shows the “falling on a slope” rule, which effectively imposes a 30^∘ angle of repose. Modified from Tucker et al. (2016).

Download

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.

Figure 6Lattice grain simulation of emptying of a silo. Light-shaded grains are stationary; darker-shaded ones are in motion. Black cells are walls (rock). Time units are indicated in seconds. From Tucker et al. (2016).

Download

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.

Figure 7Transitions representing rock-to-regolith transformation by weathering (a) and regolith disturbance (b), in which a stationary particle becomes mobile and switches position with a air cell. The illustration represents one of the six possible orientations.

Download

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 δ³.

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⁻¹ 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).

Figure 8Equilibrium topographic cross sections using only regolith particles (no rock) and a variety of disturbance frequencies (d) and time interval between base-level fall events (τ). Fast basal incision and/or infrequent disturbance lead to planar threshold hillslopes; slow basal incision and/or frequent disturbance lead to parabolic hillslopes.

Download

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.

Figure 9Dimensionless mean hillslope height, h, as a function of dimensionless disturbance rate d^′ for a range of hillslope lengths. Data points include 125 sensitivity analysis runs in which $d\in [{\mathrm{10}}^{-\mathrm{3}},{\mathrm{10}}^{-\mathrm{2.5}},{\mathrm{10}}^{-\mathrm{2}},{\mathrm{10}}^{-\mathrm{1.5}},{\mathrm{10}}^{-\mathrm{1}}]$, $\mathit{\tau }\in [{\mathrm{10}}^{\mathrm{2}},{\mathrm{10}}^{\mathrm{2.5}},{\mathrm{10}}^{\mathrm{3}},{\mathrm{10}}^{\mathrm{3.5}},{\mathrm{10}}^{\mathrm{4}}]$, and λ as shown in the legend.

Download

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=2d.

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{‾}}{\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), δ²:

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

where $\stackrel{\mathrm{‾}}{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²θ in Eq. (12) approaches unity, and using the prior relation N_a=2d, 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}{\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² yr⁻¹, 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⁻¹.

Figure 10Relationship between dimensionless diffusivity and mean gradient, from the series of 125 model runs of which a subset is shown in Fig. 8.

Download

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.

Figure 11Final equilibrium profiles from Grain Hill runs with rock and weathering. Domain size is 222 rows by 257 columns, and uplift interval ranges from 100 to 10 000 years.

Download

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.)

Figure 12Mean equilibrium gradient and regolith thickness for models with rock and weathering, as a function of d^′ and w^′. Data represent 125 runs with $d\in [{\mathrm{10}}^{-\mathrm{3}},{\mathrm{10}}^{-\mathrm{2.5}},{\mathrm{10}}^{-\mathrm{2}},{\mathrm{10}}^{-\mathrm{1.5}},{\mathrm{10}}^{-\mathrm{1}}]$ yr⁻¹, $w^{\prime }\in [{\mathrm{10}}^{\mathrm{0}},{\mathrm{10}}^{\mathrm{0.5}},{\mathrm{10}}^{\mathrm{1}},{\mathrm{10}}^{\mathrm{1.5}},{\mathrm{10}}^{\mathrm{2}}]$, and $\mathit{\tau }\in [{\mathrm{10}}^{\mathrm{2}},{\mathrm{10}}^{\mathrm{2.5}},{\mathrm{10}}^{\mathrm{3}},{\mathrm{10}}^{\mathrm{3.5}},{\mathrm{10}}^{\mathrm{4}}]$ yr. Horizontal dashed lines show model's angle of repose for regolith.

Download

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.

Figure 13Comparison between rate of material input, 1∕τ (cells ∕ year), with effective rate of weathering, wn_ra∕n_H, from 125 model runs (see text).

Download

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₀ 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{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⁻¹, 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.

Figure 14Regolith thickness versus time, as predicted by inverse-exponential theory (log growth; solid cyan curve) and the Grain Hill model with a range of ratios of disturbance rate (d) to weathering rate (w). Time (horizontal axis) is non-dimensionalized by multiplying by w.

Download

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.

Figure 15Two examples of cliff-and-rampart morphology. (a) Near Palisade, Colorado, USA, after recent rock-fall event (photo courtesy of D. Nathan Bradley and Dylan Ward). (b) Colorado Plateau, Utah, USA. Note that contact between lower rampart and subvertical slopes, both of which have formed in a gray shale unit, occurs without any apparent break in lithology.

Download

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.

Figure 16Quasi-steady model hillslope profiles created using a collapse rule, under four different combinations of d^′ and w^′. Insets show magnified views of a portion of each hillslope.

Download

Figure 17Time series showing transient erosion of a steep rock slope under a stable base-level, highlighting formation of ramp-and-cliff morphology. Simulation shows 20 000 years of slope evolution under $d=w={\mathrm{10}}^{-\mathrm{3}}$ yr⁻¹. Nominal width, assuming δ=0.1 m, is 12 m.

Download

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.

Figure 18Examples of models that include blocks. Rock (black) weathers to blocks (dark red), which can only move downward or downward plus laterally. Blocks in turn weather to regolith (light brown) (GIF-format animations of similar runs are available as an online resource; see Tucker, 2018a).

Download