Evolution of Coseismic and Post‐seismic Landsliding After the 2015 Mw 7.8 Gorkha Earthquake, Nepal

Coseismic landslides are a major hazard associated with large earthquakes in mountainous regions. Despite growing evidence for their widespread impacts and persistence, current understanding of the evolution of landsliding over time after large earthquakes, the hazard that these landslides pose, and their role in the mountain sediment cascade remains limited. To address this, we present the first systematic multi‐temporal landslide inventory to span the full rupture area of a large continental earthquake across the pre‐, co‐ and post‐seismic periods. We focus on the 3.5 years after the 2015 Mw 7.8 Gorkha earthquake in Nepal and show that throughout this period both the number and area of mapped landslides have remained higher than on the day of the earthquake itself. We document systematic upslope and northward shifts in the density of landsliding through time. Areas where landslides have persisted tend to cluster in space, but those areas that have returned to pre‐earthquake conditions are more dispersed. While both pre‐ and coseismic landslide locations tend to persist within mapped post‐earthquake inventories, a wider population of newly activated but spatially dispersed landslides has developed after the earthquake. This is particularly important for post‐earthquake recovery plans that are typically based on hazard assessments conducted immediately after the earthquake and thus do not consider the evolving landslide hazard. We show that recovery back to pre‐earthquake landsliding rates is fundamentally dependent on how that recovery is defined and measured. Clarity around this definition is particularly important for informing a comprehensive approach to post‐earthquake landslide hazard and risk.

process, however, with the result that much of our existing understanding of post-seismic landsliding comes from just two earthquakes: the 1999 M w 7.6 Chi-Chi earthquake in Taiwan, and the 2008 M w 7.9 Wenchuan earthquake in China. The Chi-Chi earthquake triggered >20,000 landslides and elevated suspended sediment concentrations in rivers draining the affected area by a factor of four relative to decadal background levels (Dadson et al., 2004). Analysis of stream gauge data and time series landslide mapping following this event suggested a post-seismic decay in suspended sediment flux (Hovius et al., 2011) and the rate of new landsliding  that returned to pre-event levels within approximately 6 years.  extended this analysis to three other shallow thrust-faulting earthquakes, demonstrating that the rate of new landslide occurrence peaked immediately after each event and then decayed to pre-event values within 1-4 years, with the recovery time tentatively suggested to be proportional to earthquake magnitude. Khan et al. (2013) and Barth et al. (2019) suggested similarly short recovery times after the 2005 M w 7.6 Kashmir and 2012 M w 7.8 Haida Gwaii earthquakes, respectively. This recovery time is a critical factor in post-seismic landslide hazard, because it determines the time period over which persistent and repeated disruption might be expected across the earthquake-affected area.
In contrast, studies of the Wenchuan earthquake have suggested a more complex response, with a wide range in inferred recovery times. The earthquake triggered >60,000 landslides over an area of >100,000 km 2 (Dai et al., 2011;Xu et al., 2014). A large and growing number of studies have assessed spatio-temporal changes in the occurrence of new landslides as well as the remobilization and transport of pre-and coseismic landslide debris. These studies have used a variety of approaches, but have tended to focus on small portions of the overall landslide-affected area. For example, Fan et al. (2018) ;Fan, Scaringi, Domènech, et al. (2019); and Fan, Scaringi, Korup, et al. (2019) compiled four post-earthquake inventories over a 463 km 2 area near the epicenter (∼0.4% of the total affected area). New post-seismic landslides made up less than 10% of each of their inventories; instead, post-earthquake landslide activity was dominated by remobilization of existing landslide scars and their deposits into debris flows. Both the proportion of landslides that were remobilized and the extent of remobilization declined over time, returning to background levels within approximately seven years. Similarly, Tang et al. (2016) compiled six inventories spanning 2005-2015 over a 172 km 2 area near the epicenter (∼0.1% of the affected area), demonstrating a 100-fold decrease in the number of active landslides in the 7 years after the earthquake. While most post-earthquake landsliding occurred in the first three years, landslide rates had still not decayed to pre-earthquake levels by the end of their study. In contrast, other studies have argued that the rates of landsliding (e.g., R. Huang & Li, 2014) or post-seismic debris flows (e.g.,  are likely to remain elevated for several decades, with C. Li et al. (2018) speculating that a full return to pre-earthquake conditions could take 50-100 years. Studies of vegetation damage and recovery, using satellite image-based vegetation proxies (e.g., Chen et al., 2020;Ni et al., 2019;Shen et al., 2020;Yang et al., 2018;Yunus et al., 2020), have also tended to infer longer recovery times, on the order of several decades, but have pointed to strong spatial heterogeneities in recovery, with time scales expected to be longer in the areas with the highest density of coseismic landslides (Ni et al., 2019).
These disparities in inferred recovery time likely stem from two separate sources: (1) real differences in recovery time between different earthquakes, which we would like to understand, and (2) differences in definition, methodology, and resolution, which we would like to correct for. Even the term "recovery" is used in a range of ways, from a decline in the rate of new landslides or debris flows to pre-earthquake values, to stabilization of coseismic landslide debris, to the re-establishment of pre-earthquake vegetation. There is little consensus on whether to focus solely on the occurrence of new post-earthquake landslides on previously unfailed hillslopes or on the continued evolution of pre-and coseismic landslide scars. For example, landslide rates following the 2005 Kashmir earthquake differ markedly depending on the specific method used to generate the inventory, with full recovery estimates ranging from <5 years (Khan et al., 2013) to >13 years (Shafique, 2020). Inventories that contain only new landslides can show a rapid post-earthquake decline in landslide rates (e.g., , but will underpredict the hazard from pre-seismic or coseismic landslides that remain active. Identifying landslides that remain active but which have not changed overall shape or size in a distinct, measurable way is problematic without detailed field surveys, often resulting in only broad definitions of activity. There is evidence, however, that the hazard posed by landslides can persist, even when the overall mapped footprint does not change. For example, in their field-based study of landslide development along two critical highways in Nepal in the three years following the 2015 M w 7.8 Gorkha earthquake, Tian et al. (2020) documented repeated activity and continued hazard across KINCEY ET AL. the majority of sites investigated. Importantly, this included landslides that had not changed in overall size but had either reactivated or had undergone a change in dominant failure mechanism, such as from dry debris cascades to wet debris flows. Increasing revegetation of a coseismic landslide was also not necessarily indicative of a cessation of hazard, with five revegetating landslides noted as still active three years after the earthquake. Multi-temporal inventories that focus solely on new or substantially altered landslide footprints may potentially therefore underestimate the overall level of landslide hazard, and potentially overestimate the rate of landscape recovery.
In addition, many studies have focused on only small portions of the total earthquake-affected area, despite growing evidence of strong spatial variations in both post-seismic landsliding and recovery (e.g., Shen et al., 2020;Yunus et al., 2020). Differences in methodology may play a role as well; airborne and satellite images are well-suited for distinguishing between vegetated and unvegetated areas, but not for determining the stability of the ground surface or the extent of landslide activity and downslope transport (Bernard et al., 2020). Indirect measures of landslide activity, such as downstream suspended-sediment discharge (e.g., Dadson et al., 2004;Wang et al., 2015;F. Zhang et al., 2019), may reflect changes in sediment supply or transfer to rivers rather than changes in landslide activity. Comparison of such direct and indirect measures relies on a comprehensive understanding of sediment transport pathways and storage within the landscape which is often lacking. Finally, it is also important to recall that recovery of landslide rates to pre-earthquake levels is not synonymous with the removal of all coseismic landslide debris from the landscape. Careful documentation of debris volumes and transport in portions of the Wenchuan earthquake rupture area has demonstrated that >80% of the material liberated by coseismic landsliding continued to reside on hillslopes or in low-order channels after a decade (Fan et al., 2018;R. Huang and Li, 2014;S. Zhang & Zhang, 2017). This suggests that the response to large earthquakes evolves over multiple overprinting time scales: relatively short-term changes in the rate of new and persistent mass wasting on the one hand, and long-term enhancement of erosion rates and sediment flux due to hillslope damage and sustained high rates of sediment supply on the other (Fan et al., 2018).
Together, these issues mean that we do not yet have a comprehensive picture of how landslide inventories evolve in space and time after a large earthquake, despite the importance of this picture for understanding dynamic landslide hazard and the associated risk. To begin to address this gap, we document coseismic and post-seismic landsliding across the area that was significantly impacted by the 2015 M w 7.8 Gorkha earthquake in Nepal over the period 2014-2018. We use a consistent mapping methodology, based on preand post-monsoon medium-resolution satellite imagery, that is designed to optimize our ability to assess change through time. Critically, we map each of the 11 time epochs independently and in full, resulting in an internally coherent set of 11 inventories that each include all visible landslides at each point in time. As a result, our approach deliberately accounts for the persistence and evolution of pre-existing landslides, rather than just the occurrence of new failures, providing an alternative perspective on the rate of post-seismic landscape recovery. This allows the first documentation of both the evolution of the total distribution of landslide activity and detailed shifts in landslide form and location resulting from the earthquake and the four subsequent monsoon seasons.

The 2015 M w 7.8 Gorkha Earthquake and Landslide Inventories
The 25 April 2015 M w 7.8 Gorkha earthquake initiated ∼80 km northwest of Kathmandu in the district of Gorkha, with rupture propagating eastwards for ∼140 km along the Main Himalayan Thrust (Avouac et al., 2015). A series of large aftershocks followed, including a M w 7.3 event on 12 May ∼75 km east-northeast of Kathmandu in Dolakha district ( Figure 1). Intense shaking triggered extensive landsliding across the east-west extent of fault rupture, with the density of landslides generally increasing toward the east (Martha et al., 2017). Multiple coseismic landslide inventories have been produced for the Gorkha event (e.g., Kargel et al., 2016;Martha et al., 2017;Meena & Tavakkoli Piralilou, 2019;Roback et al., 2018;Tiwari et al., 2017;Williams et al., 2018;Xu et al., 2016), resulting in estimates of landslide numbers ranging from <5,000 to ∼25,000. The most detailed assessment, and therefore probably the most complete inventory to date, was produced from high-resolution satellite imagery by Roback et al. (2018) and documented ∼25,000 landslides with a combined area of ∼90 km 2 , distributed across an area of 28,300 km 2 of central Nepal.
Assessments of post-seismic landslide rates and distributions are so far fairly limited for the Gorkha earthquake sequence. Tian et al. (2020) used imagery and field surveys to document the evolution of c. 30 landslides along two major highway corridors. They found that most of the landslides enlarged or remained active through the first three years after the earthquake, with frequent remobilization of pre-and coseismic landslide material into debris flows. Marc et al. (2019) suggested that rainfall-induced landslide rates in three selected catchments covering 7% of the total affected area were ∼3-6 times higher for the 2015 monsoon when compared to pre-earthquake levels. Between 70% and 80% of the landslides in their 2015 monsoon inventory were in new locations when compared with the coseismic landslides from Roback et al. (2018), suggesting that new failures occurred on hillslopes weakened by the earthquake rather than by reactivation of pre-existing landslides. By 2017, the rate of new rainfall-induced landslides in all three catchments had returned to pre-earthquake levels (Marc et al., 2019). Enhanced rates of post-seismic debris flows also appear to have been relatively short-lived, with the exhaustion of loose material mobilized from coseismic failures during the 2015 monsoon and a reduction in the number of new debris flows back to pre-earthquake levels by 2016 (Dahlquist & West, 2019). Dahlquist and West (2019) reported that only ∼2% of coseismic landslides developed into post-seismic debris flows, with the majority of failures likely to either KINCEY ET AL.   (Credit: AW3D30-JAXA). Stars show the main epicenters of the 2015 Gorkha earthquake sequence: 25 April M w 7.8 and 12 May M w 7.3. Gray polygons show areas of persistent cloud or snow cover throughout the study period, which were excluded from the analysis. Red box shows location of panel (c). (b) Spatial distribution of coseismic landslides (epoch 4) mapped across the 14 administrative districts that were most intensively affected by the earthquake. (c) Example of multitemporal landslide inventory results for an area of the upper Bhote Kosi river basin in Sindhupalchok. Note complex patterns of overlap between polygons from different inventories; this complicates spatial correlation and analysis of regional-scale evolutionary patterns, and motivates our focus on area density estimates. ALOS, Advanced Land Observing Satellite; DEM, digital elevation model. be in more stable hillslope positions (Dahlquist & West, 2019) or directly coupled with fluvial channels and so exhausted of loose sediment immediately following the earthquake (Roback et al., 2018).

Landslide Inventory Mapping
We mapped landslides from a time series of freely available medium-resolution satellite imagery (Landsat and Sentinel-2) covering the period from 2014 to 2018, and focusing on an area of 25,575 km 2 covering the 14 administrative districts that were most intensively affected by the 2015 Gorkha earthquake, as identified by the Government of Nepal ( Figure 1) This represents just over 90% of the landslide-affected area mapped by Roback et al. (2018). Mapping was divided into 11 individual epochs, including pre-and post-monsoon inventories for each individual year, as well as an additional coseismic inventory for 2015 (Table 1). Landsat 8 imagery was used for epochs 1-5 (2014-2015), with the lower resolution (30 m) multispectral bands being pan-sharpened to 15 m resolution using the panchromatic band 8. Panchromatic sharpening involves the radiometric transformation of a lower spatial resolution multiband image using a higher resolution image band, effectively increasing the spatial resolution of the original color image to provide enhanced detail. We then switched to Sentinel-2 imagery (10 m spatial resolution) when coverage became available from epoch 6 (pre-monsoon 2016) onwards. Medium-resolution imagery was chosen because it spans the full pre-and post-earthquake time period, allows for low-cost long-term monitoring over the entire earthquake-affected area, and suffers less geometric distortion compared to high resolution imagery, which is a particular problem in high-relief environments (Weiss & Walsh, 2009). For each epoch, we chose images that spanned the shortest possible time window while maintaining a consistent time series of mid-point dates ( Table 1).
Identification of individual landslides was based on visual interpretation of true color (R-G-B) and false color (NIR) composites using bands 2-3-4-5 for Landsat 8 and 2-3-4-8 for Sentinel-2. Landslides were typically identified by the spectral contrast between exposed sediment or bedrock within the failure and the surrounding, undisturbed vegetation. Mapping of identified landslides in each epoch was conducted within QGIS v2.18 by two trained mappers using a 5 × 5 km fishnet grid to guide the mapping process. All landslides were manually digitized as vector polygons. The resolution of the satellite imagery meant that it was not possible to distinguish between failure source areas and deposits (cf. Roback et al., 2018). We characterized landslides by their planimetric area and by the ratio of area to perimeter (A/P), which indicates the degree of elongation; more elongate landslides have lower A/P values than those that are more equant. We did not convert area to volume due to uncertainties in both scaling laws (e.g., Larsen et al., 2010), potential aggregation of multiple polygons (e.g., G. Li et al., 2014;, and the wide variety of landslide types triggered by the event (Kargel et al., 2016 Unlike previous multi-temporal inventory assessments (e.g., Fan et al., 2018;Fan, Scaringi, Domènech, et al., 2019;Fan, Scaringi, Korup, et al., 2019;C. Li et al., 2018;Marc et al., 2019), we mapped the full extent of all landslides visible on imagery from each epoch, treating each epoch independently. This meant that we mapped landslides based solely on visibility in the imagery and irrespective of whether they had been documented during earlier epochs, rather than only recording landslides which had newly occurred or substantially changed since the preceding epoch. This approach was selected for a number of reasons. Focusing only on new landslides has the potential to underestimate any persistent hazard posed by previously failed slopes which have not yet fully stabilized or revegetated. In contrast, mapping all visible landslides from each epoch provides a measure of the total area of impacted ground for any given date, as well as allowing the overall rates of hillslope recovery and revegetation back to pre-earthquake levels to be derived. Understanding the degree to which the landslide hazard footprint persists or alters in the years following a large earthquake is crucial for determining a realistic measure of spatially and temporally variable hazard and risk. Our approach also avoids the need to classify the extent of reactivation of each landslide based on qualitative criteria (e.g., Fan et al., 2018;C. Li et al., 2018). With the medium-resolution imagery used here, it is not possible to observe variations in reactivation or remobilization of debris within individual landslide scars. Thus, our hillslope recovery rates (in terms of landslide areal density) are expected to be slower compared to studies that examine only new or clearly reactivated landslides in each epoch.
To assess the impact of this conceptual approach on our derived recovery rates, we also applied an alternative methodology based on the analysis of individual 10 × 10 m cells rather than entire mapped polygons. In this approach, every 10 × 10 m cell within the study area was assigned a series of binary values showing whether it was a landslide ("1") or not ("0") within each of the mapping epochs. Cells that had been newly impacted relative to previous epochs, whether through the occurrence of entirely new landslides or expansion of existing landslide footprints, could then be isolated from persistent landslide cells and analyzed separately in terms of their distribution and recovery rate. Additional information relating to this alternative methodology is provided in supporting information 7, while a detailed consideration of the different conceptual approaches and the implications for assessing post-seismic hillslope recovery is discussed in Section 5.1.
A common issue with large-scale multi-temporal landslide inventory analysis is the presence of intermittent and spatially variable cloud, shadow, and snow cover, which typically precludes the direct comparison of data from affected areas. Automated masking of cloud, shadow, and snow cover was undertaken using the Fmask (Function of mask) software v4.0 (Qiu et al., 2019), with manual checking and editing of the resulting mask layers to ensure accuracy. Areas of persistent cloud or snow (i.e., those for which the ground was not visible in any of the 11 epochs) totaled 1308.32 km 2 , 5% of the total study area, and were concentrated in the High Himalaya in the north of the study area. These areas were excluded from the analysis, leaving a "seen ground" mapping extent of 24266.68 km 2 .
A second issue with multi-temporal inventories is intermittent visibility-that is, landslides which are visible in some epochs but obscured in others, due to cloud cover, shadowing, or poor image contrast. To address this issue, we developed a novel approach of forward propagation of landslide polygons from successive epochs based on a series of structured spatial overlap queries. Any landslides visible in both epoch n−1 and epoch n+1 that were obscured in epoch n, and that did not intersect a pre-existing epoch n polygon, were merged from epoch n−1 into the epoch n inventory. This process was extended to also extract landslides visible in epoch n−1 and epoch n+2 that were obscured in both epochs n and n+1. An additional spatial query identified landslides missed due to relief-based shadowing, a problem that was particularly acute on north-facing high-relief slopes in the study area. This involved selecting landslides from epoch n−1 that were within a one-pixel distance (with 15 m as the coarsest pixel resolution) of an epoch n+1 landslide but were more than one pixel distant from any landslide present in epoch n, and then merging these selected landslides into the epoch n inventory.
No landslides were propagated forward more than two epochs, equivalent to a maximum time interval of 12 months, and propagated landslides were tagged as inferred rather than directly mapped. Results were manually checked to ensure that the propagated landslides were correctly identified, and any erroneous polygons were removed. The number of landslides propagated forwards varied between epochs due to the differing extents and distributions of cloud and snow cover (supporting information 1). An average of 940 landslides were propagated into each mapping epoch, representing a mean of 8% of each inventory. The maximum number of inferred landslides was for epoch 3 (pre-earthquake 2015) when cloud cover was highest, with 2,094 landslides being propagated from the two previous epochs. This equates to 33% of the epoch 3 inventory but this is proportionally higher than later epochs, due to the relatively low numbers of landslides in the pre-earthquake inventories (see Table S1 for full inventory statistics).

Mapping Uncertainties
Errors occur in all manually mapped landslide inventories and relate to issues such as image clarity and the interpretative skills and experience of the mappers (Fan et al., 2018;Meena & Tavakkoli Piralilou, 2019). Such errors were minimized for this study through (1) initial training of the mappers by repeat delineation of landslides within a 2,500 km 2 area that spanned a range of physiographic settings and landslide densities to identify and then correct individual mapper bias (see supporting information 2), and (2) manual checking and editing of each completed inventory by a team of independent analysts. Each 5 × 5 km fishnet cell was assessed by a minimum of one mapper and by two independent analysts. This validation process involved the visual checking of each landslide polygon against the base imagery for that particular epoch. Landslides that had not been identified during the initial mapping were added and features that were determined to have been mistakenly mapped as landslides (such as quarries) were deleted from the inventory. Particular care was taken to ensure that any amalgamated landslides were segmented into individual polygons representing discrete failures, in order to provide accurate area-frequency distributions for each epoch .
Direct comparisons were made between the coseismic inventory from this current study and that produced by Roback et al. (2018) using higher-resolution satellite imagery (supporting information 3). We also quantified the influence of sensor type and spatial resolution that might result from our two imagery types on landslide area characteristics through repeat mapping of a 121 km 2 subset area using a range of base imagery (supporting information 4). This involved the repeat blind mapping of landslides visible on a Pleiades image tile from September 2015 (coincident with epoch 5) which had been resampled to different spatial resolutions (2, 10, and 15 m), as well as equivalent mapping from a concurrent Landsat 8 tile which had been pan-sharpened to 15 m resolution.

Landslide Area Density Analyses
Direct analysis of multi-temporal landslide polygon inventories is potentially complicated by small variations in polygon boundaries and overlaps between epochs, which can frustrate attempts to correlate polygons belonging to individual landslides. Without reliable correlation across epochs, it is difficult to establish the evolutionary trajectory of landsliding at a regional scale. An alternative means of monitoring temporal changes in landslide inventories is to divide the landscape into a grid of cells, and to consider the epochs in which cells first become affected by landsliding and then become inactive again once the hillslopes have healed and revegetated to the point that no potentially active landslides are visible. This framework can then be used to assess when each cell first and last experiences a landslide, the duration of activity, and the spatial extent of landslide activity at any given point in time.
We chose a cell size of 1 × 1 km as a compromise between retaining the spatial fidelity of the original mapped landslide inventories and detecting broad-scale spatio-temporal trends in landsliding across the entire affected area. For each cell, we assigned a value for (1) the first epoch in which that cell was intersected by one or more landslides, (2) the last epoch in which that cell was intersected, (3) the total number of epochs for which the cell was impacted by landsliding, and (4) the per-epoch landslide area density, expressed as a percentage of the cell area. These metrics enabled an assessment of the spatio-temporal variability in landslide area density and the distributions of density values recorded in our inventories, as well as providing insights into potential controls on the spatial distribution and timing of when particular cells became active (here termed the landslide "birth" epoch) or inactive (the landslide "death" epoch).

Physical Controls on Landsliding
A range of topographic and distance-based variables were summarized by individual landslide polygon, and each polygon was also attributed with the birth and death results associated with the specific 1 km 2 grid cell in which it was located. This approach allowed the different control variables to be analyzed for their potential influence on spatial and temporal variability in landslide birth and death. Kernel density estimations were generated for each variable and by each birth/death category, with the results then being differenced from equivalent density estimations for the overall study area.
Topographic information was sourced from the Advanced Land Observing Satellite World 3D-30m (AW3D30) data set, a freely available global digital elevation model (DEM) at 30 m spatial resolution. Although slightly coarser than the mapping resolution, this DEM represents the highest-resolution openly available data with full coverage of the study area, and has been shown to be particularly suitable for modeling of topographic derivatives in mountainous terrain (Boulton & Stokes, 2018). We focused on derivatives that have previously been shown to relate to coseismic landsliding: slope, aspect, plan curvature, profile curvature, topographic wetness index, upslope contributing area, and hillslope relief (Parker et al., 2015;Robinson et al., 2017). The relative position of a landslide along the hillslope profile was calculated by measuring the minimum distance from the polygon perimeter to the nearest ridge top, to capture the position of the head scarp, and to the nearest river, to define the position of the toe. Distance values were then normalized to 0-1, with 0 representing the river and 1 the ridge. Earthquake-specific layers were also calculated, including distance to epicenter and slope direction to epicenter. The slope direction to epicenter was only calculated relative to the 25 April mainshock epicenter in Gorkha, rather than to both the mainshock and 12 May aftershock epicenters. This was because previous studies have demonstrated that the majority of landsliding associated with the Gorkha earthquake sequence was associated with the 25 April mainshock, with only a small number of new landslides triggered by the 12 May aftershock (Roback et al., 2018).
Rainfall data covering the study period were analyzed to highlight any variation in precipitation trends that could potentially have contributed to temporal variability in landslide rates. Gridded monthly precipitation totals at a spatial resolution of 0.1° (11 × 11 km) were obtained from NASA's Global Precipitation Measurement (GPM) Integrated MultisatellitE Retrievals for GPM (IMERG) data set, covering the period from 2014 to 2018 (supporting information 5). Individual NetCDF files were converted to geotiffs and then stacked to form a multiband raster containing 60 bands representing monthly precipitation totals. Mean monthly rainfall totals for the entire study area were then calculated by averaging the cell values from each raster band.

Results
The spatial pattern of landsliding in epoch 4 (coseismic) is shown in Figure 1b, while the detailed temporal variability between inventories is shown in Figure 1c for a small section of the study area. The full inventory includes 160,332 individually mapped landslide features across the 11 epochs. Inventory data for the whole area of interest are described in Figure 2 and Table S1 in terms of number and area, in Figure 2a in terms of areal density, and compared to the timing and intensity of the monsoon in Figure S8. Because of the complexity of these patterns, we focus first on aggregate statistical descriptions of our inventories, before turning to the distribution and evolution of landslide area density estimates.

Multi-temporal Landslide Inventory Statistics
Landslide numbers are consistently low for each of the three pre-seismic inventories, averaging ∼6,400, with little discernible difference between pre-and post-monsoon inventories ( . Importantly, by the end of the study period the number of visible landslides remained 4,908 greater than were recorded in the coseismic inventory (13,684) and 12,193 greater than the pre-earthquake inventory average.
Total landslide area followed broadly the same pattern as landslide count, with consistency before the 2015 earthquake (mean = 64.63 km 2 , σ = 1.71 km 2 ) and sharp increases immediately following the earthquake (154.68 km 2 ) and the 2015 monsoon (179.03 km 2 ) ( Figure 2a). Landslide area decreased steadily in 2016-17, albeit with pronounced fluctuations between pre-and post-monsoon epochs, with pre-monsoon inventories being ∼22 km 2 (∼16%) higher on average than their subsequent post-monsoon inventories. While this appears large, this variability averages less than c. 30 m 2 /km 2 across our study area. In contrast, areas in both pre-and post-monsoon 2018 inventories were similar, and still exceeded the area recorded by the immediate post-earthquake inventory (epoch 4).
The area-frequency distribution of landslides in the coseismic inventory is markedly different to the preand post-seismic inventories (Figure 2b). Prior to the earthquake, the mean landslide area was ∼10,000 m 2 with an A/P ratio of 14, and these values were relatively invariant across the three epochs. On average, coseismic landslides were both larger and less elongate than landslides before and after the earthquake, with a mean area of 11,300 m 2 and an area/perimeter (A/P) ratio of 15.3. Following the earthquake, both average area and A/P ratio declined through post-monsoon 2016 and then stabilized at values that were substantially lower than pre-earthquake levels-indicating a systematic shift to both smaller and more elongate failures than were observed before the earthquake.
Comparison with the independent inventory compiled by Roback et al. (2018) shows that our inventory contains fewer small landslides (areas <10 3 m 2 ) and an overall shift in the area distribution toward larger failures, consistent with our use of coarser resolution imagery (supporting information 3; Figure S4). Our coseismic inventory exhibits a characteristic probability density distribution with a rollover between 1,000 and 1,500 m 2 and power-law scaling with an exponent of 2.93 beyond a size threshold of 48,000 m 2 (Table S2). In comparison, the portion of the comparable inventory compiled by Roback et al. (2018) that overlaps with our study area has a rollover between 100 and 300 m 2 and a power-law exponent of 2.74 beyond a threshold area of 18,000 m 2 . These differences are primarily due to the base imagery used for the mapping in each case, with our coseismic inventory being generated from 15 m resolution pan-sharpened Landsat 8 imagery, compared with the 0.2-0.5 m resolution Worldview imagery used by Roback et al. (2018). Repeat mapping of the 121 km 2 subset area outlined above (Section 3.2) shows a similar change in the probability density distributions when progressively coarsening image spatial resolution from 2 to 15 m (supporting information 4; Figure S6). The overall distribution is shifted toward proportionally larger failures as image resolution is degraded and smaller failures become undetectable, with a corresponding increase in the average size of mapped landslides and decrease in the overall number of landslides in the inventory ( Figure S7). Power-law scaling of these subset inventories broadly indicates that the minimum landslide size threshold for the power-law fit increases as image resolution is degraded, although with pronounced variability which is likely the result of a relatively small mapping area and inventory sample size (Table S3).

Evolution of Landslide Area Density Distribution
When averaged using a 1 × 1 km grid, the overall landslide area density increased markedly as a result of the earthquake. The per-epoch mean area density followed the same pattern as total landslide area: increasing to 0.62% after the earthquake and 0.72% after the 2015 monsoon, declining through 2016-2017, and remaining higher than the coseismic value after the 2018 monsoon ( Figure 2a). The spatial distribution of landslide area density provides a useful visualization of the pattern of pre-seismic, coseismic, and post-seismic landsliding (Figure 3a). The landscape-scale impact of the earthquake is clearly visible as a widespread increase in landslide area density at epoch 4 across the entire study area, focused in particular along a broad northwest-southeast band following the physiographic divide between the Lesser Himalaya/Middle Hills and the High Himalaya ( Figure 3a). In each of the pre-earthquake epochs, 15% of grid cells within the study area were impacted by landslides, increasing to 24% following the earthquake then to 29% following the KINCEY ET AL.
10.1029/2020JF005803 10 of 27 Pairwise differencing between individual epochs shows that changes in the spatial pattern of landslide density through time were highly spatially correlated (Figure 4a). The change between pre-seismic and coseismic inventories (epochs 3-4) clearly shows the focusing of landslide activity within the major valleys draining the High Himalaya, particularly within Gorkha, Sindupalchok, Rasuwa, Dhading, and Dolakha districts (see Figure 1 for district names). Between the earthquake and the 2015 monsoon (epochs 4-5) there was a systematic shift in landslides toward the higher-elevation regions of the northerly districts of the study area. Within those areas, however, there were also clusters of reduced landslide density, with more areas experiencing a reduction in landslide density (16% of the study area) than an increase (13%) between post-monsoon 2015 and pre-monsoon 2016 ( Figure 4b). As expected from the temporal changes in landslide area (Figure 2), area densities tended to decrease during subsequent monsoons and increase between the end of one monsoon and the start of the next. In epochs 9-10 and 10-11, we observe secondary increases in landslide density through the central and north-central parts of the study area, including the northern part of Kabhrepalanchok district-well outside of the main northwest-southeast band of intense landsliding.
Overall, landslide area density across much of the area was still markedly higher in post-monsoon 2018 than it was immediately prior to the 2015 earthquake, especially within the major river valleys that flow from the High Himalaya (Figure 5a). Over 25% of the study area had an increased density relative to pre-earthquake conditions, compared with only 7% with a decreased density (Figure 5c)  more spatially complex (Figure 5b), with approximately the same proportion of cells having experienced increased (19%) and decreased (16%) landslide densities (Figure 5c). Much of the central northwest-southeast band of intense coseismic landsliding experienced a decrease in density over this time period, with the greatest changes concentrated on the lower slopes of major river valleys. At the same time, there are substantial zones of increased post-seismic landsliding in western Sindhupalchok and northern Nuwakot, as well as the higher ground across the northern margin of the study area. A more dispersed area of increased landslide density is also present within the Lesser Himalaya in the southeast of the study area, stretching across Okhaldunga, Ramechhap and northern Sindhuli districts (Figure 5b). Qualitative observations from districts recording pronounced post-seismic changes demonstrate a range of commonly observed behaviors, including the occurrence of new post-seismic landslides, expansion of pre-existing landslides through both retrogression and downslope channelized runout, and recovery of coseismic landslide scars; some of these behaviors can be seen in the detailed changes visible in Figure 1c.

Spatial and Temporal Distribution of Landslide Birth and Death
Categorizing grid cells by landslide "birth" and "death" results in nine discrete classes ( Figure 6). For example, there were three classes associated with landslides that were triggered by the 2015 earthquake ("coseismic birth"), distinguished based on whether the landslides were no longer visible after the coseismic inventory ("coseismic death"), disappeared during one of the post-seismic epochs ("pre-2018 monsoon death") or were still visible in the final mapping epoch ("post-2018 monsoon survival"). The distribution of these classes is spatially correlated, with clustering of areas that have recovered since the earthquake and those where new post-earthquake landslides have occurred (Figure 6a). Grid cells in which the landslide hazard persisted to the end of the study period-that is, areas where landslides were still visible in the post-monsoon 2018 inventory-dominate the distribution. Those cells where landslides were born before the earthquake are distributed throughout the study area, although somewhat clustered in the high-relief areas to the north and also the foothills to the south and east of the Kathmandu valley. In contrast, cells where landslides were born during the earthquake are largely restricted to the main northwest-southeast band of intense landsliding between the two epicenters. Cells where landslides were born after the earthquake are particularly concentrated in the valleys peripheral to the coseismic distribution, with a notable shift northward toward areas of higher elevation and the upstream reaches of the drainage network.
To assess the extent of clustering, we considered cells in each category in turn, and examined the distribution of categories within the eight cells adjacent to each cell. Clustering is noticeably greater for the three classes containing landslides which persisted to the end of the study period than for those which recovered before that time (Figure 6b). For example, 34% of cells adjacent to those which were activated in the earthquake and which still contained mapped landslides in post-monsoon 2018 are from the same class, compared with only 21% for cells which were coseismically activated but which had recovered prior to post-monsoon 2018.
In general terms, cells which recovered at any point during the time series were more spatially isolated and tended to occur around the periphery of clusters of cells containing persistent landsliding.
More grid cells were initially activated during the three pre-earthquake mapping epochs (17% of the total number) than were activated by the earthquake (9%). Importantly, this observation indicates that, although the density of landsliding caused by the earthquake was much higher (Figure 3), the spatial area impacted was comparatively confined (Figure 6c) were first impacted by landsliding during a post-earthquake epoch, with 73% of these still containing one or more landslides by post-monsoon 2018.
For the majority of epochs, the proportion of grid cells newly affected relative to the preceding epoch was low (1%-2%), with the exception of the earthquake (10%), the 2015 monsoon (6%), and pre-monsoon 2018 (3%) (Figure 6d). About 11% of grid cells in the post-monsoon 2018 epoch had been newly affected relative to the earthquake distribution, indicating that the spatial pattern of landsliding had altered steadily but progressively since the earthquake.

Topographic Controls on Regional-Scale Landslide Birth and Death
Of the 10 topographic derivatives tested, four variables showed clear skill in explaining landslide birth and death: elevation, slope angle, slope direction to epicenter, and normalized distance from river to ridgeline. The other topographic variables showed either no consistent relationships or no clear difference between birth and death epochs and so were omitted from further analysis. We focus here on understanding potential controls on patterns of landslide activation and recovery post-earthquake, and we note that more 10.1029/2020JF005803

Elevation
Landslides within grid cells that were first activated during pre-earthquake epochs preferentially occurred at elevations between 2,500 and 5,000 m, representing the elevation range of the intermediate hillslopes covering much of the High Himalaya to the north of the study area (Figure 7a). In contrast, the distribution of coseismically activated cells was shifted toward lower elevations between 1,500 and 4,000 m. This corresponds with the band of intense seismic shaking caused by the 2015 earthquake and correlates with other studies of the coseismic inventory (e.g., Roback et al., 2018). The elevation distribution for landslides occurring within cells that were first activated after the earthquake shows a shift back toward higher elevations, with a modal peak in the density difference at 3,000 m. This shift toward higher elevations is also reflected in the northward movement of high post-seismic landslide area density (Figures 4-6). In all three birth categories, landslides occurred much less often at lower (<1,500 m) and higher (>5,000 m) elevation ranges. For context, 46% of the mapping area has an elevation <1,500 m and 8% of the area has an elevation >5,000 m.
Elevation is also correlated with the timing of grid cell recovery (Figure 7b). Landslides within cells that recovered prior to the earthquake tended to preferentially occur within elevations between 1,500 and 4,500 m. This is lower than the range for pre-earthquake birth cells, suggesting that lower-elevation locations were more likely to recover. Locations which recovered after the earthquake were similarly shifted toward lower elevations, likely reflecting the elevation range at which most coseismic landslides occurred. The coseismic landslide death distribution (i.e., landslides that disappeared between the coseismic 2015 and post-monsoon 2015 epochs) is peaked at higher elevations of 3,000-4,500 m, albeit with a small sample size (n = 117 cells) compared with the other classes. Cells in which landslide activity persisted through the end of the study period typically occurred at higher elevations (mostly >2,500 m), again indicating that hillslopes at higher elevations are recovering more slowly.

Slope Angle
The influence of slope angle on landslide birth and death also differs between categories. Landslides within grid cells with a pre-earthquake landslide birth preferentially occurred on slopes of between 30 and 55°, with a modal peak of 40° (Figure 7c). The coseismic birth distribution was shifted toward noticeably higher slope angles, occurring more commonly than would otherwise be expected between 35 and 65° and with a pronounced peak at 45°. Landslides in cells activated after the earthquake showed a partial return to pre-seismic conditions. Locations which recovered prior to the earthquake were typically within the same range of slope values as those where pre-seismic landslides had occurred, 30-55° (Figure 7d). Locations which recovered after the earthquake broadly reflect the slope distribution of coseismic births, but with a shift toward slightly lower slopes. In contrast, the distribution for cells in which landsliding persisted beyond post-monsoon 2018 was offset toward higher slope angles and with a peak at 42°, indicating that landslides on steeper hillslopes were generally more likely to persist whereas those on shallower slopes were likely to recover more quickly. Cells which recovered during the coseismic inventory show a very different distribution, with a peak at much lower slope angles (20°), but this curve is again defined by a relatively small sample size.

Slope Direction to Epicenter
Slope direction relative to the Gorkha earthquake epicenter is defined so that 0° represents hillslopes that directly face the epicenter, 90° represents those perpendicular, and 180° represents those that face directly away from the epicenter. Pronounced asymmetry in the density of landslides due to slope direction rela-KINCEY ET AL.

10.1029/2020JF005803
16 of 27 Figure 7. Topographic controls on landslide birth and death, displayed as difference in kernel density values relative to the density distribution from the overall study area. Results are shown for elevation (a, b), slope angle (c, d), slope direction relative to the Gorkha earthquake epicenter (e, f), and normalized distance from river to ridge (g, h), with the left column showing landslide birth categories and the right column showing landslide death. Difference data are based on the individual landslide polygons, with the density distribution for the overall study area shaded in light gray. Blue, black, and red lines show results for landslides in cells that were activated or recovered before, during, or after the earthquake, respectively. Green lines in the right column show results for landslides which persisted through the end of the study period, irrespective of birth epoch. The half-width of the kernel has been delimited in dark gray for each graph. Density values within these areas are likely to be less accurate due to edge effects during the density estimation. Further information on how to interpret this figure is provided in Section 6 of the supporting information.
tive to the seismic source has been previously demonstrated for coseismic landslide inventories (Robinson et al., 2017), due to localized variability in topographic amplification of seismic waves (Meunier et al., 2008). The distribution of values for grid cells that were activated during and after the earthquake shows that landslides occurred preferentially on hillslopes facing away from the epicenter (90-175°) (Figure 7e). The pre-earthquake distribution also shows, however, that landslides occurred preferentially at these slope direction values. Since this distribution is entirely independent of the direction of seismic shaking, the pattern is likely to instead be related to the direction of prevailing monsoonal rains in the region, which typically originate from the south-southeast and which may therefore preferentially impact upon south-facing hillslopes with relative slope directions of 90-180°.
The pre-earthquake and coseismic landslide death distributions show that hillslopes approximately perpendicular to the epicenter (75-150°) were more likely to recover than those facing more directly toward or away from the epicenter (Figure 7f). Again, because the pre-earthquake distribution is independent of the earthquake, this likely reflects the aspects on which most rainfall-triggered landslides occurred, and may also indicate that hillslopes directly facing the direction of monsoonal rains were less likely to recover than those which were partially sheltered. In contrast, grid cells which recovered after the earthquake as well as those that persisted beyond post-monsoon 2018 broadly reflect the equivalent birth distributions. This again suggests that, for this particular earthquake and study area, the slopes most affected by the earthquake are also those which are most susceptible to rainfall-triggered landsliding due to the predominant direction of monsoonal rains.

Normalized Distance from River to Ridgeline
Landslides are often observed to occur more frequently toward the top and bottom of a normalized hillslope profile. Seismic shaking is amplified by topography, triggering a concentration in landsliding near ridge tops (Meunier et al., 2008). In contrast, locations close to a river can also result in larger landslide frequencies due to high pore fluid pressure and undercutting of hillslopes (Densmore & Hovius, 2000). Landslides in grid cells activated prior to the Gorkha earthquake showed a clear tendency toward lower hillslope positions close to rivers and were underrepresented toward ridge tops (Figure 7g). In contrast, coseismic landslides occurred more uniformly along the entire hillslope profile, with the exception of the extreme ridge top locations. Landslides in cells activated after the earthquake show a similar distribution to the coseismic inventory, although shifted slightly toward the pre-seismic pattern.
Although pre-earthquake landslides were more likely to occur at lower hillslope positions, landslides within cells which recovered prior to the earthquake were typically at mid-slope locations (Figure 7h). Landslides within cells which recovered after the earthquake broadly followed the coseismic and post-seismic birth distributions. Landslides in cells which had not recovered by post-monsoon 2018, however, show a clear tendency toward lower hillslope positions close to the river network. This persistent activity may be due to fluvial undercutting of hillslopes (Cook et al., 2018) or downslope movement of landslide material by debris flows (Dahlquist & West, 2019;Tian et al., 2020).

Discussion
Our results provide the first comprehensive, regional-scale demonstration of how the pattern of landsliding has varied in the years after a large continental earthquake. While our findings agree in some ways with the understanding gained from studies of previous earthquakes, there are also some important differences-most notably, our observation that the number and area of potentially active landslides remained high through the first 3.5 years after the earthquake, and were nowhere close to recovering to pre-earthquake levels. Our study also documents systematic and profound shifts in the pattern of landslide activity over time, both during and after the earthquake, and at both the regional scale and the scale of individual hillslopes. Below, we discuss the spatio-temporal evolution of landsliding in more detail and speculate on some potential controls behind the patterns that we observe, before considering the implications of our results for post-earthquake hazard, response, and reconstruction efforts.

Conceptual Approach
Determining the rate of landscape recovery after a large earthquake depends fundamentally on (1) the definition of recovery and (2) the data that are used to assess it. Most previous multi-temporal inventory studies have tended to focus on either the occurrence of new landslides (e.g., Marc et al., , 2019 or new landslides plus visible reactivation of coseismic landslides (e.g., Fan et al., 2018;C. Li et al., 2018) to define the rate of post-seismic landsliding; the temporal decay of this rate to pre-earthquake levels is then used as a proxy for recovery to assumed pre-earthquake conditions of hillslope stability (e.g., . While this approach accurately characterizes the susceptibility of the landscape to new failures, it will underestimate the hazard posed by coseismic and post-seismic landslides as long as those failures are persistent within the landscape-that is, as long as they continue to move and supply sediment to downslope areas. The widespread occurrence of debris flows after large earthquakes (e.g., Dahlquist & West, 2019; illustrates the ease with which landslide debris may be remobilized from pre-existing failures. Some studies have attempted to address this shortcoming with qualitative estimates of reactivation or remobilization (e.g., Fan et al., 2018;C. Li et al., 2018). In the absence of intensive field monitoring, however, this determination is challenging, even with high-resolution imagery; furthermore, neither field monitoring nor the acquisition and analysis of high-resolution imagery are yet feasible over the scale of the entire rupture area in a large continental earthquake. An exclusive focus on new landsliding also implicitly assumes that revegetation happens rapidly enough to stabilize existing failures and impede the further transport of debris downslope. A growing number of studies have now demonstrated that the link between revegetation and reduced failure rates after the Wenchuan earthquake is complex (e.g., Shen et al., 2020;Yang et al., 2017Yang et al., , 2018, with both spatial and temporal variations in vegetation regrowth and landslide occurrence (Ni et al., 2019). This complexity is to be expected from the diverse mechanical and hydrological ways in which vegetation affects shallow hillslope stability, especially in a landscape recently impacted by an earthquake (e.g., Cohen & Schwarz, 2017;Stokes et al., 2014).
In contrast, by mapping each epoch independently, our approach implicitly assumes that, to first order, landslide activity is only effectively impeded once bare earth is no longer visible in the medium-resolution imagery. Until that point, we assume that the landslide can potentially undergo continued failure or remobilization of landslide debris and can continue to supply sediment to downslope areas. This is clearly also an oversimplification of recovery-not least because the point at which bare earth is no longer visible will depend upon the specific sensor, band combination, and spatial resolution that are used. By employing Landsat 8 and Sentinel-2 data throughout our study, we have sought to minimize these variations across the different epochs. The trade-off with this approach is a decreased ability to resolve both the smallest failures ( Figures S4 and S6) and the extent of continued activity or remobilization within existing failures (e.g., Fan et al., 2018). Our approach therefore provides a precautionary view of the potential hazard posed by existing landslides, due to the overestimation of the time period over which they may remain active, but this may be a more defensible approach when assessing risk. In sum, both approaches have advantages and disadvantages, and it is critical to bear their differences in mind when considering the results of any multi-temporal analysis. These differences point to the need to better understand the detailed biophysical relationships between what is visible in imagery, and how vegetation recovery and mass wasting evolve in the aftermath of a large earthquake.

Temporal Variability in Landslide Rates
Our results demonstrate that both the number of landslides and the total landslide area remained well above coseismic levels to the end of 2018, approximately 3.5 years after the Gorkha earthquake ( Figure 2). This indicates that the recovery period when measured in terms of the number and extent of persistent and therefore potentially active landslides present within the landscape is considerably longer than initial estimates suggested (e.g., Marc et al., 2019), and is instead likely to be on the order of 10 years or more, based on extrapolation of the landslide area results beyond our current time series. We note that studies which mapped all visible failures in the years after the Wenchuan (e.g., Chen et al., 2020;R. Huang and Li, 2014;C. Li et al., 2018;Tang et al., 2016) and Kashmir (Shafique, 2020) earthquakes have also inferred decade-scale recovery times.
Our results do match the rapid recovery inferred by other studies, however, in terms of a return to pre-earthquake landslide and debris flow rates within 1-2 years (Dahlquist & West, 2019;Marc et al., 2019), when we focus only on the changing number of newly impacted 1 km 2 grid cells (Figure 6d) or the landslide birth rate based on individual 100 m 2 landslide cells (supporting information 7). Analyzing the multi-temporal inventory data in this way indicates that the number of new landslide cells remained high in post-monsoon 2015 but then declined rapidly, with the post-monsoon 2018 level being only marginally higher than that for the post-monsoon 2014 inventory ( Figure S11). The difference in per-epoch landslide cell birth and death counts is even more pronounced, with the net change indicating that more landslide cells were disappearing than being newly created by the pre-monsoon 2016 epoch. This suggests that, when measured simply in terms of the occurrence of new landslide activity, the recovery rate following the earthquake was comparable to that obtained from previous studies that have taken a similar methodological approach, with the landscape returning close to pre-earthquake levels within just a few years. Importantly, however, this approach overlooks the persistence of existing landslide scars within the landscape, as illustrated by the high number of cells which were still mapped as landslides in post-monsoon 2018.
While the spatial and temporal patterns of post-earthquake recovery revealed by our data are highly complex, there is evidence for distinct clustering and evolution of landslide activity after the earthquake (Figure 6). The complexity of the changes that we describe reflects the numerous landslide types and sizes triggered by the earthquake (Kargel et al., 2016), which are likely to each change through time in a different manner. It is therefore reasonable to assume that the recovery of each landslide reflects both failure mechanism and setting, and so exploring recovery as a function of landslide type or mechanism would be desirable. While the nature of coseismic landslides is to some extent controlled by conditions at the time of the earthquake (Roback et al., 2018), the degree to which sediment production from coseismic landslides is fully accomplished during initial shaking must also define the future potential for continued landsliding in the period that follows . It is widely understood that some coseismic landslide types are more sensitive to reactivation under rainfall (e.g., as debris flows), but equally that some evolve via mechanisms that have a more complex relationship with controlling conditions and hence time. This includes, for example, progressive landslides initiated by the earthquake (Terzaghi, 1950), which are likely to exhibit some independence from environmental forcing and so remain difficult to associate with characteristic timescales of recovery beyond those set by site-specific conditions.
Interestingly, there remains only limited forensic research on how landslide failure mechanism or the evolution of hillslope-scale rock strength influence the trajectory of post-seismic recovery (e.g., Chen et al., 2020), beyond the exhaustion of supply to debris flows (e.g., Dahlquist & West, 2019). It seems reasonable to expect some form of asymptotic reduction in the intensity of many post-seismic landsliding processes, due for example to exhaustion of supply. Temporal changes in other processes may be quite different in form, however, which may be reflected in the landscape-scale patterns that we report. For example, we note an incremental post-seismic shift to failures at higher elevation which may tentatively reflect a retrogressive mode of landslide change; this is only mechanically possible for some types of landslide in some positions in the landscape. We also note that some processes which are incorporated into the inventory mapping, such as landslide runout and scarp stabilization, may act to cancel each other out in aggregated statistics, and so are not captured in averaged measures of landslide activity. Therefore, without a single dominant coseismic landslide type, it remains to be shown how such mechanism-dependent timescales interleave, and at what scale these can be summarized into a single event 'recovery' timescale. New post-seismic landslides are one relatively well-constrained component of this wider range of processes, but again that approach commonly aggregates landslides driven by a number of failure mechanisms. It is clear that a fuller understanding of post-earthquake landsliding needs to unpick how these various landslide types and mechanisms evolve and recover.
We have not converted total landslide area to volume because of uncertainty relating to scaling relationships (e.g., G. Li et al., 2014) as well as potential amalgamation of landslide areas  that cannot be resolved with our medium-resolution imagery. Nonetheless, the persistence of large areas of potentially active landsliding across the study area through to the end of the post-monsoon 2018 period is consistent with, although not conclusive proof of, the continued presence of large amounts of coseismic and post-seismic landslide debris. We tentatively infer from this observation that much of this material is there-fore still sequestered within the landscape. Similarly, R. Huang and Li (2014) documented that 80%-90% of erodible material remained on hillslopes 5 years after the Wenchuan earthquake, while Fan et al. (2018) suggested that more than 90% of the coseismic landslide debris produced during the Wenchuan earthquake was still contained on hillslopes and in low-order channels after 7 years.
Somewhat surprisingly, our results show slightly higher total landslide areas for the pre-monsoon inventories in 2016 and 2017, compared to their post-monsoon equivalents ( Figure 2). All else being equal, we would expect the monsoon both to trigger new landslides and to reactivate existing landslides, leading to higher total post-monsoon areas. Close examination of the inventories suggests, however, that this counter-intuitive result may be due at least in part to the visibility of landslide outlines under post-monsoon vegetation conditions, relative to pre-monsoon conditions (Figure 8). Ground visibility is greater in our medium-resolution imagery in pre-monsoon imagery, whereas landslide boundaries in the post-monsoon mapping periods are partly obscured by overhanging vegetation. This leads to persistent decreases in the mapped extent of post-monsoon landslides compared to their pre-monsoon equivalents. In support of this observation, we note that the area distributions for the pre-monsoon 2016 and 2017 inventories are shifted toward larger areas relative to the post-monsoon inventories (Figure 2b), but that the number of mapped landslides does not show the same fluctuation ( Figure 1a). Comparable seasonal variability has been noted in multi-inventory landslide studies that analyze vegetation recovery (e.g., Ni et al., 2019;Yang et al., 2018), but our results indicate that these phenological variations may also be important for inventories created from manual mapping of optical imagery. Importantly, the observed cyclicity in pre-and post-monsoon landslide area totals does not appear to correlate with interannual variability in monsoon rainfall totals (Figure S8), and remains small in proportion to the wider longer-term changes that we observe. Our results do, however, highlight a potential uncertainty in single epoch inventories, or single inventories that are mapped from imagery that spans multiple vegetation seasons.

Spatial Distribution of Landslides
Our area density results indicate that, even at the resolution of our approach (1 km 2 ), there is clear and consistent evidence for both persistence of existing pre-seismic and coseismic landslides and new post-seismic landsliding. Overall, the most intense areas of landslide activity appear to have migrated northward, from a northwest-southeast belt of intense coseismic activity to lower-order parts of the regional-scale drainage network by the end of the study period ( Figure 5). We also observe an overall shift in activity up and away from trunk streams toward higher parts of the landscape (Figure 7a Somewhat paradoxically, this large-scale shift in the locus of landslide activity occurs at the same time that debris from existing pre-seismic and (especially) coseismic landslides is remobilized and transported downslope and into the drainage network ( Figure 1c). Changes in area-perimeter ratios through time show increasing elongation of landslide form between the earthquake and post-monsoon 2016, followed by a reduction in elongation through 2018 (Figure 2b). By the end of the study period the area-perimeter ratio was still well below the average levels from the three pre-earthquake inventories, likely indicating widespread progressive runout and remobilization of coseismic failures in the years after the earthquake. This remobilization has led to persistent disruption to infrastructure, particularly the highway network, as documented by Tian et al. (2020). It is important to distinguish these distinct patterns of migration from each other, as they are occurring at different scales, and again may otherwise cancel out in wider scale averaged measures of change. We note as well that our ability to assess remobilization and downslope transport is almost certainly limited by the constrained widths of landslide and debris-flow deposits within the confines of the channel network, especially given the resolution of our imagery and concealment by vegetation at lower elevations. Thus, our estimates of downslope runout and landslide elongation must be taken as minima.
It would be reasonable to expect that the spatial footprint of landsliding triggered by the 2015 earthquake might evolve with time to one more characteristic of monsoon-triggered landsliding. By the end of the study, the area of the landscape that exceeded pre-earthquake landslide densities remained largely coincident with the full coseismic landslide footprint (Figures 4a and 5a), rather than being concentrated around the areas with the highest coseismic landslide densities. The implication is that even after 3.5 years, a wide area remained subject to continued landslide hazard, and thus that at the event scale coseismic landslide density alone is a poor predictor of the potential for persistent landsliding in the post-earthquake period.
In other words, persistent landsliding should be expected anywhere within the affected area, not just in the areas with highest coseismic landslide density-a key point for post-earthquake recovery that we explore below. Within the overall affected area, the spatial changes in landslide distribution through time are not random, and show clear patterns in terms of which areas have recovered since the earthquake and which have not ( Figure 6). Clustering was greater for areas with landslides which persisted until the end of the study period than for those which recovered prior to that time, with the latter tending to be more spatially isolated and located around the periphery of clusters containing persistent landslides (Figures 6a and 6b).
An important outcome of our study is the recognition of substantial new post-seismic landsliding, occurring in areas that did not fail during the earthquake. Almost 10% of cells across the study area show post-seismic onset of landsliding that persisted through the end of the 2018 monsoon ( Figure 6c). The clustering of these areas indicates that at this scale, the location of coseismic landsliding is not necessarily a good predictor of the location of post-seismic landsliding, which raises important considerations for hazard assessment. We note specifically that, had we chosen to focus on a smaller study area, recognition of this pattern would have been considerably more difficult. In future large continental earthquakes, it is therefore important to establish the overall spatio-temporal distribution of landsliding, rather than making wider inferences based on small subsets of the affected area.

Controls on the Evolution of Co-and Post-seismic Landsliding
The systematic spatial evolution of co-and post-seismic landsliding after the 2015 Gorkha earthquake suggests that there are some underlying physical controls on the locations that are most susceptible to reactivation and continued activity, as well as those that have recovered. Understanding these controls is critical to our ability to anticipate how landslide hazard is likely to change after a future event. Here, we compare our results to inferences from previous large continental earthquakes.
Stark spatial variations in the speed of recovery have been noted following the 2005 Kashmir and 2008 Wenchuan earthquakes (e.g., Khan et al., 2013;Ni et al., 2019;Shafique, 2020;Shen et al., 2020), but to our knowledge our study is the first to demonstrate this effect over (nearly) the entire earthquake-affected area.
A number of studies have tried to link these variations to elevation, slope angle, or aspect. For example, Khan et al. (2013) argued that landsliding rates stabilized within the first few years after the Kashmir earthquake, and that sites at lower elevations recovered faster in terms of the regrowth of vegetation and inferred cessation of activity. We observe qualitatively similar shifts after the Gorkha earthquake, although the lack of spatio-temporal information in Khan et al. (2013) means that we cannot directly compare the spatial pat-terns observed in Nepal (Figures 3 and 4) with those in Kashmir. Similarly, both Fan et al. (2018) and Yunus et al. (2020) observed that post-seismic landslides after the 2008 Wenchuan earthquake were more likely to occur at higher elevations, a trend that Yunus et al. (2020) linked to slower rates of vegetation recovery in those locations. In contrast, C. Li et al. (2018) argued that new landslides in part of the Wenchuan earthquake-affected area occurred preferentially at lower elevations except during periods of high rainfall; they tentatively linked this preference to river erosion and human activities, and distinguished between new landsliding and the progressive shift of coseismic landslide debris toward lower parts of the landscape. At present, we do not know whether the progressive shift in landslide activity toward the headwater reaches of the drainage network after Gorkha, as visible in Figure 5b, reflects large-scale vegetation change, or is due to another mechanism, such as the subsequent failure of hillslopes that were damaged in the earthquake but did not fully collapse (cf. Dadson et al., 2004), large-scale retrogression of failure upslope from existing landslide locations, or some other cause. Investigation of this spatial shift is a priority for future research.
Slope angle and aspect have also been invoked in disparate ways to explain landslide patterns after past earthquakes. Fan et al. (2018) argued that post-seismic landsliding typically occurred on lower slope angles, reflecting remobilization of landslide debris within the channel network-a pattern also noted by Hovius et al. (2011), who showed that the locus of landslide activity migrated to lower hillslope positions over time after the 1999 Chi-Chi earthquake. In contrast, both C. Li et al. (2018) and Yunus et al. (2020) documented a shift in new landslides toward progressively higher slope angles, a trend that Yunus et al. (2020) linked to slower revegetation on those steeper hillslopes. Our results provide evidence of a return in post-seismic landsliding toward similar slope angle values to those that were active before the earthquake (Figure 7c), as well as continued persistence of landslide activity on the steepest hillslopes in the landscape (Figure 7d). In terms of aspect, Fan et al. (2018) showed that post-seismic landslides typically occurred on the same southeast-facing hillslopes as the coseismic inventory. They suggested that this could be due to both climatic factors as well as the distribution of earthquake damage. This similarity is broadly reflected in our results, which demonstrate coseismic and post-seismic landsliding on hillslopes with common orientations (Figure 7e). In their study of landslide patterns in the 10 years after the Wenchuan earthquake, Yunus et al. (2020) showed that west-and northwest-facing hillslopes, in the shadow of the regional monsoon wind direction, had fewer landslides but also recovered more slowly than hillslopes at other aspects. This suggests that post-earthquake hillslope recovery will vary over multiple spatial and temporal scales depending on, among other factors, the distribution of hillslope damage, the local aspect, and the prevailing climatic conditions of the area.
A number of prior studies have also linked post-seismic landslide persistence to human activities, including road-building, construction, terracing and other agricultural activities, and mining (e.g., Khan et al., 2013;C. Li et al., 2018). Anecdotal observations from the Gorkha earthquake show that comparable links are likely to have been present in Nepal as well. This is an important phenomenon, not least because reconstruction activities are likely to be focused in the first few years after a large earthquake when landslide rates also remain high. Similarly, it is highly likely that (re)construction on ground damaged by the earthquake, or in areas previously considered safe for occupation or infrastructure, will subsequently be far more challenging. Some of the post-seismic increases in landslide density within the Lesser Himalaya and south of Kathmandu ( Figure 5b) correspond directly to major road corridors, and are almost certainly indicative of failures along either these corridors or along lesser branch roads, visible even at the regional scale.
Landslides associated with the construction of new rural roads are visible in our medium-resolution imagery (Figure 9), particularly in the pre-and post-monsoon 2018 inventories. This uptick tallies with the proliferation of rural road construction associated with the set-up of the new federal structure of government in Nepal in 2017 (Sudmeier-Rieux et al., 2019). We are not aware of any systematic attempt to document the growth of the rural road network in Nepal or its links to persistent post-earthquake landslide activity, and this remains another priority for future research. It may be the case that post-seismic landsliding will not return to pre-earthquake levels in locations undergoing widescale infrastructure development due to the consequent impacts on the landscape.

Implications for Recovery and Reconstruction Following a Large Earthquake
It is common that the scientific effort focused on understanding catastrophic earthquakes is poorly aligned with decision making and information needs on the ground (Datta et al., 2018;Williams et al., 2018). A key challenge relates to providing timely information on the likely future dynamics of landsliding triggered by earthquakes in a timeframe that dovetails with response and reconstruction. In this context, the rate and pattern of landslide activity and landscape recovery after a large earthquake are a first-order control on long-term seismic hazard, and so both synoptic and site-specific knowledge of what is likely to occur are highly valuable for both recovery and reconstruction, and future preparedness planning. We have documented the persistence of large areas of potentially active landsliding across the entire Gorkha earth-KINCEY ET AL.
10.1029/2020JF005803 23 of 27 quake-affected area, even after 3.5 years, indicating that hillslopes have not recovered as rapidly as might be inferred by studies that focus only on new failures (e.g., Marc et al., 2019). Our data also illustrate that, importantly, the nature of landslide hazard and therefore risk can significantly change through time after an earthquake. Our mapping illustrates that landslides triggered on the day of the earthquake, or even in the first few years after the earthquake, can evolve in terms of location. This dynamic risk must be reflected in messaging for earthquake-affected populations, and recovery and reconstruction efforts need to be cognizant of the potential continued changes that the landscape might experience. Our approach is in part intended to inform this type of effort in future and to feed into risk sensitive planning that accounts for rainfall-and earthquake-triggered landslide risk (e.g., Milledge et al., 2019), as well as the processes that play out in the aftermath of a large earthquake. Our approach provides a more conservative view of the persistence of landslide activity than earlier studies, which is arguably more appropriate in guiding risk averse planning that must also recognize people's everyday livelihood concerns and wider systemic risks (Oven et al., 2021). While there is an understandable desire to rebuild, doing so too quickly and ignoring the risk of post-seismic landsliding and remobilized landslide debris can have devastating impacts on people, their homes, and infrastructure.
Our results demonstrate three key messages for dynamic landslide hazard. First, for the study area as a whole, coseismic landslide occurrence-but not landslide density-is a reasonable first-order guide to the locations of persistent post-seismic landsliding over the 3.5 years of our study period (Figures 6c and 6d). Second, the coseismic landslide pattern cannot capture the nearly 10% of cells that show new and persistent post-earthquake landsliding-so the coseismic pattern should not simply be projected forward in time to anticipate the evolution of the hazard. Finally, the spatial clustering observed in our data means that areas of post-earthquake landslide activity should be expected to persist, all else being equal. These findings are useful for triaging areas of greatest risk in the immediate aftermath of a large earthquake, and for highlighting the manner in which landslide hazard is likely to evolve as disaster response and reconstruction progresses. More locally, one important caution illustrated by the dynamic nature of the hazard is that cultural memory or experience of living with landslide risk in these regions developed prior to a large earthquake may be less relevant in dealing with coseismic and post-seismic landsliding in the years after an earthquake (Alexander, 2000).
We see several promising ways in which our multi-temporal inventories could be used in order to provide a better assessment of future landslide hazard across the Gorkha earthquake-affected area in Nepal. The first relates to providing a regional-scale understanding of the continued evolution of landslide activity, as coseismic and post-seismic landslide debris is remobilized and as existing landslide complexes evolve, with a view to informing ongoing recovery and reconstruction ( Figure 1c). Further analysis would require physically based modeling of the release, entrainment and transport of landslide debris (e.g., Croissant et al., 2019), perhaps embedded within a wider consideration of the earthquake-triggered hazard chain (e.g., Fan, Scaringi, Domènech, et al., 2019;Fan, Scaringi, Korup, et al., 2019). Second, our inventories could be used to tune regional-scale post-earthquake landslide susceptibility models, in order to improve their skill in anticipating the occurrence of new failures relative to static models that do not account for the Gorkha earthquake. Finally, our data provide an assessment of the trajectory of stability for individual landslides, which adds considerably to static one-off site assessments, providing insight for ongoing mitigation of landslide risk in this region of Nepal. This is relevant to the National Reconstruction Authority's on-going National Geohazards Assessment, which seeks to support landslide-affected communities following the 2015 earthquake.

Conclusions
We have systematically mapped pre-, co-, and post-seismic landslides across the 14 administrative districts of Nepal that were most intensively affected by the 2015 Gorkha earthquake, covering about 90% of the total area affected by coseismic landsliding and continuing until the end of the 2018 monsoon. This represents, to our knowledge, the first systematic multi-temporal landslide inventory that spans nearly the entire rupture area of a large continental earthquake across both the earthquake itself and the following four monsoon seasons. Our mapping documents a sharp increase in both number and area of landslides caused by the earthquake, and further increases after the first post-earthquake monsoon in 2015. The numbers and total area of mapped landslides have fluctuated in the following years, with greater areas in the pre-monsoon inventories-likely due to seasonal changes in vegetation and visibility of landslide margins. Nevertheless, both landslide number and area remained greater than the coseismic inventory at the end of this study period, indicating a high degree of persistent and evolving landslide activity. It is therefore important to recognize that landslide hazard remains significant several years after this earthquake.
Coseismic landsliding was concentrated in a northwest-southeast band along the physiographic divide between the Lesser and Higher Himalaya. In general, coseismic landslides were larger on average than those that occurred before the earthquake, and they occurred at lower elevations, on steeper slopes, and higher up on individual hillslope profiles. The years immediately after the earthquake have seen the continued evolution of landsliding across the whole impacted area, rather than a progressive constriction around the most severely impacted locations. We see clear evidence that post-seismic landslide activity has shifted northwards, higher in the drainage network, although there is also some evidence of recovery to pre-earthquake patterns in terms of slope gradients and hillslope positions. This northward shift is superimposed upon a smaller-scale change in the pattern of landslide activity, as pre-and coseismic landslide material has been remobilized and transported downslope into the drainage network. We also see substantial evidence of post-seismic landsliding in areas that were not badly affected by coseismic impacts, both in the far northern margins of the rupture area and within parts of the Lesser Himalaya.
When amalgamated into landslide area density estimates on a 1 km 2 spatial grid, our results indicate that most areas of high landslide activity have persisted throughout the study period, irrespective of whether activity started before or during the earthquake. These areas of persistent activity are highly spatially correlated and are somewhat predictable on the basis of simple topographic metrics, which is important for understanding the temporal evolution of landslide hazard. We see evidence of local systematic shifts in landsliding, with notable reductions in landsliding in the lower slopes of major river valleys over time, and a tendency for spatial clustering of specific landslide birth and death trajectories. Further studies of post-earthquake landslide hazard could focus on the continued remobilization and transport of existing landslide debris, or on the ways in which our inventories can be used to constrain changes in landslide susceptibility.
Our study takes a distinct conceptual approach in mapping all visible areas of disturbance during each epoch, rather than concentrating solely on the occurrence of new landslides. We demonstrate that unpacking landsliding after an earthquake in this manner reveals a wide variety of mass wasting processes that may each follow a different trajectory to recovery, so that recovery warrants both careful definition and choice of how it is described. This methodology assesses directly the time period over which an existing landslide may remain active, and thus yields a more precautionary approach to generating information about hazards that is suitable for understanding the regional-scale development of landsliding.