Glassy dynamics of landscape evolution.

Soil creeps imperceptibly downhill, but also fails catastrophically to create landslides. Despite the importance of these processes as hazards and in sculpting landscapes, there is no agreed-upon model that captures the full range of behavior. Here we examine the granular origins of hillslope soil transport by discrete element method simulations and reanalysis of measurements in natural landscapes. We find creep for slopes below a critical gradient, where average particle velocity (sediment flux) increases exponentially with friction coefficient (gradient). At critical gradient there is a continuous transition to a dense-granular flow rheology. Slow earthflows and landslides thus exhibit glassy dynamics characteristic of a wide range of disordered materials; they are described by a two-phase flux equation that emerges from grain-scale friction alone. This glassy model reproduces topographic profiles of natural hillslopes, showing its promise for predicting hillslope evolution over geologic timescales.

Creep in soil-mantled hillslopes takes place due to slow rearrangements of grains that can be reasonably approximated as a viscous flow (1). At critical slopes or under certain perturbations like rain events, however, soil may fail catastrophically to create landslides (2). These processes control the erosion and form of hillslopes and the delivery of sediment to rivers (3). Despite widely varying materials and environments, hillslope soil motion falls into two distinct categories: slowly creeping “earthflows” associated with surface velocities that cover 10 orders of magnitude up to m/s and rapid landslides (including debris flows, mudflows, etc.) that are faster than m/s (4–7) (Fig. 1). Much progress has been made in mechanistic models for the latter; in particular, continuum models based on mass and momentum conservation for the granular and fluid phases are able to reproduce important aspects of soil failure and mass-movement runout (2, 8). Models for hillslope soil creep, however, lack a mechanistic underpinning. For over 50 y, a heuristic “diffusive-like law”—in which sediment flux is proportional to topographic gradient —has been used to model landscape erosion (9–11). For soil to creep at subcritical gradients, it is supposed that dilation occurs as a result of (bio)physical perturbations such as freeze/thaw/swell, rainsplash, tree throw, and burrowing animals (12–14).

Fig. 1.

Landslide and creep phenomenology. (A) Rapid landslide in San Salvador, El Salvador. Image courtesy of Wikimedia Commons/USGS. (B) Slow earthflow in Serranía del Interior, Venezuela. Image courtesy of Wikimedia Commons/Fev. (C) Ranges of surface velocities observed for various types of slow and rapid landslides. The datapoints in red, brown, magenta, and green correspond to the observations reported or documented by Cruden and Varnes (4), Hungr et al. (5), Hilley et al. (6), and Saunders and Young (7), respectively. (D) Schematic cross-section of a soil-mantled hillslope.

Hillslope soil creep has not been connected to the creep phenomenon observed in diverse amorphous and granular systems with a wide range of materials and particle shapes, including dry (15) and fluid-driven (16, 17) granular flows. These materials belong to a large class of systems—including glasses, pastes, foams, gels, and suspensions—that are known to exhibit interesting rheological properties termed “glassy dynamics” that include slow dynamics such as compaction (18), hysteresis and history dependence of the static configurations (19), and intermittency and spatially heterogeneous dynamics (20). These behaviors are thought to be a natural consequence of two properties shared by all glassy materials: structural disorder and metastability (21). In such materials, thermal motion alone is not enough to achieve complete structural relaxation, and consequently relaxation times are extremely large compared with the timescale of a typical experiment. As a result and for practical purposes, glassy materials are nonequilibrium systems with long memory (22). For the evolution of soil-mantled landscapes over geologic time, such glassy dynamics should be relevant.

Dense granular flows on inclined planes are a reasonable idealized model for hillslope soil transport, and the rheology of such flows has been thoroughly examined in simulations and experiments (23–31). Local (32, 33) and then nonlocal (34–36) constitutive relations developed for dense, steady flows establish that the effective friction is a nonlinear function of the (nondimensionalized) shear rate. We note, however, that most of these studies assume that granular piles are jammed below the angle of repose; i.e., they do not examine subcritical creep dynamics (for an exception see ref. 37). Creep involves grain-scale rearrangement—due to structural and mechanical disorder in the pack (31, 38)—that induces exponentially small but finite particle velocities below threshold (39). Although creep has been recognized for some time (15), only recently have researchers begun probing the nature of the creep to fluid transition in granular heap flows. Experiments performed with an inclined granular layer show localized and isolated events—microfailures—in the bulk at inclinations below the bulk angle of repose. As the inclination increases, microfailures occur more frequently until they coalesce to form an avalanche (40). It has been suggested that this change from creeping to flowing states is a dynamical phase transition (41) that is phenomenologically similar to a liquid–glass transition (40, 42, 43). In the athermal granular system, the transition occurs in the vicinity of a critical force (instead of the critical temperature for thermal glasses). One way to probe this transition is to define an order parameter that describes a dynamical quantity that changes dramatically across the phase transition. Recent simulations have considered the particular case of the onset of erosion of a single layer of monodisperse grains on a substrate; Yan et al. (44) and Aussillous et al. (45) showed that this can be mapped to a plastic depinning transition, a generic phase transition associated with irreversible strain and disorder (42). In their work, sediment flux (current) works as an order parameter that goes sharply to zero at a critical driving force below which flow does not occur (44, 45). In the presence of additional noise, however—which may take the form of internal disorder, sidewalls/boundary effects, or external perturbations—localized flow may occur below the critical force in amorphous systems (34, 46–48). Experiments have demonstrated that the onset of erosion in 3D granular flows is indeed continuous, accompanied by subcritical flow and creep (16), and exhibits dynamics qualitatively consistent with a liquid–glass transition (40); however, the dynamical nature of this transition has not yet been examined in such systems. In this paper we conduct 3D numerical experiments of granular heap flows to first document subthreshold creep and then demonstrate that the creep to landslide transition is continuous and is quantitatively consistent with a plastic depinning transition. We then present evidence of such glassy dynamics in the soil transport rates and topographic profiles of real hillslopes in nature.

Granular Hillslope Model

Building on the work described above, and recent success in linking idealized granular-physics models to sediment transport (49), we examine the granular origins of a creep transport equation and the transition to landsliding by developing a granular hillslope model using the discrete element method (DEM) (Fig. 2). Simulations are performed with LAMMPS (lammps.sandia.gov). In typical soil-mantled hillslopes there is a meters-deep, mobile surface layer composed of mostly unbonded particles that are sand sized and smaller; it is underlain by weathered bedrock (saprolite), composed of increasingly more bonded particles and larger rock fragments, and eventually by unaltered bedrock (50) (Fig. 1). The surface layer may flow as a sheet or as a channel, depending on confinement and material factors. We model an idealized representation of this mobile surface soil: a layer of polydisperse spheres with diameter range , average diameter , and depth , overriding an immobile and frictional bottom boundary (Fig. 2); more details of the model implementation are available in . The system is inclined at various gradients that are below, near, and above the bulk angle of repose (; ). From instantaneous particle motions, we compute vertical profiles of time- and horizontal ()-averaged downslope grain velocity, , over the duration of each model run (see for averaging procedure).

Fig. 2.

Hillslope DEM simulations. (A–C) Snapshots correspond to time (s) after inclining the granular hillslope; particle colors represent their downslope (-dir) instantaneous velocities, while black lines show time-averaged downslope velocities. The bulk (macroscopic) angle of repose for this set of simulations is . (A) corresponds to a hillslope just below onset of dense flow at the surface, where the pack is almost fully creeping. (B) is right at the transition point. (C) shows a fully developed dense granular flow in more than half of the model depth. D and E show time-averaged downslope (-dir) velocity () and local friction coefficient () profiles, respectively, for the granular hillslope at . The critical depth and critical downslope velocity at the transition to creep are indicated. The value is further used in E to determine the critical friction coefficient .

To study the dynamical behavior of our model hillslope in the phase transition framework described above, we calculate the per-grain friction coefficient , where is the per-grain stress tensor in the plane and is the average confining/normal stress (). Similar to velocity, we perform time and horizontal averaging to produce vertical profiles of the local friction coefficient for each run; since only the values averaged over time and the x direction are reported, we retain the notation to emphasize that measured friction is local in that it changes with depth. Our simulations show two distinct phenomenological behaviors for inclinations below and above . Below we observe slow particle velocities, a slow decay rate of the downslope particle velocity () with depth, and hot spots of intermittent and localized motion (51). We interpret this regime as creep (Fig. 2), where similar dynamics have been reported in experiments (15, 16, 40). As crosses we observe the abrupt emergence of a fast and continuous surface layer flowing over the creeping regime, whose velocity and thickness increase with increasing above critical friction coefficient (Fig. 2). The rapid decay of with depth and the continuous nature of the particle motions are consistent with the dense granular flow regime (15, 16). We verified that the DEM simulations reproduce the rheology observed in heap-flow experiments that are similar to our model setup (). The transition from dense-granular flow to creep is associated with a kink in the mean velocity profile, which is used to define a critical depth () and critical particle velocity () for each inclination modeled (Fig. 2). We also observe a kink in the profile of friction coefficient with depth, which occurs at the same critical depth; we infer the associated value as the critical friction associated with the creep transition, (Fig. 2), and compute this also for each inclination. We emphasize here that creeping and dense flow regimes can take place in the same column of soil, with dense flow at the top and creep at the bottom. The variation of downslope particle velocity () vs. friction coefficient is qualitatively similar for four different inclinations, below and above the bulk angle of repose (Fig. 3). These observations, and previous work (44), suggest the possibility of a generalized relation between particle velocity and local friction. Comparison of normalized average grain velocity, , vs. normalized local friction coefficient, , for all inclinations confirms this notion and reveals a striking pattern (Fig. 3). For the creep regime () simulations show an exponential flow relation, , with a transition friction coefficient of that is similar for all model runs. For the dense granular flow regime () the functional form is a power law, , where is the critical exponent, the critical point, and is the order parameter. This is in agreement with the suggestions by Fisher (52) and Chauve et al. (39) that the pinned to sliding transition is a second-order phase transition in which the order parameter obeys power-law scaling close to the critical point. A similar depinning transition with an exponential-law relation at a low driving force and a power-law relation above a critical threshold has been also reported in failure of inhomogeneous brittle materials (53). For a plastic depinning transition it is expected that the critical exponent (42), consistent with our simulation results. A precise value for the exponent cannot be determined from numerical results alone. Moreover, finite-size effects and dimensionality of our system may influence the value of the critical exponent. Nevertheless, a representative value , shown for illustration purposes (Fig. 3), is in the general range reported for colloidal and granular systems near the critical point (16, 54, 55).

Fig. 3.

General flow behavior in simulations. (A) DEM results showing local downslope velocity () as a function of local friction coefficient () for four different inclinations, below and above the bulk angle of repose. (B) DEM results showing normalized local downslope velocity () as a function of normalized local friction coefficient () for the four different inclinations shown in A. Dashed line illustrates an exponential scaling for creep regime with critical gradient and a power-law scaling for the range of large flux with a power-law exponent . Relations are not fitted to the data; they are shown only for illustrative purposes.

Hillslope Evolution as a Depinning Transition: Field Evidence

We acknowledge that our highly idealized simulations may not translate to downslope movement of heterogeneous soil in the complex natural environment. The above findings indicate, however, that hillslope soil movement may be governed—at least in part—by generic glassy dynamics associated with disordered granular systems. To search for signatures of this behavior in the field, we have collected observations of sediment flux () as a function of local hillslope gradient from four different published sources with different climatic and uplift conditions, as well as different soil types and material properties (). Unfortunately, data collection methods and their associated timescales differ among these studies, potentially contributing to noise and limiting our ability to directly compare among these different field sites. The data presented by Yoo et al. (56) are from Frog’s Hollow, a semiarid eucalyptus grassland savannah hillslope located about 80 km south–southeast from Canberra, NSW, Australia. They calculated sediment transport rates by taking soil samples along a hillslope transect, measuring the soil mass production rate and the elemental chemistry of soils and saprolite, and then using an iterative modeling process integrating over the timescales of chemical weathering. Data presented by Gabet (57) are from the Santa Ynez Valley in the tectonically active transverse ranges near Santa Barbara, CA, with a semiarid Mediterranean climate. The measurements were carried out over annual timescales using sediment traps installed on hillslopes. Data by Martin and Church (58) and Martin (59) are from the Queen Charlotte Islands, off the coast of British Columbia in Canada, with an oceanic climate and extremely frequent precipitation. The original datasets in these studies are from reports by Rood (60, 61), which are landslide inventories completed by identifying landslides on aerial photographs in the region. They cover approximately a 40-y period of activity.

In studies by Yoo et al. (56) and Gabet (57), the sediment flux was originally measured in units of mass flux for each basin. The averaged bulk density of soil in the study by Yoo et al. (56) is reported as 1,800 and in the study area by Gabet (57) is reported as 1,770 . We used these densities to convert mass sediment flux to volumetric sediment flux (Fig. 4 ). The studies by Martin (59) and Martin and Church (58) originally reported measurements of volumetric sediment flux. They presented measurements for 23 different basins: Thirteen sites (group A) have the majority of their basins composed of soft volcanic and sedimentary rocks, while the 10 remaining sites (group B) have a larger proportion of their basins composed of hard volcanic rocks and granites. Here, we analyze the data from their group A and group B hillslopes separately to calculate critical slope and flux at that critical slope. We further calculated the mean value of sediment flux measurements for each hillslope gradient class in each group; these are presented in a single sediment flux vs. hillslope gradient relationship for the studied region (Fig. 4 ). We assumed a constant error bar of 0.1 for all gradient measurements (classes) in the studies by Martin (59) and Martin and Church (58), even though such an error in measurements is not explicitly stated in their studies. All data here were extracted from figures in the cited papers and also reports cited therein. For the case of studies by Martin (59) and Martin and Church (58), the values of sediment flux and their error bars at slopes smaller than 20° are calculated based on the estimates presented in table 1 in Martin and Church (58).

Fig. 4.

Field data showing measured sediment flux vs. hillslope gradient for natural hillslopes reported previously in the literature. (A) Gabet (57). (B) Group A basins in Martin (59) and Martin and Church (58). (C) Group B basins in Martin (59). (D) Yoo et al. (56). (E) Rescaled data combined from A–D. Note that normalized flux is equivalent to normalized velocity while normalized gradient () is equivalent to normalized friction, allowing comparison with numerical results (Fig. 3). The dashed line in E shows the same bipartite flux relation as in Fig. 3 for illustration purposes: i.e., an exponential flux relation for gradients below critical gradient and a power-law relation for larger gradients.

For each study we observe a kink in the relation between flux and slope, which allows us to determine critical values of and for each field site by eye from inflection points in the plots (). These critical values are illustrated with arrows in Fig. 4 . The normalized sediment flux and the normalized hillslope gradient are then calculated as and , respectively. Note that sediment flux , where the angle brackets indicate averaging over the depth of the soil column (). While studies report values for , most are derived from surface measurements and therefore mostly reflect surface velocities of the flows, while flow depths are generally poorly constrained. Normalization removes this depth dependence, however, since . Hillslope gradient is also related to friction coefficient; assuming a naive hydrostatic and 1D behavior for a depth-averaged earth flow, , where is acceleration due to gravity. Thus, data plots of normalized flux and slope for the field data are equivalent to normalized velocity and friction presented from numerical simulations.

The critical values for slope and flux likely encode climatic, tectonic, and soil properties unique to each site (62). Plotting normalized values from all field sites as against , however, collapses the data and results in a pattern that is similar to our model results of against (Fig. 4). We take this similarity as strong evidence for a generic depinning transition, where field data are consistent with an exponential flux relation for gradients below critical gradient and a power-law relation for larger gradients. Moreover, the values inferred for critical slopes at each field site are physically meaningful; they correspond to a reasonable range of reported values for the angle of repose (or friction coefficient) of soils () (63, 64). The scatter in the data, especially above the critical gradient, can be due to several factors, including (i) limited sediment availability on natural hillslopes, because high slopes transport soil faster than it can be produced from weathered bedrock, and (ii) different measurement methods used in different studies, including variations in the detection limits for transport.

Some additional features of the data warrant mention, in terms of their physical interpretation. Hillslope creep velocities (fluxes) are measured over years to decades and thus average over event-based and seasonal fluctuations in flow speed that often occur due to precipitation and temperature effects (7, 58), bioturbation (14), and other disturbances. The timescales associated with measured velocities (fluxes) for above-critical flows are also longer than those of the individual landslides that (presumably) occur. Thus, the reported velocities (fluxes) of soil motion should be understood as the long-term average of episodic and relatively fast events and intervening periods of relatively slow motion. This is not unlike mountain uplift that results from repeated fault slip; average uplift rates are not representative of slip events, but are nonetheless meaningful for considering landscape erosion that occurs over geologic timescales (65). Another notable feature is that fluxes vary widely for slopes slightly to moderately larger than critical gradient (). We speculate that flows within this slope range may occur as either creep or landslides, depending on environmental forcing [e.g., pore pressure (2)], supply (or availability) of material, and soil thickness. As a result, the creep-to-landsliding transition in natural landscapes is much more variable than in a constant forcing situation such as our model. We suspect that averaging over suitably large (geologic) timescales, if possible, would recover a flux–slope relation that is similar to model expectations but with a more diffuse flow transition. Such an idea may be tested by examining the topographic hillslope profile produced by erosion and rock uplift over geologic timescales (10).

Modeling Landscape Evolution with a Glassy Flux Model

We propose a bipartite “glassy flux model” to represent hillslope soil transport that joins the exponential and power-law relations associated with creep and landsliding regimes, respectively, aswhere is the Heaviside step function that acts to blend the two transport regimes across the transition (39), and is a constant associated with a particular field site. Eq. implicitly assumes steady flow and therefore is applicable only for long timescales that integrate over very many flow events (). The flux Eq. is related to hillslope erosion through conservation of mass,where and are the bulk densities of sediment and rock, respectively, is the rate of landscape elevation change, and is the rock uplift rate. We use Eqs. and to model the steady-state form of a hillslope, i.e., the topography associated with a balance between uplift and erosion () that results from the new flux equation. The left boundary condition of the model hillslope is no flux, representing a drainage divide, and the right boundary condition is a fixed elevation that represents base level. We calibrate and compare model results to hillslope topography data in the Oregon Coast Range (OCR). First, we assume a value for the critical exponent for simplicity because it cannot be better constrained from the data, and then we estimate the values for critical gradient and from values of sediment flux vs. gradient reported by Roering et al. (10) from a site near Coos Bay (). Next, we model the transient evolution of hillslope topography by iteratively solving Eqs. and starting from a flat initial condition with constant uplift rate and densities determined from Roering et al. (10) for the OCR. The profile is evolved for 20 million y, roughly the time that the OCR has been uplifting (66), and we verify that the hillslope reaches a steady-state topography over this time.

We compare model results to hillslope topographic profiles extracted from aerial light detection and ranging (LiDAR) data at an OCR site 60 km west of Eugene, OR (). The glassy flux model reasonably captures elevation and gradient profiles for both short hillslopes (distance about 40 m) where gradients are mostly below critical gradient and longer hillslopes which significantly exceed the critical slope (Fig. 5). Hillslope gradient profiles show a clear kink at the critical slope value derived from the glassy flux model (Fig. 5); the corresponding angle of 26° represents a reasonable value for the transition from creep to landsliding. The explicit incorporation of landslide (dense-granular flow) dynamics allows the glassy flux model to reproduce the flattening out of hillslope profiles as they lengthen; this flattening has been previously reported and cannot be reproduced with diffusion-like flux equations (68) (). We also verified that the model reproduces observed hillslope topography in a different climatic and geologic setting in California ().

Fig. 5.

Hillslope topography of the OCR derived from publicly available airborne LiDAR data (67). (A) Regional perspective view, showing locations of two example hillslopes. (B and C) The elevation–distance (B) and gradient–distance (C) relationships for representative profiles of hillslopes 1 (black dots) and 2 (red dots) in A. Blue dashed line is the prediction of the “glassy” flux model with and . See for more examples.

Discussion

We have developed a model for hillslope soil transport that describes behaviors from creep and slow-earthflow to landsliding and fast-flow regimes. Although the DEM simulations are highly idealized, we suggest that the underlying dynamics are general. While soil creep on hillsides has been viewed as the result of external perturbations (1, 10, 14), model results show that this is not necessary. We found that the addition of perturbations, through random noise added to the locations of some grains, increased the flux magnitude in the creeping regime but did not change the functional form of the flux–slope relation (). Noise had little influence on the fast-flow regime (). We also confirmed that changing grain shape does not change the qualitative behavior in simulations, by replacing spheres with elongated particles having a 3:1 aspect ratio (). More broadly, observations provide strong evidence that the creep–landslide transition exhibits glassy dynamics that may be modeled as a plastic depinning transition at a critical normalized force. The subcritical exponential relation and the supercritical power-law scaling (Eq. ) are expected behaviors based on theoretical and experimental studies of amorphous systems (42). In other words, such behavior is the generic consequence of dynamical phase transitions in disordered materials. A dimensionless number, local normalized friction coefficient in the numerical model, , is calculated in terms of the effective tangential and normal stresses. Although the current model does not account for changes in fluid pore pressure that are known to influence downslope flow velocity (2), we suggest this effect might be viewed as a perturbation to the effective pressure terms that could be combined in the future with the framework we propose. The effective critical friction coefficients determined for the different field settings examined here fall in the range expected for soil mixtures (63, 69).

These results show how recent advances in the physics of disordered materials can be used to explain the evolution of natural landscapes over geologic timecales. The functional form of the flux equation used in this work is a specific case of a more general form , where represents mechanical internal and external noise (48). We suggest that creep, and its associated slow subcritical flow, takes place in our numerical system and in natural hillslopes due to (i) internal disorder of the particulate packing and (ii) bedrock and saprolite boundary layers that surround the mobile regolith, which continuously inject disorder that may induce creep through nonlocal effects (34, 46, 47). It is an open question for soil-mantled hillslopes whether, and under what conditions, the injection of porosity and noise from external perturbations (plants/animals, freeze/thaw/swell, etc.) produces distinctly different dynamics from the sources of disorder considered here.

Finally, we note that one of the hallmarks of granular and amorphous materials is the emergence of rate weakening in the vicinity of their dynamical phase transition (31, 70). Although the fundamental mechanisms of this phenomenon remain to be explored, we believe the picture provided here can help to understand the origins of rate- and state-dependent friction behavior that has recently been proposed to characterize slow and fast landslides (71, 72).

Materials and Methods

Details of the implementation of the DEM model are described in , and the protocol for calculation of velocity and stress profiles from DEM simulations is provided in . The influence of perturbation on the behavior of the DEM model in presented in . We compare the observations from our DEM simulations with a local granular rheology model applied to a heap-flow granular experiment in . The influences of grain shape on the slow creep regime in DEM simulation and experiments are discussed in . Implementation of the landscape evolution model is described in . Details for measurements of hillslope profiles from LiDAR data are described in . The collections of sediment flux and hillslope surface gradients in the study by Roering et al. (10) are described in .