An open dataset for landslides triggered by the 2016 Mw 7.8 Kaikōura earthquake, New Zealand

On November 14, 2016, the northeastern South Island of New Zealand was hit by the magnitude Mw 7.8 Kaikōura earthquake, which is characterized by the most complex rupturing mechanism ever recorded. The widespread landslides triggered by the earthquake make this event a great case study to revisit our current knowledge of earthquake-triggered landslides in terms of factors controlling the spatial distribution of landslides and the rapid assessment of geographic areas affected by widespread landsliding. Although the spatial and size distributions of landslides have already been investigated in the literature, a polygon-based co-seismic landslide inventory with landslide size information is still not available as of June 2021. To address this issue and leverage this large landslide event, we mapped 14,233 landslides over a total area of approximately 14,000 km2. We also identified 101 landslide dams and shared them all via an open-access repository. We examined the spatial distribution of co-seismic landslides in relation to lithologic units and seismic and morphometric characteristics. We analyzed the size statistics of these landslides in a comparative manner, by using the five largest co-seismic landslide inventories ever mapped (i.e., Chi-Chi, Denali, Wenchuan, Haiti, and Gorkha). We compared our inventory with respect to these five ones to answer the question of whether the landslides triggered by the 2016 Kaikōura earthquake are less numerous and/or share size characteristics similar to those of other strong co-seismic landslide events. Our findings show that the spatial distribution of the Kaikōura landslide event is not significantly different from those belonging to other extreme landslide events, but the average landslide size generated by the Kaikōura earthquake is relatively larger compared to some other large earthquakes (i.e., Wenchuan and Gorkha).


Introduction
The Mw 7.8 Kaikōura earthquake struck the northeastern South Island of New Zealand at 12:03 local time (11:03 UTC) on November 14, 2016. It was one of the most complex earthquakes ever recorded (Cesca et al. 2017;Hamling et al. 2017). The earthquake hypocenter has been located at 42.69° S, 173.02° E, at a depth of 15 km in a well seismic instrumented region. The earthquake ruptured a sequence of active faults, terminating offshore and markedly affecting coastal and inland areas across the northeast of the South Island (Kaiser et al. 2017). The global moment tensor solutions revealed that this multi-fault rupture earthquake event is rather complex and that the overall rupture process propagated from south to north, connecting the region between the Hikurangi subduction system of the North Island and the oblique collisional regime of the South Island (Alpine Fault) (Cesca et al. 2017). The effects of the Kaikōura earthquake on the infrastructure and the environment were extremely severe and widely distributed, with a cost assessed in between 3 and 8 billion NZD in rebuilding efforts (Kaiser et al. 2017).
The preliminary investigations of landslides triggered by the 2016 Kaikōura earthquake were conducted by Dellow et al. (2017), Rathje et al. (2017), and Jibson et al. (2018). Following these investigations, Massey et al. (2018) introduced the first comprehensive landslide inventory of the Kaikōura earthquake. They mapped 10,195 landslides over a total area of approximately 10,000 km 2 . Massey et al. (2020) upgraded this inventory by mapping 29,557 landslides overall. Massey et al. (2018) compared the Kaikōura event with some previous earthquake-induced landslide (EQIL) events that also occurred in New Zealand. They showed that the number of landslides is less than what would have been expected in comparison to the 1929 Mw 7.8 Murchison earthquake but similar to the smaller magnitude 1968 Mw 7.1 Inangahua earthquake. They indicated that this could be because of the morphological constrains of the area affacted by the earthquke and/or fault rupturing mechanism including multiple segments.
Overall, comparative studies between different EQIL events are becoming increasingly common in the geomorphological literature (Tatard and Grasso 2013;Xu 2014;Fan et al. 2018) following the number of available landslide inventories, which has significantly increased during the last two decades (Tanyaş et al. 2017;Fan et al. 2019).
The 2016 Kaikōura earthquake offers an extensive landslide event where we can test our knowledge regarding the EQILs. For instance, Allstadt et al. (2018) used the Kaikōura event as a test case to evaluate the performance of the global seismically induced landslide models (Godt et al. 2008;Nowicki et al. 2014;Nowicki Jessee et al. 2018). They showed that all the models dramatically overpredicted the hazard level for this specific case. Other than sharing the polygon-based landslide inventory, this paper shares a similar motivation and standpoint of Allstadt et al. (2018). But, here we aim to test whether landslides triggered by the 2016 Kaikōura earthquake are less numerous and/or less large in planimetric extent than those triggered by other strong earthquakes or vice versa. Following this research question, first, we created a landslide inventory for the Kaikōura earthquake. Second, we examined the characteristics of landslides in terms of lithologic units and seismic and morphometric characteristics. And ultimately, we compared the Kaikōura inventory with respect to five other co-seismic landslide inventories associated with the 1999 Chi-Chi, 2002Denali, 2008Wenchuan, 2010Haiti, and 2015 Gorkha earthquakes. We did so by examining the landslide size statistics and the overall area affected by landsliding in response to ground motion.

Tectonic setting and earthquake characteristics of the study area
The active tectonics of New Zealand is driven by convergent and transform plate boundary processes. The Pacific Plate subducts beneath the Australian Plate (Hikurangi Subduction Zone) along the eastern coast of the northern Island, while the Australian Plate subducts beneath the Pacific Plate (Puysegur Subduction Zone) along the southwestern coast of the southern island. Geodetic data indicate that the rate of the motion between these two plates gradually decreases from the north (44 mm/year) to the south (32 mm/year) (see Nicol et al. 2017;Norris and Cooper, 2001). The relative motion of these two subductions is accommodated by the Alpine Transform Fault Zone (ATFZ). This fault zone is a rightlateral strike-slip fault zone that stretches along the western coast of the southern island. The ATFZ splays at the northeastern part of the southern island and constitutes ~ 80 km wide Marlborough Fault Zone (MFZ) sheared by NE-SW striking right-lateral strikeslip faults parallel to each other. The Hope, Clarence, Awatere, and Wairau faults are major splays from south to north (see Fig. 1). The Fig. 1 Tectonic setting and landslide distribution map of the study area regional strain is partitioned among these splays. The slip the rates of these faults are 20-25, 3-5, 4-6, and 3-5 mm/year, from the south to the north (Shi et al. 2017). The slip rates of major splay faults increase to the Hope Fault which is located just north of the epicenter of the Kaikōura earthquake. The slip rate of the Hope Fault has already ruptured during the 1888 (Mw 7-73) and 1780 (Mw 7.2) earthquakes. Paleoseismological data reveals that the recurrence interval of the Hope Fault is approximately 180 to 310 years (Langridge et al. 2003), but the recurrence intervals of the Clarence Fault, Awatere Fault, and Wairau faults are ~ 1700 years, ~ 820-950 y ears, and ~ 1150-1400 years, respectively (Van Dissen and Yeats 1991;Langridge et al. 2003;Mason et al. 2006;Zachariasen et al. 2006). The temporal distribution of recurrence intervals is consistent with the spatial distribution of slip rates. These data indicate that the highest regional strain has concentrated at the epicentral area of the Kaikōura earthquake. Although the NE-SW striking geometry of the major splays is evident, the presence of secondary structures with different strikes indicates an anastomosing fault network in the area. The geometry of the rupture of the Kaikōura earthquake seems to have a similar pattern. The Jordan and Kekerengu faults' surface ruptures stretch along with the northeastern continuation of the Hope Fault but thrust faults at the eastern end of the Emu Plain Humbs, Hundalee, and Papatea faults stretch roughly N-S (Litchfield et al. 2014;Langridge et al. 2016;Wang et al. 2018) and indicate complex oblique thrust nature of the mainshock (Ellis et al. 2017;Kaiser et al. 2017). GPS and InSAR (Interferometric Synthetic Aperture Radar) data reveal that most of the dip-slip movement occurred along the northeastern part of the rupture (Du et al. 2018). These kinematic differences give rise to the large heterogeneity of the slip along individual fault planes.

Landslide mapping
An extensive landslide interpretation was carried out using sets of optical 10 m resolution Sentinel-2 satellite images for both the pre-and post-earthquake conditions (Fig. 2). In particular, we used nine pre-seismic images acquired between 13th of September and 26th of October 2016 and nine post-seismic images between 22nd of November and 15th of December 2016. We primarily used images with cloud and shadow coverage of < 3% of the entire earthquakestruck region. The landslides were identified and mapped as polygons using multi-temporal visual image interpretation based on satellite imagery and morphological elements of landslide diagnostic indicators. Overall, 14,233 individual landslides were mapped over a total area of approximately 14,000 km 2 . We also used highresolution Google Earth images to confirm if the mapped landslides were triggered by the earthquake but not by the successive rainfall events. The confirmation was done after we create the inventory because the Google Earth images showing the conditions right after the earthquake (14th of November 2016) were not initially available.
Notably, 90% of these landslides clustered in an area of about 4000 km 2 , which could be defined by a ~ 15-km buffer zone around the rupturing fault lines (Fig. 1). Shape factor, drainage disruption, and the existence of lakes were used as diagnostic elements to identify landslides that have dammed rivers (see Fig. 1 for a complete overview and Fig. 2 for a detailed example). The distribution of 101 damming landslides was identified and mapped as a subcategory of the total landslide distribution.

Controlling factors' overview
We examined the co-seismic landslides in relation to lithologic units and seismic and morphometric characteristics (summarized in Table 1). Lithologic units are associated with the shear strength of hillslope materials, and therefore, we take them into consideration as factors possibly contributing to the slope failure mechanism. We accessed the lithological map (Lynn 1985) of the study via The Land Resource Information System Portal (https:// lris. scinfo. org. nz/) of New Zealand.
To analyze the effect of ground motion, we examined two different peak ground acceleration (PGA) maps of the Kaikōura earthquake provided by the U.S. Geological Survey (USGS) ShakeMap system (Worden and Wald 2016) and Bradley et al. (2017). PGA is a seismic proxy, and specifically, the deterministic estimate of PGA provided by the USGS ShakeMap system is widely used in susceptibility analyses of co-seismic landslides (e.g., Godt et al. 2008;Nowicki et al. 2014;Tanyaş et al. 2019a). We used the most updated version of the ShakeMap products (i.e., ShakeMap Atlas v4). As an alternative to this product, we also used the grid of the predicted mean PGA generated from broadband ground-motion simulations (Bradley et al. 2017).
Distance to rupture zone is another common seismic proxy, which contains a similar spatial signal to the PGA because, overall, the intensity of seismic shaking decreases away from the rupture zone. For instance, Massey et al. (2018) already examined both the ground motion proxies and the distance to the rupture zone as part of the multivariate analyses conducted for the landslides triggered by the Kaikōura earthquake. They conclude that distance to rupture zone works better in terms of explaining the landslide distribution of the inventory they mapped. In our analyses, we also examined the spatial distribution of landslides based on the distance to fault surface rupture. For the calculation, we used the fault surface rupture provided by the GNS Science Active Faults database (https:// data. gns. cri. nz/ af/; Litchfield et al. 2018). For each grid cell, we calculated the distance from the surface rupture to create a map representing the distance to the rupture zone. This procedure was done to ensure that the minimum distance was calculated orthogonally with respect to the rupturing tectonic line.
We also included aspect, slope, and local relief as morphometric characteristics, for they are indicators of the slope geometry and of some components of the driving forces acting over the hillslope in/stability equilibrium. We used TANDEM-X (a Ter-raSAR-X add-on for digital elevation measurements with 12-m spatial resolution) to calculate these DEM derivatives. To calculate local relief, which is the maximum difference in elevation within a defined neighborhood, we scanned the study area with a circle of 1 km radius. We then computed the aspect as an additional morphometric characteristic because the rupturing direction may contribute to amplify the ground motion on slopes facing towards the same direction of the co-seismic slip (e.g., Gorum et al. 2014). To test this possibility, in addition to the aspect map, we also calculated the relative position of each grid cell with respect to lines passing through the rupture planes (i.e., relative position with respect to rupture zone, RP). We made this calculation in the same manner as the aspect, being expressed in degrees ranging from 0 to 360. Ultimately, we used a landform classification to examine if co-seismic landslides are predominantly associated with a particular landform. For instance, the upper part of the slopes, near the ridge crests, should theoretically correspond to landforms where co-seismic landslides are abundant because of topographic amplification (e.g., Rizzitano et al. 2014). We obtained landform classes using the r.geomorphon module (Jasiewicz and Stepinski 2013) available in GRASS GIS (Neteler and Mitasova 2013). This algorithm automatically identifies and returns 10 landform classes (flat, summit, ridge, shoulder, spur, slope, hollow, footslope, valley, and depression) on the basis of pattern recognition.
To initially prompt considerations on the spatial distribution pattern of co-seismic landslides, we calculated the landslide point density, by using the ArcGIS kernel density function. Specifically, we used a radius of 5 km to calculate the point density per square kilometer. We then identified five classes with equal sizes based on a quantile classification.
Ultimately, we also processed the PGA signal by creating equally spaced bins of 0.1 g and 0.2 g. We used this information as our reference and examined the variations of each other landscape and landslide characteristics per PGA bin. We chose to do so because, within the epicentral area, the landslide response may be saturated. As a result, any other predisposing factor may not stand out, living the ground motion to dominate the landslide scenario. Conversely, slicing the PGA map into equal bins offers a better perspective where the landslide characteristics can be analyzed specifically for portions of the landscape that underwent the same level of ground shaking.

Fig. 2
Overview of the landslide mapping procedure. The first column (a, d) reports the pre-Kaikōura situations, followed by the second column (b, e) depicting the same area after the Kaikōura earthquake. The third column (c, f) shows a sample of the mapped landslides. The two different rows simply report a detail captured from two different areas also indicated in Fig. 1

Landslide size characteristics
We examined the size statistics of co-seismic landslides via both probability and frequency density distribution of landslide planimetric areas (Malamud et al. 2004). We plotted the corresponding probability and frequency densities in a 2D plane described by landslide-area-bins and the corresponding non-cumulative probability/frequency density values. We calculated two size statistical measures (power-law exponent, β, and the landslide event magnitude scale, mLS) using the frequency-density distribution of landslides, which is argued to exhibit power-law scaling (Guzzetti et al. 2002;Malamud et al. 2004;Tanyaş et al. 2019b).
The first size statistical measure we computed is the β, which refers to the slope of the power-law tail. It is used not only for quantifying the total denudation caused by landsliding (Hovius et al. 1997) but also for the landslide hazard (Guzzetti et al. 2005).  show that the β of globally available co-seismic landslide inventories varies from 2.0 to 3.6, with a mean value of 2.5. The second size statistical measure we computed is the mLS, which is a proxy for the energy released during landsliding and strongly correlated with the total landslide area. To calculate mLS, we used the method proposed by , which represents an updated version of the original approach suggested by Malamud et al. (2004).

Comparison between similar landslide events
We examined if the area affected by landslides, the mLS, and other landslide size characteristics produced by the Kaikōura earthquake were less than what has been observed for earthquakes with similar magnitude. We did so by following two approaches. First, we tested the method proposed by Tanyaş and Lombardo (2019) to predict the area affected by landslides. The authors develop this method focusing on non-flat areas, which are morphologically assumed to be susceptible to landsliding (i.e., slope ≤ 5° and local relief ≤ 100 m) and examined landslide-affected areas with respect to PGA. They perform these analyses by examining 20 co-seismic inventories, out of which they found the 0.12-g contour to be the minimum PGA containing at least 90% of the landslides. In their work, Tanyas and Lombardo demonstrated that knowing the percentage of the area potentially susceptible to landslides solely on the basis of relief and slope steepness can indicate the PGA level within which the vast majority of landslides are confined into. In this work, we applied this approach because its success or failure could show if the examined landslide event associated with the Kaikōura earthquake is an anomaly or it is consistent with other examples from the literature.
Second, we analyzed landslide size characteristics, comparing them with respect to five other co-seismic landslide inventories (Chi-Chi, Denali, Wenchuan, Haiti, and Gorkha) ( Table 2). These   (Schmitt et al. 2017;Tanyaş et al. 2017). Also, we compared our findings with two other co-seismic landslide inventories compiled for the Kaikōura earthquake (Massey et al. 2018(Massey et al. , 2020. The polygon-based inventories are not available in digital form as of the 23rd of June 2021. And, only landslide source area centroids lacking landslide size information were shared for the inventories compiled by Massey et al. (2018Massey et al. ( , 2020. However, the landslide probability density is reported in the two articles Massey and co-authors published. Therefore, to keep our comparison consistent with their work, we digitized the curves they plotted against the landslide source area presented in Fig. 3b of Massey et al. (2020). This procedure may have induced some degree of uncertainty, but it was the only available solution to ensure that the co-seismic landslide inventory we mapped for Kaikōura could be compared to those compiled by other researchers. For the five co-seismic landslide events gathered from the literature, we examined seismic and morphologic characteristics of those sites and compared them with the Kaikōura case. We used (1) local relief, (2) area affected by co-seismic landslides, (3) average landslide area, and (4) total landslide area.
We use PGA bins as our reference layer to examine the relation between landslide occurrences and all other factors under consideration. In the analyses, we added a preprocessing step where we masked out flat areas (i.e., slope ≤ 5° and local relief ≤ 100 m, Tanyaş and Lombardo, 2019). As a result, for the identified sites, we examined the abovementioned landscape characteristics.
First, we calculated local relief following the same method described in "Controlling factors' overview." Second, we examined the area affected by co-seismic landslides, which is a measure of the extent of failed hillslopes. To calculate this property, we took the earthquake rupture zone as the origin and calculated the area of non-flat zones in a cumulative manner within each PGA bin where landslides were present. Third, we examined the average landslide size within each PGA bin to investigate their variation away from the rupture zone. Ultimately, we calculated the cumulative total During the analyses, we divided the study areas into two halves and perform all investigations both for the footwall and for the hanging wall separately, to better assess the landslide distribution. Generally, this is particularly important when studying landslides in thrust fault systems, because most of the failures tend to concentrate on hanging walls (e.g., Sato et al. 2007;Gorum et al. 2011;Tatard and Grasso 2013) because of the greater ground acceleration (Oglesby et al. 2000). Therefore, if we did not treat the footwall and hanging wall separately, some of the examined factors (e.g., distance to rupture zone) that we expressed in terms of mean values could end up being exceedingly smoothed when considering the entire study site.

Factors controlling landslide distribution
Most of the landslides triggered by the event were disrupted rock, debris, and soil falls and slides with a limited number of rock avalanches clustered along the Hapuku river, and lateral spreads, especially on road edges and river banks Massey et al. 2018Massey et al. , 2020. The spatial distribution of landslides shows a high concentration around the central sector of the study area. Specifically, landslides cluster around the Manakau, Upper Kowhai, Papatea Faults, and Jordan Thrust ( Fig. 3a and b), where steep hillslopes characterize most the local landscape features (Fig. 3c). Among the two PGA estimates, USGS's product provides relatively low values (~ 0.5-0.7 g) around this high landslide density zone compared to the one provided by Bradley et al. (2017) (~ 0.7-1.6 g) (Fig. 3d and e). High landslide concentrations are also observed around the North Leader Fault Zone and the northern part of the Humps Fault Zone, but for this area, values from the USGS's PGA map (~ 0.4-0.8 g) are overall larger than the Bradley et al.'s (2017) (~ 0.7-1.6 g). Overall, the intensity of ground motion is positively correlated with slope failures unless some other morphometric factors counterplay the shaking effect (e.g., Tanyaş and Lombardo 2019). This implies that the match between high ground shaking and high landslide density zones offers interesting characteristics to further investigate. In fact, 60% of the landslide population lies within 5 km distance from the fault surface rupture and 90% are located within a 10 km buffer zone, although landslides were triggered even more than 55 km far from the rupture zone (Fig. 3a). However, the spatial distribution of landslides is always controlled by multiple factors (morphometric, seismic and geotechnical features) and their relative interactions. Therefore, strong ground shaking alone may not explain the whole spatial distribution. The deterministic estimates of ground motion provided by the USGS ShakeMap system show that the maximum PGA value in the Kaikōura earthquake reached 1.04 g, whereas the lowest PGA value associated with landsliding was 0.09 g (Fig. 3d). The same dataset also shows that the PGA contours of 0.5 g and 0.7 g encompass 75% of landslide population. Based on the broadband simulations provided by Bradley et al. (2017), maximum and minimum PGA values associated with landslides are 1.59 g and 0.04 g, respectively (Fig. 3e). Bradley et al.'s (2017) data also shows that landslides are more equally distributed among PGA bins (Fig. 3e).
Taking aside the difference between PGA maps, high landslide densities are specifically observed around the junctions of Humps and Leader faults as well as Papatea and other faults (Fig. 3a). This may indicate that not only the ground motion but also the multiple fault rupturing mechanism might have influenced slope failures and cause high landslide densities.
We also examined the variation in slope steepness and local relief associated to landslide locations throughout PGA bins. Based on Bradley et al.'s (2017) dataset, the highest values for both of the morphometric variables are observed over the areas exposed to PGA larger than 0.8 g ( Fig. 4a and b). Within this zone, the median slope steepness varies between 35 and 40° and then from 0.8 g to lower ground shaking areas, it gradually decreases to ~ 15° (Fig. 4a). Local relief shows some fluctuations within the area experienced PGA > 0.8 g. For instance, the area bounded by 1.2 g and 1.4 g has median local relief of ~ 560 m, whereas for the other PGA bins (i.e., PGA > 0.8 g), the median values changes from 680 to 760 m (Fig. 4b). However, the area bounded by 1.2 g and 1.4 g constitutes less than 2% the area affected by landslides and 2.2% of the landslide population (Fig. 3e). Therefore, these observations show that relatively steeper morphological features observed in the area were affected by stronger ground shaking (e.g., PGA > 0.8 g). As a result, overall landslide sizes also gradually decrease from high to low ground shaking zones. Specifically, the largest landslide size in terms of median values is ~ 7000 m 2 within the area exposed to PGA > 1.4 g (Fig. 4c).
Landslides also appeared to cluster within specific landform types. Overall, 85% of landslides occur on slope, ridge, and spur (i.e., subordinate ridge; see Fig. 5a). This confirms the general observation that co-seismic landslides mostly occur on the upper portion of the topographic profile due to topographic amplification (e.g., Rizzitano et al. 2014). The landslide size distribution also indicates that large landslides most commonly initiated at summits and ridges (Fig. 5b). Moreover, summits exposed to higher PGA refer to hillslopes where we observe the largest landslides (Fig. 5b). Lithologically, the study area is quite uniform. More than 70% of the entire area is covered by greywacke (see Fig. 6). Consequently, 78% of co-seismic landslides occur in greywacke. The second most common (7%) lithologic unit associated with landslides is limestone. Limestone also refers to the unit where the average landslide size is at its peak (~ 12,000 m 2 ).
The aspect map of the study area also reveals a relation between rupturing direction and the exposition of failed hillslopes (Fig. 7). On one hand, the rupturing initiated from the SW part of the study area and propagated towards NE direction (Wang et al. 2018). On the other hand, co-seismic landslides are mostly seen on east-facing hillslopes. However, the dominant aspect of the examined area also Error bars indicate ± 1σ of average landslide size. Gw, greywacke; Al, undifferentiated floodplain alluvium; Ms, mudstone (weakly indurated); Cw, weakly consolidated conglomerate; Hs, sandstone (strongly indurated); Ss, sandstone (weakly indurated); St2, schist; Ar, argillite; In, ancient volcanoes, minor intrusives (dikes and sills); Lo, loess; Ls, limestone; St1, semi-schist; Um, ultramafics; Wb, windblown sand; Gn, plutonics; Vo, lavas and welded ignimbrites; Tb, pyroclastics (ash and lapilli); Pt, peat points out southeast to east. To further investigate this observation, we also examined the distribution of aspects separately based on relative position with respect to rupture zone (RP). This analysis shows that co-seismic landslides are still abundant on hillslopes facing northeast, east, and southeast (Fig. 7), although in some sectors, the dominant landscape does not dictate these directions (e.g., west, northwest, southwest and southeast sectors). Also, in some sectors, multiple high concentration zones exist for the aspect of the landscape and, yet, landslides still facing east and southeast directions (e.g., north and northeast sectors). Therefore, the aspect of the failed slope is most likely associated with the rupturing direction (e.g., Del Gaudio and Wasowski 2007).

Comparison of the Kaikōura landslides with previously occurred seismic landslide event
To predict the area affected by landslides triggered by the Kaikōura earthquake, we examined the area bounded by the 0.12-g contour as indicated in "Comparison between similar landslide events." We then calculated the extent of landslide susceptible area. Figure 8a shows that on the basis of the PGA maps provided by the USGS ShakeMap system (Worden andWald 2016) andBradley et al. (2017), 44% of the area affected by the Kaikōura earthquake is morphologically susceptible to landslides (i.e., areal coverage of landslide-susceptible areas, ACLSA = 44%). Second, we used this measure to predict the PGA value covering of the majority of landslide population (blue and pink contour lines in Fig. 8a). Both contour lines correspond to the 0.36-g PGA (Fig. 8b), which in turn encompasses 97% and 88% of the total landslide population, based on USGS ShakeMap system's and Bradley et al.'s (2017) products, respectively. This shows that the extent of landslide affected area is overall consistent with the observation from the literature. In fact, the method proposed by Tanyaş and Lombardo (2019) aims to delineate the area where at least 90% of the landslide population occurred. The fact that the two estimate satisfy or are close to satisfy the 90% requirement further implies that the limits of landslide affected area are predominantly controlled by morphologic and seismic characteristics similar to other co-seismic landslide events, previously examined in the literature (list labeled in Fig. 8b).
We also compared our co-seismic landslide inventory with two previously complied inventories for the same earthquake (Massey et al. 2018(Massey et al. , 2020. For this comparison, we plotted the probability density distributions of all inventories. Figure 9a shows that Massey et al. (2018Massey et al. ( , 2020's inventories are rich in small size landslides. This is a natural consequence of the higher resolution images (e.g., 0.3 m ground resolution orthorectified air photographs) they used to map landslides. Massey et al. (2020) provided these trends on the basis of landslide source areas, whereas in our inventory, source and deposit of landslides are not differentiated. As a result, we used the area of entire landslide polygons to create the probability density curve. Therefore, this is most likely another reason why our inventory includes a smaller number of small landslides compared to Massey et al. (2020).
The power-law exponent we calculated for our inventory is − 2.9, whereas it is − 1.9 and − 2.1 in Massey et al.'s (2018Massey et al.'s ( , 2020 inventories. The same reasons we mentioned above are valid to explain the difference between power-law exponents. Also, mapping preferences could influence the power-law exponents, because estimation of a power-law exponent is quite sensitive to the mapping criteria (Tanyaş et al. 2019b).
We also plotted frequency-area distribution (FAD) curves for five additional co-seismic landslide inventories (Chi-Chi, Denali, Haiti, Wenchuan and Gorkha), each one of these being constituted by polygons representing the whole landslide body (i.e., both source and deposit). We then compared these FADs with respect to the inventory we mapped for the Kaikōura earthquake (Fig. 9b). Results show that the power-law exponent of our inventory (β = − 2.9) is similar to the ones estimated for the Gorkha (β = − 2.8) and Wenchuan (β = − 3.1) cases. Also, we calculated the landslide-event magnitude for each inventory. Our findings show that the Kaikōura landslide event (mLS = 4.9) has a similar magnitude to the Denali (mLS = 4.9), Chi-Chi (mLS = 5.1), and Gorkha (mLS = 4.6) ones.
Moreover, we compared the Kaikōura EQIL inventory to the same five co-seismic landslide inventories, this time focusing on landscape characteristics. Specifically, we examined the spatial variance of local relief, cumulative landslide affected area, average landslide area, and cumulative landslide area with relation to PGA bins, for both hanging and foot walls (Fig. 10). We did this because in the Kaikōura earthquake, all the landslides are essentially located within the hanging wall, if we consider the megathrust causing the earthquake. To split the study area into two parts corresponding to the hanging and foot walls, we used the fault surface ruptures provided in Litchfield et al. (2018). Figure 10a shows how local relief varies away from the rupture zone, across each PGA bin. On the hanging wall, where more than 60% of landslides occurred, the average local relief is relatively lower (< ~ 500 m) than most other cases, where PGA < 0.7 g. For this zone, the only case where the local relief is smoother than it is in Kaikōura corresponds to the landscape affected by the Haiti (< ~ 400 m) earthquake. Also, within this zone, local relief in Wenchuan and Gorkha (~ 500-1000 m) appears to be significantly higher than the Kaikōura case. However, in the Kaikōura case, the topography becomes steeper for the areas exposed to higher ground shaking (PGA > 0.7 g) and the local relief reaches comparable levels to those observed in the Gorkha case.
As for the footwall, in the Kaikōura site, the local relief is still steeper than Chi-Chi, Denali, and Haiti cases overall. For the area that experienced PGA > 0.8 g, the local relief values significantly fluctuate and in some zones it becomes even steeper than the Wenchuan site, especially for the ground shaking confined by the 0.8-g PGA contour. In fact, this PGA interval also refers to the zone where most co-seismic landslides were triggered (Fig. 3d). In turn, we can state that the Kaikōura site has a relatively steeper topography than Haiti and Denali but smoother than the other cases in terms of local relief (Fig. 10b). Figure 10c shows the variation in the cumulative areas affected by earthquakes through PGA bins. It also indicates that Kaikōura and Haiti are the only earthquakes that have shaken the landscape with greater PGA than 0.15 g. The Wenchuan earthquake then follows as the second strongest one (PGA max = ~ 1.2 g). In both Kaikōura and Haiti, the highest ground shaking zone (e.g., PGA max > ~ 1.2 g) is only effective inside the epicentral region (the innermost ~ 100-150 km 2 ) and the affected area gets narrower away from this zone. Overall, the total area affected by Kaikōura (~ 13,000 km 2 ) earthquake is compatible with Haiti (~ 11,000 km 2 ) and Chi-Chi (~ 12,000 km 2 ) earthquakes, whereas in the Wenchuan, Gorkha, and Denali cases, the areas affected by landslides are larger (Fig. 10d).
The average landslide size in the Kaikōura earthquake (~ 7000 m 2 ) is mostly similar, but overall relatively larger than the size associated to the landslides triggered by the Wenchuan earthquake (~ 6000 m 2 ) in both hanging wall and footwall (Fig. 10e). Among the examined cases, the Denali earthquake appears as an outlier with significantly larger average landslide sizes (~ 30,000 m 2 ).

Fig. 10
Comparison between co-seismic landslide inventories. The yellow star schematically represents the rupture zone and separates the plots representing foot (left column) and hanging (right column) walls' conditions. Each row reports a different parameter, separated for foot and hanging walls, whereas the central panels report the combination of both. Row a and panel b show the local relief; row c and panel d show the cumulative area affected by the earthquake; row e and panel f show the average landslide area; row g and panel h show the cumulative landslide area. Notably, the considered earthquakes may have different rupturing mechanism; these are summarized in the legend at the bottom for thrust, strike-slip and complex mechanisms including both Ultimately, we compared the cumulative landslide areas of inventories. Figure 10g shows that the cumulative landslide area estimated for the Kaikōura case is similar to the Denali in both hanging and foot walls, and overall, it is the fourth largest landslide event (100 km 2 ) in terms of total landslide area, after Wenchuan, Chi-Chi, and Denali (Fig. 10h).

Discussion and conclusions
This study presented a new landslide inventory for the 2016 Kaikōura earthquake. Two other landslide inventories for the same earthquake had already been introduced in a successive manner as version 1.0 and version 2.0 (Massey et al. 2018(Massey et al. , 2020. However, except for the landslide source area centroids lacking landslide size information (Massey et al. 2018(Massey et al. , 2020, the full inventory has not been made publicly available as of June 2021. Therefore, to share this inventory and to examine the distribution of landslides and their size statistics, we created our version of the co-seismic landslide inventory of the Kaikōura earthquake. We used optical Sentinel-2 satellite images with 10 m spatial resolution for landslides mapping. Because of this, we might have missed some small landslides or mapped them as part of relatively larger landslide bodies. In particular, the probability density plot of our landslide inventory (Fig. 9a) shows that the rollover point (i.e., the point below which probability density values decrease for smaller landslides) is around 1000 m 2 . This implies that 1000 m 2 corresponds approximately to the smallest landslide area for which we can assume to have completely mapped landslides' size in our inventory (Tanyaş et al. 2019b). Notably, this is larger than the rollover point (100 m 2 ) reported for the two other landslide inventories described by Massey et al. (2018Massey et al. ( , 2020 and also mapped for the Kaikōura earthquake (Fig. 9a). This difference may primarily reside in the fact that Massey and co-authors mapped the source area separately from the depositional one.
The completeness of an EQIL inventory can be examined by a semi-quantitative method called Completeness Index (Tanyaş and Lombardo 2020). This method suggests examining a landslide inventory in relation to two quantitative parameters, addressing two qualitative questions: -Is there a minimum size threshold used for mapping landslides (i.e., is the rollover size smaller than 500 m 2 )? -Are landslides mapped for the entire landslide-affected area or just for a subregion?
Based on this method, our inventory can be classified with "moderate completeness" because we examined the entire area affected by landslides but the rollover size is larger than 500 m 2 (Tanyaş and Lombardo 2020).
The compiled landslide inventory allowed us to (1) examine the factors controlling the landslide distribution and (2) make a comparison with other co-seismic landslide inventories. Consequently, we tested whether the Kaikōura landslide event is smaller or less hazardous overall than other events associated with strong earthquakes. Our findings show that the distribution of landslide is mainly controlled by the coupling effect of seismic shaking and steepness of topography, as one can expect. Therefore, the Kaikōura landslide event was not exceptional, nor in terms of landslide magnitude nor in terms of area affected by landslides. This being said, we should also stress that our inventory has moderate completeness and we might have missed some small landslides. Therefore, the number of landslides we identified could be less than the actual number of landslides. However, the identified landslide event magnitude and the area affected by landslides should be reflecting the actual event as neither of them fully depends on small landslides.
Our findings show that the areal boundary of the landslide affected area was successfully predicted (Fig. 8) based on the method developed by previous observations available in the literature (Tanyaş and Lombardo 2019). In fact, this supports the argument that the landslide event is not exceptional case. Also, to expand on the landslide magnitude and area affected statement made above, we report below comparative metrics with other coseismic landslide events. For instance, the Wenchuan event has the highest mLS (6.1) and both topographic steepness (Fig. 10b) and area affected by seismic shaking (Fig. 10d) are larger in the Wenchuan case compared to the Kaikōura (mLS = 4.9). The Chi-Chi event has the second-largest mLS (5.1), which is quite similar to the Kaikōura event in terms of mLS. In fact, the overall steepness of the terrain and boundary of the seismically affected area is also quite similar to the Kaikōura event (Fig. 10). The Denali earthquake has the same mLS (4.9) as our case. Actually, our study area presents a slightly steeper topography and yet the area affected by the Denali earthquake (~ 48,000 km 2 ) is approximately four times larger than the Kaikōura (~ 14,000 km 2 ). Therefore, the larger extension of the area affected by the Denali earthquake explains why its mLS is equal to the Kaikōura case despite the relatively smooth topography of the Denali area. The Haiti earthquake has the smallest mLS (4.3), and as we can expect, it presents both the smoothest topography and the smallest areal extent affected by earthquake (Fig. 10). The only case that could not be fully explained by this simple comparison rationale is the 2015 Gorkha earthquake (mLS = 4.6). The Gorkha case has one of the steepest terrain conditions among the six cases (Fig. 10b). Also, the area affected by the earthquake is (~ 57,000 km 2 ) significantly larger than the Kaikōura (Fig. 10d). However, the maximum PGA values recorded in the Gorkha earthquake (PGA max = ~ 0.9 g) is lower than the Kaikōura earthquake (PGA max = ~ 1.5 g) and this could be the reason for having a slightly smaller mLS than the Kaikōura earthquake.
In terms of landslide size statistics, the average landslide size in the Kaikōura earthquake is even larger than in the Wenchuan case. The same situation is also valid for the Chi-Chi and the Denali earthquakes; their average landslide sizes are larger than the ones in the Wenchuan. This could be caused by the characteristics of ground motion controlling the size of landslides. For instance, Jibson and Tanyaş (2020) show that increasing periods and duration of seismic ground motion are the physical factors driving increased landslide sizes for large earthquakes. In this regard, they show that, for instance, the median duration of the Denali earthquake is longer than the Wenchuan earthquake.
Notably, comparison between different landslide events is usually no trivial task because each landslide event has its own eventspecific conditions characterized by different seismic, morphologic, climatic, and anthropogenic factors. For example, the 2008 Wenchuan and the 2015 Gorkha earthquakes are considered two events Landslides 19 & (2022) 1417 that occurred in similar geomorphic (Kargel et al. 2016) and seismotectonic settings having comparable earthquake magnitudes (Goda et al. 2015). Nevertheless, as we also presented above, the resultant landslide events are quite different in terms of total landslide area and the number of landslides . This implies that even if there are some similarities between events, the distribution and characteristics of co-seismic landslides still vary depending on, for instance, existing of a surface rupture (Xu 2014) or asperity zones (Gorum et al. 2014), different faulting style, fault geometry (Tatard and Grasso 2013;Gorum et al. 2014;Gorum and Carranza 2015), and characteristics of ground shaking associated with directivity (e.g., Del Gaudio and Wasowski 2007;Barth et al. 2020) or topographic amplification (Sepúlveda et al. 2005). This also means that our analyses are not adequate to present the relationship between the complex rupturing mechanism and spatial/size distribution of landslides because we do not have detailed characteristics of the ground motion such as frequency and duration content of seismic waves (Jibson et al. 2004;Jibson and Tanyaş 2020). Additional studies considering physics-based ground-motion simulations in 3D geologic structures are needed to explore the link between the rupturing mechanism and spatial/size distribution of landslides.