The impact of earthquakes on orogen-scale exhumation

Individual, large thrusting earthquakes can cause hundreds to thousands of years of exhumation in a geologically instantaneous moment through landslide generation. The bedrock landslides generated are important weathering agents through the conversion of bedrock into mobile regolith. Despite this, orogen-scale records of surface uplift and exhumation, whether sedimentary or geochemical, contain little to no evidence of individual large earthquakes. We examine how earthquakes and landslides influence exhumation and surface uplift rates with a zero-dimensional numerical model, supported by observations from the 2008 M_w 7.9 Wenchuan earthquake. We also simulate the concentration of cosmogenic radionuclides within the model domain, so we can examine the timescales over which earthquake-driven changes in exhumation can be measured. Our model uses empirically constrained relationships between seismic energy release, weathering, and landsliding volumes to show that large earthquakes generate the most surface uplift, despite causing lowering of the bedrock surface. Our model suggests that when earthquakes are the dominant rock uplift process in an orogen, rapid surface uplift can occur when regolith, which limits bedrock weathering, is preserved on the mountain range. After a large earthquake, there is a lowering in concentrations of ¹⁰Be in regolith leaving the orogen, but the concentrations return to the long-term average within 10³ years. The timescale of the seismically induced cosmogenic nuclide concentration signal is shorter than the averaging time of most thermochronometers (>10³ years). However, our model suggests that the short-term stochastic feedbacks between weathering and exhumation produce measurable increases in cosmogenically measured exhumation rates which can be linked to earthquakes.

1 Introduction

Surface uplift of a mountain range is controlled by the balance of additive uplift processes and removal of material by surface, typically fluvial, processes. Earthquakes produce rock uplift and equally importantly generate exhumation via landsliding (Avouac, 2007; Keefer, 1994; Marc et al., 2016a). Landsliding events can scour deeply into bedrock, causing many hundreds or thousands of years of exhumation of the bedrock surface in a geologically instantaneous moment (Fig. 1) (Li et al., 2014; Marc et al., 2019; Parker et al., 2011; Stolle et al., 2017). Existing mass balances on single (Parker et al., 2011) or sequences (Li et al., 2014, 2019; Marc et al., 2016b) of earthquakes demonstrate that landslide volumes of large thrust earthquakes are comparable to, and may exceed, rock uplift. However, earthquakes are not the only rock uplift process in mountain belts; aseismic mechanisms such as viscous and elastic crustal deformation (Meade, 2010; Simpson, 2015), lithospheric delamination (Hales et al., 2005; Molnar et al., 1993), and isostatic rock uplift (Molnar, 2012; Molnar et al., 2015) can also contribute. The total time-averaged mass balance, and any net surface uplift, is also affected by erosion between earthquakes (Hovius et al., 2011; Marc et al., 2019; Yanites et al., 2010). Therefore, the contribution of earthquakes to the generation of long-term surface uplift of mountains remains poorly constrained (England and Molnar, 1990; Li et al., 2019).

Figure 1 Landsliding can generate large volumes of sediment but only the largest consistently transport sediment off the hillslope. The volume of bedrock eroded is also dependent on the thickness of regolith on the hillslope before the landslide occurs. Panel (A) shows the 2016 Aranayake landslide in Sri Lanka which was triggered by intense rainfall in heavily weathered soils. The failure occurred along the regolith–bedrock interface and very little bedrock was eroded. Panel (B) shows that in bedrock-dominated areas, such as the Usoi Dam in Tajikistan (produced by the 1911 earthquake), significant volumes of regolith can be generated by landsliding but little measurable erosion as the regolith produced has remained in the catchment as a landslide dam for over 100 years. The 2008 Wenchuan earthquake produced over 60 000 landslides (C), many of which have significantly smaller transport lengths than the hillslope (highlighted in red) and so cannot be easily removed from the catchment. For these reasons, we define landsliding as weathering rather than as erosional agents. Panel (A) is sourced from https://blogs.agu.org/landslideblog (last access: 3 July 2020), (B) is from https://www.mergili.at/worldimages/picture.php?/9330 (last access: 3 July 2020), (C) is sourced from Google Earth, using imagery provided by Landsat/Copernicus and Maxar Technologies, 2019 and draped over a digital terrain model. © Google Earth.

Despite the importance of earthquake-triggered landslides in the total erosion budget of mountain belts, there is little evidence of large earthquakes in most sedimentary or exhumation records. Increased rates of sedimentation and changes in sedimentary characteristics, linked to the abundance of loose sediment, have been identified in lakes and reservoirs immediately proximal to large faults (Howarth et al., 2012; Zhang et al., 2019). However, these pulses in sedimentation can be difficult to separate from extreme hydrological events (Zhang et al., 2019). Steep slopes within mountain belts, as well as typical observations of plentiful exposed bedrock on those slopes, have been used to support the idea that earthquake-generated sediment is rapidly removed from orogenic belts (Dingle et al., 2018; Li et al., 2014; Niemi et al., 2005; Parker et al., 2011). However, observations collected after earthquakes suggests that regolith, transportable sediment produced by landsliding or weathering, remains in low-order catchments (Pearce and Watson, 1986) (Fig. 1a), single large landslide deposits (Korup and Clague, 2009; Stolle et al., 2017), landslide dams (Fan et al., 2012; Ouimet et al., 2007) (Fig. 1b), and valley fills (Blöthe and Korup, 2013) for up to 10⁴ years. An estimated 44 % of the sediment produced in the Himalaya is stored in the long term in some form before it exits the mountain range, demonstrating the importance of storage and sediment recycling to orogen-scale sediment fluxes (Blöthe and Korup, 2013). Similar findings are found in studies on post-earthquake sediment fluxes: in New Zealand, up to 75 % of the sediment produced by the 1929 Murchison earthquake remained in catchments 50 years after the earthquake (Pearce and Watson, 1986), while in Taiwan between 92 % and 99 % of the sediment produced by the Chi-Chi earthquake remained in catchments 8 years after the earthquake (Hovius et al., 2011). The slow removal time of earthquake-derived sediment suggests residence times could exceed recurrence times of earthquakes in tectonically active mountain ranges. If significant volumes of sediment remain in the landscape for times longer than the recurrence period of earthquakes, large earthquakes could contribute significantly more to surface uplift than currently assumed (Li et al., 2014; Marc et al., 2016b; Parker et al., 2011). Constraining the long-term contribution of earthquakes to a mountainous orogen is impossible without fully constraining the rate at which coseismic landslide deposits are removed from the mountain range and the contribution of non-seismic erosion and uplift. The rock uplift minus the coseismic landsliding produced by earthquakes is commonly assumed to be equal to surface uplift, which directly assumes there is no sediment cover in the mountain range, and thus all landslides are in bedrock (England and Molnar, 1990; Li et al., 2014, 2019; Marc et al., 2016b). If sediment remains in orogens for extended periods of time, it is likely it will be remobilised, reducing the erosion and increasing the surface uplift of an earthquake. Remobilisation of coseismic sediment can also occur between earthquakes, altering sediment fluxes and estimates of erosion rates from cosmogenic nuclides (Andermann et al., 2012; Dingle et al., 2018; Yanites et al., 2009). Constraining any aseismic uplift and how it affects coseismic uplift is also be important for understanding whether earthquakes build mountain ranges (Hovius et al., 2011; Li et al., 2014; Parker et al., 2011; Royden et al., 1997). Observations of the 2008 Wenchuan earthquake and its aftermath give us insight into the export of landslide-generated regolith from mountain belts. Over 60 000 landslides were produced by the earthquake, with a total volume of approximately 3 km³ (Li et al., 2014) (Fig. 1c). In mountain ranges, there are significantly more small landslide deposits than large ones, such that the magnitude and frequency of these volumes follows a power law (Malamud et al., 2004a; Marc et al., 2016a; Stark et al., 2001). Only the largest, or most mobile, landslides are most likely to deposit sediment directly into the channel; for the Wenchuan earthquake, less than half of the total volume of the regolith produced is connected to the channel network immediately after the earthquake (Li et al., 2016) (Fig. 1c). Mapping of the landslide deposits reveals much of the landslide material produced by the earthquake remains on the hillslope and in small-order channels, and many deposits are undistributed by erosional processes (Fan et al., 2018; Zhang et al., 2016). When looking at an orogen as a whole, landsliding only transports sediment a very small distance, and hence we define landsliding as a weathering rather than erosional process.

The rate at which landslide deposits can be evacuated from a catchment is typically related to the capacity of the fluvial system to transport the influx of sediment (Croissant et al., 2017; Yanites et al., 2010). If a significant proportion of sediment produced by an earthquake is not able to be eroded by the fluvial system, the residence time of the sediment is likely to be increased due to the need for stochastic hillslope processes, such as debris flows, to remobilise the sediment into the channel before it can be exported (Bennett et al., 2014; Croissant et al., 2019; Hovius et al., 2011). The timing of debris flow triggering is related to the interaction of storms and slope hydrology, which cannot be easily predicted, while the volume remobilised by debris flows is primarily controlled by the non-linear process of sediment entrainment during run out (Horton et al., 2019; Iverson, 2000). The rapid stabilisation of landslides after earthquake, combined with the stochastic triggering of debris flows, could be one reason why large pulses of sediment associated with single earthquakes are rare or absent from downstream sinks (Fan et al., 2018; Zhang et al., 2019). Rather than a sink receiving a large impulse of sediment, which can be easily recorded via a change in average grain size or sedimentation rate, the rate change is instead smeared or shredded by stochastic sediment transport processes like flooding or debris flows (Jerolmack and Paola, 2010). The averaging time for different measures of erosion rate (e.g. cosmogenic vs. thermochronometric) may strongly affect the probability of measuring a single earthquake. If the recording time of erosion (10³ years for cosmogenic radio nuclides or 10⁴⁻⁵ years for thermochronometry methods) is similar to or less than the return period of large earthquakes, then any difference between short- and long-term erosion rates could be due to the influence of earthquakes (Kirchner et al., 2001; Ouimet, 2010). By investigating the variation of erosion rates with varying timescales, or with coseismic landslide density (Niemi et al., 2005; West et al., 2014), we may be able to identify the long-term impact of earthquakes on orogen-scale erosion rates.

In this paper, we use a zero-dimensional volume balance model to explicitly track earthquake generated sediment through time in a hypothetical orogen based upon the Longmen Shan region. We use the tracked sediment thickness in order to understand how earthquake-triggered landslides (EQTLs) affect exhumation and surface uplift at orogenic scales. Our model co-varies the amount of aseismic uplift in the orogen, imposed earthquake magnitude–frequency relationships, and both the timing and maximum magnitude of earthquakes, under multiple possible uplift regimes, in order to fully investigate the role of earthquakes in orogen-scale volume budgets. We then use these scenarios to investigate whether earthquakes can be identified from erosion or exhumation records over different timescales and by modelling cosmogenic radionuclide concentrations of sediment leaving the orogen. Finally, we test our hypothesis that the variation in recorded erosion rates in a landscape can be an indicator of seismic activity using a global database of cosmogenic radionuclide derived erosion rates.

Definitions of terms

This paper is precise in its use of terminology around the changes in elevation of the various surfaces we discuss. These follow standard modern definitions (England and Molnar, 1990), but as we are explicitly measuring the generation of regolith and the movement of two surfaces (the topographic surface and the bedrock surface), the potential for ambiguity requires us to clearly define these terms:

• Erosion is the transport away of material from a site and thus a change in the topography. In our context, it describes the full evacuation of regolith material from the model domain, i.e. removal of sediment from the entire mountain range.

• Weathering is the in situ conversion of rock into regolith. In our model, rock must be converted to regolith before it can be transported. We explicitly separate the role of EQTLs generating regolith by weathering bedrock from their role as (inefficient) eroders of bedrock and regolith – i.e. we separate their role in producing loose material from their role in transporting that material. EQTLs on average occur on hillslopes near ridge tops, typically with transport lengths less than hillslope lengths, with only the largest impinging on the fluvial system (Li et al., 2016).

• Regolith is the mobile transportable layer of sediment at the surface. In the model, regolith can be created by two distinct weathering mechanisms: landslides cutting into bedrock to create transportable debris, and soil production by physiochemical processes.

• Uplift is the increase in elevation of a material or surface in an absolute frame of reference. We distinguish rock uplift, topographic surface uplift, and bedrock surface uplift. Rock uplift is the expected increase in the bedrock surface before considering erosion, it is either produced by coseismic or aseismic means. Bedrock surface uplift is the vertical distance moved by the bedrock surface after erosion. Topographic surface uplift is the total vertical distance moved after considering changes in the bedrock surface and the thickness of the regolith layer.

• Exhumation is the approach of the rock mass towards the topographic surface and/or the bedrock surface, in the frame of reference of that surface. Our model tracks both surfaces; therefore, it is possible to have a rock uplift event that causes exhumation relative to the bedrock surface without exhumation relative to the topographic surface, by thickening the regolith.

• Denudation is used almost as a synonym for exhumation, but where the frame of reference is the rock mass or the bedrock surface, and the topographic surface moves towards it.

2 Methods

2.1 Zero-dimensional volume valance model

Here, we present a generalised zero-dimensional mountain volume balance model which we use to test the impact of regolith storage on bedrock surface uplift and exhumation. In the absence of sufficient empirical evidence on the long-term spatial distributions of rock uplift, exhumation, and regolith volumes in mountain ranges, we simulate these parameters by treating the evolution of a landscape as a series of dimensionless seismic volume balances.

In our model, we define the change in topographic surface elevation (S_T) with time as

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{T}}}{\mathrm{d}t}}=U-E,\end{array}$$

where U (units of length/time) is the thickness of rock entering the orogen during a time step and resulting in rock uplift (calculated as the volume entering the orogen across the area of the model per unit time), while E is the rate of regolith removed (length/time) from the topographic surface and thus is the long-term erosion rate of the orogen. Rock can enter the orogen in two ways: either via shortening and thickening of the crust during coseismic thrusting earthquakes (U_co) or through one of a number of aseismic uplift mechanisms (U_as), such as lower crustal flow (Royden et al., 1997) or lithospheric delamination (Hales et al., 2005). Coseismic deformation (Bonilla et al., 1984; Wells and Coppersmith, 1994), and hence the volume of rock uplifted, scales as a function of earthquake magnitude (Li et al., 2014; Marc et al., 2016b). The addition of mass to a column of crust by thickening will produce an isostatic compensation which will reduce the overall surface uplift response to rock uplift. In our model, we apply a simple compensation based upon the relative densities of the crust and mantle to account for the isostatic response (Densmore et al., 2012; Molnar, 2012; Turcotte and Schubert, 2002). The calculated response is applied immediately to the volume balance and the surface uplift. We do not consider interseismic strain as an uplifting mechanism in this model due its limited contribution to permanent surface deformation (Avouac, 2007). There is ongoing debate to how much surface uplift can be attributed to aseismic vs. coseismic sources and how they interact (Hubbard and Shaw, 2009; Royden et al., 1997). Acknowledging the complexity of the debate, we simplify the aseismic component of uplift and generalise it as the proportion of uplift that cannot be accounted for by U_co. Hence, topographic surface uplift can be represented as

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{T}}}{\mathrm{d}t}}=U_{\mathrm{co}}+U_{\mathrm{as}}-E,\end{array}$$

where the ratio between coseismic uplift rate (U_co) and aseismic uplift rate (U_as) is defined by the term α which represents the proportion of the total uplift rate that is caused by aseismic uplift, such that $(\mathrm{1}-\mathit{\alpha })U=U_{\mathrm{co}}$ and αU=U_as.

In our zero-dimensional model, the thickness of regolith (R) removed from the surface of the orogen is defined as

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{CLRP}}{\mathrm{d}t}}+{\displaystyle \frac{W}{\mathrm{d}t}}-E,\end{array}$$

where CLRP (coseismic landslide regolith production) is the average thickness (length/time) of weathered material generated by coseismic landslides, all of which is assumed to be transportable, and W (length/time) is the thickness of rock weathering caused by all other mechanisms (simply the thickness of material removed from the bedrock when there is no regolith cover). W is included to ensure that erosion can continue even when regolith is not present. In our model, weathering does not occur when there is a covering of regolith as the background weathering rates for our study site in the Longmen Shan region are unknown. This way of including weathering in our model allows it to be an emerging property rather than a fixed rate. The rate of elevation change (length/time) of the bedrock surface (S_B) can now be described as

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{B}}}{\mathrm{d}t}}=U_{\mathrm{co}}+U_{\mathrm{as}}-{\displaystyle \frac{\mathrm{CLRP}}{\mathrm{d}t}}-{\displaystyle \frac{W}{\mathrm{d}t}}.\end{array}$$

We can now define our rate of surface uplift again as a combination of the bedrock surface elevation and the regolith thickness:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{T}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{S}_{\mathrm{B}}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}}.\end{array}$$

The model represents the average topographic surface uplift, regolith generation, and bedrock surface lowering for the area (A) of coseismic displacement of the largest possible earthquake for a fault found within a mountain belt. The length of the modelled area (L) is set by the length of the surface rupture on the fault that generates the maximum earthquake, while the width is the distance to the estimated line of zero strain based upon the dip of the modelled fault (θ) and the focus depth (D) (Li et al., 2014).

$$\begin{array}{}\text{(6)}& A=L\cdot \left({\displaystyle \frac{D}{\mathit{\theta }}}\right)\end{array}$$

As surface uplift rates for mountain ranges are hard to define (England and Molnar, 1990), we set the model to a flux steady state, where U is set to equal the long-term erosion rate (E). For each time step in the model, an earthquake with M_w>5 is randomly chosen from a power law distribution, and the coseismic rock uplift volume associated with this earthquake is calculated using an empirical scaling relationship (Li et al., 2014) between magnitude and rock uplift volume:

$$\begin{array}{}\text{(7)}& \log \left(V_{\mathrm{u}}\right)=\mathrm{1.06}(\pm \mathrm{0.22})M_{\mathrm{w}}-\mathrm{8.40}(\pm \mathrm{1.44}),\end{array}$$

where V_u is the volume of rock uplift generated by an earthquake of magnitude M_w. This volume is scaled by α and divided by the model area A to calculate U_co. M_w 5 earthquakes are the smallest that regularly produce coseismic landsliding so represent the smallest earthquakes of interest to our study (Marc et al., 2016b). We use an optimising algorithm to fit the uplift produced by Eq. (7) to ensure the model remains in a flux steady state. The algorithm uses the uncertainty within Eq. (7) to fit the model parameters so that the average uplift produced during a time step is equal to the long-term erosion rate. The use of the aseismic uplift scaling parameter α has the effect of increasing the time-averaged rock uplift of time steps of small earthquakes and decreasing the rock uplift of large earthquakes.

Regolith is generated in the model based on calculations of bedrock lowering by CLRP and weathering by other mechanisms. We use the scaling between landslide volume (V_l) and earthquake magnitude proposed by Malamud et al. (2004b) (Fig. 2a)

$$\begin{array}{}\text{(8)}& \log V_{\mathrm{l}}=\mathrm{1.42}{M}_{\mathrm{w}}-\mathrm{11.26}(\pm \mathrm{0.52}),\end{array}$$

V_l is converted to a depth of landsliding by dividing by the area of the model space (A). The area of the landscape affected by landsliding of the largest earthquakes is greater than the model space, so a scaling is applied based on the landslide density of the Wenchuan earthquake. Alternative models of coseismic landslide volume as a function of seismic moment (Marc et al., 2016a; Robinson et al., 2016) cannot easily be scaled into a zero-dimensional model space due to their reliance upon earthquake source depth and landscape metrics. These models describe the relationship between earthquake magnitude and total landslide volume as a curve around a hinge magnitude. The shaking produced by an earthquake correlates with the length and width of its surface rupture; however, the width (depth) of the rupture is limited by the thickness of elastic crust. At a maximum magnitude (∼M_w 6.75), the scaling between earthquake magnitude and shaking changes resulting in a curved relationship between total landslide volume and earthquake magnitude. This has the effect of reducing the importance of large earthquakes in the surface uplift balance. As our chosen model does not include this threshold, the larger earthquakes of this model will be more erosive than other models. All empirical models relating coseismic landslide volume and earthquake magnitude have large uncertainties in them. We acknowledge these uncertainties by applying a normal distribution of error using the uncertainty bounds on the landsliding volume produced by each earthquake (Fig. 2a), reducing the difference between the different models. The total new regolith generated by coseismic landslides is then calculated as the difference between the average depth of landsliding and the average thickness of the regolith.

$$\begin{array}{}\text{(9)}& \mathrm{CLRP}=\left({\displaystyle \frac{V_{\mathrm{l}}}{A}}-R\right)\end{array}$$

Figure 2 The average depth of regolith produced by an earthquake is impacted by the earthquake magnitude and the thickness of regolith that is on the hillslope before the earthquake occurs. (a) Scaling of landslide volume with magnitude from Malamud et al. (2004b) and the average residence time of landslide sediment in the Longmen Shan region based upon an erosion rate of 0.62 mm yr⁻¹. The two lines represent the minimum and maximum volumes of landsliding generated, within the bounds on Eq. (8). (b) Variability of regolith production, expressed as volume per area, with existing depth of regolith on the hillslope, for four representative earthquake magnitudes. Coloured areas represent the variability of the landslide volume produced by an earthquake, randomly sampling within the bounds of Eq. (8). These coseismic landslide regolith production functions (CLRPFs) emerge from the model rather than being set in advance, and the variability at each magnitude is driven by noise inherent in the relationship between magnitude and landslide size (Eq. 8). The inset shows the probability distribution function for regolith thickness across the whole model run, integrating the effects of the CLRPF through time. A small but non-zero spatially averaged modal regolith thickness persists, but significantly larger thicknesses regularly occur. (c) Variability of surface elevation through time for model runs with identical earthquake sequences but with varying additional proportions of aseismic uplift. The assumption of steady state prevents any long-term permanent uplift, and so all variations in surface elevation are driven by the sequence of earthquakes and changes in regolith thickness. Increasing the aseismic contribution to uplift reduces the uplift of large earthquakes, resulting in much less variable surface uplift and therefore exhumation.

2.2 Model implementation: Longmen Shan

We test our model using the Longmen Shan region due to the wealth of studies of the area both prior to and after the 2008 Wenchuan earthquake. These studies allow for the small number of parameters in our model to be well constrained. The Longmen Shan region marks the eastern margin of the Tibetan Plateau and the western edge of the Sichuan Basin (G. Li et al., 2017). Hillslopes are at their threshold steepness with pervasive bedrock and limited channel storage (Li et al., 2014; Ouimet et al., 2009; Parker et al., 2011). The front of the range is dissected by three parallel dextral-thrust oblique-slip thrust faults, two of which, the Yingxiu-Beichuan and Pengguan faults, ruptured during the 2008 Wenchuan earthquake (Densmore et al., 2010; Liu-Zeng et al., 2009). Prior to this earthquake, geodetic measurements recorded limited shortening rates suggesting a possible role for lower crustal flow in driving surface uplift (Clark et al., 2005; Kirby et al., 2003; Royden et al., 1997). However, significant shortening associated with the Wenchuan earthquake supports an important, possibly exclusively coseismic, surface uplift element (Hubbard and Shaw, 2009).

To apply the model to the Longmen Shan region, we use a power-law relationship of earthquake frequency and magnitude derived from historical earthquake records (Z. Li et al., 2017):

$$\begin{array}{}\text{(10)}& N>M_{\mathrm{w}}=\mathrm{3.93}-\mathrm{0.91}{M}_{\mathrm{w}}.\end{array}$$

N is the number of earthquakes that occurs above a certain magnitude in a year. The smallest earthquake we model is a M_w 5, which occurs on average every 5 years. This relationship gives a return time of 1816 years for earthquakes of the same magnitude as the Wenchuan earthquake. Other studies (Densmore et al., 2007; G. Li et al., 2014, 2017) have proposed a return time of anywhere between 500 and 4000 years for a M_w 7.9 earthquake. We use the uncertainty in the frequency of Wenchuan-sized earthquakes to vary Eq. (9) to test the impact of earthquake frequency on exhumation and surface uplift. We use an apatite fission track-derived exhumation rate of $\mathrm{0.62}(+\mathrm{0.14}-\mathrm{0.08})$ mm yr⁻¹ (G. Li et al., 2017) to represent the long-term sediment flux from the orogen (E). The model area is set by Eq. (6) using parameters derived from the Wenchuan earthquake. The length (L) is the surface rupture of the Wenchuan earthquake (240 km), the focal depth (D) was between 14 and 18 km deep, and the dip angle ranged between 40 and 90^∘ (Li et al., 2014). These parameters yield an area (A) of 2600 km². We run the model for 25 Myr to allow multiple analyses over various timescales and vary α between 0 and 1.

2.3 Exhumation calculations within the model

We calculate exhumation for 2 kyr intervals, which is the average time our model takes to exhume 1.2 m of rock through the rock–regolith interface given our assumed value of E. This depth is chosen as it broadly representing the recording timescale for cosmogenic radionuclides (Gosse and Philips, 2001). Exhumation is calculated as the difference between rock uplift (U_co+U_as) and bedrock surface uplift (S_B) over the recording time. We randomly choose 10 000 2 kyr samples from each 25 Myr run to produce a distribution of exhumation rates. We investigate the change in exhumation rates due to different proportions of coseismic and aseismic rock uplift and varying earthquake frequency and maximum earthquake magnitudes.

2.4 Cosmogenic radionuclide calculations

We also calculate the cosmogenic ¹⁰Be flux out of the model through time. For each time step, we add ¹⁰Be atoms to the system using published production rates and attenuation lengths to simulate the depth profile of cosmogenic concentrations (Balco et al., 2008; Braucher et al., 2011; Granger and Muzikar, 2001).

$$\begin{array}{}\text{(11)}& P_{\left(z\right)}=P_{\mathrm{0}}{e}^{-z(\mathit{\rho }/\mathrm{\Lambda })}\end{array}$$

The production rate (P) of ¹⁰Be decreases exponentially with depth below the topographic surface, z, from the production rate at the topographic surface (P₀) based upon the density of the bedrock (ρ) and the attenuation length (Λ). The production rates (atoms g⁻¹ yr⁻¹) and attenuation length (g cm⁻²) depend on the radiation being modelled. In our model, we simulate spallation (production rate of 5.784 atoms g⁻¹ yr⁻¹, attenuation length of 160 g cm⁻²), fast muons (production rate of 0.0418 atoms g⁻¹ yr⁻¹, attenuation length of 4320 g cm⁻²), and slow muons (production rate of 0.014 atoms g⁻¹ yr⁻¹, attenuation length of 1500 g cm⁻²), and combine them to give a total concentration at depth intervals set by the long-term erosion rate (Braucher et al., 2011; Granger and Muzikar, 2001). When an earthquake generates regolith, the top depth intervals are mixed, and the constant erosion rate is applied to remove regolith from the surface. After a spin-up time of 10 kyr, the model tracks concentration of ¹⁰Be in the sediment leaving the model. The spin-up time is the time required for the concentrations to reach a stable value which is perturbed by earthquakes. As erosion in the model is constant, and set to the long-term exhumation rate, any change in concentration represents the effect of stochastic EQTLs on exhumation.

3 Results and discussion

3.1 Coseismic landslide regolith production

Within our model, regolith generated by the largest earthquakes can remain on hillslopes for ∼1000 years (Fig. 2a). The average thickness of new regolith that is produced in an earthquake (expressed as volume per area; i.e. a depth) is a strong function of both the pre-existing depth of regolith prior to new failures, the magnitude of the earthquake, and stochastic differences in the volume of landslides for a given earthquake magnitude (Fig. 2b). The primary control on the total landslide volume produced by an earthquake magnitude is the strength of the shaking, with topography and rock strength as secondary factors (Marc et al., 2016a; Valagussa et al., 2019). If shaking can produce similar volumes of landsliding regardless of how much bedrock or regolith is on the hillslopes, landslide deposits in a mountain range with widespread regolith will contain less fresh bedrock, as regolith will make up a greater proportion of material mobilised by the earthquake. As regolith makes up a greater volume of the landslide deposits, less bedrock weathering occurs. This effect will be particularly powerful in areas where large earthquakes occur frequently in similar locations. However, the distribution of earthquake magnitudes exerts a stronger influence on total regolith production through time, as more frequent small earthquakes can only ever weather small depths of regolith from the bedrock (Fig. 2b). The decline of coseismic landslide regolith production (CLRP) with existing regolith thickness is reminiscent of soil production functions described for soil mantled landscapes (Heimsath et al., 1997), and by analogy, we term the non-linear relationship between regolith production rate and the average depth of weathering by landslides seen here a “coseismic landslide regolith production function” (CLRPF). However, unlike a “traditional” soil production function, two elements of stochasticity are inherent to a CLRPF. One reflects the role of shielding of the bedrock surface (S_B) from lowering when the regolith layer is thicker than the average depth of the generated landslides (Larsen et al., 2010), and that thickness is dependent on the past history of landsliding in the model. The other reflects the inherent randomness in the size and distribution of the landslides that occur in response to an earthquake of a given magnitude, i.e. within Eq. (8).

As expected, total seismic regolith production is dominated by the largest earthquakes, which produce the largest mean landslide volumes (Malamud et al., 2004b) (Fig. 2b). Summing through time, earthquakes produce 42 % of the total regolith generated by the model; M_w>7 earthquakes account for ∼65 % of the total earthquake-generated regolith and thus 27 % of the total regolith production. However, because smaller earthquakes occur often but produce little regolith, allowing the layer to thin, the time-averaged and spatially averaged regolith layer is predominantly thin – the modal thickness is only 0.02 m, and the mean is 0.03 m (Fig. 2b, inset). Thus, the model shows that although mountain belt regolith cover appears thin almost all the time, at times it can be thick enough to severely affect the short-term exhumation rates of the mountain range. The Longmen Shan region is primarily classified as a bedrock landscape with little storage, but significant volumes of sediment remain in the mountain range after the earthquake, in a similar way to that simulated by the model (Fan et al., 2018; Zhang et al., 2016). Variability in surface uplift through time (Fig. 2c) is affected by the pre-earthquake regolith thickness and therefore the sequence of earthquakes which occur before it. Where large earthquakes are closely spaced in time, pre-existing regolith can limit weathering of the bedrock surface, encouraging uplift of the topographic surface. In catchments close to active faults, the bedrock is likely to be heavily fractured and the shaking is more intense, producing larger landslides with greater densities (Marc et al., 2016a; Meunier et al., 2013; Valagussa et al., 2019). If the regolith is not fully removed from these catchments in between earthquakes, it is possible that the CLRPF may encourage greater surface uplift. In our model, the regolith produced by earthquakes is spread evenly across the landscape, which does not occur in reality. Even in the most impacted catchments in Wenchuan, landslide density is rarely above 10 % km⁻², suggesting that remobilisation of landslide regolith on the hillslope may be rare unless the regolith can remain in the catchment for multiple earthquake cycles (Dai et al., 2011; Marc et al., 2015; Parker et al., 2011, 2015). The remobilisation of previous coseismic landslide regolith is, likely to be a local effect mainly impacting catchments close to active fault belts. Ultimately, this interaction between surface uplift and regolith depth is controlled by (1) the time between earthquakes, (2) the magnitude of the previous earthquake, and (3) the rate of regolith removal. The closer together, in both time and space, large earthquakes occur and the slower regolith is removed from hillslopes, the greater the impact of the CLRPF on the surface uplift of a mountain range.

3.2 Regolith generation and volume budgets of earthquakes

Our model demonstrates that the contribution of earthquakes to the uplift and weathering budgets of mountains varies with earthquake magnitude, frequency, and the relative contribution of aseismic weathering and erosion, i.e. erosion or weathering not directly related to earthquakes such as rainfall-triggered landsliding or chemical weathering. When coseismic rock uplift is the dominant uplift mechanism, we can classify four distinct landscape response styles at different earthquake magnitudes (Fig. 3a and c, zones 1–4). On average, an earthquake of magnitude M_w<5.6 lowers the topographic surface (S_T). Here, the erosion of regolith (E) out of the model space is greater than the rock uplift produced by the earthquakes (zone 1 in Fig. 3a). For earthquakes with magnitudes $\mathrm{5.6}<M_{\mathrm{w}}<\mathrm{7.6}$, the bedrock surface (S_B) rises because the rock uplift rate is greater than the typical regolith generation rate (zones 2 and 3). The regolith thins because the change in the regolith thickness due to CLRP is less than the erosion rate (E) out of the model (zone 2). Conversely, the change in regolith thickness due to CLRP exceeds erosion rate in earthquakes with magnitudes M_w>6.4, so the regolith layer increases in thickness (zones 3 and 4). The largest earthquakes with M_w>7.6 lower the bedrock surface due to the large volumes of regolith produced (zone 4). However, much of the regolith is not removed before the next earthquake, resulting in a net topographic surface uplift, primarily due to thickening of the regolith. For a theoretical purely aseismic uplift scenario, where earthquakes produce landslides but do not create rock uplift (earthquakes produce only horizontal motion, this scenario is not realistic but acts as an extreme end member), earthquakes with M_w>6.5 would cause bedrock surface lowering, while smaller earthquakes would permit bedrock surface uplift due to low CLRP. Earthquakes with M_w>6.5 produce a thick layer of regolith which can persist until the next earthquake, limiting bedrock surface weathering and resulting in net uplift of the bedrock surface.

Figure 3 Interplay of changes to the modelled rock uplift rate, the topographic and bedrock surface uplift, and the resulting regolith thickness through time, classified according to earthquake magnitude. The total bedrock and topographic uplift produced by each earthquake magnitude through the model run is summed up and divided by the run time to produce a rate. Panel (a) represents a run with 100 % coseismic uplift, while panel (b) is purely aseismic. Each time a recorded surface intersects the horizontal axis, we separate the chart into a zone which is further described in the text and in panel (c).

Figure 4 Kernel density plots (bandwidth of 0.01) of exhumation and denudation rates in various scenarios. Dashed lines indicate the position of the mean ± the standard deviation of the distribution. Each curve is made up of 10 000 samples taken from their respective model runs. (a) Exhumation rates in different uplift regimes. (b) Denudation rates while varying the maximum earthquake magnitude in a run; the run with a maximum magnitude 5 has only earthquakes of a magnitude 5. (c) Denudation rates while varying the frequency of Wenchuan size earthquakes from every 500 to every 4000 years.

Figure 5 Variability of (a) topographic surface uplift and (b) the recorded concentration of cosmogenic nuclides leaving the orogen after a representative magnitude 8 earthquake within the model run. Red lines are the mean for the elapsed time since the earthquake, while the grey lines are the real-time concentrations.

3.3 Earthquakes and exhumation

We explore how earthquakes affect exhumation through direct calculation of exhumation of rock at the bedrock surface (S_B) relative to the topographic surface. There is very little (∼1 %) variability in erosion rates in model runs with maximum earthquake magnitudes of M_w<5.0. However, the introduction of stochastic weathering of many tens of cm of the bedrock surface by earthquakes with M_w>7 introduces variability in exhumation. When large earthquakes are present, exhumation rates have a standard deviation of 0.055–0.081 mm yr⁻¹ (9 %–14 % of the long-term exhumation rate) and a range of 0.77–0.89 mm yr⁻¹, with the lower figures reflecting a greater contribution of aseismic uplift (Fig. 4a and b). Decreasing the return times of Wenchuan-like earthquakes from 4000 to 500 years produces more variable distributions of exhumation rates, with the standard deviation of exhumation rate increasing from 0.044 mm yr⁻¹ (7 % of long-term average) to 0.076 mm yr⁻¹ (12 % of the long-term average) (Fig. 4c). Taken together, our model results suggest that up to 14 % of the variability in a sample of exhumation rates from a single geographical region could be associated with the time since the last large (M_w>7.0) earthquake. However, this variation may only be seen in exhumation or surface uplift records with recording times of less than 1000 years (Fig. 5a). Pre-Wenchuan earthquake measurements of cosmogenically derived erosion rates are between 40 % and 60 % lower than low-temperature thermochronometrically derived exhumation rates (Ouimet, 2010). Stochastic exhumation of low-concentration bedrock by EQTLs may explain some of that difference.

Cosmogenic radionuclides provide a record of potential earthquake-driven changes in exhumation because they have a relatively short averaging time that is close to the frequency of large earthquakes in many mountain belts. Our modelling results demonstrate the scale of stochastic variability in surface uplift and exhumation. We extend this analysis by simulating cosmogenic concentration in the model to estimate the potential impact of a large earthquake on both the cosmogenic concentration through time and the distribution of cosmogenic concentrations that are likely to be measured. We assume each earthquake thoroughly mixes the regolith down to its average weathering depth. The cosmogenic analysis (Fig. 5b) shows that immediately after a large earthquake, mixing of low-concentration bedrock material with higher concentration regolith lowers the concentration of radionuclides exiting the model. Regolith exiting the mountain range has a lower cosmogenic nuclide concentration for 200–300 years after the earthquake; after this period of low concentrations, there is a peak of concentration higher than the long-term average before a rapid return to the long-term average concentration (Fig. 5b). In the case of a representative magnitude M_w∼8 earthquake, the concentration falls initially by around one-third. However, the process of mixing also increases the concentration of nuclides close to the regolith–bedrock interface compared to the values before mixing, so that as the regolith is slowly eroded through time the lower half releases concentrations greater than the long-term average. Therefore, in landscapes with frequent M_w>7 earthquakes and regular long-term storage of regolith, it is possible to record more variable cosmogenic concentrations than might be expected, including positive as well as negative excursions from the long-term mean (Fig. 5b).

Figure 6 Reanalysis of detrital cosmogenic radionuclide derived denudation rates for mountain belts around the world compiled by Harel et al. (2016), in the context of peak ground acceleration and tectonic environment. (a) A box plot indicating the median (central orange line), quartiles (end of box) and the range (“the whiskers”) of denudation rates in the analysed localities ordered by their average seismicity (in brackets), defined as the maximum peak ground acceleration (PGA) of a 475-year return period. Points indicate values outside the range, ±1.5 times the interquartile range. (b) Standard deviation (SD) of denudation rates for each mountain belt against seismicity represented by the 475-year return PGA (m s⁻²). (c) Standard deviation of denudation rates for each mountain belt compared to the standard deviation of steepness indexes. The areas of highest variability are found in steep, tectonically active mountain ranges.

These modelling results provide testable predictions of the exhumation-related changes to cosmogenic concentration caused by large earthquakes. Hence, we look to examine whether the predicted variability might represent some of the variability associated with cosmogenic erosion rates in seismically active areas using a global dataset compiled by Harel et al. (2016). Harel et al. (2016) collated detrital cosmogenic ¹⁰Be concentrations for 59 geographical areas, separated into areas of similar climates, and recalculated the erosion rates using consistent production rate and shielding corrections. We limit our sampling to those sites with >18 measurements and basins larger than 10⁵ m² to limit sampling bias (Dingle et al., 2018; Niemi et al., 2005). We compare the probability density distribution of erosion rates from within those geographic regions to seismicity, as represented by the 475-year return peak ground acceleration (PGA) derived from a global seismic hazard map (Giardini et al., 1999; Harel et al., 2016) (Fig. 6a). Due to the size of the geographical areas, there may be multiple seismic hazard levels recorded; we simply use the mean value to classify the area. The use of a single number to characterise a large area can underestimate the potential PGA. While a single number may not accurately describe the entire area, it does allow us to compare the variability of denudation rates with seismic hazard. We crudely classify the regions as dominantly coseismic or dominantly aseismic: regions with thrust faults (identified from a literature review) and erosion rates greater than the median are deemed coseismic, while slowly eroding regions with no thrust faults are aseismic. “Mixed” regions are those that do not fall under either of these classifications. Coseismic landscapes, as expected, have higher average cosmogenically determined erosion rates with higher standard deviations (Kruskal–Wallis H-test statistic of 14.14; p value of 0.00017) (Fig. 6a and b). Relief is a major control on erosion rates, with steeper catchments having higher erosion rates than shallower ones (Montgomery and Brandon, 2002). The most seismically active mountain ranges are also among the most varied in relief as they have some of the steepest catchments in the world. Therefore, we need to test whether variability in erosion rates is more closely related to variation in catchment steepness or the seismicity of the area. Within the database compiled by Harel et al. (2016), they include a normalised channel steepness index which we use to compare the impact of catchment steepness on erosion rate. The channel steepness index equation is defined by

$$\begin{array}{}\text{(12)}& M_x={\left({\displaystyle \frac{E}{K{A}_{\mathrm{0}}^m}}\right)}^{\frac{\mathrm{1}}{n}},\end{array}$$

where M_x is the steepness index, E is the erosion rate, K is the erodibility coefficient, A₀ is a reference area of 1 m², and m and n are empirical constants. The index is normalised by assuming fixed values for m and n (Harel et al., 2016). We would expect that in areas with highly variable steepness indexes the erosion rates are also more variable. We find that while areas with higher seismicity have more variable erosion rates, the variation in erosion rates correlates much more closely with the variation in steepness index (Fig. 6b and c). Steeper basins in tectonically active mountain ranges are more susceptible to coseismic landsliding (G. Li et al., 2017; Marc et al., 2016a) and thus will have more variable denudation rates through time, depending on the residence time of the landslide regolith and the frequency of earthquakes, than shallower basins. Landscapes with more variation in basin relief could enhance the temporal perturbations in denudation rates produced by earthquakes, but the contribution of tectonics to long-term variation is difficult to isolate.

We also explore the averaging time required to reach the mean exhumation rate in model runs with different contributions from seismic and aseismic uplift. After a large earthquake that produces tens of centimetres of instantaneous weathering of S_B, exhumation rates measured with different averaging times converge to the long-term mean rate within hundreds to thousands of years (Fig. 5a). Bedrock surface exhumation rates are impacted by both surface uplift and lowering rates, so as a result the timescale of the perturbation is impacted by the dominant form of uplift in the mountain range. As a result, the more dominant coseismic uplift is in a landscape, the longer the recording time needs to be before a reliable exhumation record can be made. Landscapes with more frequent earthquakes have more variable exhumation rates which require longer averaging times to achieve accurate measurements of the long-term average exhumation rate. Regardless of the frequency of earthquakes in a mountain range, the events of the greatest magnitude remain uncommon while being the dominant contributors to weathering. Hence, the relationship between exhumation rate and averaging time is consistent with the Sadler effect that has been described for sedimentary systems (Schumer and Jerolmack, 2009). Unlike sedimentary systems, where it is possible to measure sedimentation rates from timescales of seconds to millions of years, there are few measures of exhumation that we can use and many of these have long averaging times. Thermochronometric methods average across timescales of 10⁵–10⁷ years, much longer than the recurrence times of individual earthquakes. There is a possibility that cosmogenic radionuclide analysis records individual earthquakes, where earthquake-triggered landslides weather bedrock that has a low cosmogenic concentration (Wang et al., 2017; West et al., 2014), although enhanced erosion during and immediately after an earthquake complicates the cosmogenic signal in practice.

Despite representing close to half of the weathering flux of mountain belts, stochastic earthquakes still remain substantively missing from our models of the development of mountain belts. The modelling here demonstrates that stochastic uplift and exhumation by large earthquakes is likely to be averaged away across the timescale of most thermochronometers, with the variability in uplift and exhumation representing around 15 % of the average exhumation rate. Even when a large earthquake occurs at the edge of a mountain belt, as has occurred in the Wenchuan region, the variable exhumation signal is further shredded by sediment transport by floods and debris flows, such that even sinks that are within 40 km of the epicentre show limited evidence for large earthquakes (Zhang et al., 2019). This result along with our model demonstrates the importance of understanding the processes by which landslide sediment are mobilised out of catchments and the time taken for these processes. Without improved understanding of the cascading nature of sediment transport from catchments, it is unlikely we will be able to identify earthquakes other than within smaller basins or sinks immediately adjacent to the epicentre (Howarth et al., 2012). Large basins (greater than 10 000 km²) have been shown to be large enough to average out the perturbations in cosmogenic radionuclide concentrations caused by large landsliding events (Dingle et al., 2018; Marc et al., 2019). While the largest basins are able to offer reliable estimates of long-term erosion rates, the detrital cosmogenic nuclide concentrations from smaller basins will be more affected by bedrock landsliding caused by earthquakes. Therefore, smaller basins could be suitable targets to recognise variations in erosion rates due to earthquakes. Our model suggests that the impact of a large earthquake is not necessarily big enough to perturb the denudation rate of an orogen for the whole of a cosmogenic nuclide recording time. The combination of averaging times, shredding, and the relative contributions of large earthquakes to long-term exhumation rates helps us to understand the lack of clear signatures for single earthquakes in sedimentary or exhumation records.

4 Conclusions and implications

Our simulations show that the regolith generated by large earthquakes can reduce the rate of weathering and exhumation of rock due to its potentially long residence time on hillslopes. Reducing exhumation rates also increases the uplift rate of the bedrock surface, but these effects are small when compared to the role of the magnitude and frequency of earthquakes. These results demonstrate that background tectonic processes and rates are the dominant control on surface uplift, while the more important role for large earthquakes is their control on weathering. Small earthquakes contribute very little to both uplift and weathering resulting in below-average rock exhumation being recorded if a large earthquake does not occur during the averaging time of the exhumation record. While large earthquakes produce higher-than-average rock exhumation rates, the slow removal of regolith from the orogen reduces the magnitude and timescale of the signal. The relatively long timescales of exhumation records prevent the recording of orogen-scale variation in exhumation due to a single earthquake. A better understanding of the controls on bedrock weathering by earthquake-triggered landslides is required to identify signals of earthquakes in sedimentary records. Higher-resolution exhumation records and the growing recognition of the complex nature of exporting landslide sediment from mountainous catchments will help to explore this problem.