Performance and limits of a shallow-water model for landslide-generated tsunamis: from laboratory experiments to simulations of flank collapses at Montagne Pelée (Martinique)

Geophys. J. Int. (2023) 233, 796–825 https://doi.org/10.1093/gji/ggac482 Advance Access publication 2022 December 06 GJI Heat Flow and Volcanology Performance and limits of a shallow-water model for landslide-generated tsunamis: from laboratory experiments to simulations of flank collapses at Montagne Pel´ ee (Martinique) P. Poulain ,1A. Le Friant,1A. Mangeney,1S. Viroulet,1,2E. Fernandez-Nieto,3M. Castro Diaz,4M. Peruzzetto ,1,5G. Grandjean,5F. Bouchut,6R. Pedreros5and J-C. Komorowski 1 1Universit´ eParisCit ´ e, Institut de Physique du Globe de Paris, departements of volcanic systems and sismology, CNRS, F-75005 Paris, France. E-mail: poulainpab[email protected] 2Institut de M´ ecanique des Fluides de Toulouse (IMFT) - Universit´ e de Toulouse, CNRS, 31400 Toulouse, France 3Dpto. Matematica Aplicada I, ETS Arquitectura - Universidad de Sevilla, Avda. Reina Mercedes S/N, 41012 Sevilla, Spain 4Departamento de Analisis Matematico, Estad´ ıstica e Investigaci´ on Operativa y Matem´ atica Aplicada. Universidad de Malaga, Facultad de Ciencias, Campus Teatinos S/N, 29081 M´ alaga, Spain. 5BRGM (French Geological Survey), Orl´ eans, France 6Laboratoire d’Analyse et de Math´ ematiques Appliqu´ ees (UMR 8050), CNRS, Univ. Gustave Eiffel, UPEC, F-77454, Marne-la-Vall´ ee, France Accepted 2022 November 28. Received 2022 November 28; in original form 2022 February 8 SUMMARY We investigate the dynamics and deposits of granular flows and the amplitude of landslidegenerated water waves using the HySEA depth-averaged shallow-water numerical model, both at laboratory and field scales. We evaluate the different sources of error by quantitatively comparing the simulations with (i) new laboratory experiments of granular collapses in different conditions (dry, immersed, dry flow entering water) and slope angles and (ii) numerical simulations made with the SHALTOP code that describes topography effects better than most depth-averaged landslide-tsunami models. For laboratory configurations, representing the limits of the shallow-water approximation in such models, we show that topography and non-hydrostatic effects are crucial. When topography effects are accounted for empirically—by artificially increasing the friction coefficient and performing non-hydrostatic simulations—the model is able to reproduce the granular mass deposit and the waves recorded at gauges located at a distance of more than two to three times the characteristic dimension of the slide with an error ranging from 1 to 25 per cent depending on the scenario, without any further calibration. Taking into account this error estimate, we simulate landslides that occurred on Montagne Pel´ ee volcano, Martinique, Lesser Antilles as well as the generated waves. Multiple collapse simulations support the assumption that large flank collapses on Montagne Pel´ ee likely occurred in several successive subevents. This result has a strong impact on the amplitude of the generated waves and thus on the associated hazards. In the context of the ongoing seismic volcanic unrest at Montagne Pel´ ee volcano, we calculate the debris avalanche and associated tsunamis for two potential flank-collapse scenarios. Key words: Numerical modelling; Tsunamis; Submarine landslides; Volcanic arc processes; Volcanic hazards and risks. 1 INTRODUCTION Volcano flank collapses are recurrent processes that are part of the evolution of a volcanic edifice (Siebert 1984;McGuire1996). These gravity induced mass movements and debris avalanches can trigger tsunamis (e.g. Kelfoun et al.2010; Abadie et al.2012, Gylfadottir et al. 2017;Barrettet al.2019; Williams et al.2019; Perttu et al.2020;Pariset al.2020). The dynamics of these natural granular flows and thus the generation of the tsunami is still an open issue (e.g. Yavari-Ramshe & Ataie Ashtiani 2016). Numerical simulations provide a unique tool to quantify these processes (see Yavari-Ramshe et al.2015;Behrenset al.2021, for a review). However, significant uncertainties have to be taken into account in the simulation of landslide-generated tsunamis regarding (i) the position and characteristics of the released mass, (ii) the estimation of sources probabilities (Lane et al.2016; Løvholt et al.2020), (iii) 796 C The Author(s) 2022. Published by Oxford University Press on behalf of The Royal Astronomical Society. Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 Model for landslide-generated tsunamis 797 the rheological laws and parameters to be used in the simulations and (iv) the inevitable approximations made in the physically based model. Each of these uncertainties is a source of active research. Uncertainties regarding the initial conditions, the material composition (de Haas et al.2015), the rheological laws and associated parameters (Delannay et al.2017; Peruzzetto et al.2020) can be addressed by running multiple simulations (Løvholt et al.2020), requiring huge computational efforts. Simplified depth-averaged models are thus often used for field-scale studies even if they are less accurate than full 3-D models (Gittings 1992; Abadie et al.2008; Horrillo et al.2013;Maet al.2013). Simplified models are based on the so-called thin-layer approximation for the landslide and shallow-water or long-wave approximation for the tsunami: the landslide thickness, respectively water depth, is assumed to be small in comparison with the landslide downslope extension, respectively, tsunami wavelength. Since the pioneering work of Savage & Hutter (1989), depth-averaged landslide-tsunami models have seen significant developments (Gittings 1992; Mangeney et al.2000; Fernandez-Nieto et al.2008; Abadie et al.2010; Kelfoun 2011;Kelfoun&Vargas2015). There are still however a number of mathematical, mechanical and numerical shortcomings in these models that, when not appropriately addressed, may induce significant error in the simulation. In this work, we focus on part of these modelling challenges and on quantifying the effect of their accurate implementation. We know there is presently no model for landslide-generated tsunamis that accurately deals with all the aspects described below at a reasonable computational cost and no precise idea of the errors made when neglecting some of them. Despite strong differences between models, most studies conclude that their model well reproduces laboratory experiments or field data. We believe that a key point for future model improvement is to show, quantify and analyse uncertainties to reduce bias in modelling. A major modelling challenge, that we will further call topography effects, concerns the accurate derivation of thin-layer depth-averaged equations for granular flows for a general 2-D topography b(x,y) (Peruzzetto et al.2021). This derivation induces two related difficulties: (1) complex terms linked to topography curvature have to be taken into account and (2) the thin-layer approximation and depth-averaging should be done in a reference frame following the topography. The first rigorous derivations of depth-averaged flow equations over arbitrary topographies, with a curvature small enough to be compatible with the thin-layer approximation, were provided by Bouchut et al.(2003) and Bouchut & Westdickenberg (2004)fora1-Dand 2-D topography, respectively, and solved numerically by Mangeney-Castelneau et al.(2005) and Mangeney et al.(2007a) in the SHALTOP code. In the resulting equations, two curvature terms appear: one ensures that the flow velocity remains tangent to the topography (as assumed for thin-layer models) and the other appears when frictional rheologies are used. To our knowledge, only a few landslide models accurately account for these effects. Only the curvature term related to friction is included in the RAMMS code (Christen et al.2010; Fischer et al.2012) while in FOAM, Rauter & Tukovic (2018) and Rauter et al.(2018) take curvature terms into account implicitly by solving numerically for the bottom pressure. As there is no model for landslide-generated tsunamis accounting for curvature effects, we will use SHALTOP here to (i) estimate the error in landslide dynamics and deposits estimated by models simulating both the landslide and the tsunami without properly accounting for topography effects and (ii) calibrate the models to artificially reduce this error. A related challenge, poorly addressed in the literature, is the different coordinate systems that have to be considered to describe the landslide and the tsunami. Indeed, to accurately derive the thin-layer or shallow-water approximation, it is necessary to measure the thickness of the material involved in the direction roughly normal to the surface of the considered thin flow. Thus, the shallow approximation should be imposed in a Cartesian horizontal/vertical frame for the tsunami part and in a coordinate system adapted to the topography variations for the landslide part as stated above. Very few models consider two different coordinates systems (Ma et al.2015;Delgado-Sanchez et al. 2020) and none account for topography curvature. While the Volcflow model (Kelfoun et al.2010) uses a topography-related reference frame for both the tsunami and landslide parts, the D-Claw model (George et al.2017) uses Cartesian horizontal/vertical coordinates for both parts, as in the HySEA code used here. The introduction of non-hydrostatic terms in the equations [Garres-Diaz et al.2020; eq. (4.9) in Bouchut et al. 2016a] helps reduce the error related to inappropriate coordinate systems (Delinger & Iverson 2004; Zhang et al.2021a; Garres-Diaz et al. 2020). But additional considerations must also be taken into account in the modelling. The third challenge is to describe grain-fluid interactions and dilatancy in landslides. Depth-averaged thin-layer landslide models have recently been improved to partly take into account these interactions (Pitman & Lee 2005, Pudasaini 2012; Mergili et al.2017), in particular dilatancy and associated pore fluid pressure effects (Iverson & George 2014; Bouchut et al.2016b,2017; Pudasaini & Mergili 2019). However, these models contain several parameters that are difficult to estimate for field applications and, in their present state, do not accurately describe the curvature of complex topographies. In the simulation of tsunamis generated by landslides, the moving mass is commonly considered as a simplified grain–fluid mixture, where fluid effects are reduced to buoyancy and drag between the granular layer and the water layer, without dilatancy effects (Fernandez-Nieto et al.2008; Yavari-Ramshe et al.2015;Maet al.2015, Zhang et al.2021b). The main advantage of these models is to involve only a few parameters. This is in particular the case for the HySEA model (Macias et al.2021b) that will be used in this work. The fourth challenge concerns the coupling strategy between the landslide flow and tsunami generation. A first approach considers the combination of the landslide and water as a single layer containing a mixture (the grain-water mass is considered as an effective media; George et al.2017) or a multiphase fluid (the grain and fluid equations are both solved; Horrillo et al.2013; Grilli et al.2014; Tappin et al.2014). A second strategy is to couple models dealing with the different processes: granular flow along the slope, tsunami generation and tsunami propagation (Abadie et al.2012;Kimet al.2019). A third approach, used in our study, considers the granular flow and the water as two separate layers that are coupled through friction between the two layers as in HySEA (Fernandez-Nieto et al.2008; Macias et al.2021b). Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 798 P. Poulain Finally, different level of approximation can be made for the landslide and tsunami description. Lots of models solve shallow-water (or thin-layer) depth-averaged equations, assuming hydrostatic pressure for both the granular and water layers (Jiang & LeBlond 1992; Heinrich et al.1998; Mangeney et al.2000; Majd & Sanders 2014; Fine et al.2003,2005; Fernandez-Nieto et al.2008; Yavari-Ramshe & Ataie-Ashtiani 2015). These models suffer from the lack of dispersion in tsunami propagation, potentially leading to significant bias. Indeed, wave dispersion can be significant (Glimsdal et al.2013;Maet al.2015) when the water wavelength is smaller than or about the same order of magnitude as the water depth (Yavari-Ramshe & Ataie-Ashtiani 2016). For more than two decades, more advanced depth-averaged models to describe water wave propagation have been developed based on Boussinesq-type equations (non-hydrostatic pressure), that are weakly or fully dispersive (Kirby et al.1998; Zhou et al.2011;Shiet al.2012; Popinet 2015;Kirbyet al.2022), as is the case for the HySEA model (Macias et al.2021b). Such depth-averaged non-hydrostatic models are also valid only for a sufficiently small ratio of the water depth to the tsunami wavelength, a limitation inadequately discussed in the literature (Macias et al.2021a). The error when using HySEA and the performance and shortcomings of these types of models will be estimated here by comparing simulations with (i) results from the landslide model SHALTOP better describing topography effects but not the water wave and (ii) eight laboratory experiments for dry, immersed and dry entering water initial configurations. These experiments are conducted in what can be considered to be the worst situations (strong slope variations, non-shallowness of part of the flow, conditions beyond the limits of the Boussinesq one-layer model, etc.), in order to quantify maximum errors. After estimating the error, we will apply the code to past volcanic flank collapses in Martinique for which we have geological data. Once calibrated on these events, we will use this model to preliminary estimate the impact of potential flank collapses related to the increasing volcanic activity of Montagne Pel´ ee, Martinique, France (OVSM-IPGP 2019-2021). We will analyse these results in light of our study of the model limitations, in order to better understand the influence of the volume and morphology of the released mass, multiple-stage events, and the presence of a simplified erodible bed on the run-out distance and the emplacement of the deposits as well as on the generated water waves and their propagation. In the next section, we first describe the HySEA and SHALTOP models. In Section 3, we present and compare the laboratory experiments and associated simulations of granular flows using the two models. In Section 4, we investigate Montagne Pel´ ee flank collapses simulated with HySEA and discuss the results in Section 5. 2 NUMERICAL MODELS To simulate tsunami generated by landslides, we use the HySEA code describing a grain–fluid layer flowing over a fixed bottom z(x,y), with a water layer on top of it (Macias et al.2021b). In HySEA, the shallow-water and thin-layer approximations as well as the flow averaging are done in a Cartesian horizontal/vertical reference frame for both the granular and the water layers. As described above, this assumption induces error that increases with the slope angle when it is applied to describe granular flow (Delgado-Sanchez et al. 2020). To quantify such error related to topography effects, we use the SHALTOP code describing dry granular flows over complex topographies, extended to account for buoyancy effects (see Section 2.3 for the granular rheology). In this way, we can quantify how much the friction parameters should be artificially increased in HySEA to roughly reproduce the deposit simulated when topography effects are appropriately taken into account. Due to depth-averaging, these models assume that the flow is either completely mobile or static throughout its depth, which is not always the case, especially during the starting and the stopping phases (e.g. Ionescu et al.2015; Martin et al.2017)orifmaterialisentrainedor eroded (Iverson & Ouyang 2015). To approximate erosion processes in Section 5.2 for field application, we followed Moretti et al.(2012) by assuming that the static mass on the slope is just added to the flowing mass as the latter passes over the erodible bed. The velocity is then averaged over the whole mass depth (initially static mass and flowing mass). This oversimplified approach made it possible to quite well reproduce the dynamics of natural erosive flows when compared to seismic data. It is obviously less accurate than models describing both the static and flowing part in the granular avalanche (Fernandez-Nieto et al.2016). In HySEA and SHALTOP as in most landslide models, the water and the granular mass are assumed to be incompressible and homogeneous. Therefore, natural complexities such as density variations due to the expansion or contraction of the material and their impact on pore fluid pressure, segregation and fragmentation processes, or the incorporation of air and/or water are not taken into account, potentially leading to significant bias in the simulation (see Delannay et al.2017, for a review and Rauter et al.2022). 2.1 HYSEA numerical model The two-layer HySEA code is a 2-D extension of the model proposed by Fernandez-Nieto et al.(2008), but using Cartesian coordinates and accounting for dispersion effects (i.e. non-hydrostatic) during water wave propagation. HySEA has been developed by the EDANYA group (Castro et al.2005; Macias et al.2015) and is based on an efficient hybrid finite-volume-finite-difference numerical scheme on GPU architectures (Macias et al.2021a). From the depth-averaged equations, the models solves the unknowns (h1,u1)and(h2,u2), representing the vertical height and horizontal velocity of the water and granular layer, respectively, averaged in the vertical z-direction. The velocities u1and u2are thus 2-D vectors in the horizontal (x,y) plane. Two additional unknowns are involved: the non-hydrostatic pressure pand the vertical velocity w1, leading to Green–Naghdi equations (Green & Naghdi 1976) written as a non-hydrostatic pressure system. The mass and momentum balance equations Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 Model for landslide-generated tsunamis 799 for the water and granular layers are: ∂th1+∇·(h1u1)=0, ∂tu1+u1·∇u1+∇(g(b+h1+h2)+p)+p h1 ∇h1+3 2(b+h2)=Sf,1, ∂tw1+u1·∇w1=3 2 p h1 , ∇·u1+w1−u1·∇(b+h2) h1/2=0, ∂th2+∇·(h2u2)=0, ∂tu2+u2·∇u2+∇ (g(b+h2)+g(1−)h1)=Sf,2(1) with Sf,1=mfh2 h2+rh1u2−u1(u2−u1)if h2>0, −gn2 h4/3 1 u1u1if h2=0.(2) Sf,2=−(1−)mf h1 h2+rh1 u2−u1(u2−u1)−gμu2 u2,(3) where b(x,y) defines the topography and gravity acceleration gis reduced by a factor: κ=1−ρ1/ρ2(4) with ρ1the density of water, ρ2the density of the granular material and r=ρ1/ρ2. The first line of Sf,1and Sf,2represent an empirical friction term between the water and granular layer involving the square of their relative velocity and a drag coefficient mf. When the water is directly in contact with the bottom—second line of eq. (2), a similar friction between the water and the bottom is assumed involving the Manning coefficient n. These friction terms are used in HySEA and other landside-tsunami models and the detailed study of these friction terms is beyond the scope of this paper. The 3 first equations of system (1) represent the mass, horizontal and vertical momentum balance equations for the water while the forth equation is related to the water incompressibility condition. The two last equations of system (1) represent the mass and horizontal momentum balance equations for the granular layer for which the pressure is assumed to be hydrostatic and a Coulomb-type friction law is used with a the friction coefficient μ(see Section 2.3). In system (1), ∇is the gradient vector in the horizontal (x,y) plane. When hydrostatic pressure is assumed for the water layer, eq. (1) are reduced to: ∂th1+∇·(h1u1)=0, ∂tu1+u1·∇u1+∇(g(b+h1+h2)) =Sf,1, ∂th2+∇·(h2u2)=0, ∂tu2+u2·∇u2+∇(g(b+h2+(1−κ)h1)) =Sf,2.(5) The non-hydrostatic and hydrostatic versions of HySEA will be used and compared in the following sections for laboratory and field-scale simulations. It is well known that hydrostatic models (i.e. shallow-water) predict an accurate phase velocity only for values of kH 1, where His the water depth and k=2π λis the wavenumber [in practice when H λ<0.05 (Le M´ ehaut´ e1976), i.e. for kH <0.1π]. Macias et al.(2021a) showed that the HySEA multilayer non-hydrostatic model makes it possible to extend the validity range of the phase velocity but, for one-layer models (i.e. depth-averaged), the dispersive relation deviates from full potential theory for kH >2.5 (see their table 1). In particular, for small landslides occurring far from the shore and at great depth, it is difficult to accurately simulate the first waves emerging on the top of the landslide. However, waves calculated far from this non-valid region are still accurate. Multilayer models strongly improve these results but are out of the scope of this paper. Let us now continue to simplify the model to show how it can be reduced to the extended version of the SHALTOP model (Section 2.2) that will be used to quantify the error related to topography effects in HySEA. If we assume that the wavelength of the free surface perturbations surface are long compared to the water depth: ∇(b+h1+h2)=0,(6) the velocity equation for the granular layer can be simplified using ∇h1=−∇(b+h2). We then obtain the reduced system: ∂th1+∇·(h1u1)=0, ∂tu1+u1·∇u1+∇(g(b+h1+h2)) =Sf,1, ∂th2+∇·(h2u2)=0, ∂tu2+u2·∇u2+∇(gκ(b+h2)) =Sf,2. (7) Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 800 P. Poulain Finally, if we neglect shear friction between the water and the granular layer (Sf,1and Sf,2), the evolution of the granular layer is decoupled from the evolution of the fluid layer. Therefore, the evolution of the granular layer can be obtained from the solution of the system: ∂th2+∇·(h2u2)=0, ∂tu2+u2·∇u2+∇(κg(b+h2)) =−κgμu2 ||u2||.(8) This corresponds to eqs (13)and(14) implemented in SHALTOP (Section 2.2) when substituting κgfor g,h2for hand u2for uwith the simplifications c=1,s=0,H=0 in SHALTOP (see next section). When the friction between the water layer and the avalanche is neglected, then, the only effect the water plays is a role in the acceleration terms—divide eq. (8)byκ. 2.2 The SHALTOP numerical model The code SHALTOP solves the depth-averaged equations for dry granular flows on a complex 2-D topography b(x,y) (Bouchut et al.2003; Bouchut & Westdickenberg 2004; Mangeney-Castelnau et al.2005; Mangeney et al.2007a). Contrary to HySEA, the granular flow is assumed to be thin in the direction normal to the topography compared to its downslope extension. The equations are depth-averaged along this normal direction. SHALTOP has been shown to reproduce experimental and natural granular flows, with the limitations inherent to the thin-layer approximation (e.g. Favreau et al.2010; Moretti et al.2012; Peruzzetto et al.2019). Peruzzetto et al.(2021) showed that curvature effects due to a complex topography can drastically change the dynamics and deposit of natural granular flows. They are taken into account in SHALTOP through the curvature tensor of the topography b(x,y): H=∂2b ∂x2 ∂2b ∂x∂y ∂2b ∂x∂y ∂2b ∂y2.(9) The normal to the topography at each point is calculated through a 3-D unit upward normal vector [we added an arrow on vectors that are not in the (x,y)plane]: En =−∇b 1+||∇b||2,1 1+||∇b||2=(−s,c),(10) where the scalar c=cos(θ) is the cosine of the angle between the vertical direction and the normal En. The flow is described by: h≥0,u,(11) where his the thickness of the avalanche in the direction normal to the topography, u=(ud,ut)(. dfor downslope and .tfor transverse) is a parametrization of the velocity defined in order to simplify the equations with topography, and uis the velocity norm (see Mangeney-Castelnau et al.2005, for more details). In the horizontal Cartesian coordinate formulation, the mass and momentum balance equations model can be expressed as: ∂t(h/c)+∇·(hu)=0 (13) ∂tu+cu.∇u+1 cId −sst∇(g(hc +b)) =−1 cutHus+1 cstHuu−gμcu c2u2+(s·u)2)1+utHu gc + ,(14) where his the granular thickness in the direction normal to the topography. The subscript +stands for the positive part, x+=max(0,x)—see Bouchut and Westdickenberg (2004) and Mangeney-Castelnau et al.(2005) for more details. The signature of the thin-layer approximation related to the topography can be observed in the presence of cand sin the equations, and, in the right-hand side of eq. (14), in the first two curvature terms and in the last term representing the Coulomb-type friction law. SHALTOP has been designed to simulate dry granular flows. However, through the use of reduced gravity, we model buoyancy even if viscous drag effects are not included. Then, replacing gravity gby κgmakes it possible to obtain the submarine flow equations when the drag friction that couples the avalanche and the water column is neglected [see eq. (8) of the simplified HySEA equations in Section 2.1]. Note that in SHALTOP, κis taken to be equal to eq. (4) for the immersed and dry flow entering water cases, even though in the last case it should have been equal to 1 before the material enters water. This leads to a slower velocity of the aerial part of the collapse for the dry flow entering water case. 2.3 Friction law In their seminal work, Savage & Hutter (1989) used a Mohr–Coulomb internal rheology with a constant Coulomb basal friction coefficient μ. However, this simple law does not allow reproduction of the steady uniform flows observed in the laboratory for a given range of slope angles. Based on these experiments involving spherical glass beads flowing on a rough bed, Pouliquen (1999) showed that: Fr =u ghcos (θ)=βh hstop (θ),(15) Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 Model for landslide-generated tsunamis 801 where uand hare the mean velocity and thickness of the flow, Fr is the granular flow Froude number, hstop(θ) is the thickness of the deposit left after a steady uniform flow on a slope inclined at an angle θand the β=0.136 is a constant for glass beads. This observation was used by Pouliquen (1999) to propose an empirical friction law valid for steady uniform flows, with a friction coefficient depending on the thickness and velocity of the flow. Pouliquen & Forterre (2002) extended this law based on experimental measurement of the thickness hstop(θ) and of the angle θstart(h) at which a layer of thickness hstarts to flow. This law was assumed to be valid for the whole range of possible thicknesses and Froude numbers even though the experimental data only concerned steady and uniform flows. The friction law involves three friction angles δ1,δ2and δ3, the constant βand a fitting parameter Lthat is of the same order of magnitude as the grain diameter. The angle δ3 corresponds to the asymptote of the curve θstart(h). Depending on the value of the Froude number, the flow is assumed to be in a dynamic (Fr >β), intermediate (0 <Fr <β) or static regime (Fr =0) and the friction coefficient can be written as: μ(h,Fr)=⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ μ1+μ2−μ1 1+hβ LFr ,Fr ≥β, Fr βξμ1+μ2−μ1 1+h L −μ3−μ2−μ1 1+h L+μ3+μ2−μ1 1+h L ,0<Fr <β, min μ3+μ2−μ1 1+h L ,|tanθi−∇h|,Fr =0, (16) where μ1=tan δ1,μ2=tan δ2,μ3=tan δ3,iis the unit vector in the x-direction and Lis the typical grain diameter. Note that the friction law in the intermediate regime is given by a power law extrapolation between the friction laws in the static and dynamic friction regimes. Simulation results have been shown to be poorly affected by the value of the empirical parameter ξas soon as ξ<10−2(Pouliquen & Forterre 2002; Mangeney et al.2007a). As in these studies, we choose ξ=10−3. Several studies have focused on these friction laws and possible extensions to reproduce laboratory experiments on erosion/deposition waves (Mangeney et al.2007b; Edwards & Gray 2014; Edwards et al. 2017,2019;Russellet al.2019) or self-channelling flows and levee formation (Mangeney et al.2007a; Rocha et al.2019) but the behaviour of granular material in the intermediate regime is still largely unknown. This μ(I) type rheology has been first derived from dry granular flows and further extended to account for grain-fluid effects by changing the inertial number or even by introducing dilatancy (e.g. Pailha & Pouliquen 2009). Introducing these grain fluid effects is beyond the scope of this paper. In this study, we use the Pouliquen & Forterre (2002) friction law in HySEA and SHALTOP to simulate subaerial, immersed and dry entering water granular flows observed at the laboratory scale (Section 4) and field scale (Section 5). 3 EXPERIMENTAL DATA Landslide-generated tsunamis have been benchmarked. At the laboratory scale, benchmarking mainly focused on comparison between simulated and recorded water waves on gauges (Kirby et al.2022). Indeed, experimental data rarely provide both the detailed time evolution of the granular mass and the water waves. For example, only pictures of some granular flows are provided in the NTHMP Landslide Tsunami Benchmarking Workshop (Kirby et al.2022) even though laboratory experiments on water waves generated by granular flows are more and more common (e.g. granular column collapses on a horizontal bed by Robbe-Saule et al.2021, Sarlin et al.2021 or gas-fluidized granular collapse by Bouguoin et al.2020,2021). Furthermore, benchmark exercises involve a few experiments in different settings and allow the fitting of model parameters to reproduce the water waves. The fitted parameters are found to be different for each experiment, with values that are not always physically meaningful, such as for example the low friction coefficients used by Macias et al.(2021b) in the Pouliquen and Forterre friction law. To quantify the error in HySEA simulations of granular flows with independently measured friction coefficients (i.e. without calibration), we process a series of laboratory experiments on granular masses released on sloping to horizontal beds performed by Viroulet et al.(2014, 2019) for water gauges in entering water and dry experiments; unpublished data for immersed experiments and analysis of dynamics for all three configurations). These new experiments provide a unique data set of water waves generated by granular flows together with the time change of the granular mass profiles, conducted in the same experimental set-up for different slope angles, initial conditions (subaerial, immersed and dry flow entering water) and grain diameters. These experiments represent demanding conditions for the models because (i) the shallow-water approximation is not always satisfied for the granular mass (especially at the beginning), (ii) the wave generation area has characteristics out of the range of validity of the dispersive model and (iii) the bottom slope is far from horizontal and therefore the horizontal/vertical reference frame is not appropriate to describe the granular flow. As a result, the error in HySEA results when simulating these experiments can be considered in the upper range, making it possible to clearly highlight the limits of such models. Quantification of this error is essential for analysing the results obtained when applying the model to simulate field-scale landslides and tsunamis where characteristic dimensions can be, at some locations and/or times, outside of the limits of model assumptions. 3.1 Experimental setup The experimental set up consisted of a 2.20 m long, 0.20 m wide and 0.40 m high wave tank with a reservoir filled with 2 kg of glass beads of diameter d=4 mm and density ρ=2500 kg m−3that is initially released on the slope (Fig. 1); the set-up is similar to the one described in Viroulet et al.(2014) where details can be found. Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 802 P. Poulain 40cm 40cm 220cm Water free surface Gauge 2 Gauge1 Lifting gate Glass beads 30cm 45cm h h h = 0.15 m Hi Li Camera (b) (b) 1i 1 2 Figure 1. Experimental set-up: (a) initial set-up with relative position of the gauges (1, 2 and 3) and initial mass. (b) During the experiment showing the thickness of both the granular layer and water column and the position of the camera. Several experiments were performed with different slope angles. For each slope angle, depending on the initial water depth, the granular collapse was either dry, fully immersed or dry entering water (Fig. 2). A layer of the same material as the collapsing beads was glued on the slope and on the bottom of the tank. A vertical gate was used to maintain the initial granular materials of maximum thickness Hi=0.15 m from the surface of the beads glued on the bottom of the tank (Fig. 1). This gate was opened vertically using a system of mass and pulley. The velocity of the gate removal was 1.2 m s–1, lasting about 0.11 s for a granular flow duration of 1.3–2.2 s. In this study, we focus on experiments obtained on two slope angles (35◦and 45◦) with three initial water heights above the surface of the particles glued on the bottom of the tank (0, 0.15 and 0.30 m). An additional experiment was studied for the dry case with d=1.5 mm. The evolution of the granular collapses and the generated waves was recorded using a high-speed camera at a hundred frames per second and we extracted their profiles by image processing every 0.1 s. The error related to the manual picking (each measurement is repeated three times), image distortion, set-up, etc. is estimated to be about 1.5d≈6mmford=4 mm. The time evolution of the water free surface during the propagation was also recorded by several wave gauges. Table 1summarizes the main characteristics of these experiments. 3.2 Granular flow The dynamics and deposit as well as the generated water waves obviously depend on the configuration (dry, immersed or dry flow entering water) as observed for example on the 35◦slope angle (Figs 2and 3a–c). Note that the gate affects the shape of the released mass that is entrained by the gate during the initial instants when it is opened (Fig. 3) as was also observed in granular collapse experiments on inclined planes (see e.g. figs 6a and c of Mangeney et al.2010). In the dry case, the granular mass reaches the bottom of the tank about 0.4 s after the release (Fig. 2c). The granular flow is relatively thin on the slope (2–3 cm). After reaching the bottom, the flow slows down (Fig. 2c) and stops at about 40 cm from the gate (Figs 2dand3a). In the immersed case (Fig. 2e), the flow is slower and thicker (around 5 cm) because of strong interactions with water (Figs 2f–g) and reaches the slope break 0.6 s after the opening of the gate. While flowing down the slope, the flow front becomes thicker than the rest of the flow (Fig. 2g) because of dilatation effects (Section 3.3). The front deposit is located 5 cm further than in the dry case (Figs 2hand3b). The time-series of the granular profiles (Fig. 3) show that the flow is faster on the slope in dry conditions than in the immersed case because of Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 Model for landslide-generated tsunamis 803 a b c d e f g h i j k l 0.2 s 0.4 s 0 s 2.0 s 0.4 s 0 s 2.0 s 0.2 s 0.4 s 0 s 2.0 s 0.2 s Dry ( = 35°) Immersed ( = 35°) Entering Water ( = 35°) 0 cm153045 0 cm153045 0 cm153045 0 cm153045 0 cm153045 0 cm153045 0 cm153045 0 cm153045 0 cm153045 0 cm153045 0 cm153045 0 cm153045 Figure 2. Snapshots of three different experiments (Table 1) showing the time evolution of the granular collapse in dry (a–d) (Experiment 1 in Table 1), immersed (e–h) (Experiment 2 in Table 1) and dry flow entering water (i–l) (Experiment 3 in Table 1) configurations. The experiments are performed on a bed with a constant slope θ=35◦followed by a horizontal plane. reduced gravity for underwater flow, whereas the deposits extend further in the immersed condition because of less deceleration at the end of the flow. This may be due to the strong contraction of the granular mass (see Section 3.3) in the immersed case as it reaches the flat bottom, inducing an increase in pore pressure, thus reducing the granular friction force (e.g. Rondon et al.2011; Pailha & Pouliquen 2009; Bouchut et al.2016b,2017). The deposit shapes are also different with more mass located near the front in the immersed case. Finally, for the case of the dry flow entering water (Fig. 2i), the front is strongly decelerated and rounder as the mass enters the water (Fig. 2j) but it then accelerates sufficiently to become similar to the fully immersed case (Fig. 3). The flow velocity is still however lower than in the dry case. During its propagation down the slope, the granular front is thick (3–5 cm) and the flow surface has a complex shape (Fig. 2k). When the flow reaches the bottom, the flow decelerates more rapidly than in the immersed case, leaving a significant part of the deposit upslope (Fig. 2l). This larger deceleration in the dry mass entering water case may be due to the smaller dilatation of the flow during its spreading (Section 3.2.2) leading to subsequent smaller compaction, inducing a smaller increase in pore fluid pressure and thus a higher friction force on the grains (Bouchut et al.2016a,b). The shorter run-out distance (i.e. front position of the deposit) of the dry mass entering water case compared to the dry case may be due to the friction with water, especially when the mass enters it. For a 45◦slope angle, the overall behaviour is similar to the 35◦experiments (Figs A1 and 3d–f). As the slope angle increases, the surface of the granular mass becomes more irregular (Figs A1g, k compared to Figs 2g, k). At 45◦, the flow down the slope is faster for the case of a mass entering water than for the fully immersed case, both still slower than for the dry case (Fig. A2). On the 45◦slope, the travelling distance on the horizontal plane is slightly smaller than at 35◦, despite higher inertial effects, and less material remains on the slope. This decrease in run-out for 45◦compared to 35◦could be related to the longer sloping part that compensates for the higher gravitational force and higher pressure gradients related to a higher initial height Hiat 45◦(Fig. A2). The variation of the slope angle (within this range) does not greatly affect the run-out distance since the abrupt change of the topography towards the horizontal plane greatly slows the flow and is the main effect controlling the mass stopping. As for the 35◦slope, the run-out distance is the largest in the immersed case, then in the dry case, while the shortest run-out is observed in the dry mass entering water case. Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 804 P. Poulain Table 1. Summary of the experimental results and their numerical simulations with θthe slope angle and h1ithe initial water height (in cm), Hiand Lithe height and length of the initial mass (in cm), and the friction angles used to perform the simulations. Friction angles and run-out distance obtained with HySEA fitted on SHALTOP are written in brackets. We used the following parameters: Experiments 1–6 used 4 mm beads, experiments 7 and 8 used 1.5 mm beads. For experiments 1, 4, 7, 8 the water level is 0. For experiments 2, 5, the water level is 30 cm. For experiments 3, 6, the water level is 15 cm. For experiments 1-3 and 7 the slope angle is 35◦, for experiment 4– and 8 the slope angle is 45◦. For experiments 1–6 the non-fitted friction angles (δ1,δ2,δ3)usedinthe simulations are δ1=22.1◦,δ2=31.8◦,δ3=23.3◦. For experiment 7 and 8 the non-fitted friction angles (δ1,δ2,δ3) used in the simulations are δ1=24.5◦,δ2 =34.2◦,δ3=25.7◦. Fitted angles used in HySEA simulation of experiment 1–3 are δ1=27.1◦,δ2=36.8◦,δ3=28.3◦, in experiment 4-6 the fitted angles are δ1=29.1◦,δ2=38.8◦,δ3=30.3◦. Fitted angles used in HySEA simulation of experiment 7 are δ1=29.5◦,δ2=39.2◦,δ3=30.7◦,andδ1=31.5◦,δ2= 41.2◦,δ3=32.7◦for simulation of experiment 8. Initial mass height Hi, length Li (cm) Run-out distance from the gate (cm) Difference of the run-out distance between experiment and simulations (cm) Time-shift for gauges t(s) Maximum amplitudes of the generated waves at Gauge 1 (cm) Hydrostatic Non hydrostatic Experiment 1 (35◦Dry) Experiment Hi=10.3 Li=13.1 39.25 - - - HySEA Non-fitted (fitted) 52.32 (43.36) 13.07 (4.11) -- - SHALTOP 39.68 0.43 - - Experiment 2 (35◦Immersed) Experiment Hi=11.0 Li=13.2 44.94 - - 0.32 / –0.65 HySEA Non-fitted (fitted) 52.96 (43.68) 8.02 (1.26) 0.15 0.7 / –0.25 0.59 / –0.57 SHALTOP 37.44 7.5 - - Experiment 3 (35◦Entering Water) Experiment Hi=10 Li=13.53 34.7 - - 1.53 / –1 HySEA Non-fitted (fitted) 80 (62.7) 45.3 (28) 0.385 2.88 / –0.08 1.32 / –0.3 SHALTOP 37.4 2.7 - - Experiment 4 (45◦Dry) Experiment Hi=11.4 Li=10.78 38.6 - - - HySEA Non-fitted (fitted) 56.8 (43.36) 18.2 (4.76) -- - SHALTOP 37.1 1.5 - Experiment 5 (45◦Immersed) Experiment Hi=11.5 Li=10.94 41.84 - - 0.65 / –0.88 HySEA Non-fitted (fitted) 57.76 (44.0) 15.92 (2.16) 0.18 0.95 / –0.34 0.59 / –0.87 SHALTOP 37.12 4.72 - - Experiment 6 (45◦Entering Water) Experiment Hi=12.0 Li=10.88 33.78 - - 2.1 / –1.4 HySEA Non-fitted (fitted) 92.64 (72.32) 58.86 (38.54) 0.37 3.37 / –0.32 1.49 / –0.46 SHALTOP 37.12 3.34 - - Experiment 7 (35◦Dry) Experiment Hi=10.32 Li =14.1 38.56 - - HySEA Non-fitted (fitted) 48.32 (39.84) 9.76 (1.28) -- - SHALTOP 35.84 2.72 - - Experiment 8 (45◦Dry) Experiment Hi=11.8 Li=11.7 39.84 - - HySEA Non-fitted (fitted) 52.48 (40.16) 12.64 (0.32) -- - SHALTOP 36.0 3.84 - - 3.3 Dilatation of the granular material Even though most of the numerical models assume that the granular layer is incompressible, it is well known that granular flows experience both dilatation and compaction (e.g. Roux & Radjai 1997; Rondon et al.2011; Pailha & Pouliquen 2009; Bouchut et al.2016a,b, 2021). This compressibility impacts the shape and velocity of the flow and thus may be a source of error when not accounted for in numerical models, especially for granular flows under water (Rondon et al.2011; Bouchut et al.2021; Rauter 2021). We quantify these effects here to better interpret the comparison between incompressible codes and experimental results and therefore discriminate the origins of the error. Dilatation and compaction of the whole mass are estimated here by measuring, on the snapshots of the experiments, the evolution of the area of the longitudinal profile of the flowing and deposited mass Awith respect to the area of the profile of the initial mass (before collapse) Ai. At each time, we calculate the relative difference (A–Ai)/Ai. Fig. 4shows the time evolution of the dilatation/compaction of the mass for Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 Model for landslide-generated tsunamis 811 0.6 0.4 0.2 -0.2 0.8 -0.4 -0.6 -0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 0.6 0.4 0.2 -0.2 0.8 -0.4 -0.6 -0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Gauge 1 Gauge 2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 -1 -0.5 0 0.5 1 1.5 2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 -1 -0.5 0 0.5 1 1.5 2 Gauge 1 Gauge 2 Amplitude (cm)Amplitude (cm) Amplitude (cm) Amplitude (cm) Time (s) Time (s) Time (s) Time (s) 0.6 0.4 0.2 -0.2 0.8 -0.4 -0.6 -0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 0.6 0.4 0.2 -0.2 0.8 -0.4 -0.6 -0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Experimental data HySEA water elevation mf = 0.1 m-1 HySEA water elevation mf = 0 m-1 HySEA water elevation mf = 1 m-1 HySEA water elevation mf = 0.5 m-1 depth (m)depth (m) distance (m) distance (m) depth (m) depth (m) distance (m) distance (m) 0.4 s 3.0 s 0.4 s 3.0 s HySEA avalanche mf = 0.1 m-1 HySEA avalanche mf = 0 m-1 HySEA avalanche mf = 1 m-1 HySEA avalanche mf = 0.5 m-1 0.1 0.05 0.4 0.3 0.2 0.1 0.1 0.05 0.4 0.3 0.2 0.1 Experiment error Experiment avalanche 0.6 0.5 0.4 0.3 0.2 0.1 0.1 0.05 0.6 0.5 0.3 0.2 0.1 0.1 0.05 0.4 (a) (b) (c) (d) (e) (f) (g) (h) = 0.15 s = 0.15 s = 0.37 s = 0.37 s Immersed ( = 35°) Entering Water ( = 35°) Immersed ( = 35°) Entering Water ( = 35°) Figure 7. Impact of the avalanche/water friction coefficient (mf) on the evolution of the free surface from the experiments and the simulations at gauges 1 (a, c) and 2 (b, d) for the non-hydrostatic HySEA simulations with a fitted friction coefficient for a 35◦slope configuration. Simulation waves have been synched in time to the experiment. The time-shift tis represented by the grey arrow. Impact of the avalanche/water friction coefficient (mf) on the evolution of an avalanche on a 35◦slope for immersed (e,f) (Experiment 2 in Table 1) and dry flow entering water (g, h) (Experiment 3 in Table 1) configurations. The blue lines are the evolution of the water free surface of the experiments. The red, green, light blue and purple lines represent the evolution of the water free surface of the simulations. HySEA one-layer model are not always in their validity range (Macias et al.2021a). This is in particular the case in the generation area where small wavelengths are generated (Figs A6f–h and j–k, A8f–h and j–k). However, our simulations show that theses waves are still confined in the generation area during the simulation and only slightly pollute the longer wavelength waves recorded further away at the gauges. Indeed, non-hydrostatic simulations reproduce relatively well the waveforms at the gauges despite these spurious waves trapped in the generation area. Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 812 P. Poulain Consequently, the comparison between simulations and laboratory experiments demonstrates the limitation of the codes and makes it possible to quantify, in simple situations, the error in the calculated granular flow and deposits and in the generated water wave. With an estimate of the error, it is possible to better interpret the simulation of geological events, in particular in view of quantifying the associated hazards. Let us now compare the hydrostatic and non-hydrostatic approximation of the water layer in the natural context of the Montagne Pel´ ee volcano (Martinique). 5 APPLICATION TO A GEOLOGICAL CASE: INSIGHT INTO MONTAGNE PEL ´ EE LANDSLIDES AND GENERATED WAVES 5.1 Geological description and data preparation 5.1.1 Flank-collapse history and unrest at Montagne Pel´ ee volcano Martinique is located in the central part of the Lesser Antilles Arc. Volcanic activity on Martinique covers a period from the early Miocene to present. Montagne Pel´ ee volcano in northern part of Martinique is an explosive-type active volcano that has experienced numerous magmatic and phreatic eruptions including the 1902–1905 and 1929–1932. During the last 130 Ma years, three flank collapses reshaped the western part of the Montagne Pel´ ee volcano (Fig. 8a) with volumes estimated from 1.8 km3for the last event at 36 ka BP (DAD3) to 25 km3for the oldest event (Le Friant et al.2003a; Boudon et al.2007; Le Friant et al.2015,2019, Brunet et al.2016,Solaroet al.2020). The volcanic debris-avalanche component of these landslide deposits is restricted to proximal areas and tends to stop at the break in slope in the bathymetry whereas the distal part of the landslide deposit, called SLD (Submarine Landslide Deposit), results from widespread failures of pre-existing low topographic gradient seafloor sediments (Le Friant et al.2013,2015,2019). In Fig. 8(a), DAD3 is the debris avalanche deposit related to the last flank-collapse event (36 ka) whereas the morphological bulge (the three blue lines) result from the accumulation of the two older debris avalanche deposits related to the two older flank-collapse events (Brunet et al.2016;LeFriantet al.2015,2019). The dynamics of these debris avalanches and their deposits remain unclear. Brunet et al.(2017) discussed the landslide dynamic by modelling the emplacement of DAD3 using SHALTOP and HySEA models. They performed sensitivity tests with SHALTOP and HySEA and concluded that large flank collapses (about 20 km3) probably occurred in several events with successive volumes smaller than 5 km3entering the sea. In their simulations, they did not took into account the presence of an erodible bed, the precise reconstruction of the topography and the variation of the friction angle with collapsed volumes. The last magmatic eruption produced a large lava dome inside an explosive crater that was excavated in the lava dome from the 1902–1905 eruption. The 1929–1932 lava dome was built on the border of the summit caldera (Crat` ere de l’Etang Sec). Its prominent talus slope extends to the south several hundred metres below the rim. The activity of Montagne Pel´ ee volcano decreased significantly after the end of the last eruption in 1932 (Fig. 8b). Fumarole activity has stopped since 1970. However, since April 2019, volcanic seismic activity has been increasing. In November 2020, the OVSM-IPGP recorded a new marked increase in the seismicity and the appearance of long-period earthquakes and volcanic tremors (OVSM-IPGP 2020-2021). Observation at the end of December 2020 of a zone of degraded vegetation on the southeastern flank of Mount Pel´ ee has suggested that diffuse CO2of magmatic origin may be degassing. These observations reflect an increase in the activity of the hydrothermal system of Montagne Pel´ ee. However, the OVSM-IPGP has not recorded any evidence of large-scale deformation or high-temperature sulfur gas emission. In this context of reactivation, the past flank-collapse history of Montagne Pel´ ee, that is its structure with the last dome of 1929 located within and overlapping the rim of the Etang Sec caldera, with a steep non-buttressed talus slope rising nearly 800 m above the surrounding topography, could favour instability, as observed on other volcanoes (Harnett et al.2019). Such failure could depressurize the magmato-hydrothermal system, trigger laterally directed explosions (blasts) and generate the emplacement of potentially tsunamigenic debris avalanches. To conduct a preliminary assessment of this problem with a numerical simulation of failure, we first perform sensitivity tests on the collapsed volumes (by adapting the friction angles used with volumes) based on the current geological knowledge of Montagne Pel´ ee volcano. We also test the effects of an erodible substratum. These sensitivity tests allow us to provide a range of run-out distances and sea surface elevations as a function of the different collapsed volumes. We then investigate two (non-exhaustive) deterministic scenarios of potential collapse (1929 Dome scenario and Dome slope scenario) related to the current heightened seismicity and volcanic unrest (Fig. 8). 5.1.2 Collapse scenarios and data preparation The stability of Montagne Pel´ ee volcano is difficult to assess because of the lack of geotechnical data. To perform sensitivity tests and calibrate the model, we use the geometric and geological constraints of the three large flank-collapses identified on Montagne Pel´ ee volcano to define the scars and volumes of the collapse scenarios (Le Friant et al.2003a,b; Brunet et al.2016,2017). To reconstruct pre- and post-topography of Montagne Pel´ ee volcano, we use the topography provided by IGN and the bathymetry collected during the AGUADOMAR expedition (Deplus et al.2001). The topography and bathymetry were combined to build a 100 m resolution Digital Terrain Model (DTM). The DTM was modified to reconstruct the topography before the flank collapses by replacing the submarine deposits with a smooth slope. This DTM was then modified to reconstruct the pre- and post-topography for different volumes ranging from 1.8 to 16.2 km3(Fig. 9a). Note that in this Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 Model for landslide-generated tsunamis 813 (a) (b) Figure 8. (a) Terrestrial digital elevation model (from IGN) with Montagne Pel´ ee volcano (red star) and bathymetric map from the AGUADOMAR (1999) and CARAVAL (2002) cruises, adapted from Brunet et al.(2017). The three horseshoe-shaped structures are represented on land (Le Friant et al.2003a,b). The extents of submarine landslide deposits are plotted with black lines from Brunet et al.(2016). The debris avalanche deposit DAD3 related to the last flank collapse event is in blue and the three blue contour lines represent the bulge associated with older debris avalanche deposits (associated with older flank collapses). The three diamonds represent the position of the gauges used to monitor the sea surface elevation generated by the debris avalanches. Inset: Lesser Antilles Arc map on the left-hand side, with Martinique Island and both the active arc (red dotted line) and older arc (orange dashed line), (figure adapted from Brunet et al.(2017)). (b) Terrestrial digital elevation model (from IGN) showing the Montagne Pel´ ee volcano and the old flank-collapse structures in purple. The two lava domes of 1929 and 1932 are indicated (pink and blue) inside the ‘Etang Sec’ caldera (yellow line). The two potential flank-collapse scenarios related to the reactivation of the volcano are indicated with horseshoe-shaped structures in red and orange. The areas of ongoing seismic activity are contoured in blue and orange and the degraded vegetation area is contoured in yellow. Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 814 P. Poulain 1400 1300 1200 1100 1000 900 800 0 200 400 600 800 1000 1200 1400 1600 700 800 900 1000 1100 1200 1300 0 200 400 600 800 1000 1200 1400 Distance (m) Heigth (m) Distance (m) Heigth (m) 1929 dome scenario Dome slope scenario V V 1400 1200 1000 800 600 400 200 0 -200 -400 -600 -800 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 V V V V Distance (m) Heigth (m) C1 C1’ C3 C3’ C4 C4’ C2 C2’ C2 C2’ 0 5 10 15 20 25 30 -2500 -2000 -1500 -1000 -500 0 500 1000 1500 0 5 10 15 20 25 30 35 40 Slope angle (°) Heigth (m) Distance (km) mean slope angle 3.7° 14.1° 9.2° 21.4° 15.4° (a) (b) (c) (d) Figure 9. (a) Collapse geometries of scenarios with different volumes for the sensitivity tests. (b) Slope variations along a profile crossing the volcano and submarine flanks. The mean slope angle is shown with dotted lines averaged along sections between five main slope breaks. (c) Geometry of the 1929 Dome scenario. (d) Geometry of the Dome slope scenario. study, the reconstruction of the first event result in an estimated volume of 16.2 km3for the largest event, which is a bit smaller than the estimation of Le Friant et al.(2003a). The reconstruction of the pre-event topography was done using the isolines from the area of the volcano outside the collapse scar. Four different flank collapse scenarios involving volumes of (a) 1.8 km3, (b) 16.2 km3and (c) 1.8 km3subdivided into three successive subevents of 0.6 km3, and (d) 16.2 km3subdivided into three successive subevents of 6.9 km3, 5.9 km3and 3.7 km3, have been investigated. Finally, to provide a first assessment of the potential instability of the edifice in the context of the unrest at Montagne Pel´ ee since 2019, we focus our attention on the 1929–1932 lava dome located within the summit caldera. Given the current unrest conditions related to the location of volcano-tectonic seismicity and the zone of damaged vegetation, we consider that the scenario involving the partial collapse of the 1929–1932 dome is the most critical to investigate and to model (Fig. 8b). We consider two scenarios. The first scenario concerns the 1929–1932 lava dome (1929 Dome scenario) with an estimated volume of 57 ×106m3. The second scenario concerns the collapse of part of the dome slope (Dome slope scenario, estimated volume 9.8 ×106m3; Figs 8b, 9c and d). We use a 30 m resolution DTM provided by Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 Model for landslide-generated tsunamis 815 Table 2. Main characteristics, friction angles and results of the simulated scenarios for volcano landslides at Montagne Pel´ ee. SHALTOP Volu m e (km3) Run-out distance (km) Offshore distance (km) Covered surface (km2) Deposit maximum thickness (m) Friction angles δ1,δ2,δ3(◦) Single events - 1.8 39.4 33.4 184.2 29.3 2, 12, 4 HySEA Volume (km3) Run-out distance (km) Offshore distance (km) Covered surface (km2) Deposit maximum thickness (m) Friction angles δ1,δ2,δ3(◦) Single events - 0.01 2.5 0 3.4 8.8 4.7, 14.7, 6.7 - 0.05 10.3 4.5 22.5 48 3.7, 13.7, 5.7 - 0.6 26.5 21.4 122.5 29 2.4, 12.4, 4.4 - 1.2 31.3 26 168.6 30.1 2.1, 12.1, 4.1 - 1.8 38.6 32.9 343.4 28.3 2, 12, 4 - 6.9 39.3 33.3 522.4 51.2 1.4, 11.4, 3.4 - 12.8 46.5 41.3 782.1 54.9 1.1, 11.1, 3.1 - 16.2 50.4 44.4 1105.7 53.3 0.9, 10.9, 2.9 Successive subevents (erodible deposits) 1.8 km3subevent 1 0.6 26.5 21.4 122.5 29 2.4, 12.4, 4.4 subevent 2 0.6 27.7 22.5 148.1 33.3 2.4, 12.4, 4.4 subevent 3 0.6 29.2 23.5 183.3 44.6 2.4, 12.4, 4.4 16.2 km3subevent 1 6.9 39.3 33.3 522.4 51.2 1.4, 11.4, 3.4 subevent 2 5.9 41.3 35.3 673 57.7 1.4, 11.4, 3.4 subevent 3 3.7 41.4 35.4 688.3 69.9 1.6, 11.6, 3.6 Successive subevents (non-erodible deposits) 1.8 km3subevent 1 0.6 26.5 21.4 122.5 29 2.4, 12.4, 4.4 subevent 2 0.6 27 21.5 133.5 43 2.4, 12.4, 4.4 subevent 3 0.6 28 22.3 155.2 53 2.4, 12.4, 4.4 16.2 km3subevent 1 6.9 39.3 33.3 522.4 51.2 1.4, 11.4, 3.4 subevent 2 5.9 39.5 33.5 564 71.7 1.4, 11.4, 3.4 subevent 3 3.7 39.6 33.6 578.9 91.3 1.6, 11.6, 3.6 Martinique unrest scenario Dome slope scenario 0.0098 5 0 4 27.9 4.7, 14.7, 6.7 1929 Dome scenario 0.057 13.9 8.2 33.9 57 3.7, 13.7, 5.7 IGN, the Litto3D data (IGN-SHOM) and the HOMONIM DTM provided by SHOM to reconstruct the topography and the bathymetry and perform numerical simulation. 5.2 Sensitivity tests on volumes, number of events and erodible substratum 5.2.1 Debris avalanches The values of the friction angle to be used is still the subject of debate. Lucas’s empirical law derived from dry debris avalanches (Lucas et al.2011) has yielded a good fit for some volcanic contexts (Zhao et al.2015) but was not able to reproduce the deposit in Martinique (Brunet et al.2017) or in Guadeloupe (Perruzzeto et al.2019). In these last two cases, a smaller friction coefficient was needed, suggesting that the observed greater mobility of the volcanic debris-avalanches could be explained by the presence of meteoric water, pressurized hydrothermal fluids and highly hydrothermally altered material in the sliding mass. To determine the friction angles to be used in this study, we start from the Lucas’s empirical relationship (Lucas et al.2011): tan (δ)=V−0.0774 (17) and combine it with the results of Brunet et al.(2017), who have calibrated the angles for the 1.8 km3volume collapse of Montagne Pel´ ee with geological constraints. For the 1.8 km3, Lucas’s empirical law suggests that δ=11◦. However, Brunet et al.(2017) show that the deposits are best reproduced using δ=6◦for a Coulomb’s type law, probably because of the presence of water. We therefore modify Lucas’s empirical law by introducing a constant αto fit the result of Brunet et al.(2017): tan (δ)=αV−0.0774 (18) with α=0.54. We then deduced the δangles for each volume reported in Table 2. Given that Brunet et al.(2017) reproduced the 1.8 km3collapse of Montagne Pel´ ee with δ=6◦for the Coulomb friction law and with δ1=2◦,δ 2=12◦and δ3=4◦for the Pouliquen and Forterre friction law, Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 816 P. Poulain SHALTOP HySEA (a) (b) 5 10 15 20 25 30 35 0 010203040 Distance (km) Distance (km) 5 10 15 20 25 30 35 0 010203040 Distance (km) Distance (km) 5 10 15 20 25 30 35 0 010203040 Distance (km) Distance (km) HySEA (d) 5 10 15 20 25 30 35 0 010203040 Distance (km) Distance (km) HySEA (c) 0 5 10 15 20 25 30 Thickness (m) 5 10 15 20 25 30 35 0 010203040 Distance (km) Distance (km) HySEA (e) 0 5 10 15 20 25 30 Thickness (m) 1 = 2°, 2 = 12°, 3 = 4° 1 = 2°, 2 = 12°, 3 = 4° mf = 0 1 = 2°, 2 = 12°, 3 = 4° mf = 1x10 n = 0.02 m-1/3 s 1 = 2°, 2 = 12°, 3 = 4° mf = 1.5x10 n = 0.02 m-1/3 s 1 = 2°, 2 = 12°, 3 = 4° mf = 2x10 n = 0.02 m-1/3 s V = 1.8 km V = 1.8 km V = 1.8 km V = 1.8 km V = 1.8 km n = 0 m-1/3 s 14°44.7’N 61°18’W 14°44.7’N 61°18’W 14°44.7’N 61°18’W 14°44.7’N 61°18’W 14°44.7’N 61°18’W Figure 10. Final deposits calculated for the 1.8 km3flank collapse with (a) the SHALTOP model and (b–e) the HySEA model. mfrepresents the drag coefficient between the water and the granular avalanche and nis the Manning coefficient. Blue lines represent DAD3 observed on the submarine flank of the volcano. Dashed white lines represent the result of the simulation conducted by Brunet et al.(2017). we calculate δ1from equation (18) by imposing: δ1=δ−4◦(19) and calculate δ2and δ3using the formula defined in Section 4.1. All the angles used in the series of simulations are reported in Table 2. In this section, we first test the influence of one single event versus multiple-stage events and the influence of a potential erodible substratum using the friction angles defined previously. Single slide events To compare HySEA and SHALTOP for the 1.8 km3flank collapse, we use the following angles: δ1=2◦,δ 2=12◦and δ3=4◦,and we first set to zero the friction between the water and granular layer and the Manning coefficient in HySEA (Figs 10aand b). SHALTOP gives a more accurate representation of the deposit while the run-out distance using HySEA is longer. As seen from the experiment, when the slope is smaller than 15–20◦, (here the average slope is about 7.6◦, Fig. 9b), we do not have to fit the friction angle in HySEA to get closer to the deposit obtained with SHALTOP. Here, we then just introduce a drag in HySEA between the water and granular layer with a coefficient (mf) and the Manning coefficient (n). We used a Manning coefficient n=0.02 m−1/3 s since our laboratory-scale study showed that it barely affects the water wave and even less the granular flow (Fig. A11). We test different values of the drag coefficient mf=1×10−5, 1.5 ×10−5,2×10−5m−1(Figs 10c, d and e). Note that the drag coefficient mfhas dimension of m−1so that its value at the field scale should be about 10−3smaller than its the laboratory-scale value. Including a drag coefficient of 1.5×10−5m−1and a Manning coefficient of 0.02 m−1/3 s, we obtain the best deposit with a run-out distance similar to the observed deposit (Fig. 10d). We then compare the flow dynamics between SHALTOP (δ1=2◦,δ 2=12◦and δ3=4◦)andHySEA(δ1=2◦,δ 2=12◦and δ3=4◦with a friction coefficient and a Manning coefficient of 0.02 and 1.5×10−5), respectively. With SHALTOP, most of the volume of the slide initially flows rapidly downslope in the horseshoe structure, then to the coastline and offshore (Fig. A12a). A relatively small part of the slide spreads to the north whereas the main part of the slide bounces off the south side of a pre-existing erosive channel. The maximal run-out distance, reached after 300 s, is 33.4 km from the coastline (Fig. A12c). When the front stops, some material is still moving at the back of the front. This material then accumulates in the front and centre of the deposits and a tail develops behind the main mass at between 300 and 600 s (Figs A12c and d). As expected, the results are similar even though the flow simulated with SHALTOP is slightly faster (Figs A12a–h). The maximum thickness of material reaches 29.3 and 28.3 m for SHALTOP and HySEA, respectively, which is coherent with observations from Le Friant et al.(2003a, b). Results from numerical simulations using SHALTOP and HySEA confirms that HySEA can be used as it gives a good run-out estimation of the flow. Note however that SHALTOP better reproduces the detailed morphology of the deposit with HySEA displaying a larger spreading of the deposits on the steep submarine slopes (Figs A12e–h). We use HySEA in the following section. To investigate additional scenarios other than the 1.8 km3event, we use HySEA to calculate the run-out distances for 7 volumes (ranging from small flank collapse volumes to the largest volume reconstructed on Montagne Pel´ ee volcano): V=10 ×106m3,V=50 ×106m3, V=0.6km 3,V=1.2km 3,V=6.9km 3,V=12.8km 3,V=16.2km 3,(Fig.11). The debris avalanche related to the smallest event Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 Model for landslide-generated tsunamis 817 10 20 30 40 50 60 30 40 50 60 70 80 90 Distance (km) Distance (km) 3 successive slides V1 = 3.7 km V2 = 5.9 km V3 = 6.9 km 3 successive slides V1 = 3.7 km V2 = 5.9 km V3 = 6.9 km 2 successive slides V1 = 3.7 km V2 = 5.9 km 2 successive slides V1 = 3.7 km V2 = 5.9 km 3 successive slides V1 = 0.6 km V2 = 0.6 km V3 = 0.6 km 3 successive slides V1 = 0.6 km V2 = 0.6 km V3 = 0.6 km 2 successive slides V1 = 0.6 km V2 = 0.6 km 2 successive slides V1 = 0.6 km V2 = 0.6 km 1 slide V = 16.2 km 1 slide V = 12.8 km 1 slide V = 6.9 km 1 slide V = 1.8 km 1 slide V = 1.2 km 1 slide V = 0.6 km 1 slide V = 50x10 m 1 slide V = 10x10 m 14°49’N 14°27.5’N 61°5’W 61°27’W Figure 11. Final deposits pathways and extensions simulated with the HySEA model for all the sensitivity tests on Montagne Pel´ ee. (10 ×106m3) does not reach the sea. The debris avalanche associated with the 50 ×106m3enters the sea and propagates over a distance of 2 km after the coastline. The 0.6 km3event deposit stopped at 21.4 km after the coastline. The dynamics of the 6.9 km3event are similar to those of the 16.2 km3event presented afterwards but with a smaller run-out distance of 33.3 km after the coastline and covering an area of 522.4 km2. Finally, a maximum run-out distance of 41.3 km after the coastline is obtained for a flank collapse of 12.8 km3. Regarding the flow dynamics of the debris avalanche associated with the largest volume (16.2 km3), most of the material of the slide flows rapidly downslope to the southwest in the horseshoe structure and then to the coastline and offshore (Figs A13a and b). The avalanche continues to flow, creating one main location of deposition (purple arrow in Fig. A13d) located at the front of the flow. A part of the material stays within the structure. The debris avalanche reaches a maximum run-out distance of 44.4 km from the coastline (Table 2) and a maximum deposit thickness of about 53.3 m. These results are not in agreement with the observations as the bulge (three blue lines in Fig. 8a), related to the debris-avalanche deposits associated with the event, is located 35 km from the coastline. One slide event versus three subslide events To test the effect of a flank-collapse occurring in successive retrogressive subevents as suggested by Brunet et al.(2017), the 1.8 and 16.2 km3volumes were each divided in three smaller volumes that slide one after the other, but the values of the friction angle were modified for each volume, contrary to the previous work of Brunet et al.(2017;Table2). First, for the 1.8 km3volume, we investigate the influence of the nature of the substratum (erodible or not) on the flow dynamics and deposit. Two sets of simulations were performed considering each previous deposit either as an erodible substratum or as a solid topography for the following event. The parameters used are summarized in Table 2. For the erodible case, Fig. A14 shows: (i) a yellow dashed line delimiting the extent of the deposit resulting from the first subevent simulation (0.6 km3) with a run-out distance of 21.4 km from the coastline after 3000 s, (ii) a green dotted line delimiting the extent of the deposit resulting from the second subevent simulation (0.6 km3) with a run-out distance of 22.5 km from the coastline after 3000 s and (iii) a coloured area representing the deposit resulting from the third subevent of 0.6 km3after 250 s. Figs A14a and b show that the third debris avalanche spreads above the previous deposit forming an erosion/deposition wave as observed for laboratory-scale dry granular flow experiments and simulations (Pouliquen & Forterre 2002; Mangeney et al.2007b; Edwards et al.2017; Viroulet et al.2019) and for field-scale simulations (Moretti et al.2012, their fig. 2d). The wave propagates beyond the front of the previous deposit, increasing the run-out distance (23.5 km from the coastline; Figs A14a, c and e). For the non-erodible substratum situation, the dynamics are different and no erosion/deposition waves are observed for the second or third subevents. The material just accumulates on the top of the thickest part (Figs A14b, d and f) and the run-out distances are smaller. Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 818 P. Poulain (a) (b) stneve-busevisseccuseerhTtneveenO 5 10 15 20 25 30 35 10 20 30 40 5 10 15 20 25 30 35 10 20 30 40 Time = 3000 s Time = 3000 s Distance (km) Distance (km) Distance (km) Distance (km) 00 2 (c) (d) 40 60 80 40 60 80 10 20 30 40 50 60 10 20 30 40 50 60 Distance (km) Distance (km) Distance (km) Distance (km) 5 0 10 15 20 25 30 10 0 20 30 40 50 60 1 = 2°, 2 = 12°, 3 = 4° mf = 1.5x10-5 n = 0.02 m-1/3 s 1 = 2.4°, 2 = 12.4°, 3 = 4.4° mf = 1.5x10-5 n = 0.02 m-1/3 s 1 = 0.9°, 2 = 10.9°, 3 = 2.9° mf = 1.5x10-5 n = 0.02 m-1/3 s V1 : 1 = 1.4°, 2 = 11.4°, 3 = 3.4° mf = 1.5x10-5 n = 0.02 m-1/3 s V2 :1 = 1.4°, 2 = 11.4°, 3 = 3.4° V3 :1 = 2.6°, 2 = 11.6°, 3 = 3.6° V = 1.8 km V = 1.8 km (0.6 + 0.6 + 0.6) V = 16.2 km V = 16.2 km (6.9 + 5.9 + 3.7) 14°44.7’N 61°18’W 14°44.7’N 61°18’W 14°49’N 61°6’W 14°49’N 61°6’W Figure 12. Final simulated deposits of the 1.8 km3collapse scenario occurring in (a) one single event or (b) three successive subevents. Blue lines represent the observed DAD3 associated with the last flank collapse event of Montagne Pel´ ee. Dashed white lines represent the result of the simulation conducted by Brunet et al.(2017). Final deposit of the 16.2 km3collapse scenario occurring in (c) one event or (d) three successive subevents (6.9 km3, 5.9 and 3.7 km3, with erodible deposits). Blue contour lines represent the observed bulge limits indicating the maximum run-out distance of the 16.2 km3flank collapse events at Montagne Pel´ ee (Brunet et al.2017). In both cases (erodible or non-erodible), the final run-out distance for the 1.8 km3flank-collapse occurring in three smaller subevents is shorter than the observed deposit (blue line in Figs A14e, f and Table 2) suggesting that this event occurred in one single event (Figs 12a and b). We then investigate, for the erodible case, the effect of successive sub-events for the maximum volume of 16.2 km3that we divided into three volumes (6.9, 5.9, and 3.7 km3). When considering one single event, the maximum thickness of material is located at the front of the deposit (53.3 m) and the final run-out distance is longer (44.4 km; Fig. 12c). The final run-out distance is smaller for three successive subevents (35.4 km versus 44.4 km; Fig. 12d). The three successive subevents lead to a maximum thickness of deposit located near the front (69.9 m). The dynamics are similar for the successive subevents 2 and 3 with the material spreading across the previous deposits, not really increasing the run-out distance. The best fit between the observations of the bulge and the results of the numerical simulations is obtained for three successive subevents. These results tend to support the hypothesis, for the 16.2 km3scenario,of a large event that occurred in smaller subevents at Montagne Pel´ ee. Therefore, the 6.9 km3seems to be a maximum volume for any subevent, but it is possible that the 16.2 km3flank collapse is composed of more than three even smaller subevents. 5.2.2 Generated waves In this section, we simulate the wave generation and propagation with the hydrostatic and non-hydrostatic versions of the HySEA model for the five different flank collapse scenarios of Montagne Pel´ ee involving volumes of 50 ×106m3, 1.8 km3, 1.8 km3subdivided into three successive sub-events, 16.2 km3, and 16.2 km3subdivided into three successive subevents. Below, we detail the results of the first three scenarios and indicate systematically the values of the sea surface elevation obtained with both the hydrostatic and non-hydrostatic models as follows: value with hydrostatic model/value with non-hydrostatic model. Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 Model for landslide-generated tsunamis 819 Figure 13. (a) Comparison of the sea surface elevation at three gauges calculated with the hydrostatic and non-hydrostatic versions of the HySEA model, for the case of a 1.8 km3collapse. The gauges location are shown in Fig. 8(a). First, the time-series of the sea surface elevation at three locations (see positions in Fig. 8a) of sea surface elevations generated by the 1.8 km3volume occurring as a single event are shown in Fig. 13. Gauge 1 (blue), in front of the flank-collapse structure, show a maximum amplitude of 104.2 m/70.4 m for the first generated wave (Fig. 13) whereas the second (red) and third (green) gauges show a first wave of 39 m/22.6 m and 18 m/7.7 m, respectively. A second wave with an amplitude of 30 m/35.9 m is estimated at the red gauge. Using the hydrostatic model, the first wave is 32 per cent higher at gauge 1, 42 per cent higher at gauge 2 and 57 per cent higher at gauge 3, than when using the non-hydrostatic model. Using the hydrostatic and the non-hydrostatic versions of HySEA gives similar results for the avalanche part as observed in the laboratory-scale study (Section 4). Then, waves generated by three successive flank collapses of 0.6 km3were simulated numerically. The first subevent generated a wave of 55 m/28.9 m on the first gauge (blue), 14.5 m/7.4 m on the second gauge (red) and 4.6 m/3 m on the third gauge (green) (Fig. A15a). The second subevent resulted in larger waves (Fig. A15a). For the last 0.6 km3subevent, the three gauges suggest waves of amplitude similar to that of the second subevent (Fig. A15a). The increasing amplitude of the waves generated by the successive events is related to their increasing velocity as the successive events start to flow higher on the slope. Note that, for the case of an erodible substratum, the wave generated by the last subevent generates a stronger wave than on a rigid substratum. The wave also displays a slightly different waveform as for instance on gauge 2 in Fig. A16. Fig. A15(b) shows the maximum water elevation at the three gauges for the 50 ×106m3collapse. The gauge (blue gauge) located in front of the flank-collapse (3 km from the coast) recorded a maximum amplitude of 8 m/4.3 m for the first generated wave at 105 s (Fig. A15b). The amplitude of the first wave significantly decreases during its propagation to reach only 1 m/0.27 m and 50 cm/20 cm at the second (red) (located close to Bellefontaine) and third gauge (green) (located 20 km from the coast), respectively (Fig. A15b). The maximum amplitude and wavelength of the generated waves for all the scenarios as a function of the volume of the event are represented in Fig. 14 for volumes varying from 50 ×106m3to 16.2 km3for hydrostatic and non-hydrostatic cases. We observed that, for the case of flank collapses on La Montagne Pel´ ee volcano, the wave amplitude increases logarithmically as a function of increasing volume whereas the wavelength increases as a function of volume according to a power law. Typically, the calculated wavelengths λvary from 1 km to more than 30 km. If we apply the validity criteria of the HySEA hydrostatic model (kH <0.314) and of the HySEA non-hydrostatic one-layer model (kH <2.5), we observe (Fig. A17) that for the shortest waves (λ=1 km), the hydrostatic assumption is valid only within a narrow band close to coastline (H<55 m) and the non-hydrostatic assumption is valid for waves propagating at depths less than 400 m. As the water depth in this case is always less than 3 km, the HySEA non-hydrostatic model is always valid for the generation/propagation of waves with wavelengths λ>7.5 km in the calculation domain. As the avalanche deposits are located at depths less than 2.7 km, the non-hydrostatic model is always valid for the generation/propagation of waves with wavelengths λ >6.8 km. For intermediate (15 km) and long wavelengths (30 km), the hydrostatic assumption is valid for wave propagation at depths below 850 and 1700 m, respectively. 5.3 Potential collapse scenarios related to ongoing volcano-tectonic activity Taking into account the ongoing unrest at Montagne Pel´ ee volcano since 2019, we considered two potential instability and collapse scenarios: the collapse of the 1929 dome with a volume of 57 ×106m3called the ‘1929 Dome scenario’ and the collapse of part of the steep talus pile of the 1929 dome with a volume of 10 ×106m3called the ‘Dome slope scenario’. Numerical simulations were performed for both scenarios using the previously calculated friction angles (Section 4.2.1). From these simulations, we estimate maximum run-out distances Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024 820 P. Poulain Figure 14. Maximum amplitude and wavelength of the generated waves as a function of the volume of the collapse scenarios simulated for sensitivity tests at Montagne Pel´ ee volcano. (a) (b) Figure 15. Final deposits calculated with HySEA for the (a) 1929 Dome collapse scenario and (b) Dome slope scenario. Older flank-collapse structures are in grey with the rims of the different structures numbered from 1 to 5. for the deposits as displayed in Figs 15(a) and (b) for simulations with δ1=3.7◦,δ 2=13.7◦,δ3=5.7◦and δ1=4.7◦,δ 2=14.7◦,δ3=6.7◦, respectively. In both cases (Figs 15a and b), the main part of the modelled debris avalanche flows southwest towards St Pierre. The main part of the flow is channeled between two older morphological structures that already channeled the block and ash flow deposits of the 1902–1905 and 1929–1932 eruptions, filling the Rivi` ere Blanche valley (Le Friant et al.2003a). A discontinuity (4 in Fig. 15) is formed by a cliff facing northwest on the left bank of Rivi` ere S` eche. The rim (2 in Fig. 15a) of the structure is a continuous cliff, facing southeast on the right bank of Rivi` ere Claire. For the 57 ×106m3event, part of the debris avalanche flows to the northwest, north of rim 3. The avalanche from the 1929 Dome scenario reaches 14 km (and thus 8 km from the coastline). In contrast, the avalanche from the Dome slope scenario reaches only 5.3 km from the volcano’s summit, without reaching the coastline. The maximum thickness of the deposits is 30 m (1929 Dome scenario) and 27.9 m (Dome slope scenario) in valleys when channeled. The area covered by the deposits ranges from 33.9 km2for the 1929 Dome scenario and 4 km2for the Dome slope scenario. These preliminary simulations show that in both cases, an instability related to the ongoing seismic-volcanic crisis could have consequences for the population. The friction angles we used here empirically account for the presence of water as discussed in Section 5.2.1. The recent observation of damaged vegetation on the flank of the volcano suggests that the hydrothermal system is active. As on other volcanoes, both phreatic and magmatic historical eruptions at Montagne Pel´ ee trigger the ascent and outflow of aquifer water and pressurized hydrothermal fluids in the host rock thereby increasing fluid saturation in the potentially unstable areas as well as reducing friction along low-strength listric and basal layers. In addition to the flow path of debris avalanches, the large simulated volumes of material deposited high upstream in the valleys will act as sediment sources for the genesis of intense debris flows (lahars) during torrential rainfall, often for an extended period of time following the initial collapse event (e.g. Peruzzetto et al.2022). They will also act as sediment barriers allowing the formation of perched water ponds that could generate, upon dam rupture, highly damaging and mobile water waves that Downloaded from https://academic.oup.com/gji/article/233/2/796/6880158 by Biblioteca de la Universidad de Sevilla user on 03 April 2024