As input data, we primarily used the global EQILs catalogue managed by the US Geological Survey (USGS), which serves as a centralized repository that adopts a standardized format (Schmitt et al., 2017; Tanyas et al., 2017). We considered polygonal inventories and excluded datasets with fewer than 50 mapped landslides, as we observed that these often show limited spatial coverage and variability, thus compromising the robustness of method. By applying this threshold, we obtained a total of 36 polygonal inventories from the USGS catalogue. Additionally, we incorporated 4 inventories mapped by ourselves (Ferrario, 2022; Ferrario et al., 2024), bringing the total to 40 inventories, which refers to 33 earthquakes. For 6 of these earthquakes, multiple inventories were available; we utilized these case histories to evaluate the different sources of uncertainties in our workflow.

Table 1 presents the main characteristics of the inventories, listed chronologically from oldest to most recent. We considered the polygons regardless of the type of movement; this information could be available for some case studies, but the ESI-07 categorization is based on landslide volume only. Figure 1 illustrates the geographical distribution. Most of the inventories are from tropical regions, either in Asia or Latin America, although the dataset also contains a few inventories outside the tropics. This distribution is significant, as landslides in tropical regions remove dense vegetation, making remote detection more efficient. In contrast, in continental, arid, or anthropized areas, sparser vegetation and smaller landslides make satellite-based mapping more challenging. The number of earthquakes and associated landslides, categorized according to magnitude class or earthquake rupture kinematics, is presented in Figure 1. Moment magnitude range between Mw 5.3 and 9.1, with the majority of the events falling within the Mw 6.0 – 8.0 range. Regarding earthquake rupture kinematics, 50% of the selected inventories correspond to thrust earthquakes, while 42% are associated with strike-slip earthquakes.

Summing up the number of landslides for each inventory, a total of 430,996 is reached. Of these, 197,481 belongs to the Wenchuan inventory compiled by Xu et al. (2014a), and 69,606 landslides are from the inventory of Li et al. (2014), together accounting for 62% of the total. The selected inventories vary significantly in size, spanning several orders of magnitude in terms of landslide number. The smallest inventory includes 52 landslides from the Cappellades, Costa Rica earthquake (Ruiz et al., 2020), while the largest comprises a nearly 200,000 landslides from the Wenchuan, China earthquake (Xu et al., 2014a).

Figure 2 illustrates the methodological workflow adopted in this study. It comprises five sequential steps (A-E), starting with data collection (step A), described in the previous section. Steps B-E are described in the following subsections: each step produces an output that can serve as input for further analyses.

The core of our analysis is Step C of the workflow (Figure 2), which involves the calculation of LAP, LND and ESI-07 values for each 1-km² grid cell.

Figure 3 presents the obtained grid maps for the 2009 Cinchona (Costa Rica, ID 16 in Table 1) earthquake. In this case, the investigated area encompasses 502 cells, and the maps display the spatial distribution of ESI-07, LAP and LND values. The highest values are concentrated in the center of the area and the distribution is fairly similar, with lower values observed toward the periphery of the investigated area. The plots in Figures 3D, E depict the relative frequency of LAP and LND values, categorized according to ESI-07 classes. For lower intensities (ESI-07 VI and VII) LAP and LND exhibit low to medium values; ESI-07 VIII covers a broader range of LAP and LND values. Intensities IX and ≥ X show a different pattern: LAP is constrained to high values only (LAP >4%), whereas LND encompasses nearly all the classes (i.e., about 20% of the cells with LND = 2 are categorized as ESI-07 ≥ X). A comparable spatial distribution of LAP, LND and ESI-07 intensities was observed in all the inventories analyzed, confirming the general trend observed in the inventory of Cinchona earthquake. In the Supplementary Material, we provide the plots for other case histories.

To evaluate the consistency and repeatability of our procedure, we analyze each inventory individually. Our working hypothesis was that ESI-07 values should be positively correlated with either LAP or LND values. We consider this hypothesis to be reasonable, however it has never been tested systematically. We computed the median LAP or LND value for each ESI-07 degree. Table 4 presents the data for the 40 inventories, while the plots of Figure 4 include 33 case histories. When multiple inventories were available for a single earthquake, we retained only the most complete dataset (i.e., with the highest landslide count and largest investigated area), to avoid double counting.

All the analyzed inventories provide consistent trends: as LAP or LND increase, so does the ESI-07 intensity. LAP values (Figure 4A) show a steady increase over the entire range of intensities, moving from 0.06% at ESI-07 VI to LAP 24% at ESI-07 ≥ X. Additionally, a low dispersion is observed for ESI-07 VI to VIII, while a higher dispersion is observed for ESI-07 IX and ≥ X.

Regarding LND (Figure 4B), a positive correlation is observed for ESI-07 VI to IX, with the highest density at ESI-07 IX (median LND 8 landslides/km²), after which it decreases to 4 for ESI-07 ≥ X. The dispersion is greatest at intensity IX, with remarkable differences among individual inventories, especially for the Gansu earthquake (ID 19 in Table 1). The lower LND values for intensity ≥ X likely reflect a “physical” limitation: the largest landslides (volume > 10⁶ m³) correspond to an area of approximately 130,000 m² according to our baseline A-V relation. This constraint limits the number of such large landslides that can fit within each cell of the 1 km² grid, thus fundamentally limiting the resulting LND values.

By analyzing the case histories, we obtained consistent results over a wide range of magnitudes and territorial settings. We then proceed to Step D of the workflow, deriving empirical regressions linking ESI-07 values to either LAP or LND.

We considered the median of the 33 earthquakes and tested various functional forms. Following the observation that LND values increase in the ESI-07 range VI to IX, while ESI-07 ≥ X has a lower LND value, we decided to fit the data over the range of intensities VI to IX. For LAP values, data are fitted over the entire range of intensities, i.e., ESI-07 VI to ≥ X. Both power-law and logarithmic forms consistently gave the best fitting results; for LAP values, power-law model performed slightly better than the logarithmic law. We selected the power-law functional form for both LAP and LND, because we want to prioritize the LAP dataset, which we deem as better suitable to be used to obtain ESI-07 values; the equations have the general formula: E S I 07 = a ·   x b (2)

Where x stands for LAP or LND, a and b are fitting coefficients.

Figure 5 and Table 5 summarize the obtained curves and fitting coefficients. We consider LAP as a better descriptor: indeed, LND depends more strongly on the quality of input data and is a more local metric, being more influenced by the topography, geo-mechanical characteristics and mapping style of the user who generated the inventory. One common issue in landslide inventories is the amalgamation of multiple landslides into a single polygon (Marc and Hovius, 2015); which can significantly affect the computation of LND, but does not affect LAP values. Nevertheless, in the case of point inventories, LND is the only suitable option. Therefore, we consider it appropriate to provide empirical equations for LND as well.

The number and dimension of triggered landslides are influenced by several factors, related to terrain (e.g., elevation, slope), as well as climatic, geological and seismological parameters. A plethora of studies have analyzed susceptibility factors, even focusing on LAP or LND metrics (e.g., Chang et al., 2021; Dai et al., 2023; Shao et al., 2023). Here we limit our analysis to seismological and tectonic factors., specifically moment magnitude and earthquake kinematics, to evaluate the possibility of deriving empirical regressions dependent on these two factors. We grouped the landslide inventories according to either magnitude or kinematics, looking for eventual differences in the obtained results. One reasonable hypothesis that we may test is whether a higher magnitude, which corresponds to a greater energy release, leads to higher LAP and LND values. Similarly, different earthquake kinematics may vary in their effectiveness in triggering landslides. Fault geometry, the type of fault movement and the presence or absence of surface faulting indeed all influence the distribution and dimension of the area affected by coseismic landslides (Gorum et al., 2014; Xu, 2014).

The boxplots of Figure 6 present the distribution of LAP and LND values according to moment magnitude, divided into 3 classes. As already observed, LAP values exhibit a much narrower distribution than LND. No clear dependence on magnitude is found for LAP values, as the 3 magnitude classes considered here have very similar median LAP values and overall distribution. LND showed lower values for stronger earthquakes (Mw class 7.0–7.9) with respect to weaker earthquakes. The result is not unexpected when considering the general scenario: stronger earthquakes tend to mobilize substantial masses of material, resulting in fewer but more extensive landslides, a key factor that reduces the overall LND value.

Figure 7 presents the distribution according to earthquake kinematics. A robust evaluation can be conducted for thrust and strike-slip earthquakes (15 case histories each), while normal earthquakes are not adequately represented (2 cases only). LAP again provide more informative results than LND. Thrust events tend to have slightly lower LAP values, especially for intensity classes IX and ≥ X.

Overall, it is evident that magnitude and kinematics exert a limited influence on LAP and LND values. We notice a clear step in LAP values, especially from ESI-07 intensity VII to intensity ≥ X, indicating that intensity classes are characterized by increasing LAP values. This confirms that LAP provides the most robust relationship.

The ESI-07 intensity value assigned to each landslide fundamentally depends on two key factors: the area of the polygon and the choice of an area-volume scaling relation. Both elements represent sources of epistemic uncertainty, as they are dependent on the mapping procedure and the methodological choices made by the users. We investigate these aspects taking advantage of the six earthquakes where multiple landslide inventories are available.

The inventories may have been realized by different research groups (e.g., Harp and Jibson (1995); Townsend et al. (2020) for the 1994 Northridge earthquake) or represent successive refinements of an initial inventory (e.g., Basharat et al., 2014; 2016 for the 2005 Kashmir earthquake). In some cases, the number of mapped landslides and the dimensions of the investigated area are similar (e.g., Sepulveda et al. (2010) and Gorum et al. (2017) for the 2007 Aisen earthquake), while in other cases they are sensibly different (e.g., Gorum et al. (2013); Harp et al. (2016) for the 2010 Haiti earthquake). The inventories are produced with different methods (manual vs. semi-automatic mapping) and using pre- and post-event images with variable resolution.

We applied five different A-V relations (Table 2), selected based on their derivation from global databases (Larsen et al., 2010; Guzzetti et al., 2009), their focus on earthquake-triggered landslides or their relevance to study areas examined in this work (Xu et al., 2016). For each A-V relation, we calculate the frequency of landslide polygons belonging to each ESI-07 intensity degree. Data are presented in Figure 8: red lines refer to the most complete inventory (labeled with “a” in Table 1), while green and blue lines refer to inventories labeled as “b” or “c”, respectively.

Several considerations arise from Figure 8: the most striking observation is that input data (i.e., landslide inventories) are the first-order element in driving the ESI-07 assessment, while the selected A-V relation plays a secondary role. Indeed, in almost all the plots (with the exception of the Aisen and partly of the Kashmir events) the lines are clustered by color, representing the different inventories. The critical role exerted by different inventories on derived products was investigated by Bornaetxea et al. (2023) in a susceptibility study in Spain, where they demonstrated that altering the input landslide inventory leads to significantly different results.

The A-V relation by Xu et al. (2016) and Larsen et al. (2010), soil deviate from the other equations, especially for ESI-07 intensity VI. For the Northridge, Haiti and Nepal events a sensible difference in the number of small landslides (ESI-07 intensity VI) is observed between the inventories, pointing to a different degree of completeness.

Our results support the arguments made by Ferrario (2022) and Sridharan et al. (2023), who claim that the choice of a given A-V relation has limited implications on the final output. We underscore the importance of thoroughly delineating coseismic landslides and selecting the most appropriate A-V scaling relation, which is consistent with the local terrain and geologic conditions. We also recall that ESI-07 classes are broad in terms of volume (see Table 3), which helps smooth out overall variability. If the inventories are instead used for purposes other than intensity assessment, for instance to calculate the amount of mobilized sediments or denudation rates, a well-constrained A-V relation is a crucial input.

Landslides are not the only type of earthquake environmental effect and thus an evaluation based solely on this effect may be misleading. Recognizing this, our study aims to find a link to bridge the gap between two distinct yet complementary technical-scientific communities: a community focused on the construction of landslide inventories and the analysis of predisposing/triggering factors for slope movements, and another community dedicated to the individuation of EEEs and their intensity assessment using the ESI-07 scale. Both approaches ultimately aim to improve the general understanding of seismic hazard and to generate crucial information and instruments for mitigation purposes. By integrating these methodologies, a more comprehensive characterization of earthquake impacts could be achieved, accounting for the diverse environmental, geological and geomorphological settings where different EEEs occur. To evaluate the reliability and applicability of the ESI-07 methodology developed in this study, we compared our result with independent data, searching for those events that were already analyzed in terms of ESI-07 scale. A recent review of the earthquakes analyzed with the ESI-07 scale (Ferrario et al., 2022) included a dataset of more than 150 events with an available ESI-07 epicentral intensity (I0) estimate. A total of 15 earthquakes from this review overlap with our dataset, as shown in Table 6. The data were extracted from Ferrario et al. (2022), supplemented with the most recent studies of Sridharan et al. (2023); ID 25 in Table 1), Naik et al. (2023); ID 28 in Table 1, in this case isoseismals are from an unpublished thesis, see statements and declarations section) and Ferrario et al., 2024 (ID 29–31).

In the last column of Table 6 we report the highest ESI-07 value assigned in our study (i.e., deriving from the analysis of slope movements only) together with the percentage of cells with such value. We notice that the ESI-07 I0 is generally higher than the estimate derived from slope movements; this fact is not surprising, since i) intensities XI and XII cannot be assigned from slope movements and ii) widespread environmental effects other than slope movements may have been triggered. Additionally, we analyzed slope movements on emerged lands, while some of the earthquakes had their epicenter offshore (e.g., ID 18 Tohoku-Oki, which triggered a devastating tsunami). In a few cases (ID 10, 30 and 31) landslides returned to an ESI-07 value slightly greater than I0; nevertheless, the percentage of cells with such anomalous outcome is less than 0.5%. We remark that ESI-07 I0 is estimated based on the amount of permanent ground deformation (surface faulting and tectonic uplift/subsidence) or on the dimensions of the area affected by secondary effects (Michetti et al., 2007); thus, the analysis of slope movements alone provides no clue for ESI-07 I0 estimation.

We selected 4 earthquakes where reliable ESI-07 isoseismals are available and attempted to make a careful comparison of the macroseismic field proposed in the literature with the grid results obtained in our study. The selected events include Guatemala (ID 1), Denali (ID 8), Wenchuan (ID 14) and Palu (ID 28); therefore, we exclude events with only one or two ESI-07 isoseismals (e.g., Davao, ID 29–31) or an incomplete macroseismic field due to large offshore areas (e.g., Tohoku, ID 18).

We superimposed the ESI-07 isoseismals to the grids obtained in our study and we calculated the difference between the two for each grid element. It is worth mentioning that all the 4 investigated earthquakes generated extensive surface faulting, which resulted in ESI-07 isoseismals elongated in the direction of the seismogenic source (Figure 9). The histograms show the difference between ESI-07 from isoseismals and ESI-07 derived in this study from slope movements; positive values indicate that the estimate from isoseismals is higher. The histogram patterns for the four investigated earthquakes appear to be similar, with a frequency peak at 1 degree of difference, except for the Denali event, where the distribution is slightly shifted toward higher values. This analysis must be intended as a first approximation, since the territory covered by the isoseismals and the area investigated with landslides may differ: for instance, the Guatemala inventory covers only a small portion of the isoseismals (Figure 9B). Moreover, the grid analysis results in a “salt and pepper” pattern, with high variations between cells close to each other, while isoseismals by definition tend to smooth out such high-frequency variations and provide a more homogeneous description of earthquake effects. Finally, we recall that the maximum ESI-07 value obtainable from landslides is ≥ X; thus, all the cells with ESI-07 isoseismal XI and XII have a minimum “difference” value of 1 and 2, respectively. A total of 1815 cells lie within isoseismal XI, and 924 cells within isoseismal XII, constituting 21% of the analyzed grid elements; we believe that the overall distribution of the values is not heavily affected by the inability to fully compare the ESI-07 values from isoseismal lines and the analysis of slope movements.

Despite the above-mentioned limitations, here we document the suitability of this comparison and claim that additional insight could be gained from the analysis of a wider number of case histories, including earthquakes that did not generate surface faulting. Indeed, some studies have already pointed out a correlation between landslide concentration and distance from the seismic rupture (Chen et al., 2022; Zhao et al., 2023).

Our results support the notion proposed by Ota et al. (2009) and Silva et al. (2013), who first attempted to an ESI-07 evaluation via a box-counting approach. Other attempts were realized on the 2008 Wenchuan (China) earthquake and on a Mw 6.6 earthquake occurred in Luding County (China). Xu et al. (2013) provided a correlation among LND or LAP and intensity degrees on the Chinese macroseismic scale, while Yang et al. (2023) determined a so-called macro-epicenter based on the maximum value of landslide area density and defined numerical thresholds for a qualitative zoning of coseismic landslide intensity.