The elevation of continental interiors over time is demonstrably variable. A major part of change in elevation within the continental interior is likely driven by density changes within the upper mantle and by global mantle convection. For example, upper-mantle flow has been invoked as the cause of Neogene uplift of the interior Rocky Mountains and Colorado Plateau, warping and tilting sediment transport slopes that link to the widespread deposition of gravel units within the Great Plains. These geomorphic and sedimentologic features, however, can also be generated by an increase in runoff, since erosion will promote change in elevation due to isostatic compensation and the loading of the lithosphere by the deposition of sediment. To explore the consequences of change in topography and climate, we use a general length-dependent diffusive sediment transport law to model both erosion and deposition that includes the concentrative effects of river systems. The simplicity of the approach means that we can collapse sediment transport to one dimension and couple erosion and deposition with plate flexure. We find that for a landscape that is gently tilted (slope of order of 10–3), a change in runoff has a minor effect on transport gradient, as sediment transport and associated flexural response maintain topography at a similar elevation. However, there can be a significant change in depositional style when the degree of tilt is altered by, for example, a local change in upper-mantle density. An increase in buoyancy within the upper mantle, which increases slopes, leads to a transient reduction in grain sizes deposited at a fixed location. This behavior is due to a temporary retreat of the zone of erosion into the catchment and a transient increase in accommodation space relative to sediment supply. A reduction in tilt has the opposite effect, the older deposits are eroded, and the erosion-deposition transition rapidly moves down system. There is convincing evidence that the formation of thin and laterally extensive conglomeratic units of the Great Plains was due to a reduced rate of subsidence. Based on the results of our model, we suggest that the deposition of widespread conglomeratic units within continental interiors is generally a consequence of a reduction in slope, as the dynamic support for regions of high topography is reduced.

Sediment accumulation within the continental interior and at passive margins is unsteady and nonuniform, as highlighted by changes in sedimentary facies, the caliber of deposits (e.g., gravels, sands, and silts), and by change in sediment accumulation rates. In some cases, these changes can be linked to clear tectonic or climatic events that affected regional topography or sediment flux, while decoding the reason for change in other records is not so clear (Armitage et al., 2011). For example, the Zambezi Delta succession records an increase in sediment delivery from the Oligocene until the Quaternary (Walford et al., 2005). This change in sediment accumulation has been associated with: the onset of extension in the East African Rift; regional uplift and tilting due to a deep thermal anomaly; and drainage reorganization of the Zambezi River catchment (Walford et al., 2005).

In North America, the New Jersey margin experienced a notable increase in sediment accumulation during the Miocene (Poag and Sevon, 1989; Mountain et al., 2009). The reasons for this increase in sediment accumulation are less clear, but there is evidence of coeval rejuvenated erosion within the Appalachian catchments (Gallen et al., 2013; Boettcher and Milliken, 1994), which could be related to regional uplift driven by mantle flow (e.g., Spasojevic et al., 2008). However, an increase in surface runoff during the Miocene due to regional climate change could also have led to increased erosion and sediment delivery to the New Jersey margin. There is also ongoing debate about the mechanism of deposition of widespread gravel units during the Miocene–Pliocene within the Great Plains of the United States: whether these units signify a change in erosion and deposition due to climatic shifts (Wobus et al., 2010; Tucker and van der Beek, 2013), are linked to long-wavelength tilting of the continental interior (McMillan et al., 2006; Duller et al., 2012), or represent a change in threshold slope due to an autogenic change in runoff (Engelder and Pelletier, 2013).

The widespread deposition of a coarse conglomerate unit within the Spanish Pyrenees during the Paleocene-Eocene transition is temporally linked to an increase in surface runoff driven by an abrupt change in climate (Schmitz and Pujalte, 2007; Armitage et al., 2011; Manners et al., 2013). This raises the possibility that similar deposits of laterally extensive gravel sheets are a signature of change in runoff. However, in northwestern North America, the deposition of gravel units throughout the Cretaceous and Cenozoic has been causally linked to changes in patterns of uplift and subsidence (Heller et al., 2013). Going back further into the geological history of North America, change in the Sloss sequences (Sloss, 1963) within the North American continental platform have been linked to large-scale tilting and subsidence due to large-scale mantle anomalies, which are likely a consequence of subduction (Mitrovica et al., 1989; Coakley and Gurnis, 1995; Burgess and Gurnis, 1995).

The key issue that hinders our attempts to accurately decode the cause of observed increases in denudation and sediment accumulation is that both can be a function of change in climate, and the same change in denudation and accumulation could be caused by tectonically or buoyancy-induced changes in surface uplift.

Erosion of bedrock by flowing water is driven by detachment of rock when river systems are incising and there is no alluvial cover. The CONUS soil data set of Miller and White (1998) for the United States would suggest that large regions of the continental interior are effectively covered in a transportable regolith. We could therefore infer that on a gross scale, erosion is not governed by bedrock detachment, but the transport of this regolith. Pelletier (2011) used the CONUS data set of Miller and White (1998) to make a preliminary estimate of the relative importance of erosion by bedrock detachment and sediment/regolith transport. This work proposed that erosion becomes increasingly limited by the transport of sediment as relief increases (Pelletier, 2011). This inference is contentious, however, because the depth to bedrock also reduces as elevation increases (Miller and White, 1998). Furthermore, bedrock incision is not uniform through time, as sediment cover will inhibit or drive incision (Sklar and Dietrich, 2001). This leaves the open question: Can we assume that erosion at a large spatial and temporal scale is limited by the transport of sediment?

Over days to thousands of years, it is arguable that individual events, storms, and rapidly fluctuating climate must be considered, as the alluvial cover will change the way in which erosion operates (e.g., Lague, 2010). Deposition of widespread conglomerate units such as the Miocene Ogallala Formation within the Great Plains, and change in sediment accumulation at a passive margin, however, occur over durations of more than a million years (e.g., Walford et al., 2005; Cather et al., 2012). A package of stratigraphy that represents a million or more years of deposition holds information about thousands of storm events, as sediment or as time gaps. In essence, when viewed over geological time scales, it is arguable that the multiple individual storm events become averaged out (Paola et al., 1992). Furthermore, over such long time scales, there is evidence that sediment transport across sedimentary systems is buffered to periodic changes in runoff, but it is sensitive to shifts to new climatic or tectonic regimes that last for millions of years (Métivier and Gaudemer, 1999; Blum and Tornqvist, 2000; Castelltort and van den Driessche, 2003). We will therefore focus on exploring how the coupled system responds to a single shift in surface runoff and upper-mantle density.

Our aim is to model erosion and deposition along the length of an ancient sediment routing system such as the Ogallala Formation, which starts at the Laramie Range, Wyoming, and spreads out onto the Nebraskan Great Plains. Given that the majority of such a low-slope sedimentary system traverses a landscape that is covered in transportable sediment or regolith, as for example in the present-day Great Plains (Miller and White, 1998), we will assume that erosion and deposition are controlled by the transport of sediment. We build on the work of Smith and Bretherton (1972), Flemings and Jordan (1989), and Paola et al. (1992) to develop a coupled model of sediment transport with lithosphere flexure. In the early 1970s, Terence Smith and Francis Bretherton published a mathematical framework showing how a model of sediment flux–dependent erosion can re-create realistic morphologies. Two decades later, an early model that coupled deposition and flexure was published where sediment transport was assumed to be a linear function of slope (Flemings and Jordan, 1989). Three years later, Chris Paola published a derivation for the conservation of mass for transport within braided alluvial channels and alluvial fans. The derived diffusion equation is referred to as Exner’s equation of conservation of bed sediment, after Felix Exner (e.g., Exner, 1920; Paola et al., 1992). In deriving the equations for the change in elevation, it was shown that the diffusion coefficient is a function of runoff (Paola et al., 1992).

In this article, we will first present the model equations for sediment transport, grain-size fining, and flexure of the lithosphere. We wish to explore how a landscape will respond to a change in regional topography due to upper-mantle flow, and to a change in surface runoff. To model a change in topography driven by a density change in the upper mantle, we introduce a positive (upward) load on the elastically defined lithosphere. This is a simplified representation of the dynamic support that is believed to be responsible for anomalously elevated mountainous regions of the continental interiors, such as the eastern Rocky Mountains in Wyoming and Colorado (Karlstrom et al., 2012).

We run three different experiments on the coupled system. First, we create an elevated region by introducing a permanent density anomaly below the elastically defined lithosphere. Second, we increase and decrease surface runoff within the model domain by increasing the imposed precipitation rate after a 10 m.y. period of model evolution. Third, we increase and decrease the density of the anomaly driving topographic change in the model after a 10 m.y. period of model evolution. We chose a 10 m.y. initial duration on the basis that it is of a similar order of magnitude to the observed periods between change in sediment accumulation within the continental interior and at the continental margins (e.g., Cather et al., 2012). The results of this study are compared to the record of sediment accumulation across the Great Plains during Miocene–Pliocene times.

2.1. Sediment Transport

Following Dietrich et al. (2003), we begin with a simple idealized landscape composed of bedrock, thickness η (units of m), and a surface layer, regolith, of thickness h (units of m; see Fig. 1). This landscape is forced externally through uplift or subsidence, U (units of m yr–1). Bedrock is transferred into regolith at a rate, P (m yr–1), and regolith is transported across the system with a sediment flux, qs (m2 yr–1). Assuming that the density of regolith produced and transported is equal to the bedrock, within this simple system the rate of change in bedrock thickness is
and the rate of change in regolith thickness is
It then follows that the rate of change in landscape elevation is the sum of the two rates of change:
To solve for the change in surface elevation, we must make a further assumption. One assumption is that the thickness of the regolith remains roughly constant through geological time, ∂th ∼ 0. This is the equivalent to assuming that any newly generated regolith is instantaneously transported out of the model domain. This leads to the rate of change in landscape elevation being
The production of regolith becomes a key consideration when modeling regions where rivers are incising into bedrock and also for exploring soil production and weathering. Yet from studies of the thickness of present-day regolith across the United States (Miller and White, 1998) and assuming that this is representative of most regions of low relief within the continents, it is plausible that there is a supply of transportable regolith that is readily available along the majority of the sediment routing system. We make such an assumption, and to solve for the change in landscape, we carry through the summation in Equation 3, assuming the density of transported material is equal to that deposited:
Equation 5 is a form of the Exner continuity for mass. Using this type of continuity equation and following the derivation of Paola et al. (1992), the change in elevation with time can then be solved by a diffusive equation of the form
where the sediment flux is a function of local slope, and the diffusion coefficient ν is dependent on the water flux, qw (m2 yr–1; all symbols are listed in Table 1). Assuming that bed load is transported following the empirical Meyer-Peter–Müller transport laws (Meyer-Peter and Müller, 1948), then the diffusion coefficient is given by (Paola et al., 1992),
where the constant A = 1 for the case of a meandering river in an alluvial plain, and A = 0.15 for a braided river. Cf = 0.01 is a dimensionless drag coefficient, C0 = 0.7 is the volume concentration of sediment in the bed, and s = 2.7 is the sediment specific gravity. Using these values, the diffusion coefficient is given by ν = 0.10qw in the meandering case and ν = 0.67qw for the braided case. Using a grain-size–dependent critical Shields stress rather than the Meyer-Peter–Müller bed-load transport law, a similar form of the diffusion coefficient (Eq. 7) was arrived at by Marr et al. (2000) for gravel and sand, giving values of the order of ν = 0.1qw for gravel and ν = 1.0qw for sand. Assuming a catchment length of 100 km and a precipitation rate of 1 m yr–1, the diffusion coefficient at the catchment outlet can be estimated to be 104–105 m2 yr–1.

In classic models of foreland basin stratigraphic development, it has been assumed that this diffusion coefficient is constant with length down the system, with ν = 100–5000 m2 yr–1 (Flemings and Jordan, 1989; Sinclair et al., 1991). Depending on the boundary conditions, this simple relationship would lead to rounded or convex-up one-dimensional (1-D) profiles (Métivier, 1999). Taking a constant diffusion coefficient that is uniform in space is appealing, given its simplicity, yet it does not allow for a change in erosion and deposition due to change in water flux down the length of the catchment.

Here, we formulate a model for sediment transport that allows for the runoff to increase down the length of the system, as greater quantities of water will be captured within the fluvial network. Following the work of Smith and Bretherton (1972), we assume that sediment flux is a function of both the slope and the surface runoff,
where c is a constant value that is similar to the product of the constants in Equation 7. This relationship states that sediment transport is the sum of a part that is a constant function of slope and a second term that accounts for the increasing water flow down system. This relationship for sediment flux is similar to that proposed by Smith and Bretherton (1972); however, we have assumed that there is no power relationship between sediment flux and water flux, and between sediment flux and slope, i.e., qsqwn (∂xz)m, where n = 1 and m = 1 . The exponents n and m are dependent on which model of bed-load transport is thought representative of the large-scale transport; for example, Einstein-Brown bed-load transport laws give n = 2 and m = 2, or Meyer-Peter–Müller bed-load transport laws give n ∼1 and m ∼1 (Smith and Bretherton, 1972; Wobus et al., 2010). We have chosen to set n = 1 and m = 1, not because we believe the Meyer-Peter–Müller bed-load transport laws to be more representative of long-term sediment transport, but to make the system equation simpler, such that sediment flux is simply a function of slope and water flux.
The water flux is found assuming a spatially uniform distribution of precipitation rate, α, and then calculating the downhill flow path distance (see Smith and Bretherton, 1972; Simpson and Schlunneger, 2003),

The combination of Equations 8 and 9 leads to a dominantly diffusive equation, where the effective diffusion coefficient is a function of space and increases down system.

In these equations for sediment transport, we do not make the distinction between the region of the landscape that is eroding and that which is undergoing deposition. Equations 8 and 9 have been used to model the erosion of upland catchments for 1-D models of normal fault-bounded catchment-fan systems (Densmore et al., 2007; Armitage et al., 2011, 2013) and two-dimensional (2-D) models of wedge-shape initial topographies (Simpson and Schlunneger, 2003). In the 1-D catchment-fan models, deposition within the sedimentary fan was calculated from a geometric mass balance where there is a continuity of slope between the catchment and fan, at the fan head. However, these equations can equally apply to the alluvial plain and alluvial fans (e.g., Paola et al., 1992; Granjeon and Joseph, 1999). Therefore, to model both erosion and deposition within the continental interior, where there is typically a transportable regolith, we propose that this transport-limited setup is reasonable.

2.2. Flexure

The vertical displacement of the continental lithosphere can be described by the displacement of an elastic to visco-elastic layer, depending on the time scale of observation (Watts et al., 1982; Willett et al., 1985). To understand the first-order response to loading and unloading upon the distribution of sediment, we treat the upper continental lithosphere as an elastic beam of a preset elastic thickness, because the relaxation time scale of continental lithosphere is of the order of a few thousand years (Mitrovica and Forte, 1997), and foreland basin architecture can be adequately reproduced by a model of elastic flexure (Flemings and Jordan, 1990). The displacement, w, of the beam due to loading is then simply given by the fourth-order differential equation,
where Df is the flexural rigidity, pi is the positive (upward) load imposed to simulated uplift due to a density anomaly in the mantle, ρm is the mantle density, ρfill is the density of material eroded and deposited, g is the acceleration due to gravity, and Δz is the change in elevation due to erosion and deposition. We have assumed that the density of what is eroded is equal to that which is deposited. This is therefore assuming no mass is lost to the system through a process such as chemical weathering and the removal of minerals to the ocean. This is therefore an upper estimate for the deposited load exerted on the lithosphere. The flexural rigidity is given by
where E is Young’s modulus, νP is Poisson’s ratio, and Te is the effective elastic thickness. We keep the elastic thickness constant within our model simulations because the changes in elevation due to erosion and deposition are small relative to the elastic thickness.

2.3. Coupling Sediment Transport and Flexure

Erosion and deposition are coupled to flexure by first solving for the changes in displacement due to an imposed load, pi (Eq. 10). The shape of the load is a rectangle, 100 km wide and 50 km thick, positioned at the center of the 2000-km-long domain, and it is described as a density anomaly of between 50 and 100 kg m–3, which forces the lithosphere upward. In sections 3.1 and 3.2, the imposed load is held constant at 50 kg m–3. In section 3.3, we either decrease or increase the magnitude of the load: It is initially100 kg m–3, and after 10 m.y. of model evolution, it is reduced to 50 kg m–3, or it is initially 50 kg m–3, and after 10 m.y. of model evolution, it is increased to 100 kg m–3. The initial elevation, z, is given by the displacement due to the imposed load (Fig. 2). Topography is of an elevated region flanked by depressions due to the flexural response of the imposed load. These flanking basins are deeper for lower lithosphere elastic thickness (Figs. 2B and 2C). The load subsequently changes as mass is redistributed by sediment transport (Eq. 8), where the load due to the sediment transport is ρfillgΔz. The elevated central region is eroded, and deposition occurs within the flanking basins. For internal sinks within the model, water flows down to the lowest point in both the positive and negative x directions. It is then assumed that the water leaves the system, and the eroded or deposited surface is left behind. The removal of this water could be rationalized as transport of that water out of the plane of the 1-D model.

Sediment transport is solved for using a finite-element numerical model. Equations 8 and 9 are made dimensionless, where
and L is the system length. Equations 8 and 9 become
and
where

We use the initial topography or topography from the previous time step, z, to solve Equation 14 for the unknown dimensionless water flux, forumla We solve for the water flux by numerically integrating Equation 14. With the water flux calculated, the length-dependent diffusion coefficient, forumla, can be calculated for each element within the model domain. Now Equation 13 can be solved. The change in topography and displacement is solved implicitly using a finite-element method with linear weighting functions. The time step is 1000 yr on a 1-D grid of 2000 elements.

2.4. Down-System Grain-Size Distribution

We are interested in looking at the way in which change in uplift rate and runoff will alter the stratigraphy in basin successions. We wish to explore the grain-size variation in the deposit, yet the sediment flux is assumed to not be a function of grain size. To a first order, this assumption is likely reasonable (Paola et al., 1992), although in reality, grain size will affect the shear stress required to transport grains down system (e.g., Dade and Friend, 1998). The diffusion coefficient on the Exner equation of mass conservation has been estimated to increase by an order of magnitude once all gravel has been deposited (Marr et al., 2000). However, to maintain the conceptual simplicity of our model, this effect will be ignored.

To explore how stratigraphic grain size may reflect the change in forcing of the system, we use a model of selective deposition by mass (Fedele and Paola, 2007). This model states that greater rates of grain-size fining take place where more deposition occurs. The model solution takes an input grain-size distribution and apportions this input distribution down the basin. For gravel, the grain-size distribution forumla is given by
where forumla is the dimensionless down-system length of deposition, forumla0 is the mean input grain size, φ0 is the variance of the input grain distribution Cv = 0.25 (Armitage et al., 2011), Cg = 0.7 (Duller et al., 2010), and forumla is the spatial transformation forumla of given by (Paola and Seal, 1995)
where forumla is the dimensionless distribution of deposition down system, and forumla is the equivalent down-system distribution of sediment flux. The model of self-similar grain-size fining, Equation 16, implicitly assumes that sediment grain-size distributions are normal. Therefore, the characteristic mean grain size is given by
and the variance is
where D50 is the median grain size, and D84 is the 84th percentile of the distribution of grain sizes. We chose D50 = 40 mm and D84 = 70 mm for our model runs.

3.1. Model Parameters and Model Evolution in the Absence of Change

To understand how the basic system behaves, a model with no change in density anomaly (50 kg m–3) and fixed precipitation rates of either 1 or 2 m yr–1 is presented in Figure 3. We explore the model behavior for two elastic thicknesses, Te = 20 and Te = 80 km. Models of erosion and deposition coupled to lithosphere flexure assumed a diffusion coefficient, ν, of 100–5000 m2 yr–1 (Flemings and Jordan, 1989; Sinclair et al., 1991). In our model, the diffusivity is also a function of down-system length (Eq. 8). To model sediment transport, we use the lower value from previous models and assume κ = 100 m2 yr–1. For the water flux–dependent part, following the derivation of Paola et al. (1992), c in Equation 8 can be related to the channel sinuosity, the drag coefficient on the base of the flow, sediment density, and sediment concentration in the bed (see Eq. 7). We will explore the model behavior with c = 0.1 to c = 0.01, similar to the values derived by Paola et al. (1992) and Marr et al. (2000). Steady state within this model would manifest in constant values of sediment flux across the system and, for these parameter values, we find that it takes more than 50 m.y. for a steady state in sediment flux to be achieved (Fig. 3). Furthermore, after 50 m.y. of erosion, there remains elevated topography at the center of the model domain (Fig. 3).

The longevity of the model landscape is not surprising, given that solutions to the sediment transport equation for the erosion of a simple 10 km ramp-like initial condition approach a steady state after 5 m.y. in the absence of flexure (Armitage et al., 2013). In Armitage et al. (2013), we showed that for the erosion of a ramp of fixed length with fixed boundary elevation, the time to reach steady state scales as forumla, which from Equation 15 can be rewritten as τ ∝ κ/(cα)2. The time scale to reach steady state has a more complex dependence on system length through the solution to a set of simultaneous equations built up of Bessel functions (see Armitage et al., 2013). For example, the time scale associated with a 100-km-long landscape is approximately two times longer than for the equivalent 10-km-long landscape. However, increasing c, which relates the sediment transport to the surface runoff, by one order of magnitude reduces the response time by approximately two orders of magnitude.

In the coupled models, we find that sediment flux decreases substantially during the first 10 m.y. for both values of c (Figs. 3B and 3D). If c = 0.01, a 50% reduction in sediment flux takes 10.5 m.y. when the elastic thickness, Te, is 20 km and 21.9 m.y. when Te is 80 km (Fig. 3B). The equivalent reduction in sediment flux for c = 0.1 takes 1.7 and 2.8 m.y., respectively (Fig. 3D). The long-term trend is of a gradual decline in sediment flux after the initial more rapid reduction in sediment flux. When c = 0.01, an 80% reduction in sediment flux across the model domain takes 50.5 m.y. in the case that Te is 20 km (Fig. 3B). If Te is 80 km, an 80% reduction in sediment flux is not achieved within the 100 m.y. model run. For the case where c = 0.1, the equivalent reduction in sediment flux takes 7.2 and 28.6 m.y., respectively. Doubling the precipitation rate from 1 to 2 m yr–1 increases the response time by a factor of two and creates a long-lived increase in sediment flux relative to the model where precipitation rate is 1 m yr–1 (Figs. 3B and 3D). The characteristic time scale for the denudation of continental landmass is estimated to be on the order of 25 m.y. (Pinet and Souriau, 1988). This suggests a lower value of c may be more appropriate. Furthermore, a lower value for c creates topography that has a low concavity, which is appropriate for exploring landscape evolution across continental interiors.

3.2. Signals Due to a Change in Surface Runoff

To explore how the idealized system evolves, we introduced a smooth transition in precipitation rates after 10 m.y. of model evolution (Fig. 4A). We explore the response for two elastic thicknesses, Te = 20 and Te = 80 km (solid and dashed lines, respectively, in Fig. 4). We also explore the model response for c = 0.1 and 0.01 (see Eq. 8).

3.2.1. Sediment Flux across the Model Domain

For all model runs, there is an increase in sediment flux following an increase in precipitation rates (Fig. 4). Sediment flux then slowly reduces; however, given the long response time of the coupled model, within the 10 m.y. period, the system does not recover to the pre-perturbed state (Figs. 4B and 4C). For the case where c = 0.01, the magnitude of the sediment flux response is larger than if c = 0.1 (Figs. 4B and 4C). This is because when c = 0.1, the landscape has lower gradients due to the effectiveness of sediment transport at removing mass (Fig. 3C). A reduction in precipitation rates leads to a reduction in sediment flux, with the system similarly shifting to a new prolonged state of gradually reducing sediment flux (Figs. 4D and 4E). Again, when c is larger, the magnitude of the response is lower due to the overall reduction in slope following increased erosion.

The long-lived state of increased/decreased sediment flux described here is caused by the interplay between the change in topography driven by the mantle density anomaly and the change in load due to erosion and deposition. As precipitation rate is changed, erosion and flanking deposition change. This redistribution of mass continuously modifies the distribution of rock uplift and acts to maintain similar gradients across the 1-D landscape (Fig. 5). The result is that sediment flux remains at an elevated value for long periods of time.

The change in precipitation rate is being imposed before the system has achieved a steady state. The coupled model takes more than 50 m.y. to achieve a significant reduction in sediment flux (Fig. 3). In previous models of erosion driven by similar transport laws, but where uplift is imposed as a vertical velocity rather than an instantaneous adjustment to a buoyant load upon an isostatically compensated lithosphere, a single step increase in precipitation rate generated a spike in sediment flux that recovered to the same steady-state sediment flux prior to the change (Armitage et al., 2011). Increasing precipitation for that same model but before the model had recovered to steady state resulted in a reduced signal in terms of sediment flux out of the catchment that nonetheless recovered to the same steady-state signal of sediment flux out of the catchment (Armitage et al., 2013). In the model we present in this paper, increasing precipitation leads to a long recovery time, and even after 20 m.y., the sediment flux signal will not attain similar values for different precipitation rates (Fig. 3). Such a long recovery of sediment flux is not predicted by previous models (Armitage et al., 2011, 2013).

3.2.2. Stratigraphy

Throughout the full 20 m.y. period for all models, there is a general trend of progradation (Figs. 6 and 7). This progradation is a function of the system slowly evolving toward a steady output of sediment flux (Figs. 3 and 4). The response recorded within stratigraphy of an increase in precipitation rate is difficult to observe without a close inspection of Figure 6 or 7.

For the case where c = 0.01, upon the increase in precipitation rates, there is a thickening over time of depositional units and a gradual increase in depositional length above the background rate at 10 m.y., associated with the increase in sediment delivery (Figs. 6B, 6C, 6G, and 6H). The lack of a strong signal of stratigraphic progradation, such as that predicted within the model of Armitage et al. (2011), is because accommodation keeps pace with sediment supply. This is because, within the coupled model, the flexural response is a combination of: (1) the load of the sediment, (2) the erosion, and (3) the imposed density anomaly, which work together to increase accommodation space generation such that there is no significant change in the rate of grain-size fining. The stratigraphic response to a reduction in precipitation rates is a thinning of the stratigraphic units, but likewise there is no strong response within the granulometry, other than a minor reduction in the rate of progradation at 10 m.y. (Figs. 6A, 6B, 6D, and 6E).

For the case where c = 0.1, there is an increase in transport of sediment due to the flow of water, qw. This causes increased erosion and deposition such that the landscape becomes quite flat within 20 m.y. (maximum slope of 2.6 × 10–4 at 20 m.y.; Figs. 3C and 7). We plotted the stratigraphic record for an elastic thickness of 20 km only in Figure 7; as for a 80 km elastic thickness, the resulting record shows the same trend but with thinner deposits. The stratigraphic evidence for change in precipitation rate is migration of the erosion-deposition transition (Fig. 7). Otherwise, there is little evidence within the stratigraphy of the change in precipitation.

For all of the models at young ages, sediment is deposited outside of the main basin due to wider flexural bulges that are shallow. These wider deposits are very thin and are eventually buried beneath the main basins. Sediment is also transported into the basins from both directions, leading to a deposit of coarser material along the far edge of the main basin. This deposition is analogous to that expected within the hanging wall of a fault-controlled basin.

3.3. Signals Due to Change in Upper-Mantle Buoyancy

We introduced an increase or decrease in topography by changing the magnitude of the density anomaly in the upper lithosphere (Fig. 8A), which causes tilting of the landscape. Again, we modeled two elastic thicknesses and two values of c (see Eqs. 8 and 10).

3.3.1. Sediment Flux across the Model Domain

Increasing the magnitude of the density anomaly that maintains the elevated central region in the model domain from 50 to 100 kg m–3, for both values of c (0.01 and 0.1), causes sediment flux to initially increase and then decrease as the system re-equilibrates (Figs. 8B and 8C). This behavior is similar to that predicted for an increase in runoff, where after the initial perturbation, the interplay between load and topography causes a gradual reduction in sediment flux as the system evolves. The key difference between a change in density anomaly and runoff is that slopes change significantly when the density anomaly is changed (Fig. 9).

For the case where c = 0.01, the response to a reduction in buoyancy in the upper mantle is a reduction in sediment flux (Fig. 8D). When the elastic thickness is 20 km, there is an additional short-lived (∼1 m.y.) relative drop in sediment flux (Fig. 8D, solid line). This minimum is due to a transient period where the imposed reduction in density driving topographic change creates a platform-like elevated region, or plateau, in the center of the model domain (Fig. 10A, 11 m.y.). Sediment flux increases briefly as the edges of this plateau are eroded at between 10 and 11 m.y. The model then evolves toward a landscape with a central peak of reduced elevation after 20 m.y. (Fig. 10A).

For the case where c = 0.1, a reduction in buoyancy creates an increase in sediment flux off the eroding regions within the model domain (Fig. 8E). This is because when the effect of water flux on sediment transport is larger, the reduction in the density anomaly from 100 to 50 kg m–3 creates a central depression (Fig. 10B). As topography lowers due to the reduced magnitude of the density anomaly, the flexural bulge due to the flanking basins becomes responsible for the generation of the greatest elevation at a model distance of 750 and 1250 km (Fig. 10B, 11 m.y.). The whole structure then inverts, as central deposition and flanking erosion act to lower topography until the landscape is almost flat (maximum slope of 1.7 × 10–4 at 20 m.y.; Fig. 10B). This inversion of topography is responsible for the increase in sediment flux when the density anomaly is reduced (Fig. 8E).

3.3.2. Stratigraphy

The stratigraphic response to a change in topography due to the mantle density anomaly is quite different from that to a change in runoff. For an increase in relief due to an increase in the buoyancy, a retrogradation of the depositional system and of grain size is recorded (Figs. 11B, 11C, 11G, 11H, 12B, 12C, 12G, and 12H). This phase of retrogradation and increase in grain-size fining is transient, however, resulting from a temporary reduction in sediment supply relative to accommodation space. Accommodation space increases due to the flexure of the lithosphere as a consequence of increased upper-mantle buoyancy. As can be seen within the chronostratigraphic diagram (Figs. 11 and 12), the depositional front has an up-system trajectory at 10 m.y., and then deposition gradually migrates back down system, with an associated lengthening of the depositional system. This signal is stronger when c = 0.1 (Fig. 12).

A reduction in the density anomaly that provides the buoyancy-driven support of the landscape creates a very different response, which is strongly dependent on the strength of sediment transport as a function of water flow (Figs. 11A, 11B, 11E, 11F, 12A, 12B, 12E, and 12F). For the case where c = 0.01, a reduction in buoyancy causes the depositional system to prograde (Figs. 11A, 11B, 11E, and 11F). This progradation is a consequence of a reduction in elevation, which causes a forward migration of the depositional front. This pushes the coarse deposits forward. In the case of a 20 km elastic thickness, a short-duration (∼1 m.y.) coarse unit is deposited upon the central elevated region and subsequently eroded away (Figs. 11E and 11F).

If the effect of runoff is stronger, c = 0.1, then there is inversion, which is to say that the centrally elevated region becomes a depocenter (Figs. 12A, 12B, 12E, and 12F). This is due to the density anomaly being of too small a magnitude to maintain the central elevation. Instead, the flanking deposits create the highest elevation due to flexure of the lithosphere at a distance of roughly 750 and 1250 km (Figs. 10B, 12A, and 12E). Deposition then switches into the central basin, and the flanking basins become abandoned (Figs. 12A, 12B, 12E, and 12F).

We have presented a model for the transport of sediment to calculate the change in topography across a 2000 km region for a change in relief driven by change in density within a 100-km-wide by 50-km-thick region in the upper mantle. Parameter values for the model of sediment transport are justifiable based on up-scaling empirical bed-load transport laws and are similar to previous models that couple deposition and lithosphere flexure (e.g., Flemings and Jordan, 1989; Sinclair et al., 1991; Paola et al., 1992). These parameter values lead to a topography that has a low concavity and hence a low relief, which is appropriate for exploring processes in the continental interior away from large mountain belts. The spatially changing load due to erosion and deposition alters topography as the modeled purely elastic lithosphere adjusts isostatically. This model suggests the following:

  • (1) In our model, the amount of material transported by the flow of water is controlled by the parameter c (Eq. 8). The value of c can be estimated from the basic properties of alluvial sediment transport and is roughly between 1 and 0.01 (Paola et al., 1992; Marr et al., 2000). We have explored the lower end of this range and found that for c = 0.1–0.01, the response time is between 1 and 20 m.y. (Fig. 3). For c = 0.1, landscape becomes relatively flat after 20 m.y., and the magnitude of change in sediment flux, following a change in precipitation or buoyancy in the upper mantle, is small (Figs. 4, 8, and 10). If the we increase c to 1, as suggested for sand transport (Marr et al., 2000), then topography that is generated by buoyancy within the upper mantle would be eroded down over a shorter period of time, <1 m.y., as the response time is inversely proportional to c, τ ∝ κ/(cα)2 (see section 3.1).

    In the case where the transport of sediment due to water flow is weak, c = 0.01, then the time scale of response to change is close to that estimated for continental denudation, 25 m.y. (Pinet and Souriau, 1988). It is likely that c is not a fixed value in space or time as sediment gets moved and the distribution of gravel and sand changes; however, our modeling study suggests that on a gross scale, if c > 0.1, landscape may be effectively beveled off, while if c < 0.1, elevated regions will remain elevated for more than 100 m.y. (Fig. 3). Therefore, for modeling long-term sediment transport, we suggest that c ∼0.01 is more reasonable.

  • (2) A change in runoff due to an increase or decrease in precipitation rate causes an increase or decrease in sediment flux out of the region of erosion (Fig. 4). The change in erosion and deposition affects the surface load, which facilitates isostatic uplift that keeps pace with the denudation. The result is that a change in runoff causes only a minor change in catchment slope (Fig. 5). The depositional system thickens and gradually lengthens as the system evolves, and, consequently, grain-size fining does not vary significantly (Figs. 6 and 7).

  • (3) An increase in relief driven by an increase in upper-mantle buoyancy causes an increase in sediment flux (Figs. 8B and 8C). This increase in sediment flux is accompanied by a transient phase of retrogradation of the depositional system, shown as a vertical reduction in grain size (Figs. 11C, 11D, 11G, 11H, 12C, 12D, 12G, and 12H). The transient retrogradation phase is due to the increase in the rate of accommodation space generation relative to sediment flux. Accommodation space is generated in the model by the flexural bulge that flanks the uplifted central region (Figs. 2 and 3). The instantaneous flexural response means that for greater contrast in density, the elevation of the central region of the model and the amount of subsidence in the flanking basins become greater. This change increases the accommodation space for sediment at a pace that is more rapid compared to the rate of increase in sediment delivery to the basin.

  • (4) A reduction in the buoyancy that maintains the elevated central region causes marked progradation of the depositional front if the transport of sediment due to runoff is low, c = 0.01 in Eq. 8 (Figs. 11A, 11B, 11E, and 11F). Progradation is due to the drop in topography creating a reduction in slope, so that the region of positive curvature migrates outward. Elevated topography is due to the imposed buoyancy and also due to the flexural response of the flanking load. The combination of the buoyancy and flexure causes the central elevated region to widen as deposition migrates outward. If sediment transport due to water flow is high, c = 0.1, then the reduced buoyancy cannot maintain the central high, and the system inverts, creating a central basin where there was previously an elevated region (Figs. 10B, 12A, 12B, 12E, and 12F). The flexural response to the flanking basins creates two elevated regions on either side of a central depocenter. This central basin then starts to fill as the topography flattens.

These numerical experiments demonstrate that the stratigraphic signature of change in tilt of the continental interior due to mantle flow is delicately controlled by the strength of erosion and deposition due to sediment transport, and mediated by the lithospheric response.

4.1. Comparison with Previous Transport-Limited Models

Earlier models that use a similar approach for sediment transport but with a different mechanism for creating change in topography have predicted both similar and different potential records of sediment accumulation. The response of this model to a change in relief due to a change in the density anomaly that drives topographic change has similarities to the previous short normal fault–controlled sedimentary fan development models of Paola et al. (1992), Densmore et al. (2007), and Armitage et al. (2011). This similarity is due to the change in fining being due to a similar shift in the ratio of sediment supply to accommodation: In section 3.3, after an increase in catchment elevation driven by increased buoyancy, accommodation space increases faster than supply. The predicted signals left in the stratigraphic record due to a change in runoff are, however, different. In the normal fault–bounded mountain catchment-fan model of Armitage et al. (2011), an increase in runoff is predicted to generate prograding conglomeratic sheet-like deposits with sediment fluxes reducing to steady-state values after a million years. However, in section 3.2, we found that for large systems, where topography change is by a flexural response to change in load from both the upper mantle and surface, a change in runoff generates only a minor signal within the granulometry accompanied by a prolonged (>10 m.y.) increase in sediment fluxes.

There are three key differences between the model developed here and the previous models of fault-bounded catchment-fans that reduce the impact of change in precipitation on the sedimentary record:

  • (1) Response times are very long in this coupled model (Fig. 3). The increased response time is due in part to the choice of parameters (κ and c in Eq. 8). The values of these parameters are based on reasonable estimates of basic physical properties of bed-load transport and are similar to previous numerical models. It is clear that increasing c or decreasing κ will decrease the initial model response time to change in precipitation. However, the long-term gradual decline in sediment flux out of the central elevated region is a function of the interplay between unloading and loading due to sediment transport. This keeps slopes elevated and allows the continued transport of material above an equilibrium value for at least 50 m.y. (Fig. 3).

  • (2) A change in precipitation does not cause a significant change in slope (Fig. 5). Increasing precipitation rates, for example, increases the load at the flanking basins and reduces the load in the central high. The instantaneous flexural response causes the eroded area to rise and the deposited regions to sink. This feedback between the removal of mass and the rebound of the surface topography keeps slopes roughly similar as precipitation is increased. Therefore, sediment flux increases and does not reduce rapidly as is the case for models where flexure of the lithosphere and isostasy were ignored (e.g., Armitage et al., 2011).

  • (3) Accommodation space is not fixed or controlled by tectonic faulting. This allows the basin to decrease and increase in size as the sediment loading changes. The result is that if sediment delivery to the basin is increased, then the accommodation space will likewise increase. Thus, a change in precipitation has very little impact on the stratigraphic record in terms of vertical granulometry (Figs. 6 and 7).

Within the construct of our model, the deposition of coarser grains in the form of a temporarily uniform far-traveled conglomeratic sheet only occurs upon a reduction in the mantle density anomaly driving tilt of the surface (Fig. 11). Such behavior is similar to that proposed by Heller et al. (1988). In Heller et al. (1988), it was suggested that widespread/far-traveled conglomerate units are deposited as mountain building ends. This hypothesis was immediately questioned based on depositional ages within the northwest Himalaya (Burbank et al., 1988), where coarser deposits prograde further down system with time and are related to rejuvenated uplift within the axial zone of the Himalaya. What we find is that an increase in elevation due to upper-mantle buoyancy does produce a down-system migration of larger grain sizes, but this gradual progradation is a symptom of the system evolving toward a steady state, rather than a direct signal of an increase in topography within the eroding landscape. We propose that coarse and laterally extensive gravel deposits are most likely a result of a reduction in catchment uplift and basin subsidence, or a reduction in tilt.

4.2. Late Cenozoic Erosion and Deposition in Southwestern North America

During the latest Cretaceous to Paleocene, a number of far-traveled conglomerate units were formed in southeastern North America. These deposits lie above a disconformity and travel down the length of the basin (Heller et al., 2013). Subsidence analysis of the basins that contain these deposits would locate deposition occurring after a period of subsidence increase (see Heller et al., 2013, their Fig. 4). In other words, deposition of gravel units occurs once subsidence decreases. In the more recent geological past, there is evidence for two periods of change in denudation within the southwestern United States during the late Cenozoic, which might be related to the deposition of far-traveled conglomeratic units and could be a function of change in climatic conditions or due to long-wavelength tilting due to mantle flow:

  • (1) During the late Oligocene to early Miocene (27–15 Ma), widespread erosion and fluvial incision occurred from the Colorado Plateau to the southern Great Plains and central and western Texas (Chapin, 2008; Flowers et al., 2008; Cather et al., 2012). This period of denudation was synchronous with extension and footwall uplift within the Basin and Range and Rio Grande Rift, and it ended with the widespread deposition of coarse-grained units of the Ogallala, Bidahochi, and Fence Lake Formations (Cather et al., 2012). The late Oligocene to early Miocene also corresponds to a period of significant change in ocean circulation, with the closure of the Tethys Ocean between Europe and Africa, and increased deep-water formation through the Faroes-Scotland Ridge (see review by Chapin, 2008).

  • (2) During the late Miocene to Pliocene (6–3 Ma), there is clear incision of the Ogallala deposits and evidence of tilting of the pre-incised surface (McMillan et al., 2002; Duller et al., 2012). The incision of the Ogallala deposits may have been driven by increased surface runoff (Wobus et al., 2010), and the tilt may have subsequently been a consequence of that erosion as the lithosphere isostatically compensated for the change in surface load. However, it could be that the tilt was due to change in the density structure of the mantle and crust, associated with warmer mantle and ignimbrite eruptions along the Jemez lineament (Wisniewski and Pazzaglia, 2002; Nereson et al., 2013), and the period of uplift led to incision of the Ogallala deposits.

We will apply our model to the Ogallala Formation, to explore if the change in depositional slope and change in gravel units deposited are more likely due to change in runoff or tilting of the continental interior. The Ogallala Formation can be split into four units of between 50 and 100 m thickness and 250–300 km length, each of duration of 4 m.y. Assuming each unit has a cross-sectional area in the shape of a triangle, the sediment flux required to deposit the Ogallala is between 1.5 and 3.8 m2 yr–1. This magnitude of sediment flux is consistent with that generated by our model of sediment transport where c = 0.01 (Figs. 3B and 3D).

The late Cenozoic Ogallala Formation is perhaps thicker than the earlier Cretaceous–Paleocene deposits that are associated with a reduction in subsidence (Heller et al., 2013). Yet it is plausible that their deposition marks the tail end of regional uplift within the Rio Grande Rift zone and the Colorado Plateau. The Ogallala deposits (18–6 Ma; Swinehart et al., 1985) are potentially correlated with volcanism and northward propagation of the Rio Grande Rift zone (e.g., McMillan et al., 2002). However, their deposition postdates peak volcanism within Colorado and New Mexico (ca. 35 Ma; McMillan et al., 2000) and the onset of extension within the Rio Grande itself (Chapin and Cather, 1994). Recent apatite (U-Th)/He data would suggest that extension within the northern and southern Rio Grande Rift was coeval, and there was no northward propagation of the rift zone (Landman and Flowers, 2013). It is therefore possible that regional uplift decreased after ∼15 m.y., after the main phase of extension, which is supported by thermochronometric measurements that suggest the Colorado Plateau has experienced little or no change in elevation since this time (Huntington et al., 2010). Therefore, as uplift reduced, the depositional system migrated downstream, causing incision of the upper deposits and the progradation of coarse units onto the Great Plains.

In the last 10 m.y., there has been further change in the topography of the southwest United States (Figs. 13A and 13B). Observations would suggest that the Ogallala deposition surface was tilted at or before 6 Ma, as the present-day slope of the Ogallala Formation is steeper than the reconstructed transport slope of the Ogallala Formation (Leonard, 2002; McMillan et al., 2002; Duller et al., 2012; Fig. 13A). From outcrop patterns of the Ogallala and Remsburg Ranch units, we can infer that the flexural hinge is located ∼160–200 km east of the Wyoming-Nebraska border (Swinehart and Diffendal, 1995, 1997; Duller et al., 2012). The transport slope of the Ogallala during the formation of this layer is similar to the present-day slope of the North Platte River (Duller et al., 2012; Fig. 13A). This phase of uplift may have been associated with increased temperatures in the mantle that led to uplift of the surface, and rejuvenated volcanism within the Jemez lineament (Nereson et al., 2013). This uplift may also have been related to the Aspen seismic anomaly to the south (Karlstrom et al., 2012). This 100–300-km-wide anomalous low-seismic-velocity zone within the upper mantle is associated with a low Bouguer gravity anomaly, which would suggest a buoyant crust and upper mantle support the high topography. Karlstrom et al. (2012) proposed from thermochronologic and geologic data that regional exhumation accelerated starting ca. 6–10 Ma, particularly within regions like that above the Aspen low-velocity zone. This would suggest that Neogene mantle convection has driven long-wavelength surface deformation and tilting over the past 10 m.y. (Karlstrom et al., 2012).

It is estimated that the surface underwent up to 600 m of increased elevation to the west (Leonard, 2002; McMillan et al., 2002; Duller et al., 2012), and the surface was then incised prior to deposition of the Broadwater Formation (Duller et al., 2012). If we assume that the lithosphere has an elastic thickness of 80 km, which gives a flexural rigidity of 4.55 × 1024 Nm and is comparable with estimates of Leonard (2002), and that c = 0.01 to be consistent with estimated sediment flux, then a decrease in buoyant support of topography by reducing the density anomaly from 100 to 50 kg m–3 generates a slope change of roughly 4 × 10–3 within 3 m.y. of model evolution (Fig. 13C). Such a change in slope is comparable to that estimated to have had occurred between 6 and 3.7 Ma (Duller et al., 2012). The model suggests that a significant change in slope due to an increase in surface runoff is, however, not possible (Fig. 5).

Based on our idealized model of sediment transport coupled to an isostatically compensated lithosphere, we propose that the observed late Cenozoic pattern of deposition within the Great Plains was primarily driven by change in catchment uplift as the eastern Rocky Mountains were gently tilted by density changes within the mantle. Over the last 30 m.y., two pulses of transient uplift have left behind the Ogallala Formation and its subsequent incision. The counter hypothesis, that this history in deposition is a consequence of change in runoff, is not consistent with our model of sediment transport coupled to lithosphere flexure.

We have developed a model of sediment transport coupled to an isostatically compensated lithosphere. It is contentious to suggest that the whole of the continental interior is covered in transportable regolith; however, given that sediment covers large proportions of the United States (Miller and White, 1998), we suggest that this transport-limited model is appropriate for modeling gross change in deposition across the continental interior. Based on this model, we propose the following:

  • (1) In the absence of tectonically controlled basin formation, change in runoff does not have a strong signature in the stratigraphic record.

  • (2) Change in topography due to change in uplift rates causes signals of progradation or retrogradation for a reduction or increase in uplift.

When we apply these model results to the history of the southwestern United States, we propose that the deposition of gravel conglomerate deposits that span millions of years, such as the Ogallala Formation, is due to a drop in catchment uplift and basin subsidence. For the case of the Ogallala Formation, this drop in catchment elevation was likely due to the ending of the Cenozoic (27–15 Ma) period of extension in the Rio Grande Rift zone. Uplift within the eastern Rocky Mountains decreased at roughly 15 Ma, and, consequently, deposition migrated outward onto the Great Plains and southward, creating the Ogallala Formation.

More recent (>6 Ma) uplift and incision of the Ogallala are more likely due to a recent phase of uplift of the Great Plains associated with a change in density within the upper mantle. There is strong evidence for buoyancy within the upper mantle below regions such as Aspen, Colorado, which is a likely source of such buoyancy-driven support of high topography. Such a change in topography due to mantle buoyancy is more likely to have promoted incision of previous deposits, as our model would suggest that erosion due to increased runoff has little effect on topography and transport slope.

This work was funded by a Marie Curie research fellowship (project number 272669) and has continued support through a Royal Astronomical Society Research Fellowship to John Armitage. We would like to thank Peter Burgess and Philip Allen for their comments on this manuscript. We would also like to thank Chris Paola, Greg Tucker, and Lithosphere Editor Eric Kirby for very constructive and thorough reviews.