Pulsed carbon export from mountains by earthquake-triggered landslides explored in a reduced-complexity model

In mountain ranges, earthquakes can trigger widespread landsliding and mobilize large amounts of organic carbon by eroding soil and vegetation from hillslopes. Following a major earthquake, the landslide-mobilized organic carbon can be exported from river catchments by physical sediment transport processes or stored within the landscape where it may be degraded by heterotrophic respiration. The competition between these physical and biogeochemical processes governs a net transfer of carbon between the atmosphere and sedimentary organic matter, yet their relative importance following a large landslide-triggering earthquake remains poorly constrained. Here, we propose a model framework to quantify the post-seismic redistribution of soil-derived organic carbon. The approach combines predictions based on empirical observations of co-seismic sediment mobilization with a description of the physical and biogeochemical processes involved after an earthquake. Earthquake-triggered landslide populations are generated by randomly sampling a landslide area distribution, a proportion of which is initially connected to the fluvial network. Initially disconnected landslide deposits are transported downslope and connected to rivers at a constant velocity in the post-seismic period. Disconnected landslide deposits lose organic carbon by heterotrophic oxidation, while connected deposits lose organic carbon synchronously by both oxidation and river export. The modeling approach is numerically efficient and allows us to explore a large range of parameter values that exert a control on the fate of organic carbon in the upland erosional system. We explore the role of the climatic context (in terms of mean annual runoff and runoff variability) and rates of organic matter degradation using single pool and multi-pool models. Our results highlight the fact that the redistribution of organic carbon is strongly controlled by the annual runoff and the extent of landslide connection, but less so by the choice of organic matter degradation model. In the context of mountain ranges typical of the southwestern Pacific region, we find that model configurations allow more than 90 % of the landslide-mobilized carbon to be exported from mountain catchments. A simulation of earthquake cycles suggests efficient transfer of organic carbon out of a mountain range during the first decade of the post-seismic period. Pulsed erosion of organic matter by earthquake-triggered landslides is therefore an effective process to promote carbon sequestration in sedimentary deposits over thousands of years.

1 Introduction

Erosion of carbon from the terrestrial biosphere results in an important transfer of carbon from land to oceans by rivers (Galy et al., 2015; Stallard, 1998). Some of the eroded organic matter escapes degradation and can be buried in sediments (Burdige, 2005; Galy et al., 2007), contributing to a long-term sequestration of atmospheric carbon dioxide (CO₂) (Hilton and West, 2020). To constrain how the export of organic carbon (OC) from land may vary over space and time and force the carbon cycle, we need to understand the factors that drive its erosion from the landscape and transfer by rivers (Berhe et al., 2018; Hilton, 2017). In mountain ranges located at convergent plate boundaries, large earthquakes (M_w>6) can trigger widespread landsliding (Keefer, 1984; Tanyaş et al., 2017), which harvests OC by mobilizing soil and vegetation from hillslopes (Allen et al., 1999; Garwood et al., 1979). Earthquakes and other events which trigger widespread disturbance of forests, such as storms, have been considered a CO₂ source due to degradation of the eroded organic matter by heterotrophic respiration (Chambers et al., 2007; Zeng et al., 2009). However, recent work has highlighted the fact that landslide-mobilized OC can be exported from mountain catchments by rivers (Hilton et al., 2008a; Wang et al., 2016) and that this export can be sustained over thousands of years, as recorded in lake deposits (Howarth et al., 2012; Frith et al., 2018). When coupled to the regrowth of the forest and soil carbon stocks at the site of erosion, this drives net CO₂ drawdown from the atmosphere (Berhe et al., 2007; Stallard, 1998). Despite this recognition, it remains a challenge to constrain how variability in climate (e.g., annual runoff) and biogeochemical processes (such as respiration) controls the fate of landslide-mobilized OC and the timescales over which their impacts are felt (Ramos Scharrón et al., 2012; Restrepo et al., 2009). We require better constraints on the role of climatic and biogeochemical processes over timeframes that span multiple seismic cycles. To achieve this, we require a theoretical framework to explore the competition between physical and biogeochemical processes from the site of erosion to the catchment outlet.

The export of OC from mountain river catchments is linked to that of fine clastic sediment (Galy et al., 2015; Hilton et al., 2012). Following widespread earthquake-triggered landsliding, fine sediment can be evacuated as fluvial suspended load within a decade in wet, subtropical settings such as Taiwan (Dadson et al., 2004; Hovius et al., 2011) or over several decades in more arid locations (Tolorza et al., 2019; Wang et al., 2015). The total landslide-derived sediment volume, including coarse material which can be mobilized by debris flows and transported as bed load, may take decades to centuries to export (Croissant et al., 2017, 2019b; Fan et al., 2018; Yanites et al., 2010). This variability in the temporal dynamics of post-seismic sediment transfer highlights the role of factors such as the connectivity of sediment to the fluvial network and the rate of transport during the post-seismic phase (Wang et al., 2015). The timescales of sediment routing are key to understanding the fate of the landslide-eroded OC. Any delay in transport provides time for OC to be degraded in the landscape (Berhe et al., 2007; Ramos Scharrón et al., 2012). This could happen, for example, if OC is exposed to an O₂-rich atmosphere or surface waters while stored in landslide deposits or colluvium, or it could happen during transient storage of eroded materials in fluvial bars.

In this study, we introduce a reduced-complexity modeling approach to quantify the fate of OC eroded from soils after a widespread landslide-triggering event, such as a major earthquake. This is also relevant to better understand how organic matter may provide a novel tool of erosion provenance in sedimentary records (Wang et al., 2020). Our aim is to explore the climatic, tectonic and biogeochemical controls on the fate of OC over decadal to centennial timescales. To do this, we describe the evolution of OC in landslide deposits by two principal processes: (i) the physical transport of OC as suspended load by fluvial export, which can act as a potential CO₂ sink; and (ii) the oxidation of OC by heterotrophic respiration and subsequent CO₂ release. We explore several approaches to quantify organic matter degradation using different biochemical models described in the literature. We conceptualize our model based on the observed processes and boundary conditions operating in the western Southern Alps of New Zealand. A key concept is that of landslide connection to the fluvial network and the potential evolution of that connection through time following a large earthquake. We consider timescales from the occurrence of a single landslide to several earthquake cycles and explore the climatic boundary conditions (mean annual runoff and runoff variability) which impact the fate of landslide-mobilized OC.

2 Methods

We use a reduced-complexity model approach to quantify the fate of landslide-mobilized OC considering the case of earthquake-triggered landslides. The model Quakos is described at length by Croissant et al. (2019b). Herein, we detail the components involved in this study and present the biochemical models used to explore the fate of eroded OC. The model considers the erosion and transport of fine-grained, soil-derived OC with a diameter that is nominally ∼1 mm or finer. This material contributes to the suspended sediment load of mountain rivers (Clark et al., 2017; Hatten et al., 2012; Hilton et al., 2008a) and can be exported to larger river systems (Bouchez et al., 2014) and sedimentary deposits (Galy et al., 2007). At present we cannot simulate the movement of coarse OC supplied by erosion of living biomass (and any associated coarse woody debris), mainly due to the lack of a robust transport law (Wohl, 2011). Recent work suggests that this coarser woody material could be a large component of erosional OC budgets (Mohr et al., 2017; Turowski et al., 2016; West et al., 2011). Indeed, it may be transported long distances and preserved in marine sediments (Lee et al., 2019). We also do not consider rock-derived OC in the model framework. While the oxidation of rock OC may be substantial in the landscape (Hemingway et al., 2018; Horan et al., 2017), the degradation rates of this material remain poorly constrained (Chang and Berner, 1999) and here we do not attempt to model its contribution to the net CO₂ budget of erosion (Hilton and West, 2020). As such, the model results described herein do not capture the full impact of erosion on the fate of OC, but they provide a framework that could be adapted in the future to include these additional carbon transfers.

2.1 Study area

A conceptual framework is provided by a hypothetical earthquake on the Alpine Fault in the western Southern Alps, New Zealand (Fig. 1). The topography of the west coast of the Southern Alps is obtained from the SRTM3 digital elevation model (DEM). This location is selected due to prior work documenting the role of landsliding in sediment and OC mobilization (Hilton et al., 2008a; Hovius et al., 1997; Korup et al., 2010) and the role of large earthquakes in long-term sediment and OC records in lakes (Frith et al., 2018; Howarth et al., 2012, 2014). The Southern Alps formed as a result of oblique plate convergence between the Pacific and Australian plates at a rate of 39.7 mm yr⁻¹ (DeMets et al., 2010). The Alpine Fault accommodates up to 80 % of the plate convergence (Norris and Cooper, 2007) in earthquakes of magnitude >7.5, with the last major earthquake occurring in 1717 CE and an average recurrence time of 263±68 years as determined from paleoseismic studies (Howarth et al., 2016, 2018). The Southern Alps form a barrier to westerly winds that leads to high precipitation rates of up to 13 m yr⁻¹ along the west coast (Tait and Zheng, 2007). Landscapes in the Southern Alps are characterized by steep hillslopes with modal gradients of approximately 35^∘ (Korup et al., 2010), and hillslope erosion is dominated by landsliding in the current aseismic period (Hovius et al., 1997).

Figure 1 (a) An example of a modeled distribution of earthquake-triggered landslides from the model Quakos. The view is centered on the Whataroa catchment on the west coast of the South Island of New Zealand, and the dot size is proportional to the landslide volume. Red lines outline some of the west coast catchments that are accounted for in the simulations. The white dots are landslides generated by Quakos but not accounted for in our simulations. (b) Illustration of the two possible connectivity statuses and their connection velocity (u_con), distance from channel d, and time for connection to the channel (t_con) for each individual landslide of volume V_ls.

The high rainfall sustains temperate rainforest on steep slopes at elevations ≤800 m, while shrubs, herbs and grassland persist up to and above the regional snowline at ∼1250 m (Reif and Allen, 1988). The carbon stocks of aboveground biomass in the western Southern Alps have been estimated at 17 500±5500 tC km⁻² based on local forest plot inventories (Hilton et al., 2011). The OC stocks of soils in the Southern Alps have been recently estimated at $\sim \mathrm{13}\mathrm{000}\pm \mathrm{4300}$ tC km⁻² (Harvey, 2019). This estimate was derived from measured OC concentrations in soil profiles collected from 52 sites in the mountain range. The OC stocks are estimated to the point of refusal of the soil auger profiles (average thickness 0.42±0.09 m) and comprise soil OC in the upper organic-rich horizons (mean thickness 0.10±0.01 m) and in the mineral-dominated B and C horizons beneath. The average value is slightly lower than previous estimates, which were inferred from sparse data (Frith et al., 2018; Hilton et al., 2011; Tonkin and Basher, 2001).

2.2 Generation of earthquake-triggered landslide clusters

Populations of landslides triggered by large earthquakes are generated in the Quakos model, with an empirical relationship between the spatial patterns of the simulated peak ground acceleration (PGA) and landslide density (see Croissant et al., 2019b, for further details). The resulting landslide density map allows for the quantification of the total area of landsliding (A_ls, tot) and determines the spatial extent of mass wasting. The co-seismic landslides are described by several components: their geometrical properties of their scar area (A_ls, m²) and deposit volume (V_ls, m³), their spatial distribution within the landscape, and their connectivity to the drainage network. Based on previous empirical studies (Malamud et al., 2004; Tanyaş et al., 2017), the full distribution of landslide areas can be described by a three-parameter inverse gamma function of the form

$$\begin{array}{}\text{(1)}& \mathrm{PDF}\left(A_{\mathrm{ls}}\right)={\displaystyle \frac{\mathrm{1}}{a\mathrm{\Gamma }\left(a\right)}}{\left[{\displaystyle \frac{a}{A_{\mathrm{ls}}-s}}\right]}^{\mathit{\rho }+\mathrm{1}}\exp \left[{\displaystyle \frac{-a}{A_{\mathrm{ls}}-s}}\right],\end{array}$$

where a is the position of the probability density function (PDF) maximum, s controls the rollover for small landslides, ρ is a positive exponent controlling the slope of the tail and Γ is the gamma function. The parameterization of Eq. (1) has been done using values obtained for a west coast landslide inventory (Hovius et al., 1997). In Quakos, a landslide population is generated by sampling the PDF until the sum of the landslide areas reaches the value of A_ls, tot. The landslides are spatially distributed within catchments following a landslide density map that depends on the PGA patterns triggered by the earthquake (Meunier et al., 2008). Landslide area is converted to volume with an empirical law of the form

$$\begin{array}{}\text{(2)}& V_{\mathrm{ls}}=\mathit{\alpha }{A}_{\mathrm{ls}}^{\mathit{\gamma }},\end{array}$$

where the empirical constants α and γ are set to 0.05 and 1.5 (Hovius et al., 1997; Larsen et al., 2010) to reflect the deep-seated bedrock landslides that are likely to be the dominant mechanism of landslide generation during large earthquakes in New Zealand (Larsen et al., 2010). V_ls is a measure of the total volume of sediment mobilized by landslides. However, in the subsequent modeling we consider only the fine-grained portion of this landslide volume that can be transported as suspended load (V_ls, fine), here set to 50 % of V_ls, as an average value considering empirical estimations varying between 10 % and 90 % (Dadson et al., 2003). The landslide volume is distributed in a single model cell.

Each landslide is assigned a probability that it is initially connected to the fluvial network (χ), which is a function of the drainage area of landslide pixels following an empirical relationship of the form (e.g., 2008, Wenchuan – Li et al., 2016; 2015, Gorkha – Roback et al., 2018):

$$\begin{array}{}\text{(3)}& \mathit{\chi }=\mathrm{c}{A}_{\mathrm{ls}}^{\mathit{\mu }},\end{array}$$

where c and μ are empirical constants. Equation (3) is used to assign the initial connectivity status of all landslides with an area smaller than 10⁶ m², while larger landslides are assumed to be initially connected (Croissant et al., 2019b). Using c=0.86 and μ=0.345, only 12 % of the landslides are initially connected to the drainage network (Fig. 1). However, because the largest landslides are more likely to be initially connected, between 40 % and 60 % of the landslide-derived sediment volume is initially available. We assume that sediment derived from initially connected landslides is immediately available for fluvial export toward the catchment outlet.

To account for processes which mobilize sediment from landslide deposits that are not initially connected to the river network (including soil creep, overland flow, shallow landsliding and debris flows), which are known to act (Dietrich et al., 2003; Roering et al., 2001) but are otherwise difficult to parametrize, our approach allows landslide deposits on hillslopes to progressively reach the river network at a rate defined as a constant “connection velocity” u_con (m yr⁻¹) (Croissant et al., 2019). The path taken by the landslides to reach the closest river is computed using a steepest descent algorithm, and their time of connection (t_con) is determined by dividing the distance to the river network by u_con. The river network is extracted from the DEM using a single flow algorithm and by setting a critical drainage area of 0.5 km².

To quantify the erosion of soil OC, the area of each individual landslide scars is combined with the local soil organic carbon stock (C_stock, tC m⁻²) determined in this study to quantify the total mass of eroded soil OC (tC):

$$\begin{array}{}\text{(4)}& M_{\mathrm{oc},\mathrm{tot}}=\sum _{\mathrm{0}}^N{A}_{\mathrm{ls},n}{\mathrm{C}}_{\mathrm{stock},n},\end{array}$$

where N is the total number of earthquake-triggered landslides. Simulations account for the uncertainties related to the OC content of soils. As outlined at the opening of Sect. 2, the transport of the aboveground biomass is not considered here as we assume it is mobilized as large woody debris, whose transport dynamics remain poorly constrained (Wohl, 2011).

2.3 Along-stream river transport capacity

In the model framework, once a landslide deposit reaches the fluvial network, fine sediment and OC are transported at a rate set by the local transport capacity of the river (Wang et al., 2015). Here, we only consider the evacuation of the fine sediment and OC that are transported as suspended load, with a nominal grain size <1 mm (see Sect. 2). First, we consider the event-based sediment discharge of the suspended load (Q_s, m³ s⁻¹) that is commonly empirically described as a power-law function of water discharge (Q, m³ s⁻¹):

$$\begin{array}{}\text{(5)}& Q_{\mathrm{s}}=k_{\mathrm{s}}{Q}^{\mathit{\epsilon }},\end{array}$$

where k_s is a constant (s m⁻³) and ε is a positive exponent. Here, the water discharge is a linear function of the local drainage area and surface runoff. The constants k_s and ε can be determined from measurements across a range of discharges. For our study location, we use gauging data from the Hokitika, Haast, Whataroa and Poerua rivers, and we find values of k_s=300 and ε=2. To encompass the contribution of all discharge events to sediment transport, the annual long-term sediment flux is computed as the integral of the convolution between sediment flux (Eq. 5) and the probability density function (PDF) of daily discharge events (PDF(Q)) (DiBiase and Whipple, 2011; Lague, 2014; Lague et al., 2005):

$$\begin{array}{}\text{(6)}& \stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}=\underset{\mathrm{0}}{\overset{Q_{\max }}{\int }}{Q}_{\mathrm{s}}\mathrm{PDF}\left(Q\right),\end{array}$$

where Q_max is the maximum discharge in the range. The PDF of daily discharges is described by an inverse gamma function of the form

$$\begin{array}{}\text{(7)}& \mathrm{PDF}\left(Q\right)={\displaystyle \frac{k^{k+\mathrm{1}}}{\mathrm{\Gamma }\left(k+\mathrm{1}\right)}}\exp \left({\displaystyle \frac{-k}{Q/\stackrel{\mathrm{‾}}{Q}}}\right){\left(Q/\stackrel{\mathrm{‾}}{Q}\right)}^{-\left(\mathrm{2}+k\right)},\end{array}$$

where k is a positive constant expressing the variability of the discharge events. A low value of k corresponds to higher discharge variability.

For each landslide, the evolution of the mass of fine sediment remaining in the landslide deposit while fluvial transport removes mass is then described by

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{ls}}^{\mathrm{riv}}}{\mathrm{d}t}}=-\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}.\end{array}$$

Equation (8) is solved to get the temporal evolution of the fine sediment mass remaining in the deposit as

$$\begin{array}{}\text{(9)}& M_{\mathrm{ls},t}^{\mathrm{riv}}=M_{\mathrm{ls},\mathrm{0}}-\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}t=\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}\left(t_{\mathrm{0}}-t\right),\end{array}$$

where M_ls, 0 is the initial mass of fine sediment mobilized by landsliding and $t_{\mathrm{0}}=\frac{M_{\mathrm{ls},\mathrm{0}}}{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}$ the time necessary to export the entirety of the fine sediment.

The mass of fine sediment exported by the river at any point in time is given by

$$\begin{array}{}\text{(10)}& m_{\mathrm{ls},t}^{\mathrm{riv}}={\stackrel{\mathrm{‾}}{Q}}_{\mathrm{s}}\mathrm{d}t,\end{array}$$

with dt being the time step over which the computation is completed.

2.4 Sediment transport model assumptions

We use a reduced-complexity approach to describe the physical processes in play, which allows us to explore the large spatial and temporal scales over which earthquake-triggered landslides impact landscapes. However, this means that the transport of fine sediment and particulate OC is subject to several assumptions. First, by considering only the transport of fine, soil-derived particles as suspended load, each landslide volume is partitioned between fine grains (that can be transported in river suspension) and coarse grains (which are transported by saltation and other bed-load processes). Second, sediment transport is treated as a detachment-limited regime, and entrained particles are never redeposited within the catchment. This assumption is likely to be reasonable for the Southern Alps, which are characterized by steep and short rivers that promote efficient suspended sediment transport (Korup et al., 2010). We also assume that the fine sediment is exported at fluvial transport capacity (Eq. 8). We do not consider in situ stabilization of landslide materials, regrowth of vegetation on those materials or subsequent remobilization by landsliding. Therefore, all landslide-derived fine grains that were initially mobilized will be eventually exported by rivers at the timescale controlled by the local transport capacity. In addition, the climatic context (mean annual runoff and runoff variability) that controls the fluvial transport capacity is assumed to be constant over the duration of the simulations. Finally, any other event that might be triggering landslides during the post-seismic phase (rainstorms, earthquakes of lower magnitude, etc.) are not considered.

2.5 Evolution of OC after landslide mobilization

Our model postulates that the fate of landslide-mobilized OC is controlled by its physical export by rivers versus its biogeochemical oxidation that can release CO₂ back to the atmosphere. The long-term sink of CO₂ by erosion comes about because a fraction of the mobilized OC is stored in sediments, while the vegetation and soil OC are renewed at the site of erosion (Berhe et al., 2007; Stallard, 1998). If eroded materials remain on hillslopes, the OC contained in the landslide deposits can be subjected to heterotrophic respiration and OC oxidation, either by soil fauna present in the original forest soil and translocated to the landslide deposit or by colonization by new fauna communities. When landslides are connected to the river network, the model framework considers the OC remaining in the deposit to experience oxidation and physical transport simultaneously. OC oxidation can be generally described by either discrete organic matter pools with a specific rate of degradation or by continuum models that seek to explore a range of oxidation rates (e.g., Arndt et al., 2013; Manzoni et al., 2009). In this section, we describe the different models of organic matter degradation used in our approach.

2.5.1 Oxidation models

There are no models that specifically describe the degradation of organic carbon within landslide deposits, but there are likely to be close analogies with those that are applied to soils and/or marine sediment (e.g., Arndt et al., 2013; Manzoni et al., 2009). The simplest model of the degradation of organic matter represents the carbon stock as one compartment, the so-called “single pool” model (Berner, 1964). This pool is described by a mass of OC contained in the soil layer mobilized by a landslide that loses carbon (during respiration the OC loss would be to CO₂) via first-order kinetics at a rate described by a single oxidation constant, k_ox (% yr⁻¹). For a situation with no further OC inputs, such as a single landslide, the OC loss of the deposit over time can be summarized as

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{oc}}^{\mathrm{ox}}}{\mathrm{d}t}}=-k_{\mathrm{ox}}{M}_{\mathrm{oc}}^{\mathrm{ox}},\end{array}$$

where $M_{\mathrm{oc}}^{\mathrm{ox}}$ is the mass of OC remaining in the landslide deposit.

An alternative is to describe the mass of OC by multiple compartments with different degradation rates and oxidation constants, referred to as the “multi-pool model”. Here, the pools are treated in parallel, and the user defines the distribution of the initial OC mass between different pools. Each pool is described by an individual oxidation constant, which could conceptually relate to the reactivity of different organic compounds (e.g., Minderman, 1968) and/or the association of organic matter with different mineral phases (e.g., Hemingway et al., 2019; Mayer, 1994). This approach attempts to capture the observation that as organic matter degrades and ages, its apparent reactivity and rate of degradation decline (e.g., Middelburg, 1989), with the first pool containing the youngest most reactive OC and the second and third pools capturing organic matter that has been retained for longer periods of time.

In our study, a single pool model may be appropriate if the residence times of landslide deposits in the catchment are short (annual timescales) so that it is only necessary to constrain the fate of the most reactive OC (Trumbore, 2000). This simple approach was applied by Wang et al. (2016) to draw conclusions on the fate of OC eroded by landslides during the Wenchuan earthquake. However, we know that at least some landslide deposits are likely to reside in the landscape for decades or longer (Clark et al., 2016; Fan et al., 2018), so a multi-pool approach may be more appropriate. In addition, landslides can mobilize the entire soil profile, including OC with different reactivity as a function of depth. For instance, the organic-rich O and A layers are highly reactive, whereas the mineral soil may have lower turnover rates (Tate et al., 1995) due to different organic compound reactivity and/or mineral–OC association (Hemingway et al., 2019). In addition, a multi-pool model may also better capture processes operating within a landslide deposit. Burial of organic matter within the sediment pile, coupled with potential waterlogging as water flow paths focus runoff from landslide scars into deposits (Emberson et al., 2016a; Lo et al., 2012), could mean that OC on the surface of deposits is more quickly degraded than OC at depth.

To constrain appropriate values of k_ox, we note that in New Zealand temperate forests the surface soil OC is estimated to turnover over every ∼20 years, which in a single pool model would correspond to k_ox∼5 % yr⁻¹ and k_ox values between 1.2 % yr⁻¹ and 0.5 % yr⁻¹ for the total soil carbon (Tate et al., 1995). Thus, for the single pool model, we explore k_ox values from 0.1 % yr⁻¹ to 2 % yr⁻¹ to reflect the total soil layer. For the multi-pool model, we assign three pools associated with k_ox of 2 % yr⁻¹, 1 % yr⁻¹ and 0.5 % yr⁻¹ to reflect the local empirical constraint on soil OC degradation (Tate et al., 1995).

2.5.2 Quantification of OC redistribution

In our modeling approach, the temporal evolution of the initial mass of mobilized OC is controlled by two processes: biochemical degradation of the OC and fluvial export. We first present the equations controlling these two processes when they operate independently and then present when they operate synchronously.

The OC biochemical degradation is described by a first-order oxidation reaction law obtained by solving Eq. (11) (Blair and Aller, 2012; Trumbore, 2000; Wang et al., 2016):

$$\begin{array}{}\text{(12)}& M_{\mathrm{oc},t}^{\mathrm{ox}}=M_{\mathrm{oc},\mathrm{0}}{e}^{-k_{\mathrm{ox}}t},\end{array}$$

where $M_{\mathrm{oc},t}^{\mathrm{ox}}$ is the mass of OC remaining in the landslide deposit considering only loss of C by biogeochemical degradation (assumed production and release of CO₂), M_oc, 0 is the initial mass of OC mobilized by the landslide, and k_ox is the first-order reaction kinetic constant for OC oxidation (Fig. 2a).

Figure 2 The fate of landslide-mobilized OC considered in terms of the proportion of the mass of OC remaining within a single landslide deposit as a function of time. (a) An example of a connected landslide for which the OC decrease is due to physical export and oxidation acting simultaneously (blue lines). The impact of the evacuation time is shown here by different values of t₀ (in years). The case of OC oxidation alone, with no physical export, is represented by the red line. (b) An example of a landslide deposit that is initially disconnected from the river network and that connects after a connection time, t_con=30 years, and is exported over a duration of t₀=50 years. The dotted ride line represents the OC oxidation alone.

The mass of OC being oxidized at any point in time is described as

$$\begin{array}{}\text{(13)}& m_{\mathrm{oc},t}^{\mathrm{ox}}=k_{\mathrm{ox}}{M}_{\mathrm{oc},t}^{\mathrm{ox}}=k_{\mathrm{ox}}{M}_{\mathrm{oc},\mathrm{0}}{e}^{-k_{\mathrm{ox}}t}\mathrm{d}t.\end{array}$$

The mobilized OC can also be physically exported from the deposit by fluvial transport. We assume that the fine particulate OC has a transport behavior which is similar to that of clastic suspended sediments, as demonstrated in a number of studies of turbulent mountain rivers (Clark et al., 2017; Hatten et al., 2012; Hilton, 2017; Smith et al., 2013). Therefore, we consider the fluvial transport of the particulate OC mobilized by landslides and of the fine sediment to be proportional. In the case in which oxidation is not an active process, the export of OC by rivers would therefore be described as

$$\begin{array}{}\text{(14)}& m_{\mathrm{oc},t}^{\mathrm{riv}}={\stackrel{\mathrm{‾}}{Q}}_{\mathrm{s}}{\displaystyle \frac{M_{\mathrm{oc},t}}{M_{\mathrm{ls},t}}}={\displaystyle \frac{M_{\mathrm{oc},t}^{\mathrm{riv}}}{t_{\mathrm{0}}-t}}.\end{array}$$

The subsequent evolution of the OC remaining within the landslide deposit, considering only loss of C by erosion and river export of OC, is described as

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{oc}}^{\mathrm{riv}}}{\mathrm{d}t}}=-m_{\mathrm{oc},t}^{\mathrm{riv}},\end{array}$$

where M_riv is the mass of OC exported by rivers and M_ls is the mass of fine sediment mobilized by the landslide. Solving Eq. (15), we determine the temporal evolution of the mass of OC exported by rivers to be

$$\begin{array}{}\text{(16)}& M_{\mathrm{oc},t}^{\mathrm{riv}}=M_{\mathrm{oc},\mathrm{0}}\left(\mathrm{1}-{\displaystyle \frac{t}{t_{\mathrm{0}}}}\right).\end{array}$$

When we combine both C loss by biogeochemical degradation (Eq. 12) and the OC loss by fluvial transport (Eq. 16), we can describe the coupled loss of the landslide-mobilized OC from a deposit (Fig. 2a):

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{oc}}^{\mathrm{riv}+\mathrm{ox}}}{\mathrm{d}t}}={\displaystyle \frac{-M_{\mathrm{oc}}^{\mathrm{riv}+\mathrm{ox}}}{t_{\mathrm{0}}-t}}-k_{\mathrm{ox}}{M}_{\mathrm{oc}}^{\mathrm{riv}+\mathrm{ox}},\text{(18)}& {\displaystyle }{M}_{\mathrm{oc},t}^{\mathrm{riv}+\mathrm{ox}}=M_{\mathrm{oc},\mathrm{0}}\left(\mathrm{1}-{\displaystyle \frac{t}{t_{\mathrm{0}}}}\right)e^{-k_{\mathrm{ox}}t}.\end{array}$$

To summarize, the evolution of the mass of OC present in a landslide deposit depends on its connectivity status. When a landslide deposit is not actively eroded by a river, the OC evolution is governed by Eq. (11). Once the landslide is connected to the drainage network, the OC evolution follows Eq. (17) (Fig. 2b). Thus, the overall fate of landslide-mobilized OC is controlled by two timescales: (i) the connection time of a landslide to the fluvial network, t_con, and (ii) the export time of landslide fine sediment, t₀.

During the unconnected phase, the initial OC mass decreases only through oxidation, and the overall mass of organic matter that has been oxidized $\left(M_{\mathrm{oc},\mathrm{uncon}}^{\mathrm{ox}}\right)$ is estimated by integrating Eq. (13):

$$\begin{array}{}\text{(19)}& M_{\mathrm{oc},\mathrm{uncon}}^{\mathrm{ox}}=M_{\mathrm{oc},\mathrm{0}}{k}_{\mathrm{ox}}\underset{\mathrm{0}}{\overset{t_{\mathrm{con}}}{\int }}{e}^{-k_{\mathrm{ox}}t}\mathrm{d}t.\end{array}$$

In the case of a connected landslide, the organic matter contained within the deposit evolves as a function of the rate of fluvial transport and oxidation. We compute the total mass of OC exported by rivers $\left(M_{\mathrm{oc}}^{\mathrm{riv}}\right)$ and the oxidized OC mass $\left(M_{\mathrm{oc}}^{\mathrm{ox}}\right)$ to quantify the OC redistribution using the following equations. However, to account for the OC loss during the unconnected phase, we must substitute the initial OC mass (M_oc, 0) by the OC mass of remaining at the end of the unconnected phase $\left(M_{\mathrm{oc},\mathrm{0}}-M_{\mathrm{oc},\mathrm{uncon}}^{\mathrm{ox}}\right)$ in Eq. (18).

The $M_{\mathrm{oc}}^{\mathrm{riv}}$ is calculated by integrating Eq. (14) combined with Eq. (18):

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{M}_{\mathrm{oc}}^{\mathrm{riv}}& =\underset{t_{\mathrm{con}}}{\overset{t_{\mathrm{0}}}{\int }}{\displaystyle \frac{M_{\mathrm{oc},t}^{\mathrm{riv}+\mathrm{ox}}}{t_{\mathrm{0}}-t}}\mathrm{d}t\\ & ={\displaystyle \frac{\left(M_{\mathrm{oc},\mathrm{0}}-M_{\mathrm{oc},\mathrm{uncon}}^{\mathrm{ox}}\right)}{t_{\mathrm{0}}}}\underset{t_{\mathrm{con}}}{\overset{t_{\mathrm{0}}}{\int }}{e}^{-k_{\mathrm{ox}}t}\mathrm{d}t.\end{array}\end{array}$$

The $M_{\mathrm{oc}}^{\mathrm{ox}}$ is calculated by integrating Eq. (13) combined with Eq. (18):

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{M}_{\mathrm{oc}}^{\mathrm{ox}}& =M_{\mathrm{oc},\mathrm{uncon}}^{\mathrm{ox}}+\underset{t_{\mathrm{con}}}{\overset{t_{\mathrm{0}}}{\int }}{k}_{\mathrm{ox}}{M}_{\mathrm{oc},t}^{\mathrm{riv}+\mathrm{ox}}\mathrm{d}t=M_{\mathrm{oc},\mathrm{uncon}}^{\mathrm{ox}}\\ & +\left(M_{\mathrm{oc},\mathrm{0}}-M_{\mathrm{oc},\mathrm{uncon}}^{\mathrm{ox}}\right)k_{\mathrm{ox}}\underset{t_{\mathrm{con}}}{\overset{t_{\mathrm{0}}}{\int }}\left(\mathrm{1}-{\displaystyle \frac{t}{t_{\mathrm{0}}}}\right)e^{-k_{\mathrm{ox}}t}\mathrm{d}t.\end{array}\end{array}$$

In the case in which the timescale of the study (i.e., one or several seismic cycles) is shorter than t₀, Eqs. (20) and (21) are integrated over the considered timescale duration. In this case the amount of OC remaining on hillslopes is obtained by

$$\begin{array}{}\text{(22)}& M_{\mathrm{oc}}^{\mathrm{land}}=M_{\mathrm{oc},\mathrm{0}}-M_{\mathrm{oc}}^{\mathrm{riv}}-M_{\mathrm{oc}}^{\mathrm{ox}}.\end{array}$$

Equations (12) to (22) describe the theoretical framework that we use to quantify the OC redistribution during post-seismic periods. They are applied to several different oxidation models that we describe hereafter. The evolution of the OC stock in each pool is managed by Eq. (12) (unconnected phase) and Eq. (18) (connection phase) for each landslide deposit.

3 OC evolution of a single landslide deposit

Before tackling the widespread mobilization of OC by a landslide population, we first focus on theoretical predictions of OC evolution from a single landslide. This section focuses on the effect of the different oxidation models and their associated oxidation constants (k_ox), as well as on the role of the connection time (t_con) and export time (t₀) in the partitioning of mobilized OC between oxidation and fluvial export (Fig. 3).

Figure 3 The proportion of landside-mobilized OC that is exported by rivers for a single landslide and different OC oxidation models, plotted as a function of the total export time (t₀). (a) Results from the single pool model using different values of the oxidation rate constant: k_ox=0.5 % yr⁻¹, 1 % yr⁻¹ and 2 % yr⁻¹. Three different times of connection to the river network are also shown (t₀=0, 10 and 50 years). (b) Results from the multi-pool model. Here the landslide is initially connected to the river (t_con=0 years). The different lines represent model outcomes with different proportions of OC in the three modeled pools, as indicated by the greyscale text. Note that a larger fraction of OC in the high-reactivity pool leads to a decrease in the proportion of OC that can be exported by the river network.

The first case (Fig. 3a) shows the evolution of the post-seismic OC content of a landslide deposit using a single pool oxidation model for different values of t₀, t_con and k_ox. The general tendency is that the proportion of OC exported by rivers decreases with t₀, as longer export times enable more oxidation to take place. The latter effect is strengthened when the landslide deposit undergoes a prior phase of disconnection from the river network. For instance, in the case k_ox=2 % yr⁻¹, when the landslide is directly connected to the river network, the quantity of the eroded OC transported by rivers is ∼99 % when the export time is very short (t₀=1 year), and it decreases to 43 % when the export time is long (t₀=100 years). In contrast, when the deposit experiences a disconnection phase of 50 years, these proportions decrease to 36 % (t₀=1 year) and 16 % (t₀=100 years). The choice of the oxidation constant value also plays a role in the fate of OC: for t₀=100 years the proportion of OC exported by rivers is 79 % for k_ox=0.5 % yr⁻¹, whereas it is only 43 % for k_ox=2 % yr⁻¹.

The multi-pool oxidation model reveals similar patterns as the single pool model, with the mass of OC exported by the fluvial network decreasing with increasing values of t₀ (Fig. 3b). However, the inclusion of OC pools with lower reactivity results in more OC being exported by fluvial transport. When the low-reactivity pool is the largest, up to 75 % of OC is exported by rivers even for long export times (t₀=100 years). In contrast, with a large high-reactivity pool, only 50 % of OC is preserved long enough to be exported by rivers with the same value of t₀.

Drawing on observations from soil profiles from New Zealand (Harvey, 2019; Tonkin and Basher, 2001), we conceptualize soil profiles from the west coast catchments as comprising (1) a highly reactive O layer, which can constitute the first pool and represents ∼10 % of the total soil mass; (2) the A and E soil layers that constitute the second pool and represent ∼20 % of the total mass; and (3) the last pool that contains the low-reactivity organic matter that is present in the deeper B and C layers and represents the remaining ∼70 % of the total soil mass. Therefore, this soil configuration, with a large slowly reacting pool, is likely to preserve the OC in landslide deposits for long timescales and therefore promote a carbon sink.

4 Timescales of erosion as a CO₂ source or sink

The single landslide scenario can be used to explore the timescales over which erosion acts as a CO₂ source or sink. For instance, immediately following a landslide-triggering event, degradation of the organic matter will release CO₂ to the atmosphere (Zeng et al., 2009). Conceptual models suggest that this impacts the rate of CO₂ drawdown via vegetation regrowth and soil re-establishment (Restrepo et al., 2009). However, over longer timescales, the OC transported by rivers has the potential to contribute to net storage, and once the forest has regrown, landslide-driven erosion represents a net CO₂ sink (Berhe et al., 2007).

To place some quantitative constraint on these timescales for the first time, we explore the fate of landslide-mobilized soil OC with the single pool degradation model (k_ox=2 % yr⁻¹) for a range of fine sediment export times, t₀ and landslide connection times t_con. The processes acting as a CO₂ sink are soil re-establishment on landslide scars and fluvial transport, under the assumption that the OC experiences minor oxidation during fluvial transport (Scheingross et al., 2019). To describe the evolution of soil OC on hillslopes following a landslide event, we assume that it tracks the regrowth of forest biomass. This is supported by the observation that surface soil horizons from the western Southern Alps have a high radiocarbon content (with a “bomb” ¹⁴C signature) that suggests rapid OC turnover on timescales of less than 50 years (Frith et al., 2018; Horan et al., 2017). As such, the evolution of the C stock of the OC-rich surface soil horizons (to 0.1 m of depth) may quickly track any regrowth of vegetation that supplies the OC inputs.

To model aboveground biomass (AGB, in Mg ha⁻¹) growth through time, we call on an empirical model defined for a set of global measurements (McMahon et al., 2010):

$$\begin{array}{}\text{(23)}& \mathrm{AGB}={\mathit{\beta }}_{\mathrm{1}}\left({\displaystyle \frac{t}{t+\mathit{\theta }}}\right),\end{array}$$

where β₁ is the asymptotic maximum biomass that a stand can achieve, t is time and θ is the age at half-saturation of the function. This model is generally used to describe the recovery of AGB following a disturbance. Based on Eq. (23), the full re-establishment of pre-event OC stocks is likely to take ∼150–200 years for aboveground biomass (McMahon et al., 2010). If we define this as the timescale over which carbon stocks approach pre-disturbance levels, t_recovery, we can explore the net CO₂ exchange associated with landslide-mobilized OC.

These calculations demonstrate how landslide erosion can be viewed as both a CO₂ source and sink (Wang et al., 2016) depending on the timescales of consideration. Without fluvial transport, the competition between aboveground biomass growth and recovery of the surface organic soil horizons versus OC degradation controls the presence and duration of a transient period of CO₂ release and CO₂ drawdown (Fig. 4a). A faster rate of soil re-establishment allows the C stock on hillsides to build up faster than the oxidation of eroded soil OC, and the erosion event initially looks like a small sink. In contrast, slower soil recovery is associated with a period of CO₂ release, whose extent depends on the value of θ. However, by the nature of the model setup, the net effect over the duration of a seismic cycle (i.e., ∼270 years) is no net CO₂ drawdown; the landslide-mobilized OC is oxidized, while the landslide scar builds up an equivalent new OC stock.

Figure 4 The modeled transient C budget of eroded organic carbon (partitioning between a CO₂ source or sink) over a seismic cycle with a duration of 270 years. Red lines are the CO₂ source from oxidation of the mobilized OC. Green lines indicate CO₂ sequestration by the re-establishment of the soil OC stocks in the landscape, which is based on a model of aboveground biomass accumulation. (a) The case in which only OC oxidation is accounted for. The parameters controlling the shape of the soil OC stock model are explored and determine whether landsliding is a net C source (solid black line) or net sink (solid dashed line). (b) The case in which OC oxidation and fluvial export are both active. Here, the parameters controlling the physical export of sediment (t₀ and t_con) are explored. Solid lines correspond to t₀=100 years and t_con=50 years, resulting in a short-term C source but long-term sequestration. Dashed lines correspond to t₀=20 years and t_con=10 years, resulting in net sequestration over the entire seismic cycle.

By including river sediment transport, the model produces a net long-term sink of CO₂ (Fig. 4b). In this case, a CO₂ source period only exists if the connection phase is long enough to allow OC oxidation to play a significant role (Fig. 4b). For instance, the initial CO₂ source lasts 85 years in the case of t_con=50 years but lasts <10 years when t_con=10 years. The size of the longer-term CO₂ sink depends on the river transport efficiency, t₀, and on the connection time t_con.

In summary, if t_recovery is ∼200 years, it suggests that landslide erosion results in a short-term transient source of CO₂ for less than ∼50 years and that the longer-term CO₂ sink persists over timescales >200 years. It is important to note that the AGB model applied here does not capture the formation of deeper, mineral-associated soil OC, which may take longer periods of time and could contribute importantly to soil OC stocks (Harvey, 2019). Nevertheless, all the cases explored with our approach show that, over the ∼270-year duration of a seismic cycle in this setting, landslides always act as a carbon sink.

5 The fate of OC mobilized by a landslide population

The view from a single landslide deposit allows us to assess the fate of eroded OC for different configurations of landslide connectivity, fluvial transport efficiency and rates of OC oxidation (Figs. 2 and 3). However, for a distribution of landslides triggered by an earthquake, we expect these to vary due to the spatial variability of river transport capacity (i.e., upstream drainage area upstream of a landslide deposit), the timescale of deposit connection, and varied landslide areas and volumes. In this section, we explore what this complexity means for the fate of mobilized OC considering a large landslide population triggered by a single earthquake.

Using the Quakos framework (Croissant et al., 2019b), we model a M_w 7.9 earthquake on the Alpine Fault that is consistent with paleoseismic reconstructions (Howarth et al., 2012, 2018). Our modeling approach predicts that ∼ 58 000 landslides would be triggered by the earthquake, impacting an area of ∼11 000 km² and mobilizing ∼1.2 km³ of fine sediment. Overall, the mass of OC mobilized from soil organic matter for this scenario is ∼3.3 MtC. From this starting point, this section examines the role of various parameters in the post-seismic redistribution of OC. We denote the cumulative timescale necessary to fully evacuate the landslide-mobilized OC as $t_{\mathrm{tot}}=t_{\mathrm{con}}+t_{\mathrm{0}}$. The OC evolution is tracked over one seismic cycle with a duration of 270 years.

5.1 The role of landslide connectivity

The spatial variation of the local transport rate ($\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}$), combined with a landslide volume (V_ls) population, leads to the emergence of a distribution of landslide export times t₀ and of distances between deposits and the fluvial network. This in turn controls the connection time t_con. Here, we explore the fate of OC as a function of the connection velocity of landslides to the river network (u_con), which also controls the distribution of t_con.

For each individual landslide generated by Quakos, we calculate the proportion of OC exported by rivers as opposed to that oxidized to CO₂ as a function of $t_{\mathrm{tot}}=t_{\mathrm{con}}+t_{\mathrm{0}}$. The general pattern that emerges is that the proportion of OC exported by rivers decreases with t_tot, ranging from 100 % of OC exported by rivers for values of t_tot<1 year to 0 % for t_tot>300 to 10 000 years depending on the value of t_con (Fig. 5). This approach allows us to assess an envelope of possible values of the proportion of OC exported by rivers as a function of t_tot, defined by two cases: (i) when t_con is large compared to the export time (red line in Fig. 5) and (ii) when the landslides are directly connected to the rivers when t_con=0 (blue line in Fig. 5). The scatter in the data emerges from the competition between the two timescales controlling the fate of the landslide-mobilized OC. A dominance of t_con promotes the oxidation (and/or storage on hillslopes) of the organic matter, while a lower value of t₀ leads to a higher proportion of OC that is exported by the fluvial network. Thus, for a fixed climate (mean annual runoff of 1 m yr⁻¹ and runoff variability k=1 in Eq. 7) the quantity of OC that is exported by the fluvial network is sensitive to the chosen connection velocity of the landslides to the rivers (Fig. 5).

Figure 5 The proportion of landslide-mobilized OC that is exported by rivers for a population of co-seismic landslides, plotted as a function of the total time the landslides remain in the catchment. All results are for a single pool model of OC oxidation and different values of the connection velocity: (a) u_con=100 m yr⁻¹, (b) u_con=1 m yr⁻¹. Each dot is an individual landslide characterized by a combination of its connection time (t_con) and fluvial export time (t₀), and its position on the x axis is $t_{\mathrm{tot}}=t_{\mathrm{con}}+t_{\mathrm{0}}$. The color of each dot represents the ratio between the export time and the connection time. The blue line is the percentage of OC exported by rivers when t_con=0. In this case, the x axis corresponds to t₀ only. The red line represents the lower limit of the percentage of OC that could be exported by rivers for a given time. In this case, t₀=0 and the x axis corresponds to t_con only. The histograms on the right of each plot represent the distribution of the percentage of OC exported by rivers from the scatter plot.

We next investigate the impacts of the connection velocity. The fate of carbon is considered to be either returned to the atmosphere as CO₂, retained in the landscape (in landslide deposits) or exported in river sediments (potential longer-term OC burial). Not surprisingly, the results for the single pool OC degradation model show that a lower connection velocity and a higher oxidation constant produce a larger proportion of C that has been oxidized and released in the atmosphere (Fig. 6). Interestingly, scenarios that show a significant storage of OC at the end of a seismic cycle only emerge for a low oxidation constant coupled to a slow connection velocity (Fig. 6). The multi-pool degradation model outputs do not differ significantly from the single pool results. However, predicted OC storage in the mountain landscape is larger when the low-reactivity pool is the largest (Fig. 7). For both oxidation models, high connection velocities (>100 m yr⁻¹) systematically promote the river transport of OC, with between 80 % and 98 % of OC being exported out of the catchments. The physical meaning of the connection velocity for the case of fine sediment is discussed in Sect. 7.1.

Figure 6 Histograms showing the fate of the landslide-mobilized OC at the end of a seismic cycle with a duration of 270 years partitioned between fluvial export, oxidation and landslide deposit storage. These results have been computed with the single pool model for different connection velocities and different oxidation constants: (a) k_ox=2 % yr⁻¹; (b) k_ox=1 % yr⁻¹; (c) k_ox=0.5 % yr⁻¹; (d) k_ox=0.1 % yr⁻¹.

Figure 7 Histograms showing the fate of the landslide-mobilized OC at the end of a seismic cycle with a duration of 270 years partitioned between fluvial export, oxidation and hillslope storage. These results have been computed for different connection velocities with a multi-pool model composed of three pools with oxidation constants k_ox=2 % yr⁻¹, 1 % yr⁻¹ and 0.5 % yr⁻¹. Panels (a)–(d) show model outcomes with different proportions of OC in the three modeled pools, as indicated by the greyscale text: (a) equal distribution between pools, (b) dominance of the low-reactivity pool, (c) dominance of the intermediate-reactivity pool and (d) dominance of the high-reactivity pool.

By considering a full landslide population, we can also consider the effect of landslide size on the fate of OC. Our results suggest that intermediate landslides (i.e., with an area ranging from 10³ to 10⁵ m²) are responsible for most of the export of river carbon (Fig. 8). While smaller landslides (A_ls<10³ m²) are removed efficiently by rivers, they mobilize a lower initial quantity of OC and are less important in the total carbon budget. In contrast, large landslides (A_ls>10⁵ m²) individually mobilize large quantities of OC but have longer export times as t₀ is a function of V_ls, fine, potentially reducing their impact on the carbon export by fluvial transport.

Figure 8 Quantification of the export of OC by river transport for different classes of landslide area (box plots) and for different oxidation constants of the single pool model (a, b, c: k_ox=0.5 % yr⁻¹; d, e, f: k_ox=2 % yr⁻¹). The yellow diamonds represent the contribution of each landslide area class to the total OC exported by rivers. The boxes show the 25th–75th percentiles, the whiskers show the 10th–90th percentiles, the dot shows the median and the horizontal bar shows the mean.

5.2 The role of climate in landslide OC evacuation

In our model, the mean annual runoff and runoff variability determine the sediment transport rate (Eq. 6 that depends on Eq. 7) and the river export of OC (Hilton, 2017; Wang et al., 2019). The fate of landslide-mobilized OC is sensitive to the export time, which is a ratio between the volume of fine sediment available for transport and the sediment transport rate. Here, we assess how the fate of OC is moderated by these climatic parameters. For the mean annual runoff, we chose a range of values that span the observed values of different mountain ranges (e.g., Bookhagen and Burbank, 2010; Hicks et al., 2011). For the runoff variability, k varies between 0.5, a value implying a high recurrence of large daily runoff events typical of Taiwan (Lague et al., 2005), and 4, a value representative of temperate regions that rarely experience flooding events (Croissant et al., 2019a).

The role of runoff and runoff variability is explored for a slow landslide connection velocity (u_con= 1 m yr⁻¹), ranging to the “full connectivity” case. Generally, the proportion of OC exported by rivers increases for higher mean annual runoff and higher runoff variability (Fig. 9). For mean annual runoff values below 1 m yr⁻¹, the proportion of OC export by rivers is not significantly affected by runoff variability. In general, the results show more sensitivity to the mean annual runoff than to runoff variability. This emerges from the parameterization of the long-term fluvial transport capacity (Eqs. 6 and 7). The suspended sediment transport law (Eq. 5) does not have a threshold, which limits the role of large runoff events in the long-term transport capacity. Therefore, the long-term transport capacity is controlled by intermediate runoff events, which are set by the mean annual runoff value (see further explanation in Croissant et al., 2019a).

The connection velocity has a major impact on the modeled fate of OC. However, climate strongly moderates this response. For instance, the full connectivity scenario shows that the proportion of the landslide-mobilized OC that is exported by rivers ranges from 30 % for low runoff and low discharge variability to ∼95 % for high runoff and variability (Fig. 9c). These model outputs predict that most feasible combinations of the climate parameters lead to a dominance of the fluvial export of OC. The low connection velocity produces a specific response (Fig. 9a). This is because the amount of exported material is limited to only the initially connected landslides (e.g., ∼30 % in these simulations). We discuss the broader implications of these results in Sect. 7.

It is important to note that in these simulations we do not attempt to capture any role of climate in the build-up of the soil OC stock in the landscape. Mean annual temperature and rainfall act to moderate the production of OC and its preservation in soil (e.g., Carvalhais et al., 2014). In addition, these variables control the development of mineral soil (e.g., Mudd and Yoo, 2010), and mineral–OC interactions are critical to organic matter stabilization (Hemingway et al., 2019). As such, a wetter climate might act to enhance erosion and export of landslide-mobilized OC (Fig. 9) but also increase its export from hillslopes by controlling soil OC stocks.

Figure 9 Proportion of landslide-mobilized OC that has been exported by rivers over one seismic cycle as a function of mean annual runoff and runoff variability. Three scenarios of landslide connection velocity are shown: (a) u_con=1 m yr⁻¹, (b) u_con=10 m yr⁻¹ and (c) full connectivity. The lines represent contours of the proportion of OC exported by rivers and are incremented every 10 %, but the dashed line represents the 50 % contour.

6 Calculating long-term fluxes – application to the Southern Alps

In the previous section, OC partitioning between the different reservoirs was expressed as percentages rather than absolute values to illustrate the flexibility of the approach. In this section, we apply the model to the New Zealand case for which we can constrain parameters using empirical datasets. In addition, we consider an upscaling approach to quantify the fate of landslide-mobilized OC over several seismic cycles. We use the single pool oxidation model for simplicity (k_ox=1 % yr⁻¹) and refer readers to Sect. 3 for further discussion of the choice of degradation models. The mean annual runoff is set at 5 m yr⁻¹ and the runoff variability to k=1 to mimic the west coast of the Southern Alps (Croissant et al., 2017). The tectonic scenario consists of a series of M_w 7.9 earthquakes occurring on the Alpine Fault every 270 years. We also explore three different landslide connectivity scenarios ranging from full connectivity to u_con=1 m yr⁻¹. Based on soil OC samples acquired in this study area, we also make the assumption that the OC is distributed equally on all slopes (Harvey, 2019).

Our approach models the pulsed, large-scale mobilization of OC from earthquake-triggered landslides (Fig. 10). The OC can be partitioned between that exported by rivers and that released to the atmosphere as CO₂ by oxidation. For every scenario, earthquakes are identified by a pulse of OC exported by rivers that reaches a maximum in the year following the seismic event. The peak is larger for the full connectivity scenario, which delivers a larger and immediate supply of sediment to the drainage network. After each earthquake, the exported OC decreases and tends towards zero at the end of each seismic cycle. The model shows the vast majority of OC export by rivers can happen during the first decade of the post-seismic period. This is consistent with observations of carbon accumulation from a lake core downstream of a small catchment draining the Alpine Fault (Frith et al., 2018) and organic geochemistry tracers of sediment provenance (Wang et al., 2020). In the model, the fast initial fluvial evacuation of OC is due to a high proportion of fine sediment being initially connected to the fluvial network, combined with a climate that promotes high sediment transport rates. In the full connectivity scenario, the total volume of landslide-derived sediment is supplied to the rivers during the co-seismic phase. The decrease in the volume of landslide-derived sediment is only controlled by the distribution of t₀ of the landslides, as no new material is introduced in the fluvial network through time. In contrast, the two other scenarios include a progressive connection of new landslides to the river network with time and therefore present higher OC export rates than the full connectivity case after ∼10 years.

Figure 10 The modeled fate of landslide-mobilized OC over several seismic cycles based on parameters set by the Alpine Fault, western Southern Alps, New Zealand, with the single pool model and k_ox=1 % yr⁻¹. Each panel shows the fluxes of OC exported by rivers (in blue) and the release of CO₂ to the atmosphere from OC oxidation (in red) for different connection velocity scenarios: (a) full connectivity, (b) u_con=10 m yr⁻¹, (c) u_con=1 m yr⁻¹. The envelope for each plot accounts for the uncertainties in the soil OC content.

The amount of OC that is predicted to be oxidized is several orders of magnitude lower than the physical export of OC. However, oxidation becomes the dominant process after the initial large wave of fluvial export in the first decade that follows each earthquake (Fig. 10). In addition, oxidation fluxes are larger when the connection velocity decreases and river export becomes less efficient at removing OC.

Over one seismic cycle, the model approach suggests that the vast majority of the OC mobilized by earthquake-triggered landslides is exhausted from the landslide deposits for every scenario; i.e., there is no net accumulation of fine OC in the landslide deposits during the time span of 270 years. Once again, the landslide connection velocity influences the fate of OC, with ∼90 % of the OC exported by rivers for the full connectivity scenario, while this value decreases to ∼28 % for the slow connection velocity case. It is worth noting that without any sediment transport, the single pool degradation model predicts that ∼89 % of the organic carbon would have been oxidized and returned to the atmosphere and the 11 % of remaining OC stored in landslide deposits over timescales longer than the seismic cycle.

7 Discussion

7.1 Comparison to previous work

Previous work has sought to quantify the provenance and flux of particulate OC following large storm events based on river gauging station data (Clark et al., 2017; Hilton et al., 2008b), remote sensing approaches (Clark et al., 2016; Ramos Scharrón et al., 2012; West et al., 2011) or lake stratigraphic records (Frith et al., 2018; Wang et al., 2020). These studies provide insights on the rates and processes of OC transport to the oceans (Galy et al., 2015; Hilton, 2017). However, they are generally unable to constrain the roles of OC oxidation in the landscape (Fig. 2) and aboveground biomass regrowth (Fig. 4) in the OC budget of a landscape recovering from these damaging events. In our study, we propose a theoretical description of these processes using a simplified numerical approach.

The export of OC during the post-seismic period depends strongly on the fine sediment export rate (Hovius et al., 2011; Tolorza et al., 2019; Wang et al., 2015, 2016). A few examples have shown that, in wet climates, suspended load fluxes after a large earthquake are characterized by a rapid increase directly after the seismic event, which is sustained for less than a decade, before returning to background levels (Hovius et al., 2011). This behavior is reproduced by our simulations (Fig. 10). However, this does not imply that all the fine sediment has been evacuated out of the epicentral area. For instance, abundant landslide deposits located in headwater areas (Ramos Scharrón et al., 2012), or fine particles trapped under coarser grains in the deposits, would slow down the evacuation of sediment once the easily accessible particles have been entrained. In that regard, the scenarios that we propose in Sect. 6 show possible trajectories of the partitioning between oxidation and fluvial export, with oxidation being the dominant process after an initial pulse of fluvial evacuation of OC. However, we note that the initial pulse is present in all simulations, independent of the connection velocity, and implies that in all cases fluvial export is the dominant fate of OC mobilized by landslides in a setting like the Southern Alps. These findings are consistent with the available measurements from the 2008 Wenchuan earthquake (Wang et al., 2016) and fit the pattern of organic matter accumulation seen in lake deposits that record four past Alpine Fault earthquakes (Frith et al., 2018).

Other research has sought to assess the role of landslide erosion as a CO₂ source or sink. For instance, Ramos Scharrón et al. (2012) proposed a post-landsliding OC evolution model that accounts for the competition between carbon gains and losses. Their modeling approach differs from ours because the processes that result in a CO₂ sink only include soil formation, aboveground biomass growth on landslide scars and OC storage in landslide deposits. In their approach, fluvial transport contributes to oxidation as part of the OC source term. Despite these important differences, their results also suggest that widespread mass wasting is a CO₂ sink over long timescales in settings where landslide deposits have long residence times (centuries), which would correspond to our low connection velocity scenario. Our results differ from theirs in that we model connection of landslides to the mountain river network, leading to scenarios in which up to 100 % of the OC mobilized by landslides is removed by rivers rather than oxidized in the landscape.

7.2 Extreme erosion events and the global OC cycle

The mobilization of OC by events that trigger widespread landsliding can produce large pulses of carbon export from the biosphere (Clark et al., 2016; Frith et al., 2018; Wang et al., 2016; West et al., 2011). These carbon fluxes are large enough to raise questions about the role of large erosive events associated with earthquakes and storms in the long-term carbon cycle (Hilton and West, 2020). In this study, we propose a framework to examine the fate of OC eroded during a large earthquake and isolate how the climatic context plays a role in promoting fluvial transport. We do not attempt to quantify any subsequent degradation of OC in the marine realm or within sedimentary deposits. In the context of “short” mountain ranges along active continental margins with sediment transport distances of <100 km and rivers that are coupled to deep-ocean sedimentary sinks (e.g., Mountjoy et al., 2018), it has been shown that 70 % or more of the OC exported in river sediments is buried in offshore sedimentary deposits and is therefore considered a long-term sink (Blair and Aller, 2012; Galy et al., 2015; Kao et al., 2014).

Our modeling exercise has shown that earthquakes can contribute to long-term biospheric carbon sinks in the Southern Alps case, with minimal within-catchment oxidation of eroded OC (Fig. 10). The model is broadly relevant to other mountain catchments in the western Pacific: Taiwan and Papua New Guinea (Dadson et al., 2005; Ferguson et al., 2011; Hicks et al., 2011; Hilton et al., 2008b). However, these settings also share geomorphic, tectonic and climatic attributes with islands of the Caribbean (e.g., Allemand et al., 2014) and of Central America (e.g., Ramos Scharrón et al., 2012), as well as some catchments draining the western American continent (Leithold et al., 2006). In these conditions, the combination of short sediment travel distance and limited potential for fluvial sediment storage may promote efficient carbon transport from land through river systems (Scheingross et al., 2019). However, a complete understanding of the role of earthquakes in the carbon cycle must consider large floodplains, such as the river systems that drain the Himalayas and the Andes. For instance, Galy et al. (2008) showed evidence for the replacement of Himalayan biospheric OC in the Ganges floodplain, and particulate OC eroded from the Andes is thought to suffer a similar fate during transport through the Amazon floodplain (Bouchez et al., 2014; Ponton et al., 2014). Drier climates may also contrast with the findings we describe because the OC stocks in soils may be lower, and landslide debris may reside in catchments for much longer periods of time.

Our model has been applied to a single case study with a specific tectonic and climatic regime (Fig. 10). It thus represents a first step in understanding how widespread landslide-triggering events impact regional to global carbon fluxes. Ideally, future work can explore the diversity of mountain ranges that present different vegetation, soil cover, topography and climate, with a specific consideration of the recurrence time of earthquakes (Fan et al., 2018). Additionally, a complete carbon budget (Hilton and West, 2020) could be assessed if other long-term CO₂ sinks via silicate mineral weathering that occurs in landslide deposits are considered (Emberson et al., 2017), alongside CO₂ release from sulfide mineral oxidation (Emberson et al., 2016a, 2016b) and oxidation of rock OC (Hilton et al., 2014).

7.3 Limitations of the approach

In this paper, we couple the reduced-complexity model Quakos to organic matter degradation models to assess the fate of OC following a large landslide-triggering event. As with any modeling approach, our model contains some limitations that are discussed here.

Prediction of the landslide pattern for an event that has not taken place during the period of modern instrumental records is a difficult endeavor. We were able to constrain the shape of the landslide distribution using observations of rainfall-triggered landslides in the Southern Alps (Hovius et al., 1997). Nevertheless, uncertainties remain as to the landslide spatial density that would emerge from a M_w∼8 earthquake on the Alpine Fault (Robinson and Davies, 2013). While this parameter does not significantly affect the amount of OC redistributed in the different reservoirs (Fig. 6), it would greatly impact the OC fluxes and therefore their comparison to background fluxes. Our computations show that for the “full connectivity” scenario, maximum post-seismic OC fluxes derived from the fluvial export would correspond to ∼2 times (for a peak density of 4 %) to ∼6.5 times (for a peak density of 8 %) the background fluxes computed from Hilton et al. (2008a).

One of the main limitations of our modeling approach is in the assumption of a connection velocity u_con. As discussed in Croissant et al. (2019b), the lack of a theoretical framework describing the evolution of non-cohesive material sitting on hillslopes after a landslide-triggering event makes the description of this process difficult. For fine particles, there are several key processes that need to be considered across a range of scales. These include those operating at the grain scale (e.g., creep, rain splash) to those operating over length scales of meters (e.g., dry raveling, runoff-driven erosion) to hundreds of meters (landsliding, debris flow). The availability of fine sediment is also likely to change over time, for instance through the initial removal of easily transportable debris (Wang et al., 2015; West et al., 2014; Zhang et al., 2019). To move more finer material, subsequent mass movements may be needed to expose it within deposits and/or within mountain river channels (Fan et al., 2019). It is tempting to compare the values of connection velocity chosen here (1 to 100 m yr⁻¹, or $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{8}}$ to $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹) to the velocity of water flow sediment export via a range of processes (e.g., DiBiase et al., 2017). In this case, the chosen connection velocities are very slow. However, they seek to describe the net motion of a large sediment pile that integrates a range of processes that occur discretely in space and time.

Numerous physically based and empirically based models have been proposed to explore and predict hillslope transport processes (e.g., Aksoy and Kavvas, 2005; DiBiase et al., 2017; Tucker and Bradley, 2010). Future work could seek to incorporate more sophisticated transport laws into the framework we describe here. This would also allow the links between landslide remobilization and climate to be explored, as there is likely to be a strong link between movement of existing landslide debris and events with high rainfall intensity and a high river stage (Marc et al., 2015). Useful insights could be gained using morphodynamic models with a realistic grain size distribution (Fan et al., 2019). An additional source of uncertainty relates to the dynamics of particulate OC transport versus clastic sediment. While the general link between the two phases has been demonstrated (Galy et al., 2015; Hilton et al., 2012), there is also evidence that discrete clasts of woody debris can be sorted by flowing water (Hilton et al., 2015; Lee et al., 2019), and its transport behavior is different from the clastic load due to a lower density and different particle shape (Turowski et al., 2016).

The model approach assumes that the majority of fine sediment is available for transport, i.e., located at the surface or near the subsurface of landslide deposits. In reality, fine sediments will be distributed within landslide deposits, and some are therefore shielded from entrainment. To export this buried fine sediment, removal of some coarser fraction would have to occur, most likely during high-intensity rainfall events that have the capacity to export the coarse fraction. This feature would have the effect of prolonging the storage time of fine sediment and OC on hillslopes. Accounting for this process would require prior knowledge of the distribution of fine sediment within the deposit, which is still largely unknown (Fan et al., 2019), and a more sophisticated sediment transport model to account for the export of the coarser fraction of the sediment (Croissant et al., 2017).

In our simulations, we simplify how the degradation of OC proceeds following its mobilization of soil from hillslopes. First, it is possible that the landsliding process sorts and mixes the organic matter in such a way to promote oxidation of some of the materials but to protect some from degradation. For instance, observations in the field suggest that large woody debris may concentrate near the surface and toe of the deposit after a landslide. Instead, finer material may mix deeper in the deposit (Hilton et al., 2008b). If the landslide deposit has low porosity and/or receives runoff focused from the landslide scar (Emberson et al., 2016a), waterlogging may promote preservation of this material under anoxic conditions. This would act to enhance the proportion of OC exported by erosion processes. Second, we do not account for OC contained within sedimentary bedrock. In the western Southern Alps, the OC content of rocks is low (∼0.15 %), but oxidative weathering is thought to be enhanced by high erosion rates (Hilton et al., 2014; Horan et al., 2017). In the future, improved constraints on the reactivity of rock-derived OC (Hemingway et al., 2018) and models that link the production of fine clastic sediment to weathering reactions (Carretier et al., 2018) could be used to explore the role of large landslide populations in this mechanism for CO₂ release.

Finally, we assume there is no OC oxidation during river transport. Our New Zealand case study is focused on short rivers with maximum distances from headwaters to oceans of less than 100 km. This assumption is supported by recent work by Scheingross et al. (2019), who showed that within-river oxidation would be minimal, with between 0 and 10 % of mobilized OC being oxidized during transport over 1000 km. We also note that in the simulations, intermediate discharge events dominate the long-term transport rates (Croissant et al., 2019a), meaning that most of the OC would be exported by flows that do not escape the river channel and thus not go overbank onto the floodplains. This dynamic could be altered if a sediment bed-load wave significantly modifies the riverbed elevation (e.g., Hancox et al., 2005; Korup, 2004). Finally, we recognize that in the Amazon River floodplain, high-altitude-derived OC from mountain catchments is thought to be lost and/or overprinted by lowland OC (e.g., Bouchez et al., 2014; Ponton et al., 2014). This means the geomorphic configuration of mountain catchments impacted by earthquake events could play an additional role in the fate of the mobilized OC.

8 Conclusions

In this study, we propose a new theoretical framework to quantify the fate of soil-derived OC mobilized by earthquake-triggered landslides. The approach combines an empirical model to compute a landslide population triggered by large earthquakes and a suspended sediment transport law with different models of OC degradation. Overall, this model allows us to quantify the fate of landslide-mobilized OC in terms of its potential contribution to atmospheric CO₂ (by respiration of OC), its transient storage in landslide deposits and its contribution to particulate OC in rivers (a potential CO₂ sink). The model is computationally efficient, and these features can be explored for a range of climatic forcing (in terms of mean annual runoff and runoff variability) and OC degradation rates over multiple earthquake cycles.

At the scale of a single landslide, we find that the fate of OC is strongly controlled by the physical processes acting on the deposit. These include its initial connectivity to the river network and its potential to connect through time. This dynamic landslide connectivity is modeled using a “connection velocity” term, which seeks to capture a range of processes that may redistribute fine clastic sediment and OC from landslide deposits. The connection velocity and local river transport rate combine to set the time available for OC oxidation. The oxidation constants also play a significant role in the ultimate fate of OC, and a multi-pool degradation model can capture the persistence of OC throughout a seismic cycle of hundreds of years.

At the scale of an entire earthquake-triggered landslide population, the connection velocity and the mean annual runoff control the redistribution of OC. A wet climate and fast connection of fine material to the fluvial network both promote the riverine export of OC. Importantly, the type and complexity of the biogeochemical degradation model do not significantly affect the results. We apply the model to the case study of the western Southern Alps, New Zealand. The simulations suggest that much of the OC eroded from mountain ranges typical of the southwestern Pacific would be efficiently transferred to the ocean during the first 1–10 years after a large earthquake. In this context, depending upon the fate of OC downstream, earthquakes are likely to promote net carbon sequestration. An extended application of this methodology tailored to other mountain ranges in different tectonic and climatic contexts would allow for a more precise determination of the role of earthquakes and widespread landslide events in the regional and global carbon cycles.