Beyond 2D inventories : synoptic 3D landslide volume calculation from repeat LiDAR data

Efficient and robust landslide mapping and volume estimation is essential to rapidly infer landslide spatial distribution, to quantify the role of triggering events on landscape changes and to assess direct and secondary landslide-related geomorphic hazards. Many efforts have been made during the last decades to develop landslide areal mapping methods, based on 2D satellite or aerial images, and to constrain empirical volume-area (V-A) allowing in turn to offer indirect estimates of landslide 10 volume. Despite these efforts, some major issues remain including the uncertainty of the V-A scaling, landslide amalgamation and the under-detection of reactivated landslides. To address these issues, we propose a new semi-automatic 3D point cloud differencing method to detect geomorphic changes, obtain robust landslide inventories and directly measure the volume and geometric properties of landslides. This method is based on the M3C2 algorithm and was applied to a multi-temporal airborne LiDAR dataset of the Kaikoura region, New Zealand, following the Mw 7.8 earthquake of 14 November 2016. 15 We demonstrate that 3D point cloud differencing offers a greater sensitivity to detect small changes than a classical difference of DEMs (digital elevation models). In a small 5 km2 area, prone to landslide reactivation and amalgamation, where a previous study identified 27 landslides, our method is able to detect 1431 landslide sources and 853 deposits with a total volume of 908,055 ± 215,640 m and 1,008,626 ± 172,745 m, respectively. This high number of landslides is set by the ability of our method to detect subtle changes and therefore small landslides with a carefully constrained lower limit of 20 m2 (90% with 20 A<300 m2). Moreover, the analysis of landslide geometric properties reveals the absence of a rollover in the landslide area distribution, which is a feature classically described in the literature. This result suggests that the rollover behaviour previously observed is due to an under detection of small landslides. Reactivated landslides represent 27.2 % of the total landslide source area and 29.9 ± 12.8 % of the total volume. Reactivated landslides are located in areas where landslide mapping methods based on 2D images are assumed to perform poorly due to the weak contrast in texture and colour between the two epochs. Our result 25 therefore suggests that the number, area and volume of landslides can be significantly under-estimated by these methods. To our knowledge, this is the first approach to create a regional landslide inventory map from 3D point cloud differencing. Our results call for a more systematic use of high-resolution 3D topographic data to assess the impact of extreme events on topographic changes in regions prone to landsliding and to infer the geometric scaling properties of landslides. 30 https://doi.org/10.5194/esurf-2020-73 Preprint. Discussion started: 15 September 2020 c © Author(s) 2020. CC BY 4.0 License.


Abstract.
Efficient and robust landslide mapping and volume estimation is essential to rapidly infer landslide spatial distribution, to quantify the role of triggering events on landscape changes and to assess direct and secondary landslide-related geomorphic hazards. Many efforts have been made during the last decades to develop landslide areal mapping methods, based on 2D satellite or aerial images, and to constrain empirical volume-area (V-A) allowing in turn to offer indirect estimates of landslide 10 volume. Despite these efforts, some major issues remain including the uncertainty of the V-A scaling, landslide amalgamation and the under-detection of reactivated landslides. To address these issues, we propose a new semi-automatic 3D point cloud differencing method to detect geomorphic changes, obtain robust landslide inventories and directly measure the volume and geometric properties of landslides. This method is based on the M3C2 algorithm and was applied to a multi-temporal airborne LiDAR dataset of the Kaikoura region, New Zealand, following the Mw 7.8 earthquake of 14 November 2016. 15 We demonstrate that 3D point cloud differencing offers a greater sensitivity to detect small changes than a classical difference of DEMs (digital elevation models). In a small 5 km² area, prone to landslide reactivation and amalgamation, where a previous study identified 27 landslides, our method is able to detect 1431 landslide sources and 853 deposits with a total volume of 908,055 ± 215,640 m 3 and 1,008,626 ± 172,745 m 3 , respectively. This high number of landslides is set by the ability of our method to detect subtle changes and therefore small landslides with a carefully constrained lower limit of 20 m² (90% with 20 A<300 m²). Moreover, the analysis of landslide geometric properties reveals the absence of a rollover in the landslide area distribution, which is a feature classically described in the literature. This result suggests that the rollover behaviour previously observed is due to an under detection of small landslides. Reactivated landslides represent 27.2 % of the total landslide source area and 29.9 ± 12.8 % of the total volume. Reactivated landslides are located in areas where landslide mapping methods based on 2D images are assumed to perform poorly due to the weak contrast in texture and colour between the two epochs. Our result 25 therefore suggests that the number, area and volume of landslides can be significantly under-estimated by these methods. To our knowledge, this is the first approach to create a regional landslide inventory map from 3D point cloud differencing. Our results call for a more systematic use of high-resolution 3D topographic data to assess the impact of extreme events on topographic changes in regions prone to landsliding and to infer the geometric scaling properties of landslides.

Introduction
In mountain areas, extreme events such as large earthquakes and typhoons can trigger thousands of landslides. Landslides are a key agent of hillslope and landscape erosion (Keefer, 1994;Malamud et al., 2004) and represent a significant hazard for local populations. Efficient and rapid mapping of landslides is required to robustly infer their spatial distribution, their total volume, the induced landscape changes and the associated direct and secondary hazards (Guzzetti et al., 2012;Hovius et al., 1997;35 Marc et al., 2016;Parker et al., 2011). Following a triggering event, total landslide volume over a regional scale is classically determined in two steps: (i) individual landslide mapping using 2D satellite or aerial images (e.g., Behling et al., 2014;Fan et al, 2019;Guzetti et al., 2012;Li et al., 2014;Martha et al., 2010;Massey et al., 2018;Parker et al. 2011) and (ii) indirect volume estimation using a volume-area relationship (e.g. Simonett, 1967;Larsen et al., 2010):

=
(1) 40 with V and A the volume and area of individual landslides, a prefactor and a scaling exponent ranging between 1.3 and 1.6 (e.g., Larsen et al., 2010, Massey et al., 2020. A first source of error comes from the uncertainty on the values of and which tend to be site specific and process specific (e.g. shallow versus bedrock landsliding). This uncertainty could lead to an order magnitude of difference in total estimated volume given the non-linearity of eq. (1) (Larsen et al., 2010). Two other sources of error arise from the detectability of 45 individual landslides themselves and the ability to accurately measure the distribution of landslide areas, due to landslide amalgamation and under-detection of reactivated landslides. Landslide amalgamation can produce up to 200 % error in the total volume estimation (Li et al., 2014;Marc and Hovius, 2015) and occurs because of landslide spatial clustering or incorrect mapping due, for instance, to automatic processing. Indeed, automatic landslide mapping (Behling et al., 2014;Marc et al., 2019;Martha et al., 2010;Pradhan et al., 2016) relies on the difference in texture, color and spectral properties such as NDVI 50 (Normalized difference vegetation index) of multispectral 2D images between pre-and post-landslide images, assuming that landslides lead to vegetation removal or significant texture change. During this process, difficulties in automatic segmentation of landslide sources can result in incorrect estimate of individual landslide area, which propagates into a much larger estimate of volume owing to the non-linearity of eq. (1). Manual mapping and automatic algorithms based on geometrical and topographical inconsistencies can enable to reduce the amalgamation effect on landslide volume estimation (Marc and Hovius, 55 2015), but it remains a source of error due to the inherent spatial clustering of landslides and the overlapping of landslide deposits and sources. Under-detection of reactivated landslides occurs because the spectral signature of images is not sufficiently altered by a new failure. This phenomenon is particularly important in areas with thin soils and sparse or total lack of vegetation (Behling et al., 2014). Guzetti et al. (2009) show that the proportion of landslide volume mobilized by reactivations can reach 62% for an individual snowmelt or rainfall events. Yet, the level of landslide reactivation in a given 60 inventory remains generally largely unknown. The delimitation of reactivated landslides is therefore critical to robustly infer total landslide volume. To better detect reactivated landslides, Behling et al. (2014Behling et al. ( , 2016) developed a method using temporal NDVI-trajectories which describe temporal footprints of vegetation changes but cannot fully address complex cases when texture is not significantly changing such as bedrock landsliding on bare rock hillslopes.
Addressing these three sources of uncertainty -volume-area scaling uncertainty, landslide amalgamation and the under-65 detection of landslide reactivation -requires to explore new approaches to obtain and analyse landslide inventories. In the last decade, the increasing availability of multi-temporal high resolution 3D point cloud data and digital elevation models (DEM), based on aerial or satellite photogrammetry and Light Detection and Ranging (LiDAR), has opened the possibility to better quantify landside volume and displacement (Bull et al., 2010;Mouyen et al., 2019;Okyay et al., 2019;Passalacqua et al., 2015;Ventura et al., 2011). The most commonly used technique is the difference of DEM (DoD) which consists in computing 70 the vertical elevation differences between two DEMs of different time (Corsini et al., 2009;Giordan et al., 2013;Mora et al., 2018;Wheaton et al., 2010). Even though this method is fast and works properly on planar surface, a vertical difference can be prone to strong errors when used to quantify changes on vertical or very steep surfaces where landsliding typically occurs (e.g., Lague et al., 2013). The "Multiscale model-to-model cloud comparison" (M3C2) algorithm implemented by Lague et al. (2013) rather considers a direct 3D point cloud comparison. This algorithm has three main advantages over a DoD: (i) it 75 operates directly on 3D point clouds, avoiding a phase of DEM creation that is conducive to a loss of resolution imposed by the pixel size and potential data interpolation, (ii) it computes 3D distances along the normal direction of the topographic surface, allowing to better capture subtle changes on steep surfaces, (iii) it computes a spatially variable confidence interval that accounts for surface roughness, point density and uncertainties in data registration. Applicable to any type of 3D data to measure the orthogonal distance between two point clouds, this approach has generally been used for Terrestrial Lidar and 80 UAV photogrammetry over sub-kilometer scales. In the context of landsliding, it has been used to infer the displacement and volume of individual landslides, using point clouds obtained by UAV photogrammetry (e.g., Esposito et al., 2017;Stumpf et al., 2015), as well as for rockfall studies (Benjamin et al.2016;Williams et al., 2018). Yet, to our knowledge, the systematic detection and segmentation of hundreds of landslides, at a regional scale, from 3D point cloud have not yet been attempted.
Here, we produce an inventory map of landslide topographic changes using a semi-automatic 3D point cloud differencing 85 method based on M3C2 and applied to multi-temporal airborne LiDAR data. We investigate the potential of this method to robustly infer landslide volumes in a region prone to landslide reactivation and amalgamation. We applied our method to a complex topography located near Kaikoura, New Zealand, where the 2016 Mw 7.8 earthquake triggered nearly 30,000 landslides over a 10,000 km² area (Massey et al., 2020). We choose a 5 km² area characterized by a high landslide spatial density along the Conway segment of the Hope fault, inactive during the earthquake, where pre-and post-earthquake LiDAR 90 were available and investigate the proportion of reactivated landslide volume on the total budget (Fig. 1). This area has a variety of vegetation cover (e.g. trees, shrubs, bare bedrock) and pre-existing landslides, and typically represents a challenge for conventional 2D landslide mapping. We illustrate the benefits of working directly on 3D data to generate landslide inventories, and discuss the methodological advantages to operate directly on point clouds with M3C2 compared to DoD in terms of detection accuracy and error for total landslide volume. 95 The paper is organized as followed: first, the LiDAR dataset is presented followed by a description of the 3D point cloud differencing method. Second, results of the geomorphic changes detection and landslides individualization method applied to the study area are presented. Detected landslides properties are then investigated in terms of area and volume, specifically for reactivated and new landslides. Finally, current limitations of the method are discussed as well as knowledge gained on the importance of landslide reactivation on co-seismic landslide inventory budget. 100 Figure 1: Maps of the regional context, location and visualization of the study area. (a) Regional map of Kaikoura with the location of the 2016 Mw 7.8 earthquake, associated active faults and the study area. (b) Orthophotos focused on the study area dated before and after the earthquake with the 5 km² LiDAR dataset extent used in this paper (all imageries are available at https://data.linz.govt.nz/set/4702-nz-aerial-imagery/, Aerial survey 2017).

Data description
In this study, we compare two 3D point clouds obtained from airborne LiDAR data before and after the November 14 2016 Kaikoura earthquake (Table. 1). Both airborne LiDAR surveys were acquired during summer. Pre-earthquake LiDAR data represents a combination of six flights performed from March 13, 2014 to March 20, 2014 for a resulting ground point density of 3.8 ± 2.1 pts/m². The vertical accuracy of this dataset has been estimated at 0.068 m to 0.165 m from check points located 110 on highways (Dolan and Rhodes, 2016). The post-earthquake LiDAR survey took place rapidly after the earthquake from December 3, 2016 to January 6, 2017 for an average ground point density of 11.5 ± 6.8 pts/m². Vertical accuracy of this dataset https://doi.org/10.5194/esurf-2020-73 Preprint. Discussion started: 15 September 2020 c Author(s) 2020. CC BY 4.0 License. assessed on ground classified points is 0.04 m (Aerial survey, 2017). The difference in acquisition dates represents a period of 2 years and 8 months. The vegetation of both datasets was removed using the classification provided by the downloaded data to keep only ground points. In addition to LiDAR data, orthophotos were used to visually validate the detection of landslides 115 from the 3D approach for new landslide sources. The pre-earthquake orthophoto was obtained on January 24 2015 (available at https://data.linz.govt.nz/layer/52602-canterbury-03m-rural-aerial-photos-2014-2015) and the post-earthquake one on December 15 2016. The resolutions are 0.3 and 0.2 m, respectively. The method developed here to detect landslides consists in a 3D point cloud differencing between two epochs using the M3C2 algorithm (Lague et al., 2013) available in the Cloudcompare software (Cloudcompare v2.11, 2020). This algorithm estimates 125 orthogonal distances along the surface normal directly on 3D point clouds without the need for surface interpolation or gridding. While M3C2 can be applied on all points, the algorithm can use an accessory point cloud, called core points, of arbitrary geometry which we impose in our case to be a regular grid with constant horizontal spacing generated by the rasterization of one of the two clouds. All the M3C2 calculations will be done in 3D using the raw point clouds, but the results will be "stored" on the core points. The use of a regular grid of core points has four advantages: (i) a regular sampling of the 130 results which allows to compute robust statistics of changes unbiased by spatial variations in point density; (ii) it facilitates the volume calculation and the uncertainty assessment on total volume; (iii) it can be directly reused with 2D GIS as a raster; and https://doi.org/10.5194/esurf-2020-73 Preprint. Discussion started: 15 September 2020 c Author(s) 2020. CC BY 4.0 License.
(iv) it speeds up calculations, although in the proposed workflow, computation time is not an issue and can be done on a regular laptop.
The first step of M3C2 consists in the calculation of a 3D surface normal for each core point at a scale D (called the normal 135 scale) by fitting a plane to the points of the first dataset located within a radius of size D/2 of the core point. Once the normal vectors are defined, the local distance between the two clouds is computed for each core point as the distance of the average positions of the two point clouds at a scale d (projection scale). This is done by defining a cylinder of radius d/2, oriented along the normal with a maximum length pmax. Distances are not computed if no intercept is found in the second point cloud. M3C2 also provides uncertainty on the computed distance at 95% of confidence based on local roughness, point density and 140 registration error as follow: where LoD95% is the Level of Detection, σ1(d) et σ2(d) are the detrended roughness of each cloud at scale d measured along the normal direction, n1 and n2 are the number of points and reg is the co-registration error between the two epochs, assumed 145 spatially uniform and isotropic in our case, but which could be spatially variable and anisotropic (James et al., 2017). M3C2 has the option to compute the distance vertically which bypasses the normal calculation, and we use this option several time in the workflow. We use the abbreviation vertical-M3C2 in that case and 3D-M3C2 otherwise.

Parameters selection and 3D point cloud differencing performance
In this section, we explain how to select the appropriate normal scale D and projection scale d to detect landslides using M3C2. 150 The normal scale D should be large enough to encompass enough points for a robust calculation, and smooth out small scale point cloud roughness that result in normal orientation flickering and overestimation of the distance between surfaces (Lague et al., 2013). However, D should be small enough to track the large scale variations in hillslope geometry. By studying roughness properties of various natural surface, Lague et al. (2013) proposed a criterion for which the ratio of the normal scale and the surface roughness measured at the same scale should be larger than about 25. We thus set D as the minimum scale for 155 which a majority of core points verify this criteria (details of the analysis can be found in section S1.1 in the supplements). As roughness is a scale and point density dependent measure, we explore a range of normal scales for the pre-earthquake dataset which has the lowest point density. We found that D ~ 10 m represents a threshold scale below which the number of core points verifying the roughness criteria significantly drops.
The projection scale d should be chosen such that it is large enough to compute robust statistics using enough points, but small 160 enough to avoid spatial smoothing of the distance measurement. Following Lague et al. (2013), M3C2 computes eq. (2) only if 5 points are included in the cylinder of radius d/2 for each cloud. In our case, the pre-EQ data with the lowest point density will in practice set the value of d. We ran 3D-M3C2 with a normal scale D = 10 m and d varying from 1 to 40 meters, and https://doi.org/10.5194/esurf-2020-73 Preprint. Discussion started: 15 September 2020 c Author(s) 2020. CC BY 4.0 License.
found that for d=5 m, eq. (2) can be computed with at least 5 data points for 95 % of the surface. The remaining 5 % of points typically correspond to area of low ground point density below dense vegetation. It was not deemed interesting to increase d 165 above 5 m for all points as it would deteriorate the ability to detect small landslides. However, to be able to generate M3C2 confidence intervals for as many points as possible, in particular on steep slopes below vegetation, we use a second pass of M3C2 with d=10 m using the core points for which no confidence interval was calculated at d=5 m.
For this study case, the lower ground point density of the pre-EQ lidar data controls D and d. If the pre-EQ dataset had a point density similar to the post-EQ data, values of D = 4 m and d = 3.5 m could have been used, improving the spatial resolution 170 of the change detection. Finally, the maximum cylinder length pmax was set to 30 m as it allowed to compute the maximum change observed in the study area. This is generally obtained by trial and error. Setting pmax too large, increases significantly computation time and may result in two different surfaces of the same point cloud being averaged (e.g., near very steep divides or in narrow gorges).

3D Landslide mapping workflow 175
Our 3D landslide mapping workflow is divided in four main steps (Fig. 2).

Registration
To detect geomorphic changes and landslides, the two datasets need to be co-registered as closely as possible and any large-185 scale tectonic deformation need to be adjusted. No significant intra-dataset deformation was detected, but the produced LiDAR datasets did have a vertical shift larger than 1 m. To correct for it, a grid of core points is first created by rasterizing the dataset https://doi.org/10.5194/esurf-2020-73 Preprint. Discussion started: 15 September 2020 c Author(s) 2020. CC BY 4.0 License.
with the largest point density. Then, a vertical-M3C2 calculation is performed and the mode of the resulting distribution is used to adjust the two datasets by a vertical shift of 1.36 m. This approach is valid only when the fraction of the surface affected by landsliding is small. A subsequent 3D-M3C2 calculation is performed to obtain a preliminary map of geomorphic change. 190 By manually selecting a threshold of change, one can identify areas that are assumed to be stable. In our case, this corresponds to regions displaying a change smaller than 0.6 m. An Iterative Closest Point (ICP) algorithm (Besl and McKay, 1992) is performed on the stable areas, and the obtained rigid transformation is applied to the entire post-earthquake point cloud to align it with the pre-earthquake one. At this stage, the two datasets are considered optimally registered for the stable areas but with an unknown registration error reg. To estimate this error, the M3C2 algorithm is re-applied on the registered data, and 195 the registration error is estimated as the standard deviation of the 3D-M3C2 distances measured on stable areas. At this stage, a 3D map of topographic change is available, but the statistically significant geomorphic change and individual landslides have not been isolated.
The registration error is defined as the standard deviation of the 3D-M3C2 distances measured between the pre-and post-event point clouds, considering only stable areas that are uniformly distributed in the studied region. These were manually selected 200 and represent 23% of the area (Fig. 3). The mean 3D-M3C2 distance is -0.01 m, showing that there is almost no bias in the registration. 95% of the measured distances are within a range of -0.33 to 0.35 m (Fig. 3). The standard deviation is 0.17 m, and we thus set reg to 0.17 m. According to eq. (2), LoD95% cannot be smaller than 0.33 m in the ideal case of negligible roughness surface. The relatively large registration error in this dataset is related to errors related to flight line imperfect alignment, some classification residual errors (e.g. low vegetation classified as ground) and residual errors related to the ICP 205 ( Fig. 3).

Geomorphic change detection 210
The registration error reg is then used in a first 3D-M3C2, using the pre-determined projection scale d=5m, to estimate the spatially variable LoD95% according to eq. (2). Then, a second 3D-M3C2 is performed at a larger projection scale (d=10m) https://doi.org/10.5194/esurf-2020-73 Preprint. Discussion started: 15 September 2020 c Author(s) 2020. CC BY 4.0 License. only on core points for which a confidence interval could not be estimated due to the insufficient points at d=5m. These points generally correspond to ground points under canopy on steep slopes, and represented typically 5 % of the core points over the entire area, but ~ 15-20 % of steep slopes prone to landsliding. The statistically significant geomorphic changes at the 95% 215 confidence interval are then obtained by considering core points for which the 3D-M3C2 distance is larger than LoD95%. This process highlights changes associated to all geomorphic processes, including landsliding, but also erosion and deposition processes in the fluvial domain. Detected changes located in the river, and specifically related to river dynamics are removed manually to select only landslides sources and deposits. This phase is performed by manual visual inspection to prevent landslide deposited at the bottom of valleys to be removed. 220

Landslide source and deposit segmentation
Prior to landslide segmentation, a vertical-M3C2 is performed on the remaining erosion and deposition points located on hillslopes in order to estimate landslide volume afterwards. This operation will be described in the next section. Core points with significant change are segmented by a connected component analysis (Fig.2.c). This technique segments a point cloud into compact sub-clouds based on a minimum distance threshold Dm and a minimum number of points per Np sub-clouds 225 (Lumia et al., 1983). A minimum distance of 2 m was chosen as it corresponds to the scale where both the amalgamation effect and the over-segmentation are limited (section S2.2.2 in the supplements). Due to the large M3C2 cylinder length pmax required to detect deep landslides, small artefacts with a large distance value, near pmax, occur. To ensure that the smallest detected landslides are not affected by these artefacts, the landslide mapping workflow was performed on the subsample versions of the densest point cloud (section 3.2) with a minimum number of points Np set to 1 during the segmentation step. All detected 230 changes are thus artefacts. The artefact area distribution was then compared with the landslide area distribution using the preand post-earthquake data (details in section 2.2.1 in the supplements). We finally impose a minimum surface of 20 m², corresponding to 20 points, as it corresponds to the smallest area where artefacts are limited. The final dataset represents all the individualized landslide sources and deposit zones detected.

Landslide volume estimation 235
While 3D normal computation is optimal to detect geomorphic changes, it is not suitable for volume estimation which requires to consider normals with parallel directions for a given landslide. As shown in figure 2.a, considering 3D normals can lead to "shadow zones", due to surface roughness, which would result in a biased volume estimate. Therefore, distances and in turn volumes are computed by using a vertical-M3C2 on a grid of core points corresponding to the significant changes ( Fig.2.b).
As the core points are regularly spaced by 1 m, the landslide volume is simply the sum of the vertical-M3C2 distances estimated 240 from the individualized landslides. While the distance uncertainty predicted by the vertical-M3C2 could be used as the volume uncertainty, it significantly overpredicts the true distance uncertainty due to non-optimal normal orientation for the estimation of point cloud roughness on steep slopes (i.e., the roughness is not the detrended roughness). For each landslide source and deposit, we thus compute the volume uncertainty from the sum of the 3D-M3C2 uncertainty measured at each core point, not https://doi.org/10.5194/esurf-2020-73 Preprint. Discussion started: 15 September 2020 c Author(s) 2020. CC BY 4.0 License. the vertical-M3C2 uncertainty. The volume uncertainty is specific to each landslide sources and deposits and depends on the 245 local surface properties such as roughness, the number of point considered and the global registration error, but not on the volume itself.

Geomorphic change and landslide detection inventory
The map of the 3D-M3C2 distances prior to statistically significant change analysis and segmentation highlights erosion (i.e. 250 negative 3D distances) and deposition areas (i.e. positive 3D distances) located on hillslopes and in the river (Fig. 4a). This provide a rare insight into topographic changes following large earthquakes. Most of the detected changes on hillslopes correspond to mass movements such as landslides and rockfalls with variable sizes, shapes and a high portion of them occurred on a previous unstable bare-rock zone in the west part of the study area. Small patch of erosion ranging between 10 m² to 100 m² occur over the entire study area. Their large number illustrate how difficult it would be to manually extract them. Most of 255 the deposit areas are on hillslopes while some the deposits of large landslides have reached the river. The erosion and deposition of sediments by fluvial processes can also be observed in the river.
The area extent of significant changes, where the absolute amplitude of change is greater than LoD95%, represents 17.5 % of the study area (Fig 4.b). Most points associated to stable areas or artificial changes are correctly filtered by the local confidence interval calculation. In particular, points located under the canopy, or associated to a low point density, or wrongly classified 260 as ground while being vegetation points, leading to locally high values of roughness, result into a higher LoD95% and therefore requires a larger topographic change to be detected as significant. After the manual removal of changes in the fluvial domain related to fluvial processes, the minimum 3D-M3C2 distance (or minimum LoD95%) on significant change areas is 0.34 ± 0.33 m and the maximum is 29.9 ± 0.61 m, both corresponding to erosion areas ( Fig.4.b).
The point cloud of significant changes was segmented to identify the landslide sources (i.e. net erosion) and deposits (i.e. net 265 deposition). The final landslide inventory, contains a total of 1431 sources and 853 deposits (Fig. 4c). The difference in the number occurs because many sources share the same deposits that are concentrated at the toe of hillslopes. For sources, the mean absolute vertical-M3C2 distance is 2.11 m, the standard deviation 2.62 m and the maximum absolute value 22.8 ± 1.17 m. For deposits, the mean absolute vertical-M3C2 distance is 2.37 m, the standard deviation 3.08 m and the maximum absolute value maximum 29.9 ± 0.61 m. 270

Landslide area and volume analysis
Recall that the subsequent analysis is based on a grid of core points with 1 m spacing. For each individual landslide, the area A was obtained by computing the number of core points inside the sources region. This represents the vertically projected area, to be consistent with the existing literature based on 2D studies of landslide statistics. The area of detected landslides ranges 280 from 20 to 42,650 m² for sources and from 20 to 33,513 m² for deposits, and the total source and deposit areas are 438,124 and 376,363 m², respectively. The area distribution of landslide sources is computed as follow (Hovius et al., 1997;Malamud et al., 2004): where p(A) is the probability density of a given area within a substantial landslide inventory, NLT is the total number of landslides and A is the landslide source area. δNLT corresponds to the number of landslides with areas between A and A + δA.
The landslide area bin widths δA are equal in logarithmic space.
The area distribution of landslides obeys a power-law scaling relationship consistent with previous studies (e.g., Hovius et al., 1997;Malamud et al., 2004), with an exponent c = -1.79 ± 0.03 (Fig.5a). We do not observe a rollover at small landslide areas 290 which is considered a characteristic of landslide distribution (Guzzetti et al., 2002;Malamud et al., 2004;Malamud and Turcotte, 1999). Varying the parameters of the segmentation does not alter this result (section S2.2.2), nor is a rollover visible at lower surface area if we reduce the minimum landslide size to 10 m² (section S2.2.1). This behaviour differs from the one observed from the landslide area distribution from Massey et al. (2020) in the Kaikoura region.
Landslide volume V was measured with a vertical-M3C2 on landslide areas detected with a 3D-M3C2. The resulting individual 295 landslide volume ranges from 0.75 ± 9.57 m 3 to 171,175 ± 18,593 m 3 for source areas, with a total of 908,055 ± 215,640 m 3 , and from 0.75 ± 17.5 m 3 to 154,599 ± 15,188 m 3 for deposits, with a total of 1,008,626 ± 172,745 m 3 . The uncertainty on total volume estimation represents 23.7 % for sources and 17.1 % for deposits. The volume distribution of the landslides sources was defined using equation (3) replacing A by the volume V and also exhibit a typical negative power-law scaling (Fig. 5.a) of the form: ( ) = . The exponent of the power-law relationship is e = -1.71 ± 0.04. In contrast to landslide source areas 300 distribution, a roll-over is visible on the landslide volume distribution around 5 to 10 m 3 . Considering that the minimum LoD95% in 3D is 0.34 m, and that the minimum landslide area is 20 m², the minimum volume that we can confidently should be 6.8 m 3 , a value consistent with the observed rollover. 35 landslides are indeed smaller than 6.8 m 3 in our inventory. They correspond to peculiar cases of very small landslides in which negative 3D distances close to the LoD95% are positive when measured vertically and thus reduce the apparent volume of eroded material. Volume estimation in these cases should be done 305 locally in 3D, but requires further development not deemed necessary given that it is restricted to the smallest and shallowest landslides of the inventory.
we determine V-A scaling coefficients using two methods: by fitting a linear model (1) on log-transformed data and (2) on 310 log-binned data. While the first method leads to a V-A relationship best describing the volume of each landslide, the second one is not affected by the varying number of landslides in each range of landslide area and leads to a V-A relationship that best matches the total landslide volume. Using the first approach, we find a volume-area scaling exponent of = 1.18 ± 0.01 and an intercept log = −0.42 ± 0.02 m 64 with a determination coefficient 2 = 0.99 (Fig.5.c). Using the second method, we find = 1.20 ± 0.02, an intercept log = −0.41 ± 0.07 m 0.6 and a determination coefficient 2 = 0.99. We also obtain a 315 good correlation R² of 0.86 and 0.82 with the Larsen et al. (2010) relationships derived from soil landslides and from mixed soil landslides and bedrock landslides, respectively (Table 2). However, R² of 0.61 and 0.72 are obtained when considering the parameters of the V-A relationships, derived by Massey et al. (2020), based on the 8,442 cleaned landslides or only soildominated ones, respectively, of the Kaikoura region. At first order, the V-A relationships we obtained are thus consistent with previous studies. Yet, the V-A scaling relationship obtained with log-binned data best predicts the volume directly measured 320 from difference of LiDAR point clouds (Table 2). If the relationships from Larsen et al. (2010) and Massey et al. (2020) were applied to our landslide area inventory, the total volume would vary from 0.541x10 6 m 3 to 1.347x10 6 m 3 (Table 2), compared to 0.908x10 6 ± 0.215x10 6 m 3 that we estimate directly. This highlights the overarching sensitivity of the total volume of eroded material to the V-A relationship (Li et al., 2014;Marc and Hovius, 2015). The closest evaluation of the total volume is based  Table 2. Grey lines are depth vs area relationship for different mean depths.

Reactivated landslides and new failures
Because the 3D measurement approach only depends on local topographic change, we evaluate the fraction of reactivated landslides in the population that would have been hard or impossible to detect with 2D imagery based on texture change 345 ( Fig.6.a). We hereinafter distinguish between new and reactivated landslides, by considering that reactivated landslides occur on bare rock areas in the pre-earthquake imagery or on vegetated areas with limited texture and colour contrast between the two epochs. These definitions were chosen following classical approaches for landslide detection based on vegetation cover analysis (Behling et al., 2014;Marc et al., 2019;Martha et al., 2010;Massey et al., 2018;Pradhan et al., 2016).
Most of the detected landslides are new failures with a total area of 318,726 m² and volume of 636,359 ± 163,496 m 3 (Fig.6.a). 350 We find that reactivated landslides have a total area 119,398 m² and volume of 271,695 ± 52,144 m 3 , which represents 27.2 % and 29.9 ± 12.8 % of the total landslide area and volume, respectively (Table 3). Figure 6b illustrates in detail an area experiencing active rock avalanches and debris flow in the pre-earthquake imagery, for which the change in contrast and shadows is likely too complex to detect a topographic change from a texture based analysis. On the contrary, the 3D change detection shows that landslide erosion is pervasive in this sector, and corresponds indeed, to the largest landslide detected by 355 our approach. The difficulty to detect small and reactivated landslides is illustrated by plotting the version 1.0 of the landslide inventory from Massey et al. (2018) in our study area ( Fig.6.a), which is a database manually validated, and constantly updated.
Only 27 landslide sources were initially mapped and none in reactivated zones, where we found 1431 landslides for which 27.2% of the total area are reactivated.   Table 3: Area and associated volume of the considered new failure and reactivation zones.

Discussion
The aim of this paper is to investigate the potential of methods based on 3D point cloud differencing to provide a landslide inventory map at a region scale from LiDAR data, a total landslide volume estimate and to overcome issues such as landslide 370 amalgamation effects on total estimated landslide volume, under-detection of reactivated landslides in 2D imagery analysis as well as limitations of the DoD approach on steep slopes. Here we first discuss the 3D workflow we have developed in comparison to traditional DoD approach, and then discuss the benefits of 3D change detection for landslide inventories creation compared to 2D imagery landslide detection.

3D point cloud differencing and landslide detection 375
3D point cloud differencing methods have already been applied in previous studies to detect geomorphic changes on single landslides using point clouds obtained by photogrammetry with drone-based images (Esposito et al., 2017;Stumpf et al., 2014Stumpf et al., , 2015. Larger-scale approaches have been attempted to evaluate the impact of a given landslide (Bossi et al., 2015). However, the DoD, based on gridded data, remains the dominant approach to evaluate topographic changes in the context of glacier dynamics, fluvial dynamics or tectonic deformation analysis (e.g. Passalacqua et al., 2015 for a review). To our knowledge, 380 the systematic detection and segmentation of hundreds of landslides from 3D point cloud have not yet been attempted.

Vertical versus 3D change detection capability, and the M3C2 algorithm
The importance of detecting change in 3D, as opposed to vertically in steep slopes can be illustrated by a simple exercise, similar to section 3.2, in which we use two random sub-sampled versions of the post-event point cloud that represent exactly the same surface without any registration error, have similar point density but different sampling of the surface. We apply a 385 vertical-M3C2 and a 3D-M3C2 to the two point clouds, and the maps of change and distance distributions are shown on figure   7a. Typical of change measurement methods on rough surfaces with random point sampling (e.g., Lague et al., 2013), a nonnull distance is often measured even though the two point clouds are samples of exactly the same surface. The distribution of measured distances is centred near zero, with a mean of −2 .10 −4 and 10 −4 m, for the vertical and 3D approach respectively. However, the 3D approach results in a standard deviation, σ=0.05 m, four times smaller than using a vertical differencing, σ= 390 0.20 m. The map of distance shows that vertical differencing systematically results in much higher distances on steep slopes than the 3D approach, while they both yield similar, low distances, on horizontal surfaces.
We thus find that 3D point cloud differencing offers a greater sensitivity to detect changes compared to classical vertical DoD.
This difference is particularly important as it propagates into a lower level of detection and uncertainty on volume calculation.
Using the M3C2 algorithm in 3D (Lague et al., 2013) also offers the benefit of accounting for spatially variable point density 395 and roughness in estimating a distance uncertainty for each core point, that can subsequently be used in volume uncertainty calculation. For instance, this approach leads to a reduced detection sensitivity in vegetated areas due to a lower ground point density and potentially higher roughness due to vegetation misclassification. By using a regular grid of core points as in Wagner et al. (2017), our workflow combines the benefits of working directly with the raw unorganized 3D data, as opposed to DoD where the relationship with the underlying higher point density data is lost, while producing a result with regular sampling that 400 can easily be used for unbiased spatial statistics, volume calculation and easy integration into 2D GIS software. Compared to DoD, if an interpolation is needed, it is performed on the results rather than on the original DEM which can lead to uncontrolled error budget management.

Tectonic internal deformation, data quality and point clouds registration
One of the most critical part of any 3D point cloud processing method is the co-registration of the two point clouds, in particular in the context of co-seismic landsliding. With LiDAR data, the registration error will generally set the minimum detectable 410 change on bare planar surfaces. Here, a rigid translation has been applied on the entire area using an ICP algorithm (Besl and McKay, 1992), de facto assuming that internal deformation during the earthquake was negligible. After applying a vertical displacement of 1.36 m, we did not observe a systematic horizontal shift of the difference map either north or south of the Hope fault. We thus conclude that the internal deformation was below the typical registration error in our study area. For larger studied regions with internal deformation and in the absence of a 3D co-seismic deformation model that could be applied to 415 the post-EQ point cloud (e.g., Massey et al., 2020), our workflow should be applied in a piecewise manner with boundaries corresponding to the main identified faults or deformation zones. For landslide inventories following climatic events, the application to very large dataset should be straightforward as no internal deformation is expected. Similarly, we also noted an internal flight line height mismatch of 0.05-0.12 m in the pre-EQ survey, that was difficult to correct after data delivery and generated some apparent large scale low amplitude topographic change (Fig. 3). This highlights the need for detailed quality 420 control (e.g., by applying M3C2 on overlapping lines) to ensure the highest accuracy possible of the LiDAR data.

Landslide segmentation
Another important aspect of the method is the segmentation procedure to individualize sources and deposits. We performed a sensitivity analysis to evaluate the impact of the minimum distance between sub-clouds Dm, on the number of segmented landslides and related geometric characteristics (section 3.2, Suppl. Material S.2.2.2). We show that below 4 m, the choice of 425 Dm has little impact on the pdf of area, pdf of volume, and the volume-area relationship. However, above 4 m, amalgamation starts to significantly alter the landslide geometry statistics. A similar analysis should be performed for each new dataset to evaluate the best segmentation scale. The connected component segmentation is a simple and rapid way to individualize landslides but given the complexity of the 3D dataset, and in particular the very large range of landslide sizes, inevitably exhibits some drawbacks and is subject to improvement. For instance, landslides occurring on both side of a same collapsed 430 divide are considered as one landslide if they are close enough (Fig. 8.a). More advanced segmentation approaches accounting for normal direction, divide organization and 3D depth maps of amalgamated scars would be needed to improve the segmentation and get more robust results on very large datasets. We note however that these issues do not affect the total landslide volume calculation.

Landslide topographic changes and feature tracking 435
Our approach focuses on landslide detection and volume calculation. The workflow we have designed is not suitable for deformation measurement based on feature tracking (Passalacqua et al., 2015). Except for a few landslides with limited displacement in which point cloud features could be potentially tracked, the severity of landsliding and the long runout of many landslides preclude any attempt in tracking features. Our approach may miss translational landslides on planar hillslopes for which topographic change occurs in the direction of the surface normal. Given that hillslopes are generally not perfectly 440 planar, any significant translation parallel to the hillslope will generate a topographic change, especially in the source area. As shown in fig. 8b, our approach detects translational landslides, although it cannot compute the deformation parallel to the hillslopes. The only element that could be easily tracked in this inventory are the barycenter of the source and associated deposit of each landslide, to explore runout dynamics, but we have not investigated this option yet. Over the studied area of ~5 km 2 , 1431 landslide sources and 853 landslide deposits were detected with the 3D point cloud processing workflow. This relatively large number of landslides is mostly associated to the large number of small landslides (< 100 m², n = 977) that were detected thanks to the resolution of the data. The scaling of the pdf of volume of sources, e = -1.71, indicates a slight tendency for the overall eroded volume to be dominated by the largest landslide (171,175 m 3 , that is 18.8 % of the total volume). The uncertainty on total landslide volume, 17.1% to 23.7 % for deposits and sources, respectively, 455 might appear large, but is based on a conservative 95% confidence interval that we use throughout our analysis. It is mostly dominated by the registration error (reg = 0.17 m) and by the lower point cloud density of the pre-earthquake LiDAR data ( are not statistically different. The larger volume of deposit is however consistent with rock decompaction during landsliding.
We also note that debris deposits form more concentrated and thicker patches at the toe of hillslopes, and are thus more 460 systematically above the detection threshold. Very shallow rockfalls might not be detected and accounted for in the source volume. Hence, we expect to detect more of the population of landslide deposits than of the population of sources.

Distribution of landslide volume and area: power-law behaviour
We obtain a range of landslide area over 3 to 4 orders of magnitude (20 to 42,650 m²) for which we constructed the pdf of area and volume. In landslide analysis, the pdf of landslide area represents the basis to estimate large-scale landslide erosion (Larsen 465 et al. 2010). This distribution has generally a negative power law behaviour for landslide areas larger than a given threshold and displays a rollover for smaller landslides (Fan et al. 2019;Guzzetti et al., 2002;Malamud et al., 2004;Malamud and Turcotte, 1999;Stark and Hovius, 2001). In this study, we find that the exponent c of the power-law for the landslide area distribution is -1.79 ± 0.03 (Fig.5a). This is roughly consistent with the exponents obtained over the entire Kaikoura coseismic landslide inventory of -1.88 (NLT = 10,195;Massey et al., 2018) and more recently of -2.10 (NLT = 29,557; Massey et al., 2020). 470 We present here, one of the first landslide volume distribution derived directly from 3D topographic data, rather than inferred from the combination of detection of the landslide area distribution on 2D data and an estimated volume-area relationship which is much more difficult to precisely estimate. Our direct measurements show that the landslide volume distribution indeed obeys a power-law relationship with an exponent e = -1.71 ± 0.04, consistent with exponents estimated in previous studies 1.0 ≤ e ≤ 1.9 and 1.5 ≤ e ≤ 1.9 for soil landslides (Brunetti et al., 2009;Malamud et al., 2004). Further analyses, similar to ours are 475 necessary to get a better handle on the landslide volume distribution, a critical information with respect to risk analysis and landslide erosion calculation.

Rollover in the distribution of landslide area
Several hypothesis have been put forward to explain the rollover behaviour at low landslide area, including the transition to a cohesion dominated regime reducing the likelihood of rupture (Frattini and Crosta, 2013;Jeandet et al., 2019;Stark and 480 Guzzetti, 2009), a cohesion gradient with depth (Frattini and Crosta, 2013), landslide amalgamation (Tanyas et al., 2019) or the under detection of small landslides (Hovius et al., 1997;Stark and Hovius, 2001). The landslide area distribution in our study does not display a rollover. We have proposed a method to rigorously evaluate the likelihood of detecting spurious landslides due to different data sampling of a rough surface at small landslide areas (Fig S.2.2.1). We set a conservative lower bound of 20 m² above which we are confident that all detected changes are true topographic change. We have also checked 485 that the segmentation distance does not impact the absence of a rollover in our data (Fig S2.2.2 in the supplementary material).
Given that Massey et al. (2020) reports a rollover at ~ 50 m² for the Kaikoura earthquake landslide inventory based on 2D optical image analysis ( Fig.5.a), our results supports the idea that the rollover behaviour observed in previous studies is likely caused by an under detection of small landslides, even with high resolution imagery (Hovius et al., 1997;Stark and Hovius, 2001). If this hypothesis is correct, the number of landslides potentially missed in previous studies can be important. If we 490 consider the power law fitting statistics for landslide area distribution from Massey et al. (2020), the number of landslides between 20 m² ≤ A < 500 m² potentially missed would be around 169,000. This represents 92% of total landslides that would not be considered. While the under detection of small landslides would not greatly affect the total landslide volume estimation, it could have consequences for our understanding of natural hazards mechanics. This promising result possibly highlights the advantages of using LiDAR data combined to our 3D differencing workflow with low level of change detection to generate 495 more accurate and complete landslide inventory datasets.

Landslide volume-area relationship
The landslide volume-area (V-A) scaling relationships obtained in this study are close to the one of Larsen et al. (2010) for soil landslides. This is consistent with the fact that 50% of the landslide thicknesses are lower than 1 m, showing that most of our inventory is relevant to shallow landsliding. The V-A scaling relationship of Massey et al. (2020) for soil landslides gives 500 the best estimation of total landslide volume but overestimate the volume of small landslides. The differences in total volume predicted by our two V-A scaling relationships show that estimates of landslide volume deduced from such relationships greatly depend on the method used to fit data. Our results suggest that fitting model on log-binned data gives a better result in total landslide volume estimation. However, measuring the volume directly from topographic data overcome the issue of choosing a peculiar V-A relationship. 505

Toward a limitation of amalgamation and reactivation biases on landslide volume estimation
The amalgamation effect is a classical issue for 2D landslide mapping and volume assessment (e.g., Li et al., 2014;Marc and Hovius, 2015b), which leads to a higher number of large landslides than should be expected, and a significant overestimation https://doi.org/10.5194/esurf-2020-73 Preprint. Discussion started: 15 September 2020 c Author(s) 2020. CC BY 4.0 License. of the total volume of landslides when using a non-linear V-A relationship. A simple solution to this last problem consists in directly measuring the total landslide volume by comparing pre-and post-earthquake topographic data (e.g., LiDAR) with a 510 DoD (Massey et al., 2020) or with the method described in this paper. While our simple segmentation approach cannot resolve the amalgamation of individual landslides perfectly, the total volume associated to sources and deposit can be robustly estimated independantly of a V-A relationship.
The detection of reactivated landslides from 2D optical imagery remain challenging due to to the weak contrast of vegetation between pre-and post-earthquake periods, or even the absence of vegetation. Recent methods have been developped to detect 515 reactivated landslide based on NDVI-trajectory induced by differences in revegetation rate (Behling et al., 2014(Behling et al., , 2016 but the impact of possibly missed reactivated landslide on total volume estimation is still poorly understood (Guzzetti et al., 2009).
2D or 3D topographic differencing methods are insensitive to the lack of texture variation and can resolve the issue of reactivation. Because vegetation barely develops on very steep slopes and cliffs, the 3D differencing approach oriented perpendicular to the local topographic slope, as opposed to vertical differencing, is critical in detecting subtle changes. Hence, 520 our 3D approach detects landslide occurring in steep areas with poor vegetation cover (Fig.6.b) that would have otherwise been missed with 2D optical imagery approaches, or incorrectly detected with vertical differencing. In this study area, the proportion of reactivated landslide area, 27.2%, is lower than new landslide failures (Table.3), and most of the large reactivated landslides that we find have not been included in the initial mapping of 27 landslides by Massey et al. (2018). Assuming that reactivated landslides cannot be detected by classical methods, the volume potentially lacking represents 29.9 ± 12.8% of the 525 total estimated landslide volume, including some of the largest failures of the studied area. Our study area was chosen based on LiDAR data availability, and may contain a particularly high proportion of reactivated landslides due to the presence of actively eroding bare bedrock hillslopes compared to the size of the study area. We expect this proportion to significantly vary when considering other landscapes with potentially varying proportion of vegetation cover, vegetation density and type (e.g. grass, shrubs, trees), lithology and ground shaking intensity. Nonetheless, our finding represents a first approach to the issue 530 of considering reactivated landslides in total landslide volume estimates, and our results indicates that this should not be neglected at least in regions dominated by a low or absent vegetation cover.
To evaluate the difference of volume that would have been estimated from traditional methods impacted by under-detection of reactivated landslides, we apply the Massey et al. (2020) V-A relationship for all landslides only to the new failures detected in our inventory. As our inventory has very likely a much lower detection level than optical methods (see 5.2.3), we only 535 consider new failures with a minimum area of 50 m² (typical of the rollover observed in Massey et al. 2020). This amount to a 2D traditional processing of the study area. This traditional approach would predict a total volume of 616,308 m 3 , significantly lower than the volume measured directly by our approach. This highlights the need to more systematically generalize the use of 3D data to improve the creation of robust landslide inventories and generate more accurate estimate of the total volume of sediment produced by earthquakes or climatic events. 540

Conclusion
In this paper, we introduced a new workflow for semi-automated landslide detection and volume estimation using 3D differencing based on high resolution topographic point cloud data. This method uses the M3C2 algorithm developed by Lague et al. (2013) for accurate change detection based on the 3D distance normal to the local surface as well as a vertical-M3C2 for volume calculation of landslide sources and deposits once a 3D connected component segmentation procedure has been applied 545 to individualize landslides. Spatially variable uncertainties on distance and volume are provided by the calculation and used in the workflow to evaluate if a change is statistically significant or not, and for volume uncertainty estimation. We provide various tests and recipes to estimate the registration error and choose the parameters of the M3C2 algorithm as function of the point cloud density to ensure the lowest level of change detection, and the best resolution of the 3D map of change. Applied to a 5 km² area located in the Kaikoura region in New Zealand with pre-and post-earthquake LiDAR, we showed that: 550  A level of 3D change detection at 95% confidence of 0.34 m can be reached with airborne LiDAR data, which is largely set by the registration error. Because it operates on raw data, M3C2 accounts for sub-pixel characteristics such as point density and roughness that are not accounted for when working on DEMs, and results in more robust statistics when it comes to evaluate if a change is significant or not. 3D point cloud differencing is critical on steep slopes, and allows to decrease the level of change detection compared to the traditional DoD. 555  Adding elevation information in landslide detection removes the amalgamation effect on the total landslide volume estimation by directly measuring it in 3D rather than considering an ad hoc V-A relationship. Amalgamation in 3D is still a potential issue when exploring individual landslide area and volume statistics given the simplistic segmentation approach that we have used, however our approach has the benefits of more systematically capturing small landslides than traditional approaches based on 2D imagery with manual or automatic landslide mapping. 560  Landslide reactivation on surfaces lacking a significant vegetation cover is classically missed with 2D imagery processing due to the lack of texture or spectral change. 3D processing fully accounts for reactivated landslides. In our study area, 29.9 ± 12.8 % of the total volume was due to landslide reactivation, highlighting that in areas with a mixture of vegetated and non-vegetated steep slopes, 2D approaches can significantly underestimate the number and volume of landslides. 565  As this method provides direct 3D measurement, landslide geometry properties can be explored and tested such as the V-A relationship, landslide area and volume distribution and others. Our results are broadly consistent with the V-A relationship scaling parameters determined by Larsen et al. (2010) and Massey et al. (2020) for soil landslides, with a scaling exponent of 1.20. The largest and deepest landslides deviate significantly from this trend, but they are too few in our database to confidently infer a scaling relationship for these. 570  No rollover is observed in the landslide area distribution down to 20 m², our conservative resolution limit. Inventories based on 2D images analysis following the Kaikoura EQ typically observe a rollover at 50 m² (Massey et al., 2020). This lend credit to the hypothesis that the rollover systematically observed in landslide area distributions generated from 2D images is related to an under detection of small landslide.
Our 3D processing workflow is a first step towards harnessing the full potential of repeat 3D high resolution topographic 575 surveys to automatically create complete and accurate landslide inventories that are critically needed to improve landslide science and managing the cascade of hazards following large earthquakes or storm events. Current bottlenecks to apply this workflow over larger scales, beyond the availability of 3D data itself, are the registration of pre-and post-EQ data when complex co-seismic deformation patterns occur, and limitations of the segmentation method in high landslide density areas.
While airborne LiDAR is best suited to vegetated environments and currently results in the best precision compared to aerial 580 or spatial photogrammetry, the workflow operates for any kind of 3D data.

Code availability
The code producing the landslide inventory in this study is available as a jupyter notebook form at https://github.com/Thomas-Brd/3D_landslide_detection, and is also archived in Zenodo: http://doi.org/10.5281/zenodo.4010806.

Data availability 585
LiDAR data used in this study can be found at: http://doi.org/10.5281/zenodo.4011629 The final landslide source and deposit information supporting the findings of this paper can be found at https://github.com/Thomas-Brd/3D_landslide_detection, or in Zenodo: http://doi.org/10.5281/zenodo.4010806