Capturing the footprints of ground motion in the spatial distribution of rainfall-induced landslides

The coupled effect of earthquakes and rainfall is rarely investigated in landslide susceptibility assessments although it could be crucial to predict landslide occurrences. This is even more critical in the context of early warning systems and especially in cases of extreme precipitation regimes in post-seismic conditions, where the rock masses are already damaged due to the ground shaking. Here, we investigate this concept by accounting for the legacy of seismic ground shaking in rainfall-induced landslide (RFIL) scenarios. We do this to identify whether ground shaking plays a role in the susceptibility to post-seismic rainfall-induced landslides and to identify whether this legacy effect persists through time. With this motivation, we use binary logistic regression and examine time series of landslides associated with four earthquakes occurred in Indonesia: 2012 Sulawesi (Mw = 6.3), 2016 Reuleut (Mw = 6.5), 2017 Kasiguncu (Mw = 6.6) and 2018 Palu (Mw = 7.5) earthquakes. The dataset includes one co-seismic and three post-seismic landslide inventories for each earthquake. We use the peak ground acceleration map of the last strongest earthquake in each case as a predisposing factor of landslides representing the effect of ground shaking. We observe that, at least for the study areas under consideration and in a probabilistic context, the earthquake legacy contributes to increase the post-seismic RFIL susceptibility. This positive contribution decays through time. Specifically, we observe that ground motion is a significant predisposing factor controlling the spatial distribution of RFIL in the post-seismic period 110 days after an earthquake. We also show that this effect dissipates within 3 years at most.


Introduction
The conditions promoting the occurrence of a landslide are governed by various climatic, morphologic, geomorphologic, geotechnical, seismic and anthropic factors and their complex interactions (Budimir et al. 2015;Reichenbach et al. 2018). These causative factors are categorized as predisposing conditions and triggering factors (e.g. IAEG 2001; Tanyaş et al. 2019a;Fan et al. 2019). Predisposing conditions (e.g. weathering, morphology) typically refer to slowly changing processes which tend to keep the slope in a marginally stable state (IAEG 2001). Landslide triggering factors (i.e. rainfall, earthquake, human activity, snowmelt, volcanic processes), on the other hand, generally refer to external stresses that cause an immediate response in terms of slope stability (Crosta et al. 2012). Therefore, triggering factors might be most crucial in terms of rapid assessment of landslide hazard because they mostly dictate the timing of landsliding.
To capture the coupled effects of earthquakes and rainfall, first we need to refer to the concepts of triggering and predisposing factors. There are two possible scenarios ( Fig.  1): scenario A where an earthquake triggers landslides after a rainfall event(s), and scenario B where rock masses that are already disturbed by ground shaking subsequent rainfall event triggers landslides.
In scenario A, the preconditioning effect of rainfall followed by (or during) an earthquake is dependent on rainfall intensity and duration, which may cause an increase in pore pressure and a reduction in the ratio of resisting forces to driving forces (factor of safety, FoS) Faris and Wang 2014;Zhang et al. 2019;Martino et al. 2020).
In scenario B, the preconditioning effect of seismic shaking followed by a rainfall event could be related to multiple earthquakes, which may have previously and permanently modified the force balance in a given hillslope (Saba et al. 2010;Fan et al. 2018). Because of the multiple earthquakes responsible for the damage given to rock masses, the accumulated effect of seismic shaking in a given area may add another level of complexity that is difficult to account for.
There are several studies showing how earthquakes increase landslide susceptibility and play a major role in RFIL events in post-seismic periods (Lin et al. 2006;Saba et al. 2010;Tang et al. 2016;Yang et al. 2017;Fan et al. 2018;Tian et al. 2020;Kincey et al. 2020). Specifically, for the 2008 Wenchuan earthquake, the post-seismic landslide evolution has been examined in several articles (e.g. Chen et al. 2020a;Tang et al. 2016;Xiong et al. 2020;. For instance, Fan et al. (2018) assess postseismic slope failures for a subset of the area affected by landslides since 2008 and until 2015. They report an elevated landslide susceptibility after the Wenchuan earthquake, with the majority of landslides being observed as remobilized coseismic failures. They examine the elevated susceptibility by monitoring the variation in landslide rate, which is a term referring to the number of landslides mapped over a period of time. They show that the landslide rate gradually decreases and returns to its pre-earthquake level within approximately seven years. Fan et al. (2018) also note that the recovery of the landslide rate to pre-earthquake periods may be longer in areas affected by stronger earthquakes. Chen et al. 2020a) examine another subset of the area affected by co-seismic landslides for the Wenchuan earthquake, in the period between 2008 and 2018. And, they report that hillslopes and landslide deposits largely stabilize 10 years after the earthquake. Besides, several studies suggest a possible link between these long periods where we observe an elevated landslide susceptibility and large co-seismic landslide deposits (e.g. Chen et al. 2020a;Fan et al. 2018;Xiong et al. 2020;Yunus et al. 2020) estimated to sum up to 2.8 km 3 specifically for the Wenchuan case (Li et al. 2014). As regards the factors controlling the length of the period characterized by the elevated landslide susceptibility, Fan et al. (2018) refer to possible contributions of various processes such as progressive strengthening of the deposits due to grain coarsening (Hu et al. 2018) and, in particular, revegetation (Chen et al. 2020b;Yunus et al. 2020). In fact, studies examining the same concept associated with vegetation recovery argue that landslide recovery has not been fully completed yet (Xiong et al. 2020;Yunus et al. 2020), and it may take up to 25 years in total (Chen et al. 2020b).
The overview above summarizes the recent literature on post-seismic landslide recovery time, taking Wenchuan simply as a reference. Notably, differences still exist in the literature in terms of recovery times in other geographic areas and even on the exact definition of recovery time itself (Kincey et al. 2020). In fact, if on the one hand, the concept of recovery time has found a large support across the whole geoscientific community. On the other hand, scientific debates are still taking place on the required time window for the recovery to take place. Taking aside the differences in opinion among research groups, one key element has certainly been agreed upon, and this corresponds to the physics behind this process. Specifically, Ambraseys and Srbulov (1995), already a few decades ago, summarized the stages of landslide genesis in co-seismic and two post-seismic phases. Co-seismically, populations of landslides trigger across a given landscape as a function of ground motion and its duration, together with the geometry of the slope, and the undrained strength of the material mobilised during the earthquake. In the first post-seismic phase, landslides can trigger in response to rainfall if the Fig. 1 Schematic diagram showing two possible scenarios ending up with earthquakeinduced landslides (EQIL, scenario A) and rainfall-induced landslides (RFIL, scenario B). t 0 and t 1 indicate the chronologic order of events in two scenarios residual undrained shear strength (disturbed by the previous ground shaking) on the slip surface is less than that required to maintain static equilibrium. The second post-seismic stage is governed by creep and consolidation processes, together with destabilising hydrostatic forces. These govern the landslide pattern of occurrence whenever the given landscape is exposed to high or extreme precipitation. These post-seismic phases clearly emerge in a recent study published by Tian et al. (2020), where the authors summarize the current literature of post-seismic landslide evolution and determine two primary drivers: the amount of co-seismic source material and the precipitation pattern. They also note that stronger and more numerous earthquake aftershocks can prolong the recovery time.
Nevertheless, the abovementioned studies mostly focus on landslide rates observed at pre-, co-and post-seismic periods but not the predisposing effect of ground shaking in a multivariate scheme implemented to model post-seismic landslides. There is no proposed method to explicitly account for the legacy of previous earthquakes as a predisposing factor in RFIL susceptibility assessment. However, some proxies (e.g. distance to fault, distance to earthquake epicenter) are suggested to use in landslide susceptibility assessments to capture part of this effect (e.g. Gallen et al. 2015;Kritikos et al. 2015;Massey et al. 2018;Parker et al. 2015). For instance, Kirschbaum and Stanley (2018) propose a predictive model for RFIL where the landslide susceptibility map created including a variable regarding the seismicity. In particular, they use distance to fault lines among the conditioning factor. They then combine the landslide susceptibility map with the landslide triggering rainfall threshold. Furthermore, Quesada-Román et al. (2019) focus on capturing the predisposing effect of earthquakes in RFIL susceptibility assessments via numerical methods and use distance to the epicenter as a proxy. Specifically, they run a logistic regression model for landslides triggered by the 2016 Hurricane Otto, Costa Rica, 4 months after the 2016 Bijagua earthquake, which affected the same site. They show that a higher landslide density is observed close to the epicentral area where intense precipitation was also recorded. However, either distance to fault line or distance to epicenter has some limitations because these approaches do not consider the distribution of ground shaking, which can vary irrespective of distance from the fault line/epicenter.
A more comprehensive alternative to distance fault lineament/epicenter is represented by ground motion parameters (GMPs). GMPs (e.g. peak ground acceleration, PGA; Modified Mercalli Intensity, MMI) are commonly used in predictive model developed for EQIL (Nowicki et al. 2014;Nowicki Jessee et al. 2018;Tanyaş et al. 2019a).
This study aims to capture the role of ground motion as a predisposing factor in a landslide susceptibility assessment for RFIL in Indonesia. We map multi-temporal landslide inventories covering both co-and post-seismic phases associated with four earthquakes occurred in Indonesia: 2012 Sulawesi (M w = 6.3), 2016 Reuleut (M w = 6.5), 2017 Kasiguncu (M w = 6.6) and 2018 Palu (M w = 7.5) earthquakes. We hypothesize that the peak ground acceleration (PGA) map of the last strongest earthquake can partially explain the spatial distribution of landslides triggered by rainfall after an earthquake. To test this hypothesis in the validity domain of the areas under study, we opt to test the PGA contribution in RFIL susceptibility models built by using a binary logistic regression (BLR), which is the most common statistical model used in landslide susceptibility assessments (e.g. Reichenbach et al. 2018). Then, we run two separate tests to address the question mentioned above.
The first one involves fitting a BLR model to examine the relative contributions of both PGA map and two morphometric covariates (i.e. slope and distance to stream) in landslide susceptibility assessments carried out for co-and post-seismic multi-temporal landslide inventories. In doing so, we monitor how the contribution of PGA (i.e. regression coefficient of PGA) changes through time from coseismic to post-seismic periods.
The second test consists of fitting a BLR model for a case where ground motion and rainfall data are contextually available. In this case, the susceptibility model features morphometric properties, the PGA of the last strongest earthquake, together with rainfall proxies (i.e. daily accumulated and 7day antecedent precipitation) associated with RFIL. This operation ensures not only that we can estimate the PGA contribution to the estimated probability but also that we can compare it to the contribution of the rainfall proxies.
We should note that all our analyses are representative only for the areal boundaries encompassing multi-temporal inventories, which are mapped for a subset of the total area affected by co-seismic landslides. Notably, the characteristics of landslides may vary spatially, and therefore, the contribution of PGA that we assess in the susceptibility analyses is representative for the boundaries of examined areas.

Materials and study areas
We examine two seismically active sites (study areas A and B) from Indonesia (Fig. 2) and create multi-temporal landslide inventories via visual interpretation of PlanetScope (3-5 m), rapid eye (5 m) images acquired from Planet Labs (Planet Team 2017) and high-resolution Google Earth scenes.
To create the multi-temporal landslide inventories, we examine all available satellite images and choose the largest available cloud-free regions, for both sites. All the multi-temporal images we use for mapping convey the real landslide distribution over time during co-and post-seismic phases. Notably, the inventories are not created following a fixed temporal resolution. We map as many inventories as the imagery availability allowed. In each inventory, landslides that have previously occurred are eliminated, and only new failures are included. We systematically examine the satellite scenes through visual observation and map landslides as polygons. We delineate landslide source and depositional areas as a part of the same polygon.
The study cases we examine are different from most of the previously investigated situations in the geomorphological literature where the effect of major earthquakes (7.9 > M w > 7.0) (e.g. 1999Chi-Chi, 2005Kashmir, 2008Wenchuan, 2015 Gorkha earthquakes) is tested against a large sample of coseismic landslides (Lin et al. 2006;Saba et al. 2010;Fan et al. Fig. 2 This figure shows the areal extent of the study areas and the number of landslides/unstable slope units for the examined multi-temporal inventories respectively for a, b study area A and c, d study area B. The epicentral locations of earthquakes (M w > 5) occurred since 2010 are indicated by stars in panels a and c. Timelines in y axes of panels b and d are arbitrarily spaced. The elevation legend given in panel a and the legend indicating the date of the earthquakes given in panel c is common for both panels a and c 2018; Kincey et al. 2020). For those earthquakes mentioned above, USGS ShakeMap reports maximum PGA values ranging from 0.83 to 1.36 g (PGA max,chi-chi = 0.83 g, PGA max,Kashmir = 1.36 g, PGA max,Wenchuan = 1.14 g and PGA max,Gorkha = 0.87 g) (Worden and Wald 2016). Here, we examine two sites where ground motion originated from strong earthquakes (6.9 > M w > 6.0) with one exception, which is the 2018 Palu earthquake. Overall, the maximum PGA values are lower than other cases examined in the literature (PGA max,Sigli = 0.20 g, PGA max,Reuleut = 0.60 g, PGA max,Sulawesi = 0.32 g, PGA max,Kasiguncu = 0.25 g and PGA max,Palu = 0.85 g) (Worden and Wald 2016) (Table 1).

Study area A
Study area A was primarily hit by two earthquakes of magnitude greater than 6 (US Geological Survey 2017) since 2010: 21th January 2013 Sigli (M w = 6.1) and 6th December 2016 Reuleut (M w = 6.5) earthquakes (Fig. 2a). We scan an area of 1356 km 2 to create multi-temporal inventories. For the 2013 Sigli earthquake, we map only one post-seismic landslide inventory using satellite scenes acquired on 27th July 2016, about 3 years after the earthquake. This means that landslide data with temporal depth is not available for the Sigli earthquake. Therefore, we cannot include this post-seismic landslide inventory because it will not support the analyses through time. For the 2016 Reuleut earthquake, which occurred along a strike-slip fault, we map one co-seismic and three post-seismic landslide inventories (Fig. 2b). The maximum PGA estimated for this earthquake is 0.06 g, whereas PGA values range between 0.05 and 0.47 g in the study area (Fig. 3). We also observe that in all inventories, landslide sizes are relatively small and the average size of co-seismic landslides (6227m 2 ) is larger than the post-seismic ones (1181, 4849 and 4257 m 2 ) (Fig. 4a).
Intermediate, basic volcanic and mixed sedimentary rocks appear as the dominant lithologic units in which landslides are triggered (Sayre et al. 2014). The spatial distributions of landslides associated with each inventory are presented in Fig. 3, and details of the landslide inventories are also given in Table 1.
Study area A shows a rare situation where we have more post-seismic landslides than co-seismic ones. Specifically, the Reuleut earthquake triggered only 60 co-seismic landslides we interpreted as shallow translational slides, whereas the post-seismic inventory compiled 110 days after the earthquake contains 742 rainfall-triggered landslides. We consider it as a rare case because, in the literature, the peak landslide rate is mostly associated to co-seismic landslide event (Guzofski et al. 2007;Saba et al. 2010;Tang et al. 2016). In study area A, we also see that, in all post-seismic inventories, less than 1% of the landslide population is associated with reactivated co-seismic landslides, and the rest is characterized by new landslides.

Study area B
Study area B (Fig. 2c) was affected by three major earthquakes of magnitude greater than 6 (US Geological Survey 2017) since 2010: 18th August 2012 Sulawesi (Indonesia, M w = 6.3), 29th May 2017 Kasiguncu (Indonesia, M w = 6.6) and 28th September 2018 Palu (Indonesia, M w = 7.5). We use the co-seismic landslide inventories of these earthquakes, which were already mapped and examined by Lombardo and Tanyas (2020). We expand the dataset and also map post-seismic landslide inventories. To map the landslides, we examine an area of 1078 km 2 , where metamorphic and acid plutonic rocks are the dominant lithologic units (Sayre et al. 2014).
Since we map landslides over the same areal extent, we inevitably examine a different subset of the earthquake affected area for each earthquake in terms of level of ground shaking. For instance, in the Sulawesi earthquake (strike-slip), the examined area crosses the epicentral area, whereas in the Kasiguncu case (normal fault), it does not cover the area we observe the highest ground shaking (Fig. 5). Specifically, the maximum PGA estimate is 0.25 g in the Kasiguncu earthquake, whereas PGA values range between 0.03 and 0.16 g in the examined area. As for the Palu earthquake, which occurred along left-lateral strike-slip fault with north-northwestward strike (Socquet et al. 2019), the study area still covers a wide range of PGA values changing between 0.11 and 0.70 g ( Fig. 5 and Table 1).
It is important to note that the first post-seismic Sulawesi inventory is created approximately a year after the earthquake, and the inventory includes mostly co-seismic landslides, with a minor amount due to RFIL. Therefore, we consider it as a coseismic landslide inventory, although some RFIL noise cannot be excluded.
In total, we create 18 inventories of co-and post-seismic landslides (the latter being triggered by different rainfall events) associated with each earthquake in study area B (Fig. 2d). However, some post-seismic landslide inventories include only a few landslides. To increase the sample of landslides and avoid large uncertainties in the statistical analyses, we therefore aggregate some of the landslide inventories. In turn, we work with three post-seismic inventories for each earthquake (Table 1). The details of the aggregation are given in "Methods," and the spatial distribution of aggregated landslides is presented in Fig. 5.
Overall, the co-seismic landslide inventories include 520, 386 and 725 landslides triggered by the Sulawesi, Kasiguncu and Palu earthquakes, respectively. In each case, the majority of landslides are characterized as shallow slides. Also, in each case, the percentage of post-seismic landslide that appears to have interacted with previous failures is less than 5%. In other word, the majority of post-seismic landslides are new failures. Similar to the Reuleut case, also in these multi-temporal inventories, landslide size area relatively small and overall, co-seismic landslides are larger than their post-seismic counterparts (Fig. 4 b, c and d).

Methods
We structure our analyses in a three-stepped procedure summarized in Fig. 6. In step 1, we identify a suitable landscape partitions and pre-process morphometric, seismic variables and rainfall proxies to organize the dataset required to run a susceptibility model for each study area. In step 2, as part of susceptibility analyses, we examine ground shaking as a predisposing factor of landslides and investigate its relevance in each model trained by using the available inventories. As a result, we monitor how the PGA role in each model changes through time. In step 3, we focus on study area A and examine the coupled effect of PGA and rainfall proxies on the spatial distribution of a RFIL inventory. In particular, we use the first post-seismic inventory associated with the Reuleut earthquake (M w = 6.5) because it appears as a rare event where many RFIL occur on a site although only a few landslides are triggered by the earthquake. Our rationale is that the higher post-seismic landslide rate (compared to its co-seismic counterpart) may be due to the legacy effect of the Reuleut earthquake. If this is the case, then the PGA of the Reuleut earthquake should be able to explain part of the spatial dependence in a susceptibility model built with RFIL. Step 1 We first mask the flat regions, which are not prone to landsliding, to increase the accuracy of the landslide susceptibility model (Kritikos et al. 2015;Tanyaş et al. 2019b). We use the methods developed by Alvioli et al. (2018) to automatically define and remove flat areas. Specifically, we use GRASS GIS (Neteler and Mitasova 2013) r.geomorphon script of Jasiewicz and Stepinski (2013) to identify various landform classes. This algorithm calculates landforms and associated geometry using pattern recognition. The algorithm self-adapts to identify the most suitable spatial scale at each location and checks the visibility of the neighborhood to assign one of the terrestrial forms. Following the identification of pixels, we use the method developed by Alvioli et al. (2018), which gets rid of the sparse "flat" pixel. The algorithm starts from the pixels classified as "flat" by r.geomorphons and shrinks the borders of the flat raster map by a few pixels and then grows it again; the procedure is repeated until sparse pixels disappear. As a result of this procedure, we mask flat regions and exclude them for the rest of analyses. We divide the study areas into mapping units. This is a crucial step in any susceptibility assessment because the chosen mapping unit determines how dependent and independent variables are represented in space and are used to prepare the training and validation subsets for susceptibility modelling (Rossi and Reichenbach 2016). Also, the single mapping units are the geographic objects for which the probability of landslide occurrence is estimated.
Among the mapping units proposed, grid cells and slope units (SUs) are the most common terrain partitioning methods available in the literature (Reichenbach et al. 2018). We choose the SU, because they internally reflect similar hydrological and geomorphological conditions and are considered a well suited terrain subdivision for landslide susceptibility modelling (Carrara 1988;Guzzetti et al. 2006;Alvioli et al. 2016). SUs subdivide the terrain between streamlines and ridges under the constrain of slope and aspect within-unit homogeneity (Alvioli et al. 2016). For the automatic partitioning of a landscape into SUs, we use r.slopeunits, an open source software developed by Alvioli et al. (2016).
We use only a few independent variables to limit crossmodel differences due to variable interactions, allowing us to study in detail how the seismic shaking effect changes over time. Specifically, we use slope steepness and distance to stream as morphometric variables and peak ground acceleration (PGA) as a proxy for ground shaking. Slope steepness controls the ratio of resisting forces to driving forces and is notoriously related to the occurrence of landslides. Slope steepness is the most frequently used covariates of the entire landslide susceptibility literature (Reichenbach et al. 2018), and it also appeared statistically significant in 95% of all Black lines indicate co-seismic landslides inventories whereas red, blue and green lines refer to first, second and third postseismic landslide inventories, respectively landslide logistic regression studies (Budimir et al. 2015). Distance to stream is another important covariate particularly for rainfall-triggered landslides and used as a proxy reflecting hydrogeological stresses affecting hillslope stability (Budimir et al. 2015;Reichenbach et al. 2018). Overall, hydrogeological conditions are assumed to be less favourable towards the river channel due to the concentration of the groundwater flow and the destabilizing effect of river incision that contributes to slope instability (Reichenbach et al. 2018).
This covariate set will be used only to support the analyses in step 2 (Fig. 6), in which rainfall proxies will not be considered. Therefore, we will disregard the potential contribution of  (Worden and Wald 2016). Acquisition dates of pre-and post-scenes used to map landslides are indicated at the top of each panel rainfall patterns due to our lack of knowledge on which specific rainfall event(s) may be responsible for the slope failures. On the other hand, in step 3, we focus on a particular RFIL event(s) where we can estimate the most likely precipitation proxies triggering landslides (Fig. 6). We will elaborate this point in step 3 below.
To express the presence/absence of the landslide distribution over each study area at the SU level, we consider unstable conditions for SU covered by landslide polygons for a surface greater than 2%. This threshold value has been applied in several studies (Guzzetti et al. 2006;Galli et al. 2008;Rossi et al. 2010) to limit the mapping inaccuracy when digitizing the landslide inventory. This means that, if more than 2% of a SU is intersected by landslide(s), the SU is unstable, and in the susceptibility model, it refers to the presence condition.
Some of the post-seismic inventories contain only a few landslides (Fig. 2 b and d), which implies that some models may be affected by large uncertainties. For instance, we have seven post-seismic landslide inventories for the 2017 Kasiguncu earthquake, and the first two of them include 52 and 44 SUs (Fig. 2d) intersected by landslides (i.e. unstable SUs), whereas latter have two, three and one unstable SUs (Fig. 2d). To create the third post-seismic landslide inventory, we thus aggregate those inventories with only a few samples. We do so without interrupting the temporal order of inventories. The aggregated landslide inventories are presented in Fig. 5, where the acquisition dates of preand post-scenes used to map landslides are also indicated. Fig. 6 Workflow of the applied three-stepped procedure. Notably, each analytical step is associated with 1000 bootstrap simulations to retrieve the sampling distribution and confidence intervals of each parameter.
Step 2 indicates the exploratory analyses carried out in areas A and B. These are aimed at retrieving the temporal evolution of the covariates' effects. And, the temporal evolution of the performance as we measure the effect of the PGA.
Step 3 indicates the analyses carried out in area A, where we have run a consistent modelling protocol to the one presented in step 2, but this time simultaneously including ground motion and precipitation indices as covariates We summarize the characteristics of the covariates' distribution within each SU using mean and standard deviation for each of the selected independent variables Tanyaş et al. 2019a). Due to the coarse resolution of the precipitation proxies, we summarize the rainfall distribution per mapping unit uniquely by using the mean. Prior to any modelling step, we initially standardize each independent variable by mean zero-unit variance. This can be achieved by subtracting the mean value of each covariate (centring) and divide by the standard deviation. This operation ensures that the estimated regression coefficients will also be in the same unitless scale, making their effects on the final susceptibility model comparable (Camilo et al. 2017;Lombardo and Mai 2018).
Step 2 We use the binary logistic regression (BLR) method to assess the contribution of PGA in landslide susceptibility analysis conducted for each inventory. BLR is the most common approach used in the geoscientific literature to predict where landslides may occur (e.g. Budimir et al. 2015;Reichenbach et al. 2018). This statistical method is a multi-variate regression used when the target variable (Y) is expressed by two classes (i.e. presence/absence or stable/unstable slopes or 0/1). In landslide susceptibility studies, the aim is to model the conditional probability p(Y = 1| X n ), or in brief P(x), that Y is positive given a set of n covariates (X n ). As for the covariates, they can be both numerical and categorical in nature.
A BLR is denoted as follows: where β 0 is the global intercept and β n is the vector of regression coefficients associated with each covariate X n , and the results are probabilities confined between 0 and 1. The same equation can be re-written as follows: Where the left term is referred to as logit, later denoted as η or logit link function, and the solution is sought by estimating β 0 and β n . The probability for the ith slope unit in the study area can be calculated knowing the observed class, y i and the associated vector of covariates, x i , via the likelihood function. For the ith slope unit, the probability of the slope unit to be unstable or stable is either p, if y i = 1, or 1 − p, if y i = 0. The likelihood can be then written as: The logarithm of the likelihood or log-likelihood turns the above products into sums, as follows: which is then maximized by differentiating the loglikelihood with respect to the parameters, by setting the derivatives equal to zero. Alternatively, one can minimize the negative log-likelihood, which is identical to solve Eq. 4, but numerically easier to handle.
In the BLR scheme summarized above, we use two types of covariates: (1) a time-invariant morphometric set (slope μ , slope σ , Dist2Stream μ , Dist2Stream σ ) and (2) time-variant ground motion parameter (PGA μ , PGA σ ) to infer the functional relation with respect to the stable/unstable condition. And we develop only explanatory models to make inference in a way that can support expert-based interpretation of the process at hand. Therefore, we use the whole number of SU observations and associated covariate values and do not make crossvalidation. This is because we aim to understand instability processes with respect to ground motion and precipitation stresses but not to develop a predictive model. A predictive model is used to forecast future events, by calibrating it over a subset of the available information and validating it over the remaining cases.
However, by fitting the whole data, we obtain single parameter estimates, neglecting the model uncertainty. To retrieve the uncertainty associated with each estimate we present in this manuscript, we then implement an additional bootstrap step. Bootstrapping is a simulation-based technique, for which data is resampled with replacement (e.g. ) each time generating a new dataset generated from the distribution of the original one. This offers the chance to fit numerous times the given explanatory model, therefore retrieving the sampling distribution and confidence intervals for otherwise single parameters. In this work, we used the R (R-Team 2014) package boot (Canty 2002) to retrieve the distribution of each regression coefficient. To evaluate the overall modelling performance, we use receiver operating characteristic curves and their integrated area under the curve (ROC and AUC, respectively; Hosmer and Lemeshow (2000). For clarity, we remind here the reader that the ROC curves are constructed in a plane defined between the true positive rate (TPR) and false positive rate (FPR). TPR and it can be calculated from any confusion matrix as the ratio between true positives (TP) over the sum of TP and false negatives (FN) (see, Rahmati et al. 2019). FPR can be calculated from any confusion matrix as the ratio between false positives (FP) over the sum of FP and true negatives (TN) (see, Rahmati et al. 2019).
Overall, we implement 1000 bootstrap replicates. From each model we built, we store the information related to the regression coefficients. The coefficients are unitless and therefore are comparable among the covariates because they are in the same scale and have the same covariate sets across the models considered in step 2 (Fig. 6).
To evaluate the modelling performance, we use the receiver operating characteristics (ROC) curve. This curve is built by computing confusion matrices for specific cutoff values used to binarize the estimated probability. In other words, for a specific probability cutoff, one can calculate the proportion of true positives and true negatives. Then by using a large number of cutoffs, one can plot the pair of coordinates in a 2D plane described by 1-specificity (or FP/FP + TN) and sensitivity (or TP/TP + FN) (TP true positive, TN true negative, FP false positive and FN false negative; details can be found in Fawcett 2006). Because building a ROC curve requires the use of a large number of cutoff, the ROC is often referred to as one of the most efficient and cutoff-independent metric for statistical classifiers (Hosmer and Lemeshow 2000). The integral of the ROC curve or the AUC is then used to rank classification performance (see, Hosmer and Lemeshow 2000).
We also use run alternative models by excluding the PGA signal from each model we present. As a result, we aim at assessing whether the inclusion of the ground motion contributes to explain RFIL landslides (also known as jackknife test, e.g. Lombardo et al. 2016). And if so, how its inclusion contributes through time. Specifically, for each jackknife test, we calculate the difference in AUCs for each model, built with and without the PGA layers. We recall here that all the analyses contain a bootstrap step. To ensure an equal comparison between models with and without PGA, the 1000 resampled datasets are consistent in both cases.
Step 3 We extend the analyses by fitting a BLR model specifically for the post-Reuleut earthquake landslide inventory. In this case, we consider rainfall-related proxies in addition to morphometric and seismic factors (Fig. 6), aiming at testing if the signature of the ground motion effect can still be retrieved from RFILs. Since the first post-Reuleut RFIL inventory was mapped 110 days after the earthquake, we cannot point to a specific rainfall trigger for this RFIL inventory. Therefore, we use two rainfall proxies to represent potential triggering scenarios: daily accumulated and 7-day antecedent precipitation from IMERG data. Out of the large number of rainfall pattern realizations within 110 days, we opted to isolate the most likely candidate to have triggered landslides in two analytical steps.
1. -we examine the 20-year time series (1th January 2000 to 31st March 2020) of daily accumulated precipitation expressed as one average value representative of each area under consideration. Using the same dataset, we also generate a time series of 7-day antecedent precipitation. We then derive the distributions of daily and antecedent precipitation, extracting from each of the two all the values above the 95th percentile. For each extreme event, we resample the IMERG product to the same resolution of the PGA grid (approximately 1-km resolution) via an inverse distance weighted interpolator (Watson and Philip 1985); this being done for spatial consistency between rainfall and ground shaking proxies. 2. -we introduce the rainfall events extracted in the previous step as independent variables in a post-seismic RFIL susceptibility model. Specifically, we do this for the first post-seismic RFIL inventory associated with Reuleut earthquake. However, adding every possible realization of these extreme rainfall events is not suitable because some events may not be related to the RFIL; therefore, they can act as noise in the model. Moreover, these rainfall events could cause issues related to the presence of multi-collinearity (when two or more covariates are linearly related; see Amato et al. 2019) or high dimensionality of the covariate space (when the number of covariates is very large or even larger than the number of observations; see Castro Camilo et al. 2017). These issues can be addressed by implementing various regularization techniques, and among those, the least absolute shrinkage and selection operator (LASSO; Tibshirani 1996) has proven to be a valid tool (Camilo et al. 2017), which is why we chose it to reject the non-contributing rainfall events mentioned above.
LASSO constrains the number of covariates by adding a penalty term referred to as L 1 -norm, which corresponds to the sum of the absolute coefficients. The penalty acts on the likelihood shown in Eq. 5: where ℓ * is the new likelihood, ∑ N n¼1 jβ n j is the L 1 -norm and λ is introduced to balance the two terms. This procedure forces to zero the regression coefficients that have a negligible contribution to the model. This requires the term λ to be estimated. Its domain is confined between zero and infinity. More specifically, when λ = 0, then all the regression coefficients remain unchanged whereas, when λ → ∞, then all the regression coefficients are shrunk to zero. The parameter λ is commonly retrieved via cross-validation routines. For instance, the R (R-Team 2014) package glmnet (Friedman et al. 2009) by default examines 100 λ values, each one included in a 10-fold crossvalidation scheme. As for the metric the shrinkage is compared to, we have again selected the AUC (Hosmer and Lemeshow 2000). From the selected covariate subset and to maintain consistency with respect to the analyses carried throughout the manuscript, we then run 1000 bootstrap replicates to retrieve the sampling distribution of the ground motion and rainfall regression coefficients, together with the morphometric ones. This procedure is graphically summarized in step 2 (Fig. 6), and we also remind the readers that all the estimated coefficients are expressed in the same scale to ensure the reciprocal comparison.

Results
For both sites, we created SUs with median SU size of 0.46 km 2 (study area A) and 0.48 km 2 (study area B) (Fig. 7). These are the two respective spatial partitions we use to run a BLR for each inventory using six covariates (β Dist2Stμ , β Dist2Stσ , β Slopeμ , β Slopeσ , β PGAμ , β PGAσ ). Figure 8 shows the regression coefficients of examined covariates obtained for each of the modelled inventories. The boxplots are obtained from 1000 bootstrap replicates and represent the sampling distribution of the relative contribution of each covariate to the probability of landslide occurrence. Overall, slope steepness and PGA have positive regression coefficients for all co-seismic inventories. To examine how the regression coefficient of PGA changes through time, we also plotted the coefficient obtained for each temporal inventory as well as the corresponding AUC value (Fig. 9). Same as above, the sampling distribution of each parameter is retrieved by bootstrapping.
In study area A, for the co-and post-Reuleut susceptibility models, we used the Reuleut PGA map. The results indicate that the PGA has a positive weight on classifying a given slope unit as "landslide" instead of "non-landslide" given the choice of predictors (Fig. 9a). We observed a gradual decay in both regression coefficients (Fig. 9a) and AUC values (Fig.  9b) from the 1st to 3rd post-seismic inventories, indicating that the positive contribution of seismic shaking decreases over 2 years (from 14th December 2016 to 5th January 2019). Specifically, within 2 years, the median regression coefficient of PGA decreases from~0.6 to~0.2.
For the co-Reuleut landslide inventory, the regression coefficient of PGA and the AUC of the corresponding susceptibility model indicate a pattern we do not observe in the study area B (Fig. 9). In particular, both the AUC and the PGA regression coefficient of co-Reuleut inventory are lower than their counterparts calculated for the first post-Reuleut landslide inventory (Fig. 9a). Notably, the uncertainty-here calculated by bootstrapping for 1000 replicates-around the median of the PGA regression coefficients also changes through time, which is due to the difference in sample sizes (or to the fewer number of unstable SUs per temporal inventory).
As for the modelling performance, for each model, we calculated the bootstrapped AUC distribution, whose median values range from 0.64 to 0.72 overall (Fig. 9b). Specifically, the median AUC estimated for and from the co-seismic to the 2nd post-seismic landslide inventory is above 0.6. Conversely, almost the entire range of AUC calculated for the 3rd post-seismic landslide inventory is below 0.6, which is the limit of acceptable modelling performance (Hosmer and Lemeshow 2000). The goal of this work is not necessarily to obtain the best modelling performance, which could be achieved by including additional covariates, but rather to capture the relative contribution of PGA through time. In this regard, on the one hand, we did not add the rainfall signal because the examined post-seismic landslides were triggered by unknown rainfall events over a large time span. On the other hand, we did not include additional terrain properties to avoid the effect of variable interactions. By variable interaction, we refer to the fact that, in a multi-variate scheme, the Fig. 7 Overview of slope units generated for a study area A and b study area B. The size characteristics of the slope unit partition are reported in the boxplots, for each corresponding study area effect of one variable may depend on the value of another variable. So, by assuming a much larger covariate set than the one we used, our interpretation should have also accounted for inter-dependencies among covariates, which we avoided As a result, the AUC values show reasonable modelling performance for a model simply built on 6 covariates ( Fig.  9). Results also show that the lowest AUC value is obtained for the 3rd post-seismic RFIL inventory where we also identified the lowest positive weight on landslide occurrences associated with PGA layer. The statistically significant drop in AUC value from 2nd (AUC = 0.70) to 3rd (AUC = 0.64) may imply that the contribution of PGA impact may have decayed over time (lower β values) and that, in general, the model is less able to explain the unstable SU distribution over space.
For study area B, we again examined three post-seismic landslide inventories per each earthquake (Fig 9 c and d). By observing the post-Sulawesi susceptibility models, which are all mapped approximately 3 years after the Sulawesi earthquake, the mean PGA regression coefficient appears to be either very low and positive or negative for all the post-seismic inventories (Fig. 9c). This shows that the positive weight of PGA is low for RFIL susceptibility model conducted 3 years after an earthquake. As for the modelling performance, we identified a gradual decrease in AUC values from co-seismic to post-seismic periods (Fig. 9d) as we observed in study area A.
A similar situation is shown for the Palu earthquakes ( Fig.  9 c and d), where the susceptibility model of co-seismic landslide inventory shows a positive regression coefficient that decreases with time. In this case, the post-seismic landslide inventories are mapped within a year from the Palu earthquake. Here, although we observed a decay, the regression coefficients are still positive and relative higher than the Sulawesi case.
Our observation out of three cases presented above (i.e. Reuleut, Sulawesi and Palu) are consistent with each other. Up to 3 years, the PGA in each model contributes to increase the probability of landslide occurrence. However, the Kasiguncu case reveals a slightly different pattern in terms of variation in regression coefficient of PGA. In particular, similar to other cases, the regression coefficient of PGA is positive for the co-seismic phase and decays in post-seismic period (Fig. 9c). Nevertheless, the Kasiguncu is different from other cases since the median regression coefficient of PGA associated with the second post-seismic inventory became negative within 4 months after the earthquake. This is the most rapid decay among the examined cases. Also, the regression coefficient in the third post-seismic inventory, which is mapped a year after an earthquake, is relatively higher than its former post-seismic counterparts, and this is not consistent with other observations. We will elaborate this observation from an interpretative standpoint in "Discussion." Aside from these observations on regression coefficient of PGA, the AUC values show a gradual decay in each case examined in study area B (Fig. 8d).
To further elaborate on the role of PGA, we also run alternative models excluding the PGA itself (hence keeping only β Dist2Stμ , β Dist2Stσ , β Slopeμ , β Slopeσ ). We do so by using exactly the same sampling strategy to create 1000 bootstrap replicates. As a result, we could calculate the performance difference between the initial models we presented above and the alternative models. Figure 10 shows that the differences are always positive, which implies that the ground motion, when included, captures some spatial dependence in the landslide distribution which is otherwise unaccounted for in the new set of alternative models without PGA. As expected, relatively larger differences are associated to the models corresponding to co-seismic landslide inventories. The Reuleut case is the only exception, though this is not entirely surprising because the regression coefficients of PGA we calculated for the Reuleut inventories revealed the same variation from co-seismic to post-seismic phases.
In the next phase of the analyses, we focus on the first post-Reuleut RFIL inventory for which the estimates are shown in Fig. 9a highlights that the PGA is a positive contributing factor. More specifically, we run a BLR-based susceptibility model for this inventory, including all covariates we used in step 2 plus a time series of rainfall proxies (daily accumulated and 7-day antecedent precipitation) (Fig. 6). Each element of the time series and associated proxies has been identified in Fig. 11 a and b as an extreme event compared to the last 20 years of precipitation within the same season.
The red peaks in Fig. 11 show that a large number of extreme rainfall events could be responsible for the first post-Reuleut RFIL inventory. These rainfall events are graphically summarized in Fig. 12. The plots show the spatial distribution of normalized landslide density (i.e. number of landslides per km 2 rescaled from zero to one) (Fig. 12a), PGA (Fig.12b), the identified daily accumulated (Fig. 12c-i) and 7-day antecedent (Fig.  12j-n) rainfall. The figure shows some degree of similarity between the spatial distribution of RFIL (Fig. 12a) and the PGA (Fig. 12b), with a decreasing pattern from northeast to southwest. Conversely, there is no explicit agreement between the landslide distribution and any rainfall pattern.
Beyond the visual comparisons, to statistically isolate the most likely triggering events, we initially perform a LASSO variable selection. We assume that irrelevant rainfall patterns or events that are not responsible for the landslide initiation should not be selected by LASSO. We tested 100 λ values to shrink the parameter space, running a 10-fold cross-validation for each λ and storing the associated AUC values. Figure 13 a graphically summarizes this information. Two possible best λ choices are reported in the figure, corresponding to the vertical dashed lines. The dashed line to the left corresponds to the most conservative choice where a relatively large number of covariates is still allowed to keep the AUC performance stable or even to improve it. The dashed line to the right corresponds to the most penalized model with acceptable AUC performance with respect to the initial one featuring all the covariates. The corresponding sets of covariates are fourteen and six, respectively. From the second limit, the AUC rapidly decays. Here, we stress something peculiar out of the two covariate sets. The most conservative model selects 14 covariates out of which 8 (magenta rhombi) are rainfall proxies. These are associated with a negative regression coefficient (Fig. 13b). This is an unrealistic effect. The rainfall is the primary cause of the instability process; therefore, it is expected to be positive. Hence, a negative regression coefficient, even if statistically significant, should suggest that the specific rainfall covariate does not have any realistic or interpretable meaning. Therefore, our expert choice would be to remove the eight rainfall proxies with an unexplainable role. Interestingly, the second λ, which is more aggressive in constraining the parameter space, selects six covariates (blue rhombi) by removing the eight rainfall proxies. As a result, we consider the second λ to be the most reasonable out of the two. And, on the basis of the six selected covariates, we have run an additional set of analyses featuring 1000 bootstrapped replicates whose distributions of regression coefficients are plotted in Fig. 13c.
Results show that PGA μ has the largest contribution to the model (Fig. 13c) and that the most likely rainfall trigger corresponded to the 4th January 2017 event, which appears to be the third contributor out of the six, after slope σ .

Discussion
In this study, we primarily used the regression coefficients of ground shaking (i.e. PGA) to explore the relevance of the earthquake legacy effect on RFIL susceptibility assessments. We extent this consideration with a time-variant approach and examine how the ground shaking effect changes through time.

Temporal evolution of earthquake legacy effect
Our findings are meant to open up an interesting scientific discussion for earthquake legacy has not been explicitly investigated in the context of multi-temporal landslide susceptibility models. In the geoscientific community, the legacy of previously occurred earthquakes is already reported to be a factor which increases the landslide susceptibility of any area during the post-seismic phase (e.g. Tang et al. 2016;Yang et al. 2017;Kincey et al., 2021). However, there is no agreement on how long this elevated susceptibility will be maintained after an earthquake. Tian et al. (2020) indicate two main factors controlling this time period: the amount of co-seismic source material and precipitation pattern. In other words, they suggest that the elevated susceptibility could be nullified in a relatively short period of time if there is not much co-seismic source material but a strong precipitation pattern. Fig. 11 Precipitation regime represented by a daily accumulated precipitation and b 7-day antecedent precipitation for the area affected by the 2016 Reuleut earthquake. Yellow stars show the date of the earthquake. Vertical dashed black lines indicate the dates of the satellite imagery used to map RFIL. The mean and 95% confidence interval of daily and antecedent precipitation are calculated from a 20-year time series and are shown by black line and grey-shaded area, respectively. Red lines indicate the time period that precipitation is higher than the historic 95th percentile This is exactly the case for the first post-seismic landslide inventory of the Reuleut earthquake, where 7 extreme daily and 5 extreme 7-day antecedent precipitation events have hit the area after the Reuleut (Fig. 11). Here, we would like to point out that our LASSO implementation allowed us to isolate the most likely trigger out of the 12 possible proxies mentioned above.
It is important to highlight once more that the geoscience community has not yet found a reference time window after an earthquake for which the ground motion increases the landslide susceptibility. In this work, we cannot explicitly define such a time window, especially in a globally valid context. Our findings are certainly localized and should be framed in the context of probabilistic models. Specifically, our findings are representative for a tropical environment where large coseismic landslide deposits do not exist. Therefore, the contribution of PGA layer could be different through time in another environmental setting (e.g. Wenchuan), for instance, if large amount of co-seismic materials deposited on hillslopes and/or precipitation rate is relatively lower and not persistent compared to the area we examined. Moreover, we examined only a subset of the total area affected by earthquakes, and thus, our findings are representative for the areal boundaries encompassing multi-temporal inventories. But, they also point Fig. 12 Spatial distribution of RFIL on the site hit by the 2016 Reuleut earthquake overlaid by a normalized landslide density, b PGA map of the 2016 Reuleut earthquake, c-i daily accumulated precipitation maps for the days where precipitation amount is relatively high (see Fig. 11a) and j-n 7-day antecedent precipitation maps for the periods where precipitation amount is relatively high (see Fig. 11b) out at an interesting assumption that earthquake legacy may still play a significant role in the spatial distribution of RFIL. Overall, we can highlight two main interesting observations. First, the findings obtained in this study area may suggest a temporal scale where no effects can be inferred from the ground motion to RFIL inventories. In none of the cases, we capture a positive PGA effect on the susceptibility model (the mean regression coefficients appear to be negatives) 3 years after earthquakes. This could be a result of the probabilistic framework, and future studies with a longer time series may provide further evidence to support this hypothesis, both within the study area as well as in other geographic contexts. Also, examining the whole area instead of a subset of area affected by earthquake could provide a better insight into temporal evolution of earthquake legacy effect.
Second, our LASSO selection extracted the closest extreme rainfall to the earthquake. In addition to this, the PGA coefficient decreases through time (see Fig. 10 a and c). This might indicate that, after the first extreme rainfall stress, the landscape may release the majority of the "available" landslides conditioned by the shaking disturbance. In other words, during the initial stages when the ground motion effect on the susceptibility is at its highest after the earthquake, disturbed hillslopes can more easily fail, following the spatial footprint of the seismic shaking. As a result, the following rainfall events may trigger landslides in other slopes (different spatial distribution), leaving behind any sign of ground motion related patterns.
The aforementioned interpretation could also explain why the first post-seismic inventory associated with the Reuleut earthquake triggered more landslides than the earthquake itself. As a matter of fact, seismic shaking could cause reduction in soil and rock mass strength parameters in the near vicinity of area affected by an earthquake. This reduction does not necessarily cause failures on some hillslopes but decrease FoS in way that subsequent rainfall event(s) could cause instabilities. Given this explanation, we could assume that the extreme rainfall event occurred on 4th of January 2017 has triggered landslides on hillslopes which was already disturbed by 6th December 2016 Reuleut earthquake.
Notably, among the examined cases, the regression coefficients of PGA calculated for the post-Kasiguncu landslide inventories show a variation in time that is not fully consistent with three other cases (Reuleut, Sulawesi and Palu) (Fig. 10). Specifically, the mean regression coefficient of PGA is around zero for the second post-seismic inventory. This should not be surprising looking at the intensity of ground shaking observed in other cases. For instance, in the Reuleut earthquake, the maximum PGA inside the boundary of the study area is 0.50 g (Fig. 3), and the regression coefficient gradually decreases and became closer to zero within two years (Fig. 10). Also, in the Palu earthquake, the maximum PGA is 0.68 (Fig.  5), and PGA is still significant 1 year after an earthquake, whereas in the Kasiguncu earthquake, the PGA counterpart is 0.10 g (Fig. 5), which is much smaller than two other cases. Therefore, the earthquake legacy could also disappear soon after an earthquake. If the damaged hillslopes are not so widespread and the damaged ones already failed in the first few rainfall events right after an earthquake, the earthquake effects could be nullified rapidly. Nevertheless, even considering the Fig. 13 Panel a shows the results of the LASSO-penalization at varying λ, summarized in terms of AUC values plotted along the y axis and the number of resulting variables as a second axis to the top. The red squares correspond to the mean AUC values estimated in a 10-fold CV scheme for each λ. The error bar around corresponds to the 95% CI of the AUC distribution across CV replicates. The first vertical bar to the left corresponds to the best LASSO model (14 covariates out of 18), and the right vertical bar shows the models at which the AUC starts to decay quite fast because of the lack of covariate information (6 covariates out of 18). Panel b reports the LASSO-penalized regression coefficients for the best model shown in panel a. Panel c shows the range of regression coefficients calculated using 1000 bootstrap replicates. explanation given above, the distribution of the PGA coefficients estimated for the third post-seismic landslide inventory appears as anomaly (Fig. 10). This anomaly could be caused by the limitation of the dataset in terms of its areal coverage. Since we do not have cloud-free satellite scenes to create the multi-temporal inventories in a way to encompass the entire area affected by landslides, we ended up mapping only a part of it. This inevitably limited our capacity to monitor the evolution of landslides over time. Therefore, the mentioned anomaly might be a result of a similarity between the pattern of rainfall event(s) occurred in between second and third postseismic landslides. Moreover, the time gap between these two inventories is 1 year (27th September 2017-26th September 2018), and regrettably, clear identification of the rainfall event(s) triggered these landslides is extremely complex if possible at all, for such a long time window.

Further research directions
Overall, and in the validity domain of the study areas we considered, we probabilistically showed that the PGA map of the last strongest earthquake may be an informative predisposing factor for RFIL susceptibility models to be built after an earthquake. And yet, we do not have enough observations to retrieve a generic and applicable rule to do so within a specific time frame. Our work investigates the boundary conditions of the validity of this legacy effect, where we statistically show its presence at least for a minimum of 4 months and its absence within 3 years. However, for larger earthquakes such as Wenchuan or Gorkha, the persistence of elevated susceptibility conditions is longer due to the strong level of disturbance.
This observation needs to be checked, and our hypothesis further demonstrated in other studies to provide additional evidence of the temporal changes of the ground motion legacy, in various environmental settings. In fact, multi-temporal EQIL and RFIL inventories in different environmental settings are particularly important to provide a better vision regarding the earthquake legacy effect and its general validity. Also, even if the legacy effect is a concept upon which the community agrees and we were able to numerically capture it here, several additional questions still need to be addressed. For instance, a better constraint on the temporal decay of the earthquake legacy still needs to be defined. Here, we retrieved its persistence at a relatively coarse temporal resolution, which may need to be further increased to investigate the phenomenon even further. And, as mentioned before, the legacy signal we have retrieved here may also change from a landscape and/ or climatic setting to another.
It is worth noting that we could not disregard some possible sources of uncertainty in the data we used. For instance, a more robust analysis and possible interpretation could have been drawn mapping the whole landslide-affected area. In another source of uncertainty, this time of numerical nature could be due to the downscaling step we added to match the rainfall and PGA spatial resolution. This step can certainly have played a role in the model, but we expect it to have been smoothed when aggregating at the SU level. These sources of uncertainty may have propagated in the model and therefore biased our interpretation. However, we stress here that we have done our absolute best collecting the data and including a bootstrap simulation step to at least account for the model uncertainty.
In future studies, we will prioritize sites containing a large amount of co-seismic landslide deposits on hillslopes, in a high-relief mountainous environment where precipitation rates are low and strongly seasonal. This should also provide an alternative situation where the earthquake legacy may statistically behave in a different way, leading to different and contrasting interpretations. This is certainly an interesting research topic, but it also requires multi-temporal landslide inventories to be available with high temporal resolution.
Overall, the topic on earthquake legacy effect from earthquakes to subsequent RFIL is still in its infancy. This is the case because our knowledge is mostly limited by a few multitemporal landslide inventories. However, in the last two decades, the geoscientific community started to focus on this topic, due to the increasing availability of multi-temporal landslide inventories. In fact, access to this information is fundamental, and as the community progresses in collating other inventories and the landscape evolution through time, much more research could be developed. For instance, one of the most important research directions would be understanding the legacy effect from a mechanical standpoint, not only at the slope scale but also in relation to large populations of landslides. To do this, an even greater effort will need to be put into place to geotechnically characterize the subsurface. In other words, future improvements may be reached for instance by selecting a smaller study site and examining the mechanical properties of a landscape where ground motion and rainfall discharges are responsible for the reduction in the strength of hillslope materials.

Conclusions
In this study, we focus on a rarely investigated concept in landslide susceptibility assessments that looks at the coupled effect of earthquakes and rainfall in triggering landslides . We approach the concept by initially considering seismic shaking as a predisposing factor over two time series of landslide inventories, both triggered by earthquakes and rainfall. Subsequently, we focus on a specific case aiming to capture the legacy effect of ground motion on RFIL susceptibility and compare it to the precipitation effect. The analyses have all been carried out in a BLR framework to study the distribution of covariate effects. These have been estimated via 1000 bootstrap simulations. We also run an additional LASSO-penalized framework to select the best set of explanatory variables for RFIL.
Our findings suggest that the signal of ground shaking evolves through time in post-seismic periods. Specifically, we show that the PGA map of the last strongest earthquake can be used as a predisposing factor for RFIL occurring after the earthquake (i.e. up to 3 years) within the two study areas. This seismic proxy captures part of the spatial dependence in post-seismic RFILs, and it does so with a decreasing capacity through time. In fact, on the basis of the four cases we examined, the signal completely disappears. In other words, the mean or median regression coefficient in each case starts from a largely positive value and decreases down to near-zero values after approximately 3 years from the initial earthquake. Furthermore, we observe that the coupled effect of earthquakes and rainfall, if co-existing in a susceptibility model, improves the accuracy of the susceptibility estimates for subsequent RFIL, induced by extreme precipitation discharges. Regarding the latter, we statistically notice that the effect of extreme precipitation and the ground motion legacy is particularly relevant if it occurs shortly after the earthquake for we assume that the mobilization of weakened slopes is promoted in this stage of the landscape evolution. Actually, for the specific case we examined, we statistically observed that the ground motion legacy explains the RFIL distribution more than the actual rainfall trigger (the PGA mean and/or median regression coefficient is positive and larger than the rainfall proxy). This consideration is also statistically significant as the vast majority of the PGA coefficient distribution is larger and outside the 95% confidence interval of the rainfall coefficient distribution.
It should be noted that both the cases we examined are located in Indonesia, which are in a tropical environment. Therefore, our conclusions still need to be verified in other contexts. Future work could demonstrate similar ground motion and rainfall combined effects, which could represent a further advancement in regional/global near-real time statistically and physically based susceptibility modelling, or early warning systems. In fact, currently, none of these fundamental applications for planning and preparedness contextually feature the potential triggering effect of seismic shaking and precipitation. Therefore, current models lack the ability to explain the residual spatial dependence that the weakening effect of the ground motion may exert onto a given landscape, thus preparing it to potentially fail with a larger landslide rate than the expected. However, to achieve such tasks, the geotechnical characteristics of hillslope materials will need to be continuously monitored. Such monitoring context will make it possible to assess the actual effects of ground motion and its legacy effects on slope stability as the slope materials are exposed to rainfall discharges through time.
Availability of data and material The multi-temporal landslide inventories used in this study are available through NASA Landslide Viewer (https://landslides.nasa.gov).
Code availability Not applicable.

Declarations
Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.