Sensitivity of rift tectonics to global variability in the efficiency of river erosion.

SignificanceThe efficiency of erosion in leveling relief mainly depends on climate and strength of exposed rocks. However, whether erosion is sufficiently efficient to influence the architecture of a tectonic plate boundary remains a topic of debate. Here, we analyze continental rift landscapes reworked by river incision to assess a globally representative range of fluvial erosion efficiency. We then simulate crustal extension exposed to surface processes acting within this documented range. We find that more efficient erosion favors the growth of half-grabens over horsts, which can explain contrasting tectonic styles across the Basin and Range province and the East African Rift. This suggests that variability in Earth's geological structures partly reflects variability in hydrological conditions and associated surface processes.

Deformation at tectonic plate boundaries is driven by a combination of far-field tractions and gravitational forces that arise from changes in crust and lithosphere thickness, temperature, and surface topography (1, 2). Topography is, in turn, reworked by erosion and sedimentation, prompting the development of models that coupled tectonics and surface processes (3–12). As numerical simulations improved and became able to resolve the development of narrow fault zones in the upper crust (Fig. 1), a consensus formed around the fact that leveling topography through efficient surface processes helps focus brittle deformation onto fewer faults, leading to narrower orogens (4, 8) or larger offsets on rift-bounding faults (5–7, 9, 13–15). Several outstanding questions remain, however, beginning with whether the efficiency of surface processes varies sufficiently in space and time to modulate the localization of brittle deformation. Most modeling studies to date have tested the effect of a wide range of landscape evolution parameters on a theoretical spectrum from no erosion/sedimentation to complete leveling of relief, but few have focused on the sensitivity of tectonic style to erosional efficiencies within a range documented by geologic observations. Moreover, it is unclear whether natural landscapes can encode a record of surface processes within this natural range of variability. This has proven especially challenging in orogens, where deformation involves multiple simultaneously active thrusts and reflects force and mass balances on a scale much larger than that of individual faults (1, 16, 17).

Fig. 1.

Tectonic and surface processes relevant to half-graben growth. (A) Shuttle Radar Topography Mission topography of the Lemhi range (Basin and Range, United States) and adjacent sedimentary basin shaped by slip on a southwest-dipping master normal fault (52). masl: meters above sea level. (B) Two-dimensional numerical simulation (run H65) of half-graben growth in visco–elastoplastic crust coupled to sedimentary infilling in subsiding areas (yellow) and fluvial incision in uplifting areas. In this snapshot, the master fault (F1; region where nonrecoverable strain ε exceeds critical weakening strain ε) has accumulated 6 km of throw in the 16-km-thick brittle upper crust. The white line shows the 420 °C isotherm, which effectively marks the brittle–ductile transition. The three gray bands initially lay flat at 5-, 10-, and 15-km depths and now show internal deformation of the fault-bounded blocks.

These challenges can be addressed by focusing on half-graben settings within continental rift zones. These relatively simple, regional-scale (<200-km) structures form when a large fraction of extensional strain is taken up by a single master normal fault (18–21). Slip on the master fault upwarps and downwarps the footwall and hanging wall blocks, respectively, thereby shaping a small mountain range and basin (Fig. 1). While some asymmetric basins can grow as half-graben structures for well over 10 My (e.g., refs. 18–22), theoretical work suggests that a master fault cannot grow indefinitely, as the buildup of flexural and topographic stresses in the adjacent blocks eventually triggers the formation of a new crustal-scale fault (23–25). Such relocalization of strain has been documented in the field both on the hanging wall (26) and on the footwall side (27).

The maximum amount of slip that can be accommodated by a half-graben master fault depends on the strength of the brittle lithosphere and can be enhanced by suppressing relief growth with infinitely efficient surface processes (7). Most half-grabens undergo surficial mass redistribution through glacial and fluvial processes that strip material from the uplifting footwall and deposit sediments onto the subsiding hanging wall (28). However, surface processes in natural systems are never efficient enough to completely level the relief caused by active rock uplift. Sediments typically fill most of the accommodation space that develops in the hanging wall, such that the topographic depression remains on the order of approximately hundreds of meters deep, sometimes covered by a lake. On the other hand, footwall relief can vary greatly depending on the efficiency of erosional agents, which is itself modulated by climatic conditions and exposed lithology (29–31). Our approach is, therefore, to assess whether geomorphologically plausible erosional efficiencies can affect the maximum amount of extension taken up by a single half-graben structure and influence its subsequent evolution. We do so by focusing on settings primarily reworked by river incision (32) and observables such as the vertical component of fault offset, hereafter referred to as fault throw (Fig. 1), which help quantify the strain accommodated by individual half-grabens.

Coupled Simulations of Rifting and Surface Processes

We simulate the coupled evolution of a half-graben system undergoing surface processes in a 2.5-dimensional model domain (Fig. 1), ignoring any thermal or mechanical contributions from magmatic processes that often accompany rifting ( has details). A 30-km-thick crustal layer overlying the mantle is modeled as a visco-elastoplastic continuum in the vertical (x, z) plane subjected to a far-field horizontal extension rate of 1 mm/y (33, 34). The coldest, shallowest portion of the crust is effectively elastoplastic and can spontaneously form faults through softening of its cohesion and friction over a characteristic amount of plastic deformation ε = 0.5 when the Mohr–Coulomb criterion is met (25). To localize strain onto an initial master fault (fault F1), a narrow 55°-dipping weak zone is initially seeded in the upper crust (Fig. 1). Deeper, hotter regions of the domain deform viscously according to experimentally determined plagioclase (crust) and olivine (mantle) flow laws (35, 36). The strong temperature dependence of these rheologies enables us to investigate a wide range of crustal strength profiles by changing the geotherm through a spatially uniform crustal radiogenic heat production rate, which we vary between 0.35 and 1.4 μW/m3 (). The warmest geotherm we consider produces a thin (∼11-km) faulted upper crust overlying a low-viscosity lower crust and upper mantle. The coldest geotherm produces a 22-km-thick upper crust; a strong, viscous lower crust; and a viscous upper mantle withstanding stresses ≥10 MPa to as deep as ∼40 km. These end-member cases, hereafter referred to as “weak” and “strong” crust, sample the range of strengths inferred from the depth distribution of seismicity across rift systems () but exclude strength profiles that favor the formation of very large offset detachment faults and associated core complexes (e.g., a very thin brittle upper crust transitioning abruptly to a low-viscosity lower crust) (2, 24, 25). We note that because of their regional focus, our simulations emphasize processes related to crustal thinning and cannot fully capture long-wavelength uplift induced by mantle upwelling.

Deformation of the domain’s free surface is coupled with a river incision model in the horizontal (x, y) plane following the approach in ref. 14 and the implementation in ref. 37. The efficiency of river incision is varied through the coefficient of erosion K of a stream power law with slope and drainage area exponents of 1 and 0.5, respectively (38, 39). To enable comparison with natural landscapes that may have formed under different tectonic rates, we hereafter refer to the simulations’ dimensionless erosional efficiency defined as (analogous to 1/N in ref. 40 and typically between ∼0.5 and ∼5; see below), where A0 is a reference area arbitrarily set to 106 m2 and U0 is the fastest rock uplift rate in the landscape, typically occurring at the master fault. Subsiding areas are instantaneously filled with low-density sediment such that their elevation remains at zero altitude, except in a set of reference simulations with low erosional efficiency.

Simulations with varied erosional efficiency and crustal strength (Figs. 2 and 3) were run until a new master fault F2 develops spontaneously (i.e., a new crustal-scale shear band accumulates a plastic strain greater than or equal to ε) (Fig. 2). Slip on fault F1 leads to elastoplastic flexure of the hanging wall and footwall blocks, where secondary faults (plastic strain less than ε) rapidly develop (Fig. 1). This effect, along with rapid fault rotation toward gentler dips (41), leads to an effective throw rate (the vertical component of the fault slip rate) of ∼0.7 mm/y on F1. Footwall uplift and basin subsidence are asymmetric, with ∼36% of the throw rate contributing to range uplift and ∼64% contributing to sedimentary basin development in the early stages of deformation. We measured a maximum footwall uplift rate of 0.25 (±0.05) mm/y on F1, during the early stages of fault growth (i.e., first ∼3 km of throw) across all of our simulations. We, therefore, used U0 = 0.25 mm/y to convert our models’ coefficients of erosion into EEs. As river networks develop on F1’s rising footwall, erosion tends to approximately balance rock uplift, and mean footwall relief ceases to grow (Fig. 4). This dynamic equilibrium is attained with a lesser fault throw when the coefficient of erosion is increased and corresponds to a more subdued footwall relief. After a few kilometers of slip on F1, a secondary fault has typically accumulated enough strain (greater than or equal to ε) to be considered a new master fault F2, at which point F2 starts creating local relief away from the initial half-graben (Fig. 2). We note that F1 continues to slip when F2 becomes active, but past this point, it can no longer be characterized as a half-graben master fault (Fig. 2).

Fig. 2.

Localization of a new master fault modulated by crustal strength and erosional efficiency. Snapshots of half-graben growth models taken after a new master fault F2 has broken spontaneously. Total throw accommodated on the initial master fault F1 is indicated. Symbols and colors match those in Fig. 1. Panels correspond to different thicknesses of the brittle upper crust (19 km in panels A and B; 16 km in panels C and D) and erosional efficiency (EE = 0.4, A and C; EE = 4, B and D). Panels correspond to simulations (A) F70, (B) F60, (C) H70, and (D) H60 as listed in .

Fig. 3.

Maximum vertical offset on half-graben structures and their subsequent tectonic evolution modulated by erosional efficiency and crustal strength. Throw accumulated on the master fault F1 when a new master fault F2 breaks in all numerical simulations. Squares, triangles, and circles indicate that new master faults broke in the hanging wall of F1, in the footwall of F1, or on both sides of F1, respectively. Labels 1B and 2A–D refer to runs displayed in corresponding panels of Figs. 1 and 2. (Left) These points are from simulations with no sedimentary infilling. Documented range (0.19 to 5.87) includes data from nine half-grabens and three horsts.

Fig. 4.

Documenting the erosional modulation of rift structure. (A) Half-graben relief as a function of fault throw. Black curves show the relief of F1's footwall in simulations with a brittle-layer thickness of 14 km (runs M50, M55, M60, M65, and M70) and varied erosional efficiency. Dots indicate when F1 can no longer be considered the dominant fault (F2 reaches ε). Dashed blue and red curves correspond to simulations with EE = 1.26 and brittle-layer thicknesses of 22 km (run C65) and 11 km (run V65), respectively. Horizontal bars correspond to nine half-grabens in the Basin and Range province, Rio Grande Rift, Taupo Rift, and East African Rift (blue, black, and red colors indicate brittle-layer thickness of >20 km, between 10 and 20 km, and of <10 km, respectively). (B) River elevation vs. upstream distance corrected for variations in drainage area and uplift rate from the nine half-grabens shown in A. A footwall range in topographic equilibrium plots along a line of slope = 1/EE. Green lines indicate the EE values used in our numerical simulations. (C) Vertical offset on half-graben master faults vs. erosional efficiency, with the same color code as in A. Colored bands show the maximum throw on F1 as a half-graben master fault in simulations of varied brittle–ductile transition depths. (D) Transition between the hanging wall snapping (squares) and the footwall snapping (triangles) regimes in numerical simulations (black curve). Hollow rectangles correspond to half-grabens from A–C. Filled triangles indicate horst structures, and the filled rectangle is a possible example of hanging wall snapping (north basin of Lake Malawi).

The maximum throw that can be accommodated by F1 as a half-graben master fault increases with EE at a rate strongly modulated by crustal strength (Fig. 3). In the case of weak crust, a 100-fold increase in EE increases maximum throw by a factor of ∼2.3, while in strong crust, the increase is only by ∼1.4. The location of F2 also appears jointly influenced by erosional efficiency and strength (Fig. 3). Simulations with more efficient erosion and strong crust tend to localize F2 in the hanging wall of F1 with the same dip direction (“hanging wall snapping” regime in ref. 25) (Figs. 2 and 3), while simulations with inefficient erosion and weak crust localize new synthetic (F2) and antithetic faults (F2’) in the footwall of F1 (“footwall snapping” regime) (Figs. 2 and 3). F2’ typically ends up dominating F2 in its topographic and structural expression, forming a new asymmetric basin that delineates a horst block (Fig. 2 ). Finally, we ran a set of reference simulations with low EE (0.4) and no sedimentary infill (Fig. 3, Left). In these “empty basin” cases, the initial half-graben never exceeds throws of 4 to 5 km, with little influence from the strength of the crust. Filling the basin while keeping EE low increases this maximum throw to 4.5 to 7.5 km, with a more pronounced effect in stronger crust.

In the theoretical framework of refs. 23 and 24, sustaining slip on F1 requires energy supplied by the work of far-field tectonic forces. This energy is dissipated by frictional resistance on F1 and plastic yielding of the adjacent blocks, while a fraction of it is converted to elastic strain energy and gravitational potential energy (GPE) associated with topographic growth and Moho uplift (2, 7, 41, 42). The amount of force required for F1 to slip thus increases with increasing throw as flexural stresses increase and as footwall relief develops. F2 breaks in intact crust when doing so requires a lesser force than continuing to slip on F1 (6, 25). When erosion and sedimentation roughly balance rock uplift and subsidence, respectively, topographic GPE ceases to increase with increasing fault slip, and relief growth no longer contributes to the force increase. The rate of increase in force, therefore, diminishes with greater fault slip, and the threshold to break F2 is attained later (7). In strong crust, the main driver of the force increase is flexure (24, 25), and the GPE associated with topography buildup is a lesser fraction of the total energy cost of sustaining fault slip. Its modulation by surface processes, thus, has a lesser effect on the evolution of the system. Conversely, the deformation of weak crust is more sensitive to the alleviation of topographic GPE (7).

Assessing the Efficiency of River Erosion

For erosion to influence natural rift tectonics, the mechanisms revealed by our simulations must be active over a range of erosional efficiency that is representative of Earth conditions. This range can be constrained by analyzing the elevation profile of rivers draining half-graben footwalls (Fig. 1). According to the stream power model (38, 39), the river incision rate at some upstream distance x from base level equals the erosion coefficient K(x) multiplied by the upstream drainage area A(x) and local slope S(x), respectively, raised to powers m ∼ 0.5 and n ∼ 1. Thus, when a dynamic equilibrium is attained between erosion and uplift U(x), the local slope can be expressed as S = (U/KA)1/. To first order, river slope increases with increasing upstream distance because A(x)—a proxy for river flow rate—sharply decreases as one approaches the drainage divide. Integrating slope upstream thus yields a concave-up river profile (40, 43). Because A(x) is readily measured in digital elevation models, river concavity not explained by changes in drainage area may be attributed to spatial variability in the uplift field and/or the coefficient of erosion. This idea forms the basis of our methodology for estimating EE ( has details).

For simplicity, we assume that K is spatially uniform at the scale of a half-graben footwall and write the spatially variable uplift rate field as U(x) = U0 × U(x). U0 is the uplift rate at the fault, and U(x) is a dimensionless shape function that represents the decay in uplift away from the master fault. In a half-graben footwall, U*(x) is well described by a thin plate flexure model. Such models involve a characteristic decay length that can be straightforwardly estimated from topographic analyses (e.g., ref. 18). Specifically, we fit U*(x) to the strike-averaged, normalized footwall topography that dips away from the fault because its shape plausibly encodes the decaying shape of the rock uplift field. This cannot be said for the steep, fault-facing scarp of the footwall range, the base level of which is essentially pinned at the fault, where rock uplift is fastest. For this reason, we restrict our analysis to footwall rivers that drain away from the master fault. Armed with an informed estimate of U*(x) and direct measurements of A(x) (but no knowledge of U0), we calculate a transformed upstream distance, corrected for changes in uplift rate and drainage area: . This variable change from x to χ is useful because the slope of river elevation plotted against χ happens to be exactly 1/EE (Fig. 4). Erosional efficiency can, therefore, be assessed without knowledge of local tectonic rates or climatic conditions, provided U(x) is well approximated. Benchmarking our method against the synthetic landscapes produced by our simulations suggests that this is the case and that the decay of footwall relief is an appropriate proxy for the decaying uplift pattern.

We apply our method to nine half-graben footwalls from four rift systems (), yielding the river elevation–χ plots shown in Fig. 4. Scatter in the data leads to considerable overlap across sites; however, clear EE signatures exist (Fig. 4 ). For example, the weaker lithology and wetter climate of the Paeroa Range in the Taupo Rift of New Zealand results in a much higher EE (2.81; with lowest/highest estimates: 1.99/3.02) than the Santa Rosa border fault in the Basin and Range (EE = 0.62; 0.51/1.11), despite a slip rate that is approximately three times faster (32). Our estimates of EE are contingent on whether topographic equilibrium has been reached in these settings. This is most plausible in Basin and Range half-grabens, where total fault throw (>4 km) greatly exceeds mean footwall relief (∼1 km) (Fig. 4) (29). A conservative estimate of documented EE variability restricted to Basin and Range faults is 0.62 to 2.66, while a broader estimate encompassing all sites is 0.51 to 5.87 (Fig. 4). This range of EE in our simulations brackets the observed trend of footwall relief vs. fault throw (Fig. 4). Specifically, the models and the data are consistent with a rapid increase in relief over the first ∼2 km of throw followed by stabilization at 800 to 1,200 m. The decrease in relief past a peak height in our models reflects the fact that as strain begins to accumulate on secondary faults (Fig. 1), the throw rate on F1 decays, and EE accordingly increases.

Sensitivity of Tectonic Styles to Erosional Efficiency

For the documented range of EE, our simulations provide a reasonable upper bound on half-graben fault throws (Fig. 4). Models with thin brittle crust yield maximum throws of 5 to 7 km, compatible with the measured throws of Basin and Range half-grabens (Fig. 4). This upper limit is sensitive to EE (e.g., increasing EE from 1.26 to 4 results in an ∼33% increase in maximum half-graben throw in 14-km-thick upper crust) (Fig. 3). By contrast, models with a strong crust result in maximum half-graben throws of 7 to 8 km (Fig. 4), not strongly affected by EEs between one and four. This is consistent with the estimated throw of ∼7 km on the Livingstone fault, which borders the north basin of Lake Malawi (44), has an EE ∼ 3, and breaks through a ∼35-km-thick seismogenic layer (45) (Fig. 4). Interestingly, the seismic moment release rate of the Livingstone fault is roughly matched by that of its intrabasin faults (45). This could be a signature of incipient or ongoing strain migration into the hanging wall of a fault that has reached its maximum throw (44), with the important caveat that seismic observations cover too short a time span to provide conclusive evidence. At any rate, basin-ward migration of strain in the north basin of Lake Malawi would be consistent with hanging wall snapping behavior in our simulations (Fig. 2), which has been documented in several other rift settings (26). Lake Malawi features the thickest brittle crust and the largest erosional efficiency in our compilation, both of which favor the hanging wall snapping regime outlined by our parametric study (Fig. 4).

Instances of footwall snapping are best documented in horst structures, which are often asymmetric and feature a dominant border fault (e.g., refs. 27 and 46), suggesting that horst-bounding antithetic faults initiate sequentially rather than simultaneously. We complement our half-graben dataset with river profile analyses from three horst blocks (), assuming that the spatial pattern of rock uplift is controlled by the dominant fault. The Rwenzori horst, which stands out as a structural and topographic anomaly along the western branch of the East African Rift (46), is characterized by a thick brittle crust (25 to 30 km) (47) and the lowest EE in our compilation (0.28; 0.19/0.46)—possibly reflecting the combination of a rapid rock uplift rate and strong gneissic lithologies. By contrast, horst blocks such as the Ruby Mountains from the Basin and Range have a documented history of strain relocalization on antithetic faults (27) and likely developed in brittle crust about half as thick as the Rwenzori horst. We measure significantly greater EEs in horst blocks from the Basin and Range (0.72 to 2.60), consistent with the regime diagram outlined by our simulations. In summary, while footwall snapping takes place across all documented erosional efficiencies in weaker crust, it requires very inefficient erosion to occur in stronger crust (Fig. 4).

Conclusions and Perspectives

Overall, the natural variability that we document in erosional efficiency reflects variability in the strength of exposed lithologies, climatic conditions, and tectonic deformation rates. Erosion will be less efficient over a faster-slipping fault system exhuming stronger lithologies, and larger offsets will be necessary for erosion to balance uplift. The result will be greater footwall relief, resulting in larger topographic stresses that prompt the formation of new faults in the footwall, possibly turning half-grabens into horsts. It is of course possible that a range of other processes may have impacted the evolution of each of our study sites on a case by case basis. Magmatic intrusions or underplating, for example, likely influenced deformation and uplift in the Rwenzori Mountains/Edward basin region (48) and the volcanically active Taupo Rift. Our modeling also ignores structural inheritance, which can profoundly influence the architecture of rifts. Interestingly, our results also open the possibility of lithological inheritance (i.e., modulating the tectonic evolution of a system by exposing more or less erodible lithologies).

While our simulations emphasize the influence of footwall erosion, they also suggest an important role of sediment deposition in the hanging wall basin, which can increase the maximum throw of a half-graben by as much as 50% compared with an empty-basin end-member (Fig. 3). Because our 2.5-dimensional modeling approach precludes any along-strike variation in tectonic uplift, our simulations cannot develop fault tips or evolving segments, which form distinct sedimentary basins along the axis of continental rifts (49). Our approach cannot account for sediment transfer between these basins, which typically supplies more sediment to an individual basin than its associated footwall range can contribute (50, 51). Our assumption of complete infilling of subsiding areas, however, captures this effect to first order. We did not explore the role of partial infilling of hanging wall basins, which is relevant for underfilled endorheic basins disconnected from along-axis sediment transport. We expect that such systems may experience greater initial subsidence of the valley floor, increasing the GPE cost of half-graben growth, as seen in our simulations with empty basins (Fig. 3, Left). Interestingly, progressive rift growth and segmentation, as well as volcanic activity and climate change, can alter the sediment routing system over time, draining the catchments of overfilled basins into underfilled ones (51). While our purposefully broad modeling approach cannot tackle each of these rift-specific effects, our findings broadly imply that naturally occurring variability in lithological and climatic conditions exerts a detectable modulation on the tectonic makeup of Earth’s plate boundaries.