## Abstract

The paleoslope estimation method uses a threshold-shear-stress criterion, together with field-based measurements of median grain size and channel depth in alluvial gravel deposits, to calculate the threshold paleoslopes of alluvial sedimentary basins. Threshold paleoslopes are the minimum slopes that would have been necessary to transport sediment in those basins. In some applications of this method, inferred threshold paleoslopes are sufficiently steeper than modern slopes that large-magnitude tectonic tilting must have occurred in order for sediments to have been transported to their present locations. In this paper, we argue that autogenic cycles of channelized fluvial and sheet flow in alluvial sedimentary basins result in spatial and temporal variations in the threshold slope of gravel transport that can, under certain conditions, cause gravel to prograde out to distances much longer than previously thought possible based on paleoslope estimation theory (i.e., several hundred kilometers or more from a source region). We test this hypothesis using numerical models for two types of sedimentary basins: (1) an isolated sedimentary basin with a prescribed source of sediment from upstream, and (2) a basin dynamically coupled to a postorogenic mountain belt. In the models, threshold slopes for entrainment are varied stochastically through time with an amplitude equal to that inferred from an analysis of channel geometry data from modern rivers. Our models show that when local threshold slope values vary stochastically and sediment supply is relatively low compared to transport capacity, alluvial gravels can persistently prograde at slopes far below the threshold slopes predicted by paleoslope estimation theory. The result holds whether the stochastic changes in threshold slope are autocorrelated along the entire channel profile or occur in localized sections of the channel profile. As such, our model results suggest that long-runout gravels do not require steep regional slopes in order for transport to occur. We conclude that the minimum progradational slopes of fluvial sedimentary basins adjacent to postorogenic mountain ranges are functions of both the mean and coefficient of variation of water-flow depths as well as the density and texture of the bed sediment.

## INTRODUCTION

*S*

_{c}is the threshold channel slope,

*D*

_{50}is the median grain size of the mobile bed sediments,

*H*is the flow depth required to entrain the gravel, and ψ (calculated to be 0.094 when the input grain size is the

*D*

_{50}) is a global constant for all gravel-bed channels with noncohesive banks. Equation 1 assumes that flow is quasi-steady, bed forms are not present, and channel banks are noncohesive. Gravel-bed channels with noncohesive banks modify their cross-sectional geometries (e.g., channel width-to-depth ratio) to produce a bank shear stress that is nearly equal to the threshold of entrainment (Parker, 1978; Andrews, 1984); hence, channel width can be eliminated from the analysis that leads to Equation 1. When this relationship was tested against a data set of modern gravel-bed rivers compiled by Church and Rood (1983), scatter was observed in the ratio between threshold slopes calculated with Equation 1 and observed channel slopes (i.e.,

*S*

_{c}/

*S*

_{obs}). To within one standard deviation of the geometric mean, predicted threshold slopes for channels with noncohesive banks were up to a factor of 1.72 both greater and lesser than the observed slopes (Paola and Mohrig, 1996).

Equation 1 has recently been used to constrain the paleotopography of gravel-dominated alluvial deposits in the North American Cordillera (McMillan et al., 2002; Heller et al., 2003). These deposits are considered to be long-runout gravels because they extend several hundred kilometers from the bedrock drainage basins from which they originated. Heller et al. (2003) calculated paleoslopes based on measured median grain sizes and proxies for mean channel depth within outcrops of the Shinarump Conglomerate of the Upper Triassic Chinle Formation, Lower Cretaceous conglomerate units (e.g., Buckhorn Member of the Cedar Mountain Formation in central Utah and the Cloverly Formation in Wyoming), and the Tertiary Ogallala Group. These long-runout gravels were deposited on top of finer-grained and previously flat-lying sedimentary units. The underlying units were interpreted to be flat-lying during deposition because the depositional facies within the stratigraphic units are indicative of low transport energy. In the absence of tectonically driven tilting, gravel deposition must have initially occurred in the proximal portion of the basin until the upper surface of the clastic wedge aggraded to threshold slope, and gravel could then be transported farther into the basin. The threshold slopes calculated from Equation 1 were significantly greater than the measured slopes of the upper surface of the clastic wedge (i.e., the grade of the fluvial system). As such, Heller et al. (2003) inferred that pre- or syndepositional tectonic tilting must have taken place in order to produce the slopes necessary to entrain sediment of the textures observed in these deposits.

A key assumption of the paleoslope estimation method, i.e., using Equation 1 to calculate threshold paleoslopes, is that the flow depth inferred from the thickness of the fining-upward sequence in outcrop is equivalent to the minimum effective flow depth during gravel transport. Several authors have pointed out that the fining-upward sequence thickness can either underestimate (i.e., if they are partially truncated by erosion) or overestimate the reach-averaged flow depth during gravel transport (Paola and Mohrig, 1996). Although this is a source of error for the paleoslope estimation method, this assumption allows mean data measured from outcrops to be used in place of flow depth *H* in Equation 1. However, flow depth varies in space and time in ways that could affect the accuracy of the paleoslope estimation method.

For example, autogenic cycles of channelized fluvial and sheet flow have long been observed in laboratory flume studies of fans, including those driven by constant water and sediment supply (e.g., Schumm, 1977; Whipple et al., 1998; van Dijk et al., 2008) (Fig. 1). In the first published description of this behavior, Schumm (1977) described the elements of an autogenic cycle as follows: “First, water issued from the canyon and spread in a sheet at the fan head…. This expected runoff configuration was interrupted periodically by fan-head trenching… Following development of a fan-head trench, the source area channel was rejuvenated, and increased sediment delivery aided the backfilling of the trench and restoration of the smooth fan-head configuration.” Whipple et al. (1998) showed that autogenic cycles of channelized fluvial and sheet flow are associated with order-of-magnitude changes in the effective width of the flow (e.g., see fig. 10 *in*Whipple et al., 1998). The continuity equation for discharge, together with Equation 1, implies that order-of-magnitude variations in the effective width of flow result in comparably large changes in flow depth and threshold slopes of entrainment. As such, autogenic cycles on fans can result in significant temporal variations in the threshold slopes of transport.

Autogenic cycles of channelized fluvial and sheet flow are difficult to study in the field given the nearly ubiquitous imprint of late Quaternary climatic changes on landscapes (e.g., van Dijk et al., 2008). Whether they are argued to be autogenic (Schumm and Hadley, 1957; Patton and Schumm, 1975; Bull, 1997; Pelletier and DeLong, 2004) or allogenic (Knox, 1983; Bull, 1991) in nature, however, cycles of cutting and filling in alluvial channels that result in large changes in the effective width of flow occur in many actively aggrading fluvial systems (i.e., those dominated by distributary channel networks). The cutting/entrenchment phase of such cycles is often accompanied by local channel narrowing as unit stream power, incision rate, and rate of channel narrowing operate in a strong positive feedback. Conversely, channels in the depositional phase of this cycle commonly widen as they aggrade, both in nature (e.g., Bull, 1997; Pelletier and DeLong, 2004) and in experiments (e.g., Sheets et al., 2002). Natural channels can oscillate between these two modes through time, raising the possibility that they might transport the coarsest load only during episodes of active channel incision when the channel is far from an equilibrium or self-forming geometry. In the southwestern United States, the typical time scale for the cutting and filling of large, valley-floor channels is ∼1000 yr (Waters and Haynes, 2001; Mann and Meltzer, 2007). Although the cyclic nature of channelized fluvial and sheet flow in piedmont alluvial channels and channel fans has been most well studied in the southwestern United States, such cycles are not unique to arid fluvial systems, as the flume experiments described previously demonstrate. Such models have neither climate nor weather, yet they nevertheless exhibit substantial deviations from equilibrium or self-forming channel geometries. Meandering channels also exhibit spatial and temporal oscillations in which lateral meandering and deposition occur preferentially in relatively low-gradient sections of the valley floor, which then builds a steep alluvial toe slope that facilitates entrenchment and narrowing downstream (Jones and Harper, 1998). Given these considerations, it is reasonable to expect that threshold slopes for entrainment in natural rivers vary stochastically in time. Episodes of periodic entrenchment triggered by climatic changes or autogenic dynamics could result in time periods of lower-than-average threshold slopes. Moreover, given that the narrowing phase of autogenic cycles is associated with channel entrenchment and the reworking of sediments downstream, and that the widening phase is associated with aggradation (see schematic diagram in Fig. 1B), it is possible that the stratigraphic record preferentially preserves channels and channel fans with unusually wide and shallow flows.

If the threshold slope of entrainment varies stochastically through time, then the transport of gravel may be controlled not only by a mean threshold condition (governed by mean grain diameter and channel depth) but also the variability through time. In this paper, we test the hypothesis that gravel may be able to prograde episodically at channel slopes substantially lower than the minimum value predicted by paleoslope estimation theory, albeit for geologically brief periods of time when the threshold slope for transport is anomalously low due to channel entrenchment. In this paper, we test the hypothesis that gravel transport is a function of both the mean and coefficient of variation in threshold slope within a numerical modeling framework. Specifically, we pose the following question: Given a sufficient period of time, can gravels prograde several hundred kilometers from their source regions at slopes significantly lower than the threshold slope predicted by Equation 1 given realistic variations in the threshold slope of entrainment in space and/or time?

## NUMERICAL MODEL DESCRIPTION

### Concept and Scope of Numerical Modeling

We used two types of numerical models in order to test our hypothesis that gravels can be transported over long distances at regional slopes lower than the threshold slope predicted by the paleoslope estimation method. The first model, herein called the isolated sedimentary basin model, represents a gravel wedge prograding along a horizontal plane in response to a sediment source at a fixed position at the upstream end of the basin. Deposition and sediment reworking along a channel profile are simulated using a numerical model that includes mass conservation combined with a slope-dependent transport relationship that contains a grain-size–dependent threshold value. This first model allows us to assess the effects of climate change and autogenetic processes alone on gravel progradation, because sediment supply is held constant, and tectonic movements are not considered. The second model, herein called the dynamically coupled postorogenic foreland basin model, represents a gravel wedge prograding within a foreland basin in response to erosional unloading of an adjacent postorogenic mountain belt. This model allows us to observe the effects of feedbacks between erosion in the mountain belt and the evolution of the adjacent sedimentary basin.

The motivation for applying two types of models is that we want to isolate the effects of stochastic variation in threshold slopes on gravel transport in our isolated sedimentary basin model independent of tectonics and then compare the results to those of a model that includes tectonics by coupling a basin and its adjacent mountain belt. The dynamically coupled postorogenic foreland basin model allows us to determine whether or not feedbacks between sediment supply to the basin and the base level of erosion for the mountain belt (controlled by transport efficiency of sediments across the basin) act in concert to control the rates and slopes at which a gravel wedge can prograde within a basin. In a numerical modeling study, Paola et al. (1992) showed that varying sediment supply and accommodation rates sinusoidally through time had significant effects on the regional slope of an aggrading alluvial fan. The three long-runout gravels described by Heller et al. (2003) were each associated with a particular period of tectonic activity during which sediment supply and accommodation were not constant. For this reason, we include a dynamically coupled model of a mountain belt interacting with an adjacent foreland basin in addition to a model of an isolated basin with a prescribed sediment supply.

*C*

_{v}is the coefficient of variation, σ

_{Sc}is the standard deviation of threshold slope, and is the mean of threshold slope. As the

*C*

_{v}of the threshold slope distribution increases, the recurrence interval between periods of threshold slopes significantly lower than average decreases. To estimate an appropriate value for

*C*

_{v}, the mean and standard deviation of the natural logarithm of the observed threshold slopes, i.e., μ and σ, were calculated from the Church and Rood (1983) data set and related to σ

_{Sc}and using the following expressions (Turcotte, 1997):

*C*

_{v}is more appropriate than the

*C*

_{v}for channel slope as an estimate of the variability of threshold slopes of individual gravel-bed rivers through time because mean flow depth and median grain size are far more variable between rivers (as in the Church and Rood [1983] data set) than within an individual river to which the paleoslope estimation method might be applied. A threshold slope of 0.0014 is predicted by Equation 1 for a gravel-bed river with an average flow depth of 1.68 m and a median bed-material grain size of 25 mm. We use this threshold slope as our representative value in the simulations described later herein because it is within the range of threshold slopes calculated for the Ogallala Group, a type-example long-runout gravel deposit. When a

*C*

_{v}value of 1.66 and a mean threshold slope of 0.0014 are used to generate a lognormal distribution, ∼70% of threshold slopes are less than the mean (Fig. 2B).

### Sediment Transport Component of the Numerical Model

*q*

_{s}is the volumetric unit sediment flux (m

^{2}/s),

*k*

_{g}is the transport coefficient (m

^{2}/yr),

*S*is the local channel slope (unitless),

*h*is the channel-bed elevation (m),

*t*is time (yr), and

*x*is distance along the channel from the drainage basin headwaters (m). The transport coefficient

*k*

_{g}can be estimated from an empirical relationship that depends on discharge and river type (Paola et al., 1992). Paola et al. (1992) estimated transport coefficients of ∼1.0 × 10

^{4}and 7.0 × 10

^{4}(m

^{2}/yr) for braided and meandering rivers, respectively, draining catchments of length equal to 100 km and a mean annual precipitation of 1 m. We chose a value of 1.0 × 10

^{4}(m

^{2}/yr) for our numerical experiments, which corresponds to braided, gravel-bed streams with relatively small drainage areas.

Threshold slope (*S*_{c}) values were sampled from a lognormal distribution with a mean of 1.4 × 10^{−3}. The coefficient of variation was varied between experiments but nominally had a value of 1.66 (i.e., the value calculated for naturally occurring rivers from the Church and Rood [1983] data set). Once the threshold slope was determined for a given point along the channel longitudinal profile, that value was held constant for the prescribed sampling period, which we defined as the length of time that the channel was characterized by a particular threshold slope. The sampling period value may have a significant effect on both types of models considered in this paper because it determines whether the simulated alluvial channels remain in a state of low threshold slope for long periods of time or whether channels rapidly jump back and forth between states of enhanced or limited gravel entrainment.

The time scales of autogenetic and allogenic cut-and-fill cycles documented in natural river systems (e.g., Waters and Haynes, 2001; Mann and Meltzer, 2007; Harvey et al., 2011) can be used as a guide for choosing sampling periods for both models. Autogenic processes that lead to changes in alluvial channel position through time (i.e., such as channel meandering and avulsion) have received much attention in the literature (Miall, 1996). However, very few studies have focused on autogenic processes that lead to alluvial channel incision. In the southwestern United States, for example, the recurrence interval of late Holocene arroyo cut-and-fill cycles for the San Pedro River, which drains a region of 1000 km^{2}, has been well constrained (Waters and Haynes, 2001). Based on the radiocarbon dates at the base of each scour, arroyos were incised and filled within 500–4000 yr, with a predominant period of ∼1000 yr. Similar trends were observed for Holocene cut-and-fill cycles for the Wildhorse Arroyo and Archuleta Creek of the Upper Cimarron River drainage basin in New Mexico (Mann and Meltzer, 2007). Radiocarbon and optically stimulated luminescence (OSL) dating of arroyo cycles in Buckskin Wash, Utah, revealed that at least four cut-and-fill cycles have occurred since 3 ka, and thus, the predominant period was between 500–1000 yr (Harvey et al., 2011). Although the periodicity of cut-and-fill cycles in the southwestern United States would be different than that of channels located in other climates, we infer that 1000 yr is a reasonable value for the threshold slope sampling period. As such, we chose to hold threshold slopes constant over a period of 1000 yr, which is equal to the time step for both the isolated sedimentary basin and dynamically coupled postorogenic foreland basin models.

Another important issue involving sediment transport within our models is whether the threshold slope should be uniform or varied along the channel profile during each sampling period. Observations of alluvial-fan and valley-floor channels reveal oscillations or instabilities in channel width and slope that appear to have wavelengths on the order of 100 m to 10 km (Pelletier and DeLong, 2004). Therefore, it is possible to have discrete sections of the river profile with lower threshold slopes compared to the majority of the river profile. If this is the case, then gravel transport should primarily take place within these localized regions of increased slope and low width-to-depth ratios. This would suggest that models aimed at understanding the role of autogenetic oscillations in river systems should use spatially non-autocorrelated threshold slopes sampled at a relatively small spatial wavelength. However, channels also respond to climatic and autogenic forcing factors that trigger incision/aggradation along the entire river. Such cases would suggest that models aimed at understanding the role of allogenetic events in river systems should use threshold slope values that vary in time but are highly spatially autocorrelated along the channel longitudinal profile. An example of such an allogenetic event would be entrenchment along an entire river profile due to a major decrease in the ratio of sediment to water from upstream. In this paper, we consider both spatially autocorrelated and spatially non-autocorrelated variations in threshold slope values. When threshold slope values were spatially non-autocorrelated in the numerical models, the threshold slope was sampled from a lognormal distribution and held constant over the sampling period for each alluvial channel pixel representing 1 km. When changes in threshold slope were spatially autocorrelated, a single threshold slope was sampled and applied along the entire alluvial channel profile for the given sampling period.

### Further Description of the Dynamically Coupled Postorogenic Foreland Basin Model

Our dynamically coupled postorogenic foreland basin model is a 2-D model that solves for the longitudinal profile of a channel evolving due to bedrock erosion and flexural-isostatic response to erosional unloading of a mountain belt in addition to sediment transport in an adjacent foreland basin subject to stochastic changes in threshold slopes of entrainment. A 2-D model is sufficient for the purposes of this study because at least some examples of long-runout gravels (e.g., the Ogallala Group conglomerates) are effectively 2-D systems, i.e., gravels were primarily transported perpendicular to a long, linear mountain belt. While the model is 2-D, it should be emphasized that variations in threshold slope implicitly include variations in channel geometry in the third dimension, i.e., channel width. The mountain belt topography at the beginning of each simulation represents the profile for a high-elevation plateau (similar to that of the modern central Andes or Himalaya).

*k*

_{e}is the bedrock erodibility (yr

^{–1}for

*n*= 1.0 and

*m*= 0.5),

*A*is the drainage basin area (m

^{2}), and

*m*and

*n*are constants that determine the dependence of local erosion on discharge and channel slope. Evidence from theory and field studies predict that the ratio of

*m*/

*n*is close to 0.5 (Whipple and Tucker, 1999). Although the slope exponent

*n*can range between 0.66 and 2.0, depending on the relationship between slope and stream power (shear stress), we chose to make the stream power linearly proportional to the slope (

*n*= 1.0 and

*m*= 0.5), consistent with the assumption of many other studies (e.g., Kirby and Whipple, 2001; Snyder et al., 2000). The range in bedrock erodibility (

*k*

_{e}) can be up to five orders of magnitude (10

^{−2}to 10

^{−7}), depending on lithology and the values of

*m*and

*n*(Stock and Montgomery, 1999). We chose a value of 10

^{−6}yr

^{−1}in order to scale the erosion rate to reasonable values (1 mm/yr or less) during the early portion of the simulations, when the channel slopes in the bedrock portion of the model are highest. Alluvial deposition occurs within the bedrock drainage basin when erosion exceeds transport capacity, and, therefore, bedrock that is temporarily buried by alluvium does not erode until the overlying sediment is removed. In this way, the transition between the bedrock and alluvium is free to migrate laterally in response to the coupled evolution of the mountain belt–foreland basin system in the model.

*w*is the deflection of the crust (m),

*D*is the flexural rigidity (Nm), ρ

_{m}is the density of the mantle (kg/m

^{3}), ρ

_{s}is the density of the mountain crust or foreland sediment (kg/m

^{3}),

*g*is gravity (m/s

^{2}), and

*L*(

*x*) is the topographic load (kg/ms

^{2}). A Fourier transform method was used to solve this PDE (partial differential equation) at each time step in the model following Watts (2001). Appropriate flexural rigidities for a high-elevation plateau can range between 2.4 × 10

^{23}and 3.0 × 10

^{24}Nm, when Poisson’s ratio is 0.25 and the Young’s modulus is on the order of 10

^{11}Pa (Stewart and Watts, 1997; Jordan and Watts, 2005). We chose to apply a flexural parameter of 150 km. The flexural parameter is defined as the following (Turcotte and Schubert, 1992):

^{23}Nm, i.e., within the range of values for mountain belts with high plateaus. In the dynamically coupled postorogenic foreland basin model, Equation 8 was applied twice, first to determine the initial flexural profile beneath the topographic load of the initial profile of the mountain belt at time zero, and second to solve for rock uplift in response to erosional unloading of the mountain belt at each time step. The flexural-isostatic component of the model was validated using comparison with analytic solutions for a line load. A mountain belt with an average height of 2.8 km and a width of 550 km triggers a forebulge that is 300–400 m high, a foredeep that is 290 km wide, and maximum depth of 3 km when the rigidity is on the order of 6.0 × 10

^{23}Nm. An issue that arose during the study is: Should the forebulge be a significant topographic barrier between the foredeep and backbulge basins? A topographically significant forebulge would prevent gravel progradation out of the foredeep until the alluvial systems could aggrade to the forebulge spill point. Along the central Andes, however, there is no topographic expression of the forebulge despite the low elevations observed 200–300 km from the thrust front in Bolivia and northern Argentina (i.e., 150–500 m above sea level [a.s.l.]; Horton and DeCelles, 1997). Conversely, the predicted forebulge region based on 2-D flexural modeling in the Himalayan foreland basin is elevated ∼300–400 m above the Ganges plain (Bilham et al., 2003). Interpreting this relief as exclusively the result of forebulge uplift is complicated by the presence of preexisting topography in the Central Indian Plateau. Although the Himalayan foreland may be an exception, topographically expressed forebulges are not commonly observed within modern terrestrial foreland basins, most likely due to erosion and/or dynamic subsidence (DeCelles and Giles, 1996; Catuneanu et al., 2000). In order for the model to predict a forebulge amplitude that is consistent with modern analogs (i.e., a maximum of 150–300 m a.s.l.), we beveled the forebulge crest at the beginning of the model down to the elevation of the initial fluvial profile. However, this initial topographic adjustment did not modify the flexural response (i.e., rock uplift and/or subsidence that occurs due to flexure of the lithosphere as a result of erosional unloading of the mountain belt) of the forebulge during the simulation, which is an important control for sediment accommodation and depositional basin regional slope, to postorogenic denudation of the adjacent mountain belt.

### Numerical Model Implementation

Parameters for both the isolated sedimentary basin and dynamically coupled postorogenic foreland basin models are reported in Table 1. The isolated sedimentary basin model simulated 10 m.y. of gravel progradation using a time step of 1 k.y. and a pixel spacing of 1 km. The mean threshold slope was held constant at a value of 0.0014, and the threshold slope sampled from the lognormal distribution was applied to the entire channel profile. Simulations with a range of *C*_{v} values and a range of ratios of the sediment supply to the transport coefficient, *q*_{s}/*k*_{g}, were conducted to test the sensitivity of gravel progradation to variations in these parameters. Sediment flux into the basin was held constant at a rate of 6 m^{2}/yr for the simulations with variable *C*_{v}, consistent with a *q*_{s}/*k*_{g} value of ∼0.00075. We chose this value for sediment supply because it is equivalent to the average sediment supply coming out of the drainage basin during the period of highest sediment flux for each of the dynamically coupled postorogenic foreland basin model simulations discussed in the Results section. Simulations for each *C*_{v} value were repeated 20 times, and the reported results are an average of those model runs. Following the simulations with variable *C*_{v} and *q*_{s}/*k*_{g}, we also tested the sensitivity of gravel progradation rates and slopes to the presence or absence of spatial autocorrelation in threshold slope along the channel profile.

The dynamically coupled postorogenic foreland basin model simulated 20 m.y. of gravel progradation using a time step of 1 k.y. and pixel size of 10 km. The initial topographic profile for the bedrock drainage basin was an actively uplifting mountain belt with frontal slopes of 0.011 and a low-relief plateau interior. In the first simulation of this model type, threshold slopes were sampled from a lognormal distribution with a *C*_{v} of 1.66 and a mean value of 0.0014, i.e., equivalent to that of the Church and Rood (1983) data set and, therefore, representative of natural variations of slopes of channels with similar mean flow depths and grain sizes. Changes in threshold slope were spatially autocorrelated along the channel profile. Simulations with spatially non-autocorrelated changes in threshold slope were not considered here due to the results obtained when comparing the isolated sedimentary basin model simulations with uniform or nonuniform changes in threshold slopes, which we address in the Discussion section. Topographic surfaces were sampled for display at 0.0, 0.25, 1.0, and 20.0 m.y. of simulation. We also report on additional simulations that were run to test the effect of varying *C*_{v} and the sampling period on the evolution of postorogenic gravel progradation.

## NUMERICAL MODELING RESULTS

### Results for the Isolated Sedimentary Basin Model

#### Effect of C_{v} on Gravel Progradation

Time-series results for regional slope and the position of the gravel front generated by the isolated sedimentary basin model with spatially autocorrelated changes in threshold slope are shown in Figure 3. Regional slopes (defined herein as the ratio of the relief of the depositional basin to the length of the basin) decrease systematically with increasing *C*_{v}. When the *C*_{v} value is greater than or equal to 1.0, the gravel wedge progrades at regional slopes that are less than the mean threshold slope (i.e., 0.0014). Thus, the model results suggest that natural rivers with a *C*_{v} of 1.66 and the channel conditions considered here can prograde at regional slopes significantly lower (i.e., a factor of 2 in this case) than that predicted by the paleoslope estimation method.

The location of the gravel front increases as a power-law function of time with an exponent that increases from ∼0.5 to 0.8 as the *C*_{v} increases from 1 to 8 (Fig. 3B). Even prior to 1.0 m.y., clastic-wedge progradation in the model is rapid, and gravels begin to reach downstream distances that are considered “long-runout,” i.e., greater than tens of kilometers (Heller et al., 2003), when the *C*_{v} is high. Beyond 1 m.y., progradation rates decrease to long-term average rates of between 30 and 55 km/m.y., with the precise value dependent on *C*_{v}. Increasing the *C*_{v} value by a factor of 4 increases the total progradation by less than a factor of 2 after 10 m.y. In such cases, gravel-bed channels with a *C*_{v} value comparable to natural rivers, i.e., 1.66, are able to transport gravel up to ∼400 km away from their source regions given sufficient geologic time (i.e., 10 m.y.).

#### Effect of q_{s}/k_{g} on Gravel Progradation

Numerical modeling studies have shown that alluvial basin slopes are also affected by sediment supply rates and sediment transport rates (Paola et al., 1992). Therefore, we conducted simulations to determine the variation in the behavior of this model with respect to variations in *q*_{s}/*k*_{g} (Fig. 4). We chose not to test the sensitivity of gravel progradation to *q*_{s} and *k*_{g} individually because the regional slope of a clastic wedge is sensitive to the ratio between these two parameters rather than their individual values. The regional slopes shown in Figure 4 are the long-term regional slopes illustrated in Figure 3. When the *C*_{v} value is equal to 0, the resulting regional slopes are always greater than the mean threshold slope (i.e., 0.0014) calculated using Equation 1. When threshold slopes are allowed to vary stochastically in time, long-term regional slopes decrease below the mean threshold slopes when *q*_{s}/*k*_{g} is less than 0.00125. Values of *q*_{s}/*k*_{g} above ∼0.00125 correspond to cases in which sediment supply to the basin exceeds the ability of the gravel-bed channels to convey that sediment at the mean threshold slope. As the *q*_{s}/*k*_{g} ratio increases above 0.01, the long-term regional slope increases as a power-law of *q*_{s}/*k*_{g}. In such cases, the gravel channels aggrade with little dependence on the value of *C*_{v}. Conversely, values of *q*_{s}/*k*_{g} less than ∼0.00125 correspond to cases in which the transport capacity at the calculated mean threshold slope approaches the sediment supply. The *q*_{s}/*k*_{g} value used in the analysis shown in Figure 3 is located in this region of the plot. In this situation, the *C*_{v} value has a significant effect on the long-term regional slope. The relationship between regional slope and *q*_{s}/*k*_{g} implies that sedimentary basins located adjacent to actively uplifting bedrock drainages may not be able to produce long-runout gravels with low regional slopes early on while sediment supplies are high. However, once active uplifting ends, the regional slopes of the gravel-bed rivers can lower in response to infrequent periods of gravel transport.

#### Effect of Local versus Large-Scale (i.e., Autocorrelated) Variations in Threshold Slopes on Gravel Progradation

Simulations with spatially autocorrelated and non-autocorrelated changes in threshold slope were conducted to test the sensitivity of gravel progradation to the presence/absence of spatial autocorrelations in threshold slope along the channel profile. We found that the sedimentary basin regional slopes were identical in the two cases, i.e., the results for spatially non-autocorrelated threshold slopes lie directly on top of the results obtained with spatially autocorrelated changes illustrated in Figure 4. In the case of spatially autocorrelated threshold slopes, there are time periods of high effective transport rates (when threshold slopes are low across the entire basin). One might expect such cases to have higher effective transport rates compared to cases with spatially non-autocorrelated threshold slopes, but the results of our model do not support this expectation. In the case of spatially non-autocorrelated variations in threshold slope, gravel is transported in small episodic steps that are localized in space but occur commonly through geologic time. In the case of spatially autocorrelated variations in threshold slope, gravel is transported long distances in brief periods of unusually low threshold slope, but then the basin experiences long periods with no transport anywhere until the entire basin can again achieve a transport condition characterized by a low threshold slope (e.g., a phase of active entrenchment). Therefore, the selection of either spatially autocorrelated or non-autocorrelated changes in threshold slope does not significantly affect the results of either the isolated sedimentary basin or dynamically coupled postorogenic foreland basin model, provided that sampling periods are small, i.e., on the order of millennia, compared to the duration of the model.

### Results for the Dynamically Coupled Postorogenic Foreland Basin Model

In addition to climatically or autogenically driven changes in threshold slope, tectonically driven changes in sediment supply and basin accommodation might affect both the regional slope and rates at which gravels are transported. Three types of results are reported for the dynamically coupled postorogenic foreland basin model: (1) a description of the ways in which bedrock drainage relief, rock uplift rates, sediment supply rates, and regional slope evolve during a simulation with a *C*_{v} value that is representative of natural rivers (i.e., 1.66); (2) a comparison of regional slope results for simulations with varying *C*_{v} values; and (3) comparison of regional slope results for simulations with varying sampling periods.

#### Effect of C_{v} and the Coupling of Sediment Supply and Accommodation Space Creation on Gravel Progradation

The first 250 k.y. of the simulation with a *C*_{v} value of 1.66 is marked by knickpoint retreat toward the plateau divide and an increase in basin-average erosion rates (Fig. 5A). Steep slopes and large contributing drainage areas (i.e., a proxy for discharge) lead to peak bedrock incision rates of 0.7 mm/yr near the mountain front. These high erosion rates near the mountain front begin to decrease slopes, and thus erosion rates, over the first 250 k.y. as the main bedrock channel adjusts to slower (i.e., postorogenic) rock uplift rates. Overall, erosional adjustment of the drainage network over this time period does not significantly change maximum erosion rates because the topography after 250 k.y. of simulation still closely resembles the original topography at the end of active deformation.

High basin-averaged erosion rates over the first 250 k.y. following the cessation of active deformation lead to isostatic rebound, gravel input into the depositional basin, and gravel front progradation. Erosional unloading of the mountain belt causes between 40 and 80 m of rock uplift within the bedrock drainage basin over the first 250 k.y. Bedrock incision near the topographic divide is an order of magnitude less than incision at the mountain front in this time period, and, as a result, peak elevations increase by 40 m at the topographic divide. Overall, however, the mean elevation of the mountain belt is less than its original value, as it must be, because only passive (isostatic) uplift is occurring. Isostatic rebound in the adjacent depositional basin increases proximal alluvial slopes because rock uplift exponentially decreases away from the mountain front. High basin-average erosion rates during the first 250 k.y. also lead to an increase in gravel flux into the depositional basin. Gravel flux reaches a maximum value of 10 m^{2}/yr and an average of ∼5–6 m^{2}/yr at this time. As a consequence of increasing sediment supply and proximal channel slopes, the gravel front rapidly progrades ∼150 km (i.e., a time-averaged rate of 0.6 m/yr) from the mountain front in the first 250 k.y.

After 1 m.y. of postorogenic time, basin-averaged erosion rates continue to decrease as the increase in peak elevations and bedrock channel slopes near the topographic divide is accompanied by a decrease in channel slopes near the mountain front (Fig. 5A). In contrast to the initial 250 k.y., the highest erosion rates (0.6 mm/yr) are concentrated within the intermediate portion of the drainage basin instead of at the mountain front. Reduced slopes and temporary sediment storage on the time scale of 100,000 yr both act to limit erosion rates within the lower portion of the drainage system. Overall, the basin-averaged erosion rates at this point in time approach rates expected for postorogenic denudation (Matmon et al., 2003; Reiners et al., 2003). Erosion in the drainage basin leads to ∼100 m and 140 m of surface uplift at the mountain front and topographic divide over the 0.75 m.y. period. Isostatic rebound in the mountain belt is accompanied by subsidence in the forebulge region. As the surface topography lowers in the forebulge region and sediment flux remains high, the distal toe of gravel deposition migrates 100 km farther from the mountain front.

Following 20 m.y. of postorogenic denudation, basin-averaged erosion rates are on the order of 0.01 (mm/yr), as bedrock incision is concentrated near the topographic divide where steep slopes still persist (Fig. 5A). Both peak elevation and erosion rates are comparable to values characteristic for postorogenic mountain belts that have only been eroding for tens of millions of years. Approximately 610 and 3780 m of exhumation occur at the mountain front and topographic divide, respectively, over the 19 m.y. period between 1 m.y. and 20 m.y. following the start of the simulation. Isostatic rebound at the mountain front due to erosional unloading of the mountain belt continues to steepen proximal foreland slopes. The distal toe of gravel deposition migrates an additional 300 km into the foreland basin in this time interval. Gravel progradation persists despite the factor of 3 decrease in gravel flux from the bedrock drainage system compared to the 0.25–1.0 m.y. time interval (Fig. 5B).

Changes in rock uplift and sediment supply rates have an effect on regional slope of the foreland basin. Figure 5B shows the time-dependent behavior of the regional slope. The most rapid change in regional slope (i.e., decreasing from values of 8.0 × 10^{−4} to values of ∼5.0 × 10^{−4}) occurs early in the simulation (prior to 200 k.y.). Initially, increasing gravel flux from the bedrock drainage and nonuniform rock uplift rates near the mountain front lead to rapid gravel migration into the foredeep. This initial period is followed by a regional slope increase that lasts until ∼1.5 m.y., when the toe of the gravel wedge reaches the distal foredeep and mean gravel flux from the mountain front begins to decrease. A topographically significant forebulge would impede gravel progradation beyond the foredeep as the distal toe is forced to aggrade to a forebulge spill point. In this case, however, the forebulge is eroded by the fluvial system to be in grade with the foredeep channel. Beyond 1.5 m.y., the gravel wedge toe is located in the forebulge to backbulge region where rock uplift rates are extremely low, and thus, the decrease in regional slope with time is predominantly controlled by decreasing gravel fluxes at the mountain front and the stochastic variation in threshold slopes.

Regional slopes are nearly linear in log-log space beyond 8–9 m.y. of progradation, and thus represent a power-law relationship between regional slope and time (Fig. 5B). It can also be argued that the decrease in regional slope is nearly linear as early as 2 m.y. following the end of active uplift. When power-law curves are fit to the data in Figure 5B through linear regression of the log-transformed data between 2 and 8 m.y. and between 8 and 20 m.y., the resulting power-law exponents are −0.11 and −0.34, respectively (with correlation coefficients of approximately −0.982 and −0.999, respectively). The differences in exponent values may be related to changes in sediment supply and foreland slopes near the toe of the gravel wedge.

#### Effect of Varying the Value of C_{v} on Gravel Progradation

Stochastic changes in threshold slope are the main mechanism for driving gravel progradation below the mean threshold slope, and thus they have a significant effect on the long-term decrease in the foreland basin regional slope with time. The variation in postorogenic regional slope as a function of the coefficient of variation of threshold slope (i.e., ranging between 0.0 and 3.0) and time is summarized in Figure 6. Heller et al. (2003) proposed that the gravel wedge must achieve and maintain a mean threshold slope to prograde gravel. Figures 5 and 6, however, show a decrease in regional slopes through time that suggests that gravel-dominated foreland basins do not maintain the predicted mean threshold slope (i.e., 0.0014 for our models) as time progresses. When there is no variation in threshold slope, the long-term regional slope approaches the mean threshold slope predicted by Equation 1 for the given parameters. However, the regional slope may not decrease to the mean threshold slope until the sediment input from the drainage basin approaches the transport capacity as demonstrated in Figure 4. When the *C*_{v} is 0.5, ∼5.5 m.y. of postorogenic gravel transport is required to reduce regional slopes below the calculated mean threshold slope value. Regional slopes never increase above the calculated mean threshold slope when the *C*_{v} is greater or equal to 1.0, and thus, the results suggest that the rate of gravel progradation for a given grain size far out into the foreland basin is strongly dependent on *C*_{v} and time, in addition to the value of the transport coefficient. When power-law curves are fit to the data, a negative correlation exists between the *C*_{v} value and the power-law exponent for slope decrease after 7 m.y. The relationship between the *C*_{v} and power-law exponent is best approximated by a logarithmic function (correlation coefficient value of 0.9869). The negative relationship between *C*_{v} and the power-law exponent occurs because higher variability in threshold slopes increases the transport capacity of gravel-bed channels early on, when sediment fluxes are high (i.e., between 0.25 and 1 m.y.), and prevents significant aggradation. As such, regional slopes are already low at the beginning of the final phase of regional slope decrease.

#### Effect of the Sampling Period on Gravel Progradation

Thus far, we have shown how regional slopes respond to changes in the threshold slopes held constant over sampling periods of 1000 yr. However, gravel-bed channels that maintain threshold slopes over longer time scales may evolve differently than channels that vary threshold slopes over shorter time scales. Simulations using the dynamically coupled postorogenic foreland basin model were run for different sampling periods to determine the dependence of gravel progradation on the sampling period duration (Fig. 7). We defined sampling period as the length of time over which the threshold slope for gravel transport was held constant before a new threshold slope was picked from the lognormal distribution. The curves illustrated in Figure 7 represent an average of trials conducted for sampling periods of 500 yr, 1 k.y., 10 k.y., and 100 k.y.

The solid line in Figure 7 represents the sampling period applied to the simulation with a *C*_{v} of 1.66. The decrease in long-term regional slope (i.e., the trend in regional slope beyond 10 m.y.) averaged over many trials remains unchanged as the sampling period is decreased by a factor of two or is increased by one to two orders of magnitude. By 10 m.y., each of the sampling period simulations had experienced enough changes in threshold slope that the regional slope values resulting from the different sampling intervals converged. At shorter time scales, regional slopes for simulations with 10–100 k.y. sampling periods were greater by a factor of two than the slopes for simulations with 0.5–1 k.y. sampling periods. However, this trend may disappear when regional slope is averaged over more trials. Overall, the first-order trend in regional slope, i.e., an increase in regional slope between 0.3 and 1 m.y., followed by a long-term decrease below mean threshold slope, occurred in each simulation. As such, the first-order behavior of the results for either model is not significantly biased by the size of the sampling period.

## DISCUSSION

Results from the simulations presented here suggest that temporal variability in the threshold slope of entrainment has a significant effect on the regional slopes at which gravels prograde in fluvial sedimentary basins, especially in cases of relatively low sediment supply. When gravel progrades with a constant sediment supply in the absence of tectonic influences, the regional slope reaches a constant value, dependent on both the mean and variability (quantified by *C*_{v}) of the threshold slope and the relative sediment supply (quantified by *q*_{s}/*k*_{g}), after only ∼1 m.y. According to our model results, gravel-bed channels that experience greater temporal variation in threshold slope values develop lower regional slopes than those with relatively constant threshold slope values. This behavior is independent of whether the changes in threshold slope are localized (i.e., autogenically forced instabilities in the channel cross section) or are autocorrelated along the entire river profile (i.e., climatically forced entrenchment). If *q*_{s}/*k*_{g} is sufficiently low and/or the *C*_{v} value is sufficiently high, the regional slope can be significantly lower than the value predicted by the paleoslope estimation method. Conversely, when sediment flux is high, channel slopes are controlled by sediment supply, and the “gravel pumping” mechanism of episodic transport at regional slopes below the predictions of the paleoslope estimation method does not apply. In such cases, the channel is overwhelmed with sediment and aggrades to a slope well above the threshold for transport. It is only when sediment supply is low relative to transport capacity that channel reaches can spend significant time periods at or below the threshold for transport. In our model, sediment is stored during these points in space and time, only to be transported downstream when low threshold-slope conditions return.

During postorogenic denudation of a coupled orogen-sedimentary basin system, the regional slope of a gravel-dominated basin decreases in a power-law relationship with time following an initial adjustment period of ∼1 m.y. The time scale of the adjustment period depends on how rapidly the sediment supply begins to decrease following the end of active uplift. A decrease in regional slope initially occurs because channels are adjusted to conditions of higher sediment supply associated with orogenic uplift. However, the long-term regional decrease results from a combination of a decrease in sediment supply as the orogen loses relief and a stochastic variation in threshold slope. The duration and magnitude of the long-term regional slope decrease are dependent on *C*_{v}, *k*_{g}, and the initial relief of the orogen. A power-law relationship with time indicates that the postorogenic foreland basin regional slope can become much lower than the mean threshold slope predicted by paleoslope estimation theory if the *q*_{s}/*k*_{g} ratio is sufficiently small. In the absence of major axial rivers, preexisting topography, and additional tectonic events that affect the postorogenic foreland basin, gravel units will continue to prograde away from the mountain belt through time over a period of more than 20 m.y.

The results for both models suggest that gravels can prograde hundreds of kilometers at low regional slopes over a few million years when sample periods are on the order of millennia. Although 1000 yr is an appropriate (mid-ranged) value for the sampling period, the time scales of cut-and-fill cycles can vary between 500 and 4000 yr. This motivates the question: As channels spend longer periods of time at a given threshold slope (i.e., longer sampling periods), must the long-term progradation of gravels happen more rapidly or at significantly lower slopes than channels that more frequently change their threshold slopes? We infer from our model results that the long-term progradation rates (>10 m.y.) and foreland behavior are approximately time invariant for the tested sampling periods (Fig. 7). As long as the duration of the gravel transport is much larger than the sampling period, then the long-term gravel progradation will occur at similar rates and regional slopes. At shorter time scales, autogenic or climatic processes that lead to changes in flow depth, and thus threshold slopes, every 0.1–1 k.y. are more efficient at prograding gravel than changes in flow depth on time scales of 10–100 k.y., when the *C*_{v} value is representative of modern rivers (i.e., 1.66).

### Evidence for Channel Instability in Ancient Long-Runout Gravel Deposits

If variations in channel width-to-depth ratios over geologically brief time scales (i.e., ∼0.001 to ∼1.0 m.y.) occur during gravel transport, then the internal stratigraphic architecture of the gravel deposits should reflect this behavior. Cycles of entrenchment and fill, which cause periods of low threshold slopes over millennial time scales, have been directly observed and inferred from Pleistocene and Holocene alluvial-fan and valley-floor river outcrops. Although climatic conditions for the Quaternary may be quite different from other periods of time, cut-and-fill cycles are not limited to the Quaternary. For example, DeCelles et al. (1991) interpreted a series of fifth-order lithosomes and fifth-order surfaces (as defined by Miall, 1995) within alluvial-fan deposits from the Paleocene Beartooth Conglomerate in Montana as a cycle of entrenchment and backfilling. Holbrook (2001) observed a continuum of channel-form bounding surfaces within the Cretaceous Muddy Sandstone of southeastern Colorado that represent river scour at different scales. Nested valleys or bundles of channel belts truncated and bound by higher-order surfaces were present throughout the formation. These nested valleys were interpreted to be analogous to modern cut-and-fill deposits that form over time scales of 10^{3} yr. Clear evidence for cycles of entrenchment and periods of low deposition adjacent to the active channel can be found within the internal stratigraphic architecture of the Tertiary Ogallala Group observed in Nebraska. Valley-fill deposits of the Ogallala Group are separated by major erosional surfaces and contain Miocene fossil assemblages of distinctly different ages, which suggest that multiple cycles of cut and fill occurred after the initial incision of the main paleovalleys (Diffendal, 1982). In some outcrops, the younger deposits cut completely through to basal Ogallala and into the underlying mid-Tertiary stratigraphy. Farther south in Texas and New Mexico, cycles of deposition and pedogenesis can be seen within the fine sand to coarse silt lithofacies of the Ogallala Group. In some outcrops, eolian deposits that contain cycles of deposition and pedogenesis cap fluvial sections in the paleovalley fill and paleo-uplands (Gustavson and Holliday, 1999). In the Petrified Forest National Park of northern Arizona, multiple cycles of valley cut and fill have been identified in the lower part of the Chinle Formation (Hasiotis et al., 2001). Lenticular channel deposits of the Shinarump Member cut into contemporaneous overbank deposits that show slight to moderate pedogenesis. Based on the observation of the western U.S. long-runout gravels and other evidence for entrenchment and fill cycles contained within the rock record, it is realistic to invoke stochastic variations in threshold slope to drive gravel transport at low gradients over long distances.

### Conditions Necessary for Producing Long-Runout Gravels

Our modeling results suggest that ancient orogens that have undergone postorogenic denudation (e.g., sufficient relief and appropriate drainage lithologies for generating gravel) should produce long-runout gravels provided that sufficiently resistant lithologies exist in the system for gravel clasts to be preserved without being abraded to sand and provided that sufficient time elapses for the sediment to be episodically reworked down the fan in the episodic or “stepping-stone” manner proposed in this paper. One explanation for the relative rarity of long-runout gravels in the geologic records is the inability of most lithologic types to withstand abrasion over transport distances of hundreds of kilometers (Kodama, 1994; Lewin and Brewer, 2002). Although channel conditions may be suitable for transporting gravel out to the distal foreland basin, gravel may nevertheless be broken down into sand-sized particles within 50–100 km from the source. This hypothesis is consistent with the fact that many long-runout gravels (e.g., the Shinarump and Lower Cretaceous gravels) are dominated by resistive lithologies. For example, ∼24% and 72% of point counts within samples from the Buckhorn Conglomerate were composed of chert nodules and monocrystalline quartz (Currie, 1998). Chert and quartzite are also prevalent within the Lakota, Cloverly, and Shinarump conglomerates (Blakey and Gubitosa, 1984; Zaleha et al., 2001). If time were the only factor in long-runout gravel occurrence, then we would expect to see long-runout gravels in the foreland basin of the central Andes. Toward the eastern edge of the central Andes, for example, Tertiary synorogenic sediments are being exhumed within the Subandean zone, which is the active fold-and-thrust belt today. Within the Tertiary units, the Yecua and Tariquia Formations, which comprise approximately half the thickness of the Tertiary units, are mainly composed of sandstones and mudstones, with some conglomerate units present (Uba et al., 2005). Consequently, studies of the large megafans that drain this region in southern Bolivia and northern Argentina have observed limited gravel to abundant sand and mud along the entire megafan profile (Horton and DeCelles, 2001). Conglomerates should be absent from the rock record if material that is prone to erode into finer grain sizes is being exhumed in the bedrock drainage regions that are the predominant supply of sediment to the depositional basin.

## CONCLUSIONS

The evolution of gravel transport within foreland basins is significantly influenced by the temporal variation in threshold slopes related to autogenic cycles of channelized flow and sheet flow. When the gravel supply is held constant and subsidence is negligible, the regional slope of a prograding gravel wedge reaches a long-term value that is dependent on sediment supply, expressed as *q*_{s}/*k*_{g}, and both the mean and coefficient of variation of the threshold slope of entrainment. When the *C*_{v} value is less than one, the effective regional slope remains equal to or greater than the mean threshold slope, and gravel progradation is minimized. However, when the *C*_{v} value is greater than one and the basin sediment supply is low, gravel can be transported at regional slopes that are well beneath threshold slopes predicted by the paleoslope estimation method. As such, caution should be exercised when evaluating regional slopes measured from the rock record in a tectonic context.

The results of our dynamically coupled postorogenic foreland basin model show that the decrease of foreland basin regional slopes closely follows a power-law relationship with time, although variability in sediment supply and rock uplift cause significant scatter in foreland basin slopes through time. This implies that gravel will continue to prograde in the foreland basin following the end of active uplift for a significant period of time. The power-law regional slope decrease with time observed in the model does not appear to depend on the sampling period. The spatial autocorrelation of threshold slope change also has little effect on the first-order behavior of regional slope decrease, and, thus, either climatically or autogenically driven changes in channel geometry are theoretically capable of transporting gravel at regional slopes lower than that predicted by paleoslope estimation theory. The three described long-runout gravels of the western United States contain lithologies that are highly resistive to grain-size reduction through abrasion, which we infer to also be a necessary condition for the transport of gravel long distances. Both macroscale and internal geometries within the Lower Cretaceous and Ogallala gravels support the possibility of significant variation in channel width-to-depth ratios through time.

We thank ExxonMobil for partial funding of this project. We thank Eric Kirby, Chris Paola, and two anonymous reviewers for critical but constructive reviews that significantly improved the manuscript.