Unified landslide hazard assessment using hurdle models: a case study in the Island of Dominica

Climatically-induced natural hazards are a threat to communities. They can cause life losses and heavy damage to infrastructure, and due to climate change, they have become increasingly frequent. This is especially the case in tropical regions, where major hurricanes have consistently appeared in recent history. Such events induce damage due to the high wind speed they carry, and the high intensity/duration of rainfall they discharge can further induce a chain of hydro-morphological hazards in the form of widespread debris slides/flows. The way the scientific community has developed preparatory steps to mitigate the potential damage of these hydro-morphological threats includes assessing where they are likely to manifest across a given landscape. This concept is referred to as susceptibility, and it is commonly achieved by implementing binary classifiers to estimate probabilities of landslide occurrences. However, predicting where landslides can occur may not be sufficient information, for it fails to convey how large landslides may be. This work proposes using a flexible Bernoulli-log-Gaussian hurdle model to simultaneously model landslide occurrence and size per areal unit. Covariate and spatial information are introduced using a generalised additive modelling framework. To cope with the high spatial resolution of the data, our model uses a Markovian representation of the Matérn covariance function based on the stochastic partial differential equation approach. Assuming Gaussian priors, our model can be integrated into the class of latent Gaussian models, for which inference is conveniently performed based on the integrated nested Laplace approximation method. We use our modelling approach in Dominica, where hurricane Maria (September 2017) induced thousands of shallow flow-like landslides passing over the island. Our results show that we can not only estimate where landslides may occur and how large they may be, but we can also combine this information in a unified landslide hazard model, which is the first of its kind.


Introduction
Dominica is situated within the Atlantic hurricane belt and, as such, is highly vulnerable to high-intensity weather events, as evidenced by its long extreme-events history, stretching as far back as the Great hurricane of 1780. In addition, the island's infrastructure and populations are concentrated along the coastal areas, particularly in the south and west, where their situation and building conditions hardly help withstand one extreme natural hazard after another. In this paper, we focus on the aftermath of hurricane Maria, which originated from a wave leaving the west coast of Africa on 12th September 2017. The wave moved westward across the Atlantic, creating deep convection, and consolidated into a tropical depression 580 nautical miles east of Barbados on the 16th September (Pasch et al. 2018;Fobert et al. 2021). By the evening of the 18th, Maria had intensified and hit the Island of Dominica twelve hours later as a category five hurricane (Pasch et al. 2018), causing a total of 9960 landslides on the island, 8,576 of which were classified as debris slides (Van Westen and Zhang 2018). The infrastructure and transport sectors were the main avenues of the USD 930.9 million in damages incurred by Maria, with the agriculture and tourism sectors not far behind (ACAPS 2018). This prompted a plea by Dominica's prime minister for the funds to make Dominica a fully climate-resilient nation, determined to protect their island (Gibbens 2019).
The Government of the Commonwealth of Dominica, in conjunction with the Caribbean Disaster Emergency Management Agency and the Caribbean Development Bank, commissioned a post-disaster needs assessment from hurricane Maria in order to estimate the total damage, the damages per sector, and to identify recovery needs (A-CAPS 2018). The reconstruction of the destroyed infrastructure was funded by a loan from the World Bank and was a part of the project ''Enhancing Resilient Reconstruction in Dominica''. The project promoted the idea that Dominica could limit the damage from natural hazards by improving the uptake of resilient building practices, aiming to accelerate short-term recovery and strengthen long-term resilience to climate-related risks. Landslide hazard assessment is critically valuable to this programme. It can help define land-use capability, detect areas where intervention is needed to stabilise slopes, and identify appropriate mitigation measures.
To help tackle the above challenges, this article builds a joint probability model to map the Island of Dominica by susceptibility and size of landslides. Our model is able to detect unstable areas that could potentially host large landslides, which is key to identifying regions in need of hazard mitigation. The model belongs to the class of hurdle models, i.e., it is a two-part model that specifies a Bernoulli likelihood for landslide occurrence and a Gaussian likelihood for landslide log-sizes given a positive occurrence. The observations (presence/absence for the Bernoulli likelihood and log-sizes for the Gaussian likelihood) are assumed to be conditionally independent given a Gaussian latent process that drives the trends, dependencies and nonstationary patterns observed in the data. The latent process is characterised by a linear predictor with a flexible additive structure allowing the incorporation of covariates and spatial effects. Specifically, we rely on generalised additive models (GAMs) to flexibly model the covariates' influence using fixed and random effects, which are also known as linear and non-linear effects in terms of their influence. Spatial dependence between locations is characterised using a Gaussian Process with Matérn covariance structure. The Matérn family of covariance functions is widely used in spatial statistics due to its flexible local behaviour and interesting theoretical properties (Guttorp and Gneiting 2006;Stein 2012). To allow for fast inference, we use the stochastic partial differential equation (SPDE;Lindgren et al. 2011) approach that provides accurate Markovian representations of the Matérn covariance. Under a Bayesian framework, we assume relatively weak but highly interpretable Gaussian priors for all the model components and hyperparameters involved and fit our model using the integrated nested Laplace approximation (INLA; Rue et al. 2009). INLA uses the Laplace approximation to compute posterior distributions of interest. It is developed for the vast class of latent Gaussian models and implemented in the R-INLA library (Bivand et al. 2015). For related contributions, see Opitz et al. (2018), Castro-Camilo et al. (2019 and Lombardo et al. (2021).
Throughout this paper, we will interchangeably refer to susceptibility or probability of occurrence (see Lombardo et al. 2020;Titti et al. 2021). As for the notion of hazard, this further extends the concept of susceptibility by looking at measurable quantities of hazardous processes, such as size and frequency (Guzzetti et al. 1999). However, the literature on statistical applications in landslide science has never presented a model where the susceptibility is combined with the expected size of landslides, which we propose here for the first time.
The remainder of the paper is organised as follows. Section 2 introduces the variables considered for modelling. Section 3 outlines our methodology, describing the latent Gaussian model framework and the two likelihoods that constitute our hurdle model. We also provide details on prior distributions and inference using INLA. In Sect. 4 we summarise our results. Conclusions and a critical discussion are provided in Sect. 5.

Data description
After hurricane Maria, the Global Facility for Disaster Reduction and Recovery sponsored the CHARIM (Caribbean Handbook on Risk Information Management; link here) project, intending to build a comprehensive natural hazard assessment. As part of this assessment, a team from the University of Twente (NL) mapped a large-scale landslide inventory for the whole island (see Fig. 1a), primarily using five scenes of Pléiades satellite images with a resolution of 0.5m, dated on the 23rd September and 5th October 2017, made available through UNITAR-UNOSAT (https://www.unitar.org/). Additionally, a series of Digital Globe Images were used after they were gathered for the Google Crisis Response. All the images were visually inspected by image interpretation experts. As a result, landslides were mapped as polygons, separating scarp, transport and accumulation areas, and classifying the landslides according to their types.
The inventory features 9960 landslides, with a cumulative planimetric area of 11.4 km 2 , covering 1.5 % of the Island of Dominica. To aggregate the landslide information over space, we opted for a Slope Unit partition of Dominica (see Fig. 1b). A Slope Unit (SU, hereafter) encompasses the geographic space between ridges and streamlines. In other words, they are half-basins of a given order (Amato et al. 2019) extracted by maximising areas with homogeneous aspects and are particularly suited to support landslide models for they approximate the morphodynamic behaviour of landslides well (Carrara et al. 1991). To partition the Island of Dominica into SUs, we initially used r.slopeunits, an open-source software developed by Alvioli et al. (2016). We have later refined them through manual editing in GIS, to obtain a total of 3335 SUs.
We use a two-step procedure to build our hurdle model at the SU scale. First, we assign to each SU the sum of all planimetric landslide areas intersecting the SU itself. Then, to build a binary dataset, we assign a landslide presence to SUs with positive landslide aggregated areas. Conversely, we assign a landslide absence to SUs unaffected by slope failures. Finally, to build a dataset of landslide sizes per SU, we extract SUs with positive landslide aggregated areas.
Covariates available for analysis are detailed in Table 1 (geographical) and Table 3 (geological) in the Appendix. They are a mixture of geographical and geological characteristics and constitute the morphology of Dominica's landscape. All the geographical variables are continuous, while the geological and land-use types are represented as proportions of each SU; for example, a SU could be 50% forest, 20% bare and 30% quarry.

Methods
As mentioned in Sect. 1, we aim to build a hurdle model to detect areas of higher landslide susceptibility and sizes. Our approach models the probability of observing a landslide in a SU using a Bernoulli distribution. Given that a landslide was observed, it also describes landslide log-sizes using a Gaussian distribution. We stress again that we assume that the observations (presence/absence for the Bernoulli likelihood and log-sizes for the Gaussian likelihood) are conditionally independent given a latent Gaussian structure that drives the trends, dependencies and nonstationary patterns observed in the data. Both parts of our hurdle model belong to the class of latent Gaussian models.
This section describes the general framework of latent Gaussian models and the two likelihood models. We also provide details regarding prior specification and inference based on INLA.

Latent Gaussian models
Gaussian Latent Models (GLMs) are a vast and flexible class of models that are well suited for modelling spatial data . They admit a hierarchical representation, where observations are assumed to be conditionally independent given a latent field and a set of hyperparameters. Specifically, let y be the vector of Panels c and d report the geological and land-use maps, respectively. Acronyms definition for panels c and d can be found in Table 3 Stochastic Environmental Research and Risk Assessment (2022) 36:2071-20842073 observations with components yðsÞ, s 2 S, where S is the study region, in our case the collection of all SUs. Let x with components xðsÞ, s 2 S be a latent Gaussian random field and let h 1 be a set of hyperparameters. Then, the latent Gaussian field can be hierarchically formulated as follows:

ð3:1Þ
For our hurdle modelling approach, /ðÁÞ is either the Bernoulli probability mass function for landslide susceptibility or the Gaussian density for landslide log-sizes. The vector h 1 contains hyperparameters associated with the likelihood model. So, for instance, if /ðÁÞ is the Gaussian density, then h 1 is equal to the Gaussian precision, which is the inverse of the standard deviation; see Sect. 3.2 for more details on the likelihood models. The latent Gaussian random field x has the role of describing the trends and underlying dependence in the data. It is assumed to have a Gaussian distribution with mean vector l x ðh 2 Þ and precision matrix (the inverse of the covariance matrix) Q x ðh 2 Þ, both of which are controlled by the vector of hyperparameters h 2 , which account for variability and length or strength of dependence.
The above hierarchical representation can be further simplified by assuming that yðsÞ only depends on a linear predictor gðsÞ. The linear predictor for our hurdle model has an additive structure with respect to some fixed and random effects and a term that accounts for the spatial dependence between SUs. Specifically, our linear predictor can be written as where a is an intercept and ðw 1 ðsÞ; . . .; w M ðsÞÞ > are a subset of the covariates detailed in Tables 1 and 3 . . .; f K ðÁÞg are random (or non-linear) effects defined in terms of a set of covariates ðz 1 ðsÞ; . . .; z K ðsÞÞ > . For our hurdle model, the specific form of the functions f k ðÁÞ is that of a Gaussian random walk of order 1 (RW1), defined over a binned version of our covariates. Specifically, for any continuous covariate The precision parameter s 0;k controls the strength of dependence among neighbouring covariate bins, or in other words, the level of smoothness of the random walk. Finally, the term uðsÞ in (3.2) is a zero-mean Gaussian field with a stationary Matérn covariance function (Matérn (2021), we set m ¼ 1 to reflect the rather rough nature of the underlying physical process. The link between the representations in (3.2) and (3.1) is as follows: if g ¼ ðgðs 1 Þ; . . .; gðs jSj ÞÞ > and we assume that ða; b > Þ > has a Gaussian prior, then the joint distribution of x ¼ ðg > ; a; b > ; f > Þ > is Gaussian and yields the latent field x in the hierarchical representation 3.1; see, e.g., Rue et al. (2017) for more details.

The likelihood models
We model landslide susceptibility with a Bernoulli distribution. Specifically, using the notation from Sect. 3.1, we have that yðsÞ ¼ O L ðsÞ 2 f0; 1g and /ðyðsÞ j g Bern ðsÞÞ BernðpðsÞÞ, where pðsÞ ¼ PrfO L ðsÞ ¼ 1g; see Fig. 2a for the spatial distribution of O L ðsÞ. Note there are no hyperparameters for the Bernoulli likelihood, i.e., the vector h 1 is empty. The probability pðsÞ is related to the linear predictor gðsÞ through the logit link, so pðsÞ ¼ expfg Bern ðsÞg= ð1 þ expfg Bern ðsÞgÞ.
Landslide size distribution is positively skewed, with extremely large values elongating the right tail. In cases like this, it is standard practice to use a monotonous transformation such as the natural logarithm to obtain a roughly symmetric, Gaussian-like distribution (see Lombardo et al. 2021). Given that a landslide has occurred in the SU s 2 S, we model its log-size using a Gaussian distribution. Using the notation from Sect. 3.1, we have that yðsÞ ¼ logfA L ðsÞg j O L ðsÞ ¼ 1 and /ðyðsÞ j g Gauss ðsÞ; h 1 Þ N ðg Gauss ðsÞ; s À1 Þ, where h 1 ¼ s. Figure 2b shows the spatial distribution of logfA L ðsÞg. The linear predictor is linked to the Gaussian mean via the identity link. The global precision hyperparameter s (reciprocal of the standard deviation) determines the concentration of all values yðsÞ around their mean gðsÞ, s 2 S. The linear predictors g Bern ðsÞ and g Gauss ðsÞ follow the general equation (3.2). Nonetheless, their specific form depends on the influence of the covariates in Tables 1 and  3 over the landslide susceptibilities and sizes. Therefore, we conduct variable selection using a stepwise forward procedure for most of the covariates in both parts of our hurdle model. This procedure was based on numerical techniques such as the Deviance Information Criterion and the Watanabe-Akaike information criterion (DIC and WAIC, respectively; Meyer 2014; Gelman et al. 2014), and graphical methods such as the probability integral transform (PIT; Gneiting et al. 2007) and fitted versus observed plots. Due to their definition or interpretability, some covariates were not tested this way. Instead, they were included or excluded based on expert opinion. Specifically, expert advice was considered for land-use types and lithology types (included linearly as they are represented as a proportion of a SU and therefore, their sum is constrained to 1); mean Eastness and mean Northness (included linearly as they are complementary measurements representing the sine and cosine of the aspect of a SU); mean slope (included non-linearly based on previous analysis, see, e.g., Tanyaş et al. 2022); and the standard deviations of all covariates (included linearly since nonlinear standard deviations lack reasonable interpretation). Fig. 2 a Observed landslide presence/absence data; b observed landslide planimetric area, aggregated at the SU level as the sum of all landslides; c Triangulation mesh (gray) with SU centroids (blue) used to discretise Dominica Island and fit the spatial effect uðsÞ in (3.2). The inner boundary (red) delimits the island, whereas the extension to the outer boundary (black) avoids possible boundary effects All the covariates related to the perimeter were excluded due to collinearity issues. The selected covariates and the specific form they enter into the linear predictor equation for each model are detailed in Table 2.

Prior specification
Here we describe the form of the third layer of the hierarchical representation in (3.1), i.e., the prior distributions for the likelihood hyperparameters h 1 ¼ s (the Gaussian precision) and the hyperparameters of the model components in the linear predictor, h 2 ¼ ða; b; s 0 ; r 2 ; rÞ > . Here, a is a vector of length two that contains the intercept for both likelihood models, b is a vector of length 48 that contains the fixed effects for both likelihood models and s 0 is a vector of length ten that contains the RW1 precision parameters for both likelihood models; see Table 2. The vectors r 2 and r have length 2 and contain, respectively, the marginal variances and range parameters of the Matérn covariance of uðsÞ in (3.2) for both likelihoods.
Non-informative priors are a common choice when little expert knowledge is available. We use this approach to define prior distributions over the fixed effects and the intercepts. Specifically, we chose a zero-mean Gaussian prior with a precision of 0.001 for all fixed effects and intercepts. Prior information with different strengths can also be defined using the penalised complexity (PC) prior approach ). This procedure penalises excessively complex models by placing an exponential prior on a distance to a simpler base model, which helps to stabilise the estimation. Priors then shrink model components toward their base models, thus preventing over-fitting. For the precision s of the log-size observations, we set a weakly informative prior such that the probability of observing a standard deviation (1=s) larger than the empirical standard deviation for the response is 0.01. For the precision parameters of the random walks of order 1, we set weak prior distributions where the probability that the standard deviations ð1=s 0;k Þ corresponding to SU A , SU MD and TWI are greater than 5, 0.5 and 0.1, respectively is 0.01. For the precision parameters of the remaining random walks of order 1, we set relatively weak prior distributions where the probability that the standard deviations of the random walks are greater than the empirical standard deviation of the response is 0.01. Due to our lack of prior knowledge, we argue that this choice made sense since, a priori, the effects on the linear predictor can only be interpreted relative to the likelihood and the implicit scaling in the likelihood. For both Matérn range parameters, we guided our selection using the empirical variogram and set a prior distribution where the probability that the range is smaller than 25km is 50%. Finally, for the Matérn marginal variances, we set a prior distribution where the probability that the variance is larger than 0.25 is 50%.

Inference with INLA
In a Bayesian framework, the interest lies in the joint posterior distribution of unknown parameters and hyperparameters. In the context of GLMs (see (3.1)), this joint posterior distribution can be written as  INLA uses numeric techniques to approximate the integrals in (3.5), with an additional step to approximate the marginal posterior distribution for the spatial effect uðsÞ in (3.2). Recall that uðsÞ is a zero-mean Gaussian field with a stationary Matérn covariance function. Therefore, computation of the posterior marginal distribution requires the inversion of an n Â n matrix, which becomes computationally demanding for large n. To speed up this step and make inference feasible, INLA uses the stochastic partial differential equation (SPDE) approach of Simpson et al. (2017) to approximate the Matérn covariance function with numerically convenient sparse precision matrices. Specifically, the SPDE approach constructs an approximation to the precision matrix (defined as the inverse of the covariance matrix) and finds an SPDE whose solutions have desired covariance and precision structures. Lindgren et al. (2011) show how to find these solutions by discretising the study region with a triangulation mesh and representing the stochastic process as a sum of basis functions multiplied by coefficients. They also show that these coefficients form a Gaussian Markov random field, for which methods for fast computation of precision matrices already exist (see, e.g., Krainski et al. 2018 for more details on the SPDE approach). The mesh construction is key to representing the spatial process and can affect the accuracy and speed of the estimation. This step is similar to the choice of integration points on a numerical integration method (Krainski et al. 2018); therefore, we studied the sensitivity of our estimates over changes in the triangulation mesh. We found that a mesh with 4, 600 nodes is a good compromise between the accuracy and speed of the algorithm (see Fig. 2c). INLA and the SPDE approach are conveniently implemented in the R-INLA software (Bivand et al. 2015), and we use them to obtain fast and accurate approximations of the posterior distributions of all the parameters and hyperparameters involved in our models.

Application to landslide hazard assessment in the Island of Dominica
This section presents the results of our hurdle model in terms of the statistical findings and the implications for landslide hazard assessment. For the Gaussian model, out of 19 covariates used linearly, eight appeared to be significant, which indicates that the model is 95% certain of the role (either positive or negative) of the given covariate with respect to the landslide log-size. We can see from the plot in Fig. 3   . Figure 4 displays the posterior means and corresponding 95% credible bands of the random effects for both models. We can see the highly non-linear influence of most of these covariates on landslide log-sizes and occurrences. For example, looking at the plot for landslide susceptibility, we can see that SU A has a moderate non-linear effect, with a positive effect peaking at approximately 0:67 km 2 . ELl displays a slight concave effect, with a positive effect on landslide susceptibility between 0 and 0.4 m=1000. We can also see that PRCl has a relatively linear and mild effect over landslide occurrences, while SLOl is among the most significant non-linear effects with very narrow credible intervals. As expected, steeper slopes are more at risk.

Hurdle model: main results
For landslide sizes, we can see that SU A has a moderate non-linear effect, with a negative effect for smaller areas and an increasing and eventually positive effect for larger areas. SU MD seems to have a mild concave and then convex effect on landslide log-size before becoming relatively linear. D2Sl is relatively linear when considering the credible bands, with negative effects at small distances. ELl has an overall negative effect on landslide size except for the range 0-0.4 m=1000. SLOl steepness is again among the most significant effects with very narrow credible intervals. As expected, larger slopes-those above 30°-are more at risk. Finally, TWIl has a convex relation with landslide log-size. Figure 5 show the posterior mean estimates at the mesh nodes of the spatial fields for both models. We can notice that the spatial field covers almost the entirety of the island for the Bernoulli model, except for the north-eastern coastlines. This is well in line with observations of Maria passing over Dominica and discharging such an amount of rainfall that most of the shallow material draping over the island has become prone to fail. The same cannot be observed in the case of the spatial field estimated for the Gaussian model. In this case, large landslides appear to share a high degree of spatial dependence in the centre of the island.
To assess the ability of the Bernoulli model to classify landslides correctly, the left panel of Fig. 6 shows the receiver operating characteristic (ROC) curve. It is constructed by plotting the true positive rate (TPR, also called sensitivity) versus the false positive rate (FPR, also calculated as 1-specificity). The TPR boils down to the number of unstable SUs that have been correctly classified divided by the total number of unstable SUs. As for the FPR, this measure is calculated as the number of misclassified stable SUs divided by the total number of stable SUs. The input to the Bernoulli model has a dichotomous nature, whereas the output is a probability. Therefore, to match the output to the input, one needs to use a probability threshold. The way a ROC curve is then constructed is done by selecting a large number of these thresholds and storing the TPR and FPR at each iteration. The greater the area under the ROC curve (AUROC), the better the model and its classification abilities (Zou et al. 2005). The Bernoulli model has an AUROC of 0.927, so we can safely affirm that the model does an excellent job distinguishing between both classes.
To assess the goodness of fit of the Gaussian model, the middle and right panels of Fig. 6 show a histogram of the probability integral transform (PIT) values and a plot of observed versus predicted values, respectively. PITs are commonly used to assess model calibration, i.e., the statistical consistency between the predictive distribution and the observations (Gneiting et al. 2007). If a model is wellcalibrated, then the observations should be indistinguishable from a random draw from the model. For a large number of observations, the PITs histogram serves as a tool to empirically check for uniformity. As expected, large landslide sizes seem to be underestimated by the Gaussian model. Nonetheless, the PITs in Fig. 6 are not too far from the way a histogram of uniform numbers might look, and it seems fair to assume that our model is decently calibrated. The observed versus predicted values plot shows a relatively good performance for moderate landslide sizes, although there seems to be a greater bias for landslide sizes below the first quartile.

Interpretation of the covariates' role
The advantage of a statistical model over other data-driven approaches is that the association between dependent and independent variables can be clearly interpreted. This is particularly useful for examining the geomorphological reasonability of a given model. Below, we will present a few examples of our interpretation for those covariates that behave close to our assumptions from the mechanical perspective. In this case, the interpretation should cover a dual aspect of our hurdle model, both from the landslide occurrence and size perspectives. As previously mentioned, Fig. 3 reports the posterior distribution of the regression coefficients estimated for the Bernoulli and the Gaussian models. There, among the lithological classes, ignimbrites (IGO, IPL, IGY) positively contribute to increasing the mean probability of landslide occurrence, with only IPL being significant among the three. From an interpretative standpoint, we recall that ignimbrite is a pyroclastic flow deposit largely comprised of pumice with subordinated ashes. Thus, a positive contribution to the landslide occurrence is reasonable because its mineral structure is prone to weathering, and it is well known to promote slope instabilities (e.g., Chigira and Yokoyama 2005). An   Fig. 3 Posterior means (dots) of fixed linear effects (except the intercept) with 95% credible intervals (vertical segments) for the Bernoulli and Gaussian models. The horizontal grey dashed lines indicate no contribution to the landslide occurrence or size, respectively analogous situation can be seen in gravel and alluvium materials (RGA). These are unconsolidated deposits, inherently susceptible to slope instabilities anytime the landscape evolution has set them to drape over steep topographies. Interestingly enough, most of the lithotypes selected for the Bernoulli model do not appear in the Gaussian case. There, Older Pleistocene volcanics are also reported with a significant and negative contribution to the estimated landslide sizes. Irrespective of the Bernoulli or Gaussian framework, no land-use class plays a role in explaining the landslide occurrence or size.
Another point in common between the Bernoulli and the Gaussian models is the role of SU MD=A . This covariate expresses how elongated a given SU is. The larger the value, the more stretched the SU appears, whereas the smaller the value, the more rounded the SU is. Thus, the posterior mean negative value reported for both models may indicate that narrow SUs are not only less prone to fail, but also the lesser availability of material does not allow for large landslides to be generated.
As for the non-linear effects shown in Fig. 4, there we can observe two noteworthy behaviours. SU A appears to influence both the Bernoulli and Gaussian models with a negative effect for very small SUs, while the contribution becomes increasingly positive for larger SUs. Similarly to the elongation/roundness (SU MD=A ) effect mentioned above, this can be interpreted as if larger SUs can host Posterior means (dark blue curves) of random effects with 95% pointwise credible intervals (light blue shaded area) for the Bernoulli and Gaussian models more and larger landslides compared to very small SUs. There, the conditions required for a failing mechanism to occur are much more unlikely to manifest simply because of the smaller extent. The other interesting effect corresponds to SLOl. In both models, the slope steepness behaves as a sigmoidal function. This is likely because no shallow flow-like landslide can occur at low steepness values. Conversely, at medium steepness values, we experience a sudden increase in the effect of this covariate, which is to be expected, up to the point where we reach an asymptotic level or even a decrease. A decrease can be interpreted as if very steep topographic conditions cannot host soil, which is washed away by normal erosional processes. Thus, the absence or near-absence of soil implies that no landslide can manifest or that, at best, a small one will mobilise the thin detrital layer draping over the stable bedrock.

Landslide hazard components
Figure 7 presents our findings in terms of posterior means and width of the 95% credible intervals (CI) for landslide susceptibility and sizes (in log-scale) across the island. There, large portions of the areas characterised by higher landslide probability (central north and southern coastlines) coincide with areas where larger landslide sizes are expected. Conversely, the north-eastern coast seems consistently small in landslide size and probability of occurrence.
From the uncertainty plots, we can see that areas with higher occurrence probability (0.86 to 1) have overall small 95 % CI width, implying that these areas are well estimated with relatively small uncertainty. We can also see that areas with moderate occurrence probability (0.47 to 0.86) display higher uncertainty, while areas with the smallest

Unified landslide hazard
The Bernoulli model addresses the island's susceptibility under extreme conditions, such as those induced by Maria. It estimates how prone a given location is to host slope failures, but it is blind to how large landslides may become once they initiate and propagate downhill. To compensate this limitation, the Gaussian model provides information on the expected size of landslides per slope unit, again under extreme meteorological stress. But once more, the model is blind to which slope was effectively prone to fail. So far, we have implemented and presented the results of the two models independently from each other. This framework is already more than enough to satisfy the requirements of the most accepted landslide hazard definition (see Guzzetti et al. 1999), taking aside the temporal dimension. Nevertheless, much more can be done to provide for the first time a unified version of the landslide hazard for the purely spatial context. Below we explain how we can provide a new and unified (data-driven-specific) landslide hazard assessment by combining the two elements of our proposed hurdle model. Our motivation to provide a unified framework stemmed from the need to provide end-users with spatially distributed information regarding how likely each slope would be to release specific landslide sizes if expected to be unstable under analogous extreme conditions to those brought by Maria. Small landslides should carry the least hazard; thus, they should be of little interest. Conversely, as the landslide size increases, the expected hazard should proportionally follow. So, to quantify the hazard for moderate and relatively large landslides within the range of observed landslide sizes, we use the two parts of our hurdle model to compute the probability of observing landslide sizes above the empirical 50%, 75%, 90% and 95% quantiles. Specifically, we compute: for u ¼F À1 logðA L Þ ð0:5Þ;F À1 logðA L Þ ð0:75Þ;F À1 logðA L Þ ð0:90Þ and F À1 logðA L Þ ð0:95Þ, whereF logðA L Þ is the empirical cumulative distribution function of landslide log-sizes. Note that this procedure can be performed for any landslide size u of interest. We here focus on the empirical 50%, 75%, 90% and 95% quantiles to illustrate how our model can be used to predict the probability of exceeding medium to large landslides. The exceedance probabilities PrðlogfA L ðsÞg [ uÞ in (4.1) and their uncertainties can be easily computed using posterior samples from both the Gaussian and Bernoulli models. Specifically, we follow a Monte Carlo (MC) procedure and generate N ¼ 500 samples from the posterior predictive distributions (PPDs) of the Gaussian and Bernoulli models. The PPDs account for uncertainty in the data and the model fitting. We then compute empirical estimates of the probabilities in the right side of (4.1) to obtain one estimate of PrðlogfA L ðsÞg [ uÞ. Finally, we replicate this procedure M ¼ 1000 times to obtain MC mean estimates and confidence intervals for PrðlogfA L ðsÞg [ uÞ, computed as the 2.5% and 97.5% quantiles of the MC estimates. Figure 8 shows plots of the exceedance probabilities and the width of the associated 95% confidence intervals for the four quantiles detailed above. It is interesting to observe how the exceedance probability of the landslide sizes changes from one hazard map to the other. As one should expect, the slopes prone to release median landslide sizes are quite numerous, and their number decreases towards larger landslides. To briefly touch on risk perspectives (although not explicitly integrated into our hurdle model), landslides greater than 90% of the landslide area distribution are particularly likely to occur in the southernmost sector of the island. There, the village of Berekua appears to be potentially vulnerable. The same situation can be seen slightly northwestward for the much larger settlement of Roseau. We stress here that this type of consideration would not be possible in the simple binary case, where the corresponding susceptibility map highlights most of the island as unstable (see Fig. 7a). Note that the confidence intervals in Fig. 8 measure the accuracy of the Monte Carlo approximation and should not be interpreted as a measure of the dispersion in the posterior predictive distribution of the exceedance probabilities.

Discussion and conclusion
A complete landslide hazard assessment should address three components: where, when (or how frequently) and how large landslides may be. The three components mentioned above have always been addressed separately in the scientific literature produced so far (for non-physicallybased models). In this work, we combined two of the three, leaving the temporal characteristic out of the scope because of the limited access to long time series of landslide inventories triggered in response to hurricanes within the Island of Dominica. Despite this limitation, the model we propose represents a substantial improvement with respect to traditional presence/absence susceptibility models because those are blind to the actual threat that a landslide may pose, depending on its size. We consider this achievement a stepping stone for further experimentation in the hope that one day this unified hazard line of research may impact international guidelines for disaster risk reduction. In this work, we focused on the landslide inventory generated by Maria because of the large sample size and the availability of data on a wide variety of geographical and geological variables, with no missing values. However, further effort can be made to implement an extension in space-time of our hurdle framework. In this line, we are already investigating the combination of the hurricane Maria inventory with the landslide inventory mapped after hurricane Erika in 2015.
Although the lack of more temporal replicates inhibits the applicability of complex spatio-temporal models, the statistical characterisation of both processes could provide an exciting perspective and additional means of validation on landslide occurrence and size probabilities across the island. Furthermore, studying the similarities and differences in the values for landslide size and landslide occurrence could be a promising start towards an even more comprehensive hazard map of the island.
Another potential improvement would be carried out by extending this hurdle model toward different modes of slope failures and propagation. Currently, we only model shallow flow-like landslides that either started as debris slides or flows and that generally evolved into debris flows due to the high water content of the moving mass. However, extensions of our model could be implemented to distinguish various landslide classes, including deep-seated ones, which have a completely different failure mechanisms and propagation behaviour.
Ultimately, we also envision the integration of actual rainfall data. Unfortunately, in the specific case of Maria, no reliable rainfall data is available. However, other inventories in data-rich conditions could help integrate rainfall into the modelling strategy, which in turn could enable testing for a near-real-time application of our hurdle model.
To promote repeatability and reproducibility, the codes and data to fit our hurdle model are freely available from https://github.com/BryceErin/ULHAdominica. We hope that this will help users utilise and extend the hurdle modelling approach to other geographic contexts, landslide types and triggers.