Landslide-Tsurrogate v1.0: a computationally efficient framework for probabilistic tsunami hazard assessment applied to Mayotte (France)

Landslide-Tsurrogate v1.0 is an open-source Python and MATLAB tool that helps scientists quickly estimate the tsunami hazards generated by submarine landslides. Instead of running thousands of heavy deterministic numerical simulations, the software builds surrogate models that reproduce the main results with a fraction of the computational cost. The method relies on a mathematical approach called generalized polynomial chaos expansion, which efficiently explores how uncertain landslide parameters affect tsunami generation. Users can perform sensitivity analyses, identify the most influential parameters, and quantify the variability of possible outcomes. The tool includes a Jupyter Notebook User Manual and interactive MATLAB and Jupyter Notebook interfaces, making it easy to understand the methodology, set up the surrogate simulations and visualize the results. The Landslide-Tsurrogate v1.0 model's performance is demonstrated through a real-world test case involving five zones in Mayotte (France). For this application, the surrogate models achieve convergence with only 135 deterministic simulations per zone and produce probabilistic results in less than 2 s within the user-friendly interfaces used on a basic laptop, demonstrating the computational efficiency of the approach. Beyond this example, the framework can be applied to any coastal region prone to submarine landslides. By combining physical modeling, statistical analysis, and user-friendly design, Landslide-Tsurrogate v1.0 enables faster and more transparent probabilistic tsunami hazard assessments.

1 Introduction

In the last decade, surrogate model approaches – providing fast-running approximations of complex dynamical models – have been gaining significant traction in the geoscience community, particularly within the Destination Earth framework which promotes the creation of Digital Twins. For example, Formaggia et al. (2013) built a surrogate model of basin-scale geochemical compaction, Wang et al. (2016) studied acoustic uncertainty predictions, Sraj et al. (2014) estimated the wind drag parameter forcing an ocean model, Giraldi et al. (2017) documented the propagation of earthquake ocean floor displacement uncertainty to tsunami wave parameters and Bulthuis et al. (2019) used a surrogate model to quantify the uncertainty of the multi-centennial response of the Antarctic ice sheet to climate change. Compared to traditional sampling techniques such as Monte Carlo simulations, surrogate models are markedly more efficient for uncertainty quantification. Once constructed for a particular application or geographic domain, their use requires very limited further computational resources (Knio and Le Maître, 2006; Najm et al., 2009). Furthermore, surrogate models are inherently flexible: they can be constructed based on sparse or incomplete prior knowledge and tailored to situations where observational and computational resources are limited. For instance, uniform prior distributions can be adopted in the absence of comprehensive observational data, and sparse grids can be used to reduce the number of model evaluations required during training.

Given these advantages, surrogate models are particularly well-suited for Probabilistic Tsunami Forecast (PTF) in regions where observational data are scarce and computational capacities are limited. Unlike traditional high-fidelity simulations, which can be computationally intensive, surrogate models offer a more efficient framework to approximate tsunami dynamics across a wide range of potential source scenarios. These models enable rapid propagation of uncertainties associated with various tsunami generation mechanisms – be it seismic, volcanic, landslide, or atmospheric – to estimate the probability of coastal hazard metrics exceeding critical thresholds pertinent to infrastructure planning, emergency response, and long-term risk mitigation.

Recent studies have demonstrated the efficacy of surrogate modeling in diverse probabilistic tsunami modeling contexts. For instance, Denamiel et al. (2019, 2020, 2021) developed a surrogate model approach to quantify meteotsunami hazards in the Adriatic Sea, addressing challenges posed by limited observational data and the coarse resolution of atmospheric models. In the Cascadia Subduction Zone, Guillas et al. (2018) introduced a functional emulator to approximate high-resolution tsunami simulations, facilitating efficient uncertainty quantification in tsunami hazard assessments for the Pacific Northwest. Additionally, Salmanidou et al. (2021) constructed statistical emulators, enabling high-resolution probabilistic tsunami predictions in northern Cascadia. The integration of machine learning techniques into surrogate modeling has further expanded its applicability and efficiency. For example, Brecht et al. (2025) utilized physics-informed neural networks to model tsunami inundation, offering a meshless approach that simplifies the modeling process while maintaining accuracy. Such advancements highlight the growing potential of surrogate models to serve as foundational tools for tsunami hazard assessment, especially in contexts where data limitations and computational constraints hinder the application of fully deterministic modeling strategies. Collectively, these examples underscore the versatility and effectiveness of surrogate models in probabilistic tsunami modeling, providing a computationally viable means to assess tsunami hazards across various scenarios and regions. This versatility is especially valuable in areas where tsunamis can cause damage to infrastructure and pose serious risks to the population.

A compelling example of such a setting is landslide-induced tsunamis which differ from their seismic counterparts by (i) their extremely short arrival times since the source is often located near the coast, (ii) their potentially larger initial wave heights near the source, and (iii) the absence of clear precursors detectable by traditional Early Warning Systems (EWS; Roger et al., 2024; Lemoine et al., 2020). This combination of rapid onset and high impact underscores the need for efficient hazard assessment tools capable of supporting local preparedness and emergency planning.

In this study, a surrogate modeling framework specifically designed to support probabilistic tsunami forecasting (PTF) workflows for landslide-generated tsunamis is presented: the Landslide-Tsurrogate v1.0 model. The framework combines sparse grid sampling (Smolyak, 1963; Gerstner and Griebel, 1998; Constantine et al., 2012) with generalized Polynomial Chaos Expansions (gPCE; Xiu and Karniadakis, 2002; Soize and Ghanem, 2004) in order to approximate the response of computationally expensive tsunami models while requiring only a limited number of deterministic simulations and simultaneously providing direct access to sensitivity indices and uncertainty quantification metrics. Consequently, gPCE-based surrogate modeling offers an efficient alternative to machine learning approaches that often require extensive training data and may provide limited interpretability regarding the influence of uncertain input parameters.

Importantly, this surrogate modeling strategy is not tied to a specific tsunami generation mechanism. The gPCE approach approximates the response of the numerical tsunami model as a function of a set of uncertain input parameters, and is therefore applicable to any tsunami source that can be parameterized with a finite set of inputs. The main adaptation required for different tsunami sources lies in the definition of these input variables and their probability distributions. The framework can therefore be readily adapted to other tsunami sources provided that the dominant uncertain parameters controlling tsunami source generation can be identified. For example, in the case of submarine landslides considered in this study, the stochastic variables describe the physical and geometrical properties of the landslide (e.g., location, volume, etc.). In the case of sub-aerial landslides entering the ocean, if the model does not describe the sub-aerial flow, the parameters could also include the landslide mass entering the water, the entry velocity, the impact geometry, or the impact location. Similarly, for seismic tsunamis, the stochastic variables may represent earthquake source parameters such as fault slip, rupture length, or hypocentral location. Finally, for the forecast of atmospheric tsunamis (or meteotsunamis), the methodology has already been applied by defining the parameters of pressure disturbances (e.g., amplitude, direction, speed, period, etc.) driving these events in the Adriatic Sea (Denamiel et al., 2019, 2020).

Figure 1 Location of the study area. (a, b) General maps. The rectangle in (a) gives the extent of the map in (b), and the rectangle in (b) the extent of map in (c). (c) Mayotte lagoon geomorphology including Grande Terre, Petite Terre and Piton locations. Map data: ESRI World Physical Map and © Google 2026. Source: (a, b) ESRI | Powered by Esri and (c) Imagery © Google 2026, Map data © Google 2026.

Here, this broadly applicable approach is applied in the Mayotte region, where since the beginning of the seismo-volcano crisis in 2018, one particular concern was the question of a tsunami triggered by submarine landslides linked to intense seismic activity. Mayotte is a volcanic territory located in the northern Mozambique Channel of the Indian Ocean. This archipelago, composed of the main islands of Grande Terre and Petite Terre, is characterized by a shallow submarine shelf that encloses a large lagoon (Fig. 1). Since May 2018, the island of Mayotte has registered intense seismic activities related to the birth of a large new submarine volcano 50 km offshore Petite Terre (Cesca et al., 2020; Lemoine et al., 2020; Feuillet et al., 2021). The intense seismo-volcanic crisis that has affected Mayotte since 2018, the location of earthquakes near the steep slopes surrounding the lagoon around the islands, and the construction of a new volcanic structure (Lemoine et al., 2020; Feuillet et al., 2021; Mercury et al., 2022) may trigger submarine instabilities offshore East of Mayotte. Such collapses could trigger submarine landslide-generated tsunamis with serious consequences for the coastal population as shown by recent simulations (Roger, 2019; Lemoine et al., 2020; Poulain et al., 2022; Marbœuf et al., 2025). This application in Mayotte serves both to validate the surrogate model approach and to highlight its operational potential for hazard-prone regions with limited resources.

The article is structured as follows. To ensure transparency and reproducibility, Sect. 2 presents the full mathematical framework of Landslide-Tsurrogate v1.0, which forms the scientific and computational basis of the model. Section 3 and Appendix A outline the implementation of the framework, from the selection and representation of input uncertainties to their integration within a complete tsunami hazard assessment workflow. These sections aim to provide sufficient detail for readers to produce surrogate models in different contexts. Section 4 and Appendix B then demonstrate the use of Landslide-Tsurrogate v1.0 through an application to submarine landslide-generated tsunamis in Mayotte, serving as a real-world validation of the approach. Finally, Sect. 5 presents the limits of the gPCE-based surrogate modeling approach and Sect. 6 summarizes the main conclusions and discusses potential future developments.

2 Surrogate Modeling Approach

The primary objective of this study is to develop a model capable of delivering faster-than-real-time probabilistic forecasts of tsunami hazards triggered by submarine landslides. In this context, the uncertainty associated with submarine landslide parameters needs to be fully explored through the evaluation of large ensembles of landslide–tsunami simulations. Running such ensembles with complex deterministic numerical models can be computationally prohibitive, which motivates the development of surrogate modeling approaches capable of rapidly approximating the system response across the parameter space. To address this challenge, Landslide-Tsurrogate v1.0 employs a computationally efficient surrogate modeling strategy based on generalized Polynomial Chaos Expansions (gPCE; Xiu and Karniadakis, 2002; Soize and Ghanem, 2004). Hereafter, to illustrate the mathematical concepts, Fig. 2 provides a schematic overview of the main components of the surrogate modeling approach.

Figure 2 Schematic overview of the main mathematical concepts used in Landslide-Tsurrogate v1.0: (a) definition of the problem including quantity of interest, stochastic variables and total order of the gPCE; (b) multi-indices, (c) Legendre polynomials, (d) orthogonal polynomial bases, (e) sets of abscissas and (f) number of deterministic simulations needed, depending on the total order of the gPCE; (g) comparison of the quantity of interest with the surrogate model results depending on the total order of the gPCE.

2.1 Theoretical Framework

A key advantage of the gPCE method lies in its suitability for uncertainty propagation, as statistical moments and sensitivity indices can be directly derived from the gPCE coefficients. The choice of gPCE as the mathematical framework for Landslide-Tsurrogate v1.0 is grounded in the assumption by Ernst et al. (2012) that the gPCE representations of scalar tsunami hazard quantities – such as maximum tsunami elevation, maximum tsunami speed, and time of arrival, hereafter collectively referred to as ξ_max – converge. Specifically, the mean-square error between the gPCE approximation ξ_gPCE,p and the hazard quantity ξ_max tends toward zero as the expansion order p increases. This assumption is justified provided that ξ_max has finite variance; a condition that is generally reasonable from a physical standpoint. Figure 2 illustrates this assumption for a theoretical non-physical modeled quantity of interest ξ_max (no unit) and an expansion order $p\in \mathrm{0},\mathrm{\dots },\mathrm{6}$.

Additionally, one of the main difficulties in assessing submarine landslide-tsunami hazards lies in generating plausible landslide scenarios. This challenge stems from scarce data on critical factors such as the volume, shape, location, and frequency of tsunamigenic landslides, the availability of high-resolution digital elevation models, and incomplete knowledge of landslide mechanics and material properties (Roger et al., 2024; Løvholt et al., 2020). Consequently, stochastic variables required to generate submarine landslide scenarios may vary in number r between regions and have various prior distributions (i.e., normal, uniform, etc.). However, due to the lack of observations, the stochastic variables are often assumed to follow uniform distributions – i.e., within the range of the uniform distributions, the generation of any submarine landslide is equally probable. For a set of elementary events (ω∈Ω) and $i\in \mathrm{1},\mathrm{2},\mathrm{\dots },r$, the stochastic variable prior distributions are thus given by:

$$\begin{array}{}\text{(1)}& R_i:\left\{\begin{array}{l}\mathrm{\Omega }\to \mathbb{R}\\ \mathit{\omega }\to R_i\left(\mathit{\omega }\right)\end{array}\right.,R_i\left(\mathit{\omega }\right)\sim \mathcal{U}\left(\right[a_i,b_i\left]\right)\end{array}$$

Figure 2a shows an example of Landslide-Tsurrogate v1.0 employing r=3 stochastic input variables, each defined as uniformly distributed random variables. In this illustrative case, no specific physical meaning or units are assigned to these variables.

Importantly, the stochastic variables must be mutually independent, and the bounds of each uniform distribution must be selected to encompass all conditions known to trigger landslide-generated tsunamis within a given region. The use of uniform distributions and the assumption of independence are both grounded in the maximum entropy principle. According to this principle, when limited information is available, the distribution with the highest entropy is preferred, as it is based on the fewest assumptions (Jaynes, 1957). Nevertheless, adopting this approach implies that many of the randomly generated submarine landslide scenarios may, in fact, lack the capacity to produce significant tsunamis in the area under investigation.

2.2 Generalized Polynomial Chaos Expansions

2.2.1 General Framework

The gPCE mathematical framework is based on polynomial expansions that decompose into deterministic coefficients ξ_α and multivariate orthogonal polynomial bases ${\mathrm{\Psi }}_{\mathit{\alpha }}(R_{\mathrm{1}},R_{\mathrm{2}},\mathrm{\dots },R_r)$ depending on the multi-index α. The multi-index $\mathit{\alpha }=[{\mathit{\alpha }}_{\mathrm{1}},{\mathit{\alpha }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\alpha }}_r]$ provides the polynomial degree α_i associated with the stochastic variable $R_i,i\in \mathrm{1},\mathrm{2},\mathrm{\dots },r$. Figure 2b illustrates the derivation of the multi-index α for r=3 and $p\in \mathrm{0},\mathrm{2},\mathrm{4},\mathrm{6}$.

The optimal expansion – minimizing the error for a maximum total order p such as $\left|\mathit{\alpha }\right|=\sum _{i=\mathrm{1}}^r{\mathit{\alpha }}_i\le p$ – is then given by the following truncated gPCE:

$$\begin{array}{}\text{(2)}& {\mathit{\xi }}_{\text{gPCE},p}(R_{\mathrm{1}},R_{\mathrm{2}},\mathrm{\dots },R_r)=\sum _{\left|\mathit{\alpha }\right|\le p}{\mathit{\xi }}_{\mathit{\alpha }}{\mathrm{\Psi }}_{\mathit{\alpha }}(R_{\mathrm{1}},R_{\mathrm{2}},\mathrm{\dots },R_r)\end{array}$$

In the context of Landslide-Tsurrogate v1.0, the multivariate orthogonal polynomial basis is expressed as the product of the Legendre polynomials ${\mathrm{Le}}_{{\mathit{\alpha }}_i},i\in \mathrm{1},\mathrm{2},\mathrm{\dots },r$ which are orthogonal with respect to uniform distributions – i.e., $\mathbb{E}[{\mathrm{\Psi }}_{\mathit{\alpha }},{\mathrm{\Psi }}_{\mathit{\beta }}]={\mathit{\delta }}_{\mathit{\alpha },\mathit{\beta }}\Vert {\mathrm{\Psi }}_{\mathit{\alpha }}{\Vert }^{\mathrm{2}}$ with δ_α,β the Kronecker symbol – and are given by:

$$\begin{array}{}\text{(3)}& {\mathrm{\Psi }}_{\mathit{\alpha }}(R_{\mathrm{1}},R_{\mathrm{2}},\mathrm{\dots },R_r)=\prod _{i=\mathrm{1}}^r{\mathrm{Le}}_{{\mathit{\alpha }}_i}\left[{\displaystyle \frac{(\mathrm{2}{R}_i-a_i-b_i)}{(b_i-a_i)}}\right]\end{array}$$

Figure 2c, d, illustrate the values of the Legendre polynomials ${\mathrm{Le}}_{{\mathit{\alpha }}_i}$ for ${\mathit{\alpha }}_i\in \mathrm{0},\mathrm{1},\mathrm{\dots },\mathrm{6}$ and the orthogonal polynomial bases ${\mathrm{\Psi }}_{\mathit{\alpha }}(R_{\mathrm{1}},R_{\mathrm{2}},R_{\mathrm{3}})$ for $p\in \mathrm{0},\mathrm{2},\mathrm{4},\mathrm{6}$.

The deterministic coefficients ξ_α are the projection of the landslide-generated tsunami hazards represented by ξ_max onto each polynomial basis of the gPCE:

$$\begin{array}{}\text{(4)}& {\mathit{\xi }}_{\mathit{\alpha }}={\displaystyle \frac{\mathrm{1}}{\Vert {\mathrm{\Psi }}_{\mathit{\alpha }}{\Vert }^{\mathrm{2}}}}\mathbb{E}\left[{\mathit{\xi }}_{\text{max}}{\mathrm{\Psi }}_{\mathit{\alpha }}\right(R_{\mathrm{1}},R_{\mathrm{2}},\mathrm{\dots },R_r\left)\right]\end{array}$$

However, in the absence of an analytical formulation of tsunami hazards, numerical integration methods are employed to approximate the expectation $\mathbb{E}\left[{\mathit{\xi }}_{\text{max}}{\mathrm{\Psi }}_{\mathit{\alpha }}\right(R_{\mathrm{1}},R_{\mathrm{2}},\mathrm{\dots },R_r\left)\right]$. One common method for computing multidimensional integrals is to express them as a series of nested one-dimensional integrals using Fubini's theorem (also known as the tensor product rule). However, this technique leads to an exponential increase in the number of required function evaluations – and thus, in the presented context, the deterministic landslide-tsunami simulations – as the number of dimensions grows, a phenomenon known as the curse of dimensionality. To address this issue, three main alternative methods have been developed to mitigate the computational burden: the Monte Carlo, Sparse Grids and Bayesian Quadrature methods.

2.2.2 Pseudo-Spectral Approximation (PSA)

Within Landslide-Tsurrogate v1.0, the deterministic coefficients ξ_α are computed using a Smolyak-based sparse-grid approach (Smolyak, 1963): the non-intrusive Pseudo-Spectral Approximation (PSA) method (Gerstner and Griebel, 1998; Constantine et al., 2012). One of the key advantages of PSA is that the error decreases rapidly as the number of needed deterministic simulations increases (Fornberg and Sloan, 1994). In this context, the truncated generalized Polynomial Chaos Expansion (gPCE) can be reformulated to yield the main equation of the Landslide-Tsurrogate v1.0 model:

$$\begin{array}{}\text{(5)}& \begin{array}{r}{\mathit{\xi }}_{\text{PSA},p}(R_{\mathrm{1}},R_{\mathrm{2}},\mathrm{\dots },R_r)=\sum _{max(\mathrm{0},p-r+\mathrm{1})\le \left|\mathit{\alpha }\right|\le p}(-\mathrm{1}{)}^{p-\left|\mathit{\alpha }\right|}\\ \left(\begin{array}{c}r-\mathrm{1}\\ p-\left|\mathit{\alpha }\right|\end{array}\right)\sum _{\mathit{\beta }}{C}_{\left|\mathit{\alpha }\right|,\mathit{\beta }}{\mathrm{\Psi }}_{\mathit{\beta }}(R_{\mathrm{1}},R_{\mathrm{2}},\mathrm{\dots },R_r)\end{array}\end{array}$$

For each α from $max(\mathrm{0},p-r+\mathrm{1})\le \left|\mathit{\alpha }\right|\le p$, the multi-index $\mathit{\beta }=[{\mathit{\beta }}_{\mathrm{1}},{\mathit{\beta }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\beta }}_r]\in \mathbf{A}$ stores all the polynomial combinations needed to calculate the tensor product $\sum _{\mathit{\beta }}{C}_{\left|\mathit{\alpha }\right|,\mathit{\beta }}{\mathrm{\Psi }}_{\mathit{\alpha }}(R_{\mathrm{1}},R_{\mathrm{2}},\mathrm{\dots },R_r)$ with:

$$\begin{array}{}\text{(6)}& C_{\left|\mathit{\alpha }\right|,\mathit{\beta }}={\displaystyle \frac{Q_{\mathit{\alpha }}^{(\mathrm{1}:r)}\left({\mathit{\xi }}_{max}{\mathrm{\Psi }}_{\mathit{\beta }}\right)}{\Vert {\mathrm{\Psi }}_{\mathit{\beta }}{\Vert }^{\mathrm{2}}}}\end{array}$$

and the multivariate quadratic rule:

$$\begin{array}{}\text{(7)}& Q_{\mathit{\alpha }}^{(\mathrm{1}:r)}\left({\mathit{\xi }}_{max}{\mathrm{\Psi }}_{\mathit{\beta }}\right)=\sum _{\mathit{\alpha }}{w}_{\mathit{\alpha }}^{(\mathrm{1}:r)}{\mathit{\xi }}_{max}\left(s_{\mathit{\alpha }}^{(\mathrm{1}:r)}\right){\mathrm{\Psi }}_{\mathit{\beta }}\left(s_{\mathit{\alpha }}^{(\mathrm{1}:r)}\right)\end{array}$$

The complete set of abscissas $s_{\mathit{\alpha }}^{(\mathrm{1}:r)}$ – i.e., the combination of parameters sampling the defined ranges of the stochastic variables – and their associated weights $w_{\mathit{\alpha }}^{(\mathrm{1}:r)}$ required for computing the multivariate quadrature rule are then obtained through the tensor products:

$$\begin{array}{}\text{(8)}& \left\{\begin{array}{l}{s}_{\mathit{\alpha }}^{(\mathrm{1}:r)}=s_{{\mathit{\alpha }}_{\mathrm{1}}}^{\left(\mathrm{1}\right)}\otimes s_{{\mathit{\alpha }}_{\mathrm{2}}}^{\left(\mathrm{2}\right)}\otimes \mathrm{\cdots }\otimes s_{{\mathit{\alpha }}_r}^{\left(r\right)}\\ w_{\mathit{\alpha }}^{(\mathrm{1}:r)}=w_{{\mathit{\alpha }}_{\mathrm{1}}}^{\left(\mathrm{1}\right)}\otimes w_{{\mathit{\alpha }}_{\mathrm{2}}}^{\left(\mathrm{2}\right)}\otimes \mathrm{\cdots }\otimes w_{{\mathit{\alpha }}_r}^{\left(r\right)}\end{array}\right.\end{array}$$

where the univariate sampling of each stochastic variable R_i is defined by the S_i abscissas $s_{{\mathit{\alpha }}_i}^{\left(i\right)}={\left\{s_k^{\left(i\right)}\right\}}_{k=\mathrm{1}}^{S_i}$ and the corresponding weights $w_{{\mathit{\alpha }}_i}^{\left(i\right)}={\left\{w_k^{\left(i\right)}\right\}}_{k=\mathrm{1}}^{S_i}$, for $i=\mathrm{1},\mathrm{2},\mathrm{\dots },r$. In practice, the unique set of abscissas extracted from $s_{\mathit{\alpha }}^{(\mathrm{1}:r)}$ forms the ensemble of input parameters that will be fed to the deterministic landslide-tsunami model. For the illustrative case, Fig. 2e shows examples of different set of unique abscissas $s_{\mathit{\alpha }}^{(\mathrm{1}:\mathrm{3})}$ selected from the distributions of the r=3 stochastic variables depending on the total order $p\in \mathrm{0},\mathrm{2},\mathrm{4},\mathrm{6}$.

2.3 Quadrature Rule

The selection of an appropriate quadrature rule for sampling the stochastic variables is critical, as it determines both the accuracy of the surrogate models and the total number of deterministic landslide-tsunami simulations to be performed. Among the most commonly used quadrature rules for uniform distributions are the Gauss–Legendre (Abramowitz and Stegun, 1965), Clenshaw–Curtis (Clenshaw and Curtis, 1960), and Gauss–Patterson (Patterson, 1968) rules.

The Gauss–Legendre quadrature is the standard n-point rule for uniform distributions, offering exact integration for any polynomial of degree up to 2n−1. However, a key limitation of the Gauss–Legendre rule is that its abscissas are not nested – i.e., the points used in the n-point rule are entirely different from those in the (n−1)-point rule, preventing reuse of previous evaluations.

Table 1 Number of abscissa S depending on the level of different nested quadrature rules.

In contrast, nested quadrature rules allow the abscissas of a lower-order rule to be included in higher-order ones, which significantly enhances computational efficiency. The Clenshaw–Curtis and Gauss–Patterson rules are the most frequently used nested quadrature rules for uniform distributions. The main difference between these two lies in their precision and growth rate: the Gauss–Patterson rule typically provides nearly double the accuracy and point growth compared to the Clenshaw–Curtis rule. In nested rules, the concept of level is introduced to denote the quadrature rule where the abscissas are nested. At level l, the number of abscissas required is 2^l+1 for Clenshaw–Curtis and ${\mathrm{2}}^{l+\mathrm{1}}-\mathrm{1}$ for Gauss–Patterson (Table 1). To make Gauss–Patterson rules more efficient – by limiting the rapid growth in the number of abscissas – the introduction of successive quadrature levels can be delayed. This delayed Gauss–Patterson rule (Holtz, 2011) uses a reduced set of abscissas (Table 1) while preserving high integration accuracy, typically at least 2l−1. Consequently, this approach significantly reduces the number of required quadrature points without sacrificing the quality of the surrogate model.

In Landslide-Tsurrogate v1.0, only the Gauss–Patterson (GP) and delayed Gauss–Patterson (DGP) quadrature rules are implemented, as they offer greater precision than the Clenshaw–Curtis rule. Despite this theoretical advantage, both Gauss-Patterson and Clenshaw–Curtis rules have demonstrated comparable performance in practice (Trefethen, 2008). However, following the findings of Holtz (2011) and Denamiel et al. (2019), the DGP rule not only reduces the number of deterministic simulations required but might also enhances the accuracy of the surrogate models. Therefore, the use of the DGP rule is strongly recommended for sampling the definition intervals of the stochastic variables in Landslide-Tsurrogate v1.0. Figure 2f presents the number of deterministic simulations N needed to create gPCE-PSA-based surrogate models with the DGP quadrature rule for the modeled quantity of interest ξ_max depending on the total order $p\in \mathrm{0},\mathrm{\dots },\mathrm{6}$.

Finally the GP and DGP rules solve univariate quadrature problem for stochastic variables following the standard distribution $Z\sim \mathcal{U}\left(\right[-\mathrm{1},\mathrm{1}\left]\right)$. The relation between standard distributions and the distributions used in Landslide-Tsurrogate v1.0 can be expressed as follows, for a set of elementary events (ω∈Ω) and $i\in \mathrm{1},\mathrm{2},\mathrm{\dots },r$:

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{l}{Z}_i:\left\{\begin{array}{l}\mathrm{\Omega }\to \mathbb{R}\\ \mathit{\omega }\to Z_i\left(\mathit{\omega }\right)\end{array}\right.,Z_i\sim \mathcal{U}\left(\right[-\mathrm{1},\mathrm{1}\left]\right)\\ \\ R_i=\left(\frac{a_i+b_i}{\mathrm{2}}\right)+\left(\frac{b_i-a_i}{\mathrm{2}}\right)Z_i,R_i\sim \mathcal{U}\left(\right[a_i,b_i\left]\right)\end{array}\right.\end{array}$$

Similarly the relationship between the abscissa and weights calculated for the standard distributions and the ones used in Landslide-Tsurrogate v1.0 can be written, for $i\in \mathrm{1},\mathrm{2},\mathrm{\dots },r$:

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{l}\left\{\begin{array}{l}{s}_{{\mathit{\alpha }}_i}^{\left(i\right)}={\left\{s_k^{\left(i\right)}\right\}}_{k=\mathrm{1}}^{S_i}\\ \\ w_{{\mathit{\alpha }}_i}^{\left(i\right)}={\left\{w_k^{\left(i\right)}\right\}}_{k=\mathrm{1}}^{S_i}\end{array}\right.,Z_i\sim \mathcal{U}\left(\right[-\mathrm{1},\mathrm{1}\left]\right)\\ \\ \left\{\begin{array}{l}\left(\frac{a_i+b_i}{\mathrm{2}}\right)+\left(\frac{b_i-a_i}{\mathrm{2}}\right)s_{{\mathit{\alpha }}_i}^{\left(i\right)}\\ \\ \left(\frac{b_i-a_i}{\mathrm{2}}\right)w_{{\mathit{\alpha }}_i}^{\left(i\right)}\end{array}\right.,R_i\sim \mathcal{U}\left(\right[a_i,b_i\left]\right)\end{array}\right.\end{array}$$

Landslide-Tsurrogate v1.0 is built upon the full system of equations (Eqs. 1 to 10) presented in the previous subsections, while the numerical simulations required to compute the quadrature rule can be generated using any suitable landslide-tsunami model. For a given region and a given set of stochastic variables, the surrogate models are built at specific locations $m\in \mathrm{1},\mathrm{2},\mathrm{\dots },M$ generally situated along the coastline. Figure 2g illustrate the behavior of a surrogate model in function of the total order $p\in \mathrm{0},\mathrm{2},\mathrm{4},\mathrm{6}$. However, the convergence of surrogate models is slower or even impossible when different sets of abscissas (i.e., input parameters of the landslide-tsunami model) are linked to the same value of ξ_max (e.g., 0 m for the maximum elevation in the inundation zone when there is no inundation). Consequently, in this study, the locations where the surrogate models are built are deep enough to remain submerged – i.e., outside of the inundation zone – ensuring that results of maximum elevation, maximum speed and time of arrival are continuous – i.e., defined with different values for all the deterministic landslide-tsunami simulations.

2.4 Convergence and Accuracy

Once generated with Landslide-Tsurrogate v1.0, two different evaluation steps are needed to demonstrate the surrogate model skill in estimating the tsunami hazards – maximum elevation, maximum speeds and time of arrival. These steps are the convergence and accuracy of the surrogate models based on the comparison between the tsunami hazards predicted with the N-member ensembles of deterministic simulations ${\left\{{\mathit{\xi }}_{\text{max}}^{m,n}\right\}}_{n=\mathrm{1}}^N$ and surrogate model results ${\left\{{\mathit{\xi }}_{\text{PSA},p}^{m,n}\right\}}_{n=\mathrm{1}}^N$ at the chosen locations $m\in \mathrm{1},\mathrm{2},\mathrm{\dots },M$ and for the total order $p\in \mathrm{1},\mathrm{2},\mathrm{\dots },P$.

Firstly, the convergence of the surrogate models is analyzed with the normalized root-mean-square errors at each total order $p\in \mathrm{1},\mathrm{2},\mathrm{\dots },P$ and for each chosen location $m\in \mathrm{1},\mathrm{2},\mathrm{\dots },M$:

$$\begin{array}{}\text{(11)}& e_p^m={\displaystyle \frac{\sum _{n=\mathrm{1}}^N{\left({\mathit{\xi }}_{max}^{m,n}-{\mathit{\xi }}_{\text{PSA},p}^{m,n}\right)}^{\mathrm{2}}}{\sum _{n=\mathrm{1}}^N{\left({\mathit{\xi }}_{max}^{m,n}-{\mathit{\xi }}_{\text{PSA},\mathrm{0}}^{m,n}\right)}^{\mathrm{2}}}}\end{array}$$

The convergence criteria is fulfilled when the errors (Eq. 11) follow $\underset{p\to P_{\mathrm{t}}}{lim}{e}_p^m=\mathrm{0}$. In this context, the total order of the gPCE selected to build the surrogate models is P_t. Figure 2g represents this convergence of a gPCE-PSA-based surrogate model ξ_PSA,p with $p\in \mathrm{0},\mathrm{2},\mathrm{4},\mathrm{6}$ towards the illustrative quantity of interest ξ_max when p increases.

Secondly, evaluating the accuracy of the surrogate models theoretically requires an ensemble of deterministic simulations that (1) are sampled using a Monte Carlo method and (2) are large enough – typically on the order of 10⁴ runs – to produce statistically robust results. In practice, however, this approach is not viable, since the primary motivation for constructing surrogate models is to reduce the computational burden associated with running a large number of deterministic simulations. Additionally, at higher total polynomial orders, reusing the same simulations used to build the surrogate risks overfitting. Within the Landslide-Tsurrogate v1.0 framework, which relies on nested sparse grids, a more practical strategy involves constructing surrogate models at a total order higher (P_t+1) than the one selected for convergence (P_t). The additional deterministic simulations ΔN generated at this higher order being independent – i.e., not used to build the surrogate models – and effectively sampling the full range of the stochastic parameter space can therefore serve as a robust basis for accuracy assessment. This is shown for the illustrative case in Fig. 2e, f where, for example, the additional set of abscissas and, hence, the additional deterministic simulations, $\mathrm{\Delta }N=+\mathrm{48}+\mathrm{72}=+\mathrm{120}$ needed to go from p=4 (in yellow) to p=6 (in red) clearly covers the entire space of definition of the stochastic variables $R_{\mathrm{1}},R_{\mathrm{2}},R_{\mathrm{3}}$.

2.5 Sensitivity Analysis

A key challenge is to determine which stochastic variable R_i for $i\in \mathrm{1},\mathrm{2},\mathrm{\dots },r$ most significantly influence the variability of landslide-generated tsunami hazards at the selected locations of interest $m\in \mathrm{1},\mathrm{2},\mathrm{\dots },M$. In other words, the goal is to identify the parameters to which the hazard estimates produced by the surrogate models are most sensitive.

This sensitivity analysis is carried out using variance-based decomposition methods (ANOVA), in particular Sobol' sensitivity indices (Sobol', 2001; Saltelli, 2002), which can be efficiently derived from the coefficients of the gPCE representation (Crestaux et al., 2009). Following the notations introduced in the previous subsections, the Sobol' indices are derived from the gPCE coefficients for which $\mathrm{1}\le \left|\mathit{\alpha }\right|\le P_{\mathrm{t}}$ and β∈A, β≠0. These gPCE coefficients can be expressed as:

$$\begin{array}{}\text{(12)}& \widehat{{\mathit{\xi }}_{\mathit{\beta }}^m}=(-\mathrm{1}{)}^{P_{\mathrm{t}}-\left|\mathit{\alpha }\right|}\left(\begin{array}{c}r-\mathrm{1}\\ P_{\mathrm{t}}-\left|\mathit{\alpha }\right|\end{array}\right)C_{\left|\mathit{\alpha }\right|,\mathit{\beta }}^m\end{array}$$

At each chosen location $m\in \mathrm{1},\mathrm{2},\mathrm{\dots },M$, the variability of the tsunami hazards D^m and the total Sobol' indices for each stochastic variable $S_{R_i}^{\mathrm{T},m}$ with $R_{i}i\in \mathrm{1},\mathrm{2},\mathrm{\dots },r$ are given by:

$$\begin{array}{}\text{(13)}& \left\{\begin{array}{l}{S}_{R_i}^{\mathrm{T},m}=\frac{D_i^{\mathrm{T},m}}{D^m}\\ \\ D^m=\text{Var}\left[{\mathit{\xi }}_{\text{PSA},Pt}^m\right]=\sum _{\begin{array}{c}\mathit{\beta }\in \mathbf{A}\\ \mathit{\beta }\ne \mathrm{0}\end{array}}{\left(\widehat{{\mathit{\xi }}_{\mathit{\beta }}^m}\right)}^{\mathrm{2}}\\ \\ D_i^{\mathrm{T},m}=\sum _{\mathit{\beta }\in {\mathbf{A}}_i^{\mathrm{T}}}{\left(\widehat{{\mathit{\xi }}_{\mathit{\beta }}^m}\right)}^{\mathrm{2}},{\mathbf{A}}_i^{\mathrm{T}}\in \left\{\mathit{\beta }\in \mathbf{A}:{\mathit{\beta }}_i>\mathrm{0}\right\}\end{array}\right.\end{array}$$

The total Sobol' indices (Eq. 13) quantify the relative contribution (ranging from 0 to 1) of each stochastic variable to the variability of landslide-generated tsunami hazard predictions at the selected locations. Since total indices account for both individual effects and all possible interactions with other variables, their sum may exceed 1 – reflecting that interactions are counted multiple times. In practice, if a total sensitivity is negligible (e.g., below 0.01), the corresponding stochastic variable can be fixed to a deterministic value without significantly impacting the resulting hazard distributions.

3 Landslide-Tsurrogate v1.0

Landslide-Tsurrogate v1.0 is a computational tool designed to support the development of gPCE-based surrogate models for regions affected by landslide-generated tsunamis. It includes a comprehensive set of functions (i.e., MATLAB and Python codes), a step-by-step User Manual provided as a Jupyter Notebook to guide users through the surrogate modeling workflow, and example user-friendly interfaces implemented in both Jupyter Notebook and MATLAB, all illustrated with the Mayotte test case (https://doi.org/10.5281/zenodo.17519408; Denamiel, 2025).

The construction of surrogate models with Landslide-Tsurrogate v1.0 follows seven building blocks or steps, as illustrated in Fig. 3. This section discusses in detail the strengths and limitations of the methodological choices implemented in Landslide-Tsurrogate v1.0. For the steps that involve user-defined inputs (marked with ^*), the various options available to users are also outlined to highlight the versatility of the tool.

Figure 3 Overview of the seven Landslide-Tsurrogate v1.0 building blocks enabling faster-than-real-time Probabilistic Tsunami Forecast (PTF) based on surrogate models. These building blocks are: (a) Step 1: User Specifications, defining the stochastic variables; (b) Step 2: Input Parameters, providing the inputs for deterministic simulations; (c) Step 3: Deterministic Simulations, using complex numerical models; (d) Step 4: Data Transfer and Formatting; (e) Step 5: Surrogate Models, created with the deterministic simulations; (f) Step 6: Evaluation, assessing surrogate model convergence and performance; and (g) Step 7: User-Friendly Interface, generating the PTF outputs based on user inputs.

1. User Specifications^*

   This step involves defining the stochastic variables used to represent the uncertain characteristics of potential submarine landslides. These variables describe the physical and geometric properties of the landslide source and can be grouped into three main categories: (i) soil parameters, such as cohesion, friction angle, and density; (ii) landslide geometry, including the thickness, length, width, and orientation of the shapes used to represent the landslide volume initially released; and (iii) initial conditions, such as the mean sea level, which may vary due to tidal fluctuations and long-term sea-level rise. In principle, the initial released mass and failure characteristics of a submarine landslide could be determined using mechanical slope stability models that incorporate measured soil properties (e.g., density, cohesion, internal friction angle), as well as hydrogeological conditions, structural discontinuities, and fault networks. However, such detailed information is rarely available for offshore environments. Consequently, many landslide properties remain poorly constrained and represent a major source of epistemic uncertainty in landslide-tsunami hazard assessment.

   Within the Landslide-Tsurrogate v1.0 framework, these uncertainties are explicitly represented through stochastic variables. Potential landslide locations, for example, may be identified using geomorphological indicators such as steep seabed slopes, sediment accumulations, previously mapped instability features or areas experiencing high seismicity. Landslide volumes, geometric parameters, and rheological properties can then be sampled within plausible ranges derived from regional geological knowledge, previous studies, or expert judgment. Similarly, mean sea level can be treated as an additional stochastic variable to capture the nonlinear interactions between tsunami propagation, tidal variability, and long-term sea-level rise.

   By representing the uncertain landslide characteristics using a limited set of stochastic variables (e.g., location, volume, friction angle, and sea level), the dimensionality of the problem can be kept manageable while still capturing the dominant sources of uncertainty affecting landslide-generated tsunamis. These variables define the stochastic parameter space used to generate the deterministic landslide-tsunami simulations required to construct the surrogate models. The stochastic space is subsequently sampled using sparse quadrature rules, allowing the surrogate models to efficiently approximate the response of the deterministic simulations across the range of possible landslide scenarios.

   Finally, it is important to note that the definition of stochastic variables and their ranges is inherently site dependent. Some uncertainties reflect natural variability in environmental conditions (aleatory uncertainty), while others arise from limited knowledge of subsurface properties and landslide triggering mechanisms (epistemic uncertainty). The proposed framework does not prescribe a unique parameterization but instead provides a flexible structure for propagating these uncertainties through the deterministic simulations.

2. Input Parameters

   This step is critical in constructing the surrogate models within the Landslide-Tsurrogate v1.0 framework. It involves generating the input parameter sets for the deterministic simulations.

   The input parameters correspond to the stochastic variables used to describe the physical and geometric properties of the landslides. Once these stochastic variables and their probability distributions have been defined in the previous step, a sparse Gauss–Patterson quadrature is constructed up to a user-defined total order P. The quadrature nodes – representing collocation points in the stochastic parameter space – are mapped from the canonical space to the physical space. Each combination of parameters defines a deterministic landslide scenario that is then simulated using a numerical landslide-tsunami model.

   The total number of deterministic simulations N is determined by three main factors: (i) the number of stochastic variables r, (ii) the selected total polynomial order P used to construct the surrogate model, and (iii) the quadrature strategy adopted (Gauss–Patterson or delayed Gauss–Patterson). In practice, increasing either the number of stochastic variables or the total order increases the number of collocation points required to accurately represent the stochastic space. The quadrature rule determines how efficiently this space is sampled.

   For example, using six stochastic variables (r=6) and a selected total order P=5, the number of required simulations is N = 10 625 with the Gauss–Patterson (GP) sparse grid and N=1889 using the delayed Gauss–Patterson (DGP) approach. The DGP method therefore provides a substantial reduction in the number of deterministic simulations while maintaining comparable surrogate accuracy.

   The choice of P represents a trade-off between surrogate model convergence/accuracy and computational cost. Higher orders improve the ability of the surrogate model to represent nonlinear interactions between parameters but require more deterministic simulations. Because each deterministic simulation involves running a full numerical tsunami model, this stage determines the computational cost of the framework.

   For this reason, surrogate models can be constructed incrementally by truncating the polynomial expansion to lower orders (e.g., P=4) and only generating additional deterministic simulations to reach the polynomial order needed for the convergence of the surrogate models (e.g., P_t=6).

3. Deterministic Simulations^*

   This step primarily involves selecting and running a numerical model to simulate the landslide-generated tsunamis based on the input parameters generated in the previous step. The surrogate models developed in Landslide-Tsurrogate v1.0 are trained on deterministic tsunami simulations, which therefore constitute a fundamental component of the overall framework. The choice of the numerical model used to generate these simulations should be guided by the physical complexity of the system considered. Depending on the characteristics of the landslide source, the bathymetry, and the relevant hydrodynamic processes, different levels of model complexity may be required. In some cases, depth-averaged shallow-water models may provide adequate accuracy for hazard assessment, while in other settings more sophisticated approaches, such as dispersive or fully three-dimensional models, may be necessary to properly represent the generation mechanism and the resulting wave field. The surrogate model does not impose constraints on the deterministic solver and can in principle emulate the outputs of any numerical model, provided that the training dataset sufficiently samples the parameter space and that the underlying simulations capture the relevant physical processes. Consequently, the reliability of the probabilistic hazard estimates ultimately depends on the adequacy of the deterministic model used to generate the training simulations.

   Among the state-of-the-art numerical models available within the submarine landslide community, several fully nonlinear and dispersive solvers have been developed with varying levels of complexity, computational cost, and intended applications. NHWAVE (Ma et al., 2012) is a three-dimensional, non-hydrostatic model particularly suited for simulating submerged landslides and the resulting tsunami generation. It provides high accuracy in capturing the wave generation and early propagation phases, but its 3D nature makes it computationally expensive, especially for large domains or probabilistic assessments. Geowave (Watts et al., 2003) is a coupled model that integrates a Boussinesq wave solver with a solid slide motion solver, allowing for flexible representation of landslide kinematics but relying on simplified assumptions for slide rheology. BING (Løvholt et al., 2005) represents the landslide as a viscous fluid and is often coupled with long-wave tsunami solvers such as COMCOT (Liu et al., 1998) or COULWAVE (Lynett et al., 2002; Lynett and Liu, 2004), striking a balance between physical realism and computational cost. Finally, Landslide-HySEA (Macías et al., 2021a, b) offers a unique advantage through its multi-layer non-hydrostatic shallow water formulation, which is especially important for capturing the vertical structure and dynamics of submarine mass movements. This multi-layer approach enables the realistic representation of both submerged and subaerial landslide processes and their interaction with the water column, leading to improved accuracy in wave generation. The HySEA suite is also GPU-accelerated, which significantly reduces computational time despite its complexity, making it particularly attractive for regional hazard assessments that require large ensembles of simulations.

   Each of these models has trade-offs: NHWAVE and Landslide-HySEA offer higher fidelity but at greater computational cost, while models like Geowave and BING are more efficient but rely on simplifying assumptions. In addition, the setup of the simulations including domain size, horizontal resolution, number of layers (e.g., MHWAVE and Landslide-HySEA), Manning friction coefficient, CFL number, dry land threshold value, etc., will greatly impact the accuracy of deterministic simulations and, hence, of the surrogate model-based PTF. The balance between numerical cost and accuracy is thus constraint by the availability of computational resources which, practically, plays a crucial role on the choice of the numerical model and simulation setup.

4. Data Transfer and Formatting^*

   The deterministic simulations required to build accurate surrogate models can generate a large volume of data. As such, transferring and formatting this dataset becomes a critical step. However, data management largely depends on the computing environment used (e.g., HPC vs. local server). Landslide-Tsurrogate v1.0 only requires the user to provide a specifically formatted file that must contain the maximum elevation, maximum speed, and time of arrival of the N tsunami waves extracted at the selected M locations, where the surrogate models are to be constructed.

5. Surrogate Models

   This is the main building block of Landslide-Tsurrogate v1.0. Here, the deterministic coefficients of the gPCE decomposition are computed based on the selected quadrature rule and the faster-than-real-time probabilistic tsunami hazard assessments are provided at the chosen M locations. A couple of important aspects of the methodological choices used to create the surrogate models should be highlighted.

   Firstly, in Landslide-Tsurrogate v1.0, the logarithm of the tsunami hazard quantities is used to build the surrogate models instead of their raw values. This offers several advantages, particularly when dealing with quantities that span multiple orders of magnitude or exhibit skewed distributions. Log-transforming these quantities tends to normalize their distribution, reduce heteroscedasticity (i.e., variability that changes with the magnitude), and suppress the influence of extreme outliers. This results in smoother response surfaces that are more amenable to polynomial approximations such as those used in generalized Polynomial Chaos Expansions (gPCE). Additionally, the log-transformation improves the numerical stability and convergence of the surrogate models, especially in regions where the hazard quantity may approach zero. However, care must be taken when applying the log-transform, as it requires all input values to be strictly positive. This condition is generally met by selecting offshore locations outside the inundation zone. Additionally, the log-transform can compress extreme values – slightly reducing the model's sensitivity to rare but intense events. Overall, applying the logarithmic transformation – used by default in Landslide-Tsurrogate v1.0 – remains a valuable and robust strategy for enhancing surrogate model performance in tsunami hazard assessments.

   Secondly, a notable characteristic of the generalized Polynomial Chaos Expansion (gPCE) approach is that the surrogate models are constructed independently at each location of interest. That is, the gPCE is built pointwise, meaning that the statistical relationships between input parameters and tsunami hazard quantities – such as maximum elevation, speed, or arrival time – are captured separately for each site without considering spatial correlations. This is in contrast to other methods such as Gaussian Process Regression (GPR; e.g., Rasmussen and Williams, 2005) or Deep Neural Networks (DNNs; e.g., Luo et al., 2025), which can exploit shared structure or dependencies across space and potentially model multiple outputs jointly. While the pointwise nature of gPCE ensures computational simplicity and interpretability (Xiu and Karniadakis, 2002; Sudret, 2008), it may limit the method's ability to account for spatial coherence in the tsunami signal, particularly in complex coastal environments where hazard responses are spatially correlated.

   In the context of landslide-generated tsunami hazard assessment, the gPCE choice made in Landslide-Tsurrogate v1.0 is particularly well-suited. The underlying physics of tsunami generation and propagation are smooth and governed by well-understood equations, making them ideal candidates for polynomial-based surrogates. Moreover, the relatively low number of stochastic input parameters – carefully defined by the user – keeps the dimensionality within the range where gPCE remains both accurate and efficient. The method's ability to provide an interpretable analytical structure, its natural convergence monitoring through sparse quadrature rules, and its rapid evaluation speed make it a robust and practical choice. While alternative surrogate approaches may offer advantages in high-dimensional or highly nonlinear contexts, gPCE provides a highly effective balance of performance, transparency, and computational efficiency for the specific demands of probabilistic tsunami hazard assessment.

6. Evaluation

   The evaluation step is a key component of Landslide-Tsurrogate v1.0, aimed at assessing the convergence (necessary to define the gPCE truncation order), accuracy, and sensitivity (needed to exclude certain stochastic variables) of the surrogate models.

   A key strength of the gPCE framework used in Landslide-Tsurrogate v1.0 lies in its ability to assess convergence at total order P_t through hierarchical polynomial orders. By comparing the surrogate model outputs at a given total order P with those at the next order P+1, especially when using the Gauss–Patterson nested quadrature rules, users can quantitatively evaluate the stability and accuracy of the expansion and define the order P_t for which the convergence is reached. This nested structure allows for computational efficiency since the nodes from the lower-order quadrature are reused in higher-order computations, minimizing redundant simulations. Moreover, this approach provides a natural mechanism for adaptive refinement, as discrepancies between orders highlight regions where the surrogate may be under performing. However, this method also has limitations. The convergence between successive orders may be slow or misleading in cases where the model response is highly non-linear or discontinuous in the stochastic space, particularly in the presence of rare but extreme events. Additionally, while nested grids improve efficiency, the cost of running additional deterministic simulations for higher orders can become prohibitive as the dimensionality increases. Nonetheless, convergence checks remain a practical and valuable diagnostic tool for ensuring surrogate model robustness within the gPCE framework.

   In Landslide-Tsurrogate v1.0, the practical approach used to assess the accuracy of gPCE surrogate models is to validate the surrogate built at total order P_t using the deterministic simulations corresponding to the quadrature nodes at order P_t+1. One key strength of this method is its efficiency: it avoids the additional computational burden of calculating all coefficients required for a full order P_t+1 surrogate, while still providing a meaningful diagnostic of the surrogate's generalization capability. This approach acts as a cross-validation scheme, testing how well the surrogate interpolates or extrapolates to unseen, but structured, parameter combinations. Moreover, because the P_t+1 quadrature includes more nodes and often better captures variability in the stochastic space, it can reveal weaknesses in the surrogate built at order P_t, especially in regions where the model response is highly sensitive. However, this methodology has its limitations. The validation points are not randomly distributed, but follow the deterministic structure of the quadrature grid, which might bias the assessment – particularly if certain regions of the input space are poorly covered by the higher-order nodes. Also, since the evaluation is based only on a subset of the full parameter space, it may overlook local deficiencies or provide an overly optimistic view of surrogate performance in regions where higher-order interactions dominate. Furthermore, this method does not capture errors related to the truncation of the polynomial chaos expansion itself, which a full surrogate at order P_t+1 might better expose. Still, when used judiciously, this validation strategy offers a valuable compromise between computational cost and model fidelity.

   Finally, using total Sobol' sensitivity indices to discard stochastic variables is a practical strategy for reducing the complexity of surrogate models in Landslide-Tsurrogate v1.0. The main strength of this approach lies in its ability to quantify the overall influence of each input variable – including interactions – on the output variance, thus offering a solid basis for identifying parameters with negligible impact. Discarding such variables can significantly lower the computational cost of the surrogate construction without sacrificing accuracy. However, a key limitation is that Sobol' indices are global measures (i.e., they reflect the properties of the surrogate models across all the scenarios generated to cover the interval of definition of the stochastic variables), and their estimation can be sensitive to sampling strategy and noise. Additionally, a variable with low total influence on the entire domain might still be critical in localized regions or under specific conditions. Hence, while total Sobol' indices provide valuable guidance, decisions to eliminate variables should be made with caution and, ideally, confirmed by localized or scenario-specific analysis.

7. User-friendly Interface^*

   To facilitate the practical use of surrogate models by a wide range of end users – including scientists, engineers, and decision-makers – developing user-friendly interfaces is essential. Several approaches can be adopted depending on the target audience and deployment environment. A Python-based GUI offers flexibility, integration with scientific libraries (e.g., NumPy, Matplotlib), and cross-platform compatibility. It is well-suited for research and prototyping, but may require users to install dependencies and be comfortable with a scripting environment. A MATLAB-based GUI benefits from a cohesive environment and strong plotting capabilities, making it particularly attractive in academic settings where MATLAB is already in use. However, it depends on licensed software, limiting accessibility and long-term sustainability. In contrast, web-based services using frameworks like Flask (Python), Streamlit, or full-stack applications (e.g., with React or Django) can reach a broader audience and require no local installation – users only need a browser. These solutions are ideal for sharing models with stakeholders or policy users, but can be more complex to develop and maintain, especially if real-time performance or secure data handling is needed. Ultimately, the choice of interface depends on balancing accessibility, development resources, and the technical proficiency of the end users.

In practice, Landslide-Tsurrogate v1.0 is implemented following the technical description and workflow presented in Appendix A based on the seven steps (or building blocks) described above.

4 Mayotte Test Case

The Mayotte test case demonstrates the application of the Landslide-Tsurrogate v1.0 model by following the sequential steps presented as building blocks in Fig. 3 and as a workflow pipeline in Fig. A2.

4.1 Steps 1 and 2: User Specifications and Input Parameters

Mayotte's submarine geomorphology is strongly shaped by its volcanic origin. Steep submarine slopes offshore transition abruptly from shallow shelf (<9°) to flank slopes reaching 25–60° locally, segmented by deep canyons and gullies (Lemoine et al., 2020; Feuillet et al., 2021). The eastern submarine slope of Mayotte would be probably prone to slope failure because it is composed of accumulation of carbonate and volcanic sediments on steep slopes. These mass-wasting events could be triggered by a combination of factors: steep slopes, weak layers of eroded volcanic ash or clay, hydrostatic pressure from groundwater, and seismic shaking linked to ongoing seismic swarms and active volcanic degassing near submarine edifices (Lemoine et al., 2020; Feuillet et al., 2021). Sultan et al. (2023) analyzed physical properties and mineralogy and performed dynamic triaxial tests on 25 sediment cores offshore of the eastern side of Mayotte to understand the hazards related to earthquake-induced submarine liquefaction. They show that the main parameter controlling the liquefaction potential offshore of Mayotte is the presence of low-density layers with high calcite content accumulating along the slope during lowstands.

Additionally, in the context of the on-going seismo-volcanic activity, Lemoine et al. (2020), Poulain et al. (2022) and Marbœuf et al. (2025) already implemented numerical simulations of potential submarine landslides and associated tsunamis in Mayotte. First, they showed that the most impacting submarine landslide scenarios were located on a large portion of the transition between the lagoon and offshore zones, presenting the steepest slopes and could have a volume that varied greatly depending on these locations. Second, they also highlighted the uncertainties around the friction angles in the landslide rheology and their high impact on the granular flow and deposit. Finally, of all the free parameters involved in their model, they demonstrated that these friction angles have the strongest impact on the wave field, compared, for example, to the Manning coefficient and the interlayer friction coefficient. Consequently, based on these previous studies, in this first attempt to apply the Landslide-Tsurrogate v1.0 framework in Mayotte, the number of stochastic variables is limited to three: location, volume and friction angle.

Figure 4 (a) Computational domain and 30 m resolution topo-bathymetry of Mayotte. Landslides locations used to build the surrogate models along (b) the Piton (red dots), Piton Offshore (white dots) and South Piton (white triangles) zones and (c) the Mayotte South (black dots) and Mayotte Offshore (white dots) zones. Dashed white lines represent the isolines at steepest slopes.

In practice, within the Landslide-Tsurrogate framework, Mayotte submarine landslides are pragmatically located along the isolines of steepest slope (extracted from high-resolution bathymetry) close to the seismic swarm related to the seismo-volcanic activity. They serve as coherent proxies for collapse-prone zones in the absence of additional comprehensive subsurface data. As illustrated in Fig. 4, five distinct zones have been selected in the Mayotte region based on the isolines of steepest slope: Piton, Piton Offshore, Piton South, Mayotte South, and Mayotte Offshore. For the volume of the submarine landslides, volumes between 1–200 × 10⁶ m³ have been selected based on (i) the works of Lemoine et al. (2020) and Poulain et al. (2022), (ii) the bathymetric data and (iii) modeling studies of landslides in shallow marine environments (Løvholt et al., 2015). The largest volumes $V=[\mathrm{1},\mathrm{200}]$ × 10⁶ m³ are used for the Piton location while $V=[\mathrm{1},\mathrm{100}]$ 10⁶ m³ are chosen for the deeper zones of Piton Offshore and Mayotte Offshore. Smaller landslide volumes of 1–50 × 10⁶ m³ are selected for the Piton South and Mayotte South zones that are associated with lower scarps volumes. The rheology and associated friction angles that should be used in models of submarine landslides are still a debate. However, previous works calibrated the friction coefficient that could be used in the Mayotte's case (Poulain et al., 2022; Marbœuf et al., 2025). Following these studies, a conservative range of friction angles [4,12]° is used to simulate potential submarine landslides in Mayotte. Such small friction angles are also necessary to reproduce the dynamics and deposits of submarine landslides involving large volumes (Poulain et al., 2023). The origin of this small dissipation in large landslides is however a challenging question (Poulain et al., 2022). Furthermore, the selected small values are aimed to describe the worst-case scenarios as often considered in tsunami hazard assessments, since low friction angles lead to higher velocities and runout distances (Løvholt et al., 2015).

In summary, for each of the five selected zones, the Landslide-Tsurrogate v1.0 framework is implemented with three stochastic variables: the along-isoline distance from the northernmost point (i.e., the point with the maximum latitude; D), the volume of the submarine landslide (V), and the friction angle (δ). These variables are defined using uniform probability distributions based on the range of physically plausible values derived from geological and geomorphological constraints, and expert knowledge:

• Piton (12.76° S, 45.31° E):

  $$\begin{array}{}\text{(14)}& \left\{\begin{array}{l}D\sim \mathcal{U}\left(\right[\mathrm{0},\mathrm{919}\left]\mathrm{m}\right)\\ V\sim \mathcal{U}\left(\right[\mathrm{1},\mathrm{200}]\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{\mathrm{3}})\\ \mathit{\delta }\sim \mathcal{U}\left(\right[\mathrm{4},\mathrm{12}\left]\mathit{°}\right)\end{array}\right.\end{array}$$

• Piton Offshore (12.77° S, 45.31° E):

  $$\begin{array}{}\text{(15)}& \left\{\begin{array}{l}D\sim \mathcal{U}\left(\right[\mathrm{0},\mathrm{1190}\left]\mathrm{m}\right)\\ V\sim \mathcal{U}\left(\right[\mathrm{1},\mathrm{100}]\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{\mathrm{3}})\\ \mathit{\delta }\sim \mathcal{U}\left(\right[\mathrm{4},\mathrm{12}\left]\mathit{°}\right)\end{array}\right.\end{array}$$

• Piton South (12.81° S, 45.30° E):

  $$\begin{array}{}\text{(16)}& \left\{\begin{array}{l}D\sim \mathcal{U}\left(\right[\mathrm{0},\mathrm{5045}\left]\mathrm{m}\right)\\ V\sim \mathcal{U}\left(\right[\mathrm{1},\mathrm{50}]\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{\mathrm{3}})\\ \mathit{\delta }\sim \mathcal{U}\left(\right[\mathrm{4},\mathrm{12}\left]\mathit{°}\right)\end{array}\right.\end{array}$$

• Mayotte South (12.90° S, 45.26° E):

  $$\begin{array}{}\text{(17)}& \left\{\begin{array}{l}D\sim \mathcal{U}\left(\right[\mathrm{0},\mathrm{6071}\left]\mathrm{m}\right)\\ V\sim \mathcal{U}\left(\right[\mathrm{1},\mathrm{50}]\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{\mathrm{3}})\\ \mathit{\delta }\sim \mathcal{U}\left(\right[\mathrm{4},\mathrm{12}\left]\mathit{°}\right)\end{array}\right.\end{array}$$

• Mayotte Offshore (12.89° S, 45.30° E):

  $$\begin{array}{}\text{(18)}& \left\{\begin{array}{l}D\sim \mathcal{U}\left(\right[\mathrm{0},\mathrm{3568}\left]\mathrm{m}\right)\\ V\sim \mathcal{U}\left(\right[\mathrm{1},\mathrm{100}]\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{\mathrm{3}})\\ \mathit{\delta }\sim \mathcal{U}\left(\right[\mathrm{4},\mathrm{12}\left]\mathit{°}\right)\end{array}\right.\end{array}$$

For the maximum total order corresponding to P=6 and the delayed Gauss-Patterson (DGP) sparse grid sampling selected for the Mayotte test case, the number of combinations of location, volume and friction angle that define the input parameters of the deterministic simulations is N=207 per zone. This number results directly from the DGP sparse quadrature construction for r=3 stochastic variables and a polynomial expansion of total order P=6 (Eqs. 14 to 18). In practice, the surrogate models are constructed using a training dataset consisting of N=135 deterministic simulations per zone corresponding to P_t=5, while the remaining 72 simulations are used as an independent testing dataset to assess the predictive performance and robustness of the surrogate model. Consequently, a total of 1035 deterministic simulations are required to construct and validate the surrogate models across the five defined zones. The choice of P_t=5 represents a compromise between surrogate accuracy and computational cost and depends on the convergence of the surrogate models that will be further discussed in Sect. 4.3.2. Lower orders may not adequately capture nonlinear interactions between landslide parameters, while higher orders would significantly increase the number of required deterministic simulations. For tsunami applications with a limited number of stochastic variables (typically r≤5), polynomial orders P_t∈4…6 provide a suitable balance between accuracy and computational efficiency (e.g., Denamiel et al., 2019, 2020).

More generally, the number of deterministic simulations required by the Landslide-Tsurrogate framework depends primarily on the number of stochastic variables and the selected polynomial order. In practical applications, it is recommended to (i) limit the number of stochastic variables to the most influential physical parameters and (ii) select a polynomial order that captures the dominant nonlinear behavior without excessively increasing the number of deterministic simulations. The computational cost associated with these deterministic simulations and the parallelization strategy adopted for the Mayotte test case are discussed in Sect. 4.2.1.

4.2 Step 3: Deterministic Simulations

4.2.1 Deterministic Model Setup

The topo-bathymetric dataset used in this study is derived from the compilation presented by Lemoine et al. (2020), which combines multiple data sources with spatial resolutions ranging from 1 to 100 m. These datasets were merged to construct an initial Digital Elevation Model (DEM) with a horizontal resolution of 10 m. For the deterministic tsunami simulations, this DEM was resampled to a 30 m resolution by removing rows and columns from the raster grid. Keeping in mind the total number of 1035 simulations to be performed, this resolution was selected as a compromise between spatial accuracy and computational efficiency while maintaining a realistic representation of the coastal and nearshore bathymetry. The computational domain covers a 21×25.5 km area defined by the local coordinate bounds $[x_{min},x_{max}]\times [y_{min},y_{max}]=[\mathrm{520}\mathrm{075},\mathrm{541}\mathrm{075}]\times [\mathrm{8}\mathrm{570}\mathrm{890},\mathrm{8}\mathrm{596}\mathrm{390}]$ in the EPSG:4471 (RGM04/UTM zone 38S) coordinate system. The computational domain and the 30 m Mayotte's bathymetric data are displayed on Fig. 4 left panel. Middle and right panels of Fig. 4 show all the submarine landslide locations and the associated isolines of steepest slopes (dashed white lines) for the five considered zones: Piton (red dots on the middle panel), Piton Offshore (white dots on the middle panel), Piton South (white triangles on the middle panel), Mayotte South (black dots on the right panel) and Mayotte Offshore (white dots on the right panel). The coastline from the SHOM Histolitt dataset is represented in all panels by a bold black line.

All the simulations are performed with Multilayer HySEA, the layer-averaged version of the HySEA model family (Macías et al., 2021a, b) described in Appendix B. Multilayer HySEA couples the landslide dynamics and the generated water waves through a friction term between landslide and water layers. Multilayer HySEA has been shown to be well suited to simulate landslide generated tsunamis in Mayotte. Indeed, a multilayer configuration is required to properly describe the effect of the strong topography variations related to the presence of a coral reef while keeping a reasonable computational cost compared to full 3D models. This enables performing the numerous simulations required for the surrogate model with the parameters set according to Marbœuf et al. (2025). The main features of all the simulations are given by:

• Final time t_f=1800 s with CFL = 0.55

• Non-hydrostatic pressures with 4 vertical layers as recommended in Marbœuf et al. (2025) for Mayotte

• Open boundary conditions with a sponge layer of 4 cells absorbing the waves

• μ(I) rheology with a grain size diameter d=1 m, a solid volume fraction Φ=0.8 and the first friction angle δ as a stochastic variable: see Eqs. (B3) and (B4)

• Density ratio $r={\mathit{\rho }}_{\mathrm{f}}/{\mathit{\rho }}_{\mathrm{s}}=\mathrm{0.5}$ where ρ_f=1000 kg m⁻³ and ρ_s=2500 kg m⁻³

• Variable Manning coefficient ranging from 0.025 to 0.07 s m${}^{-\mathrm{1}/\mathrm{3}}$ offshore (Lemoine et al., 2020; Marbœuf et al., 2025)

• A mass/water friction of m_f=0.004 m s⁻¹

The second friction angle is calculated thanks to ${\mathit{\delta }}_{\mathrm{2}}-\mathit{\delta }=\mathrm{10}$° (Pouliquen and Forterre, 2002; Poulain et al., 2023; Marbœuf et al., 2025).

Finally, to quantify the computational performance of the Landslide-Tsurrogate v1.0 framework for the Mayotte test case, the cost of the deterministic simulations required to construct the surrogate models using Multilayer HySEA is addressed. Each deterministic simulation required approximately 50 min of runtime when executed on 1 NVIDIA GPU A100 card (40GB) of the S-CAPAD/DANTE platform at IPGP, France. Because the simulations associated with the quadrature nodes are independent, the training dataset was generated using a parallel workflow in which multiple deterministic simulations were executed simultaneously. A total of 1035 simulations were required to construct and validate the surrogate models across the five landslide zones considered in this study. When distributed across 6 GPU cards, the total wall-clock time required to generate the full dataset was approximately 5.6 d, whereas the equivalent sequential runtime would be approximately 36 d. This parallelization strategy substantially reduces the computational burden associated with constructing the surrogate models and makes the framework tractable for large parameter studies or probabilistic hazard analyses involving multiple potential landslide sources.

4.2.2 Landslide volumes

In Multilayer HySEA, the initial moving mass is defined from two input grids: (1) the topo-bathymetry of Mayotte, and (2) the same topo-bathymetry modified by adding the initial failure volume. For each of the 1035 simulations, these two files are required.

Figure 5 Illustration of three failure surfaces (a, b, c) and the corresponding mass thicknesses (d, e, f) created by the algorithm, superimposed to the topo-bathymetry at 30 m resolution: (a, d) simulation 124 from Piton Offshore (V=88.84 × 10⁶ m³), (b, e) simulation 0 from Mayotte South (V=25.5 × 10⁶ m³) with the cross section AA' displayed in Fig. 6 and (c, f) simulation 114 from Mayotte Offshore (V=71.84 × 10⁶ m³).

Figure 6 Submarine failure surface of simulation 0 from the Mayotte South zone. (a) The 30 m bathymetry with the failure surface, the box ℬ from step 1, the ellipsoidal footprint, its semi-axis l₁, l₂ and the rotation angle θ involved in step 3. (b) 3D visualization of the paraboloid h_e from step 3, the plane z=0 and the ellipsoidal footprint formed by the intersection. (c–d) Different steps of the algorithm plotted on the cross section AA' shown in Fig. 5.

To generate the failure surface automatically, we developed an algorithm that modifies the 10 m DEM of Mayotte based on two stochastic parameters: the along-isoline distance D, which determines the coordinates (x_c,y_c) of the failure center, and the target released volume V. The shape of the failure surface is controlled by three geometric parameters $(l_{\mathrm{1}},l_{\mathrm{2}},r)$:

• l₁ and l₂ are the semi-axes of the ellipsoidal footprint along and across the local slope,

• rl₁ controls the maximum depth (thickness) of the failure.

These parameters define a paraboloid representing the sliding mass. Figure 5 shows different failure surfaces and sliding masses generated by the algorithm. Note that the shapes of the failure surfaces differ at each location and depend on the local bathymetry.

To match a prescribed volume V, the algorithm iteratively adjusts $(l_{\mathrm{1}},l_{\mathrm{2}},r)$ by minimizing

$$f={\displaystyle \frac{|V-V_{\mathrm{c}}|}{V}},$$

where V_c is the computed volume from the generated surface. The minimization uses the Nelder-Mead method (scipy.optimize.minimize), with convergence tolerances of 0.001 on f and 0.01 on the parameters. The construction of the failure surface (illustrated in Fig. 6 for the simulation 0 of the Mayotte South zone) follows these steps:

1. A square box ℬ, centered at (x_c,y_c), defines the local computation domain (Fig. 6a).

2. A plane 𝒫 is fitted to the bathymetry inside ℬ (Fig. 6c).

3. A paraboloid thickness h_e(x,y), centered at (x_c,y_c), is created:

   $$h_e(x,y)=r{l}_{\mathrm{1}}\left[\mathrm{1}-{\left({\displaystyle \frac{x^{\prime }}{l_{\mathrm{1}}}}\right)}^{\mathrm{2}}-{\left({\displaystyle \frac{y^{\prime }}{l_{\mathrm{2}}}}\right)}^{\mathrm{2}}\right],$$

   where $(x^{\prime },y^{\prime })$ are translated and rotated coordinates (angle θ) aligned with the local slope direction (Fig. 6a, b and d).

4. The surface 𝒮 (height z_s) is defined by the points ${\mathit{x}}_{\mathrm{s}}={\mathit{x}}_{\mathrm{p}}-h_e\mathit{n}$ where x_p are the points of the plane 𝒫 and n its normal (Fig. 6e).

5. The initial bathymetry (height z_b) is locally lowered following the surface 𝒮 to form the failure surface (height z_d), i.e. $z_{\mathrm{d}}=min(z_{\mathrm{b}},z_{\mathrm{s}})$ (Fig. 6f).

6. The removed volume is computed as

   $$V_{\mathrm{c}}=\mathrm{\Delta }x\mathrm{\Delta }y\sum (z_{\mathrm{b}}-z_{\mathrm{d}}),$$

   where $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{10}$ m is the grid size.

The parameter bounds are set to $l_{\mathrm{1}},l_{\mathrm{2}}\in [\mathrm{40},\mathrm{800}]$ m and $r\in [\mathrm{0.05},\mathrm{0.2}]$, ensuring that the generated surface is smooth at 10 m resolution and not excessively deep (to avoid trapping material).

For computational efficiency, the resulting 10 m DEMs are later downsampled to 30 m resolution, with a volume discrepancy below 1.7 %.

4.3 Steps 4 and 5: Surrogate Models and Evaluation

This section aims to illustrate how the Landslide-Tsurrogate v1.0 model operates and, hence, only the results of the surrogate models implemented along the shores of Petite Terre are presented hereafter. The corresponding results for Grande Terre are not discussed in the article but are provided via the user-friendly interfaces.

4.3.1 Locations

To provide PTF along all vulnerable Petite Terre coastal areas, the surrogate models are built pointwise (i.e., with no spatial correlation between the sites) at targeted depths of 1.0–1.5 m outside of the inundation zone. These target locations are distributed along the Petite Terre coastline in areas where inundation maps are required for risk reduction, emergency preparedness, and evacuation planning. Practically, the spatial discretization used in the deterministic simulations leads to a total of M=211 coastal target points (Fig. 7). To simplify the analysis, the surrogate models are numbered and divided into four areas: east coast (1–51), south-west coast (52–113), north-west coast (114–151 and 190–211), and lagoon (152–189). Practically, the surrogate models are built using the outputs of the 1035 deterministic tsunami simulations initiated from the five previously defined submarine landslide-prone zones (Piton, Piton Offshore, Piton South, Mayotte South, and Mayotte Offshore). For each landslide scenario, tsunami propagation is simulated, and relevant hazard metrics such as maximum tsunami elevation, maximum tsunami speed and time of arrival are extracted at the coastal target points.

Figure 7 Locations and associated numbers of the 211 surrogate models built along the east, south-west, north-west coasts and the lagoon of Petite Terre. Source: Source: Imagery © Google 2026, Map data © Google 2026.

This process enabled the construction of efficient surrogate models of tsunami hazards that capture the complex, non-linear relationships between submarine landslide source characteristics and coastal impacts along the Petite Terre coastline.

4.3.2 Evaluation

The convergence, accuracy and sensitivity of the 211 surrogate models built along the Petite Terre coastline are presented in Figs. 8–10.

1. Convergence

   The convergence analysis of the 211 surrogate models built for Petite Terre (Fig. 8) provides important insights into the reliability of the generalized Polynomial Chaos Expansion (gPCE) methodology when applied to landslide-generated tsunami hazards. Across the three key tsunami variables – maximum elevation, maximum speed, and time of arrival – the results highlight both the strengths and the limitations of the surrogate modeling approach in capturing the complex dynamics of landslide-induced tsunamis in the Piton, Piton Offshore, Piton South, Mayotte South, and Mayotte Offshore zones.

   For maximum elevation, the surrogate models display robust convergence across nearly all zones and all surrogate model locations. Normalized errors ${\left\{e_p^m\right\}}_{p=\mathrm{0}}^{\mathrm{6}}$, $m\in \mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{211}$ (Eq. 11) decrease rapidly with increasing gPCE order p, with most locations m achieving normalized errors below 5 %. This demonstrates the strong capacity of the surrogate models to approximate the spatial distribution of maximum wave heights, even in dynamically complex environments. The convergence of maximum speed is less uniform and somewhat slower. While some zones – particularly Piton and Piton Offshore – achieve acceptable error levels at higher gPCE orders, localized bands of moderate errors persist, particularly for surrogate models number 25 and 27 which are located in a protected bay. This heterogeneity may be attributed to small-scale variations in bathymetry and shoreline geometry that strongly influence flow velocities (e.g., via refraction and reflection). This uneven convergence indicates that while surrogate models provide a good first approximation of flow speed, site-specific refinements may be required for accurate hazard characterization in sensitive locations.

   The time of arrival emerges as the most difficult variable to approximate using the surrogate models. With the exception of most of the east-coast locations (1–51) and the northern most locations of the north-west coast (189–211) for the three Piton zones and the Mayotte South zone (only east-coast), errors generally decrease more slowly compared to elevation and speed, with many surrogate model locations retaining high normalized errors (>60 %) even at the highest gPCE orders tested. This underscores the sensitivity of arrival time to small uncertainties in landslide initial conditions and wave propagation processes. Additionally, it should be noted that, in HySEA, the time of arrival is set to −1 when the wave height remains below 0.01 m. For the construction of the surrogate models, realistic values were required; therefore, the −1 values were replaced with 30 min, corresponding to the full length of the simulation. While this substitution ensured consistency in the dataset, it may have affected the generation of surrogate models in locations minimally impacted by the tsunami waves (e.g., the west-coast and the lagoon). This suggests that while surrogate models can offer useful probabilistic estimates of arrival times, their predictive skill remains limited in certain contexts, necessitating caution if used for early-warning purposes.

   Overall, although the convergence patterns are not uniformly perfect across all zones, all surrogate model locations and all variables, the results clearly demonstrate that the gPCE method achieves near-optimal performance by total order p=5. The marginal improvement observed at total order p=6 does not significantly enhance convergence, indicating that increasing the polynomial order beyond p=5 provides diminishing returns. This finding underscores that a total order of P_t=5 strikes a practical balance between computational efficiency and model accuracy, ensuring reliable surrogate model representations of tsunami hazards without incurring unnecessary complexity.

2. Accuracy

   The scatter plots in Fig. 9 assess the predictive skill or accuracy of the surrogate models built using the gPCE truncated at total order P_t=5. Importantly, the surrogate results are compared against the 72 additional deterministic simulations needed to reach the total order p=6. These deterministic simulations can be considered as independent as they are not used in building the surrogate models. This approach provides a robust validation of the surrogates by testing their ability to generalize beyond the training dataset.

   For maximum elevation, the surrogate models demonstrate excellent predictive accuracy across all five zones. The data points cluster tightly around the 1 : 1 reference line, with minimal dispersion, confirming that elevation – a relatively smooth and integrated measure of tsunami intensity – is reliably reproduced by the surrogates even in unseen scenarios. The predictive performance for maximum speed is slightly less consistent, with more noticeable scatter around the reference line, particularly in Piton South and Mayotte South. This behavior is expected, as speed is a more sensitive and localized variable influenced by fine-scale bathymetric features and flow dynamics. Despite this, the surrogate models capture the overall trends well, with the densest regions of the scatter (Fig. 9, occurrences between 20 and 50, in red) still following the 1 : 1 relationship. For time of arrival, the agreement remains strong overall, though some divergence appears for Piton Offshore, Piton South and Mayotte Offshore. As explained previously, this variability may partly reflect the methodological choice to replace non-physical HySEA values (originally set to −1 for wave heights below 0.01 m) with the full simulation length of 30 min. Such replacements could affect the accuracy at sites with negligible tsunami impact, yet the bulk of the results remain closely aligned with the deterministic benchmarks.

   In summary, the comparison demonstrates that surrogate models built at total order P_t=5 not only converge internally but also generalize effectively to independent simulations, providing reliable predictions of tsunami elevation, speed, and arrival time along the coastal areas of Petite Terre.

3. Sensitivity

   Finally, Fig. 10 presents the total sensitivity indices (Eq. 13) of the surrogate models for maximum elevation, maximum speed, and time of arrival of landslide-generated tsunamis across the 211 Petite Terre locations. The decomposition highlights the relative importance of the three uncertain input parameters – landslide location, volume, and friction angle – within each of the five study zones.

   Across all zones and most locations, landslide volume emerges as the dominant contributor to the output variance of the maximum elevation and flow speed, with its sensitivity index typically approaching or exceeding 0.75. This result is consistent with the established understanding that the magnitude of the failing mass strongly controls the intensity of landslide-generated tsunamis. Rather than representing a new physical insight, this outcome provides an important consistency check, confirming that the surrogate modeling framework correctly reproduces the expected physical sensitivities of the system.

   In contrast, the landslide location shows a more variable influence depending on the region and the considered metric, reflecting the role of local bathymetric conditions and coastal geometry in modulating tsunami propagation and coastal impacts. In the Piton South and Mayotte South zones, location accounts for a substantial fraction of the sensitivity, especially for maximum speed and time of arrival, reflecting the strong influence of coastal geometry and bathymetric gradients on tsunami propagation and impact timing. The other zones exhibit a reduced location sensitivity, where volume largely governs the results except for the time of arrival.

   Finally, The friction angle consistently contributes the least to the variance across all zones and metrics, with values generally below 0.1. Its role appears minor compared to volume and location, though it may still introduce subtle effects on flow dynamics in specific sites. Interestingly, particularly for the speed and time of arrival, total sensitivities exceed 1, indicating non-linear interaction effects among the input parameters and demonstrating that location and volume combine to drive complex arrival-time variability.

   Overall, the analysis demonstrates that while landslide volume is the primary control on tsunami hazard characteristics, the importance of location cannot be neglected, particularly in areas where coastal complexity modulates local hazard intensity. The friction angle, though secondary, still contributes non-negligibly in sensitive locations.

Figure 8 Convergence of the 211 Petite Terre surrogate models of maximum elevation, maximum speed and time of arrival of the landslide-tsunamis. For each of the Piton, Piton Offshore, Piton South, Mayotte South and Mayotte Offshore zones, the convergence is presented as the normalized error (Eq. 11) in function of the total order of the truncated gPCE decomposition.

Figure 9 Accuracy of the 211 Petite Terre surrogate models of maximum elevation, maximum velocity and time of arrival of the landslide-tsunamis. For each of the Piton, Piton Offshore, Piton South, Mayotte South and Mayotte Offshore zones, the accuracy is presented as scatter plots of the surrogate model results derived from the gPCE truncated at the total order P_t=5 (based on 135 simulations) as a function of the 72 independent deterministic simulations needed to reach a total order p=6. The dashed black line represents the 1 : 1 reference line.

Figure 10 Total sensitivity of the 211 Petite Terre surrogate models of maximum elevation, maximum speed and time of arrival of the submarine landslide-tsunamis. For each of the Piton, Piton Offshore, Piton South, Mayotte South and Mayotte Offshore zones, the total sensitivity is presented as bar plots.

This evaluation (Figs. 8–10) highlights a consistent relationship between sensitivity structure, convergence, and accuracy of the surrogate models. For maximum elevation and speed, in zones where landslide volume dominates the response (e.g., Piton, Piton Offshore, Mayotte Offshore), the reduced effective dimensionality (i.e., the fact that the surrogate model results mostly depend of the landslide volume) leads to faster convergence (Fig. 8) and higher predictive accuracy when tested against independent simulations (Fig. 9). Conversely, in zones where the influence of location is stronger – particularly for speed and time of arrival in Piton South and Mayotte South – the surrogate models exhibit greater scatter around the 1 : 1 reference line (Fig. 9) and require higher orders to stabilize (Fig. 8), reflecting the increased complexity of the tsunami response. It is noteworthy that the Piton South and Mayotte South zones encompass the largest geographical extents within the study area and are characterized by pronounced slope heterogeneity. Such geomorphological variability can increase the complexity of the surrogate model parameter space, thereby hindering model convergence. Furthermore, the substantial distances separating potential landslide sources within these zones imply that tsunami arrival times will be more sensitive to the location of the initiating event. Thus, the observed sensitivity patterns directly explain why some variables, locations, and zones achieve better convergence and accuracy than others, reinforcing the physical interpretation of the surrogate model framework.

4.3.3 Step 7: Mayotte Landslide Probabilistic Tsunami Forecast (PTF)

An important aspect of the Mayotte submarine landslide test case is the emphasis on user-friendly interfaces providing faster-than-real-time forecasts (Fig. A3) that allow both researchers and decision-makers to interact with probabilistic tsunami hazard data efficiently.

However, the physical constraints of the Mayotte setting highlight the short delay between the landslide and the impact on the coast (a few minutes in some cases). Consequently, running deterministic simulations in real time, even with rapid computational tools, may be challenged by the speed of tsunami onset. Once trained using the deterministic simulations, the Landslide-Tsurrogate v1.0 framework enables rapid exploration of the stochastic parameter space describing potential submarine landslides. In this study, the probabilistic hazard associated with 1000 stochastic realizations of landslide parameters at 211 coastal locations along Petite Terre can be evaluated in less than 2 seconds on any laptop while the same number of deterministic simulations would run in about 6 d across 6 GPU cards. This capability makes it possible to characterize the statistical variability of tsunami impacts across a large number of plausible landslide scenarios, which would be computationally impractical using only full deterministic simulations.

The resulting probabilistic simulations (Fig. 11) indicate that tsunamis generated in the Piton zone can reach the airport and nearby settlements in less than three minutes for a significant subset of the explored scenarios. This result highlights the extremely short warning times associated with nearby submarine landslide sources and illustrates how the surrogate framework enables rapid evaluation of arrival-time distributions across a wide range of possible landslide parameters. Such short warning times strongly limit the feasibility of traditional evacuation strategies and emphasize the need for complementary adaptation measures. These may include vertical evacuation shelters, reinforced structures, pre-identified safe zones, clear evacuation protocols, and automated alert systems – strategies that have been studied in Mayotte (Leone et al., 2023) and implemented in other tsunami-prone regions such as Japan and Indonesia.

Figure 11 LandslideProbabilistic Tsunami Forecast (PTF) derived from 1000 conditions (or samples) for the Petite Terre Airport area and considering submarine landslides with volumes of V = 100–200 × 10⁶ m³ located at the tip of the Piton zone. The map represents the median of the maximum tsunami elevation distributions along the Petite Terre coastline while the probabilistic density distributions are displayed based on the surrogate models of maximum elevation, maximum speed and time of arrival located along the eastern coast of the airport area. Source: Imagery © Google 2026, Map data © Google 2026.

By integrating fast, accessible computational tools with proactive adaptation measures, Landslide-Tsurrogate not only advances scientific understanding but also provides actionable solutions for emergency preparedness. This combination ensures that hazard assessments translate directly into timely and informed decisions, ultimately enhancing community resilience in the face of submarine landslide tsunamis in Mayotte.

5 Limitations

While the Landslide-Tsurrogate v1.0 framework provides an efficient approach for propagating uncertainties in landslide-generated tsunami forecasts, several limitations should be acknowledged. These limitations arise from the assumptions made in the representation of landslide sources, the deterministic numerical simulations used to construct the surrogate models, and the surrogate modeling approach itself.

5.1 Surrogate modeling approach

The Landslide-Tsurrogate v1.0 framework relies on gPCE to approximate the response of the deterministic landslide–tsunami simulations across the stochastic parameter space. While gPCE provides an efficient and mathematically rigorous surrogate modeling approach, it also introduces specific limitations that must be considered.

First, the construction of a polynomial chaos surrogate requires that the quantity of interest be defined for all deterministic simulations corresponding to the quadrature nodes. In practice, this means that all model evaluations must return finite values. However, in tsunami simulations certain output variables – particularly those related to coastal inundation – may exhibit undefined or discontinuous behavior. For example, at locations where the tsunami does not reach the coastline for some parameter combinations, quantities such as maximum elevation, time of arrival or maximum speed may be zero or undefined, while they become positive for other combinations. These wet–dry transitions introduce discontinuities in the response surface that are difficult to represent accurately with polynomial expansions.

Second, polynomial surrogate models generally perform best when the system response varies smoothly with respect to the stochastic parameters. In the context of landslide-generated tsunamis, strongly nonlinear processes such as shoreline interaction, or abrupt changes in landslide mobility can introduce localized nonlinearities or threshold effects. These behaviors may reduce the convergence rate of the polynomial expansion and lead to reduced surrogate accuracy unless sufficiently high polynomial orders are used.

Finally, inundation-related quantities of interest, such as maximum run-up or flooded area, are particularly challenging to represent within the gPCE framework because they depend on complex wetting–drying dynamics and may exhibit non-smooth spatial patterns. For this reason, the present implementation primarily focuses on tsunami wave characteristics at locations outside of the inundation zone, where the response is generally smoother and more amenable to polynomial approximation. Extending the framework to accurately represent inundation processes may require alternative surrogate modeling strategies or hybrid approaches combining gPCE with other statistical or machine learning methods.

5.2 Deterministic modeling strategy

The accuracy of the surrogate models ultimately depends on the ability of the underlying deterministic numerical model to represent the physical processes governing landslide-tsunami generation, propagation, and coastal interaction. The surrogate model does not introduce new physical information but rather approximates the response of the deterministic model across the stochastic parameter space. Consequently, any simplifications or inaccuracies in the deterministic model are directly inherited by the surrogate model.

In particular, the representation of submarine landslide dynamics may vary significantly depending on the assumptions made in the models (e.g., rigid-block motion, depth-averaged granular flow models, or fully coupled hydro-mechanical approaches). The suitability of a given model depends on the characteristics of the considered site, including bathymetric complexity, sediment properties, and the expected failure mechanisms. For some environments, simplified landslide parameterizations may provide reasonable approximations of tsunami generation, whereas other cases may require more advanced three-dimensional or multilayer models to capture complex flow behavior and landslide–water interactions. The present framework is therefore not tied to a specific deterministic model but relies on the user to select an appropriate model for the site under investigation.

Furthermore, in the Mayotte test case, submarine landslides are represented using simplified geometric parameterizations (e.g., paraboloid shapes for the initial mass released) and a limited number of stochastic variables describing their location, volume, and rheological properties. While this representation allows the dimensionality of the stochastic space to remain manageable, it necessarily simplifies the complex geometry and mechanical behavior of real submarine landslides.

Natural submarine slope failures may involve irregular geometries, heterogeneous sediment properties, progressive failure mechanisms, or multi-stage collapses. Such processes may not be fully captured by simplified parametric representations. As a result, the surrogate models should be interpreted as providing probabilistic estimates of tsunami forecasts for a range of idealized landslide scenarios rather than precise predictions of specific failure events.

5.3 Dimensionality of the stochastic parameter space

The efficiency of the surrogate modeling approach depends strongly on the number of stochastic variables used to describe the system. As the number of uncertain parameters increases, the number of deterministic simulations required to construct accurate polynomial surrogate models grows rapidly. Although sparse quadrature methods substantially reduce the number of required simulations compared with full tensor grids, the approach may still become computationally demanding when the dimensionality of the stochastic space is large.

For this reason, the framework is best suited to applications where the dominant sources of uncertainty can be represented by a relatively small number of stochastic variables (up to 6). In practice, careful parameter selection and sensitivity analysis may be necessary to identify the most influential parameters and avoid unnecessary expansion of the stochastic space.

5.4 Representation of uncertainty

The Landslide-Tsurrogate framework represents uncertainty through prior probability distributions assigned to the stochastic variables describing landslide characteristics and environmental conditions. The definition of stochastic variables and their probability distributions are inherently site dependent. Bathymetric conditions, sediment characteristics, tectonic setting, and coastal morphology all influence landslide dynamics and tsunami propagation.

However, in many practical applications, the available geological and geotechnical information may be limited, and the choice of probability distributions may rely partly on expert judgment or regional analogues. Due to these limitations, the prior distributions of all stochastic variables are assumed uniform in the current version of the Landslide-Tsurrogate framework. As a result, the probabilistic outputs of the framework reflect both natural variability and epistemic uncertainty associated with limited knowledge of the system.

Future applications could benefit from improved constraints on landslide occurrence, sediment properties, and triggering mechanisms, which would allow more realistic parameter distributions to be defined. Consequently, future developments of the Landslide-Tsurrogate framework should include the possibility to build surrogate models with non-uniform prior distributions of the stochastic variables.

In summary, the Landslide-Tsurrogate framework should be viewed as a flexible methodological approach that can be adapted to different environments. Applying the framework to other sites requires careful consideration of the local geological setting, the selection of appropriate deterministic models, and the definition of stochastic variables that adequately represent the dominant uncertainties affecting landslide-generated tsunamis.

6 Conclusions

This study addresses the critical challenge of efficiently assessing tsunami hazards generated by submarine landslides, a phenomenon characterized by high physical complexity and considerable uncertainty. Traditional Probabilistic Tsunami Forecast (PTF) approaches are computationally expensive, limiting the ability to explore the full range of uncertain landslide parameters. To overcome this limitation, Landslide-Tsurrogate v1.0, a surrogate modeling framework designed to approximate tsunami maximum elevation, maximum speed and time of arrival as functions of submarine landslide characteristics, has been developed to provide fast and reliable predictions suitable for probabilistic analyses.

The key takeaways of the Landslide-Tsurrogate v1.0 model are:

1. Surrogate modeling efficiency for probabilistic assessment of landslide-generated tsunamis

   The Landslide-Tsurrogate framework combines high-fidelity numerical simulations with generalized Polynomial Chaos Expansion (gPCE) and uniform distributions of the submarine landslide parameters to construct surrogate models. The methodology incorporates a rigorous training and validation workflow, ensuring that the surrogate accurately captures the non-linear response of the tsunami system while drastically reducing computational costs compared to full deterministic numerical simulations. However, adopting this approach implies that many of the randomly generated submarine landslide scenarios may, in fact, lack the capacity to produce significant tsunamis in the area under investigation. This provides strong motivation to better constrain uncertainty distributions, that is, to more accurately define the prior distributions of the stochastic parameters using additional data. Such refinement can be achieved through Bayesian inference approaches (e.g., von Toussaint, 2011), or more generally via model calibration and data assimilation techniques. These aspects fall outside the scope of the present study and are left for future work.

2. Sensitivity Analysis

   Within the gPCE framework, performing sensitivity analysis is trivial and demonstrates that key landslide parameters – such as volume and geographical location – dominate tsunami maximum elevation, maximum speed and time of arrival, confirming both the physical relevance of the model and the ability of gPCE to reproduce complex dependencies.

3. Real-world application

   Application to the Mayotte submarine landslide case illustrates the practical utility of Landslide-Tsurrogate v1.0 for real-world PTF studies. Using landslide scenarios informed by geological and geophysical data, the surrogate models along Petite Terre efficiently generated probabilistic maximum elevation, maximum speed and time of arrival distributions and identified regions of elevated tsunami risk. Results highlight both the expected range of tsunami characteristics and the relative influence of individual landslide parameters on hazard outcomes, providing actionable insight for risk management and early warning considerations.

4. Scalability

   The Landslide-Tsurrogate v1.0 approach can be extended to multi-slide scenarios, higher-resolution bathymetry, and real-time hazard forecasting for enhanced coastal risk management.

Overall, Landslide-Tsurrogate v1.0 demonstrates that surrogate modeling is a powerful tool for bridging the gap between high-fidelity tsunami simulations and practical PTF applications. The framework allows extensive uncertainty quantification, sensitivity analysis, and scenario exploration while maintaining computational efficiency. Future work will focus on addressing missing values generated by the deterministic landslide–tsunami models (e.g., time of arrival set to −1 in HySEA) to enable the extension of the surrogate model implementation to the inundation zone and the production of accurate inundation maps essential for risk assessment. In addition, Landslide-Tsurrogate v1.0 will be expanded to incorporate alternative prior distributions of the stochastic variables beyond the uniform case, thereby enhancing its general applicability. By enabling rapid yet reliable tsunami hazard assessment, Landslide-Tsurrogate v1.0 represents a valuable step toward improved coastal risk management in landslide-prone regions.

Appendix A: Landslide-Tsurrogate v1.0: Technical Description

Landslide-Tsurrogate v1.0 is implemented in both Python and MATLAB to ensure accessibility and flexibility for a broader range of users within the scientific community. The same simple structure (Fig. A1) is used for both the Python and MATLAB codes. To run the model, three folders are required: code, data, and results.

The code folder contains the seven user-oriented programs that each produce a single output file. They should be executed sequentially following the workflow presented in Fig. A2. While each program can run independently, it depends on input files generated by the preceding step. This modular structure is intentional, as certain steps – such as running deterministic simulations – require user intervention and can be time-consuming. Separating the programs in this way allows for greater flexibility and practicality in managing the workflow.

The data folder includes two MATLAB-format input files that the user must generate when running Landslide-Tsurrogate v1.0. The first, input_simus.mat, stores the outputs of the N deterministic simulations corresponding to the highest total order P_t+1. It contains the maximum tsunami elevation zeta_max_surf[N,nx,ny], the maximum tsunami speed velo_max_surf[N,nx,ny], and the tsunami time of arrival time_max_surf[N,nx,ny], all defined over the simulation domain with spatial dimensions [nx,ny]. The second file, surrogate_model_locations.mat, stores the indices index_coast[M] – following MATLAB conventions in the example included in the distribution of the code – corresponding to the M locations where the surrogate models are to be constructed.

Finally, the MATLAB-format output files of the seven main programs are created in the results folder.

Figure A1 Structure of the Landslide-Tsurrogate v1.0 software.

The sequential execution of the seven user-oriented programs and their associated outputs (illustrated in Fig. A2) corresponds to the first six core building blocks of the workflow while the implementation of the user-friendly interfaces is done outside of the main workflow and illustrated, for the Mayotte test case, in Fig. A3.

Figure A2 Workflow pipeline of the main programs and associated functions in Landslide-Tsurrogate v1.0 for the generation of the surrogate models. User-oriented programs are shown in Blue boxes, core building-block functions in Orange, and supporting functions in Green and Yellow. Functions or programs that require user editing are highlighted in Grey. The input/output flow between components is represented by arrows, labeled with the names of the variables saved in the corresponding files.

Figure A3 Example of Landslide Probabilistic Tsunami Forecast (PTF) user-friendly interfaces for the Mayotte test case: (a) developed in MATLAB and distributed as program files for experts or as a Windows installer (.exe) for non-experts and (b) developed in Python and distributed as Jupyter Notebook for experts or deployed as a web application through Voilà for non-experts. Sources: Esri, i-cubed, USDA, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the Gis User Community | Powered by Esri.

1. Landslide_Tsurrogate_step_1 _user_input

   The user specifies the number r of stochastic variables, the definition intervals $[a_i,b_i]i\in \mathrm{1},\mathrm{2},\mathrm{\dots }r$ for each stochastic variable, the maximum total order P_t+1, and the quadrature rule (0 for GP, 1 for DGP). This program generates the file output_users.mat.

2. Landslide_Tsurrogate_step_2 _input_parameters

   This program creates the combinations of parameters that sample the definition intervals of the stochastic variables, according to the quadrature rule selected by the user and following the theoretical framework described in Sect. 2. It calls the main function surrogate_model_gauss_patterson_param, which constructs the final nested grids using the Heiss and Winschel algorithm (Heiss and Winschel, 2008).

   In practice, the multi-index set α is generated by the function get_seq, while the multivariate abscissa (Eq. 8) are computed by the function gauss_patterson_nested_rule_nd, based on the univariate GP abscissa (Eq. 10) provided by gauss_patterson_nested_rule_1d.

   The final parameter set is given by the unique combinations of these multivariate abscissa. These are saved in the file output_param.mat and used as inputs for the deterministic landslide-generated tsunami simulations.

3. Landslide_Tsurrogate_step_4 _format_simulations

   After the deterministic simulations are completed, the user must provide the tsunami hazard variables – maximum tsunami elevation, maximum tsunami speed, and time of arrival – in the specific format required by Landslide-Tsurrogate v1.0. This program serves as a template for formatting those outputs and should be edited and adapted to match the structure and naming conventions of the user's simulation results. This step is essential to ensure compatibility with the subsequent surrogate modeling process and the results are saved in the file output_model.mat.

4. Landslide_Tsurrogate_step_5 _coefficients

   This program requires the inputs generated in the three previous steps to compute the deterministic coefficients of the gPCE expansion (Eqs. 6 and 7). It relies on the main function surrogate_model_gauss_patterson_coeff, which builds upon the same functions used in Step 2, with the addition of norm2_legendre_polynomials to compute the norms of the multivariate Legendre polynomials. This program produces the main output of Landslide-Tsurrogate v1.0 which is saved in the file output_coeff.mat.

5. Landslide_Tsurrogate_step_6 _evaluation

   This program also requires the inputs generated in all the previous steps. It is based on the main functions surrogate_model_gauss_patterson_evals and surrogate_model_gauss_patterson _sensi, and produces two essential output files used for the validation of the surrogate models: output_evals.mat and output_sensi.mat. Subsequently, the user can run the program Landslide_Tsurrogate_step_6_bis _visualization to explore the results and make informed decisions regarding the truncation of the gPCE expansion or the elimination of certain stochastic variables.

6. Landslide_Tsurrogate_step_7_PTF

   In this final step, the program uses inputs from Step 1 and Step 4 to compute the Probabilistic Tsunami Forecast (PTF) at the locations where surrogate models have been constructed. It relies on the function surrogate_model_gauss_patterson_PTF, which evaluates the surrogate models based on a user-defined input distribution of the stochastic variables. The resulting PTF data are stored in the file output_PTF.mat.

7. User-friendly interface

   The generation of user-friendly interfaces with MATLAB and Python is illustrated for the Mayotte test case (see Fig. A3). In MATLAB, the interface can be converted into a standalone executable, which eliminates the need for internet access or complex software installations. This ensures that users in remote or resource-limited environments can run simulations and visualize results immediately, making the tool practical for real-world decision-making scenarios. In Python, the user-friendly interface leverages Jupyter Notebooks and Voilà (https://voila.readthedocs.io/en/stable/index.html; last access: 31 October 2025) to offer the possibility to provide interactive web applications. Jupyter Notebooks are widely used in scientific computing to combine code, documentation, and visualization in a single environment. With Voilà, these notebooks can be converted into fully interactive web applications, allowing users to adjust input parameters, run simulations, and visualize results dynamically – all through a standard web browser. Importantly, both the MATLAB and Python approaches require no programming knowledge from the end user: the complexity of the model is fully hidden behind intuitive controls, while outputs such as probability density functions, cumulative distributions, and hazard maps are immediately accessible. These design choices illustrate the philosophy of Landslide-Tsurrogate: creating operational products from advanced computational research, ensuring that PTF can be used by emergency planners, local authorities, and other decision-makers.

To summarize, Landslide-Tsurrogate v1.0 provides a modular and transparent workflow for constructing and evaluating surrogate models tailored to landslide-generated tsunami hazards. Each user-oriented program is designed to function independently while maintaining a logical sequence that ensures reproducibility and flexibility. From the initial specification of stochastic inputs to the final generation of PTF outputs, the tool incorporates rigorous mathematical foundations and practical considerations to support a wide range of scientific applications. By offering both Python and MATLAB implementations and supporting user customization, Landslide-Tsurrogate v1.0 aims to be accessible to diverse users while maintaining high standards of technical robustness. This flexibility is exemplified in the next section through its application to the Mayotte test case, which demonstrates the end-to-end capabilities of the model in a real-world context.

Appendix B: Multilayer HySEA Model

The deterministic simulations for the Mayotte case study were performed using the Multilayer HySEA (Hyperbolic Systems and Efficient Algorithms) model, which belongs to a broader class of long-wave numerical models designed to simulate tsunamis generated by submarine landslides (Macías et al., 2021a, b). Earlier one-layer configurations of HySEA have been extensively validated against laboratory benchmarks (e.g., Poulain et al., 2023) as well as real-world tsunami events at field scale (e.g., Linares-Sánchez et al., 2011; Gonzalez-Vida et al., 2019). The HySEA modeling framework is formulated within the shallow-water approximation and explicitly represents both the motion of the granular landslide and the hydrodynamics of the water column, including their dynamic interaction through a frictional coupling at the landslide–water interface (Macías et al., 2021a, b). In the one-layer configuration, the water motion is described using depth-averaged non-hydrostatic (dispersive) equations, which can be written as follows:

$$\begin{array}{}\text{(B1)}& \left\{\begin{array}{l}{\partial }_{t}h+{\partial }_x\left(hu\right)+{\partial }_y\left(hv\right)=\mathrm{0}\\ {\partial }_t\left(hu\right)+{\partial }_x(h{u}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}+hp)+{\partial }_y\left(huv\right)\\ =(gh+p){\partial }_x(H-h_s)-S_x\\ {\partial }_t\left(hv\right)+{\partial }_x\left(huv\right)+{\partial }_y(h{v}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}+hp)\\ =(gh+p){\partial }_y(H-h_s)-S_y\\ {\partial }_t\left(hw\right)+{\partial }_x\left(hwu\right)+{\partial }_y\left(hwv\right)=p\\ \mathrm{\nabla }\cdot \mathit{u}+\frac{w_{\mathit{\eta }}-w_{\mathrm{0}}}{h}=\mathrm{0}.\end{array}\right.\end{array}$$

The model variables include the water surface height h, the horizontal water velocity vector u=(u  v), the vertical water velocity w, the non-hydrostatic pressure p. The gravity acceleration is denoted by g. The last equation comes from the incompressibility condition and involves the vertical velocity at the free surface w_η and at the bottom w₀. At the free surface, the non-hydrostatic pressure is prescribed to vanish. The term S=(S_x  S_y) in the momentum equations include all the dissipative terms (see Marbœuf et al., 2025 for a detailed formulation of the model).

Accounting for more vertical layers than the one-layer model described by Eq. (B1) improves the representation of the velocity field (Fernández-Nieto et al., 2016; Garres-Díaz et al., 2020; Escalante et al., 2023; Marbœuf et al., 2025) and makes it possible to better capture dispersive effects (Macías et al., 2021a). We use here 4 layers as the best compromise between accuracy and computational cost to simulate landslide-induced tsunamis in Mayotte (Marbœuf et al., 2025).

Note that the model does not differentiate between offshore and inland domains, and the propagation of flow over land is handled without introducing additional models, apart from the standard wet–dry treatment applied at shallow water depth (Macías et al., 2021a; Castro et al., 2005). No other model – such as a breaking wave or a vertical viscosity term (Roeber et al., 2010; Scala et al., 2024) – is applied.

A Savage-Hutter model (Savage and Hutter, 1989) describes the landslide motion and complements the system (Eq. B1) as follows:

$$\begin{array}{}\text{(B2)}& \left\{\begin{array}{l}{\partial }_t{h}_s+{\partial }_x\left(h_s{u}_s\right)+{\partial }_y\left(h_s{v}_s\right)=\mathrm{0}\\ {\partial }_t\left(h_s{u}_s\right)\\ +{\partial }_x(h_s{u}_s^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{1}-r\left)g{h}_s^{\mathrm{2}}\right)+{\partial }_y\left(h_s{u}_s{v}_s\right)\\ =(\mathrm{1}-r)g{h}_s{\partial }_{x}H-S_{s,x}\\ {\partial }_t\left(h_s{v}_s\right)\\ +{\partial }_x\left(h_s{u}_s{v}_s\right)+{\partial }_y(h_s{v}_s^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{1}-r\left)g{h}_s^{\mathrm{2}}\right)\\ =(\mathrm{1}-r)g{h}_s{\partial }_{y}H-S_{s,y},\end{array}\right.\end{array}$$

where the unknowns are the landslide thickness h_s and the horizontal velocity vector u_s=(u_s  v_s). The density ratio r is defined as $r={\mathit{\rho }}_f/{\mathit{\rho }}_s$ where ρ_f is the fluid density and ρ_s the grain phase density. The system (Eq. B2) differs from Eq. (B1) through two main features: (i) the reduced gravity term (1−r) g taking into account buoyancy effects when the landslide is underwater, and (ii) the absence of dispersive effects (non-hydrostatic pressure and vertical velocity). The terms ${\mathbf{S}}_s=\left(S_{s,x}{S}_{s,y}\right)$ are the friction forces where the friction coefficient μ(I) appears (Pouliquen and Forterre, 2002; Jop et al., 2006; Mangeney et al., 2007). The coefficient μ(I) is given by

$$\begin{array}{}\text{(B3)}& \mathit{\mu }\left(I\right)=\text{tan}\left(\mathit{\delta }\right)+{\displaystyle \frac{\text{tan}\left({\mathit{\delta }}_{\mathrm{2}}\right)-\text{tan}\left(\mathit{\delta }\right)}{\mathrm{1}+\frac{I_{\mathrm{0}}}{I}}}.\end{array}$$

The inertial number I is chosen here as in Martin et al. (2023)

$$\begin{array}{}\text{(B4)}& I={\displaystyle \frac{\frac{\mathrm{5}}{\mathrm{2}}\parallel {\mathit{u}}_s\parallel d}{h_s\sqrt{g(\mathrm{1}-r)h_s\mathrm{\Phi }}}}\end{array}$$

where the empirical value of the constant I₀=0.279 comes from Jop et al. (2006); Martin et al. (2023). The friction angles δ and δ₂, the grain diameter d and the solid volume fraction Φ are specified in Sect. 4.2.1.

A complete description of the multilayer HySEA model and its implementation can be found in Castro et al. (2011), Escalante et al. (2018), Macías et al. (2021a, b) and Marbœuf et al. (2025).