Modelling post‐earthquake cascading hazards: Changing patterns of landslide runout following the 2015 Gorkha earthquake, Nepal

Coseismic landslides represent the first stage of a broader cascading sequence of geohazards associated with high‐magnitude continental earthquakes, with the subsequent remobilization of coseismic landslide debris posing a long‐term post‐seismic legacy in mountain regions. Here, we quantify the controls on the hazard posed by landslide remobilization and debris runout, and compare the overlap between areas at risk of runout and the pattern of post‐seismic landslides and debris flows that actually occurred. Focusing on the 2015 Mw 7.8 Gorkha earthquake in Nepal, we show that the extent of the area that could be affected by debris runout remained elevated above coseismic levels 4.5 years after the event. While 150 km2 (0.6% of the study area) was directly impacted by landslides in the earthquake, an additional 614 km2 (2.5%) was left at risk from debris runout, increasing to 777 km2 (3.2%) after the 2019 monsoon. We evaluate how this area evolved by comparing modelled predictions of runout from coseismic landslides to multi‐temporal post‐seismic landslide inventories, and find that 14% (85 km2) of the total modelled potential runout area experienced landslide activity within 4.5 years after the earthquake. This value increases to 32% when modelled runout probability is thresholded, equivalent to 10 km2 of realized runout from a remaining modelled area of 32 km2. Although the proportion of the modelled runout area from coseismic landslides that remains a hazard has decreased through time, the overall runout susceptibility for the study area remains high. This indicates that runout potential is changing both spatially and temporally as a result of changes to the landslide distribution after the earthquake. These findings are particularly important for understanding evolving patterns of cascading hazards following large earthquakes, which is crucial for guiding decision‐making associated with post‐seismic recovery and reconstruction.

remain active for several years as debris is remobilized by subsequent rainfall, resulting in the expansion and extension of the original landslide extent (Fan, Domènech, et al., 2018;Kincey et al., 2021). The hazards posed by this remobilization process represent a significant, but generally unquantified, proportion of total landslide hazard experienced in post-earthquake landscapes.
Because the spatial distribution, extent, and stability of individual active landslides change significantly through time, through both the remobilization of existing landslides and the formation of new postseismic landslides (Fan, Domènech, et al., 2018;Kincey et al., 2021;Marc et al., 2019), the future hazard posed by the runout of debris must also change accordingly. Critically, thresholds for landslide triggering may also reduce, rendering previously safe areas hazardous (e.g., Dadson et al., 2004), and the most prevalent mass wasting mechanisms may shift to become more dominated by debris flows (Zhang & Zhang, 2017). Importantly, much of the coseismic landslide debris may also lie upstream of the channel network (Li et al., 2016) and so can run out over areas that have not previously experienced landslides or debris flows. Thus, the evolving landslide hazard after a large earthquake consists of at least four distinct components: (1) remobilization of coseismic debris, (2) continued failure of coseismic landslides, (3) new post-seismic landslides on previouslyunfailed hillslopes, and (4) remobilization of post-seismic debris.
Understanding the changing significance of each of these different components of the hazard chain should therefore be a fundamental aspect of any comprehensive assessment of post-seismic risk.
Previous studies assessing the post-seismic evolution of landslide hazard have focused on the analysis of multi-temporal landslide inventories to assess changes to existing failure extents, the occurrence of new post-seismic landslides, or both (e.g., Fan et al., 2019;Fan, Domènech, et al., 2018;Fan, Zhang, et al., 2018). For example, after the 2015 M w 7.8 Gorkha earthquake in Nepal, such an analysis has demonstrated that 3.5 years after the earthquake, landslides remained more numerous and covered a larger total area than they did on the day of the earthquake . Importantly, the changing post-seismic landslide distribution comprised both persistently-active coseismic landslides but also new post-seismic landslides that had developed since the earthquake; thus, the sources of landslide debris with the potential to runout must change over time. Similarly, in the years following the 2008 M w 7.9 Wenchuan earthquake in China, the primary hillslope failure mechanism shifted from landslides to debris flows, with Huang and Li (2014) documenting a landslide/debris flow ratio of 5:1 for the pre-seismic period and 1:1 for the initial 5-year post-seismic period. This increase in debris flow prevalence reflected the abundance of loose coseismic debris and the associated reduction in the hydrological triggering threshold required for mobilization (Fan, Zhang, et al., 2018;Ma et al., 2017). The post-seismic increase in debris flow occurrence reduced through time, potentially as a result of progressive exhaustion of supply (Qu, 2019;Yunus et al., 2020), grain coarsening due to the loss of fine sediment (Domènech et al., 2019), or revegetation of failure scars (Shen et al., 2020;Yang et al., 2018).
Debris-flow runout distances have also generally decreased through time, perhaps due to a reduction in sediment mobility and downslope progression of debris flow initiation positions (Fan, Zhang, et al., 2018;Zhang & Zhang, 2017). However, the time period over which increased debris mobilization persists after an earthquake appears to vary considerably between both settings and studies. In their study of post-seismic debris flows in the 2 years following the Gorkha earthquake, Dahlquist and West (2019) suggested that there was only a short-lived transient increase in debris flow rates, with the available coseismic sediment supply being largely exhausted during the first monsoon and a reduction in the number of new debris flows back to pre-seismic levels within a year. This rapid return to pre-seismic conditions was argued in part to be due to the low proportion ($2%) of coseismic landslides that actually transitioned into post-seismic debris flows (Dahlquist & West, 2019;Roback et al., 2018). In contrast, more recent and longer-term studies of post-seismic hillslope evolution following the Gorkha earthquake have indicated that landslide and debris flow hazards are still high relative to pre-seismic conditions.
For example, in their field-based study of landslide development during the period 2015-2018, Tian et al. (2020) documented repeated activity and continued hazard across the majority of the investigated sites, with a notable shift in the dominant failure mechanism towards debris flows. Similarly, the 2020 monsoon is known to have triggered extensive debris flows across large areas of central Nepal that were badly affected by the Gorkha earthquake , suggesting a persistent runout hazard legacy associated with the earthquake. This persistence has been noted after the Wenchuan earthquake as well, and has been ascribed to an abundance of coseismic sediment even years after the earthquake (Huang & Li, 2014;Zhang & Zhang, 2017), changing source area form and location Zhang et al., 2016) and spatio-temporal variability in triggering factors such as high rainfall events (Ma et al., 2017;Yunus et al., 2020).
A limitation of most current approaches to assessing changing landslide and debris-flow hazard after a large earthquake is that they have focused on mapping where landslides have already occurred, resulting in limited capacity to predict how the hazard will evolve. An alternative approach is to forecast the potential evolution of the hazard footprint using a runout model, and then to compare modelled outputs with multi-temporal inventories to assess and refine the runout model based on where the modelled hazard has been 'realized'.
Because we cannot reliably identify which landslides are likely to be remobilized at the full event scale, a precautionary approach is to assume that further runout of any landslide source remains possible.
Typical approaches are to model potential runout pathways based on manually-mapped landslides as potential source areas (Aaron et al., 2019), or to use a predicted source area distribution based on predefined variables (Kappes et al., 2011;Pastorello et al., 2017) or a threshold-based landslide susceptibility model (Melo & Zêzere, 2017;Paudel et al., 2020). Such assessments usually provide a snapshot of the hazard, and it is hard to test or update the hazard due to a lack of multi-epoch landslide inventories that describe how the landslide footprint evolves. Whilst one-off regional-scale debris flow hazard assessments are commonplace (e.g., Blais-Stevens & Behnia, 2016), the degree to which runout occurs during the initial mass movement versus during reactivation remains poorly understood, and so is difficult to account for in hazard assessments.
To address these issues, we use multi-epoch landslide inventories from the 2015 Gorkha earthquake  to assess the changing extent of potential runout from coseismic and post-seismic landslides in the 4.5 years following the earthquake. We use Flow-R, a spatially distributed empirical model for regional-scale estimation of debris runout (Horton et al., 2013), to model the spatial footprint of potential runout across the area that was most affected by the earthquake, using mapped landslides as source areas. We repeat this simulation based upon 13 separate epochs of landslide mapping at approximately 6-month intervals (2014-2019) and compare modelled runout to subsequent mapped landslides. This approach allows us to explore how changes in the spatial distribution and character of landslides influence the evolution of post-earthquake landslide runout and the degree to which landslides achieve or 'realize' the modelled runout, and to identify the controls on that realization of runout potential. We use this information to derive a synoptic overview of the evolution of post-earthquake landslide runout and to describe the regional-scale characteristics of landslide runout pathways, which together provide inputs into a more holistic understanding of the long-term hazard chain associated with earthquake-triggered landslides. and an additional coseismic inventory for 2015. Landslide footprints incorporated both source areas and deposits due to the resolution of the satellite imagery. Assigning levels of post-failure reactivation and remobilization using remotely-sensed imagery is not straightforward (e.g., Fan, Domènech, et al., 2018) and so we mapped all landslides visible in each epoch independently, irrespective of whether they were already present within a preceding inventory. This approach is necessary since mapping only new or substantially-altered landslides in each epoch would have removed persistent landslides which could potentially act as source zones for later runout from the analyses. Our approach therefore makes the assumption that bare, unvegetated ground equates to the presence of exposed rock or sediment that could be mobilized in a future event, which is justified based on field observations (Tian et al., 2020) and a precautionary approach to modelling cascading hazards across such a large spatial area .
A detailed description of our mapping approach and its implications for time series analysis of modelled runout is provided in Supporting Information Methods S1, while full details of the multitemporal landslide inventory construction and analysis are provided by Kincey et al. (2021).

| Runout modelling
To assess the potential runout from existing landslides, we used the flow path assessment of gravitational hazards at a regional scale (Flow-R) model version 2.0 (Horton et al., 2013). Flow-R is a spatially distributed empirical model designed to model runout paths across large spatial extents with minimum input data requirements. Runout paths from defined source areashere manually mapped landslidesare propagated on the basis of a spreading algorithm that controls the route and extent of the flow, and friction laws that determine the runout distance. The flow volume and mass are not directly considered, as these cannot be accurately quantified across large regions. Instead, the model is well-suited to regional susceptibility assessments where the full range of possible runout pathways from a large set of distributed sources must be considered. Flow-R has been utilized in a range of different applications, including debris flow hazard assessments in Switzerland (Horton et al., 2008), France (Kappes et al., 2012), Italy (Blahut et al., 2010) and Norway (Fischer et al., 2012), as well as F I G U R E 1 (a) Location of study area showing the coseismic landslide distribution (Epoch 4) in red. (b) Extract from the runout model (Rout_Mod) results showing maximum runout susceptibility (p(Rout)) from coseismic sources (Ls_Map; shown in red) (AW3D 5 m DEM ©JAXA, RESTEC and NTTDATA). [Color figure can be viewed at wileyonlinelibrary.com] modelling other gravitational hazards such as rockfall (Losasso et al., 2016;Michoud et al., 2012), snow avalanches (Horton et al., 2009;Jaboyedoff et al., 2012) and rock avalanches (Oppikofer et al., 2016).
Runout modelling used the 13 manually-mapped landslide inventories to sequentially define the potential source areas (Section 2.1) and a 10 m digital elevation model (DEM), resampled from a 5 m resolution Advanced Land Observing Satellite World 3D (AW3D) dataset, as the most appropriate compromise between accuracy of modelled flow paths, reduction of topographic noise, and processing time (Claessens et al., 2005;Fischer et al., 2012;Horton et al., 2013) (see Methods S2 for full details of the Flow-R model inputs and parameters). In this configuration, the model is testing the potential for runout from the full landslide footprint, including continued runout from the landslide source and the potential for remobilization of the deposit.
For our analyses, spreading of the runout path was based on a modified version of Holmgren's multiple flow direction algorithm (Holmgren, 1994) that includes an additional height factor (dh = 1 m) designed to increase the elevation of the central processing cell and so account for any height errors in the DEM (see Horton et al., 2013). Tuning of model input parameters to fit individual landslides was not appropriate given the regional scale of the study area (> 24,000 km 2 ) and the number and variability of source areas used in the modelling (> 195,000 landslides). Our parameterization of the model was therefore based on input values that have been developed for modelling runout in similar mountainous terrain by previous studies, and through initial testing of the model to generate plausible runout pathways (see Methods S4). Sensitivity analysis of the different parameters for a subset area located in Sindhupalchok district, central Nepal (Supporting Information Figure S3), showed that the Holmgren exponent value (x) is the most sensitive parameter in terms of influencing the overall modelled runout extent and susceptibility distributions. The travel angle and velocity threshold values determine the degree to which flow is permitted to continue downslope, but the results were relatively insensitive to the choice of the height factor (dh). The sensitivity analysis shows that our parameter choices represent a plausible but precautionary approach to the modelling of runout, which is appropriate given the potential implications of the work in terms of risk to life and livelihood (Methods S4). Sensitivity testing showed that variations in modelled runout extent vary by at most 43% depending on specific parameter choices, but most parameter sets (70%) are within 10% of the runout area modelled by this study (Supporting Information Table S3).
For every 10 m Â 10 m cell, the model assigned a susceptibility value p(Rout) with a value from 0 (no chance of being within a runout path) to !1 (very likely to be within a runout path). Landslide areas had a value of 1. The model accounted for cells located within multiple modelled runout pathways, from which a set of summary statistics were generated, including: maximum runout susceptibility from all upslope sources combined ( Figure 1b); the sum of all runout susceptibilities, which assesses total runout susceptibility irrespective of source; and runout count, which accumulates the number of distinct sources that could impact a cell irrespective of susceptibility. We refer to the locus of cells where p(Rout) > 0 as Rout_Mod. The sequential model runs based on the 13 multi-temporal landslide inventories allowed an assessment of how Rout_Mod changed over the 4.5 years since the 2015 earthquake.

| Time series analysis
Quantitative assessment of model performance using standard error statistics is somewhat complicated because Flow-R is a forward model, and runout from existing landslide source locations is timeand trigger-dependent. Thus, the occurrence or non-occurrence of actual runout for any given cell at a single point in time does not simply indicate model success or failure. In other words, the lack of realized runout does not show that the model failed in its prediction, because the actual occurrence of runout may still happen in the future if there is a sufficient trigger event. This study is therefore focused on an assessment of how the runout distribution and realization change through time after an earthquake, using a precautionary but plausible modelling approach.
We first assessed the net change in Rout_Mod across the study area using aggregate statistics, per-cell values, and frequency distributions of p(Rout) for each epoch in turn. All cells mapped as landslides are referred to as Ls_Map, with any cell classified as Ls_Map relating to a model source location for that epoch and so excluded from the associated list of Rout_Mod runout cells. We refer to cells in Rout_Mod that experienced landsliding in later epochs, which we term 'realization' of the modelled runout, as Rout_Real. For each epoch, each 10 m Â 10 m cell in Rout_Mod was categorized based on whether it became part of a landslide during that epoch (i.e., the cell was contained within Ls_Map and so it became part of Rout_Real), remained within the modelled runout extent (the cell remained in Rout_Mod but outside of Ls_Map), or no longer fell within the modelled runout extent (i.e., p (Rout) = 0), which typically occurred when the associated upslope landslide became revegetated and was no longer visible as bare ground. This process was repeated for each epoch, allowing the stepwise evolution of Rout_Mod and the match between this and Ls_Map to be assessed (see Methods S3 for additional description of the runout variables described above).
Since our analysis considered realized runout to include any cells within a modelled runout area that experienced landsliding in later epochs, this realization could plausibly occur as a result of both changes to pre-existing landslides and the occurrence of entirely new landslides. This approach is justified from the perspective of informing a comprehensive post-seismic hazard analysis, but it is still important to also consider what proportion of realization comes from existing landslides versus entirely new landslides. To calculate this, we used a series of spatial selection queries to identify which landslides had a direct physical intersection with earlier landslides and considered these to be pre-existing landslides. Any landslides that had no intersection with landslides from earlier epochs were classified as entirely new landslides (see Methods S10 for full methodological details).
Note that the percentage of Rout_Mod that is realized by the end of each epochthat is, the ratio of Rout_Real to Rout_Modis equivalent to the model precision. We would expect this quantity to change over time, as runout proceeds and as areas within Rout_Mod that are downslope of existing landslides become inundated with debris. Similarly, we would expect the model recallthe ratio of Rout_Real to the sum of Rout_Real and any new landslides that occur outside of the modelled runout areato also change over time, as new landslides occur outside of Rout_Mod. We thus assessed model performance using sequential precision-recall curves, which are preferable to other performance measures when Rout_Mod is a small fraction of the overall area (e.g., Saito & Rehmsmeier, 2015). We used Flow-R to calculate Rout_Mod using the coseismic landslide footprints (Epoch 4) as the runout source areas, and evaluated model precision and recall up to that point in time after each subsequent epoch. For this portion of the analysis, we treated realization as a cumulative process; that is, we considered that a cell in Rout_Mod was realized if it appeared as a landslide in a subsequent epoch, even if the cell was not mapped as a landslide in later epochs. This avoids complications caused by vegetation growth or obstruction of the ground surface in some but not all epochs (see Kincey et al., 2021).
To assess broader spatial variability in Rout_Mod, Ls_Map, and Rout_Real, we generated a set of geomorphological slope units across the study area, based on the methodology developed by Alvioli et al. (2016). This approach uses a DEM to partition a landscape into individual terrain units that are defined by hydrological and geomorphological boundaries , and has been shown to be appropriate for susceptibility assessments across large spatial areas (Domènech et al., 2020;Jacobs et al., 2020;Tanyas et al., 2019). The size distribution of output slope units is primarily determined by parameters controlling the flow accumulation thresholds and the circular variance in terrain aspect that is permitted within a single slope unit, which together define the acceptable degree of aspect homogeneity between adjacent units. Slope units were generated within Grass GIS version 7.8.4 using the minimum parameter settings from the range of values recommended by Alvioli et al. (2016) which delimit slope units to a scale that matches observed approximate hillslope length scales across the study area (see Methods S5 for full details of slope unit parameters). Each resulting slope unit (n = 13,456) was attributed with a range of landslide and runout statistics for each mapping epoch, including: the number of landslides, the total area of landsliding, the total modelled runout area, and statistics summarizing the maximum and summed runout susceptibilities (minimum, maximum, range, mean, sum, standard deviation). All of the topographic and environmental variables were also aggregated to each slope unit as tabular attributes using the same set of summary statistics as were used for the runout susceptibility values.

| Identifying locations favourable for landslide runout
Analysis of potential controls on runout evolution and hazard realization focused on eight topographic variables that have previously been shown to influence the occurrence of coseismic and post-seismic hillslope landslides: elevation, slope, aspect, normalized distance to stream channel, profile and plan curvature, upslope contributing area, and Melton ratio (defined here as per-slope unit relief divided by the square root of slope unit area) (e.g., Kincey et al., 2021;Parker et al., 2015;Robinson et al., 2017). Three event-specific variables from the Gorkha earthquake were also included: slope aspect relative to the epicentre, the Euclidean distance to the epicentral location, and the distance to the nearest mapped coseismic landslide. All control variables were analysed at 10 m resolution using the same AW3D DEM (see Methods S6 for full details of the potential control variables).
The significance of each variable was assessed in two ways: by (1) comparing per-slope unit aggregated summary statistics of each variable (see Section 2.3) to corresponding modelled runout areas and realization percentages, and (2)

| Characteristics of post-earthquake runout evolution
Evolution of landslide area (Ls_Map) in the pre-seismic, coseismic, and post-seismic phases was described by Kincey et al. (2021), so we provide only a brief summary here. The total landslide area for the three pre-earthquake epochs was consistently low, ranging between 60 and 64 km 2 , equivalent to 0.2-0.3% of the overall study area (

| Post-seismic changes in modelled coseismic runout realization
Although runout susceptibility evolved as a result of changes to the number, distribution, and form of landslides in subsequent postseismic epochs, it is still valuable to assess the degree to which the potential runout area modelled from just a coseismic landslide inventory was fulfilled in the years following the earthquake. This provides insights into timescales over which a runout assessment carried out in the immediate aftermath of an earthquake might remain relevant, as well as the degree to which the realization of the modelled runout area varies under different model susceptibility thresholds.
We focus initially on the modelled runout area Rout_Mod from the coseismic landslides (Epoch 4). When all modelled runout susceptibility values are included, 14% (85 km 2 ) of Rout_Mod cells from coseismic landslides became Rout_Real at some point in the 4.5 years after the 2015 earthquake ( Figure 5a). In other words, 86% of the modelled hazard remained 'unrealized' by the end of the study. This is a highly conservative test of the runout model, as cells with very low susceptibilities are included. We therefore apply a moving threshold to progressively examine the evolution of cells with higher values of p (Rout). As this threshold increases, whilst the area of Rout_Mod naturally decreases, Rout_Real initially increases (black line on Figure 5a), indicating increased model precision. This is greatest for a p(Rout) threshold of ≥ 0.6, at which 32% (10 km 2 ) of all modelled coseismic runout cells (32 km 2 ) became a landslide at some point during our study period. The percentage of Rout_Mod cells that were realized decreases marginally for higher thresholds, but this also leads to a pronounced decrease in the predicted runout area ( Figure 5a); for  Figure 5b), but at a decreasing rate as more of the coseismic Rout_Mod was progressively realized. Precision is maximized at p(Rout) values of 0.5-0.6. For all epochs, the model is considerably more skilled than a random classifier (Figure 5b). Note, however, that maximum model recall values decrease somewhat in later epochs, likely due to the occurrence of new landslides that do not fall within the coseismic Rout_Mod.
Analysis of the intersections between landslides from different epochs shows that 90% of the total post-seismic realization of the coseismic runout area came from landslides that were physically connected to the original coseismic landslide footprints, with only 10% originating from entirely new landslides (Methods S10). When evolving downslope connections between landslides are considered, this value increases to 92% of the realized area being related to preexisting landslides and only 8% coming from entirely new landslides.
In both cases, the proportion of realized runout that is attributable to pre-existing coseismic landslides decreases through time. For direct intersections with coseismic landslides, the values decrease from 91% in post-monsoon 2015 (Epoch 5) to 85% in post-monsoon 2019 (Epoch 13), with the equivalent figures for the evolving downslope connections being 94% and 89% ( Figure S14). For all post-seismic epochs combined, we find that 62% of the area impacted by new post-seismic landslides that are directly connected to coseismic landslide polygons is within the modelled runout extent, compared with 38% that is outside of the modelled extent Methods S10). When analysed through time, we find that the proportion of newly impacted ground that is within the modelled runout extent increases from 61% in post-monsoon 2015 (Epoch 5) to 65% in 2016 (Epochs 6 and 7), before decreasing to 60% by post-monsoon 2019 (Epoch 13) ( Figure S15).

| Runout hazard realization in a changing postseismic landscape
We next examine the realization of Rout_Mod from all coseismic and post-seismic epochs. Disaggregating the realized hazard data by epoch allows us to establish the time between when a cell is first predicted as being within a runout area (Rout_Mod) and when it is first intersected by a mapped landslide (Rout_Real), here termed the 'reali- Rout_Mod remains high for all post-seismic epochs (Figure 2). This implies that the modelled runout area is shifting spatially, reflecting changes to landslide footprints and the addition of new landslides in later epochs, meaning that the analysis of hazard realization needs to also track changes occurring in individual post-seismic epochs.
F I G U R E 7 (a) Extent of modelled runout in each of the mapping epochs that was realized by landslide occurrence within any subsequent epoch (Rout_Real), differentiated by susceptibility (p(Rout)) threshold. Note: The realization area data on the y-axis are stacked. (b) Model precision by epoch, defined as the realized runout area (Rout_Real) as a fraction of the total coseismic runout area (Rout_Mod), again differentiated by susceptibility (p(Rout)) threshold. For all epochs, the model achieves maximum precision at a p (Rout) threshold of 0.55-0.65; these points are connected for reference to show their evolution in time for both pre-and post-monsoon inventories. Coloured dots on the y-axis indicate the proportion of all Rout_Mod that was realized by the end of the study for each mapping epoch.

| Controls on changing post-earthquake runout characteristics
No clear correlations are apparent at slope unit level between the topographic variables (Section 2.4) and the percentage values of realized coseismic hazard (Rout_Real), suggesting that any potential causal relationships are likely to be manifest below slope unit scale. In contrast, at 10 m resolution, the kernel density estimate differencing of these variables shows marked variation between Rout_Mod cells for which the hazard was realized within 4.5 years after the earthquake, and cells for which it was not ( Figure 10). Modelled runout cells where the hazard became realized preferentially occurred at elevations above approximately 2750 m, and in particular between 3500 and 5000 m (Figure 10a), whereas at lower elevations (< 2750 m) Rout_Mod was markedly less likely to be realized. Similarly, realized runout preferentially occurred at slope angles > 35 , with a modal peak at 45 , but much less frequently for lower slope angles Rout_Mod cells that experienced landslides over the study period also preferentially occurred on slopes with aspects between approximately 90 and 240 , representing hillslopes facing orientations between east, south and west-southwest, with a modal peak around south-southeast (Figure 10g). Rout_Real cells on more northerly-facing slopes were noticeably less common. When aspect is adjusted to reflect hillslope direction relative to the Gorkha earthquake epicentre, the distribution shows that realized runout was preferential on slopes facing obliquely away (90 -175 ), despite these hillslopes being relatively infrequent in the study area (Figure 10h), and in contrast, Rout_Mod cells on slopes facing the epicentre generally did not become Rout_Real. Realized runout also preferentially occurred close to coseismic LS_map cells, especially for distances < 50 m (Figure 10i).
Beyond this distance, Rout_Mod cells were much less likely to experience landslides in subsequent epochs, suggesting that realization of the modelled runout extent was mainly achieved via iterative extension of coseismic landslides, rather than via the occurrence of entirely new landslides in downslope locations.

| DISCUSSION
Post-earthquake landsliding can represent a significant secondary hazard in the months and years after mountain region earthquakes. A considerable portion of this hazard arises from the remobilization of previously failed material, and so understanding the evolution of the associated runout holds the potential to locate, and therefore mitigate, this component of post-earthquake hazard. Our results provide the first regional-scale assessment of how this runout hazard changes in the years following a large continental earthquake, including the extent to which the overall potential runout area is realized and what local factors control the degree to which this occurs.

| Evolution of landslide runout following a large earthquake
Our results demonstrate that the overall area modelled to be at risk of runout from mapped landslides increased considerably after the 2015 F I G U R E 1 0 Kernel density distributions of all modelled coseismic runout cells (Rout_Mod) for selected topographic, seismic and distancebased variables. Differences in kernel density values relative to the overall distribution are shown for those cells where the coseismic runout hazard was later realized (Rout_Real; red) and for those where it was not (Rout_Mod ≠ Ls_Map; blue). Shaded areas show distributions of each variable for all coseismic runout cells (Rout_Mod) within the study area for reference. [Color figure can be viewed at wileyonlinelibrary.com] Gorkha earthquake. The area at risk of potential runout peaked after the 2015 monsoon, but importantly remained above coseismic levels through the end of the 2019 monsoon (Figure 2). This pattern of runout susceptibility fits with our understanding of the changing pattern of landslides following this earthquake, which suggests that the overall landslide footprint remained large through at least 2019 relative to coseismic levels Rosser et al., 2021). Our results are in contrast with those of Dahlquist and West (2019), who documented a rapid decline in debris-flow activity after the Gorkha earthquake and suggested that the transient increase in debris flow rates did not persist beyond 2016.
There are several possible reasons for the apparent differences in these two studies. One important consideration is that our modelling is based on all visible landslides in each epoch , as opposed to mapping only newly-occurring landslides (Dahlquist & West, 2019). Since we have no data on sediment availability or volume, our approach makes the conservative assumption that any existing landslide could potentially be a source of future runout. The peak in runout susceptibility after the 2015 monsoon therefore reflects the transition from more landslides being triggered than revegetated, to more landslides revegetating than being triggered, rather than an actual peak in landslide rate. It is certain that some of the landslides that persisted as bare ground in post-seismic inventories will in fact have been already exhausted of readily mobile sediment and so will no longer pose an immediate hazard in terms of secondary runout.
This means that our time series represents a precautionary scenario that will likely overestimate hazard persistence relative to studies that only document new or substantially-altered source areas. However, we know that the majority of coseismic landslide deposits remain in place on hillslopes for years after a large earthquake, with estimates for Wenchuan ranging between 80 to 90% of the initially mobilized material still being in situ after 6-7 years (Dai et al., 2021;Fan, Domènech, et al., 2018;Huang & Li, 2014;Märki et al., 2021;Zhang & Zhang, 2017). Even though the proportion of this material that is readily erodible will decrease through time (Domènech et al., 2019;Qu, 2019;Yunus et al., 2020), longer-term studies have demonstrated that elevated rates of debris flows can persist for decades Ni et al., 2019).
Similarly, enhanced rates of new post-seismic landsliding will persist in the years after a large earthquake, meaning that the source areas available for remobilization will also evolve (Dadson et al., 2004;Li et al., 2018;Zhang et al., 2016), and that new locations will become susceptible to runout (Zhang & Zhang, 2017). Defining the period over which both new post-seismic landslides and post-seismic runout play out needs to balance the likely period over which the earthquake legacy remains relevant for the failure of new landslides and the time for the landscape to experience forcing from a full spectrum of conditions that drive remobilization. Analysis of landsliding during the 2020 monsoon, which was recognized as intense in the area of the 2015 earthquake, suggests that the legacy of the Gorkha earthquake has persisted at least that long .
Importantly, our results also demonstrate that the overall potential runout extent can be estimated as a multiple of the mapped landslide area (Figure 3), thereby helping to constrain the potential magnitude of the runout phase of the overall cascading hazard. Modelled runout areas are on average four times larger than equivalent landslide areas for the pre-and coseismic periods, but five times larger for post-seismic inventories (Figures 2 and 3). This change in runout-to-source area ratio (Rout_Mod:Ls_Map) could reflect differences in the hillslope location and topographic characteristics of post-seismic landslides, which are often at higher elevations (Fan, Domènech, et al., 2018;Kincey et al., 2021) and on steeper slopes Yunus et al., 2020) than pre-seismic landslides, and so have greater potential for longer runout pathways. The similarity between the coseismic and pre-seismic runout-to-source ratios is likely the result of a combination of factors, including that coseismic landslides are typically larger, rounder in plan and, at least after the Gorkha earthquake, were less channelized than pre-or post-seismic landslides . As post-seismic landsliding progressively shifts towards a dominance of rainfall triggering, the ratio between Rout_Mod area and Ls_Map area should also decrease, a trend which may be occurring from post-monsoon 2016 (Epoch 7) onwards ( Figure 3).

| Realization timescales associated with modelled runout susceptibility
Our results show that 14% of all modelled coseismic runout cells, equivalent to 85 km 2 , experienced a landslide at some point within the 4.5 years following the earthquake (Figure 5a), and this proportion increases to 32% when the model output is thresholded by a susceptibility threshold that maximizes the model precision. We find that 90% of the total runout realization extent originates from changes to preexisting coseismic landslides, and only 10% from entirely new landslides within the modelled runout footprint ( Figure S14). This total realization area is equivalent to 57% of the total area affected by coseismic landsliding as a result of the 2015 Gorkha earthquake. Analysis of the rate at which coseismic Rout_Mod cells become Rout_Real shows a marked decrease through time, with 4.8% of the modelled runout area experiencing a landslide within 6 months of the earthquake, but just 0.3% of cells taking the full 4.5 years ( Figure 6). This indicates that the majority of actual runout from coseismic landslides occurred in the immediate post-seismic period, with limited ongoing expansion of this extent continuing through time. These results mirror the short-lived mobilization of debris flows from coseismic landslides observed by Dahlquist and West (2019) by post-monsoon 2015. We would expect that the rate at which new runout occurred from coseismic landslides will continue to decrease through time, perhaps as sediment supply decreases (Qu, 2019), fine material is preferentially removed (Domènech et al., 2019), and stabilizing vegetation becomes re-established (Shen et al., 2020;Yang et al., 2018). This is reflected in the decreasing proportion of realized runout that relates to preexisting landslides versus entirely new landslides for each of the postseismic epochs ( Figure S14).
The unrealized per-epoch percentage of coseismic Rout_Mod cells decreased consistently after the earthquake, from 81% in postmonsoon 2015 to just 50% by post-monsoon 2019 (Table S7) Ls_Map areas over time. It is important to note, however, that this pattern of change is based solely on the coseismic landslide inventory and does not include the potentially significant influence of new postseismic landsliding . This explains the apparent paradox that, although the Rout_Mod associated with coseismic landslides is decreasing through time ( Figure 6; Table S7), the overall Rout_Mod area remains high (Figure 2). New landslides provide an additional supply of sediment that can subsequently be remobilized into debris flows (Zhang & Zhang, 2017). Analysis of runout realization for landslides in the post-monsoon 2015 inventory shows that 12% (95 km 2 ) of the modelled runout area experienced a landslide by the end of the time series, with this proportion remaining above 10% through 2016 and 9% for 2017, before decreasing rapidly after this date (Figure 7). Over the same time period, 0.3% of the study area experienced landslides, of which 39% lay inside and 61% lay outside the modelled runout area, although these values include the occur- One important outstanding question is the degree to which the post-seismic runout realization timescale described earlier was determined by the timing of the earthquake itself. The Gorkha earthquake occurred in late April 2015, approximately 6 weeks before the start of the monsoon season and after a 6-month period of dry weather. This meant that the earthquake occurred in relatively dry hillslope conditions but immediately before the onset of the monsoon. It has been previously suggested that the antecedent moisture conditions at the time of an earthquake may influence the extent of coseismic landsliding (Marc et al., 2018). This has also been argued as one reason why some moderate earthquakes are able to generate more extensive coseismic landsliding, as in the case of the 2018 M w 6.6 Hokkaido earthquake, which occurred 1 day after Typhoon Jebi (Cui et al., 2021;Wang et al., 2019). Given this, the degree to which landslides run out in the post-seismic period is influenced by their behaviour and runout

| Controls on runout realization
When analysed at slope unit level, the spatial distribution of runout broadly reflects the distribution of coseismic landslides. Larger modelled runout extents are observed across much of the higher-relief areas to the north of the study area, although there is a notable region of low runout to source area ratios and high summed p(Rout) values coincident with the highest landslide densities, where the remaining hillslope area available to accommodate runout was limited (Figure 8).
Despite pronounced spatial heterogeneity in the distribution of realized coseismic runout (Figure 9a), the percentage of realized hazard is again higher in those districts most badly affected by the earthquake.
This indicates that slope units that experienced the highest densities of landsliding were also those in which post-seismic runout was most likely to occur, a finding that correlates with spatial variability in postseismic activity after the Wenchuan earthquake (Ma et al., 2017;Ni et al., 2019). Whilst this may appear somewhat unsurprising, evidence  Figure 10). This finding provides a significant step in the ability to anticipate locations that are more likely to experience post-seismic runout, and therefore where cascading hazard and associated risk is concentrated Tian et al., 2020;Zhang & Zhang, 2017).
Whilst the limited number of equivalent studies on the remobilization of coseismic and post-seismic landslides restricts a wider comparison of the controls on runout, our findings do correlate well with broader-scale assessments focusing on rates of new postseismic landsliding and the reactivation of coseismic landslides (e.g., Fan, Domènech, et al., 2018). For example, we find that runout associated with coseismic landslides is most likely to be realized at elevations > 2750 m (Figure 10a). This is to be expected given that our results are showing progressive (downslope) runout from existing landslides, which themselves typically cluster at ridge-top locations (Meunier et al., 2008). Nevertheless, it is known that new post-seismic landsliding also tends to occur at higher average elevations than preceding coseismic landsliding (Fan, Domènech, et al., 2018;Khan et al., 2013;Kincey et al., 2021;Yunus et al., 2020), and so the distributions will also partially reflect the inclusion of entirely new landslides at higher elevations. Our finding that hazard realization preferentially occurs on steep slope angles of > 35 (Figure 10b) is also reflected by studies focusing on both post-seismic debris flows (Dahlquist & West, 2019) and post-seismic landsliding Yunus et al., 2020).
Locations in the middle-upper portion of the hillslope profile were also found to be more likely to be realized than those lower down the profile, or at the ridge top (Figure 10c). This at least in part reflects the spatial distribution of coseismic landslides, which tend to cluster towards ridge tops due to topographic amplification of seismic shaking (Meunier et al., 2008), and the downslope runout from these areas to mid-slope positions during post-seismic epochs. Mid-slope positions are also likely to generate higher pore pressures and thicker overland flows, which are known debris flow triggers, as compared to ridge tops. This differs from patterns of new post-seismic landsliding after the Gorkha earthquake, which occurs more broadly across the entire hillslope profile but with a concentration towards lower hillslope positions , where high pore fluid pressure is also present and undercutting of hillslopes may be more prevalent (Densmore & Hovius, 2000). In addition, we find that realization is more likely to occur on slopes with a pronounced degree of plan or profile curvature, notably for either convex and concave slopes (Figures 10d,e), and that realization is most likely for cells with an upslope contributing area of between 40 and 15,000 m 2 m À1 (Figure 10f). Whilst inferring an underlying mechanism for concentrating landslides in these areas from these data alone is challenging, the landslide distribution may reflect a combination of flow accumulation within concavities, and a reduction in slope stability around convexities, which both reduce the local factor of safety. According to our observations, remobilization and runout concentrate where there are most likely to be notable changes to flow characteristics (i.e., local concavities or convexities). Low contributing areas correspond to ridge top positions with limited landslide material available for remobilization, while very high contributing areas correspond to valley floors or lower portions of the channel network where slope may be a limiting factor in remobilization.
Aspect-dependent asymmetry in the distribution of landslide activity has been previously documented for both coseismic (Meunier et al., 2008) and post-seismic (Fan, Domènech, et al., 2018;Kincey et al., 2021) landslide inventories, and here we show that this signal is also manifest in runout. Rout_Real is preferentially concentrated on hillslopes oriented between east, south and west-southwest ( Figure 10g), and on slopes facing away from the direction of the Gorkha epicentre (Figure 10h). These distributions may reflect a superimposition of directional asymmetry in the hillslope damage legacy generated by seismic shaking (Brain et al., 2017;Robinson et al., 2017), which influences where newly-failed ground may occur, and the dominant direction of prevailing monsoonal rains from the south-southeast in this part of the Himalaya, which combined may act to amplify the apparent directional preference for runout seen here (Fan, Domènech, et al., 2018). We also demonstrate that Rout_Mod cells proximal to pre-existing coseismic landslides, notably those within 50 m, are considerably more likely to be realized than those further away (Figure 10i). This finding indicates that runout is dominated by changes occurring over relatively short distances close to existing landslides, rather than more extensive runout over long distances, although the risk of more extensive runout should not be ignored. This pattern of activity concentrated near to existing landslides mirrors ideas around path dependency observed in other multitemporal landslide inventories (Samia et al., 2017).

| Integrating runout evolution into long-term hazard and risk management
Landslide-triggering earthquakes in mountain regions can leave a concerning legacy of unstable slopes and extensive landslide debris, which compounds the risks faced by populations whose priority is reconstruction Rieger, 2021). Any effort to understand these complex hazards, and their evolution in time, is therefore of value. Whilst our capacity for regional-scale evaluation of individual hillslope susceptibility to new landslides after an earthquake is currently extremely limited, a potentially knowable risk is that posed by runout from coseismic landslide deposits in the landscape. Such remobilization commonly and tragically can lead to significant losses in the aftermath of large earthquakes, as has been seen in central Nepal since 2015 .
In this research we argue that considerable gains can be made by isolating and characterizing the risks posed by the potential for runout from coseismic and post-seismic landslides, as this may represent a substantial portion of total landslide risk faced in the aftermath of a large earthquake. The need to model and map these risks is also clear: the long hillslopes, confined valley topography and complex drainage networks that are common in central Nepal often mean that upslope landslide hazards may not always be recognized, particularly where populations may have variable degrees of awareness of the hazards in the landscape around them. The behaviour of post-earthquake landslides is also likely to sit at odds with the lived experience of residents, because earthquakes result in larger, more numerous, and more active landslides than may have been present prior to an earthquake Rieger, 2021).
We demonstrate here that a considerable area of land (85 km 2 ) within the earthquake-affected area could have been recognized after the earthquake was at risk from runout, based upon a precautionary parameterization of Flow-R. We also show that the potential hazard associated with post-seismic runout may have a spatial extent up to five times that of the coseismic landslide footprint alone, and so must be a significant factor in post-earthquake land-use planning. Whilst only 14% of this full modelled coseismic runout extent was realized within the first 4.5 years, this represents a valuable demarcation of areas at risk that could be generated immediately after coseismic landslides have been either modelled (Robinson et al., 2017) or mapped (Williams et al., 2018). Importantly, our model was run at 10 m resolution over the full earthquake affected area, which is a resolution relevant to individual land holdings and buildings, allowing those potentially at risk to be identified.
We also show that choices involving small distances (< 50 m) can make a significant difference to exposure to runout of landslide debris, and so this fine-scale information is critical for the selection of safer places for reconstruction. The tendency for landslide debris to channelize during runout poses a particular set of risks for valley bottom settlements, which commonly occupy the deposits of previous debris flows as these areas are often the only habitable land; identifying which of these areas are at greater risk where choices for development and reconstruction are highly limited is therefore essential. An important further finding of our work is that the runout hazard following a large earthquake changes significantly and in a complex manner through time as the contributing landslide sources also evolve. Identifying both where and critically when it is safe to rebuild is equally complex.
We took a precautionary approach to the modelling of potential runout in this study, assuming that any landslide visible as bare ground in any mapping epoch represented a possible source of future runout material. This approach was justified given the absence of any meaningful data relating to sediment volumes or supply across the earthquake-affected area, but it does mean that some source locations that were already depleted of sediment will have been included within the analyses. Similarly, our approach to parameterizing and validating the model had to reflect the complexity and variability present within a landslide dataset that included > 190,000 source locations distributed across an area of > 24,000 km 2 . A useful avenue for future research would therefore be to assess in more detail how the ideal parameterization of Flow-R changes through both space and time as the population of coseismic and post-seismic landslides evolve, including consideration of appropriate model susceptibility thresholds for different hazard and risk scenarios. In particular, analysing overall realization and model performance once the rate of runout and landsliding has returned to pre-seismic levels would allow model parameters to be better refined and a less precautionary hazard assessment produced. Repeating such analyses for other multi-temporal post-seismic landslide datasets would also provide important information on the degree to which standard model calibration can be applied within different topographic and seismic contexts.
Another important avenue for future research is linking the evolution of runout potential to the changes that occur in debris flow source material, mechanisms and initiation locations highlighted by other studies (e.g., Fan, Zhang, et al., 2018;Zhang et al., 2014;Zhang & Zhang, 2017). Crucially, future work should also consider the impacts of periodic high-intensity rainfall events on determining the trajectory and timing of runout evolution, with the large-scale modelling approach presented here providing the opportunity to assess whether localized triggering factors can explain the spatial variability in patterns of reactivation and runout documented elsewhere (Tang et al., 2016;Yunus et al., 2020). Including rainfall data of sufficient spatial and temporal resolution as part of the analysis of multi-temporal source inventories and modelling of runout pathways would generate invaluable information on how runout from coseismic landslides is likely to evolve under particular environmental conditions. Integrating such information into the planning framework associated with postdisaster response therefore has the potential to substantially improve our ability to forecast evolving cascading hazards and manage the associated risks.

| CONCLUSIONS
Using a spatially-distributed empirical sediment runout model, we considered how the spatial extent and relative likelihood of potential runout from existing landslides changed in the 4.5 years following the 2015 Gorkha earthquake in Nepal. Our results indicate that runout from coseismic landslides represents a considerable component of the overall mountain hazard chain, with actual runout representing an area equivalent to 57% of the total area impacted by coseismic landsliding, and the modelled potential runout area being on average 4-5 times the equivalent coseismic landslide area. Although the modelled runout area from the coseismic landslide population decreased through time, the overall runout potential remained high, indicating that the runout hazard is changing as a result of the evolution of the post-seismic landslide distribution itself. This finding clearly demonstrates the importance of developing systematic multi-temporal landslide inventories and associated runout susceptibility assessments in the years after a large earthquake.
Predicting the precise timing of runout from existing landslides remains problematic in the absence of high-resolution and accurate precipitation forecast data, an in-depth understanding of antecedent and local hillslope conditions, and extensive early warning systems.
However, our results demonstrate the possibility to anticipate the spatial extent of future runout across an entire earthquake-affected area and to provide indicative timescales over which the runout is likely to occur. A comparison of modelled runout extents with subsequent mapped landslides shows that 14% of all modelled coseismic runout cells became landslideswhich we term 'realization' of the hazardat some point during the 4.5 years after the earthquake, equivalent to 85 km 2 of newly-affected ground. Limiting the modelled runout extent to higher susceptibility areas increases the model precision to 32% (10 km 2 ), meaning that the spatial location of a considerable area of potential future runout risk can be identified immediately after an earthquake using our precautionary approach to the modelling of secondary hazards.
Our analysis shows that the majority of runout realization occurs within the first 12 months after the earthquake, but that runout activity still persists after 4.5 years, reflecting both a lag in the subse- Our findings enhance understanding of the extent and timing of cascading hazards following high magnitude earthquakes in mountain regions. Such regional-scale modelling of runout susceptibility from existing landslides has the potential to refine prediction of where, and over what timescales, future runout may occur, thereby greatly improving our ability to inform how we manage long-term postseismic cascading hazard and risk. Embedding this knowledge within frameworks of disaster planning and decision-making has the potential to significantly improve the effectiveness of post-event recovery and reconstruction.