Modeling the area of co-seismic landslides via data-driven models: The Kaik ¯ oura example

substantial developments in data-driven models for landslide prediction. However, this improvement has been mostly devoted to models that estimate locations where landslides may occur, i


Introduction
The estimation of where landslides may occur in the future has dominated the geomorphological literature pertaining to data-driven applications since its first conceptualization in the early 1970s (Reichenbach et al., 2018).Almost no other data-driven modeling framework with a spatially-explicit (i.e., an explicit representation of a geographic space) connotation has been developed for the subsequent five decades other than the susceptibility (Fell et al., 2008).This concept boils down to estimating the probability of a landslide occurring in a given mapping unit under the influence of topography and other thematic landscape characteristics (Van Westen et al., 2008).An extension to this framework is present in the literature, although less prominent than the susceptibility, and it features the spatiotemporal characteristics of the landslide trigger, leading to the estimation of the hazard concept (Guzzetti et al., 1999).This extension has led to the development of important forecasting tools such as near-real-time models (Nowicki Jessee et al., 2018) and early warning systems (Intrieri et al., 2012;Kirschbaum and Stanley, 2018).In both cases, though, the model behind the respective results still targets landslide occurrence data in the form of presences or absences across a given landscape (Frattini et al., 2010).
As technology advanced, the information on unstable slopes also changed, being acquired and processed in multiple ways.For instance, at the origin of the susceptibility concept, geomorphologists were observing the landscape and labeling slopes to be likely stable or unstable on the basis of personal experience (Brabb et al., 1972).The birth of GIS facilitated the development of numerical tools, which started from simple analytical approaches such as bivariate statistics (e.g., Van Westen et al., 1997) and evolved into more complex modeling schemes where multiple variables simultaneously contribute to estimate the susceptibility (e.g., Steger et al., 2016).The latter frameworks belong to different and complementary approaches that have taken root in the landslide community.One corresponds to statistical models where model interpretability is favored at the expense of reduced flexibility and lower performance (i.e., analytical tasks).The other one corresponds to machine learning, where performance maximization is sought at the expense of interpretability (i.e., prediction task).Two common examples correspond to statistical models such as Generalized Linear Models (e.g., Castro Camilo et al., 2017) and machine learning models such as decision trees (e.g., Yeon et al., 2010) or neural networks (e.g., Wang et al., 2021).Between these two lies the Generalized Additive Model, also referred to as an interpretable machine learning technique (Steger et al., 2021).They ensure a high degree of interpretability typical of statistical architectures, but their structure allows for incorporating nonlinear effects, leading to flexible models with high performance (Lin et al., 2021).Nevertheless, because their target variable is a binary realization of landslide occurrences, these models lack the ability to return information on how large a landslide may be (Lombardo et al., 2021) or on how many coalescing landslides may initiate in a particular region (Lombardo et al., 2019).
There are three previous articles where a data-driven model is used to estimate landslide areas at regional scale.The first is Lombardo et al. (2021), which estimates the maximum and the sum of planimetric landslide areas within slope units.That model has a global connotation: it implies that knowing its validity in a worldwide context is not sufficient because this scale is not applicable to local territorial management.Alternatively, physically-based models are applied to predict landslide size (Alvioli et al., 2014), although the limited availability of geotechnical and hydrological data restricts their application to small regions.The applicability of the model proposed in Lombardo et al. (2021) still needs to be validated for site-specific conditions.Moreover, it must undergo tests in case of seismic-and rainfall-induced landslides.To this purpose, other contributions try to replicate a similar experimental setting at regional scales, focusing on co-seismic landslides (Aguilera et al., 2022) and rainfall-triggered landslides (Bryce et al., 2022).
In this work, we seek to produce a new analysis that follows the main workflow direction of the articles mentioned above but introduces untested innovations, such as topographic correction of the landslide area target variable and the validation of model performances using a spatially-explicit model validation.We selected the Kaikōura earthquake (7.8 M w , 13-11-2016), for which a suitable co-seismic landslide inventory is available (Tanyas et al., 2022).We partitioned the area affected by landslides into slope units (Alvioli et al., 2016), extracted the sum of all landslide extents falling within each mapping unit, and calculated the topographically-corrected surface areas.As for the model, we adopted a Generalized Additive Model structure under the assumption that the aggregated landslide area per slope unit followed a Log-Gaussian likelihood (cf.Lombardo et al., 2021).

Study area and co-seismic landslides
The Kaikōura earthquake struck the South Island of New Zealand on UTC 14 th of November 2016 at 11:02.This was the largest magnitude earthquake in New Zealand since 1855 and had a very complex rupturing mechanism (Hamling et al., 2017).The earthquake cascaded across a series of faults with dextral, sinistral, oblique, and reverse rupture mechanisms (Diederichs et al., 2019).Significant co-seismic surface deformations occurred across a large landscape extending up to 100 km from the epicenter (Cesca et al., 2017).The reported uplift amount reached up to 8 meters in some locations (Hamling et al., 2017).As a result, the earthquake severely damaged infrastructure and altered the environment (Kaiser et al., 2017).
Considering the steep terrain affected by considerable ground shaking, the earthquake resulted in a large number of landslides (Tanyas et al., 2022).Massey et al. (2020) reported more than 29,000 landslides triggered by the earthquake, whereas Tanyas et al. (2022), using slightly different criteria and lower resolution imagery, mapped 14,233 landslides over an area of about 14,000 km 2 .Considering the documented earthquake-induced landslide events (Tanyas ¸ et al., 2017), the Kaikōura event is one of the largest recorded in the literature.
This study examines an area affected by the 2016 Kaikōura earthquake (see Fig. 1) using the landslide inventory of Tanyas et al. (2022), which delineates landslide sources and deposit areas as single polygons.The inventory consists of various landslide types, including disrupted rock, debris, soil falls and slides.However, landslide types are not indicated in the original data; thus, our analyses are not sensitive to any specific type of landslide.

Slope Units
The use of a Slope Unit (SU) delineation in the framework of landslide predictive models dates back to Carrara (1988).The spatial extent of this mapping unit is usually coarser than the more common grid cells.The latter are regular polygonal objects that offer a simple spatial partition of any landscape, mainly by matching the gridded resolution of the Digital Elevation Model (DEM) available for the given study area.They are suited to express the spatial variability of continuous phenomena, such as temperature fields.However, landslides are discrete processes.The geoscientific community has long debated whether grid cells are suitable for modeling slope failures.Conversely, SUs are more suitable from a geomorphological perspective, although they require additional processing steps, such as the aggregation of fine-scaled landscape characteristics.Several automated tools have been proposed and shared within the geoscientific community (Alvioli et al., 2020).In this work, we selected SUs to partition the area affected by the Kaikōura earthquake to predict the cumulated area of landslides per mapping unit.Below we report the parameterization of r.slopeunits, the software we used.As for their interpretation, we refer to Alvioli et al. (2016).
• Circular variance = 0.4 • Flow accumulation threshold = 1,000,000 • Minimum Slope Unit area = 80,000 The circular variance is a parameter that controls how homogeneous the aspect criterion should be.For instance, a circular variance of 0 would result in a strict selection where a few pixels could be merged into a SU.Conversely, a circular variance of 1 would result in a flexible search where pixels with large differences in slope exposition would be merged in a single SU.As for the flow accumulation, this parameter controls the starting planimetric area for the SU partition.Any subsequent r.slopeunits iteration would reduce the SU unit extent.The minimum SU area and the cleansize control the final outcome.The first defines the average target planimetric area for a SU to be considered, and the second represents the size below which any artifact SU should be merged to the adjacent one.
The resulting SUs offered a medium resolution of the exposed Kaikōura landscape with 26,839 total SUs, whose planimetric area distribution has a mean of 500,000 m 2 and a standard deviation of 430,000 m 2 .

Covariates: landscape characteristics and ground motion data
This section illustrates the covariates we used to explain the variability of the co-seismic landslide area distribution in the affected Kaikōura region.Although there is extensive literature examining factors governing the probability of landslide occurrence, factors controlling the area of landslides in a spatial context is a relatively new concept (e.g., Lombardo et al., 2021).In this regard, we tested several variables representing morphometric, anthropogenic, and seismic factors as well as material properties (see Table 1 and Fig. 2).We used a 25-m-resolution DEM provided by The Land Resource Information Systems Portal of New Zealand (LRIS; https://lris.scinfo.org.nz) and tested some basic DEM derivatives, namely, slope steepness (Slope), northness (NN), eastness (EN), local relief (Relief), profile curvature (PRC) and planar curvature (PLC) to assess the role of morphometric variables on landslide area.Capturing the role of anthropogenic factors is challenging (e. g., Tanyas ¸ et al., 2022), but the study area is a remote territory, and the road cuts are the main features representing human influence on landsliding.Therefore, we calculated the Euclidean distance to the road network (e.g., Lepore et al., 2012) to capture a possible influence of anthropogenic factors.Specifically, we accessed the road network map of the study area via the Land Information Portal (https://data.linz.govt.nz) of New Zealand.As for the co-seismic ground shaking, we used the Peak Ground Acceleration (PGA) map of the Kaikōura earthquake provided by the U.S. Geological Survey (USGS) ShakeMap system (Worden and Wald, 2016).PGA is a seismic proxy, and specifically, the deterministic estimate of PGA provided by the USGS ShakeMap system is widely used in susceptibility analyses of co-seismic landslides (e.g., Nowicki et al., 2014).Also, we used a pedological soil thickness map of the study area (Lilburne et al., 2012) via the LRIS portal as a proxy for the shear strength of hillslope materials.Different from all the other covariates, we examined the soil thickness map as a categorical covariate because it includes four categories where soil depth is described as deep (D, >90 cm), moderately deep (MD, 40-90 cm), shallow (S, 20-40 cm) and very shallow (VS, <20 cm) as well as a category indicating no soil (NS) cover (note that a more desirable regolith thickness map does not exist).

Table 1
Covariate summary table.Each covariate listed here was later used in a dual form during the analyses.Specifically, we represented each covariate in this table through the mean and standard deviation values computed per SU.We do not list both terms in the table, but they will be denoted in the remainder of the manuscript via the suffix _mean and _stdev added to the acronyms reported in the table.

Data aggregation at the Slope Unit level
We used slope units to aggregate both the target variable, this being the topographically-corrected landslide area, and the covariates described in the previous Section.
The landslide extent calculation was based on the aggregation of the landslide area by summing up all landslide areas within each SU.Before this aggregation step, though, we applied a correction to reduce the underestimation of landslide area on steeper terrain due to the underlying planar projection.For this purpose, a trigonometric function based on a slope-steepness map was used to derive the "true" surface area of each landslide polygon in analogy to Steger et al. (2021).
Fig. 3a shows the distribution of the topographically-adjusted landslide area after the aggregation step mentioned above (sum for each SUs).Being the distribution strongly heavy-tailed, we opted to take the logarithm of the cumulative landslide area per SU Fig. 3b.In such a way, a Log-Gaussian model could be used to suitably explain the variability of these estimates (more details in Section 3.1).

Fig. 2.
Example of the covariate set used for the analyses.The soil depth map includes five classes: NS for No Soil, VS for Very Shallow, S for Shallow, MD for Moderately Deep, and D for Deep.Notably, the Dist2R map is shown as is only for graphical purposes.We constrained the information conveyed by Dist2R (to the model we will describe in Section 3) only up to a 500 m buffer around the road network.After this distance, we impose the covariate to cease to be informative.
As for the landscape characteristics, we computed the mean and standard deviation of each continuous covariate per SU (see Titti et al., 2022).Ultimately, we only extracted the dominant type per SU whenever the landscape characteristics corresponded to categorical properties, such as underlying lithology, land use, or soil thickness classes.

Modeling strategy
Below we provide a brief description of the model we adopted, the cross-validation scheme we implemented, and the metrics we used to assess how the estimated landslide areas matched the observed cases.

Generalized Additive Model
Generalized Linear Models (GLM) are statistical techniques designed to model linear relationships between a target variable and a set of predictors.A Generalized Additive Model (GAM) is a more flexible extension of a GLM.In analogy to GLMs, GAMs can handle a variety of error distributions but additionally account for nonlinear associations between the target variable and continuous predictors.This flexibility combined with high interpretability makes GAMs useful in data-driven studies.The presence of nonlinear relationships between landslide occurrence and environmental factors can be expected (e.g., landslides may less likely occur in flat and very steep terrain), while high interpretability of the modeling results is paramount for geomorphological interpretation and plausibility checking (Steger et al., 2017).GAMs with a binomial error distribution have been applied to model landslide susceptibility (Petschko et al., 2014), while Poissonian GAMs were used to model spatial landslide counts (i.e., intensities; Lombardo et al., 2020).A Log-Gaussian distribution within a Bayesian GAM forms the foundation to create the first data-driven model predicting landslide size per SUs, i.e., the maximum landslide area and the sum of landslide area (Lombardo et al., 2021).The Log-Gaussian GAM used within this study is based on the R-package "mgcv" (Wood and Augustin, 2002).This framework allowed us to model the topographically corrected log-area of co-seismic landslide areas at the SU-level (hereafter L A ) as a function of a covariate set that describes landscape characteristics and spatial ground motion properties.The nonlinear relationships were fitted while we restricted the maximum allowed flexibility of the underlying smoothing functions to a k-value of 4 (i.e., the maximum allowed degrees of freedom) to enhance model generalization and interpretability.The covariates are described in detail within Section 2.3, while their selection was based on a systematic procedure that included an iterative fitting and evaluation of different model realizations.In detail, we started with a full model and iteratively excluded covariates that did not meet the following two criteria: a covariate was only considered appropriate in case the underlying smoothing term was estimated to be significant at the five percent level (p-value < 0.05); a covariate did not enhance the model's predictive performance.
Besides handling nonlinear relationships, GAMs also allow visualizing the respective modeled associations between the target variable and the predictors.This model transparency is useful to enable interpretation and uncover implausible results (Zuur et al., 2009).In this sense, component smooth function (CSF; i.e., partial effects) plots were used to visualize the estimated covariate-response relationship.These plots enabled the interpretation of nonlinear effects on the aggregated landslide area per SU at a single covariate level while simultaneously accounting for the influence of other covariates in the model (Molnar, 2019).

Model performance
Below we summarise the cross-validation schemes adopted and the metrics used to assess how the model performed in explaining the spatial distribution of landslide areas.The last section explains how we then provided estimates of landslide areas for SUs that did not experience slope failures during the Kaikōura earthquake.

Cross-validation routines
To test the performance of our model, we selected two crossvalidation approaches.The first corresponded to a random crossvalidation scheme (RCV), where we repeatedly extracted a random subset of 90% of SUs within the study area for training our model (i.e., a training set) while the remaining data (i.e., a test set) of each repetition was used to calculate the performance metric.The random selection was constrained to select the same SU only once.Thus, the union of the ten replicates returns all the SUs constituting the whole study area.
However, a spatial process may embed some degree of internal spatial dependence, which may not be fully explained by a chosen covariate set.Conventional non-spatial random partitioning of training and test sets (e.g., RCV) may provide test statistics that do not capture the variability of model performance across sub-regions of a study site.Using RCV, over-optimistic performance scores are likely to be measured if spatial model predictions poorly match data within single sub-regions of an area.Spatially-explicit validation schemes, such as spatial cross-validation (SCV), can be used to estimate the spatial transferability of model performance scores within a study site and uncover spatially incoherent model predictions (Steger et al., 2017).SCV results can inform potential users of a given model about worst-case prediction skills in space about the spatial robustness of the general model setup.SCV is usually based on a repeated random splitting of training sets and test sets according to sub-areas of a study site.For this study, the underlying spatial partitioning approach was based on kmeans clustering (see Brenning, 2012 for a more detailed explanation).
In this work, we opted to report the model performance estimated via an RCV where the prediction skill was aided by residual clustering effects, as well as via an SCV where the estimated performance scores are usually lower, thus providing insight into the minimum prediction skill expected for sub-regions of the study site.Fig. 4 shows a few examples of the routines mentioned above.Specifically, the 10-fold RCV and the 10fold SCV have been iterated ten times to randomize the spatial cluster of slope units to be extracted.

Performance metrics
To assess how suitable our modeling framework is to reflect the observed landslide area per SU, we selected a dual approach featuring visual and numerical performance summaries for both CV schemes described above.The visual summary corresponds to a simple graph plotting the observed landslide areas against the estimated ones.As for the numerical summaries, the metrics we used were the Pearson Correlation Coefficient (R-Pearson; Schober et al., 2018) and Mean Absolute Error (MAE; Mayer and Butler, 1993).To these, we also added the Root Mean Square Error (RMSE; Kenney and Keeping, 1962) for completeness, although several contributions mention that the MAE is a better measure of deviance (Willmott and Matsuura, 2005).

Map-based landslide area prediction
In this section, we specify something of particular conceptual relevance.In traditional susceptibility models, one can and should use the presence-absence information across the whole study area (Petschko et al., 2014).However, information on landslide area is only associated with a subset of the SUs partitioning the Kaikōura landscape.Therefore, to produce maps of the predicted landslide area for the whole study area, we adopted the following procedure.Initially, we extracted the positive landslide areas to train and test our Log-Gaussian GAM.Subsequently, we implemented a simulation step using the estimated regression coefficients to solve the predictive equation in areas where landslide area information was unavailable.

Results
Below we present an interpretation of the model components, performance, and mapping results.

Model relationships
This section summarizes the estimated covariate effects responsible for explaining the spatial distribution of landslide areas per SU.
Fig. 5 offers an overview of all the nonlinear effects we included in the model.Although we allowed the regression coefficient to vary nonlinearly across each covariate domain, the implemented internal smoothness selection procedure selected certain covariates to be best represented via linear functions.This was the case for Slope_stdev, NN_mean, PRC_stdev, and PGA_mean.This implies that a unit increase in the covariate value would generate a proportional changedepending on the sign of the regression coefficientin the resulting landslide area.And, that the change would be the same irrespective of where that unit increase happened across the whole covariate spectrum.Moreover,  eight covariates deviated from linear behavior, of which two were only mildly nonlinear NN_stdev, Dist2R_stdev, whereas the remaining six showed much more nonlinearity (Slope_mean, EN_mean, Relief, PLC_mean, Dist2R_mean, and PGA_stdev).
Below we provide a brief overview of these covariate effects (from the most interesting linear to the nonlinear ones) by interpreting their marginal contribution (i.e., assuming all the other covariate contributions to be fixed).For example, we justify the positive increase of the estimated landslide area due to Slope_stdev because a rougher terrain may have larger quantities of hanging material susceptible to be mobilized due to the contextual water impoundment (Jiao et al., 2014).Similarly, the PGA_mean positively contributes to the estimated landslide area, and its linear behavior may be seen as a destabilizing effect of ground motion on the landscape (Tanyas and Lombardo, 2019).Furthermore, two covariates share similar nonlinear contributions.These are Relief and PLC_mean, both with pronounced sigmoidal behavior.The former can be interpreted with the positive contribution of gravitational potential energy, where at increasing values, the failing mass will experience increasing kinetic energy as it moves downhill, thus producing larger landslides overall (Yamada et al., 2018).As for PLC_mean, the planar curvature is known to control the convergence of granular materials and overland flows over a landscape (Ohlmacher, 2007).
Aside from covariates we allowed to behave nonlinearly while still carrying their ordinal structure, we also considered the nonlinear and categorical signal of soil thickness classes.As it stands out in Fig. 6, the signal carried by the prevalent soil depth class per SU does not produce a clear "monotonic" pattern in the estimated regression coefficients per class (i.e., landslide area increases or decreases systematically with soil depth).This is likely due to two reasons.First, the raw soil depth map we accessed is directly expressed into classes, which implies a loss in the continuous information a soil depth should be expressed into.Clearly, soil depth cannot be continuously measured over space because it would require excessive resources.Therefore, even the classes we used are the result of an interpolation routine, which may have smoothed the soil depth signal over space.A second and valid reason for the not straightforward interpretation emerging in Fig. 6 is that we aggregated the soil depth signal over the SU by choosing the majority rule.In this sense, a given SU is assigned with the soil depth label of the class with the largest areal extent.However, the majority class may not be the one responsible for the failure.

Model performance
The visual agreement between the observed and estimated landslide area among the three model routines we tested is summarized in Fig. 7.
There, one can see that the model fit produces the highest degree of agreement between the observed and estimated landslide areas.In the second panel, the RCV-predicted landslide areas closely follow the trend shown for the fit.As for the SCV results, the deviations from a perfect match between observed and estimated landslide areas appears slightly more pronounced compared with the other two cases.However, this is to be expected because a SCV essentially removes any residual dependence from a spatially distributed dataset, thus producing lower performance scores in a real-world data setting.In this sense, the match shown for the SCV can still be considered a suitable source of information for hazard assessment.
Fig. 8 complements Fig. 7 by informing on the correlation between observed and estimated landslide areas, together with the error between the two.Several authors have proposed a classification of the R-Pearson, and most of the literature on the topic would indicate values of around 0.6 to reflect a moderate (Mukaka, 2012) to strong (Corder and Foreman, 2011) correlation between observed and estimated landslide extents.Analogous considerations arise by examining the MAE and RMSE, with acceptable errors in both cross-validation schemes.Notably, the performance metrics reported in Fig. 8 confirm that the SCV returned a slightly poorer agreement compared to a purely random cross-validation scheme.

Landslide area predictive maps
Fitting a statistical model allows one to retrieve the set of regression coefficients through which one can estimate the expected values of the given target variable.At the same time, one can use the same set of regression coefficients to solve the predictive function for locations where the target variable is unknown.The latter concept boils down to what one could refer to as a statistical simulation (e.g., Luo et al., 2021) or model transferability (e.g., Petschko et al., 2014).Fig. 9 summarizes the estimates of the two cross-validation schemes at SUs for which we have L A observations, as well as SUs where we have not.The first row highlights the agreement in spatial patterns among the observed and predicted L A values, with a coherent pattern shown among the three images, albeit the prediction routines show some degree of smoothing as they transition from RCV to SCV.The strength of our modeling framework is shown in the second row of Fig. 9, where we transfer the predictive equations to the remainder of the Kaikōura landscape.

Discussion
The capacity of data-driven models to go beyond traditional susceptibility models is still in infancy stage.This analysis suggests that a Log-Gaussian GAM can reproduce the pattern and value range of landslide areas aggregated at the slope unit level.Out of the whole procedure, certain elements already support the replication of similar analyses, while others call for further improvements.Two of these elements are discussed separately in the sections below.

Supporting arguments
Landslide area correction with respect to slope steepness is not usually considered in most geoscientific contributions, with the exception of very few cases (e.g., Steger et al., 2021).In the context of a model that aims at predicting landslide area, we consider this an additional element to be added to the protocol proposed in Lombardo et al. (2021).A further addition is the use of spatial cross-validation.Lombardo et al. (2021) constrain the spatial cross-validation to be generated once.Herein we focus on a specific site, which makes it easier for us to replicate the spatial sampling, thus fully randomizing the spatial crossvalidation results, in line with what Brenning (2012) prescribes, albeit in a binary context.
The performances we retrieved suggest that it is worth extending the landslide area prediction further.Figs.7 and 8 show the extent to which our model estimates the observed landslide areas.Also, this is translated over the geographic space in Fig. 9, where the spatial patterns match, albeit the predictive routines show deviation from the original A L values as the tested cross-validation routines moved from the random context to the spatially constrained one.Further improvements may consolidate the concept and role of landslide area prediction within protocols of disaster risk reduction.For instance, we could combine the area model with a traditional susceptibility one.As things are, the traditional susceptibility framework does not formally account for the expected area of landslides once triggered in a given slope.However, even our landslide area framework is blind to whether a slope may be prone or not to fail.In turn, this means that these two tools are currently separated, and further efforts could be directed toward merging them into a single product that integrates two hazard features, namely spatial landslide probability and landslide area.For instance, one could model them separately and then combine them.In such a way, slopes that may morphologically be associated with large failures but are seen to be stable (low probability of occurrence) by the susceptibility component will result in a small hazard-proxy value.The same may happen in the case of slopes that are expected to be unstable (high probability of occurrence) but associated with very small landslides.In such a scenario, the estimated hazard proxy will also be low.On the contrary, only in situations where high susceptibility is associated with large expected landslides, one would obtain a level of such a hazard proxy that would inevitably require attention.Such a scheme may give rise to a new landslide hazard framework, providing a full spectrum of probabilistic estimates aimed at aiding decision-making processes for landslide risk reduction.An alternative we envision could include physics-based simulations typical of slope-scale engineering solutions.Assuming we could simulate landslides for single slopes, the physics-based framework would ensure retrieving intensity measures such as landslide velocities, kinetic energy, or impact force.These could then be passed to a similar model to the one presented here to predict momentum-related measures rather than landslide areas.

Opposing arguments
To critically review our landslide area model, one should take a step back and look at the model's fundamentals.The fact that it relies on a logarithmic transformation of the landslide area distribution per SU requires some consideration.From a mathematical perspective, this framework is sufficient to produce valuable predictive maps, as the logarithm is a monotonic transformation.Thus, landslide areas that are smaller compared to other SUs in the observed data will still be relatively smaller in the prediction, irrespective of whether we direclty model the landslide extent in linear or log scale.However, two negative elements affect this framework.An obvious one is that from an interpretation standpoint, one lacks the intuition of what a predicted value would indicate at the log scale.This argument could still be acceptable because of the monotonic transformation mentioned above.But,  reflecting on what this entails in terms of errors calls for potential improvements.A Gaussian likelihood implies that the model focuses on the bulk of the landslide area distribution.In other words, the mean landslide area will be suitably estimated, leaving the tails potentially misrepresented.The left tail, the side of the distribution with very small landslides, may be of lesser interest.However, a misrepresentation of the right tail, the side of the distribution that hosts very large landslides, can lead to poor decisions for the more dangerous ones.Notably, the performance we produced does not raise concerns to the point of considering our landslide area model inappropriate.However, we envision the next phase of the model development to explore more suitable likelihoods.The log-Gaussian context is appealing because of its easy implementation, and as long as the performance may stay along the lines of what we presented here, the choice of such likelihood can be justified.However, in the hope of further extending the landslide area prediction in different geographic contexts across different landslide types and triggers, we cannot exclude that the likelihood we chose so far may prove to be insufficient or lead to undesired errors away from the bulk of the distribution.In such cases, extreme-value theory in statistics may provide a better modeling framework.We already envision this direction to be the following research and development phase.Some of this development can be seen already in Yadav et al. (2022), where a number of distributions are tested in a Marked Point Process framework.

Conclusions
Data-driven landslide models have relied essentially on the same toolbox for over five decades.We believe it is time to review whether new tools could improve the static susceptibility framework and complement the information it provides.One of these elements consists of how landslide areas may enlarge after the landslide initiates, evolves, and potentially coalesces into larger areas of materials moving downhill.This information has been traditionally associated with physically-based models, and other kinematic parameters such as velocity.On the one hand, landslide kinematics cannot be modeled in detail via data-driven approaches because of the lack of observations.On the other hand, though, the landslide area information is contained in any standard landslide polygonal inventory.As a result, data-driven models can be trained to learn what environmental characteristics promote small to large landslides and translate this information into maps of the expected landslide areas.This idea is an uncharted territory within the geoscientific community, with few articles currently addressing this issue.
We consider our landslide area model a new venue of potential scientific interest, and prompt the geoscientific community to explore this framework further.To ensure dissemination, we shared data and codes in a GitHub repository, accessible at https://github.com/mmorenoz/GAM_LandslideSize.Further improvements may be directed toward fitting different statistical distributions tailored toward extremely large landslides or to extend the current spatial model toward its space-time counterpart.Potential implications may be translated into better hazard information for administrations to base their risk reduction plans.

Fig. 1 .
Fig. 1.Geographic summary of the co-seismic landslides triggered in response to the Kaikōura earthquake (panels a and b).Panel c shows an example of the slope unit delineation superimposed on an aspect map.

Fig. 3 .
Fig. 3. Distribution of the topographically-corrected landslide areas per SU.Panel (a) shows the sum of derived landslide areas per SU in a linear scale, whereas panel (b) highlights the same information in a logarithmic scale.

Fig. 4 .
Fig. 4. Geographical sketches of CV routines via five examples of the ten folds we implemented in this work.The first row shows an RCV, whereas the second row highlights the effect of a spatial constraint in the SU selection.

Fig. 5 .
Fig. 5. Summary of ordinal nonlinear effects on the aggregated landslide area per SU.

Fig. 6 .
Fig. 6.Summary of the categorical nonlinear effect of soil depth classes on the aggregated landslide area per SU.

Fig. 7 .
Fig. 7. Summary of the agreement between the observed landslide area per SU and the corresponding values estimated through a fit and two cross-validation (RCV and SCV) routines.

Fig. 8 .
Fig. 8. Pearson correlation coefficient, mean absolute error and root mean square error estimated for the purely random cross-validation and the spatial random cross-validation, respectively.

Fig. 9 .
Fig. 9.The first row of this figure highlights the details of the main area affected by landslides for which we have observations.The second row shows the whole study area without focusing on the SU for which we measured the landslide area.The first column plots the actual measurements and represents the target variable of our model.The second and third columns report the estimated landslide areas via the RCV and SCV routines.