Predicting Earthquake-Induced Landslides by Using a Stochastic Modeling Approach: A Case Study of the 2001 El Salvador Coseismic Landslides

: In January and February 2001, El Salvador was hit by two strong earthquakes that triggered thousands of landslides, causing 1259 fatalities and extensive damage. The analysis of aerial and SPOT-4 satellite images allowed us to map 6491 coseismic landslides, mainly debris slides and ﬂows that occurred in volcanic epiclastites and pyroclastites. Four different multivariate adaptive regression splines (MARS) models were produced using different predictors and landslide inventories which contain slope failures triggered by an extreme rainfall event in 2009 and those induced by the earthquakes of 2001. In a predictive analysis, three validation scenarios were employed: the ﬁrst and the second included 25% and 95% of the landslides, respectively, while the third was based on a k -fold spatial cross-validation. The results of our analysis revealed that: (i) the MARS algorithm provides reliable predictions of coseismic landslides; (ii) a better ability to predict coseismic slope failures was observed when including susceptibility to rainfall-triggered landslides as an independent variable; (iii) the best accuracy is achieved by models trained with both preparatory and trigger variables; (iv) an incomplete inventory of coseismic slope failures built just after the earthquake event can be used to identify potential locations of yet unreported landslides.


Introduction
In El Salvador, CA, landslides are among the most destructive natural processes and can cause mass fatalities and devastation in the built environment [1].They may occur either due to heavy rains (rainfall-induced landslides) or as a secondary event of an earthquake (earthquake-induced landslides).The extreme rainfall events that frequently affect the Central American region are indeed responsible for activating gravitational phenomena consisting of shallow and fast-moving flow landslides that may cause severe economic damage and fatalities [2].Additionally, seismic events can produce hundreds to thousands of gravitational phenomena which may cause more damage and victims than the earthquake itself, a situation that occurs in countries often affected by high-intensity phenomena.Indeed, among the various secondary effects of a seismic event (i.e., tsunamis, fires and liquefaction phenomena), 71.1% of deaths worldwide result from landslides [3].Landslides usually occur during or shortly after a seismic event and have different spatial distributions according to the slope gradient, lithology, earthquake magnitude and hypocentral and epicentral distances [4][5][6][7][8][9].A number of studies observed a correlation between earthquake intensity and landslide occurrence [10,11].Keefer identified a critical Richter magnitude of 4.0 and observed that the area affected by gravitational phenomena progressively grows up to 500,000 km 2 for an earthquake with a magnitude of 9.2 [10].However, this correlation may present exceptions.A landslide could indeed be initiated even by a weak shaking if a slope failure is imminent prior to the seismic event or if the material involved includes weakly cemented and/or incoherent rocks.Volcanic rocks, which are widely distributed across the country, are among these lithologies, as they are usually characterized by a higher mobility than non-volcanic rocks due to differences in their material properties such as granularity, collapsibility, and water content [12][13][14][15][16].These geological conditions can significantly increase the frequency of gravitational phenomena triggered and the areas affected by landslides when compared to areas that experienced comparable magnitudes [17,18].
Damage due to coseismic landslides is closely related to economic development and the ability to prevent and reduce disasters.Indeed, developed countries have the economic and scientific resources to plan and manage structures designed to contain blocks of rock and stop or divert both slow and rapid flows [19,20].They may also activate direct or indirect protection and alarm systems to favor a prompt and rapid general evacuation.Projecting effective defensive engineering structures and alarm systems requires knowing where landslides are more likely to occur.This information, which is also fundamental for land use planning aimed at mitigating landslide risk, is provided by landslide susceptibility maps.
Landslide susceptibility reflects the relative probability of slope failure occurrence in a given area [21,22].The mapping of landslide susceptibility can be achieved by using various methods, among which the stochastic approach is the most popular choice [23][24][25][26][27][28].This methodology is based on the hypothesis that "the past is the key to the future" [29][30][31].Therefore, future landslides are more likely to be triggered under the same environmental conditions that have produced past and present slope failures.The first step of landslide susceptibility modeling using a stochastic approach is to select a mapping unit, which usually corresponds to square grid cells derived from digital elevation models (DEMs) or to slope units limited by streams and water divides.Then, a value of the dependent and independent variables is assigned to each mapping unit.The dependent variable is usually dichotomous (0: absence of landslide; 1: presence of landslide).The values of the selected independent variables (also known as covariates or predictors), which serve as proxies for the preparatory factors that mainly control landslide occurrence in the study area, are usually derived from DEMs and from already available thematic maps, such as geological or land use maps.Finally, multivariate statistical techniques or machine learning algorithms are usually employed to model the relationships between the landslide spatial distribution and the variability of the predictors to calculate a probability value of landslide occurrence for each mapping unit.The most common statistical and machine learning methods employed for this aim are logistic regression [22,28,[32][33][34], discriminant analysis [30,35,36], support vector machines [37,38] or random forests [39].Comprehensive reviews of the most popular methods for landslide susceptibility zonation can be found in [30,[40][41][42].Landslide susceptibility zonation usually relies on independent variables that reflect only preparatory factors, assuming that the geographical distribution of the landslide inventories is not affected by the pattern of the trigger (either rainfall or earthquake), which is thus considered spatially homogeneous.However, if a homogenous trigger pattern can be assumed for single and relatively small areas, this hypothesis does not apply when modeling large areas or multiple study areas that may experience different intensities of the trigger.In other words, when landslide susceptibility is assessed at a regional scale or in different areas spread over a large territory, it is necessary to consider not only variables that express preparatory factors but also variables that represent the pattern and intensity of the trigger in relation to the return period.
This paper focuses on predicting the landslides triggered by two earthquakes of magnitudes 7.7 and 6.6 which occurred in El Salvador in January and February 2001, respectively, and triggered thousands of slope failures and caused more than 1150 fatalities.
The main aim of this study is to evaluate the ability of rainfall-and earthquake-induced landslide susceptibility models to predict the geographical distribution of the coseismic slope failures that occurred in 2001 under three calibration and validation scenarios.The first is a scenario in which most landslides are known (i.e., 75%), and the second is a scenario in which only a minor portion of the gravitational phenomena that occurred are identified (i.e., 5%).A comparable and more than acceptable predictive ability of these scenarios might suggest that an inventory of slope failures built just after the triggering event, even if incomplete, may help rapidly generate reliable landslide susceptibility maps.To further evaluate our approach, we evaluated a third scenario based on a k-fold spatial cross-validation [37,43,44].
The experiment was performed in fourteen sectors of El Salvador, accounting for a total area of 397 km 2 .Eight sectors (92 km 2 ) were strongly hit by the January earthquake, whereas most of the landslides triggered by the February earthquake occurred in the other six sectors (305 km 2 ).The modeling strategy was based on both rainfall-and earthquake-triggered landslide susceptibility models.These were based on statistical relationships between landslide inventories and multiple predisposing and triggering factors, established using multivariate analysis regression splines (MARS), which demonstrated excellent landslide predictive ability in previous papers [22,[45][46][47][48][49][50].The first model was calibrated using an inventory of thousands of landslides triggered on 7-8 November 2009 by the combined action of Hurricane Ida and the low-pressure system 96E in a sector extending 287 km 2 along the slopes of the San Vicente Volcano, El Salvador [49].The rainfall-triggered model, which was trained by eleven preparatory factors, was applied to all fourteen study areas to verify its ability to predict the position of the coseismic landslides that occurred in January and February 2001.On the other hand, the earthquake-induced landslide predictive models were calibrated using both preparatory and trigger variables, with the latter reflecting the intensity of ground shaking.
The article is organized as follows.Section 2 describes the available data and the methods employed to map the coseismic landslides and to calibrate and validate the predictive models.Section 3 describes the inventory of the coseismic landslides and shows the results of models' training and testing.The results are discussed in Section 4, and the main conclusions of the study are reported in Section 5.

Study Areas
The complex geodynamic context in which the small country of El Salvador is located is characterized by the convergence of the Cocos plate and the Caribbean plate through the Middle America Trench at a rate of about 7 cm/year [51] (Figure 1a).This active tectonic environment is responsible for several active volcanoes in the country which are aligned in the WNW-ESE direction and have frequently erupted and deposited widespread and poorly consolidated tephra in many parts of the country [13].These lithologies, in combination with acid pyroclastites and acid and basic effusive rocks, ensure that volcanic rocks represent 95% of the outcropping lithologies in the country (Figure 1b).
The combination of the geodynamic contexts responsible for medium and high earthquakes, steep slopes due to rugged volcanic ranges, incoherent rocks and a subtropical climate characterized by heavy storms, cause a high probability of the occurrence of earthquake-and rainfall-induced landslides in most of the country.The consequences of these gravitational phenomena are devastating for the country, which is the smallest country in Central America but has the highest population density (310 inhabitants/km 2 ).Indeed, the fertility of the soils around the volcanoes, which are oriented WNW-ESE due to the tectonic assessment [17], has generated a situation in which a large part of the population lives in the areas of the country most susceptible to the gravitational phenomena (Figure 1b).As a confirmation of this, El Salvador has the highest ratio of destroyed houses to affected risks in Central America [52].This work was carried out in fourteen sectors along the volcanic chain (Figure 2), which were among the areas most affected by the two intense earthquakes that occurred in the first two months of 2001.
ISPRS Int.J. Geo-Inf.2023, 12, x FOR PEER REVIEW 4 of 3 lives in the areas of the country most susceptible to the gravitational phenomena (Figure 1b).As a confirmation of this, El Salvador has the highest ratio of destroyed houses to affected risks in Central America [52].This work was carried out in fourteen sectors along the volcanic chain (Figure 2), which were among the areas most affected by the two intense earthquakes that occurred in the first two months of 2001.ISPRS Int.J. Geo-Inf.2023, 12, x FOR PEER REVIEW 4 of lives in the areas of the country most susceptible to the gravitational phenomena (Figu 1b).As a confirmation of this, El Salvador has the highest ratio of destroyed houses affected risks in Central America [52].This work was carried out in fourteen sectors alon the volcanic chain (Figure 2), which were among the areas most affected by the two inten earthquakes that occurred in the first two months of 2001.The first eight study areas, numbered from 1 to 8, were selected because of the high number of landslides triggered by the January 2001 earthquake.Sectors 1 to 3 are in the southwestern part of the country, close to the San Salvador volcano, whereas sectors 4 to 8 are located in the southeastern portion of El Salvador.Sectors 9 to 14 are placed in the southern central part of the country and were selected because they were close to the epicenter of the February 2001 earthquake and were strongly affected by coseismic landslides.Table 1 reports the extent, average and standard deviation of the altitude and slope angle and the main outcropping lithology of all the study areas.

The 2001 Earthquakes
El Salvador has been severely affected by medium-and high-intensity earthquakes that have triggered numerous and devastating landslides.Bommer et al. [17] report that the country is experiencing earthquakes from two major seismic sources.The first source of seismicity originates from deep earthquakes in the Benioff-Wadati zones of the subducted Cocos plate beneath the Caribbean plate, which are characterized by relatively high magnitude values (often with surface wave magnitude, Ms, and moment magnitude, Mw, values of >6.0), capable of involving large areas of El Salvador.Numerous seismic phenomena with this genesis were witnessed in the last century, many of which were particularly destructive and had dramatic impacts on the country.Among these phenomena were the September 1915 event, Ms 7.8, which caused destruction in the country's western area, and the Ms 7.3 earthquake in June 1982, which caused hundreds of earthquakeinduced landslides and had 49 victims.The second source of seismicity is a zone of upper-crustal earthquakes that coincides with the array of volcanoes and is orientated WNW-ESE [17].These intraplaque earthquakes occur at relatively lower hypocentral depths.Since the main Salvadoran cities are located just along the volcanoes' alignment, it is not hard to determine the catastrophic consequences of the seismic events that can occur, often with more significant damage than the seismic phenomena originating in Benioff-Wadati zones.Among the most significant events of this origin in the last century are the earthquake of 8 June 1917 (Ms 6.7), which had 101 victims, the 20 December 1936 (Ms 6.1) event, which resulted in more than 200 victims, and the 10 October 1986 event (Ms 7.5), which produced more than 1500 fatalities.
On 13 January 2001, the country was hit by a violent earthquake, originating from the first source of seismicity, at a depth of 60 km with a magnitude of 7.7 Mw (Figure 3).This event caused thousands of gravitational phenomena, injuring 5565, causing 844 fatalities and destroying over 100,000 buildings.The differences in damage and fatalities between this event and the one in 1982, which was of a similar magnitude, are related to the increase in the elements at risk.Indeed, a rapid expansion of urban areas occurred during this time span, which continued to grow mainly toward the flanks of the volcanoes, as revealed by the difference between the topographical maps of the country dating back to the year 1981 and aerial photos from 2001.
Mw, values of > 6.0), capable of involving large areas of El Salvador.Numerous seismic phenomena with this genesis were witnessed in the last century, many of which were particularly destructive and had dramatic impacts on the country.Among these phenomena were the September 1915 event, Ms 7.8, which caused destruction in the country's western area, and the Ms 7.3 earthquake in June 1982, which caused hundreds of earthquake-induced landslides and had 49 victims.The second source of seismicity is a zone of upper-crustal earthquakes that coincides with the array of volcanoes and is orientated WNW-ESE [17].These intraplaque earthquakes occur at relatively lower hypocentral depths.Since the main Salvadoran cities are located just along the volcanoes' alignment, it is not hard to determine the catastrophic consequences of the seismic events that can occur, often with more significant damage than the seismic phenomena originating in Benioff-Wadati zones.Among the most significant events of this origin in the last century are the earthquake of 8 June 1917 (Ms 6.7), which had 101 victims, the 20 December 1936 (Ms 6.1) event, which resulted in more than 200 victims, and the 10 October 1986 event (Ms 7.5), which produced more than 1500 fatalities.
On 13 January 2001, the country was hit by a violent earthquake, originating from the first source of seismicity, at a depth of 60 km with a magnitude of 7.7 Mw (Figure 3).This event caused thousands of gravitational phenomena, injuring 5565, causing 844 fatalities and destroying over 100,000 buildings.The differences in damage and fatalities between this event and the one in 1982, which was of a similar magnitude, are related to the increase in the elements at risk.Indeed, a rapid expansion of urban areas occurred during this time span, which continued to grow mainly toward the flanks of the volcanoes, as revealed by the difference between the topographical maps of the country dating back to the year 1981 and aerial photos from 2001.
Exactly one month later, on 13 February 2001, El Salvador experienced another intense earthquake, which increased the destruction experienced by the country.This event, originating from the second source of seismicity, had a Mw of 6.6, with a hypocenter placed at a depth of 10 km.This earthquake resulted in 315 fatalities and injured 3,399 people due to ground-shaking and landslides (Figure 3).Exactly one month later, on 13 February 2001, El Salvador experienced another intense earthquake, which increased the destruction experienced by the country.This event, originating from the second source of seismicity, had a Mw of 6.6, with a hypocenter placed at a depth of 10 km.This earthquake resulted in 315 fatalities and injured 3399 people due to ground-shaking and landslides (Figure 3).Nearly 45,000 homes were destroyed and 16,000 were seriously damaged.This earthquake affected areas comparable to those affected by the 1936 event, which was of the same genesis.Once again, the increased risk factors have led to more damage and more victims.
Most of the landslides triggered by the 2001 earthquakes were shallow debris flows; however, there were also numerous slides and debris and rock falls, according to Cruden and Varnes' classification of landslides [53].Moreover, both earthquakes triggered a few large and deep landslides that were highly damaging.
There are several doubts about whether these two earthquakes were related.However, some studies [17,54] reported that normal faulting subduction earthquakes in Central America tend to be followed by either large thrust events or shallow intraplate events within four years.This trend can also be observed in El Salvador, where earthquakes belonging to the second type of seismic genesis have occurred even after a month, as in the case of the earthquakes of 2001.To mention other examples of this trend, Bommer et al. [17] cited subduction earthquakes that occurred in 1915, 1932 and 1982 and were followed by intra-crustal earthquakes dated 1917, 1936 and 1986, respectively.
El Salvador is in a tropical humid climate system characterized by a dry season from November to April and a wet season from May to October in which about 75% of the annual precipitation falls.The earthquakes of 2001 occurred in the first two months of the year and thus during the dry season, as confirmed by the precipitation data recorded close to the study areas.Furthermore, the cumulative precipitation in the 12 months preceding the earthquakes of 2001 was lower than the annual average recorded at the same gauge stations in the 29 years before the earthquakes (Table 2).In addition, November 2000 was the last month that experienced rainfall, which was lower than the average, and unlike the averages of previous years, the following months until the earthquakes of 2001 were without rain.These data suggest that the water content did not play a decisive role in the genesis of gravitational phenomena due to the absence of or low water pore pressure in static conditions and the presence of a relatively high apparent cohesion in pyroclastic soils due to suction.

Mapping Strategy
Landslide inventories are valuable tools for investigating landscape susceptibility to these phenomena.When investigating the slope's response to a specific triggering event, the archive must exclude the landslides that occurred before and after that event.
In stochastic and empirical assessment of landslide susceptibility, the presence or absence of a gravitational phenomenon within a selected mapping unit is employed as dependent variable.Thus, the accuracy and degree of completeness of the inventory largely affect the predictive performance of the models [47, [55][56][57][58].
The 2001 landslide inventories available to date are incomplete, as they only include the main gravitational phenomena; therefore, the first step of this study was to prepare a complete inventory of the slope failures triggered by the two earthquake events in the study areas.To this aim, aerial photos provided by the CNR (Centro Nacional de Registros de El Salvador-Instituto Geográfico y del Catastro Nacional) were employed.Regarding the January earthquake, 136 aerial photos covering an area of 91.85 km 2 were analyzed, while 313 aerial photos acquired over an area of 305.72 km 2 were employed to identify the landslides triggered by the February earthquake (Figure 4).while 313 aerial photos acquired over an area of 305.72 km 2 were employed to identify the landslides triggered by the February earthquake (Figure 4).Recognizing and mapping landslides in tropical areas where vegetation recovers rapidly requires aerial/satellite images taken not long after the event.The employed aerial photos were taken in the days just after the earthquakes, when the displaced materials had not yet been removed or covered by vegetation.The aerial images cover the areas of the country that were the most heavily affected by the 2001 earthquakes.These correspond to sectors of the eastern and western parts of the country in the case of the January earthquake, and to the central sector in the case of the February earthquake.To exclude slope failures triggered by a previous seismic or rainy event from the analysis, we analyzed a SPOT-4 satellite image dated 27 November 2000 (Scene ID: 4 611-323 00-11-27 16:42 2M) with 10 m spatial resolution (European Space Agency).In addition, in the areas covered by the February 2001 aerial photos, SPOT-4 satellite images from 21 January (Scene ID: 1 611-323 01-01-21 16:44 1P) were also analyzed to identify those landslides triggered by the January earthquake.It is worth noting that the aerial and satellite images allowed us to detect only landslides due to slope collapses; thus, landslide-producing fracturing and tension cracks from short coseismic slides were not included in the inventories.
In this study, the landslides were recognized, and their bodies were mapped as vector polygons.Then, the highest point of each landslide crown was identified, assuming that the environmental properties of these points, known as LIPs (landslide identification points) [23,[59][60][61], reflect the pre-failure conditions that favored the initiation of the landslides.

Dependent and Independent Variables
Landslide susceptibility assessment involves dividing the study area into mapping units, which usually correspond to grid cells or slope units [29,62,63].As grid cell partition is largely and effectively employed in landslide susceptibility mapping [8, 64,65], this strategy was employed in our experiment.The study areas were resampled into 10 m pixels Recognizing and mapping landslides in tropical areas where vegetation recovers rapidly requires aerial/satellite images taken not long after the event.The employed aerial photos were taken in the days just after the earthquakes, when the displaced materials had not yet been removed or covered by vegetation.The aerial images cover the areas of the country that were the most heavily affected by the 2001 earthquakes.These correspond to sectors of the eastern and western parts of the country in the case of the January earthquake, and to the central sector in the case of the February earthquake.To exclude slope failures triggered by a previous seismic or rainy event from the analysis, we analyzed a SPOT-4 satellite image dated 27 November 2000 (Scene ID: 4 611-323 00-11-27 16:42 2M) with 10 m spatial resolution (European Space Agency).In addition, in the areas covered by the February 2001 aerial photos, SPOT-4 satellite images from 21 January (Scene ID: 1 611-323 01-01-21 16:44 1P) were also analyzed to identify those landslides triggered by the January earthquake.It is worth noting that the aerial and satellite images allowed us to detect only landslides due to slope collapses; thus, landslide-producing fracturing and tension cracks from short coseismic slides were not included in the inventories.
In this study, the landslides were recognized, and their bodies were mapped as vector polygons.Then, the highest point of each landslide crown was identified, assuming that the environmental properties of these points, known as LIPs (landslide identification points) [23,[59][60][61], reflect the pre-failure conditions that favored the initiation of the landslides.

Dependent and Independent Variables
Landslide susceptibility assessment involves dividing the study area into mapping units, which usually correspond to grid cells or slope units [29,62,63].As grid cell partition is largely and effectively employed in landslide susceptibility mapping [8, 64,65], this strategy was employed in our experiment.The study areas were resampled into 10 m pixels resulting from rescaling a 5 m cell DEM from the European Space Agency (ESA).To each of these pixels, a value for the dependent and independent variables was assigned.The cell size was selected to find a compromise between the spatial resolution and computation time needed to train and test the models [64,66], taking into account that a 10 m cell partition was successfully employed in previous studies [33,67].
In this study, the dependent variable is dichotomous, taking the value 1 when a pixel hosts one or more LIPs and the value 0 when it does not intersect a LIP or a landslide body.Independent variables were selected according to data availability and their supposed direct or indirect relationships with the spatial distribution of the slope failures that occurred in the study areas [68].The following eleven predictors, which reflect preparatory factors, were selected: Northness (N), Eastness (E), Catchment Area (C_AR), Convergence Index (C_IND), Downslope Curvature (D_CUR), Elevation (ELE), Geology (GEO), Normalized Difference Vegetation Index (NDVI), Slope (SLO), Topographic Positioning Index (TPI), Topographic Wetness Index (TWI) and Upslope Curvature (U_CUR).N and E were obtained by calculating the cosine and sine of slope aspect, respectively.The NDVI was calculated from images dated 1 January 2001 and 2 February 2001 to analyze the areas affected by the January and February 2001 events, respectively.The terrain attributes were extracted from the DEM, which is dated 2017 and is thus more recent than the 2001 earthquake events.Therefore, to avoid training and testing the models with topographic conditions altered by the slope failures, the landslide points coincided with the highest points of each crown (i.e., LIP) and the non-landslide points were randomly selected outside the landslide areas.Additionally, two trigger (i.e., depending on the earthquake characteristics) variables were also employed as covariates, namely: Peak Ground Accelerations (PGAs) and the distance from the hypocenter (DIST).The first was obtained by rasterizing the PGA maps from both earthquakes from the USGS ShakeMap website [69], while the second by producing a distance raster from the epicenter coordinates known through the USGS website.
The absence of multicollinearity among the predictors was evaluated by calculating the Variance Inflation Factor (VIF) (e.g., [32,70]) and by applying the rule of ten ( [26,71] and references therein).As all the predictors exhibited a VIF value < 10, they were used to build the landslide predictive models.
Table 3 summarizes the variables employed to model the landslide spatial distribution.

Modeling Technique
The likelihood of landslide occurrence at each 10 m grid cell of the study area was calculated using multivariate adaptive regression splines (MARS; [72]) as a modeling technique.MARS is a non-parametric regression method that is useful for handling complex data and is able to model non-linear relationships between dependent and independent variables, splitting their range into an optimized number of linear branches.The extreme values of these intervals, called knots, allow an optimal linear function known as a basis function (BF).The result is a weighted sum of terms which includes a basis function or the product of two or more of them.The MARS regression function is expressed as (1): where y is the dependent variable predicted by the function f (x), α is the constant, each h i (x) is a BF and β i its coefficients, and N is the number of functions.

Calibration and Validation Strategy
In this study, the landslide predictive models were calibrated and validated using five datasets (Figure 5; Table 4), which correspond to the fourteen study sectors (JAN, FEB and JF), to the area of the San Vicente (SV) volcano landslide inventory prepared by Mercurio et al. [49] and to the intersection between the latter and sectors 9 and 13 (FSE).

Calibration and Validation Strategy
In this study, the landslide predictive models were calibrated and validated using five datasets (Figure 5; Table 4), which correspond to the fourteen study sectors (JAN, FEB and JF), to the area of the San Vicente (SV) volcano landslide inventory prepared by Mercurio et al. [49] and to the intersection between the latter and sectors 9 and 13 (FSE).These datasets are matrices in which rows correspond to the pixels of a specific area, These datasets are matrices in which rows correspond to the pixels of a specific area, and columns report the values of the dependent and independent variables calculated for each pixel.The dataset of rainfall-triggered landslides (i.e., SV) includes only preparatory covariates, whereas those of the coseismic landslides (i.e., JAN, FEB, JF and FSE) report the values of both preparatory and trigger variables.The MARS predictive models were prepared using different calibration datasets or sets of predictors (Table 5).The models named M1 were trained with the dataset SV, which includes rainfall-induced landslides and eleven independent variables reflecting preparatory factors.The M1 model calibrated with the dataset SV was then transferred to all fourteen study sectors of El Salvador, calculating a probability (PSV) of landslide occurrence for each pixel of these areas.The PSV values were then employed as new predictors when training the M3 model.The M2, M3 and M4 models were calibrated with the earthquake-induced landslides of January and February 2001.The variables PGA and DIST were used in the M2 model as the only covariates and in the M3 model in addition to PSV.The M4 model contains both preparatory and triggering independent variables.The predictive ability of the M1 models was evaluated by comparing the PSV values and the presence or absence of coseismic landslides in the datasets JAN, FEB, JF and FSE.
The MARS models were calibrated and validated using three strategies.The first and second strategies consisted of the following steps.First, we randomly selected ten subsets from each dataset, which included all the event cells and the same number of non-event cells.The subsets of cells extracted from the SV dataset were used to train ten M1 models.Second, we extracted two calibration samples from the ten subsets of all the datasets, excluding SV, by randomly selecting 75% and 5% of the event cells, respectively.The remaining event cells (25% and 95%, respectively) were included in the validation samples.Both calibration and validation samples were completed by randomly selecting the same number of non-event cells.Thus, this second step produced twenty calibration and twenty validation samples from the datasets JAN, FEB, JF and FSE.The two ratios of calibration/validation samples were chosen in order to evaluate the predictive performance of the MARS models under two different scenarios: in the first scenario, the position of a large part (i.e., 75%) of the earthquake-induced landslides is known; the second scenario simulates the situation just after the earthquake occurrence, when only a very small part (i.e., 5%) of the landslides had been identified.The probability of landslide occurrence was calculated for each cell of the validation samples by averaging the score obtained from the ten MARS model runs (one for every calibration sample).This strategy was applied to increase the stability of the landslide predictions and mitigate the rare-events issue [22,32,80,81].
To further evaluate the reliability of the models, a k-fold spatial cross-validation (CV) strategy was also adopted [44,82].According to this approach, the models are trained for k-1 combined samples at time and validated using the remaining sample.This procedure is repeated k times.Unlike the classic k-fold cross-validation approach, a sample is extracted from each of k different areas.In this study, five (k = 5) different sampling areas were identified by applying the k-means clustering algorithm.Similarly, to the first and second calibration and validation strategy, the spatial CV was repeated ten times in order to assess the robustness of the models.
For all the validation samples, a receiver operating characteristics (ROC) curve [42,48,71,83] was prepared using the average MARS score to predict the value of the dependent variable at each pixel (0: LIP absence; 1: LIP presence).The ROC curve plots sensitivity (true positive rate) against 1-specificity (false positive rate).Then, the area under the ROC (AUC) curves were calculated, thus obtaining ten AUC values for the prediction of the datasets JAN, FEB, JF and FSE for each of the models M1 to M4 and validation scenarios.This procedure led to a total of 480 AUC values.An AUC value close to 0.5 reveals a random prediction, whereas a value of 1.0 reflects a perfect ability to discriminate between event and non-event pixels.AUC values larger than 0.7, 0.8 and 0.9 were considered as acceptable, excellent and outstanding, respectively [84].The Wilcoxon signed-rank test was employed to identify significant differences between the predictive skill performance of the models by setting the level of significance at 0.01.

Landslide Susceptibility Maps
The models M1 to M4 were used to generate a landslide susceptibility map for each of the four datasets.To this aim, we employed the MARS models calibrated with the samples including 75% of the LIPs.For each pixel of the datasets JAN, FEB, JF and FSE, a probability (P) of landslide occurrence was calculated for the models M1, M2, M3 and M4 by averaging the score of their ten repetitions.Since the models were trained on samples including the same number of LIP and non-LIP locations, P should be interpreted as the spatial relative probability of landslide occurrence.

The 2001 Earthquake-Induced Landslide Inventories
The landslide inventory of the dataset JAN was obtained by analyzing eight sets of aerial photos taken a few days after the earthquake of 13 January 2001.These images cover eight inhabited sectors, extending for 91.85 km 2 (Figure 6), which experienced the most significant damage due to the seismic shaking and the consequent landslides.These were mainly rotational slides, debris/earth flows and earth/rock falls.
The inventory consists of 997 LIPs, each hosted by a 10 m pixel.The percentage of event cells (i.e., cells hosting at least one LIP) compared to the total number of pixels (1,260,811) of the study area is equal to 0.08%, while there are 9556 cells intersecting the landslide polygons, corresponding to 0.76% of the total pixels.The JF inventory includes the LIPs and landslides caused by the January earthquake in the areas covered by the aerial images taken after the February event (sectors 9 to 14, counted from west to east, are Ilopango south and Rio Jiboa, S. Cruz Michapa, Cojutepeque, El Rosario, San Vicente and Tecoluca, respectively; red polygons of Figure 6).These landslides were identified by analyzing images from the SPOT-4 satellite dated 21 January 2001.In these sectors, with relatively low January PGA values, 123 LIPs were mapped, each placed in one 10 m cell and accounting for 0.002% of the total number of cells of the area under examination.When landslide bodies are considered, unstable pixels cover 0.2% of the area.

The 2001 Earthquake-Induced Landslide Inventories
The landslide inventory of the dataset JAN was obtained by analyzing eight sets of aerial photos taken a few days after the earthquake of 13 January 2001.These images cover eight inhabited sectors, extending for 91.85 km 2 (Figure 6), which experienced the most significant damage due to the seismic shaking and the consequent landslides.These were mainly rotational slides, debris/earth flows and earth/rock falls.The inventory consists of 997 LIPs, each hosted by a 10 m pixel.The percentage of event cells (i.e., cells hosting at least one LIP) compared to the total number of pixels (1,260,811) of the study area is equal to 0.08%, while there are 9556 cells intersecting the landslide polygons, corresponding to 0.76% of the total pixels.The JF inventory includes the LIPs and landslides caused by the January earthquake in the areas covered by the aerial images taken after the February event (sectors 9 to 14, counted from west to east, are Ilopango south and Rio Jiboa, S. Cruz Michapa, Cojutepeque, El Rosario, San Vicente and Tecoluca, respectively; red polygons of Figure 6).These landslides were identified by When the JAN and JF inventories are merged, the percentages of LIPs and landslide pixels in relation to the total area (sectors 1 to 14) are significantly reduced (0.017 and 0.30, respectively) due to the low PGA values within sectors 9 to 14.However, the landslides triggered in these sectors, despite being about one tenth of the number of the JAN landslides, cover a significantly higher number of pixels (Table 6).This can be explained considering that the steepest slopes are characterized by thick deposits of poorly consolidated late Pleistocene and Holocene Tierra Blanca rhyolitic tephras that erupted from Ilopango caldera, an area particularly susceptible to the initiation of gravitational phenomena [47,61].In fact, in that region, the January earthquake resulted in a low number of debris flows and rotational slides; however, they had a large extent.This is also demonstrated by the number of gravitational phenomena that were larger than 10,000 m 2 : two and twenty-two in the JAN and the JF datasets, respectively.However, two gravitational phenomena, which both occurred in sector 3 (Las Colinas neighborhood of Santa Tecla, west of the capital San Salvador), caused the worst damage and the largest number of victims.One of these was a complex landslide (rotational slide and earth flows) that covered an area of about 200,000 m 3 in which approximately 585 people lost their life.This landslide, which sadly become famous for the devastation it generated, was about 800 m long and 150 m wide, with a 50 m high escarpment activated on the north side of El Balsamo Ridge, which is composed of the andesitic cinders and some interbedded tephra of the El Balsamo Formation [17].The second event, located approximately 2.5 km from the first landslide, was a debris flow about 1500 m long which also began from the north side of El Balsamo Ridge and reached and destroyed a sector of the Carrera Panamericana, the most crucial road in the country which connects the capital with the western part of El Salvador.The highway interruption hindered the arrival of aid, creating enormous damage to the population.However, the absence of small landslides in sectors 9 to 14 compared to that of sectors 1 to 8 (41.42%) could also be related to the different resolution of the images used for mapping, as phenomena with an extension of less than 100 m 2 are hardly detectable in a satellite analysis with a 10 m resolution.The seismic phenomenon of 13 February triggered thousands of landslides in the country's central sector.Many landslides occurred where slopes were covered by thick and slightly consolidated deposits of the rhyolitic tephras of the Tierra Blanca Formation, which originated from eruptions from the Ilopango Caldera.
Landslide mapping was carried out using aerial photos taken a few days after the event in areas that experienced a very high PGA and suffered the most significant damage (sectors 9 to 14, Figure 7).The landslides that were also included in the JF inventory were discarded except for those that have been reactivated.The seismic phenomenon of 13 February triggered thousands of landslides in the country's central sector.Many landslides occurred where slopes were covered by thick and slightly consolidated deposits of the rhyolitic tephras of the Tierra Blanca Formation, which originated from eruptions from the Ilopango Caldera.
Landslide mapping was carried out using aerial photos taken a few days after the event in areas that experienced a very high PGA and suffered the most significant damage (sectors 9 to 14, Figure 7).The landslides that were also included in the JF inventory were discarded except for those that have been reactivated.The FEB inventory includes 5371 LIPs intersecting 5344 pixels.The percentage of event cells of the FEB dataset is 0.09%.The extent of the landslide bodies is 123,488 pixels, which accounts for 2.12% of the study area.This seismic event triggered landslides similar to those of the JAN dataset, but their frequency and size were higher (Table 7).
Table 7. Amounts and characteristics of the gravitational phenomena of landslides belonging to The FEB inventory includes 5371 LIPs intersecting 5344 pixels.The percentage of event cells of the FEB dataset is 0.09%.The extent of the landslide bodies is 123,488 pixels, which accounts for 2.12% of the study area.This seismic event triggered landslides similar to those of the JAN dataset, but their frequency and size were higher (Table 7).Among the many large-scale phenomena, large rotational flows of debris were mapped in the area close to the epicenter, corresponding to areas near the Jiboa River, which fortunately have low population densities.As a matter of fact, thousands of gravitational phenomena have been detected in the Tierra Blanca Formation (GEO5), some of which caused the damming of the river for hundreds of meters, causing the formation of an artificial lake that could have caused a high hydrogeological risk downstream.
The FSE inventory was generated by overlapping the areas of the datasets FEB and SV.The intersection extended for 107 km 2 and included the entirety of sector 13 and 12.5% of sector 9 in the FEB study areas (Figure 8).The FSE dataset contains 1587 LIPs pertaining to the landslides triggered by the February earthquake, which were mainly earth/debris flows.Among the many large-scale phenomena, large rotational flows of debris were mapped in the area close to the epicenter, corresponding to areas near the Jiboa River, which fortunately have low population densities.As a matter of fact, thousands of gravitational phenomena have been detected in the Tierra Blanca Formation (GEO5), some of which caused the damming of the river for hundreds of meters, causing the formation of an artificial lake that could have caused a high hydrogeological risk downstream.
The FSE inventory was generated by overlapping the areas of the datasets FEB and SV.The intersection extended for 107 km 2 and included the entirety of sector 13 and 12.5% of sector 9 in the FEB study areas (Figure 8).The FSE dataset contains 1587 LIPs pertaining to the landslides triggered by the February earthquake, which were mainly earth/debris flows.The conditional density plots of Figure 9 reveal empirical relationships between landslide inventories and independent variables.These plots show the relative frequency of non-LIP (0) and LIP (1) pixels within the predictor classes.A higher frequency of LIPs occurs in steeper slopes and where the upslope and downslope curvatures are convex and concave, respectively.The N and E plots show that landslide density was higher on the The conditional density plots of Figure 9 reveal empirical relationships between landslide inventories and independent variables.These plots show the relative frequency of non-LIP (0) and LIP (1) pixels within the predictor classes.A higher frequency of LIPs occurs in steeper slopes and where the upslope and downslope curvatures are convex and concave, respectively.The N and E plots show that landslide density was higher on the north-and east-facing slopes of the datasets JF and FSE, respectively.Landslides are very rare at high-elevation JF positions, while they are more numerous where the elevation is high in the JAN dataset.As expected, the density of February LIPs decreases with the distance from the earthquake epicenter; surprisingly, the opposite occurs for the LIPs triggered by the January seismic event.A relatively high landslide frequency is observed in hard rocks (JAN), very soft soils (FEB and JF) and hard soils (FSE), as well as in upper (positive TPI) and lower (negative TPI) slopes.The density of LIPs is slightly higher than the average in intermediate PGA classes of the datasets FEB, JF and FSE.

The Ida/96E 2009 Inventory
The SV inventory [49] covers an area of 287 km 2 around the flanks of the San Vicente Volcano and includes the LIPs of 5609 landslides triggered between 7 and 8 November 2009 [85] by an extreme rainfall event.This event was caused by the combined action of Hurricane Ida and the low-pressure system 96E, which produced 350 mm in 24 h in an area of about 400 km 2 located between Ilopango Lake and the San Vicente Volcano (Figure 8).The rainfall triggered numerous and large flows that caused fatalities and massive damage, with an economic loss of approximately a quarter of USD one billion [86].

Predictive Performance of the MARS Models
The inventories and variables used for calibrating and validating the M1, M2, M3 and M4 models are shown on Table 5.The predictive skill of the models, measured on the validation samples of the datasets JAN, FEB, JF and FSE, are revealed by the boxplots in Figure 10.These reflect, for each model and dataset, the variability and skewness of the ten AUC values calculated by applying the calibration and validation strategies described in Section 2.6.The skewness displays, for each group of ten AUCs, their first quartile (Q1), median and third quartiles (Q3), the minimum and maximum values within 1.5 times the interquartile range (IQR) above Q3 and below Q1 and eventually outliers outside 1.5 of the IQR.Moreover, Table 8 shows the average and standard deviation values of the AUC groups.Furthermore, the average ROC curves calculated for each group of ten model runs are shown in Figure 11.The M1 models, calibrated with all rainfall-induced landslides that occurred near the San Vicente volcano and then transferred to the study areas, exhibit average AUC values ranging from 0.65 to 0.81.Overall, the accuracy of these models is poorer than that of the The M1 models, calibrated with all rainfall-induced landslides that occurred near the San Vicente volcano and then transferred to the study areas, exhibit average AUC values ranging from 0.65 to 0.81.Overall, the accuracy of these models is poorer than that of the other models, except for what was observed in the datasets JAN and JF, where the M1 models perform similarly to the M2 models or even better when validation is carried out with 95% of the landslides or by applying the spatial CV.The average performance of the M1 models is not acceptable for the datasets FEB and FSE, whereas it is almost acceptable in the dataset JAN [84].The M1 models show the best performance in the dataset JF, achieving a mean AUC in the range of 0.79-0.81.
The M2 models, which include as predictors the variables PGA and DIST, achieved mean AUC values in the range of 0.66-0.88.Their performance is better than the performance exhibited by the M1 models except for the dataset JAN when the validation scenarios 2 (test rate = 95%) and 3 (spatial CV) are applied.In these cases, as well as in the datasets FEB and FSE (with spatial CV), the predictive skills of the M2 models are not acceptable.On the other hand, their ability to discriminate between stable and unstable cells ranges from acceptable to excellent for the datasets FSE (with Scenarios 1 and 2) and JF.
The validation of the M3 models revealed mean AUC values in the range of 0.69-0.90,significantly higher than those achieved by the M1 and M2 models in all datasets except when applying spatial CV to the JF dataset.The average AUC achieved by the M3 models, which were fitted using the covariates PGA, DIST and PSV, ranges from acceptable to excellent, with lower values revealed by the spatial CV in the datasets FEB, JF and FSE.
The M4 models were trained with all independent variables available, excluding PSV.The validation results reveal a further improvement in the predictive ability of the M4 models compared to the predictive ability exhibited by the other models.The average AUC values, which are in the range of 0.73-0.96,exceed the acceptability threshold in the datasets JAN, FEB and FSE, achieving excellent to outstanding levels when the models are applied to the JF datasets.
The best predictive ability of the four tested models is observed in the dataset JF (mean AUC: 0.79-0.96),which contains the landslides triggered in sectors 9 to 14 by the January earthquake.Conversely, the models exhibited some inaccuracy in discriminating between stable and unstable cells of the JAN and FEB datasets, showing mean AUC values in the ranges of 0.67-0.79and 0.65-0.78,respectively.
Overall, when the learning samples include 75% of the observed landslides, models M2 to M4 exhibit a significantly better accuracy (mean AUC: 0.69-0.96)than the accuracy achieved when they are calibrated using only 5% of the landslides (mean AUC: 0.68-0.88)or when applying the 5-fold spatial CV (mean AUC: 0.66-0.90).However, except for the dataset JF, the differences in mean AUC values are small and insignificant when the M2 models are fitted with 5% of the FEB landslides.Moreover, regarding the M1 models, there is no significant difference in accuracy when validation is carried out with 25% and 95% of the landslides.As expected, larger differences in AUCs are observed when comparing validation Scenarios 1 (test rate = 25%) and 3 (spatial CV), since the latter involves the transfer of the models between different areas.
Validation samples containing 95% of the landslides reveal a more stable predictive skill (SD: 0.003-0.018)than the predictive skill obtained when validating the models with the 25% of the event cells (SD: 0.007-0.042).The performance of the models is also stable when the spatial CV is carried out in JAN, FEB and FSE (SD: 0.003-0.019),but the AUC values are rather dispersed in JF (SD: 0.025-0.068),where Scenarios 1 (SD: 0.018-0.042)and 2 (SD: 0.014-0.018)also exhibit variability in their performances.
The ROC curves of Figure 11 confirm that predictive skill increases from the models M1 to M4.Since the FEB dataset contains the largest number of LIPs (5371), the ROC curves prepared for this dataset are smoother.Conversely, those calculated for JF are scattered, as they are derived from samples including a few landslides (123).

Comparision of the Landslide Susceptibility Maps
Figure 12 shows an example of the landslide susceptibility maps generated by models M1, M2, M3 and M4 for two zoom areas (Figure 5).and 2 (SD: 0.014-0.018)also exhibit variability in their performances.
The ROC curves of Figure 11 confirm that predictive skill increases from the M1 to M4.Since the FEB dataset contains the largest number of LIPs (5371), th curves prepared for this dataset are smoother.Conversely, those calculated for JF a tered, as they are derived from samples including a few landslides (123).

Comparision of the Landslide Susceptibility Maps
Figure 12 shows an example of the landslide susceptibility maps generated b els M1, M2, M3 and M4 for two zoom areas (Figure 5).As they are landslide-prone areas, MARS score values tend to be high.Overal generated by the M2 models show abrupt changes related to the pattern of PGA while the other maps are smoother except for the M3 and M4 maps of the dataset ures 13 and 14 reveal the values of the true positive rate (TPR), true negative rate false positive rate (FPR) and false negative rate (FNR) calculated for the zoom a Figure 12 and for the entire landslide susceptibility maps, respectively.These metri calculated by setting a probability cut-off value of 0.5.In the zoom areas, the maps a very high rate of true positive but also a very high FPR (Figure 13).This is parti As they are landslide-prone areas, MARS score values tend to be high.Overall, maps generated by the M2 models show abrupt changes related to the pattern of PGA values, while the other maps are smoother except for the M3 and M4 maps of the dataset JF.Figures 13 and 14 reveal the values of the true positive rate (TPR), true negative rate (TNR), false positive rate (FPR) and false negative rate (FNR) calculated for the zoom areas of Figure 12 and for the entire landslide susceptibility maps, respectively.These metrics were calculated by setting a probability cut-off value of 0.5.In the zoom areas, the maps exhibit a very high rate of true positive but also a very high FPR (Figure 13).This is particularly evident for the M2 and M3 models.On the other hand, the entire susceptibility maps show a better discrimination ability, with markedly lower values of FPR, even if a decrease in true positives is also observed (Figure 14).These histograms confirm that the best ability to discriminate landslide pixels is achieved for the dataset JF by the models M2 (TPR: 0.886; TNR: 0.750), M3 (TPR: 0.894; TNR: 0.797) and M4 (TPR: 0.976; TNR: 0.856).
evident for the M2 and M3 models.On the other hand, the entire susceptibility maps s a better discrimination ability, with markedly lower values of FPR, even if a decrea true positives is also observed (Figure 14).These histograms confirm that the best a to discriminate landslide pixels is achieved for the dataset JF by the models M2 ( 0.886; TNR: 0.750), M3 (TPR: 0.894; TNR: 0.797) and M4 (TPR: 0.976; TNR: 0.856).evident for the M2 and M3 models.On the other hand, the entire susceptibility maps s a better discrimination ability, with markedly lower values of FPR, even if a decrea true positives is also observed (Figure 14).These histograms confirm that the best ab to discriminate landslide pixels is achieved for the dataset JF by the models M2 ( 0.886; TNR: 0.750), M3 (TPR: 0.894; TNR: 0.797) and M4 (TPR: 0.976; TNR: 0.856).

Discussion
In January and February of 2001, El Salvador was hit by two intense seismic phenomena of different natures which produced thousands of gravitational phenomena.Both 2001 earthquakes triggered similar landslides, but their distribution varied because of different earthquake source parameters, such as magnitude, epicentral distance and hypocentral depth.Indeed, earthquakes with large hypocentral depths tend to produce a relatively low number of gravitational phenomena spread over large areas.In contrast, phenomena that occur at a relatively low hypocentral depth usually produce many landslides concentrated in areas close to the epicenter [10].The seismic event of January, with a hypocentral depth of 60 km, triggered 1200 gravitational phenomena across 398 km 2 , whereas 5371 landslides, spread over 305 km 2 , were induced by the February event, which had a hypocentral depth of 10 km.The difference in hypocentral depth caused different spatial distributions of the PGA and MMI (modified Mercalli intensity).In fact, their values decreased from the hypocenter to a greater distance in the January event compared to the February event.The slow decay of PGA values, which is observed for earthquakes with a deep and offshore genesis, triggers gravitational phenomena over larger areas than earthquakes with the most superficial genesis.Conversely, the latter are characterized by a higher concentration of slope failures around the hypocenter.
In this work, four different MARS models (M1 to M4) were produced using different predictor sets and landslide inventories, the latter containing the slope failures triggered by the Ida/96E event of November 2009 (M1) and those induced by the earthquakes of January and February 2001 (M2 to M4).Two different calibration and validation scenarios were employed.In the first scenario, three quarters of the observed landslides were used for calibrating the models M2 to M4, whereas the remaining quarter was included in the validation samples.This scenario is frequently used to produce and validate landslide susceptibility models [65,80,87,88].In the second scenario, the learning samples contained only 5% of the earthquake-induced landslides, using the remaining 95% of the LIPs for validation.Additionally, to further investigate the predictive skill of the models, a third validation approach based on a 5-fold spatial cross-validation (CV) was applied.
The AUC values calculated for each model, dataset, and validation scenario, which are summarized by Table 8 and the boxplots of Figure 10, reveal that the MARS algorithm provides reliable predictions of the spatial distribution of earthquake-induced landslides.Indeed, 73% of the 320 AUC values AUC that were computed by applying the two test rates (25 and 95%) exceed the threshold of acceptability (AUC > 0.7), and a quarter of them even reveal the excellent predictive skill of the models (AUC > 0.8).Regarding the spatial CV, these percentages decrease to 50% and 21%, respectively.This is related to the data-driven nature of the adopted approach, which may generate a very good ability to predict unknown events if calibration and validation samples are extracted from the same area, even when time partition is performed [61]; however, models tend to lose predictive skill when transferred between different areas.As confirmation of this, the models trained with only 5% of the event pixels but which were spatially randomly distributed over all the sectors of the datasets achieve, in general, a better performance, with 68% of the model repetitions achieving an AUC value > 0.7.
When comparing the performance of the models M2 to M4, which were assessed with the two different percentages of validation LIPs (i.e., 25 and 95%), it is worth noting that the mean AUC gap is never greater than 0.08 and is0.04 on average.These results suggest that reliable predictive models of earthquake-triggered landslides can be obtained even when only a small portion of the slope failures have been identified.Therefore, in this situation, which is typical of the first moments immediately after a triggering event, predictive models such as those used in this work can help to quickly identify slopes where unmapped mass movements may have occurred.Indeed, even after over twenty years, only one fifth of the landslides identified in this work were already known.The latter were the main landslides that had hit the inhabited centers and which had been mapped immediately after the events to aid civil protection operations.The difference in the AUC values between the first validation scenario (test rate = 25%) and the 5-fold spatial CV are slightly higher (max: 0.10; mean: 0.05), which is acceptable if we consider that spatial transfer is very challenging for statistical landslide susceptibility zonation.
A general growing trend of predictive skill from M1 to M4 can be observed except in the datasets JAN and JF, where the predictive performance of the M1 and M2 models is similar.Conversely, when landslides triggered by the February earthquake are considered (i.e., datasets FEB and FSE), the M2 models' predictive ability is clearly better than the ability achieved by the M1 models, which do not achieve acceptable accuracy.This could be explained by a weaker control of the earthquake variables (i.e., PGA and DIST) on the spatial distribution of the January landslides due to a higher hypocentral depth and the largest distance between the epicenter and observed slope failures with respect to the February events.This hypothesis is confirmed by Figure 9, which reveals that the frequency of January LIPs increases with the distance from the earthquake epicenter and shows no relationship with the PGA.On the other hand, the higher accuracy of the M3 models compared to the M2 models, which is observed in all four datasets, reveal that a rainfall-induced landslide susceptibility model, even if calibrated outside the validation area, can significantly improve the predictive ability of a model containing only trigger variables.However, the fact that the M4 models achieve the best accuracy highlights the need to include preparatory variables measured in the same area in which validation is performed.
The AUC values also show that the models' predictive skill observed in the dataset FSE is better than that achieved in the dataset FEB.Based on these results, we could infer that rainfall-triggered landslide susceptibility models help more in predicting earthquakeinduced slope failures when focusing on the same area in which the models were calibrated.
The landslide susceptibility maps of Figure 12 appear smoother as the number of predictors and/or training LIPs increases.Indeed, models with only two predictors (i.e., M2) and/or those that are fitted with few landslide pixels (i.e., M2, M3 and M4 calibrated for JF) produce more heterogenous maps.On the other hand, the M1 models, which include ten predictors and were calibrated with 5609 LIPs, and the M3 (three predictors) and M4 (twelve predictors) models fitted on JAN (997 LIPs), FEB (5371 LIPs), JF (123 LIPs) and FSE (1587 LIPs) generate smoother landslide susceptibility maps.As suggested by Steger et al. [56], landslide susceptibility maps with abrupt changes, such as those generated for JF by the M2, M3 and M4 models, are characterized by poor geomorphic plausibility and applicability and may be influenced by biased relationships between independent and dependent variables.Conversely, the other landslide susceptibility maps, which exhibit smoother prediction patterns, are more easily readable and thus might be more useful for land planning purposes, even though they derive from models with a slightly lower predictive skill.
The high number of false positives observed in the zoom areas of the landslide susceptibility maps (Figure 13) can be reasonably explained by considering that the models were fitted with a 1:1 sampling ratio of positives to negatives and that in these areas, even though they are among the most prone to landsliding, event pixels are very rare compared to the number of non-event pixels.Further research could be carried out to explore if increasing the ratio of negative to positive cases or changing the cut-off value could decrease the FPR while keeping the value of the TPR high.However, when considering entire maps (Figure 14), the number of false positives decreases significantly, although the TPR remains relatively high in most cases.The FPR and FNR values calculated in our study are comparable with those found in other studies that applied empirical models to model the spatial distribution of earthquake-induced landslides.Nowicki Jessee et al. [28], who developed a global empirical model for the near-real-time assessment of seismically induced landslides, found FPR and FNR values in the ranges of 0.06-0.45and 0.02-0.47,respectively.These values changed in relation to the predicted probability, with higher error rates calculated for the average probability classes.In our study, FPR values are in the ranges 0.31 (M1)-0.37(M2), 0.30 (M4)-0.58(M1), 0.14 (M4)-0.63(M1) and 0.27 (M4)-0.42(M1), for the datasets JAN, FEB, JF and FSE, respectively.Regarding the FNR, the observed values are markedly lower, being in the ranges 0.26 (M4)-0.38 (M1), 0.20 (M1)-0.36(M2), 0.02 (M4)-0.11(M2) and 0.25 (M4)-0.29 (M1).Some limitations must be considered when interpreting the results of our study.First, due to the lack of area photos covering larger areas, it was not possible to evaluate the presence/absence of landslides in a wider range of epicentral distances and PGA values.This would probably have led to an increase in the weight of these variables in the spatial distribution of landslides, thus determining a better predictive capacity of the models.Second, due to the resolution and spectral limits of the aerial photos that were visually analyzed in this study, we could not identify and consider permanent deformations of the slopes and fractures with little or no displacement, which are also important effects of an earthquake and may lead to severe infrastructural damage [89,90].Unfortunately, the lack of aerial photos could not be remedied by using satellite images as their spatial resolution is lower.Furthermore, due to the coarser resolution of the SPOT-4 satellite images with respect to the resolution of the aerial photos employed to map the slope failures of the other datasets, the JF dataset includes only landslides extending more than 100 m 2 .This limitation probably did not allow for the identification of smaller landslides and therefore resulted in a significantly lower number of LIPs in the JF dataset than in the other datasets.

Concluding Remarks
Based on the results of this experiment, the following main conclusions can be drawn for the areas of El Salvador (CA) that were most affected by the earthquakes of 2001.

•
Overall, the MARS algorithm provides reliable predictions of the spatial distribution of earthquake-induced landslides.

•
MARS models calibrated with rainfall-induced landslides (i.e., M1 models), even if exported from other areas, are able to predict landslides distributed over vast areas and caused by deep earthquakes, such as the one in January 2001, with an acceptable accuracy.

•
When MARS models are based only on seismic covariates (i.e., M2 models), including as a predictor the probability of occurrence of rainfall-induced landslides (i.e., M3 models), even when transferred from other areas, a significant increase inaccuracy can be observed.

•
The models containing both preparatory and trigger variables (i.e., M4 models) achieve the best accuracy in predicting the spatial distribution of earthquake-induced landslides.

•
Even when only a small part of the landslides produced by an earthquake is known, as is usually the case shortly after the event, it is possible to use them to calibrate models, such as those used in this work, that help to identify slopes where yet unreported landslides may have occurred.
The results of this experiment and the main conclusions listed above suggest that the approach presented in this work can be useful for predicting the distribution of landslides caused by earthquakes such as those that occurred on Salvadoran territory in January and February of 2001.Seismic events with similar characteristics occur quite frequently in El Salvador due to the continuous subduction of the Cocos Plate beneath the Caribbean Plate, which generates severe offshore earthquakes, resulting in intraplate phenomena in a relatively short time.Therefore, it is useful to provide the country with landslide susceptibility maps that are also calibrated with slope failures induced by extreme rain events, which occur with greater frequency than earthquake-induced landslides.Indeed, as observed in our experiment, these maps can be used as a basis to identify slopes in which a future seismic event is most likely to cause landslides.On the other hand, in the case that a high magnitude earthquake occurs and the positions of some landslides caused by the same event are known, the models based on predisposing variables can be integrated with trigger variables such as PGAs and epicentral distances in order to obtain predictive maps that may help local civil protection to identify sites potentially affected by others and unknown earthquake-induced landslides.

Figure 1 .
Figure 1.Location and geological setting of the study area.El Salvador is located in Central America where the Cocos Plate is subducted beneath the Caribbean Plate into the Middle American Trench (a); major cities and geological setting of El Salvador, mainly characterized by the volcano chain and by the large outcropping of volcanic rocks (b).

Figure 2 .
Figure 2. Location of the fourteen study areas: sectors 1-3 and 4-8 were heavily affected by the earthquake of 13 January 2001; sectors 9-14 were severely affected by the earthquake of 13 February 2001.

Figure 1 .
Figure 1.Location and geological setting of the study area.El Salvador is located in Central America, where the Cocos Plate is subducted beneath the Caribbean Plate into the Middle American Trench (a); major cities and geological setting of El Salvador, mainly characterized by the volcano chain and by the large outcropping of volcanic rocks (b).

Figure 1 .
Figure 1.Location and geological setting of the study area.El Salvador is located in Central Americ where the Cocos Plate is subducted beneath the Caribbean Plate into the Middle American Tren (a); major cities and geological setting of El Salvador, mainly characterized by the volcano chain an by the large outcropping of volcanic rocks (b).

Figure 2 .
Figure 2. Location of the fourteen study areas: sectors 1-3 and 4-8 were heavily affected by t earthquake of 13 January 2001; sectors 9-14 were severely affected by the earthquake of 13 Februa 2001.

Figure 2 .
Figure 2. Location of the fourteen study areas: sectors 1-3 and 4-8 were heavily affected by the earthquake of 13 January 2001; sectors 9-14 were severely affected by the earthquake of 13 February 2001.

Figure 3 .
Figure 3. Main earthquakes (Mw > 5.5) that affected El Salvador.The earthquakes that hit the country on 13 January and 13 February 2001 are highlighted.

Figure 3 .
Figure 3. Main earthquakes (Mw > 5.5) that affected El Salvador.The earthquakes that hit the country on 13 January and 13 February 2001 are highlighted.

Figure 4 .
Figure 4. Positions of the epicenters, paths and areas covered by the aerial images employed to map the landslides that occurred in January and February 2001.

Figure 4 .
Figure 4. Positions of the epicenters, paths and areas covered by the aerial images employed to map the landslides that occurred in January and February 2001.

Figure 5 .
Figure 5. Locations of the datasets used for the landslide susceptibility modeling.

Figure 5 .
Figure 5. Locations of the datasets used for the landslide susceptibility modeling.

Figure 6 .
Figure 6.PGA values and LIPs related to the seismic event of Mw 7.7 on 13 January 2001 in the fourteen study areas.

Figure 6 .
Figure 6.PGA values and LIPs related to the seismic event of Mw 7.7 on 13 January 2001 in the fourteen study areas.

Figure 7 .
Figure 7. PGA values of the seismic event of 13 February 2001 and LIPs of the landslides triggered in the FEB study areas.

Figure 7 .
Figure 7. PGA values of the seismic event of 13 February 2001 and LIPs of the landslides triggered in the FEB study areas.

Figure 8 .
Figure 8. LIPs of the landslides triggered by the seismic event on 13 February 2001 in the intersection between FEB and SV study areas (FSE).

Figure 8 .
Figure 8. LIPs of the landslides triggered by the seismic event on 13 February 2001 in the intersection between FEB and SV study areas (FSE).

Figure 12 .
Figure 12.Examples of the landslide susceptibility maps extracted from the zoom areas 1 JAN) and 2 (datasets FEB, JF and FSE).

FSEFigure 12 .
Figure 12.Examples of the landslide susceptibility maps extracted from the zoom areas 1 (dataset JAN) and 2 (datasets FEB, JF and FSE).

Figure 14 .Figure 13 .
Figure 14.TPR, TNR, FPR and FNR values calculated for the entire landslide susceptibility ma the datasets JAN, FEB, JF and FSE.

Figure 13 .
Figure 13.TPR, TNR, FPR and FNR values calculated for the landslide susceptibility maps i zoom areas 1 (dataset JAN) and 2 (datasets FEB, JF and FSE).

Figure 14 .Figure 14 .
Figure 14.TPR, TNR, FPR and FNR values calculated for the entire landslide susceptibility ma the datasets JAN, FEB, JF and FSE.

Table 2 .
Precipitation in the 12 months preceding the seismic events of 2001 and the average annual precipitation in 1971-1999 in the fourteen areas under consideration.

Table 3 .
Independent variables used to model rainfall-and earthquake-induced landslides.

Table 4 .
Datasets used to calibrate and validate the landslide predictive models.
SV San Vicente area (Mercurio et al., 2021 [49]) Rainfall November 2009 JAN Sectors 1 to 8 Earthquake January 2001 FEB Sectors 9 to 14 Earthquake February 2001 JF Sectors 9 to 14 Earthquake January 2001 FSE Intersections between FEB and SV areas Earthquake February 2001

Table 4 .
Datasets used to calibrate and validate the landslide predictive models.

Table 6 .
Amounts and characteristics of the landslides triggered by the January earthquake in sectors 1 to 14.

Table 7 .
Amounts and characteristics of the gravitational phenomena of landslides belonging to sectors 9 to 14, related to the February earthquake.