Slopes in steady-state soil-mantled landscapes tend to increase downslope in a way that balances local transport capacity with the sediment supplied from progressively larger source areas. Most formulations of sediment transport due to hillslope processes scale transport rate with local slope, which produces convex-up forms that are independent of the properties of the underlying lithologies. In this study, we document soil-mantled hillslopes that show variations in slope that mimic the underlying stratigraphy. We present stratigraphic and soil-thickness measurements, topographic analyses, and numerical models to demonstrate that variations in rock type can impact the forms of these soil-mantled hillslopes if hillslope transport rates scale with local slope and soil thickness. This demonstrates that hillslope forms in soil-mantled landscapes can be influenced by the underlying lithology through a coupling between the processes that produce soil from rock and those that transport this soil downslope.


This study investigates how the geologic substrate impacts landscape forms on soil-mantled hillslopes. Gilbert (1880) proposed that when weathering can keep pace with rates of erosion, a mantle of soil buffers surface morphologies from the structure of the subsurface geology. Gilbert (1909) then reasoned that convex forms characteristic of soil-mantled hillslopes balance the mass of mobile material produced along them with that transported by processes whose rates increase with slope. This model of slope-dependent hillslope transport was later formalized as a hillslope diffusion rule (geomorphic transport law, GTL; sensuDietrich et al., 2003) (Culling, 1960). Previous studies have also proposed that transport rates may depend on slope and soil thickness (Ahnert, 1976; Braun et al., 2001; Anderson, 2002; Roering, 2008; Furbish et al., 2009), and some empirical evidence supports this view (Heimsath et al., 2005; West et al., 2014). Soil thickness may reflect some properties of underlying rocks and rates of erosion (Ahnert, 1976), because rates of soil production are expected and observed to depend on soil thickness and lithology (Heimsath et al., 2005). As a result, variations in underlying rocks may impact soil thickness, and in turn, transport rates. Therefore, changes in transport rates for a given slope that may accompany lithologic contacts should be expressed in the form of soil-mantled hillslopes. If this mechanism is in operation, Gilbert’s (1880) proposition that a soil mantle shields topographic forms from underlying lithologic variations needs revisiting.

Hillslopes in the northern Gabilan Mesa, in the central California Coast Ranges (western USA), display the characteristic convex-up morphology predicted by Gilbert (1909) and existing GTLs. However, superimposed on this convexity are regular undulations in slope termed “shadow beds”—subhorizontal, soil-mantled features that traverse hillslopes and valleys (Dohrenwend, 1974) (Fig. 1). Individual shadow beds seem to correlate with the underlying stratigraphy of the Pancho Rico Formation, suggesting that geologic structure is impacting hillslope forms in this soil-mantled landscape.

We present field surveys of soil thickness, topographic analyses, and stratigraphic measurements that establish the connection between bedrock geology and the morphology of a soil-mantled hillslope. Using numerical simulations of landscapes and modeling field data, we demonstrate that shadow bedding is consistent with the quasi-equilibrium form of hillslopes traversing rocks with varying resistances to disaggregation given a soil thickness–dependent transport rule. We propose a modified rule for hillslope soil flux that casts flux as the integral of a soil velocity profile that decays exponentially with depth. This velocity profile shows similar dependence on depth as the activity of a number of processes thought to be responsible for soil transport.



Vertically measured soil thickness, Hv [L], evolves according to the local input of soil by production, which depends on the surface normal soil thickness, H [L], and the divergence of the flux of soil per unit width, forumla [L2 t–1] (e.g., Dietrich et al., 2003): 
The first term on the right-hand side describes the empirically calibrated rate of soil production. W0 [L t–1] is the bedrock lowering rate at zero soil thickness, α [L–1] is a constant describing the decay in production rate with increasing soil thickness, ρbr and ρs [M L–3] are the densities of bedrock and soil, respectively, and θ is the hillslope angle. Theoretical expectations (Ahnert, 1976) and empirical evidence (Heimsath et al., 2005) suggest that lithology controls rates of soil production through W0, while α is fairly constant across different climates and rock types, with a value of ∼2 m–1. Lithology-dependent variations in W0 require different soil thicknesses to achieve equivalent soil production rates over different lithologies. Adding rock uplift, U [L t–1], and solving for the evolution of surface elevations, z [L], gives 
We assume that downslope soil velocities, V(h) [L t–1], reach a maximum at the surface and decay exponentially with depth throughout the entirety of the soil column. This is similar to previously developed expressions (Kirkby, 1967; Anderson, 2002; Roering, 2008; Furbish et al., 2009), but ignores the possibility of an inactive region in deep soil columns: 
The surface normal depth at a point in the soil column is h [L], dc [L] is the scaling depth of the velocity profile, and V0 [L t–1] is the surface velocity. Integrating Equation 3 over H yields the flux per unit width: 

As dc increases relative to H, a greater fraction of the soil column is present in the high-velocity portion of the profile, and flux becomes increasingly sensitive to variations in thickness. If dc >> H,Equation 4 approaches a plug flow condition, qs = V0H. By normalizing the flux in Equation 4 as q* = qs/(V0H), such that q* = D* (1 – e–1/D*), we capture the sensitivity of flux to soil thickness with the value D* = dc/H. As D* increases, variations in soil thickness increasingly impact the flux at a given slope.

We incorporate Equation 4 into simulations of landscape evolution by allowing V0 to vary as a function of slope and a constant, k [L t–1]: 
for the case of thickness-dependent linear diffusion, or 
for thickness-dependent nonlinear diffusive transport (Roering, 2008). Here, Sc [] is a threshold slope, with values commonly of ∼1 (Roering et al., 1999).

Field Location and Data

We test this theory with observations from an area of the Gabilan Mesa contained within the Pliocene shallow-marine Pancho Rico Formation (Durham and Addicot, 1965). The Pancho Rico Formation comprises beds of mudstone with scattered very fine sand; fine sandstone with scattered coarse sand; and fossiliferous, pebbly conglomerate. Stratigraphy is exposed only at the heads of recently incised arroyos, where two stratigraphic columns were measured (Fig. 1; see the GSA Data Repository1 for stratigraphic data).

Soil thickness was measured in soil pits and by driving a rod down to the soil-bedrock contact. We averaged thicknesses measured within one pixel in a lidar digital elevation model (1 m), which amounts to 115 measurements from soil pits and 112 collected with the rod (Fig. 1). The soil-bedrock contact was reliably identified with the rod in some soils, but the rod occasionally penetrated easily disaggregated bedrock. Thus, we base our interpretations and modeling of data only on measurements collected in soil pits. Vertically measured thicknesses were converted to surface-normal measurements using local slope.

To compare the above theory with field data, we assume topographic steady state. Given this, hillslopes should be adjusted to transport the flux received from upslope (Gilbert, 1909), which should scale linearly with the upslope area on the landscape. This allows us to compare data from different hillslope positions according to the drainage area per contour length, a/b. We calculate a using the Dinf flow routing algorithm (Tarboton, 1997; Perron, 2010), and use the 1 m pixel dimension as b. Slope magnitude is calculated using a second-order, finite-difference kernel in the cardinal directions.

We focus on nonlinear relations between slope and flux, as hillslopes tend to become more planar downslope, with slopes below ∼0.5 (Figs. 1 and 2A). A new GTL prediction of flux is presented in each of the three rows of plots shown in Figure 2. These relationships show a nonlinear function of slope (Fig. 2A; Roering et al., 1999); the product of soil thickness and a nonlinear function of slope (representing transport due to plug flow; Fig. 2C; Heimsath et al., 2005); and the result of combining Equations 4 and 6 (Fig. 2E; e.g., Roering, 2008). Plotted against a/b are the components of each GTL that should vary smoothly as a function of a/b. Modeled values for these parameters, shown as dashed lines, were computed by least-squares minimization between observed data and a prediction based on the relevant GTL and the steady-state assumption (e.g., Roering et al., 1999; justification is provided in the Data Repository). For GTLs that depend on soil thickness, we use the mean soil thickness from our observations to develop this single prediction. In the right column of plots, we test if deviations from the modeled predictions correlate with soil thickness.

Simulating Landscapes

To illustrate the variable role lithology may play in shaping landscapes as a function of D*, we numerically integrate Equations 1 and 2 through the application of Equations 4 and 5 from a flat initial topography in one dimension. The model domain was subjected to uniform uplift except at the boundaries, whose elevations were held fixed. To approximate a situation such as that in the Gabilan Mesa, we assume that lithology controls W0 and let W0 vary with the amount of exhumation that has occurred at a point. Here, modeled topography is incising through laterally continuous oscillations between 2.5-m-thick “recessive” and “resistant” beds; we specify these as having values of W0 50% greater or less than a mean value, respectively. The cyclic values of W0 are advected upward by U, causing the model to evolve to a state in which relief oscillates around a mean value. To ensure that initial topographic conditions were eliminated, we ran the model for 4 × 106 model years, a total uplift of >10× the final basin relief.


Soil Thickness and Topography

For a given a/b, slopes are higher above thinner soils (Figs. 2A and 2B). This is seen as a negative bias in the residuals as a function of soil thickness (Fig. 2B). Bias increases for the case that a soil thickness–slope product relates soil thickness to transport rate at a given slope (Figs. 2C and 2D). The velocity profile–based prediction does not eliminate the noise present in the data (Fig. 2E). However, with a dc of 12 cm, the segregation of points at a particular value of a/b based on thickness (e.g., Figs. 2A and 2B) is no longer present. This is reflected in the residuals, which are generally centered near zero for all but the thinnest soils (Fig. 2F).


In Figure 3A, we plot our dimensionless flux, q*, as a function of D* to reference each of the simulation results shown. Variations in H arise from the variations in W0 in these simulations. As a result, there are steady-state values of D* for resistant (high D*) and recessive (low D*) layers; bars in Figure 3A show the range of these values within each simulation. We scale k in Equation 5 so that flux remains relatively constant for a given slope and soil thickness, which reduces fluctuations in relief as dc varies.

Soil thickness varies in all the simulations, but small values of D* (Fig. 3Ai) decouple soil production and transport in a way that prevents variations in soil thickness (and hence lithology) from being expressed in the modeled topography. For larger values of D*, balancing upslope sediment supply requires that V0 (and thus slope) must change to compensate the changes in H that arise at lithologic boundaries. D* values of ∼1 (Fig. 3Aii) produce weak correspondence between slope and soil thickness—an effect that becomes progressively more pronounced as D* increases (Figs. 3Aiii and 3Aiv).


Beds in the Pancho Rico stratigraphy can be traced from exposures in arroyos to laterally continuous undulations in slope on soil mantled hillslopes (Fig. 1). Soils tend to be thinner in steep portions of these shadow beds than in low-slope sections (Figs. 2A and 2B). A linear relationship between soil thickness and transport overestimates transport through thick soils, suggesting that transport does not occur at equivalent rates throughout the soil column (Figs. 2C and 2D). A prediction of flux based on Equation 3 accounts for variations in transport with depth in the soil column, and produces a balanced estimate of sediment flux along the length of the hillslope when dc is 12 cm (Fig. 2E). This prediction results in some bias in the residuals of the thinnest soils toward negative values (Fig. 2F), perhaps related to an underestimation of high slopes on gridded data or the fact that the proposed velocity profile is an oversimplification (e.g., Kirkby, 1967; Lewis, 1976). What remains to be demonstrated is that W0 varies with the lithologies present and that the proposed velocity profile (Equation 3) matches velocity profiles from the Gabilan Mesa.

Numerical modeling (Figs. 3Aii–3Aiv) demonstrates that landscape forms influenced by subsurface lithology can occur if rates of soil transport, qs, depend on soil thickness and slope (Equations 4–6) and if different lithologies require different soil thicknesses to achieve equivalent soil production rates (e.g., Ahnert, 1976; Heimsath et al., 2005). Shadow beds form in simulations of hillslopes crossing lithologic contacts, represented with values of W0, when soil is transported with the proposed velocity profile as long as H is small relative to dc (Figs. 3Aiii and 3Aiv).

Soil transport is commonly related to the disturbance of soils yielding a net downslope flux (Davis, 1892; Culling, 1960; Roering et al., 1999). The velocity profile we model could alternatively reflect a decrease in the frequency of soil-disturbing events with depth. We use the term “activity” to refer to some measure of the soil disturbance accomplished by a range of processes. Ground-squirrel burrows and drying cracks in soils suggest that burrowing and shrink-swell processes are two important agents in soil disturbance in the Gabilan Mesa. Data on the depth dependence of these and other soil disturbing processes suggest that their activity commonly declines with depth (Fig. 3B). This suggests that the upper portions of the soil consist of rapidly mobilized material, while deeper portions may be disturbed less frequently. In thin soils, the high-activity zone may include the entire mobile layer (e.g., high D*), generating a tight coupling between soil production and landscape form. In thick soils, transport may be decoupled from deep portions of the soil column, and landscape forms may be less sensitive to soil production processes.

How common is lithologic control on hillslope form? Our simulations of one-dimensional landscapes evolving by linear, depth-dependent transport indicate that the sensitivity of soil-mantled landscapes to thickness occurs when q* > ∼0.5–0.7 (Fig. 3A), which occurs when D* > ∼1–2. Given a hillslope where soil thickness has adjusted so that soil production rates match rock uplift rates, D* can be calculated as: 

In the case that α–1 and dc are similar and θ is small, the equilibrium D* only depends on U and W0. In this situation, values of D* will exceed ∼1, and fluxes will be sensitive to soil thickness, when the ratio of U to W0 is greater than ∼e–1 (∼0.4). For a value of W0 around ∼0.1 mm yr–1 (Dietrich et al., 2003), and assuming W0 and U are uncorrelated, soil mantled landscapes where uplift (or mean erosion rates) are greater than 0.04 mm yr–1 are likely to be sensitive to variations in soil thickness (and therefore lithology). This suggests that a broad range of soil-mantled hillslopes may encode some signature of the underlying rock type. However, the dominance of nonlinear hillslope transport (Roering et al., 1999) in steep terrains may allow small changes in slope to compensate for potentially large variations in W0. Ultimately, nonlinear hillslope transport may limit features such as shadow beds to landscapes with large variations in W0 and in which base-level lowering is rapid enough to produce thin soils, but slow enough (and with hillslopes short enough) to maintain slopes well below their critical values.

We thank K.D. Chadwick, S. Moon, A. Sanquini, G. Seixas, and E. Shelef for invaluable discussions and feedback on earlier versions of this manuscript, and R. Klier, J. Perkins, and G. Seixas for their efforts in the field. Reviews and comments by A. Heimsath, J. Roering, an anonymous reviewer, and J. Spotila greatly improved this manuscript. The generosity of L. Homen allowed this work to proceed. Support for Johnstone was provided by National Science Foundation grant EAR-TECT-105581. Fieldwork was funded through Stanford McGee/Levorsen and Geological Society of America Graduate Research grants to Johnstone. Acquisition and processing of airborne laser swath mapping lidar data were provided by the National Center for Airborne Laser Mapping seed grant program.

1GSA Data Repository item 2015039, additional information, discussion, and data, is available online at www.geosociety.org/pubs/ft2015.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.