Volume Characteristics of Landslides Triggered by the MW 7.8 2016 Kaikōura Earthquake, New Zealand, Derived From Digital Surface Difference Modeling

We use a mapped landslide inventory coupled with a 2‐m resolution vertical difference model covering an area of 6,875 km2 to accurately constrain landslide volume‐area relationships. We use the difference model to calculate the source volumes for landslides triggered by the MW 7.8 Kaikōura, New Zealand, earthquake of 14 November 2016. Of the 29,519 mapped landslides in the inventory, 28,394 are within the analysis area, and of these, we have calculated the volume of 17,256 source areas that are ≥90% free of debris. Of the 28,394 landslides, about 80% are classified as soil or rock avalanches and the remainder as mainly translational slides. Our results show that both the soil avalanches and the rock avalanches, ignoring their source geology, have area to volume power‐law scaling exponents (γ) of 0.921 to 1.060 and 1.040 to 1.138, respectively. These are lower than the γ values of 1.1–1.3 (for soil) and 1.3–1.6 (for rock) reported in the literature for undifferentiated landslide types. They are, however, similar to those γ values estimated from other coseismic landslide inventories. In contrast, for 50 selected rotational, translational (planar slide surfaces), or compound slides, where much of the debris remains in the source area, we found γ values range between 1.46 and 1.47, indicating that their slide surfaces were considerably deeper than those landslides classified as avalanches. This study, like previous studies on coseismic landslides, shows that soil and rock avalanches (disrupted landslides) are the dominant landslide type triggered by earthquakes and that they tend to be shallow.


Introduction
Quantifying rates of landslide erosion is fundamental to understanding how, and over what time scales, landslide debris generated by earthquakes and/or precipitation events is transported from hillslopes to the oceans. To do this, it requires accurate estimation of the location, type (i.e., failure mechanisms and movement processes), and volume of the initial landslides. Volume estimates are difficult to determine for coseismic landslides given: (1) the size of the area typically impacted; (2) the large number of landslides produced by a major earthquake, especially in mountainous regions; (3) the general lack of high-resolution pre-earthquake topographic models (and sometimes optical images) to compare with data sets captured after a major earthquake; and (4) the three-dimensional nature of the landslide slide surface, which may be obscured by debris and/or vegetation.
One approach to this challenge is to map the landslide characteristics in the field. However, for events generating many landslides, estimating the source volume for each landslide using field-based measurement techniques is unrealistically time-consuming. Such measurements are likely to be inaccurate because it is difficult to determine the original ground surface and landslide depth across the source area in the field, given the large surface changes involved. Therefore, using high-resolution pre-and post-landslide digital surface models (DSMs) to estimate landslide source volumes is a more promising method by which to produce efficient and accurate data sets. However, this is frequently constrained by the availability of high resolution pre-and post-landslide digital surface data.
In the past, many studies have relied on representing the landslide head scarp/centroid, length of scar (debris source and trail combined) and scar area as points, lines, and polygons respectively, using manual and/or automated mapping processes. However, such studies often do not differentiate between the landslide source area from the debris trail (e.g., Parker et al., 2011;Kargel et al., 2016), adding uncertainty.
Widely used statistical models that relate the landslide volume to the landslide scar plan area are typically constructed based on a limited number of data sets with variation in both the methods; they use and data quality (Table 1). In some cases, these relationships are based on measurements where pre-and post-earthquake topographic models are available and/or where geometrical and field measurements of landslide area and depth can be made (e.g., Xu et al., 2016). However, many older coseismic landslide studies are not based on such data; therefore, the resulting estimates of landslide volume can have considerable uncertainty that propagates through into subsequent analyses (e.g., Marc et al., 2016). This uncertainty also hinders our ability to investigate what factors might be controlling the volume of the landslides generated by a specific significant earthquake or precipitation event.
Global data sets that relate landslide volume to area have been produced by two main studies: Guzzetti et al. (2009) and Larsen et al. (2010); both mainly contain precipitation-induced landslides or landslides where the trigger mechanism is unknown. Guzzetti et al. (2009) compiled a data set of 677 landslides and derived an empirical power law relationship between landslide volume (V) and area (A), based on geometrical measurements of individual landslides: where α is the multiplier and γ is the scaling exponent derived from fitting a power law relationship to the compiled data. The power law relationship derived by Guzzetti et al. (2009) has a multiplier of α = 0.074 and a scaling exponent of γ = 1.450 (with a standard error of ±0.0086). They concluded that the relationship is largely geometrical, and not influenced significantly by geomorphological or mechanical properties of the soils or rocks, or the landslide types. Larsen et al. (2010) compiled a much larger data set of 4,231 landslides with geometrically derived measurements of total landslide volume and area, and their source materials are classified into soil or bedrock. They found that (1) shallow, soil-based landslides have a  (Larsen et al., 2010). Hovius et al. (1997) present an area to volume scaling exponent of γ = 1.5 for landslides in the western Southern Alps of New Zealand. Working in the same area, Korup (2005) presents a landslide area to volume scaling exponent of γ = 1.98 based on 23 large landslides; however, many of these were creeping landslides involving entire hillsides that are inferred to have been moving (and possibly still are) for hundreds if not thousands of years. Whitehouse (1983) presented an area to volume scaling exponent of γ = 1.25 for rock avalanches in the central Southern Alps. These data sets are likely to contain both earthquake-and rainfall-induced landslides (EIL and RIL). Hancox et al. (2014Hancox et al. ( , 2016 report scaling exponents of γ = 1.10-1.23 and 1.10-1.15 for those landslides-predominantly soil and rock (debris) avalanches-triggered by the M W 7.8 1929 Murchison Earthquake (Buller) and the M W 7.1 1968 Inangahua Earthquake, respectively, and estimated from field measurements of landslide depth, which would suggest that EIL have lower γ values than landslides initiated by non-earthquake triggers, possibly because such data sets are dominated by soil and rock (debris) avalanches (Hungr et al., 2014) also referred to as disrupted landslides (Keefer, 2002).
scaling exponent γ from 1.1 to 1.3; and (2) landslides that involved the failure of bedrock were deeper and hence had a larger volume for a given surface area, with γ in the range 1.3-1.6. The smaller exponent values for landslides in soil may occur because their volume is limited by soil thickness. These variations in the value of γ imply that to produce accurate estimates of landslide volume from aerial measurement, it is necessary to distinguish between soil and bedrock landslides. This is not always possible for coseismic landslides. An uncertainty associated with the landslide data set compiled by Guzzetti et al. (2009) is that they adopt total landslide areas and do not distinguish between the landslide source and debris trail. As a result, they may overestimate γ, which is why in this study, we use the landslide source area to estimate volume. These compilations also do not distinguish between earthquake and precipitation trigged landslides, and many do not differentiate γ values for distinct types of landslide, thus making it impossible to determine whether the γ value is affected by landslide trigger and type.
In this paper, we compare landslide volume-area relationships for landslides triggered by the 2016 M W 7.8 Kaikōura Earthquake derived from regional-scale, pre-and post-event surface difference modeling with those presented in other data sets listed in Table 1. Using these data, we investigate the influence that the landslide source material, geological structure, and landslide type have on the value of γ. From these results, we make some recommendations about how to use such global landslide volume-area relationships that consider landslide type and source material.

Kaikōura Earthquake Landslide Inventory Version 2.0 and Updates
In this study, we use Version 2.0 of the digital inventory of landslides triggered by the M W 7.8 Kaikōura, New Zealand, earthquake of 14 November 2016-which contains 29,557 hand-digitized landslides ( Figure 1). Version 1.0, published by Massey et al. (2018), included 10,195 mainly larger (>10,000 m 2 ) landslides, including some mapped from satellite imagery. Version 2.0 contains all landslides from Version 1.0-which were modified where necessary-plus additional landslides (mainly <10,000 m 2 ) mapped from 0.3 m resolution pre-and post-earthquake digital optical imagery. While Version 2.0 contains significantly more landslides than the Version 1.0 inventory, it does not contain all the smaller landslides, typically <500 m 2 that were triggered by the earthquake, as differentiating such small landslides from apparent "noise," e.g., shadows, vegetation changes, and other landslides, was difficult and time-consuming. The methodology adopted for mapping the landslides is detailed in Dellow et al. (2017) and Massey et al. (2018) and follows the proposed "criteria" listed in Fan et al. (2019). Field estimates of the mapped landslide polygon boundaries suggest a nominal positional error of 1-5_x00A0;m. However, this accuracy reduces for those landslides where debris remains within the source, as delineating the overlap between the source and debris boundaries is difficult, relying on the experience of the mapper and the quality of the shade and difference models derived from the DSMs and light detecting and ranging (Lidar) digital elevation models (DEMs), and the optical fixed wing and helicopter imagery used.
The earthquake-triggered landslides cover an area of 10,000 km 2 , with the majority concentrated in a smaller area of about 3,600 km 2 . Using the landslide sources contained in the Version 2.0 inventory and the methods described in Massey et al. (2018), we have recalculated (1) the landslide frequency and area and (2) the landslide area density, defined by where p(A L ) is the probability density of a given area within a near complete inventory-defined as the frequency density of landslides of a given source area bin (A L ), divided by the total number of landslides in the inventory. N LT is the total number of landslides in the inventory, and δN L is the number of landslides with areas between A L and A L + δA L . For the landslide area bins, we adopted bin widths (δA L ) that increased with increasing landslide source area (A L ), so that bin widths were equal in logarithmic space.
All landslides in the Version 2.0 inventory were hand digitized using the highest available resolution aerial photography (Aerial Surveys, 2017), captured between December 2016 and January 2017 (post-earthquake). During this period no significant rain events occurred, indicating the earthquake was the main landslide- triggering mechanism. Pre-earthquake aerial photography from January 2015 was used to ensure that the landslides visible in the post-earthquake imagery were induced by the earthquake, although there is a small chance that some may have been triggered either by aftershocks or non-earthquake triggers such as rain between the pre-and post-earthquake surveys. Both epochs of aerial imagery were captured at 0.3-m ground sample distances, defined as the distance between pixel centers measured on the ground. The original deformation models have a ground resolution of 90 m and were resampled to 2-m resolution using a natural neighbor interpolation. The large black star is the epicenter of the 14 November 2016 main earthquake. The smaller black stars are aftershocks ≥magnitude M6 that occurred in the two days after the main earthquake.

Digital Surface Models
DSMs were generated from each epoch of digital stereo aerial photographs, covering the main area affected by landslides (Figure 2), using the semi-global matching correlation methodology in Image Station Automatic Elevations-Extended software. Aerial Surveys (2017)

Estimating the Uncertainties in the Difference Models
Bootstrapping (Efron & Tibshirani, 1993) was used to estimate the statistical variance in the difference models that were derived from (1) subtracting the 2015 InSAR-adjusted DSM and (2) the 2015 unadjusted DSM, from the 2017 DSM. To do this, we randomly sampled the two difference models to derive sample data sets of 1,000 points attributed with the vertical change between surveys. Areas of the difference model representing holes in the data, landslides, and water bodies (>10,000 m 2 ) were excluded from the sampling. We repeated the sampling 999 times, to create 1,000 samples of 1,000 randomly selected points (i.e., the "bootstrap model") from both difference models ( Figure 3). While this approach allowed the regional tectonic displacement to be effectively removed from the difference model, a regional systematic vertical bias ("offset") of −0.45 m in the 2015 DSM remained after the adjustment. This vertical bias was "corrected" from the 2015 InSAR-adjusted DSM by lowering each point systematically by 0.45 m. The 2015 InSAR-adjusted and now corrected DSM was then subtracted from the 2017 DSM to generate a difference model. As a test, the bootstrapping was redone using this difference model, to calculate the point frequency distribution of the vertical elevation differences and the bootstrap mean and standard deviations ( Figure 3). The results show that the bootstrap mean is 0.00 (±1.83) m, error at one standard deviation (1σ).
To further explore the uncertainties in the DSM difference models, we compared the 2015 and 2017 DSMs with the DSMs and DEMs derived from the airborne Lidar surveys carried out in 2012 (pre-earthquake) and 2017 (post-earthquake), where they overlapped ( Figure 2). To do this, we compared the following elevation models (Lidar was resampled to 2 m):  (1) within the area of the InSAR-adjusted and -corrected difference model; (2) have source areas that are ≥90% free of debris, determined using the source and debris trail polygons; and (3) have source area length to area ratios of <1.2, as sources with ratios ≥1.2 represent wide but relatively short sources, which in most cases comprise multiple sources where the boundaries between each source cannot be clearly differentiated. Only these 17,256 landslide sources were considered appropriate for volume estimations. Volumes of these landslide sources were calculated using ArcGIS by summing the 2-m grid cells of the difference model within the mapped landslide source-area polygons to estimate the volume of change, where negative changes are assumed to imply erosion and positive change to imply deposition ( Figure 4). We used the hand-digitized landslide source polygons because the uncertainties in the DSM-caused mainly by vegetation effects-mean that accurate automated mapping of the source polygons from the change model could not be relied upon. The landslide source areas were rasterized to calculate the volumes from the difference model creating up to a 1-m buffer around the source-area polygons. Only negative changes were used to estimate volumes. To be consistent with the past literature cited in this paper, we use only the two-dimensional plan area of the landslide sources. Landslides with source areas smaller than the resolution of the difference model (2 m) were not included in this assessment.
We manually calculated the volumes of 50 selected translational (planar slide surfaces), rotational (arcuate slide surfaces), or compound (where slide surfaces are a mixture of planar and arcuate) landslides in the upper Cretaceous to Neogene limestones, siltstones, and sandstones-where the debris remained relatively intact and within the source area, thus obscuring the slide surface. Such landslides would be classified as "coherent" following the classification scheme of Keefer (2002). These calculations were undertaken by compiling cross sections through the landslides and comparing the differences in topography before and after the earthquake and the morphology of the post-earthquake debris, allowing their failure surfaces to be estimated (e.g., Jaboyedoff et al., 2015). Of these 50 landslides, 25 were visited in the field to validate the mapping and volume estimates. The field validation was challenging because the slide surfaces were obscured by debris, and so we used the landslide debris extent and morphology of the scar, tension cracks, and associated block displacements to project potential slide surfaces in three dimensions to estimate the volumes of the landslides.
For the 17,256 landslides considered in our analysis, general linear models adopting the general and weighted least squares method were fitted to the log of the source volume and area for each landslide data set, comprising (1) all 17,256 rock and debris avalanches; (2) the debris avalanches subdivided per GeolCodes 1-4 (Table 2), where GeolCodes 1-3 mainly comprise surficial materials derived from fluvial and/or mass movement processes and the in situ weathering of sandstones, limestones, and siltstone, and GeolCode 4, which mainly comprises weathered greywacke rocks; and (3) the 50 translational, rotational, or compound slides in predominantly rock.
We investigated the impact of varying the landslide source volumes as well as weighting them, on the scaling exponent (γ) of the fitted models. Weighting, based on the landslide volume, was used to adjust the contribution of individual landslides to the linear-least squares model fitting. This was done by weighting the source volumes in inverse proportion to the estimated uncertainty of their volume. To do this, we used the standard deviation of the source volume estimated for each landslide. The uncertainties in the volumes of erosion calculated for each landslide source area ≥90% free of debris were quantified by estimating the standard deviation of the erosion for a given landslide source (LSx): where SD BS is the standard deviation of the bootstrap mean, which is 1.9 m. LS N is the number of grid cells within each landslide source area multiplied by the plan area of an individual 2 m by 2 m grid cell (4 m 2 ), which is the ground resolution of the DSM. The weighting for each landslide used in the leastsquares line fitting was calculated in two steps: Step 1: calculate the variance of the logarithm of the volume erosion for each landslide source (V E:LSx ): where LS E is the volume of the landslide source (x) derived from the difference model (  Note: The geological unit descriptions are after Rattenbury et al. (2006) and Heron (2018). The bolds are preferred statistical models fitted to the data, which could be used by others to forecast landslide volumes from source areas.
Step 2: estimate the weighting for each landslide (x) by A weighted linear regression of the log (volume of erosion) on log (landslide area) was then calculated. This has two fitted parameters: the intercept, which is the logarithm of the multiplier α of equation 1; and the coefficient of log(area), which is the scaling exponent γ. The residuals (e) were calculated for each landslide by comparing the logarithms of the observed volume (Vo LSx ) and the predicted volume (Vp LSx ) calculated from the linear least squares models fitted to the data, adopting uncorrected and corrected landslide source volumes, with and without weighting:

Sensitivity Analysis
The sensitivity of the power-laws fitted to all the 17,256 landslides were investigated by identifying and then removing "cleaning" those landslide source areas where the calculated volumes may be affected by shadows, vegetation, and/or mapping effects that propagate through into the digital surface and difference models and thus the volume calculations. This was done in two steps.
Step 1  The landslide source areas used in the analyses are ≥90% free of debris and thus should fall within the bare ground bands based on the 2017 NDVI. This is true for 76% of the 17,256 sources, while the remaining 24% fall within the low to high vegetation bands. A selection of those sources within the low to high vegetation bands were checked using the 2017 imagery. It was found that they had relatively small source areas and were in areas of moderate to high vegetation making accurate mapping of their boundaries difficult.
Step 2 involved screening only those landslide source areas within NDVI bands 3-7 (−0.06 ≤ NDVI ≤ 0.4) in both the 2015 and 2017 epochs, thus resulting in a total of 11,162 "cleaned" landslide sources. These were then further analyzed by grouping them into bands based on their erosion area to source area ratios (erosion cover ratio or %). These ratios were calculated by dividing the area of erosion within the source by the area of the landslide source, using the InSAR-adjusted and -corrected difference model to derive the area of erosion within each landslide source. Ten equally spaced bands were used (0.1/10%) from 0 (0%) to 1.0 (100%). Power laws were then fitted sequentially to the source volume and area of those landslides within each erosion 10.1029/2019JF005163

Journal of Geophysical Research: Earth Surface
cover ratio: (a) band and (b) range. An erosion cover ratio of 0.5 (50%) was used to generate the final "cleaned" sources shown in Figures 5 and 9, representing 8,442 landslides, which were then subsequently subdivided per GeolCode.

Landslide Area and Volume
As expected, the landslide source areas generated by this event exhibit characteristic power-law scaling (Figures 5a and 5b) (e.g., Hovius et al., 1997;Guzzetti et al., 2009;Malamud et al., 2004;Parker et al., 2015). The power-law fitting statistics are N LT of 29,557, x min of 500 m 2 , and α of 2.10, and are similar to the original statistics based on the Version 1.0 inventory (Massey et al., 2018). About 80% of landslides triggered by the Kaikōura Earthquake in the Version 2.0 inventory can be classified as soil or rock (debris) avalanches under the classification scheme of Hungr et al. (2014), (Figure 6), which would be classified as disrupted landslides under the scheme of Keefer (2002). In a global study, Keefer (2002) found that disrupted landslides (debris avalanches) comprise about 86% of all coseismic landslides.
Of the 29,519 mapped landslides within Version 2.0 of the Kaikōura landslide inventory, 70.4% occurred within the Pahau terrane greywacke (GeolCode 4, Table 2). These mainly comprised slides, falls, and topples where the debris quickly broke down to form predominantly soil and rock (debris) avalanches and in some cases flows (Hungr et al., 2014) (Figures 6a and 6b). From field observations, their failure surfaces were either shallow, located along the boundary between the overlying completely/highly weathered rock and residual soil, talus, or colluvium ("regolith") and the underlying less weathered in situ greywacke rock or deeper within the less weathered rock. Of the 29,519 mapped landslides, 11.9% and 8.3% occurred in the upper Cretaceous (GeolCode 3) and Neogene (GeolCode 2) materials, respectively (Table 2; Figures 6c, 6e, and 6f), which are completely to highly weathered sandstones, siltstones, limestones, and associated residual soils and overlying colluvium and terrace sands, silts, and gravels. Of these (GeolCodes 2 and 3 combined), 43.8% were mainly mixed soil and rock (debris) avalanches, but most were landslides that comprised either  Figure 6c). In these materials, the rock masses tend to be more massive with highly persistent bedding planes and clay seams, which allow the development of "deeper-seated" failure surfaces, when

10.1029/2019JF005163
Journal of Geophysical Research: Earth Surface compared to those of a similar landslide type in the highly deformed and closely jointed greywacke. Substantial numbers of pre-Kaikōura Earthquake, large translational, rotational and compound slides were also mapped in these materials of which many reactivated, moving a few centimeters to a few meters during the earthquake. Of the 29,519 mapped landslides, 8.8% occurred in Quaternary soils (GeolCode 1), which typically comprise sand, silt, and gravel originally deposited as river terraces, subsequently uplifted and incised by streams and rivers. The landslides within these materials tend to occur on the steeper "incised" slopes and could be classified as soil avalanches.
Of the 29,519 mapped landslides, 28,394 had source areas that were fully within the area of the difference model (Figure 2a). Of these, 17,256 had source areas that were ≥90% free of debris and had source area length to area ratios of <1.2. These were predominantly soil avalanches or rock avalanches or a mixture of soil and rock (debris) avalanches. Of these 17,256 landslides, 7% were in Quaternary soils (GeolCode 1); 5% were in Neogene materials (GeolCode 2); 9% were in upper Cretaceous and Paleogene materials (GeolCode 3); and 79% were in lower Cretaceous Torlesse (Pahau) terrane (GeolCode 4; Table 2). The remaining 11,168 landslides within the extent of the difference model comprise a combination of soil and rock (debris) avalanche flows and translational, rotational, or compound slides and flows but where the debris remains in the source area, thus obscuring the slide surface and making volume estimation difficult.
The results showing the linear least-squares models fitted to all 17,256 of the landslide source volumes derived from the InSAR-adjusted and -corrected difference models are presented in Figure 7, with and without weighting, and their associated residuals (Equation 6) calculated from each of the fitted models are shown in Figure 8. The model that best fits the data is the one based on the InSAR-adjusted and -corrected difference model with weighting, as the residuals (Figure 8d) are normally distributed, the mode is centered either side of zero, and the standard deviation is the lowest of the four fitted models. The residuals from the models fitted to the landslide volumes estimated from the other difference models show a more lognormal distribution with larger standard deviations.
Our results show that the 17,256 predominantly debris avalanches, irrespective of their source geology and materials, have scaling exponents γ of 1.08 (±0.01 at 1σ) with no weighting and γ of 1.016 (±0.003 at 1σ) with weighting, when using source volumes estimated from the InSAR-adjusted and -corrected difference model and adopting general linear regression models fitted to the logged data (Table 2 and Figure 7). The range of γ from the model fitting to the landslide source volumes estimated from the InSAR-adjusted, but uncorrected, difference model, assuming weighting and no weighting are between γ of 1.052 and 1.022, respectively. These γ values are <1.5, which is the threshold corresponding to the self-similar behavior, when landslides have the same three-dimensional scaling of geometry, regardless of size (Klar et al., 2011;Wartman et al., 2013). These values, indicating relatively shallow landslide sources, are consistent with the observations of landslide geometries made in the field during the initial post-earthquake reconnaissance (Dellow et al., 2017) and more recent fieldwork by the authors as part of the field validation for this work (Figure 4).
The values of γ ( Figure 9, Table 2), are lower for those in GeolCodes 1-3 soils, derived from sandstones, limestones, and siltstones (min = 0.822 and max = 1.013) than those in GeolCode 4 greywacke sandstones (min = 1.040, max = 1.136), indicating that avalanches in predominantly soil tend to be smaller volume than those in predominantly rock.

Sensitivity Analysis
The volumes and areas of the cleaned landslide source areas derived from the sensitivity analyses are shown in Figures 7 and 9 and the fit statistics in Table 2. Supporting graphs and tables are presented in Figures S1-S7 and Tables S1-S3 in the supporting information that accompanies this paper. The cleaning of those landslides based on the NDVI and erosion cover ratio (Figures 7 and 9) show that it is the smaller volume landslides that have been preferentially removed because of the cleaning. These represent those landslides in Figure 4 showing little or no erosion. The power law fit statistics derived from the all the cleaned source areas are similar to those derived from the weighted fits to the uncleaned sources (Figure 7). This is also the case for those power laws fitted to the cleaned and uncleaned data per GeolCodes 1-4 ( Figure 9). However, the intercepts (α values) tend to be lower for those power laws fitted to the cleaned data, indicating that the weighting is not as influenced by the smaller volume landslides.

Journal of Geophysical Research: Earth Surface
NDVI: The results shown in the supplementary section highlight the importance of removing those sources with NDVI values correlated with shadows. The power laws fitted to the sources per NDVI range indicate that the intercepts (α values) increase and the gradients (γ values) decrease with increasing NDVI (vegetation height). Overall, the results show that the total volume of debris produced by the 17,256 debris avalanches if estimated using the power-law fit parameters (fitted to the cleaned source data) based on the NDVI bands would be overestimated for those power laws fitted to those avalanches in the lower bands (<−0.06), representing shadows and partial shadows. This is due to the shadows affecting the elevation of the DSM and the resulting difference model.
Erosion cover ratio: The results shown in the supplementary section indicate that the power laws fitted to the sources per erosion cover ratio range have little influence on the gradient (γ values) but a larger influence on the intercept (α values). The coefficients of determination (R 2 ) also increase with increasing erosion cover

Journal of Geophysical Research: Earth Surface
ratios. This is because at higher erosion cover ratios, a number of landslides with small volumes are reduced thus increasing the intercept and reducing the variance.

Discussion
When compared to the γ values reported in the literature for landslides in "soil" and "bedrock," the γ values of the models fitted to the Kaikōura soil and rock avalanches in the different geological units all plot at the lower end of the range reported by Larsen et al. (2010) for soil landslides (Figure 10a), regardless of the statistical variance resulting from the weighting or cleaning used to calculate the power-law fit parameters (α and γ) in this study.
Field observations indicate that many of the avalanche source areas in soils, sandstones, limestones, and siltstones (GeolCodes 1-3) were confined to the superficial regolith (soil) with very few examples in which the failure surfaces extended deeply into the underlying rock. The low γ values of 0.822-1.091-essentially 1.0 -also indicate that the depth of the landslide source remains relatively shallow and uniform, irrespective of how big the source area becomes. For those in GeolCode 4 (greywacke sandstones), field observations indicate that some source areas occurred in the regolith-comprising residual soil, colluvium, talus, and completely weathered rock-but most were within weathered rock. The γ values of these varied between 1.040 and 1.138, indicating that the depth of the source increases with increasing area. Field observations showed that as these landslides increased in volume, their slide surfaces transitioned from being within the shallower regolith to deeper into the underlying, highly jointed and less weathered rock. This made it difficult to

10.1029/2019JF005163
Journal of Geophysical Research: Earth Surface determine whether soil or rock was the dominant source material, as they typically comprised both. For those sources predominantly in rock, their depth appears to have been controlled by multiple intersecting and closely spaced joint blocks within the weathered greywacke. The greywacke rock masses typically comprise four or more joint sets (Read et al., 2000;Richards & Read, 2007), and therefore, the landslide failure surfaces have multiple discontinuities ("degrees of freedom") along which to propagate, without fracturing the intact rock. Near the surface, the rock masses weather and dilate along these joints, but at depth and therefore at higher normal stresses, the joints are tight and interlocked, increasing the effective strength of the rock mass. This appears to have limited the development of deeper seated failure surfaces Figure 9. Debris avalanche volume versus landslide plan area for the 17,256 landslides (red hollow circles) and the 8,442 "cleaned" landslides (blue solid dots) subdivided per source geology, where (a) GeolCode 1 is Quaternary soils of sand silt and gravel; (b) GeolCode 2 is Neogene completely to highly weathered sandstones, siltstones, and residual soils with minor colluvium; (c) GeolCode 3 is upper Cretaceous and Paleogene completely to highly weathered limestones, sandstones, siltstones, and residual soils with minor colluvium; and (d) GeolCode 4 is lower Cretaceous Torlesse (Pahau) terrane completely to highly weathered greywacke rock and overlying residual soil, talus, and colluvium. Volumes were estimated from the InSAR-adjusted and -corrected difference model.

10.1029/2019JF005163
Journal of Geophysical Research: Earth Surface Figure 10. (a) Landslide area to volume scaling exponents (γ) from the studies listed in Table 1 compared with the results from this study. The landslide inventories are classified based on the dominant landslide source material (soil or rock) and landslide type, where DA = debris (mixed soil and rock) avalanche, RA = rock avalanche, S = translational/rotational/compound slides. The values shown on the plot are the number of landslides used in each study to derive γ. Error bars represent the standard deviation of the γ value, where given. *Hancox et al. (2014*Hancox et al. ( , 2016 do not give standard deviation uncertainty estimates, and so the given ranges in γ values have been used. (b) The total volume of debris produced by the 17,256 debris avalanches estimated from the InSARadjusted and -corrected difference model and the volumes estimated from the given power-law scaling relationships, using the individual source areas of each debris avalanche. The uncertainty range of the total debris avalanche volume estimated from the difference model was calculated by adding and subtracting the bootstrap standard deviation of ±1.8 m (1σ) of the difference model, to the 2015 DSM, as systematic offsets. These new DSMs were then subtracted from the 2017 DSM to create +1σ and −1σ difference models, which were then used to estimate the +1σ and −1σ source volumes for each debris avalanche. The error bars for the volumes estimated using the power-law scaling relationships represent the total volumes estimated using the standard deviations of the powerlaw fits.

10.1029/2019JF005163
Journal of Geophysical Research: Earth Surface when compared with those forming the translational, rotational or compound slides in GeolCode 2 or 3 materials.
For the 50 rotational, translational, or compound slides in sandstones, siltstones, and limestones (GeolCodes 2 and 3), the γ values range between 1.46 and 1.47 based on fitting models to the mean, upper and lower estimates of their volumes (Figure 7b), thus indicating that such landslides have, in general, deeper failure surfaces than the debris avalanches, and field observations suggest these may be structurally controlled. Most of these failures occurred in the upper Cretaceous to Neogene sandstones, siltstones, and limestones, where movement occurs either along bedding planes or other persistent structural discontinuities, such as fault planes, bedding parallel clay seams, or through the rock mass. Even so, the γ values for these landslides are within the center of the range reported in the literature for "bedrock" landslides ( Figure 10 and Table 1).
In summary, the values of γ for the 17,256 predominantly debris avalanches initiated by the Kaikōura Earthquake are low compared with observed global and New Zealand values for mixed bedrock and soil landslides (Larsen et al., 2010; Figure 10a). They are closer to γ values reported by Hancox et al. (2014Hancox et al. ( , 2016 for the (mainly) debris avalanches initiated by the 1929 Murchison Earthquake and 1968 Inangahua Earthquake in New Zealand (Table 1), which occurred in similar materials to GeolCodes 1-3 (Table 2). Thus, like the Kaikōura coseismic landslides inventory, the Inangahua and Murchison coseismic landslide inventories must also contain a high proportion of shallow debris avalanches. The Kaikōura Earthquakeinduced translational/rotational slides in rock have higher values of γ, suggesting that landslide type is influenced by the source geology/material, which thus affects the γ value. There may also be a difference in the γ value between coseismic landslides and precipitation-induced landslides. Most of the landslides reported by Guzzetti et al. (2009) andLarsen et al. (2010) are precipitation induced. The change in stress caused by an earthquake acting on a slope differs from that caused by precipitation-induced increases in porewater pressure; the former involves an increase of shear stress, which increases toward the ridge crests as a result of topographic amplification of ground motion (e.g., Meunier et al., 2008), and the latter involves a reduction of shear resistance, which could occur at depth within a slope. Such differences are known to affect the spatial distribution of landslides (e.g., Meunier et al., 2008) but could also affect the depth of the failure surfaces within a slope.
Slope angle also appears to influence the failure mode (landslide type)-and is intrinsically linked to the geology/material_x2014;and whether the debris can evacuate the source. For example, the dominant landslides in this analysis are those in GeolCode 4 greywacke. The mean slope angle of all source areas in GeolCode 4 that are ≥90% free of debris is 49°(±9°at 1σ), indicating that such landslides occurred on relatively steep slopes. The mean slope angles of the sources in GeolCodes 1-3 (mainly soils) that are ≥90% free of debris are also relatively steep, 40°(±14°). In comparison, the mean slope angles of those sources in GeolCode 4 and GeolCodes 1-3_x2014;where the debris remains within more than 10% of the source area -are lower, 45°(±9°) and 35°(±14°) respectively, and most comprise avalanches in soil. There is no reason to believe that γ value fitted to these source volumes-if they could be accurately measured-would show any significant difference from the γ values derived from the sources that were ≥90% free of debris, given their observed geometry. The mean slope angle of the 50 translational/rotational/compound slides in GeolCodes 2 and 3 is 20°(±10°). Many of the other landslides on these less steep slopes in GeolCodes 2 and 3, where the debris remains in the source, are also translational/rotational/compound slides and not avalanches. These results would suggest that debris from landslides on steeper slopes is more likely to break down and evacuate the source, generating avalanches, while the debris from landslides on less steep slopes does not move as far, thus remaining in the source, generating less mobile avalanches and those landslides on lower angle slopes in GeolCodes 2 and 3 mainly comprise translational//rotational/compound slides, where the debris does not displace far.

Sensitivity
It is possible that the differences between the γ values reported in the literature and in this study could also be related to contrasting method and measurement inaccuracies on estimating volumes and the statistical methods used. For example, most of the volumes reported in Table 1 were not derived from differencing of highresolution pre-and post-failure DEMs, leading to uncertainties in the volume estimates that are not easy to quantify. Other potential issues are (1) how complete the mapped landslide inventories are and how for the same event the number of mapped landslides can vary greatly between events (e.g., the 2010 Haiti Earthquake: Gorum et al., (2013); Harp et al., (2016); and the 2008 Wenchuan Earthquake: Li et al. (2014); Xu et al., 2014); (2) mapping quality and amalgamation of multiple source areas into single larger areas, which can inflate, by up to a factor of three the total estimated volumes (e.g., Li et al., 2014;Marc & Hovius, 2015); and (3) the ability to map or estimate the extent of the landslide source area from the transport/runout zones in order to estimate the true source volume (e.g., Roback et al., (2018); Marc et al., 2019). This study, and those of Hancox et al. (2014Hancox et al. ( , 2016 and Xu et al. (2016), use the assumed landslide source area to define the landslide area. Other studies in Table 1 typically use the scar-the combined source and debris trail areas-which is less representative of the initial landslide source volume. Most other studies use the two-dimensional plan area of the landslide, and so using the three-dimensional "surface" area of the landslide will also affect the volume to area relationships. It is outside the scope of this paper to investigate such issues, which may be examined when other three-dimensional high-resolution topographic data sets become available. In addition, Larsen et al. (2010) use the reduced major axis regression technique to fit power-law trends to the data and not the general linear models used in this study, which could potentially contribute to the differences in the fitted parameters. Such impacts have not been explored in this study.
The effects of vegetation on the total volume of debris estimated from the 17,256 sources do not appear to have as much influence as shadows in the imagery used to generate the DSMs. The total landslide volumes estimated from the power-law fit parameters derived from the landslides within each NDVI band are consistent for those bands representing bare ground to high vegetation cover.

Total Landslide Volumes and Their Uncertainty
To investigate the impact of the reported scaling relationships-and the relatively minor changes between them-on the total estimate of landslide debris volume, we have calculated the total volume of debris produced by the 17,256 debris avalanches from the InSAR-adjusted and -corrected difference model. We have then compared this volume with those calculated from the given power-law scaling relationships derived from this study and those in the literature, using the individual source areas of each debris avalanche ( Figure 10b). The results in Figure 10b show that the total volume estimated in this study for mixed debris avalanches (in soil and rock combined) adopting the weighted power-law model fit, of 38.2 M m 3 is comparable to the total volume estimated from the cleaned data set (33.3 M m 3 ) and the global (mixed) model reported by Larsen et al. (2010), of 35.6 M m 3 , but these estimates are below the volumes from the difference model, which is 45.8 M m 3 . The total volume from the weighted power-law model fitted to the avalanches in rock from this study of 40.2 M m 3 provides the closest comparison to the total volume from the difference model. This is possibly because the landslide distribution is dominated by avalanches in GeolCode 4 greywacke rock. Most of the power-law fits, however, provide a slightly lower total volumes than the total mean debris avalanche volume estimated from the difference model, with the exception of the global bedrock (Larsen et al., 2010) (53.2 M m 3 ) and global (Guzzetti et al., 2009) (51.6 M m 3 ) power-law models, which overestimate the total volume.
The relatively large uncertainty (given as 1σ) in the total volume estimated from DSM differencing represents the true uncertainty of the volume estimated for each avalanche source area. The uncertainty range of the total debris avalanche volume estimated from the difference model was calculated by adding and subtracting the bootstrap standard deviation of ±1.8 m (1σ) of the difference model, to the 2015 DSM, as systematic offsets. These "offset" DSMs were then subtracted from the 2017 DSM to create +1σ and −1σ difference models, which were subsequently used to estimate the +1σ and −1σ source volumes for each debris avalanche. The fitted "model" uncertainty estimates potentially may underestimate the actual uncertainty range, when compared to the uncertainty associated with the volumes estimated for each source area from a difference model. The uncertainties estimated using the power law models-from the literature and this study-are much smaller but represent the uncertainty of the given model fitted to each of the data sets used. The cleaned avalanche source volumes derived from the difference model used in this study-for which we have high confidence on the volume measurements-show that for any given source area, the volume could vary over two orders of magnitude. The volume estimated from power-law scaling relationships for each landslide should ideally be represented as a range with models fitted to the, e.g., "upper" and "lower" volume estimates. If such measurement uncertainties were considered in the model fitting, it is possible that they would be similar in magnitude, or larger, than those estimated from DSM differencing used in this study.
To calculate the total volume of landslides material mobilized by the earthquake, it is important to estimate the volumes of the other 12,213 landslides, where the difference model cannot be used. Of these, about 80% (approx. 9,000) are debris avalanches, while the other 20% (approx. 3,000) are mainly rotational, translational or compound slides in rock. For comparison, we have plotted the total volume (and the estimated uncertainty) of the 50 rotational, translational, or compound slides in rock (Figure 10b), based on manual estimates of their volume. These 50 landslides represent most of the larger landslides triggered by the earthquake, and therefore, they represent a significant proportion of the total volume of the landslide material mobilized by the earthquake. The total volume of the remaining debris avalanches could be estimated from their source areas using for example, the upper and lower volume estimates, for the given source area, from creating envelopes around the cleaned data (Figure 9), and not the power-law trends fitted to the data, thus reflecting the true uncertainty of the estimated volume. The depth and therefore volume of the remaining rotational, translational, or compound slides in rock-albeit smaller in area than the 50 for which the volumes have been calculated-will be controlled by the structural geology and therefore will vary from one landslide to the next, making it difficult to estimate their volumes remotely. Estimating the volumes of such landslides will form an interesting future research topic, as such landslides could contribute significantly to the total volume of landslide material mobilized by earthquakes. However, the debris from such landslides triggered by the Kaikōura Earthquake has remained relatively intact and within the source area, with only minor remobilization post-earthquake, thus providing only small volumes of sediment to the river system.

Conclusions
About 80% of the landslides within Version 2.0 of the Kaikōura Earthquake landslide inventory can be classified as soil or rock avalanches, with the remaining landslides being mainly translational (or rotational or compound) slides. We found that the source geology and landslide type/failure mode influenced the source volumes of the landslides triggered by this earthquake. Of the 17,256 avalanches, where accurate estimates of their source volumes could be calculated, 79% were in Pahau terrane greywacke sandstones with γ values of 1.040-1.138. The larger a landslide is in area, the more likely it is that the slide plane will be deeper and within the underlying, highly jointed less weathered rock, making the differentiation between soil and rock as the dominant source material difficult. The remaining 3,671 avalanches were predominantly in colluvium, residual soils, and completely-to highly-weathered rocks derived from upper Cretaceous and Neogene sandstones, siltstones, and limestones. These yielded γ values of 0.822-1.060, suggesting that the depth of the landslide source remains relatively shallow irrespective of how big its area becomes. In contrast, for 50 selected rotational, translational, or compound slides in the same materials, where the debris remains in much of the source area, thus obscuring the slide surface, we found γ values range between 1.46 and 1.47, indicating that their slide surfaces are considerably deeper than those landslides classified as avalanches.
Regardless of the source material, all the γ values for the soil and rock avalanches used in this study were at the bottom end of the range reported in the literature for both soil and rock landslides and more similar to those γ values estimated from other coseismic landslide inventories. This suggests that the Inangahua, Murchison, and Wenchuan earthquake landslide inventories, like the Kaikōura landslide inventory, contain a high proportion of shallow avalanches. This study corroborates the results from previous studies (e.g., Keefer, 2002) on coseismic landslides, which determined that debris (soil and rock) avalanches/disrupted landslides are the dominant landslide type triggered by earthquakes. This study also corroborates the results from previous studies (e.g., Larsen et al., 2010) that differentiation between material type is critical for accurate calculation of landslide volumes from source-area geometries, but we also show that landslide type/failure mode also influences the landslide source volume and should also be considered when estimating landslide volumes from statistical landslide area to volume scaling relationships. Future research is needed to further explore the results from this study by investigating the role of slope angle, geology, and landslide failure mode on the resulting landslide volume and how the results compare to other landslide data sets derived from DSM differencing.