Thomas Croissant
email: thomas.croissant@durham.ac.uk

Philippe Steer

Dimitri Lague

Philippe Davy

Louise Jeandet

Robert G Hilton

Seismic cycles, earthquakes, landslides and sediment fluxes: linking tectonics to surface processes using a reduced-complexity model

Keywords: landslide, earthquake, river morphodynamics, landscape evolution

In tectonically active mountain ranges, landslides triggered by earthquakes mobilise large volumes of sediment that affect river dynamics. This sediment delivery can cause downstream changes in river geometry and transport capacity that affect the river efficiency to export this sediment out of the epicentre area. The subsequent propagation of landslide deposits in the fluvial network has implications for the management of hazards downstream and for the longer-term evolution of topography over multiple seismic cycles. A full understanding of the processes and time scales associated with the removal of landslide sediment by rivers following earthquakes however, is still lacking. Here, we propose a nested numerical approach to investigate the processes controlling the post-seismic sediment evacuation at the mountain range scale, informed by results from a reach scale model. First, we explore the river morphodynamic response to a landslide cascade at the reach-scale using a 2D modelling approach. The results are then used to describe empirically the evacuation of a landslide volume which avoids using a computationally extensive model in catchments which may have thousands of co-seismic landslides. Second, we propose a reduced-complexity model to quantify evacuation times of earthquake-triggered landslide clusters at the scale of a mountain range, examining the hypothetical case of a M w 7.9 earthquake and its aftershocks occurring on the Alpine Fault, New Zealand. Our approach combines an empirical description of co-seismic landslide clusters with the sediment export processes involved during the post-seismic phase. Our results show that the inter-seismic capacity of the mountain range to evacuate co-seismic sediment is critical to assess the sediment budget of large earthquakes, over one to several seismic cycles. We show that

Introduction

Earthquakes impact the landscapes of active mountain ranges by mobilizing large volumes of sediment through widespread landsliding [START_REF] Malamud | Landslides, earthquakes, and erosion[END_REF][START_REF] Ouimet | The hills came tumbling down[END_REF] (Fig. 1). These clusters of landslides can deliver sediment to the fluvial network and affect its dynamics over timescales lasting from decades to centuries [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF][START_REF] Hovius | Prolonged seismically induced erosion and the mass balance of a large earthquake[END_REF][START_REF] Wang | Controls on fluvial evacuation of sediment from earthquake-triggered landslides[END_REF][START_REF] Wang | Perturbation of fluvial sediment fluxes following the 2008 Wenchuan earthquake[END_REF][START_REF] Yanites | How rivers react to large earthquakes: Evidence from central Taiwan[END_REF]. Whereas numerous studies have focused on the coseismic response of hillslopes to large earthquakes (e.g. [START_REF] Gallen | Coseismic landslides reveal near-surface rock strength in a highrelief, tectonically active setting[END_REF][START_REF] Keefer | Landslide causeb by earthquakes[END_REF][START_REF] Keefer | Earthquake-induced landslides and their effects on alluvial fans[END_REF]Malamud et al., 2004a;Marc et al., 2016;[START_REF] Meunier | Regional patterns of earthquake-triggered landslides and their relation to ground motion[END_REF][START_REF] Meunier | Topographic site effects and the location of earthquake induced landslides[END_REF], comparatively few have investigated the post-seismic evolution of landscapes and sediment transport (e.g. [START_REF] Yanites | How rivers react to large earthquakes: Evidence from central Taiwan[END_REF][START_REF] Hovius | Prolonged seismically induced erosion and the mass balance of a large earthquake[END_REF][START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF]. This probably results from extensive datasets constraining co-seismic landsliding and sediment production [START_REF] Keefer | Earthquake-induced landslides and their effects on alluvial fans[END_REF][START_REF] Larsen | Landslide erosion controlled by hillslope material[END_REF][START_REF] Li | Seismic mountain building: Landslides associated with the 2008 Wenchuan earthquake in the context of a generalized model for earthquake volume balance[END_REF] A C C E P T E D M A N U S C R I P T

3 Malamud et al., 2004a;Marc et al., 2016;[START_REF] Tanyaş | Presentation and Analysis of a Worldwide Database of Earthquake-Induced Landslide Inventories[END_REF][START_REF] Tanyaş | An updated method for estimating landslideevent magnitude[END_REF], while it is more challenging to obtain measurements of pre-and post-seismic sediment fluxes and topographic changes over tens of years. Most of the observations relate to suspended sediment concentration data [START_REF] Dadson | Earthquake-triggered increase in sediment delivery from an active mountain belt[END_REF][START_REF] Hovius | Prolonged seismically induced erosion and the mass balance of a large earthquake[END_REF][START_REF] Wang | Controls on fluvial evacuation of sediment from earthquake-triggered landslides[END_REF][START_REF] Wang | Perturbation of fluvial sediment fluxes following the 2008 Wenchuan earthquake[END_REF]) completed with analytical and numerical models of bedload evacuation [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF][START_REF] Yanites | How rivers react to large earthquakes: Evidence from central Taiwan[END_REF].

Investigating how fluvial systems digest these abrupt and large sediment pulses is critical for sediment transfers [START_REF] Benda | Stochastic forcing of sediment supply to channel networks from landsliding and debris flow[END_REF], bedrock incision patterns in landslide-dominated mountain ranges [START_REF] Lague | Reduction of long-term bedrock incision efficiency by short-term alluvial cover intermittency[END_REF][START_REF] Yanites | The influence of sediment cover variability on long-term river incision rates: An example from the Peikang River, central Taiwan[END_REF], hydro-sedimentary hazards in alluvial fans [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF][START_REF] Robinson | Review Article: Potential geomorphic consequences of a future great (Mw Combining double low line 8.0+) Alpine Fault earthquake, South Island, New Zealand[END_REF] and even to quantify the feedbacks of surface processes on fault stress loading (Steer et al., 2014). They are also relevant to geochemical fluxes from mountain belts, particularly those which relate to soil erosion and organic carbon transfer [START_REF] Wang | Earthquake-triggered increase in biospheric carbon export from a mountain belt[END_REF] and inorganic carbon via silicate and carbonate weathering which can take place in landslide deposits (Emberson et al., 2016a;[START_REF] Jin | Seismically enhanced solute fluxes in the Yangtze River headwaters following the A . D[END_REF]. Over short-times scales (i.e. < 1000 years), the downstream propagation of sediment pulses have been studied principally at the reach scale using flume experiments and 1D numerical modelling [START_REF] Cui | Sediment pulses in mountain rivers: 2. Comparison between experiments and numerical predictions[END_REF][START_REF] Cui | Numerical Model of Sediment Pulses and Sediment-Supply Disturbances in Mountain Rivers[END_REF][START_REF] Lisle | The Dominance of Dispersion in the Evoluation of Bed Material Waves in Gravel -Bed Rivers[END_REF][START_REF] Sklar | Translation and dispersion of sediment pulses in flume experiments simulating gravel augmentation below dams[END_REF][START_REF] Sutherland | Evolution of a landslide-induced sediment wave in the Navarro River, California[END_REF]. These studies have focused, however, primarily on the end-member case of how a low amplitude sediment supply compares to the transport capacity of the river. [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF] proposed a 2D morphodynamic approach that examines high-amplitude sediment supplies compared to the transport capacity of the river. In the latter case, the role of dynamic river narrowing in accelerating the removal of landslide-driven sediments is critical. Despite these recent efforts, a full understanding of post-seismic sediment fluxes at the mountain range scale is still lacking.

Post-seismic sediment export is controlled by the sediment supply delivered by the landslide and by the transport capacity of the river receiving the deposit. The quantity of sediment transported by the river is, therefore, likely to be strongly dependent on the degree of connectivity of sources (landslides) to the fluvial network at the initial stage and through time [START_REF] Hovius | Supply and Removal of Sediment in a Landslide Dominated Mountain Belt: Central Range, Taiwan[END_REF]. Several studies have provided a quantification of the initial percentage of earthquake-triggered landslides that connect to the drainage network, ranging from 8% to nearly full connectivity [START_REF] Dadson | Earthquake-triggered increase in sediment delivery from an active mountain belt[END_REF]; Li , 2016;[START_REF] West | Mobilization and transport of coarse woody debris to the oceans triggered by an extreme tropical storm[END_REF]. Work on the temporal evolution of connectivity through time, however, remains an open question [START_REF] Zhang | Evolution of Mass Movements near Epicentre of Wenchuan Earthquake, the First Eight Years[END_REF]. Landslides connected directly to the river network can inject an almost instantaneous sediment load to the river. Landslides deposits that remain on the hillslope are likely to deliver the sediment in a more progressive manner, depending also on post-seismic storms (Fan et al., 2018). Once sediments mobilized by landslides reach the river channel, their export time is expected to depend mostly on the river geometry, discharge and on sediment grain size [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF]. Whereas the evacuation of one landslide has already received attention, no work has been dedicated to the evacuation of seismically-triggered clusters of landslides.

The distributions of landslides can statistically inform on the dynamics of connectivity and subsequent export by river transport. For instance, landslides older than several seismic cycles and persisting in the landscapes have been argued to indicate a low efficiency of sediment export [START_REF] Korup | Large landslides and their effect on sediment flux in South Westland, New Zealand[END_REF]. Such inference, however, cannot be made solely based on individual and old landslides, which represent outliers of the total cluster of landslides triggered by earthquakes or rainfall.

Understanding the triggering and export of landslide clusters over several seismic cycles is required to assess the topographic budget of large earthquakes [START_REF] Parker | Mass wasting triggered by the 2008 Wenchuan earthquake is greater than orogenic growth[END_REF], the role of aftershocks relatively to mainshocks, sediment fluxes at the range scale [START_REF] Hovius | Sediment flux from a mountain belt derived by landslide mapping[END_REF], the geochemical signature of these extreme events [START_REF] Frith | Carbon export from mountain forests enhanced by earthquaketriggered landslides over millennia[END_REF][START_REF] Wang | Earthquake-triggered increase in biospheric carbon export from a mountain belt[END_REF], or the impact and risks associated to these natural hazards [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF][START_REF] Keefer | Earthquake-induced landslides and their effects on alluvial fans[END_REF]. In this study, we develop a nested numerical approach which acts to: 1. simulate earthquakes over several seismic cycles, 2. trigger landslides across catchments, 3. assign the dynamic connectivity to the fluvial network and 4. determine the subsequent sediment transport. Our approach is deliberatively simplified to examine the challenges that emerge when investigating post-seismic sediment evacuation. The nested model integrates sediment export times defined at the reach scale, using the Eros river morphodynamic model (e.g. [START_REF] Davy | A precipiton method to calculate river hydrodynamics, with applications to flood prediction, landscape evolution models, and braiding instabilities[END_REF][START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF], in a statistical model that generates earthquakes and landslides at the mountain range scale, refered to as Quakos. The paper is divided in three sections. First, at the reach-scale, we investigate the impact of one or a series of landslides on the morphodynamic response of a river, the efficiency of sediment export and the persistence of downstream deposits. This allows us to develop a semi-analytical model of landslide export. Second, we focus on embedding the reach-scale model outcomes into the catchment-scale model. We define a reduced-complexity model that accounts for the different processes driving the sediment export of landslide debris that can be applied to clusters of landslides triggered by earthquakes. We apply this model to the case of a hypothetical M w 7.9 earthquake occurring in the Southern Alps of New Zealand. Third, we investigate the morphological impact of a series of earthquakes on landscape dynamics at the mountain range scale. We focus specifically on the roles of the dynamic connectivity of landslide deposits and runoff that modulate river transport capacity. Our results illustrate how the model parametrization impacts the number and volume of triggered landslides, the time persistence, the sediment fluxes leaving the mountain range and the evolution of landslide-size distributions. This leads us to highlight research needs and critical steps required to better understand and predict the impact of large earthquakes on landscape evolution and sediment fluxes.

Morphodynamic modeling of the export of landslide sediment

A C C E P T E D M A N U S C R I P T 6
In this section, we explore the fine-scale dynamics of sediment export of a river that is impacted by a cascade of landslides. We first describe the 2D morphodynamic model Eros that we use to quantify the evacuation of individual landslides at the reach scale [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF]Davy et al., 2017). Following [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF], we present the mechanisms controlling the downstream propagation of a single landslide deposited in a bedrock channel. We then explore the impact of a cascade of landslides on sediment export by introducing several landslides along the same river channel. In addition to revealing the morphodynamic evolution of the landslide deposits, the results from this section will be used in the Quakos study to tackle evacuation of landslide clusters over large areas.

Model description

This study is placed in the context of a bedrock river experiencing a high amplitude sediment forcing that causes perturbations of the river geometry including its width and slope. Therefore, an accurate quantification of landslide removal at the reach scale requires a model that contains the physical processes allowing for the feedbacks between river erosion, transport capacity, flow, geometry and sediment supply. Here, we use Eros [START_REF] Davy | A precipiton method to calculate river hydrodynamics, with applications to flood prediction, landscape evolution models, and braiding instabilities[END_REF], a particle-based model that is well-suited to simulate the evolution of a river subjected to large sediment supplies (Croissant et al, 2017a, b). The particles referred as "precipitons" are elementary volumes of water that move on the top of the topography and interact with it along the downstream path by entraining, transporting or depositing sediment. This model is composed of:

-A hydrodynamic model that predicts water depth and flow velocity patterns on high resolution topographies [START_REF] Davy | A precipiton method to calculate river hydrodynamics, with applications to flood prediction, landscape evolution models, and braiding instabilities[END_REF]. This model resolves the 2D shallow-water equations under the stationary assumption without inertia.

-A vertical and horizontal sediment transport and deposition model that is coupled with the hydrodynamic model. In the following, we only briefly describe the constitutive equations of sediment entrainment, transport and deposition, as a more detailed description can be found Davy et al., [2017]. The rate of sediment entrainment, 𝑒, is defined by the bedload transport law of Meyer-Peter and Muller [1948]:

A C C E P T E D M A N U S C R I P T 7 𝑒̇= 𝐸(𝜏 -𝜏 𝑐 ) 1.5 (1)
with 𝐸 a constant, 𝜏, the shear stress and 𝜏 𝑐 the critical shear stress. The rate of sediment deposition 𝑑 ̇ is a function of the sediment specific discharge 𝑞 𝑠 and transport length, 𝜉 [START_REF] Davy | Fluvial erosion/transport equation of landscape evolution models revisited[END_REF]:

𝑑 ̇= 𝑞 𝑠 𝜉 (2)
In the morphodynamic simulations, 𝜉 is set to 2 m to insure a bedload transport regime where the flow is close to at-capacity conditions in non-supply-limited cases. The model also includes horizontal sediment dynamics. The lateral erosion of neighboring cells is described by:

𝑒̇𝑙 𝑎𝑡 = 𝑘 𝑒 𝑆 𝑦 𝑒̇ (3)
with 𝑘 𝑒 a dimensionless coefficient (here set at 0.05 as in [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF]) and 𝑆 𝑦 the slope in the transverse direction. The lateral sediment deposition 𝑞 𝑠𝑙 is defined as:

𝑞 𝑠𝑙 = 𝑘 𝑑 𝑞 𝑠 𝑆 𝑦 (4)
with 𝑘 𝑑 a constant (here set to 0.5). The choice of parameter values emerges from different calibrations studies (e.g. [START_REF] Davy | A precipiton method to calculate river hydrodynamics, with applications to flood prediction, landscape evolution models, and braiding instabilities[END_REF], Croissant et al, 2017).

The model allows for self-formed channels to emerge and for instabilities, such as bars and braiding, to develop (i.e. the river flow to converge in a self-formed width) when the local sediment flux is non-linearly correlated with the river discharge. In Eros simulations presenting a simple tilted bed as initial topography and no sediment input, Eros shows that the river width that emerges scales with discharge at a power 0.5 which is similar to that measured on natural cases and in flume experiments (e.g. [START_REF] Métivier | Laboratory rivers: Lacey's law, threshold theory, and channel stability[END_REF].

Initial topography and boundary conditions

The model setup is similar to the one used in [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF]. The initial topography is a 3 km long bedrock channel. Its transport capacity is set by its width and slope (Fig. 2a). The water enters through the upstream boundary condition at a constant effective discharge. The landslide deposit is introduced near the upstream end of the bedrock channel. It has a Gaussian shape in the longitudinal direction and is described by its volume, length and median grain size. Bedrock incision

A C C E P T E D M A N U S C R I P T 8
is neglected as we assume that rates of bedrock erosion at short times scale (i.e. a seismic cycle)

would not be large enough to significantly affect the bedrock channel geometry. In the following, we also investigate the evacuation of several landslides introduced along the same channel. In these cases, the bedrock channel lengthens accordingly to accommodate each landslide based on the interlandslide distance that ranges from 0 to 1 km (Fig. 2b).

Sediment evacuation of a single landslide

Here, we investigate the export dynamics of a single landslide deposit in a bedrock river. To fully understand the mechanisms controlling landslide export, we run 80 simulations exploring the parameter space governing the bedrock transport capacity (i.e. width slope, river discharge) and landslides properties (i.e. median grain size, volume) (Fig. 2a). In a previous work, [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF] identified two end-members in terms of landslide evacuation which depends on the ratio between the landslide volume 𝑉 𝑙𝑠 and the river initial transport capacity 𝑄 𝑇 .

-For low 𝑉 𝑙𝑠 /𝑄 𝑇 , the width of the alluvial cover remains equal to the width of the bedrock river and the morphodynamic adaptation of the alluvial cover occurs mainly in slope. The -For high 𝑉 𝑙𝑠 /𝑄 𝑇 , the model predicts an acceleration of the evacuation of a large part of the landslide (50 to 70%) compared to the case where the landslide would be exported at a constant rate. This acceleration is caused by the dynamic narrowing of the alluvial river inside the landslide deposit. The remaining volume of sediment (30-50%) is removed by lateral erosion. This phase is less efficient because the lateral entrainment rate is a fraction of the rate of vertical incision (eq. 3).

Sediment evacuation of a cascade of landslides

After large earthquakes, rivers are likely to receive sediments from multiple co-seismic landslides. Here, we investigate the impact of a cascade of landslides on sediment transport and evacuation. Two scenarios are explored. In the first one, the total volume of sediment introduced in the river remains constant (i.e. 2.10 6 m 3 ) but is distributed along the stream in a cascade of 1 to 4 landslides separated by the same apex to apex distance. In the second one, we investigate the effect of the distance separating two individual landslides on the sediment evacuation dynamics. In both cases, sediment export is evaluated at the outlet of the most downstream landslide. In terms of morphodynamic evolution, the simulations with a cascade of landslide deposits (Fig. 2b and3a) share similarities with those of single deposits. During the first stage, the landslides are large enough to create partial or total landslide dams. The lakes forming in between two deposits have zero transport capacity. Regressive erosion leads to the progressive and simultaneous incision into each landslide deposits by a river narrower than the width of the bedrock channel. Whereas sediments eroded from the most downstream landslide are rapidly evacuated, similarly to the simulation with one landslide, sediments from the upstream landslide enter the downstream lake and form a prograding delta that lasts until the lake is drained. The river eventually incises vertically into both deposits until it reaches the bedrock surface and then removes the remaining volume of sediment by lateral erosion.

When multiple landslide deposits enter the channel, the cases where N ls ≥ 2 display similar export rates (Fig. 3b). Differences arise when these cases are compared with the evacuation of a single large landslide. In the latter, the landslide evacuation is more efficient until ~70% of the deposit is evacuated. The period of river increased efficiency (i.e. when the alluvial river width is less than 40% the one of the bedrock width) is also longer [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF]. We then vary the inter-landslide distance 𝑑 to 0, 250, 500 and 1000 m. The results show that the distance between successive landslides appears to play only a minor role in the efficiency of total sediment evacuation. When we examine landslides individually (by setting gauging stations downstream of each deposit), however, the downstream landslides are evacuated less efficiently than those upstream, especially during the phase of lateral erosion (see Supplementary Figure S1). This is because the high sediment fluxes from the evacuation of the upstream landslides slows down the erosion of those downstream as the river is already loaded with sediment. Apart from this effect, the time-evolution of sediment transport for individual landslides remains relatively similar when considering a cascade of landslides or a single landslide. Taken together, when viewed from the downstream gauging station, the overall sediment

A C C E P T E D M A N U S C R I P T
11 export response of a cascade of landslides is similar to that of a single landslide in these single channel experiments (Fig. 3).

We, therefore, assume in the following discussions that the removal of landslides can be considered independently from each other. This assumption represents a critical step in building an upscaling approach where thousands of landslides can be triggered simultaneously in the river stream.

A reduced-complexity model for the export of landslide sediment

Earthquakes generally trigger numerous landslides that will eventually reach the fluvial network at locations with varying local transport capacity. As such, a cluster of landslides can be described as a distribution of the ratio 𝑉 𝑙𝑠 /𝑄 𝑇 [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF] spanning the studied endmembers in the previous section. Whereas morphodynamic modelling provides useful information on the mechanisms of landslide removal at the reach scale, it is still too computationally demanding to apply it at the scale of a whole mountain range over hundreds to thousands of years and to accounting for thousands or more landslides. Therefore, we aim here to define a reduced-complexity model to describe the post-seismic evolution of landslide volume for any value of 𝑉 𝑙𝑠 /𝑄 𝑇 , informed by the results obtained with 80 simulations performed in the previous section for the single landslide scenario.

Using these simulations, we compute the time 𝑇 𝑒𝑥𝑝 needed to export 20 to 90% of the initial 𝑉 𝑙𝑠 as a function of 𝑉 𝑙𝑠 /𝑄 𝑇 at a 10% percent interval (Fig. 4a). As described in [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF] the export time computed for each interval follows a trend with 𝑉 𝑙𝑠 /𝑄 𝑇 that can be fitted with a function of the form:

𝑇 𝑒𝑥𝑝,𝑖 = 𝛿 𝑖 (𝜇 𝑖 𝑉 𝑙𝑠 𝑄 𝑇 ) [1 + ( 𝑉 𝑙𝑠 𝑄 𝑇 ) 𝜑 𝑖 ] (𝛽-1) 𝜑 𝑖 (5)
with 𝛿 and 𝜇 constants, 𝜑 a curvature parameter, 𝛽 an exponent (here fixed at 0.1) and 𝑖 corresponds to the studied percent interval. The values of these parameters are found by a least-square fitting using equation 5 on each percent interval (Fig. 4a, Table 1).

A C C E P T E D M A N U S C R I P T

12

We obtain a discrete description of the evolution of a landslide volume for any value of 𝑉 𝑙𝑠 /𝑄 𝑇 (Fig. 4b,c). The continuous description is obtained by interpolating linearly between the points.

Following this, the temporal evolution of landslide sediment evacuation can be assessed for two values of 𝑉 𝑙𝑠 /𝑄 𝑇 which represent two end-members responses (Fig. 4b,c). Simulations presenting a high 𝑉 𝑙𝑠 /𝑄 𝑇 ratio never succeed to evacuate 100% of the initial volume of sediment as a small fraction of sediment remains captured in lateral terraces. Therefore, to reconstruct fully the landslide export, we assume that the last 10% of sediment volume is exported at the same rate as that estimated for the last 20% to 10%. This may lead to a slight over-estimation of the efficiency of sediment transport during the lateral erosion phase.

This reduced-complexity method provides a computationally efficient way to appraise the post-seismic sediment evolution of landslide clusters. It also has the advantage of implicitly accounting for the evolution of the width and slope of the river. This reduced-complexity description of sediment export can, therefore, be applied at a larger spatial scale (i.e. a mountain range) provided that the distribution of 𝑉 𝑙𝑠 /𝑄 𝑇 is known. In the next section, we combine this method with a newly developed large-scale model, Quakos, which determines clusters of earthquake-triggered landslides, including the volume and the transport capacity of the river that they are connected to.

A C C E P T E D M A N U S C R I P T

Upscaling to sediment fluxes at the mountain range scale during several seismic cycles

Here, we propose a statistical approach to quantify post-seismic sediment fluxes at the scale of a mountain range. Quakos is composed of three main components that can be used independently (Fig. 5a): 1. a statistical earthquake generator; 2. a landslide generator that predicts the 2D distribution of co-seismic landslides and the volumes using empirical laws; 3. a sediment evacuation model based on the reduced-complexity method described in previous section. The values of the parameters used in this section can be found in Table 2.

Model description

Study area

In the following, we consider a hypothetical failure of the Alpine Fault in the Southern Alps, New Zealand as a case study (Fig. 5). The formation of the Southern Alps and its recent morphologic evolution has occurred under extreme tectonic and climatic forcing [START_REF] Tippett | Fission track analysis of the Late Cenozoic vertical kinematics of continental pacific crust, South Island, New Zealand[END_REF]. The rapid rate of convergence is largely accommodated by the Alpine Fault [START_REF] Norris | The obliquely-convergent plate boundary in the South Island of New Zealand : implications for ancient collision zones[END_REF] which paleo-seismic studies have shown ruptures in large earthquakes (M w > 7.5) with a recurrence time of 263 ± 68 years [START_REF] Berryman | Major Earthquakes Occur Regularly on an Isolated Plate Boundary Fault[END_REF][START_REF] Cochran | A plate boundary earthquake record from a wetland adjacent to the Alpine fault in New Zealand refines hazard estimates[END_REF][START_REF] Howarth | Lake sediments record cycles of sediment flux driven by large earthquakes on the Alpine fault[END_REF]. The` range of the Southern Alps extends ~450 km from southwest to northeast and rises from sea level to 3724 m at Mount Cook. It forms a natural barrier to western winds which leads to high rates of precipitation up to 13 m.yr -1 along the west coast [START_REF] Tait | Analysis of the Spatial Interpolation Error associated with Maps of Median Annual Climate Variables[END_REF]. Landscapes in the Southern Alps are characterized by steep hillslopes with modal values averaging at 35° [START_REF] Korup | The role of landslides in mountain range evolution[END_REF] and are prone to landsliding even during aseismic periods [START_REF] Hovius | Sediment flux from a mountain belt derived by landslide mapping[END_REF]. Thousands of landslides are expected to be triggered in the next large seismic event, potentially mobilizing ~1 km 3 of sediment (Marc et al., 2016;[START_REF] Robinson | Coseismic landsliding estimates for an Alpine Fault earthquake and the consequences for erosion of the Southern Alps, New Zealand[END_REF].

The topography of the Southern Alps is obtained from the SRTM3 digital elevation model (DEM), with a pixel size set at 100 m. Computations involving operations on the DEM have been performed using TopoToolbox_v2 [START_REF] Schwanghart | Short Communication: TopoToolbox 2 -MATLAB-based software for topographic analysis and modeling in Earth surface sciences[END_REF].

A C C E P T E D M A N U S C R I P T

Fault, earthquakes and peak ground acceleration

Here, we design a thrust fault of length F L = 400 km and width F W = 19 km with a dipping angle of 60° that approximatively mimics the geometry of the Alpine Fault [START_REF] Robinson | Coseismic landsliding estimates for an Alpine Fault earthquake and the consequences for erosion of the Southern Alps, New Zealand[END_REF]. In the Quakos model, earthquakes (mainshocks and aftershocks) are generated on this fault. When mainshocks are generated the ruptures cover the entire fault width, leading to earthquakes of magnitude M w = 7.9. The position of each mainshock is randomly sampled along the fault plane. Each mainshock triggers a series of aftershocks, which location, date and magnitude are determined using the BASS model [START_REF] Turcotte | BASS, an alternative to ETAS[END_REF]. The aftershocks series are generated following three statistical laws:

-The difference between the magnitude of the mainshock and its largest associated aftershock (∆M w = 1.25) is determined using a modified version of Båth's law [START_REF] Shcherbakov | A generalized Omori's law for earthquake aftershock decay[END_REF].

-The rate of aftershocks is submitted to a temporal decay described by a generalized form of Omori's law [START_REF] Shcherbakov | A generalized Omori's law for earthquake aftershock decay[END_REF]. Parameters from this law are the exponent p = 1.25 and offset c = 0.1 days.

-The spatial distribution of aftershocks is given by a spatial form of Omori's law [START_REF] Helmstetter | Predictability in the Epidemic-Type Aftershock Sequence model of interacting triggered seismicity[END_REF] with the exponent q = 1.35 and offset d = 4 meters.

All the defined parameter values of the BASS model are taken from [START_REF] Turcotte | BASS, an alternative to ETAS[END_REF] and are constant for all the simulations performed in this paper.

The range of simulated magnitudes is bounded by fault dimensions, for the largest magnitude, and by the spatial discretization of the fault, for the smallest magnitude. Simulated earthquakes have magnitudes ranging from 2.5 to 7.9. For each earthquake (mainshocks and aftershocks), the associated rupture length (𝑅 𝐿 ) and width (𝑅 𝑊 ) are estimated as a function of seismic moment, 𝑀 𝑂 following a consistent set of scaling laws based on empirical observations determined for strike-slip faults [START_REF] Leonard | Earthquake fault scaling: Self-consistent relating of rupture length, width, average displacement, and moment release[END_REF]:

𝑅 𝐿 = ( 𝑀 𝑂 𝜇𝐶 1 3 2 ⁄ 𝐶 2 ) 𝛽 (6) A C C E P T E D M A N U S C R I P T 15 𝑅 𝑊 = 𝐶 1 𝑅 𝐿 𝛿 (7)
where 𝜇 = 33 GPa is the shear modulus and 𝐶 2 = 3.6.10 -5 , 𝐶 1 , 𝛽 and 𝛿 are constants which value depends on the rupture length:

-𝑅 𝐿 < 5 km: 𝐶 1 = 1, 𝛽 = 1/3 and 𝛿 = 1.

-5 km < 𝑅 𝐿 < 45 km: 𝐶 1 = 15, 𝛽 = 2/5 and 𝛿 = 2/3.

-𝑅 𝐿 > 45 km: 𝐶 1 = F W , 𝛽 = 2/3 and 𝛿 = 1.

For a given earthquake, a map of the synthetic peak ground acceleration (PGA) is computed as a function of earthquake magnitude, ruptured fault geometry (rake, dip, dimension), fault mechanism, lithological controls and site effects (Fig. 5b) [START_REF] Campbell | NGA ground motion model for the geometric mean horizontal component of PGA, PGV, PGD and 5% damped linear elastic response spectra for periods ranging from 0.01 to 10 s[END_REF]. The theoretical framework that derives from this work is quite extensive (see equations 1 to 12 in [START_REF] Campbell | NGA ground motion model for the geometric mean horizontal component of PGA, PGV, PGD and 5% damped linear elastic response spectra for periods ranging from 0.01 to 10 s[END_REF]). Here, we consider the sediment depth Z 2.5 = 0 m, a S-wave velocity in the first 30 m of the crust V s,30 = 180 m.s -1 and a reference PGA at V s = 1100 m.s -1 , A 1100 = 0.10 g [START_REF] Robinson | Coseismic landsliding estimates for an Alpine Fault earthquake and the consequences for erosion of the Southern Alps, New Zealand[END_REF]. For other parameter values, we refer to Table 2 in [START_REF] Campbell | NGA ground motion model for the geometric mean horizontal component of PGA, PGV, PGD and 5% damped linear elastic response spectra for periods ranging from 0.01 to 10 s[END_REF].

Landslides triggered by earthquakes

Modelled maps of PGA are then used to infer the spatial density of triggered landslides, i.e.

the number of landslides by unit area. Previous studies suggest that the density of earthquaketriggered landslides is linearly dependent on the PGA [START_REF] Meunier | Regional patterns of earthquake-triggered landslides and their relation to ground motion[END_REF][START_REF] Yuan | Density distribution of landslides triggered by the 2008 Wenchuan earthquake and their relationships to peak ground acceleration[END_REF]:

A C C E P T E D M A N U S C R I P T 𝑃 𝑙𝑠 = 𝛼 𝑝 𝑃𝐺𝐴 -𝛽 𝑝 (8)
where 𝛼 𝑝 and 𝛽 𝑝 , are two empirical parameters that controls the maximum density of landslide and its spatial repartition, respectively. The parameter 𝛽 𝑝 is a critical value of PGA under which no or very few landslides are triggered. The values of these parameters are highly dependent on the studied case and the choice of their value is explained in the result section. Following the work of [START_REF] Meunier | Regional patterns of earthquake-triggered landslides and their relation to ground motion[END_REF], locations where the local slope is less than 20% (on a 100m DEM) are not affected by landsliding. In our case study, this mostly prevents landslides from being triggered on alluvial fans and large river valleys (Fig. 5c,d).

The area distribution of earthquake-triggered landslides is commonly given by an inverse gamma probability density function:

𝑝𝑑𝑓(𝐴 𝑙𝑠 ) = 1 𝑎𝛤(𝑎) [ 𝑎 𝐴 𝑙𝑠 -𝑠 ] 𝜌+1 𝑒𝑥𝑝 [- 𝑎 𝐴 𝑙𝑠 -𝑠 ] (9)
with 𝑎 a parameter controlling the position of the pdf maximum, 𝑠 a parameter controlling the rollover for small landslides and 𝜌 is a positive exponent controlling the tail of the pdf (Malamud et al., 2004a). In Quakos simulations, the area of each landslide belonging to a landslide cluster is determined by randomly sampling the 𝑝𝑑𝑓(𝐴 𝑙𝑠 ). The location of landslides is then determined according to the map of landslide density (Fig. 5d). At this stage, those two conditions are met and then landslides are distributed randomly on slopes above >20%. Future improvements to this approach may include methods to preferentially select parts of the landscape where landslides locate, and link this to the morphology of the local hillslopes at the ~ km 2 scale.

For each individual predicted landslide, its planform area is converted to volume 𝑉 𝑙𝑠 by using an empirical scaling law:

𝑉 𝑙𝑠 = 𝛼𝐴 𝑙𝑠 𝛾 ( 10
)
with 𝛼 and 𝛾 are set to 0.05 and 1.5 [START_REF] Hovius | Sediment flux from a mountain belt derived by landslide mapping[END_REF][START_REF] Larsen | Landslide erosion controlled by hillslope material[END_REF]. This parameterization is well-suited to infer the volume of deep-seated bedrock landslides which dominate the volume budget of a population of triggered-landslides.

Post-seismic evacuation of landslide sediment

A C C E P T E D M A N U S C R I P T

18 After Quakos has provided a distribution of earthquake-triggered landslides, the prediction of post-seismic landslide evacuation depends mostly on: 1. the rate of sediment supply to the channel network, i.e. the transfer of material from hillslopes to channels, and 2. the rate of sediment transport by the river.

Several studies have pointed out the importance of the initial and dynamic connectivity of landslides to the channel network on post-seismic sediment fluxes [START_REF] Li | Connectivity of earthquake-triggered landslides with the fluvial network: Implications for landslide sediment transport after the 2008 Wenchuan earthquake[END_REF][START_REF] Roback | The size, distribution, and mobility of landslides caused by the 2015 M w 7.8 Gorkha earthquake, Nepal[END_REF]. Some studies provide an estimation of initial connectivity ranging from 8% to 43 % and even to full connectivity [START_REF] Dadson | Earthquake-triggered increase in sediment delivery from an active mountain belt[END_REF][START_REF] Li | Connectivity of earthquake-triggered landslides with the fluvial network: Implications for landslide sediment transport after the 2008 Wenchuan earthquake[END_REF]; A. J. [START_REF] West | Mobilization and transport of coarse woody debris to the oceans triggered by an extreme tropical storm[END_REF]. Based on the data of [START_REF] Li | Connectivity of earthquake-triggered landslides with the fluvial network: Implications for landslide sediment transport after the 2008 Wenchuan earthquake[END_REF] and [START_REF] Roback | The size, distribution, and mobility of landslides caused by the 2015 M w 7.8 Gorkha earthquake, Nepal[END_REF], the initial landslide-channel connectivity 𝐶 of each landslide is determined as a function of its area 𝐴 𝑙𝑠 :

𝐶 = 𝑚𝐴 𝑙𝑠,𝑏𝑖𝑛 𝜔 ( 11
)
where 𝑚 is an empirical constant and 𝜔 an empirical exponent (see Supplementary Figure S3).

Equation ( 11) applies for landslides presenting an area lower than 1.10 6 m 2 . Above this threshold, we assume that landslides are always initially connected to the drainage network. This assumption is supported by empirical data and because larger landslides usually present a longer run out [START_REF] Lucas | Frictional velocity-weakening in landslides on Earth and on other planetary bodies[END_REF]. Here, 𝐶 gives the percentage of connected landslides in the considered bin of landslide area (𝐴 𝑙𝑠,𝑏𝑖𝑛 ) on a logarithmic scale.

To consider landslide deposits which locate away from the channel network on steep slopes, we consider that transport of loose debris is likely to occur at a velocity that likely depends on the climatic and meteorological context and local topographic properties. Unfortunately, studies which quantify the post-seismic sediment delivery from hillslopes to channels are scarce because of the difficulties of measuring it in the field or using remote sensing (Fan et al., 2018;[START_REF] Zhang | Evolution of Mass Movements near Epicentre of Wenchuan Earthquake, the First Eight Years[END_REF].

We, therefore, develop a simplified approach to account for a "dynamic connectivity" of landslides to rivers. This approach first determines the distance (𝑑) between the landslide and the closest river connection point using a steepest descent algorithm. The timing of connection is then obtained by setting a constant and arbitrary connectivity velocity (𝑢 𝑐𝑜𝑛 ) to each landslide and computed as 𝑡 𝑐𝑜𝑛 = 𝑑/𝑢 𝑐𝑜𝑛 . Once the landslide is connected to the river network, we assume that the whole volume is connected, i.e. completely available to be removed by the river. A future

A C C E P T E D M A N U S C R I P T

modification may be to link this dynamic connectivity to diffusive erosion processes which occur on steep surfaces (Roering et al., 1999). Here, for simplicity, we first explore scenarios with a constant 𝑢 𝑐𝑜𝑛 .

Once a landslide reaches the closest stream, its subsequent evacuation depends on the ratio between its volume 𝑉 𝑙𝑠 and the local river transport capacity 𝑄 𝑇 following equation ( 5). If the landslide volume is determined for each landslide using Quakos, the transport capacity needs to be computed. The along-stream transport capacity of bedrock rivers is set by its geometry (width and slope), local river discharge and sediment grain size. The bedrock river width (𝑊), slope (𝑆) and mean discharge (𝑄 ̅ ) are expressed as a function of the local drainage area (𝐴) as:

{ 𝑊 = 𝑘 𝑤𝑛 𝐴 0.5 𝑆 = 𝑘 𝑠𝑛 𝐴 -0.45 𝑄 ̅ = 𝑟̅ 𝐴 (12)
with 𝑘 𝑤𝑛 the normalized width index, 𝑘 𝑠𝑛 the normalized steepness index and 𝑟̅ mean annual runoff [START_REF] Lague | The stream power river incision model: Evidence, theory and beyond[END_REF]. Here, the critical drainage area used to extract the drainage network is equal to 0.5 km 2 .

To ensure the same context as the Eros simulations, the river transport capacity is described using an effective daily discharge (Q eff ) presenting a return time of one year which is a good compromise between frequency of occurrence and the amount of geomorphic work of such events.

The bedload transport capacity [START_REF] Meyer-Peter | Formulas for Bed-Load Transport[END_REF] is then computed as:

𝑄 𝑇 = 𝑊𝐾 (𝜌 𝑤 𝑔 ( 𝑛𝑄 𝑒𝑓𝑓 𝑊 ) 0.6 𝑆 0.7 -(𝜌 𝑠 -𝜌 𝑤 )𝑔𝜏 𝑐 * 𝐷 50 ) 1.5 (13)
with 𝜏 𝑐 * the critical Shields stress, 𝜌 𝑤 the water density, 𝜌 𝑠 the sediment density, 𝑛 the Manning friction coefficient, 𝐾 an erodability constant and 𝑔 the gravitational constant. For simplicity, the grain size distribution of the landslide is reduced to the median grain size descriptor 𝐷 50 which is chosen to be constant for all landslides.

To compute the value of Q eff , we assume that the range of daily discharges 𝑄, experienced at any point along the river, follows an inverse-gamma probability density function:

𝑝𝑑𝑓(𝑄) = 𝑘 𝑘+1 𝛤(𝑘 + 1) 𝑒𝑥𝑝 (- 𝑘 𝑄/𝑄 ̅ ) (𝑄/𝑄 ̅ ) -(2+𝑘) (14) A C C E P T E D M A N U S C R I P T 20
with 𝛤 the gamma function and 𝑘 a parameter linked to the variability of the hydrological forcing, here, fixed at k = 1, based on empirical data [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF][START_REF] Lague | Discharge, discharge variability, and the bedrock channel profile[END_REF]. This assumption is supported by empirical data, that demonstrate that the runoff of rivers located along the West Coast of New Zealand present a high variability (k = 1) (Croissant et al, 2017b). From this distribution, the return time (𝑡 𝑟 ) of a particular daily discharge can be assessed using:

𝑡 𝑟 (𝑄 𝑒𝑓𝑓 ) = 𝛤(𝑘/𝑄 𝑒𝑓𝑓 , 𝑘 + 1) -1 (15)
The value of Q eff can be computed using this equation.

Landslide triggering and sediment export over a seismic cycle 4.2.1 The volume of triggered landslides by earthquakes on the Alpine Fault

To assess Quakos outputs, we compare them to natural observations and existing analytical models [START_REF] Keefer | Earthquake-induced landslides and their effects on alluvial fans[END_REF]Odin Marc et al., 2016). We plot the total volume of landslide clusters as a function of earthquake magnitude (Fig. 7). Consistent with Marc et al, [2016], Quakos leads to a threshold magnitude, here ~4.5, under which no landslides are generated. This threshold results from the fact that landslides are only generated if a critical PGA is reached (eq. 6; see Supplementary Movie S1). Above this threshold, the total landslide volume triggered by an earthquake increases with magnitude and shows a sensitivity to the depth of earthquake nucleation. The variability of total landslide volume for earthquakes of equal magnitude results from the depth of the earthquake but also from the variability of the topography impacted by landsliding, including the proportion of the topography with local slopes greater than 20%. For low magnitudes, total landslide volume is strongly sensitive to the depth of the earthquake. Indeed, the width extent of the rupture is small compared to fault width, and the depth of the earthquake becomes the controlling factor to generate PGA above the critical value for landsliding. This mostly explains the spread in the distribution of total landslide volume as a function of magnitude for earthquakes with magnitude lower than ~6.5. Larger earthquakes are less prone to this spread in total landslide volume as the rupture width becomes closer to fault width. We also note that Quakos outputs asymptotically tend, for large magnitudes, towards the analytical model from Marc et al, [2016] for R 0 = 15 km, independently from Quakos depth of the earthquake. This results from the PGA model in Quakos which depends on rupture extent (Campbell 2008), covering the entire fault width for the largest earthquakes, and not on the depth of the earthquake. Whereas Marc et al. [2016] only consider the depth of the sources of seismic waves, and not the rupture extent. Total landslide volume is also sensitive to the parameters in equations 8 and 9, which control the number of triggered landslides for a given PGA and their volume distribution (see supplementary Figure S2).

A C C E P T E D M A N U S C R I P T 21 & Bozorgnia,
For a M w 7.9 earthquake occurring on the Alpine Fault, without considering aftershocks, the total number of landslides generated in the Southern Alps is ~17700 ± 500 landslides with a total volume of sediment 0.75 ± 0.07 km 3 (Fig. 6). Empirical constraints allow for a more accurate estimation of the characteristics of a future earthquake on the Alpine Fault.The total number of landslides, however, is of the same order of magnitude as most of the natural cases documented in [START_REF] Tanyaş | Presentation and Analysis of a Worldwide Database of Earthquake-Induced Landslide Inventories[END_REF] and the total volume is similar to the one estimated by two independent studies (Odin Marc et al., 2016; T. R. [START_REF] Robinson | Coseismic landsliding estimates for an Alpine Fault earthquake and the consequences for erosion of the Southern Alps, New Zealand[END_REF].

Sediment export over a seismic cycle

A C C E P T E D M A N U S C R I P T

22 Here, we explore the dynamic of sediment export over a single seismic cycle that follows a M w 7.9 earthquake. We only consider the river catchments of the West Coast where the landslide number is ~12000 and total volume is ~0.5 km 3 . We focus on the role of the initial and dynamic connectivity of landslides to the fluvial network in controlling post-seismic sediment evacuation.

Based on the morphodynamic modelling results, we assume that landslides located in the same river reach are evacuated independently of each other.

At the initial stage, ~41% of the total volume of landslide sediment is connected to the drainage network. Models with 𝑢 𝑐𝑜𝑛 equals to 0.1, 1 and 10 m.yr -1 lead to 43, 50 and 96% of the total volume evacuated over 263 years. The 'full connectivity' case leads to the highest rates of sediment transport with 70% of the landslide mass evacuated in less than 10 years (Fig. 7a). This value matches to order of magnitude the predictions of [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF] for high mean annual runoff and runoff variability and for which the full population of landslides was assumed fully connected. After 2000 years, most of the sediment volume is removed for 𝑢 𝑐𝑜𝑛 ≥ 1 m.yr -1 , whereas the rivers are starved of sediments because of the absence of new connected landslides for 𝑢 𝑐𝑜𝑛 = 0.1 m.yr -1 . This illustrates, that landslide connectivity is an important factor limiting sediment evacuation after the triggering of landslides by a large earthquake. In turn, landslide connectivity is critical to assess the fraction of co-seismic debris evacuated before the next earthquake occurs. Interestingly, the number of active landslides (Fig. 7d), i.e. connected landslides with remaining sediments, do not exactly follow the total landslide volume evolution (Fig. 7a). This occurs because the volume of deposits in the landscape is controlled by the largest landslides. For low values of 𝑢 𝑐𝑜𝑛 , model predictions show that a large proportion of the initial landslide population can be preserved while having evacuated a moderate to large proportion of the volume of earthquakeproduced sediment (Fig. 7b).

Connectivity also impacts the amplitude and duration of the rates of sediment remobilization (Fig. 7c). The full connectivity case presents rates that are at least one order of magnitude greater than any other model during the first years after the earthquake. This rate drops abruptly by three orders of magnitude in less than a century. On the contrary, for 𝑢 𝑐𝑜𝑛 < 10 m.yr -1 , the the rate of remobilization oscillates during the first 200 years before decreasing progressively. The different rates are controlled by the sediment delivery from hillslope to the channels (Fig. 7d). A low value of 𝑢 𝑐𝑜𝑛 ensures a The velocity of connection has also an impact of the evolution of the distribution of landslide volume (Fig. 8). We here account for the volume changes as landslide are being evacuated by river sediment export. For the full connectivity case, the temporal evolution of the probability density function (pdf) of landslide volume, pdf(V ls ), shows that landslides with volumes smaller than 10 3 m 3 disappear from the distribution after only five years as they tend to have characteristic timescales 𝑉 𝑙𝑠 /𝑄 𝑡 < 1 year. Over longer durations, the remaining landslides tend towards the largest areas and the distribution shrinks towards these largest areas. There is however no change in the scaling of the tail of the pdf(V ls ). For 𝑢 𝑐𝑜𝑛 = 10 m.yr -1 the shape of the pdf is preserved during 50 years until the largest landslides start to dominate the long-term signal (Fig. 8b). Moreover, the slope of the tail of the pdf changes with time and becomes less steep. This occurs because of the parametrization of the initial connectivity in Equation ( 11) which favors the connection of large landslides and tend to preserve small landslides. For 𝑢 𝑐𝑜𝑛 ≤ 1 m.yr -1 , and probably also for lower values of 𝑢 𝑐𝑜𝑛 , the shape of the pdf is preserved through time. This occurs because only a very limited fraction of the landslide population is actively connected to the fluvial network, and most landslides are preserved within the mountain range.

Upscaling to several seismic cycles

In this section, we extend model duration to several seismic cycles. The scenario is chosen to mimic the response observed on the Alpine fault, i.e. a temporal series of 12 M w = 7.9 mainshocks separated by recurrence period randomly sampled in the 263 +/-68 years range (Fig. 9). The production of co-seismic landslides is slightly different for each mainshocks, because ofthe stochastic way in which Quakos simulates individual landslide areas and volumes, with a total volume ~0.45 m 3 and an average total number of landslides of 13000 on the West Coast.

Each mainshock is followed by a series of aftershocks with magnitudes varying between 2.5 and 7.5, some of which mobilise an additional volume of sediment. In most cases, the contribution of aftershocks is generally lower than that of mainshocks for two reasons: 1. they mobilize sediment volumes that are lower by 1 to 5 orders of magnitude (Fig. 9a) and 2. most of landslide-triggering aftershocks are quasi-synchronous with mainshocks and, therefore, the total sediment production is dominated by the one of the mainshock. Aftershocks occurring between two mainshocks, however, can have a visible impact on the sediment production as highlighted by the M w 7.5 earthquake at ~2000 years that mobilize sediment volumes comparable to the mainshocks. Over the 12 seismic cycles, rivers are never able to export all the earthquake-mobilized sediment (Fig. 9c). Yet, they progressively reach a dynamic steady-state between new sources of sediment, because of landslide triggering during earthquakes, and sediment evacuation of all the previous generation of co-seismic landslides. The duration of the transient phase and the average volume of sediment remaining in the mountain at steady-state is controlled by 𝑢 𝑐𝑜𝑛 . For 𝑢 𝑐𝑜𝑛 ≥ 10 m.yr -1 the steady-state is reached after three seismic cycles with a sediment storage within the mountain range reaching a minimum of 10-20% of one M w 7.9 earthquake worth of sediment. Most of the landslides are evacuated within one seismic cycle reaching a minimum of 50-100 landslides remaining in the catchments (Fig. 9d). For 𝑢 𝑐𝑜𝑛 ≤ 1 m.yr -1 , the earthquake-produced sediment progressively piles-up inside the mountain range, stored on hillslopes, until an equilibrium situation is reached after 1500 years in which only a small proportion of the newly triggered population is removed. In these cases, the sediment volume that is stored within the mountain range at the end of Figure 10 shows the repartition of landslides generation at the occurrence of each mainshocks.

A C C E P T E D M A N U S C R I P T

The repartition is dominated, unsurprisingly, by the i-th generation of landslides at the time of occurrence of the i-th mainshock. Despite that, the full connectivity case (Fig. 10a) shows that only a few landslides of early generations are preserved over tens of seismic cycles. This is illustrated by the presence of two landslides of the 1-th generation at the occurrence of the 12-th mainshock. Because no dynamic connectivity occurs in this case, this persistence is explained by rare landslides characterized by very low transport capacity and very high 𝑉 𝑙𝑠 /𝑄 𝑡 . For low velocities of connection (Fig. 10b), preservation of landslides over tens of seismic cycles becomes less exceptional. This is illustrated by the wide diversity of generations of landslides preserved when the last mainshock occurs. Based on the results of Figure 9, we aim to define a representative export time of the sediment delivered by the whole population of earthquake-triggered landslides. This index is based on the two intrinsic timescales that are involved in these simulations: 1. the time, 𝑡 ℎ𝑖𝑙𝑙 , necessary to transfer landslide material from hillslopes to channels by dynamic connectivity and 2. the time, 𝑡 𝑟𝑖𝑣 , necessary to transport landslide sediment by the river network. 𝑡 𝑟𝑖𝑣 is computed as the time necessary to remove 90% of the total landslide volume for different mean annual runoff intensity under the assumption that all the landslides are connected. 𝑡 ℎ𝑖𝑙𝑙 is computed as the mean connection time of the whole landslide population:

A C C E P T E D M A N U S C R I P T

𝑡 ℎ𝑖𝑙𝑙 = ∑ 𝑉 𝑙𝑠,𝑗 𝑡 𝑐𝑜𝑛,𝑗 𝑁 𝑙𝑠 𝑗=1 ∑ 𝑉 𝑙𝑠,𝑗 𝑁 𝑙𝑠 𝑗=1 (16)
With 𝑁 𝑙𝑠 the number of landslides, 𝑉 𝑙𝑠,𝑗 and 𝑡 𝑐𝑜𝑛,𝑗 the volume and connection time of the j-th landslide. We then compute the mean export time 𝑡 𝑒𝑥𝑝 of the landslide population normalized by the average duration of the seismic cycle 𝑇 𝑠-𝑐 :

𝑡 𝑒𝑥𝑝 = 𝑡 ℎ𝑖𝑙𝑙 + 𝑡 𝑟𝑖𝑣 𝑇 𝑠-𝑐 A C C E P T E D M A N U S C R I P T 29
The mean export time of the landslide population increases when decreasing the mean annual runoff and or the connection velocity 𝑢 𝑐𝑜𝑛 (Fig. 11a). 𝑡 𝑒𝑥𝑝 is lower than one seismic cycle for 𝑢 𝑐𝑜𝑛 greater than 2 m.yr -1 , except when the mean annual runoff becomes close to 2 m.yr -1 , the critical runoff below which no sediment transport occurs. For a mean annual runoff greater than 3 m.yr -1 and 𝑢 𝑐𝑜𝑛 lower than 2 m.yr -1 , 𝑇 𝑒𝑥𝑝 is almost independent of the mean annual runoff. For this regime sediment export is limited by connectivity as 𝑇 ℎ𝑖𝑙𝑙 is greater than 𝑇 𝑟𝑖𝑣 and the rivers are lacking sediments (Fig. 11b). This highlights the importance of hillslope processes in controlling the export time of a landslide population. When 𝑢 𝑐𝑜𝑛 is greater than 10 m.yr -1 , sediment export is limited by river transport as 𝑇 ℎ𝑖𝑙𝑙 is lower than 𝑇 𝑟𝑖𝑣 and the rivers are loaded by sediments. The limit between the connectivity-and transport-limited regimes tends towards lower values of 𝑢 𝑐𝑜𝑛 when the mean annual runoff gets closer to its critical value around 2 m.yr -1 .

Discussion

Our theoretical modelling approach emphasizes the need to integrate various processes involved in the production and routing of co-seismic landslide debris across several seismic cycles, partly at the expense of individual process complexity. An important novelty is the spatially explicit nature of peak ground acceleration, co-seismic landslide and transport, and the inclusion of important channel morphodynamic feedback occurring for large landslides. Yet, our model is based on several assumptions. We still have a rather limited understanding and/or lack of empirical data for three key processes: generation of co-seismic landslide, landslide connectivity and hillslope transfer, and river morphodynamics. Here, we discuss these limitations and highlight challenges to overcome.

Generation of co-seismic landslides

Compilations of empirical datasets show that co-seismic landslide cluster properties (total volume, area and number) exhibit large variations which mainly depend on the tectonic forcing and the topographic properties of the affected mountain range [START_REF] Tanyaş | Presentation and Analysis of a Worldwide Database of Earthquake-Induced Landslide Inventories[END_REF]. Early work has shown that these properties are primarily correlated to earthquake magnitude [START_REF] Keefer | The importance of earthquake-induced landslides to long-term slope erosion and slope-failure hazards in seismically active regions[END_REF]. Other studies

A C C E P T E D M A N U S C R I P T 30
have improved this analysis by accounting for seismological characteristics and important topographic attributes to provide a seismologically-consistent model (Marc et al., 2016). No theoretical framework, however, is robust enough to predict with accuracy what would be the properties (volume and localisation) of a landslide cluster for a future event. In our modelling approach, we build on the previous work by using a PGA-dependent model that also accounts for topographic properties to generate co-seismic landslide clusters. The novelty of this approach relies on predictions of spatialized landslide distributions within a mountain range while being consistent with previous empirical data and analytical model. An accurate prediction of landslide cluster properties is important because several components of the post-seismic sediment evacuation depend on them, i.e.

the dependency of initial connectivity to landslide area and export time to landslide volume (see Supplementary Figure S2 andS4).

In our approach, Quakos generates earthquake-triggered landslide clusters with a spatial density that linearly depends on PGA [START_REF] Meunier | Regional patterns of earthquake-triggered landslides and their relation to ground motion[END_REF]. Slope, curvature, and distance from ridges are identified as other controlling parameters of landslide initiation [START_REF] Meunier | Topographic site effects and the location of earthquake induced landslides[END_REF][START_REF] Robinson | Rapid post-earthquake modelling of coseismic landslide intensity and distribution for emergency response decision support[END_REF][START_REF] Tanyas | A global slope unit-based method for the near real-time prediction of earthquake-induced landslides[END_REF]). Yet, to date, no consensus exists on a theoretical framework to accurately predict the spatial distribution of earthquake-generated landslides [START_REF] Reichenbach | A review of statistically-based landslide susceptibility models[END_REF]. A common theme of these assessments is that they rely primarily on the PGA patterns and have to be calibrated against a pre-existing partial co-seismic landslide catalogues.

These are non-existent for the Alpine Fault. Accounting for more detailed model of landslide susceptibility would have an impact on the landslide distances to stream distribution and then on the dynamic connectivity. In addition, in our modelling approach the maximum value of the landslide density is limited by the transport capacity computed at the outlet of catchments. Indeed, the sum of sediment volume exported out of connected landslides within one catchment cannot be greater that the transport capacity at the catchment outlet.

We have also assumed that the generation of co-seismic landslide is independent of the state of landslide export from previous events. This is a reasonable assumption in the transport-limited regime, but the accumulation of co-seismic material on hillslopes in the connectivity-limited regime may at some point reduce the likelihood of landsliding. Similarly, the reactivation of co-seismic

A C C E P T E D M A N U S C R I P T
31 landslide is not accounted for in our approach [START_REF] Marc | Transient changes of landslide rates after earthquakes[END_REF] and may well contribute to further increase the release of debris in the channel network.

Evacuation of post-seismic landslide

A better understanding of the mountain range response to a sudden increase in sediment supply requires new data on 1. the initial landslide connectivity, 2. the processes controlling the debris transport dynamics on hillslopes and rivers and 3. The export of river sediment, including the partitioning between the production of fine and coarse sediment and its implication on the subsequent sediment evacuation. We discuss these three main points in the following.

Initial landslide connectivity

At the landscape scale, our results highlight the role of initial connectivity on the modulation of the sediment evacuation and export fluxes during the first few years that follow earthquakes (Fig. 7b). The assignment of an initial connectivity status to each landslide relies on only two empirical relationships between landslide area and connectivity [START_REF] Li | Connectivity of earthquake-triggered landslides with the fluvial network: Implications for landslide sediment transport after the 2008 Wenchuan earthquake[END_REF][START_REF] Roback | The size, distribution, and mobility of landslides caused by the 2015 M w 7.8 Gorkha earthquake, Nepal[END_REF]. The wide range of initial connectivity found in the literature, e.g. for earthquakes occurring in mountain ranges presenting different tectonic and climatic contexts, however, calls for a more systematic study of this process [START_REF] Dadson | Earthquake-triggered increase in sediment delivery from an active mountain belt[END_REF][START_REF] Li | Connectivity of earthquake-triggered landslides with the fluvial network: Implications for landslide sediment transport after the 2008 Wenchuan earthquake[END_REF][START_REF] West | Mobilization and transport of coarse woody debris to the oceans triggered by an extreme tropical storm[END_REF]. The differences in tectonic and climatic settings lead to a variety of topographic configurations that could promote (or inhibit) initial connectivity and the subsequent sediment transfers from hillslope to channel. For instance, hillslope size, valley width, slope distribution, runoff distribution (frequency-magnitude) could all potentially impact the initial connectivity by controlling: 1. the abundance of large landslides with long runout, thus promoting a higher initial connectivity (see Supplementary Figure S2 and S3) 2. the distance required for the landslide to reach the nearest river and 3. the possibility of the river to erode vertically (narrow valleys) or horizontally the landslide deposits (wide valleys), the former being more efficient [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF].

A C C E P T E D M A N U S C R I P T 32

Dynamic landslide connectivity

Over centennial timescales, our modelling framework demonstrates that the distribution of connection times determines the efficiency at which earthquake-mobilized sediment is exported from a mountain range over a seismic cycle (Fig. 11). However, because of a lack of theoretical and empirical constraint, these results are based on a simple description of dynamic connectivity by a constant velocity which value varies on a large range encompassing two extreme behaviors. The 'fullconnectivity' scenario implies the injection of the totality of landslides immediately during the coseismic period leading to extremely high sediment fluxes. On the contrary, the lowest value of the connection velocities implies that most of the evacuated landslides are the one initially connected to the fluvial network, the remaining of the landslide population being stored on hillslopes during several seismic cycles.

Potentially, hillslope to channel material transfer could be described by more complex physical laws, such as a diffusion law [START_REF] Roering | Hillslope evolution by nonlinear creep and landsliding: An experimental study[END_REF] or transport laws relevant to colluvial processes (e.g., Lague and Davy, 2003). These laws would make the connection times dependent on the geometric properties of its path to the river (i.e. local slope) in addition of the distance to the stream. It would affect the distribution of connection times and therefore the dynamic of post-seismic landslide evacuation (see Supplementary Fig S6 andS7). However, mean climate variables (precipitation, temperature, vegetation growth, storms and co-seismic reactivation, remain difficult to include in these geomorphic laws. These limitations call for a better description of this process based on a more extensive collection of empirical data (Fan et al., 2018).

Given the large spatial scales (i.e. > 100 km) and the difficulties of measuring total sediment fluxes (i.e. suspended plus bedload) in the field, remote sensing is well-suited to track the postseismic evolution of sediment masses [START_REF] Scaioni | Remote Sensing for Landslide Investigations: An Overview of Recent Achievements and Perspectives[END_REF]. This can benefit from a variety of techniques to investigate landscape evolution in 2D (satellite imagery) and 3D (e.g., photogrammetry, laser scanning) and, thus, aim toward a detailed quantification of the system through time [START_REF] Cook | Glacial lake outburst floods as drivers of fluvial erosion in the Himalaya[END_REF]Fan et al., 2018a;[START_REF] Jaboyedoff | Use of LIDAR in landslide investigations: a review[END_REF]. An important issue over centennial timescales is a better understanding of the contribution of rare events to long-term dynamics. As remote sensing only gives a first order quantification of sediment mass transfer, a full understanding of the processes

A C C E P T E D M A N U S C R I P T

controlling debris transport on steep hillslopes requires systematic monitoring of other variables such as precipitation, local gradient and slope self-healing processes, such as vegetation growth and grain coarsening (Fan et al., 2018a;Fan et al., 2018b;[START_REF] Zhang | Evolution of Mass Movements near Epicentre of Wenchuan Earthquake, the First Eight Years[END_REF]. The numerical modelling framework here could be adapted based on such observations, and help calibrate the transport laws which allow experiments to run over multiple seismic cycles.

River morphodynamics and sediment export

The Quakos post-seismic evacuation model relies on Eros simulation results and, therefore, shares its limitations. The Eros model contains the physical components required to simulate the feedbacks between a river flow, its geometry and transport capacity, in particular when large supplies of sediment are considered [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF][START_REF] Davy | A precipiton method to calculate river hydrodynamics, with applications to flood prediction, landscape evolution models, and braiding instabilities[END_REF]. As with any numerical model, Eros contains simplifications. For instance, the description of the landslide deposit and its erosion is lacking some processes. Once incised, the landslide deposit is not subjected to gravitational forces.

Therefore, it does not respect a specific angle of repose typical of non-cohesive granular material [START_REF] Cantelli | Experiments on upstream-migrating erosional narrowing and widening of an incisional channel caused by dam removal[END_REF][START_REF] Carrigy | EXPERIMENTS ON THE ANGLES OF REPOSE OF GRANULAR MATERIALS1[END_REF]. Additionally, water seepage through the landslide deposit could also structurally weaken the deposit [START_REF] Hovius | Landslide-driven drainage network evolution in a pre-steady-state mountain belt: Finisterre Mountains, Papua New Guinea[END_REF][START_REF] Meyer | Potential for Seepage Erosion of Landslide Dam[END_REF]. Accounting for these two processes might lead to more efficient lateral erosion of the deposit during its vertical incision, or lead to landslide dam mechanical failure. Ultimately, this would promote a higher sediment removal efficiency in Quakos. Contrary to Eros, Quakos does not provide 'real' sediment fluxes as the assumption implied in the treatment of sediment evacuation is similar to a detachment-limited model, i.e. once the sediment is entrained it is exported out of the catchment (Whipple & Tucker, 1999). This neglects the role of the deposition and re-entrainment of coarse sediments and the possibility that sediment can be stored in mountain valleys and in terrace fills (Pearce & Watson, 1986;Schwanghart et al., 2016). The remobilization fluxes can be used, however, as a proxy of the expected amplitude and duration of sediment fluxes at catchment outlets when subjected to large sedimentary forcing.

To simplify computations, we use a median grain size descriptors for the sediment transport law, whereas grain size distributions are more complex, i.e. from fine particles to boulders. The partitioning between fine and coarse particles and their spatial organisation remains poorly known. It

A C C E P T E D M A N U S C R I P T

34 is a critical element to constrain because 1) different timescales are expected for the sediment export by bedload and suspended load [START_REF] Hovius | Prolonged seismically induced erosion and the mass balance of a large earthquake[END_REF][START_REF] Yanites | How rivers react to large earthquakes: Evidence from central Taiwan[END_REF]Croissant et al., 2017), 2) the proportion of coarse material has an impact on hydro-sedimentary hazards as it propagates downstream (see next section) and 3) the relative proportion of sand at the subsurface has an impact on the transport rate of coarser sediment [START_REF] Parker | Surface-based bedload transport relation for gravel rivers[END_REF][START_REF] Wilcock | Surface-based transport model for mixed-size sediment[END_REF]. Accounting for a more complete grain size distribution (GSD) in transport models is an important next step to more fully describe the processes at play in landslide evacuation [START_REF] Cook | Glacial lake outburst floods as drivers of fluvial erosion in the Himalaya[END_REF][START_REF] Finnegan | Fluvial boulder transport controls valley blocking by earthflows in the California Coast Range , USA[END_REF].

It is also relevant to compare model predictions to sedimentary archives of past responses of earthquakes in mountain ranges (e.g. [START_REF] Howarth | Lake sediments record cycles of sediment flux driven by large earthquakes on the Alpine fault[END_REF]. Linked to this issue, determining the GSD of co-seismic debris is an essential step, for which new extensive data may now be extracted from 3D remote sensing techniques (e.g. [START_REF] Steer | 3D granulometry: grain-scale shape and size distribution from point cloud dataset of river environments[END_REF].

As for the mobilisation of landslide deposits, the role of storms of varying frequency and magnitude should be better addressed. If a more complete GSD is included in a model, it should be accompanied by the frequency-magnitude distribution of floods because that will modulate the entrainment of different grainsizes relative to one another. In addition, large floods may promote the removal of landslide deposits by lateral channel migration and influence the connectivity between the river and landslide deposits. In a modelling framework, such as Quakos, the threshold mean annual runoff of 2 m.yr -1 below which no transfer occurs is a limitation of applying a constant effective discharge [START_REF] Lague | The stream power river incision model: Evidence, theory and beyond[END_REF]. The use of an effective discharge and a median grain size comes with the inconvenience that for large D 50 values, sediment transport is inhibited leading to the unrealistic scenario of infinite sediment storage, A full-stochastic modelling approach, however, would predict sediment transport below runoffs of 2 m.yr -1 , albeit at a reduced rate, because of the action of infrequent very large floods. To go beyond the limitations of the constant effective discharge approach, a new generation of numerical models, driven by stochastic hydrologic forcing, need to be developed to predict channel morphodynamics at timescales up to several thousands of years.

Upscaling

A C C E P T E D M A N U S C R I P T 35

Given the return time of large earthquakes (classically several hundred years), current empirical observations can only resolve a small fraction of the duration of a seismic cycle. Other approaches, such as geochronological dating [START_REF] Wang | Perturbation of fluvial sediment fluxes following the 2008 Wenchuan earthquake[END_REF], detrital thermochronology [START_REF] West | Dilution of 10Be in detrital quartz by earthquake-induced landslides: Implications for determining denudation rates and potential to provide insights into landslide sediment dynamics[END_REF], source-to-sink analysis [START_REF] Howarth | Lake sediments record cycles of sediment flux driven by large earthquakes on the Alpine fault[END_REF] and sediment tracking using organic and inorganic clasts [START_REF] Frith | Carbon export from mountain forests enhanced by earthquaketriggered landslides over millennia[END_REF][START_REF] Nibourel | Provenance analysis using Raman spectroscopy of carbonaceous material: A case study in the Southern Alps of New Zealand[END_REF], can offer complementary and more integrative information. Landscape evolution modelling is, therefore, a useful tool to explore the dynamics of mountain ranges over long time scales, provided that they can benefit from the aforementioned data collection.

A successful application of such a modelling approach would give insights on several fundamental questions regarding the dynamics of mountain ranges over short time scales. First, how co-seismic debris impacts spatial and temporal bedrock incision patterns is still poorly understood. A large and sudden sediment supply can have a dual effect: locally protecting the bedrock from particle impacts, while providing tools for downstream incision [START_REF] Lamb | A model for fluvial bedrock incision by impacting suspended and bed load sediment[END_REF]). In that sense, the rates of sediment supply from hillslopes to channels seem to be important. A low velocity of connection would slowly release the material with time, potentially enhancing the tool effect. On the contrary, a high increase of sediment supply could lead to alluviation of the river channels and inhibit the vertical erosion of the bedrock [START_REF] Lague | Reduction of long-term bedrock incision efficiency by short-term alluvial cover intermittency[END_REF][START_REF] Turowski | Cover effect in bedrock abrasion: A new derivation and its implications for the modeling of bedrock channel morphology[END_REF]. This calls for a new approach to bedrock river incision dynamics at the scale of the seismic cycle, one which includes elements of bank erosion and the lateral mobility of river (e.g., Cook et al., 2016) and the emergence of channel width [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF]. Secondly, several studies have hypothesised the potential destructive role of earthquakes, i.e. the erosive action of a seismic event could be superior to the material influx brought by surface uplift [START_REF] Li | Seismic mountain building: Landslides associated with the 2008 Wenchuan earthquake in the context of a generalized model for earthquake volume balance[END_REF]Marc et al., 2016;[START_REF] Parker | Mass wasting triggered by the 2008 Wenchuan earthquake is greater than orogenic growth[END_REF]. The full mass budget of earthquakes (erosion versus uplift), however, can only be estimated if the post-seismic timescale of sediment evacuation is known. The scarcity of studies on this matter means that the potential destructive effect of earthquakes remains uncertain.

Earthquake-associated hazards

A C C E P T E D M A N U S C R I P T 36

The model Quakos developed in this study integrates three types of potential hazards. Two of them describe the primary hazards related to the co-seismic phase with the prediction of PGA and landsliding patterns and as a function of magnitude and fault geometry. The third involves the hydrosedimentary hazards linked to large increases in sediment supply and removal in the drainage network. As such, Quakos can be used as a tool to identify critical areas that could be subjected to a series of hazards in co-and post-seismic periods.

Whereas co-seismic areas has generally received the highest attention, the sediment mass movement in the following years after the earthquakes is also of importance. The downstream propagation of landslide-derived coarse sediment can trigger a series of potential hydro-sedimentary hazards in alluvial fans, such as, riverbed aggradation, river avulsion, bank erosion and increased frequency of inundation [START_REF] Hancox | The october 1999 mt adams rock avalanche and subsequent landslide dam-break flood and effects in poerua river, Westland, New Zealand[END_REF][START_REF] Korup | Sediment generation and delivery from large historic landslides in the Southern Alps, New Zealand[END_REF][START_REF] Robinson | Review Article: Potential geomorphic consequences of a future great (Mw Combining double low line 8.0+) Alpine Fault earthquake, South Island, New Zealand[END_REF]. Previous work has shown that alluvial fan response to large sediment supply is highly sensitive to sediment concentration and the total volume of bedload sediment feeding the fan (Croissant et al, 2017b).

Therefore, the river's ability to remobilize sediment immediately after an earthquake controls how the fans will respond to the supply. Our results show that the initial and dynamic landslide connectivity, as well as the evacuation of sediments by rivers, control the amplitude of sediment fluxes and delivery of catchments to the alluvial plains (Fig. 7d). Every scenario considered in this work presents high fluxes during the first few years after the seismic event. During this time, the fluxes are controlled by the proportions of landslides initially connected to the rivers. After this period, the connection velocity controls the temporal persistence and level of hazards. Fast connection velocities promote short-lived but high-amplitude sediment fluxes. On the contrary, slower connection velocities see a gradual release of sediment at lower rates but over a longer period. Immediate post-seismic risk mitigation plans should, therefore, consider the initial degree of connectivity of sediment sources.

Analyses of the mechanisms controlling the velocity of landslide deposit connection are also critical to determine the risk susceptibility of a specific area.

One feature that is not considered in this work, but that has a large hazard potential, is the formation and failure of dams [START_REF] Costa | Formation and Failure of Natural Dams[END_REF]; T. R. Davies et al., 2007;Korup, 2005a).

Dam break can cause tremendous damage and geomorphic changes in their downstream pathway,

A C C E P T E D M A N U S C R I P T 37 including on alluvial fans [START_REF] Cook | Glacial lake outburst floods as drivers of fluvial erosion in the Himalaya[END_REF][START_REF] Hancox | The october 1999 mt adams rock avalanche and subsequent landslide dam-break flood and effects in poerua river, Westland, New Zealand[END_REF].The inclusion of this process in a long-term landscape evolution modelling framework, however, is still a difficult task and we often have to rely on event-based approach [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF]Davies et al., 2007;[START_REF] Davies | Persistent alluvial fanhead trenching resulting from large, infrequent sediment inputs[END_REF].

Conclusion

In this paper, we study the post-seismic dynamics of rivers in the aftermath of a large injection of sediment by landslides. We first use a 2D modelling framework at the reach-scale to develop a model of the sediment export of a population of landslides over the scale of a fault rupture.

This spatially explicit reduced-complexity approach is justified by the high computational cost of numerical simulation of river morphodynamics if applied over large areas for decades or longer.

Our 2D morphodynamic modeling shows that the overall evacuation of landslide deposits are not significantly affected by the upstream sediment feed from other landslides. We then use the Eros modeling results to develop a generic way to predict landslide evacuation for a range of 𝑉 𝑙𝑠 /𝑄 𝑡 values.

This allows us to account for the morphodynamic adaptation of rivers to the injection of landslide debris, without having to describe the width and slope evolution explicitly.

At the mountain range scale, our newly developed model Quakos is applied to infer post-seismic sediment fluxes. For the co-seismic stage, Quakos produces a distribution of earthquake-triggered landslide volumes in the landscape, informed by empirical laws. Landslide debris is removed, but only when landslides are connected to river channels. We introduce a connection velocity by which debris on hillslopes migrates to river channels over time. Our results show that the post-seismic sediment transfers are strongly modulated by dynamic connectivity of landslide sediment to the drainage network and the mobilization of co-seismic debris on hillslopes. Two regimes are identified:

-a transport-limited regime in which the connectivity velocity only affects the transient sediment evacuation during a seismic cycle. In this regime, most of the sediment is evacuated within a seismic cycle and the duration of landslide evacuation only depends on the capacity of rivers to transport the co-seismic debris. For instance, our results suggest that the Southern Alps of New Zealand are likely to be in connectivity-limited conditions, for connection velocities less than 10 m.yr -1 .

The theoretical framework developed here is a first attempt to bridge the gap between the timescales of earthquake occurrence and sediment evacuation over seismic cycles. It opens up new research directions to better understand the residence time of co-seismic debris in the landscape, and the rejuvenation of hillslope surfaces with implications for post-seismic organic carbon mobilization and transfer [START_REF] Frith | Carbon export from mountain forests enhanced by earthquaketriggered landslides over millennia[END_REF], weathering fluxes from landslide deposits [START_REF] Emberson | Oxidation of sulfides and rapid weathering in recent landslides[END_REF], hydro-sedimentary hazards [START_REF] Croissant | A precipiton-based approach to model hydro-sedimentary hazards induced by large sediment supplies in alluvial fans[END_REF] and the interpretation of sedimentary archives [START_REF] Howarth | Lake sediments record cycles of sediment flux driven by large earthquakes on the Alpine fault[END_REF].

Overall, this paper points towards the lack of information and physical description of the transfers of sediment from hillslopes to river channels and calls for intensifying empirical study, including long-term monitoring, and the development of spatially explicit numerical models accounting for stochastic forcing and multiple grain size.

Figure 1 |

 1 Figure 1 | Illustration of the geomorphic impact of landslides at different spatial scales a. Aerial image of the Hapuku river landslide (taken the 5 th Dec. 2016) triggered by the 2016 Kaikoura

Figure 2 |

 2 Figure 2 | Eros model set up. a. The case of a single landslide (in brown) evacuated in a bedrock channel (in beige). This snapshot illustrates an advanced stage of one simulation where the river incise the deposit with a reduced width. b. The case of landslide cascade deposited in the same channel.

  by the river at the rate set by the initial transport capacity of the bedrock river .

Figure 3 |

 3 Figure 3 | Morphodynamic evolution of a reach evacuating several landslides. a. Snapshots of different stages of the Eros simulation. This is the case of 2 landslides of V s = 1.10 6 m 3 separated by a

Figure 4 |

 4 Figure 4 | Predicting the dynamic evacuation of landslides using Eros results. a. Series of fitting functions applied to predict the landslide export when x% of it has been evacuated. b. and c. Comparison between our predictions (yellow dots, green line) and Eros output (grey line) for 2 values of 𝑉 𝑙𝑠 /𝑄 𝑇 .

Figure 5 |

 5 Figure 5 | Quakos workflow used to predict landsliding pattern at the mountain range scale. a. Quakos workflow. b. PGA pattern predicted by the model for a M w 8 scenario on the Alpine Fault, New Zealand. c. Map of the landslide density. d. Landslide pattern predicted by Quakos. The view is focused on catchments located on the Central Southern Alps. Yellow dots indicate landslides, the sizes are a function of landslide area. White dots outside of the considered catchments only represent the location of landslides. Notation: V ls , landslide volume.

Figure 6 |

 6 Figure 6 | Relationship between the total volume of landslides triggered by earthquakes and the magnitude, in the case of a strike-slip fault. Quakos results (circles) are coloured as a function of the depth of the earthquake compared to the surface (R 0 ). Quakos outputs are compared to two empirical models from Marc et al, [2016] and Keefer, [1999].

Figure 7 |

 7 Figure 7 | Temporal evolution of the M w 7.9 earthquake simulations with different connectivity properties. a. Temporal evolution of the landslide volume for different connection velocities. b. Temporal evolution of the number of landslides. c. Sediment remobilization fluxes. d Temporal evolution of the number of active landslides, i.e. landslides that are being connected to the drainage network and being actively evacuated by the river. Note: The grey area indicates the estimated return time of a M w 8 earthquake on the Alpine Fault.

  -constant delivery of sediment with a steady number of landslides being active during the first 200 years.

Figure 8 |

 8 Figure 8 | Probability density function of landslide volumes at different time steps. a. For the Full connectivity case. b. For the case 𝑢 𝑐𝑜𝑛 = 10 m.yr -1 . c. For the case 𝑢 𝑐𝑜𝑛 = 1 m.yr -1 .

Figure 9 |

 9 Figure 9 | Temporal upscaling over several seismic cycles. a. Times series of earthquakes generated on the faults. It is characterised by a series of mainshocks of M w = 7.9 followed by the sequences of aftershocks. The dot size and color are a function of the total volume of the landslide population (V ls,tot ). The grey dots indicate earthquakes that have not triggered any landslides. b. Evolution of the remobilization fluxes. c. Evolution of the total volume of sediment mobilized by the successive earthquakes. Grey lines represent the mean sediment thickness that would be deposited on the total area affected by landsliding. d. Evolution of the number of landslides in the mountain range.

  is equivalent to 120% (𝑢 𝑐𝑜𝑛 = 1 m.yr -1 ) to 600% (𝑢 𝑐𝑜𝑛 = 0.1 m.yr -1 ) of the average value of earthquake-triggered initial landslide volume. These cases are also characterized by the persistence of several thousands of landslides in the mountain range.

Figure 10 |

 10 Figure 10 | Landslide generations. a. For the Full connectivity case. b. For the case u con = 1 m.yr -1 . Notation: Ms: mainshocks.

Figure 11

 11 Figure 11 | a. Time necessary to evacuate a landslide cluster as a function of mean annual runoff and connection velocity. Black lines represent a 𝑡 𝑒𝑥𝑝 value of 25, 50 and 75 % of the duration of a seismic cycle. b. Relative importance of the two timescales controlling post-seismic landslide evacuation.

  -limited regime in which landsides debris moves slowly on hillslopes compared to the return time of a large earthquake. In this regime, a small fraction of landslides are exported during a seismic cycle, and co-seismic debris can progressively build up in the landscape. A dynamic equilibrium is ultimately reach after a transient period whose duration depends on the connection velocity of the landslides. At this stage, the sediment flux leaving the mountain range during a seismic cycle equals the co-seismic sediment production.

Table 1 :

 1 Value of the parameters of equation 5.

	T exp	𝜹	𝝁	𝝋	𝜷	R 2
	20%	2.1	0.03	1	0.1	0.90
	30%	4.6	0.035	1	0.1	0.88
	40%	7	0.033	1	0.1	0.90
	50%	10.8	0.033	0.95	0.1	0.97
	60%	13	0.035	0.9	0.1	0.93
	70%	18.1	0.030	0.97	0.1	0.94
	80%	26	0.031	0.85	0.1	0.92

Table 2 |

 2 Parameters used in QUAKOS

	Parameter	Notation	Values	Units
	Tectonic
	Fault length	F L	400	km
	Fault width	F W	19	km
	Fault dip	F D	60	°
	Fault rake	F R	172	°
	Earthquake rupture -Strike-slip
	Shear modulus	μ	33	GPa
	Exponent	β	0.66	/
	Constant	C 1	15	/
	Constant	C 10	F W	/
	Constant	C 11	1	/
	Constant	C 2	3.6.10 -5	/
	Landsliding
	Min. landslide area	A ls,min	50	m 2
	Max. landslide area	A ls,max	2.10 6	m 2
	pdf(A ls ) roll-over position	a	2000	m 2
	pdf(A ls ) roll-back	s	-200	m 2
	pdf(A ls ) tail exponent	ρ	1.4	/
	Volume-Area prefactor	α	0.05	/
	Volume-Area exponent	γ	1.5	/
	PGA-P ls coefficient	α p	4	/
	PGA-P ls coefficient	β p	0.5	/
	Connectivity parameter	m	0.87	/
	Connectivity exponent	ω	0.34	/

Acknowledgments

T.C., P.S. and L.J. acknowledge supports by the EROQUAKE project funded by the Agence Nationale de la Recherche (ANR-14-CE33-0005) and by the Region Bretagne. D.L. acknowledge support by CNRS/INSU/ALEAS project SEDIQUAKE. R.G.H. was supported by a NERC Standard Grant (NE/P013538/1). We thank the two reviewers whose comments and corrections helped to improve the quality of our manuscript. Discussions with Jamie Howarth and Alexander Densmore were greatly appreciated.

Code availability

The code, Eros, used to generate the morphodynamic simulations can be accessed at https://osur.univ-rennes1.fr/eros/. The source code can be obtained by making an inquiry to P.D. (philippe.davy@univ-rennes1.fr). The code, Quakos, can be obtained by making an inquiry to T.C. or P.S.

Authors contribution

All authors contributed equally to the design of the study and the writing of the paper. T.C. designed and performed the simulations and processed their results. T.C., P.S. developed the Quakos model. P.D. developed the Eros morphodynamic model.

River transport capacity