Towards a predictive multi-phase model for alpine mass movements and process cascades

Alpine mass movements can generate process cascades involving different materials including rock, ice, snow, and water. Numerical modelling is an essential tool for the quantification of natural hazards, but state-of-the-art operational models reach their limits when facing unprecedented or complex events. Here, we advance our predictive capabilities for process cascades on the basis of a three-dimensional numerical model, coupling fundamental conservation laws to finite strain elastoplasticity. Through its hybrid Eulerian-Lagrangian character, our approach naturally reproduces fractures and collisions, erosion/deposition phenomena, and multi-phase interactions, which finally grant very accurate simulations of complex dynamics. Four benchmark simulations demonstrate the physical detail of the model and its applicability to real-world full-scale events, including various materials and ranging through four orders of magnitude in volume. In the future, our model can support risk-management strategies through predictions of the impact of potentially catastrophic cascading mass movements at vulnerable sites.


Introduction
Instabilities in mountain slopes can lead to a variety of hazardous mass movements such as falls, topples, slides, spreads, flows, and avalanches (33,77). These processes populate a vast domain of possible sizes and velocities, ranging from sub-unit volumes to the largest terrestrial mass movements on the planet (up to billion of cubic meters) and from extremely slow velocities (10 −10 m s −1 ) to extremely rapid movements (up to 10 2 m s −1 ), remaining visible in the landscape for tens of thousands of years (35,50,51). Landslides can be triggered by precipitation, seismic activity, glacial debuttressing, permafrost thawing and by human activity (29,30,56). Where large mass movements intercept anthropic exposure, catastrophic consequences for livelihoods and infrastructures may occur (25,40). Some of the most disastrous landslides result from the complex interactions between different materials and cascading processes (16,66,84). In recent history many tragedies are related to such cascading events. Amongst the most notable ones, the Vajont rockslide in Italy (52,69), the Huascaran ice avalanche in Peru (17,57), the Kolka-Karmadon ice detachment in southern Russia (31), the Piz Cengalo rock avalanche in Switzerland (3,81), and the Chamoli disaster in India (70) are only some examples. Large rock falls can entrain sediments, ice, snow or water and develop into long-reaching avalanches and debris flows (32,57). Landslides can also impact water bodies and cause outburst floods by over-topping or, in the worst case, induce the partial or total destruction of the dam (69,70,84). In general, the dynamics of process cascades is controlled by the topography and by the physical properties of its constitutive materials, with water in all its phases (liquid, ice, and snow) playing a determinant role for mobility and for the transition between different processes (65,81). Drastic changes in the cryosphere increase the likelihood of large slope failures in steep glacierised and permafrost terrain, acting on both components of the related hazards (20,32,41,61). This situation, driven by climate change and by steadily increasing anthropic pressure to remote environments, requires important advances in our capability of assessing and mitigating the hazards related to mass movements and process cascades.
The hazard related to mass movements can be quantified through mapping and can be mitigated by means of integrated hazard management with the design of protective structures and alternative mitigation measures (26,62). Hazard is the result of the probability of occurrence and the expected intensity of a given event and as such, its quantification requires the assessment of both susceptibility and impact. The operative procedure for hazard assessment consists in the definition of scenarios including release zone location and dimension, avalanche volumes and a set of appropriate model parameters. Modern hazard management strategies require quantitative assessment of the magnitude of the mass movements, i.e. their potential impact on exposed areas, which is usually synthesised by spatially distributed values of flow height and depth-averaged velocity (21). These quantities can be calculated with the use of process-based numerical models. For this reason, understanding the physics behind the formation and the dynamics of gravitational mass movements is a topic of public concern and of permanent significance. Voellmy et al. (80) proposed for the first time a model for the calculation of the flow height and the velocity of snow avalanches. This milestone was achieved under the assumption of a dynamic block model, for which the gravity-driven force is compensated by the shear resistance, divided between a dry Coulomb-like friction term µ and a viscous-turbulent friction term ξ. With the advent of numerical computing, this pioneering approach has been extended to the solution of depth-averaged Eulerian models based on the incompressible equations of Saint-Venant, allowing simulations in complex terrain (10,47,60,89). These models allow to take into account multiple release areas, spatially variable friction parameters, and to simulate entrainment and deposition along the flow path. Due to their versatility and computational efficiency, these depth-averaged models represent a practical tool for researcher and practitioners for the simulation of mass movements including debris flow, snow, ice, and rock avalanches and the related hazards.
Process-based models attempt to reproduce real phenomena by capturing the fundamental physics governing them. At the same time, their level of detail must be balanced with ease of use and efficiency for practical applications. As a consequence, numerical models must adopt assumptions that limit the physical detail of the simulated processes. The common assumptions adopted by popular research and operational models are: homogeneous, continuous, and non-compressible materials; fixed bed or simple parametrizations of entrainment, shallow flow conditions with hydrostatic pressure distribution, constant vertical velocity profile, and depth-averaged equations. Therefore, state-of-the-art operational models do not account for three dimensional effects, complex erosion and deposition waves, dilatancy, interaction between multiple phases, phase and flow-regime transitions, fractures and collisions, segregation, and density variations. For the same reason, complex processes are often taken into account with the use of simplistic formulations that require the determination of parameters for which a clear physical interpretation is poorly constrained by observations if not completely absent. The inevitable loss in physical detail of the described processes leads to the nonphysical character of model parameters which thus cannot be evaluated based on state-of-the art geomechanical and rheological experiments. This shortcoming is classically compensated by statistical calibrations on the basis of field observations and historical records (2,10,53). In recent years, large efforts have been made in the direction of implementing more realistic two-and three-phase physical models, but these come with more parameters, at times even harder to calibrate (15,48,54,59). These unwanted but inevitable drawbacks of the trade-off between physical accuracy and efficiency is not only detrimental for the definition of the hazard scenarios and the related parameter values, which fundamentally limit the predictive capacity of the models, but it also hinders the interpretation of results, thus bearing important consequences for hazard management.
In comparison, three-dimensional (3D) simulations can fully resolve flow variations in all directions. Resolving greater physical details obviously comes at the cost of a longer computational time. Yet, large-scale 3D simulations have become possible with the increase in computational power. In recent years, particle-based continuum methods, including the smoothed particle hydrodynamics (SPH) method, the particle finite element method (PFEM) and the material point method (MPM), have gained increasing popularity in modelling geophysical mass flows, as they are able to easily handle large deformations and discontinuities. For instance, landslidetsunamis were extensively simulated using SPH (75,85), PFEM (23,64), and MPM (83,87). The dynamics of snow avalanches were simulated using SPH (1) and with MPM (42,43,72). Furthermore, the MPM has been used to simulate both the release mechanisms and flow of unsaturated and saturated (5,71,86) slopes as well as slab avalanche processes which involve mixed-mode brittle fracture propagation (27), anticrack nucleation (7,63) and dynamics (28). Fluid simulations with MPM benefit from recent advances related to improved grid-to-particle transfer schemes such as the Affine Particle in Cell (APIC) method which preserves both linear and angular momentum (36) In contrast to most depth-averaged approaches, these particlebased methods rely on well established theories and constitutive frameworks such as elastoplasticity and critical state soil mechanics. As such, model parameters represent physically interpretable quantities, which can be estimated on the basis of experimental tests (from the field or the laboratory) and from extensive literature studies. The detailed description of mass wasting processes in combination with physically constrained parameters creates the conditions for complex predictive simulations.
Here, we propose to advance our predictive capabilities of multi-phase gravitational mass movements and process cascades in complex alpine terrain through the development and the application of an innovative particle-based numerical modelling framework. First, we describe the numerical scheme and the constitutive relations used to simulate rock, debris, ice, snow and water. Second, in order to prove the potential of the modelling approach and to test its limitations, we perform four benchmark simulations based on well documented historical cases. The events are the 1963 Vajont rockslide and subsequent tsunami wave, the 2019 Piz Cengalo rock and ice avalanche, the 2019 Flüela Wisshorn rock and snow avalanche, and the 2020 Whymper ice and snow avalanche resulting from the hanging glacier collapse. These study cases cover a broad range of processes, sizes, and materials. We show that our approach is capable of explicitly simulating the mechanical behaviour of unstable masses and to reproduce the flow, erosion, entrainment and deposition of the simulated events. In addition, complex processes such as formation of lateral levees, rolling and erosion-deposition waves, avalanche fingering and block fragmentation naturally emerge from the model and provide relevant insights in the dynamics of the simulated process cascades. A critical evaluation of the results is done through quantitative and qualitative comparisons against field observations and previous numerical simulations. Finally, we discuss the potential of the model and its limitations and make suggestions for future research and practical applications.
A full 3D material point method with finite strain elastoplasticity The material point method (MPM) is a numerical technique used in computational mechanics to simulate continuous materials undergoing large deformations, topological changes, fractures and collisions. The underlying equations are the conservation of mass and momentum: where ρ is the density, t is the time, v is velocity, g is the gravitational acceleration and σ is the Cauchy stress tensor. Symmetry of the Cauchy stress tensor ensures conservation of angular momentum. Material-specific constitutive laws are required to close this system of conservation equations. The constitutive laws adopted in this study will be introduced in the next Section. In MPM, mass, momentum and the deformation gradient are tracked using material points (Lagrangian particles) which discretises the spatial extent of the material. In addition, a background Eulerian grid coupled to interpolation functions allows simple spatial differentiation of the stress field. Because there is no need for mesh connectivity, MPM does not suffer from the mesh distortion problem associated with the classical finite element method, and enables the description of fractures and collisions.
In this work, we rely on the finite strain elastoplastic MPM formulation proposed by Klár et al. (39). Figure 1 displays the computational procedure for a time integration step. First, the particles states, i.e. positions, masses, and velocities resulting from the previous time step, are transferred to the Eulerian grid (step 1) using a generalization of the Affine Particle-In-Cell method (APIC) method (36,37). This method allows for conservation of angular momentum in contrast to PIC (Particle-In-Cell) and FLIP methods (Fluid Implicit Particle-In-Cell) (9). Mass is conserved in MPM as the mass of each particle remains constant over time. Then, forces and updated velocities are calculated on the grid (steps 2 and 3). At this point, boundary conditions are applied on the grid velocities. A trial elastic deformation gradient is computed based on Eq. 3 (step 4) and the yield condition y ≤ 0 is tested to see if permanent deformation occurred (step 5). If such plastic deformation occurred in the trial step (y > 0), a return mapping  procedure (step 6) projects this non-admissible stress state back to the yield surface and the deformation gradient (step 7) is updated accordingly. If no plastic deformation occurred (y < 0), the trial gradient is conserved to the next steps. Finally, based on the updated deformation gradient and grid velocity, particle velocities are computed and particle positions are updated (step 8 and 9). A more exhaustive description of the numerical method is given in Jiang et al. (37) and Gaume et al. (27).

Multi-material constitutive relations and state equations
Alpine mass movements and process cascades involve multiple materials characterised by contrasting mechanical properties and different phases. In the following, we give details of the constitutive relations and the equation of state used here for their mechanical description. The mechanical behaviour of rock, debris, ice, and snow is described using finite strain elastoplasticity theory, making use of the common assumption of decomposing the deformation gradient F = F E F P into an elastic F E and plastic part F P . We further assume the Cauchy stress can be related to the elastic deformation through an isotropic hyperelastic relation, where J ≡ det(F ) and Ψ will in this work be chosen to be the Saint-Venant-Kirchhoff energy density function. Moreover, a yield criterion y(σ) ≤ 0 defines admissible stress states for an elastoplastic continuum and delimits the onset of plastic, irreversible deformations. Plastic deformations can thus only occur on the yield surface. In contrast, fluids are here described by an equation of state that relates the stress tensor to variations in pressure and density. Each modelled medium is described under the assumption of homogeneous mechanical properties.
In the following subsections, we present the yield criteria y(σ) used to describe rock, ice, saturated sediments, and snow, as well as details of the equation of state used to simulate water. All yield criteria are here expressed as a function of two stress invariants, the mean Kirchhoff stress p and the von Mises equivalent Kirchhoff stress q, defined respectively as where d is the spatial dimension of the problem, τ = Jσ is the Kirchhoff stress tensor and dev(τ ) = τ + pI is its deviatoric part. Positive pressure p > 0 corresponds to compression and negative pressure p < 0 corresponds to tension of the material.

Rock, ice and saturated debris
Geomechanical materials such as rock, ice, and debris are characterised by a pressure-dependent shear strength, a distinct asymmetry between tensile and compressive strength, and exhibit softening after failure. In order to model these mechanical characteristics in the simplest possible manner, we describe the yielding with a strain-softening Drucker-Prager model. According to the formulation proposed by Blatny et al. (7), the yield function is formulated as where γ = 2d sin ϕ/(3 − sin ϕ) can be related respectively to the material's friction angle ϕ, and γc represents the material's cohesion. The post-yield softening is controlled by the accumulated deviatoric plastic Hencky strain magnitude ε P S through an exponential function that allows us to adjust the brittleness of the material through the exponent χ, in particular, where ϕ ini and ϕ res are the initial and residual friction angles, respectively. The adopted nonassociative plastic flow rule is illustrated in Blatny et al. (7).

Snow
The mechanical properties of snow are modelled with the associative Modified Cam Clay developed in Gaume et al. (27). This model is able to capture the cohesive behaviour of snow, mixed-mode fracture, and compaction with a yield surface given by where p 0 is the pre-consolidation pressure, M is the slope of the cohesionless critical state line, and β the ratio between tensile and compressive strength. With increasing accumulation of volumetric plastic strain ε P V , the pre-consolidation pressure increases and the material hardens according to where ξ is the hardening factor and K the bulk modulus related to the Young's modulus and Poisson's ratio of the material. In the opposite case, the pre-consolidation pressure decreases, resulting in softening of the material and allowing fracture. After yielding, the model allows to describe the material's dynamics through a reduction of the slope of the critical state line to M s .

Water
Water flow is simulated with the use of the Tait equation of state, which combines computational efficiency with low density variations (6). The state equation relates the hydrostatic stress tensor to pressure, which in turn depends on density oscillations: where k and B are two constants. The density variations η allowed by the model can be approximated by the ratio between the flow speed v f and the speed of sound in the fluid c s : To limit the density variations to a maximum of 1 , we chose B = 3 × 10 8 Pa and k = 7 (6). The method can naturally handle effects such as splashing and breaking waves and is therefore well suited for the simulation of tsunami waves. This model has been recently validated by Wolper et al. (83) for the case of glacier calving and related tsunami propagation by comparison to laboratory experiments.

Benchmark simulations from historical events
For the purpose of demonstrating the potential of the model, we simulate four real-scale events characterised by complex process cascades. The selected study cases are all well documented events, whose volume ranges from tens of thousands to hundreds of millions of cubic meters. This selection includes various materials with contrasting mechanical properties such as: rock, saturated debris, glacier-ice, snow, and water. In chronological order (and by chance also by decreasing volume), the events are: the 1963 Vajont rockslide and subsequent tsunami wave, the 2017 Piz Cengalo rock avalanche and ensuing debris flow, the 2019 Flüela Wisshorn rock and snow avalanche, and the 2020 Whymper hanging-glacier collapse. For each event, we briefly introduce the historical and geographical context, then we present the specific model set up used for its simulation. Finally, we analyse results quantitatively and qualitatively, compare them against field observations and previous numerical studies and discuss the relevance of our modelling approach. Additionally, for each event, a movie of the simulation is presented in the supplementary material.  (23). In the following, we show how our framework allows to advance our modelling capabilities of events such as the Vajont rockslide by bringing together the points of strength of these previous studies, i.e. by including both landslide mechanics and its interaction with the reservoir lake with the finest resolution reached so far, and by utilising the most realistic geometry available.

Model set-up
We simulate the internal deformation of the rock mass on the slopes of Monte Toc, its consequent impact with the reservoir lake, and the dynamics of the tsunami waves in the valley and beyond the retention wall. The geometry is provided by the accurate reconstructions of Vacondio et al. (75) on the basis of previous work from Semenza, (69) and Superchi et al. (73). The  (8,19,87) and are given in Tab. 1.

Results and discussion
In the simulation, the rockslide terminates its downslope movement after 30 seconds, with a maximum velocity v max = 48 m s −1 . This result is in good agreement with historical studies and recent modelling investigations (87,88). In Fig. 2a, we show that the simulated displacement rates at point P 0 (as visible in Fig. 2b) is comparable to that obtained based on 2D MPM simulations by Yerro et al. (87) and in the range proposed by the literature (v P 0 max = 20 m s −1 ). In addition, our results closely relate to the ones of a recent study of the rockslide kinematics (88).
The quality of our simulation of the rockslide can be further evaluated by comparing the geometry of the deposited mass against field measurements and other numerical studies (Fig. 2b). Vacondio et al. (75) imposed a predefined rotational movement of the rock mass and did not simulate its mechanical behaviour. In contrast, in our study we account for the mechanics of the rockslide and its interaction with the water body. In general, the two modelling approaches provide a good fit with the field observations. Notably, we obtain the best results for section A, where the assumption of a rigid body appears to be less valid. For the other sections the results from the two studies are comparable, with a better performance of Vacondio et al. (75) close to the right orographic side of the valley. In fact, in this zone of the landslide, other studies have back-calculated a higher value of e.g. the basal friction angle (88) than the one used here and by Yerro et al. (87).
In Fig. 2c, we compare the discharge hydrograph from our study with the ones obtained by Vacondio et al. (75) and Franci et al. (22). The total value of the flooded water (i.e. the integral of the discharge function), the peak discharge and its timing and the trend of the curves are similar. However, with our approach, we obtain three discharge peaks, obtained only in pairs by the two other studies. Moreover, our modeled discharge exceed the other two in the first 20 to 40 seconds, whereas we obtain a lower discharge after the second peak until the end of the simulations. No discharge measurements are available to validate this result, but the ensemble of the modelling approaches provides a consistent input for further hydrodynamical studies about the propagation of the flood wave in the Piave valley. The quality of the simulation can be additionally evaluated on the basis of observations of the impact of the tsunami wave on land.
The maximum wave extent observed by Semenza (69) can be directly compared against the four numerical studies as shown in Fig. 2d. In the westernmost part of the basin, close to the village of Casso and where the wave over-topped the dam wall, our model shows the best results in good accordance with those of Crosta et al. (13). In addition, in our simulation the water reached a height of 945 m a.s.l. more than two hundred meters above the initial level of the lake, in accordance with historical evidence which reported damages to the houses of Casso up to 950 m a.s.l. (79). In the central part, at the edge of the area impacted directly by  . In a matter of seconds, the rock mass violently impacted the underlying Vadrec dal Cengal Glacier, and generated a jet stream of ice and water particles that hit the rocky Bügeleisen ridge line and thus got redirected towards the North (17 seconds in the supplementary movie). The impact of the rock mass on the glacier led to a loss of 0.6×10 6 m 3 of ice due to a combination of melting and erosion. In some areas, the small glacier was eroded throughout its thickness down to the bedrock (78). The mixed rock and ice avalanche continued its downslope flow, reached the steep glacier forefield and flew above a bedrock outcrop before impacting the sediment laden Clavera della Bondasca, where it decreased its velocity and finally came to rest after about 60 s. During its entire movement, the average velocity was estimated to 40-60 m s −1 . The reconstruction of the duration of the rock avalanche, its average velocity and even the reconstruction of its force history were possible due to inversion of low-and high-frequency seismic signals, measured at nearby stations (81). The entrainment of large volumes of glacier ice and wet debris (above 1×10 6 m 3 ) did not lead to a particularly large runout distance of the rock mass. The rock avalanche deposited in correspondence of the turn of the river, not far away form the deposition of the smaller 2011 event (volume of about 1.5×10 6 m 3 , see Fig. 3). However, the combination of entrained glacier ice and the ejection of water from wet sediments enhanced mobilisation of the avalanche deposit (81). Less than a minute after the deposition of the rock avalanche, a debris flow with an initial volume of about 0.5×10 6 m 3 initiated in absence of any precipitation, followed by twelve similar events within the following week and finally by two more events triggered by rainfall on August 31.

Model set-up
We simulate the first part of the process cascade including the rock mass detachment from the north-east face of the Piz Cengalo, the entrainment of glacier ice from Vadrec dal Cengal, the entrainment of debris in glacier forefield and in the Clavera della Bondasca and finally the deposition of the rock-ice avalanche. The volume of the unstable rock and the entrained ice are obtained by difference of the digital elevation models from 2018 and 2012 (data from Swisstopo and Kanton Graubünden). This methodology, valid for the calculation of the rock mass, also applies to the glacier due to the fact that this was completely eroded down to its bed (78). The sediments are assumed to be of homogeneous thickness (H sed = 10 m) over the entire area and the erosion and entrainment process are explicitly simulated. The resolution of the simulation is set to 1.5 m, resulting in a total of 16×10 6 material points. Rock, ice and debris are simulated with a Drucker-Prager model with different material properties. The phase change of glacier ice after the impact of the rock mass is simplified by imposing a lower value of the residual friction angle after the yield surface is reached. For details about the parameters refer to the Methods and Tab. 1.

Results and discussion
The unstable volume in the Piz Cengalo north-east face is the initial condition that originates the process cascade. Immediately, the mass starts displacing, causing large deformations and fractures that lead to a general collapse of the rock (Fig. 3a). Here, the fragmentation is solely controlled by the internal stress distribution in the homogeneous mass, whereas preconditioning weak surfaces are neglected. In a matter of seconds, the mass violently reaches the glacier, where the impact is explicitly simulated by taking into account the mechanical properties of the two materials, rock and ice. Most of the glacier mass is eroded, in good agreement with field observations. The tremendous amount of energy causes a part of the glacier to behave in a fluid manner and create a propagation wave that moves towards the north-east (according to the slope of the rock face) and hits the opposite rocky ridge at about t = 20 s. Directly after, the mass hits the sediments in the glacier forefield and creates a vast vertical movement spreading over the entire valley up to the opposite slopes towards the east, in agreement with video recordings of the event.
When the explosive energy that caused the mass to disperse in the air is dissipated, the rest of the moving mass is developing active dynamics. On the right orographic side, the avalanche flowed above the bedrock outcrop, strongly eroding the sediments and quickly reaching the Clavera della Bondasca. On the left orographic side, a second avalanche front follows and fills the valley with a delay of a few hundred of meters to the other front caused by the gentler gradient of the slope in this section. The avalanche dynamics is dominated by surging waves with strong three-dimensional effects. The rock mass, although strongly fragmented, still conserves large boulders, transported downslope and deposited mostly at the tail of the moving mass. The result of the deposition process is visible in the comparison of the aerial photographs of the bedrock outcrop and the Val Bondasca with our modelling. A large amount of boulders is clearly visible after the event as well as in our simulation. Those are, both in reality and in the numerical results, concentrated in an area above and below the bedrock outcrop and mostly on the left orographic side of the valley.  Shortly after one minute of flow, most of the mass is completely at rest, but some surges still develop close to the surface of the fresh deposit between the glacial terrace and the lower part of the valley. This material is mostly constituted by a mixture of glacier ice with some rock and sediments. This mechanism resembles Scenario 2 as described in the numerical study by Mergili et al. (46) for explaining the onset of the debris flow. From the results of our simulation, we suggest that this effect is naturally driven by the multi-phase behaviour of the granular avalanche and that it likely occurred in reality. However in this study, we did not consider poromechanical effects (74) which prevented us from simulating the ejection of water from the sediments and thus debris flow initiation and dynamics. Therefore, we can not conclude if this scenario would be sufficient to trigger the observed debris flow activity nor if more water from the sediments is required. Instead, we suggest that a combination of the two hypothesis is the most probable explanation of the real event.
Succeeding the surge phase, the mass comes to an almost complete rest after about 250 s. In Figure 3a, we show observations before and after the event as well as our modelling results in three areas in Val Bondasca. The three areas are identified in panel b, where a qualitative comparison of the final deposit and the simulation at t = 40s is shown. The final deposit in Val Bondasca is analysed in more detail in the right panel of Fig. 3a, where our modelling results are qualitatively compared against aerial imagery of the site before and after the event (or during the collapse of the rock mass). The 2017 events deposit not further than a kilometer downslope of the 2011 event, but this difference is enough for the rock avalanche to reach the curve of the river and drastically turn westwards. This effect, controlled by the topography of the valley, causes the sharp deposition pattern clearly visible in the photograph as well as in the rendering of the modelling results at the southern edge of the transition between the avalanche deposit and the pre-existing sediments in the river bed. This peculiar pattern can be considered the best indicator of the run-out distance of the process cascade excluding the debris flow. As evident in the visual comparison, the model is able to reproduce both the run-out distance and the geometrical pattern with great accuracy. Signs of the impact of the avalanche and significant deposit volumes are visible further up-slope, shortly before the bedrock outcrop, on the right orographic side of the eastern little ice age moraine. On the north-face of the Piz Cengalo, the simulated dynamics of the rock avalanche is overestimated. This effect is explainable by a mismatch between the modelled mass flowing towards this direction and the real event: the overestimation is due to aggregation of small rock fall events that in reality happened at different times. This part therefore should be interpreted with care, but it is essentially negligible when looking at the entire rock avalanche and process chain.
A more quantitative analysis is possible on the basis of the seismic data presented in Walter et al. (81) and displayed here in figure 3c. Here, the four gray bars indicate, from left to right, the initiation of the avalanche, (A) the impact of the rock mass with the glacier, (B) the impact of the avalanche with the Bügeleisen ridgeline, and (C) the slow-down of the mass in the Clavera della Bondasca. The upper plots show the observed relative force components, the speed and the seismicity based on the measurements from the Swiss Seismological Service. The four lower graphs show the results of the modelling. We plot the total kinetic energy for the rock (orange), ice (blue), and sediment (green) masses alone and the total value (in pink) against time. The rock mass displays two peaks in energy, the first of which corresponds to the initial falling movement from the Piz Cengalo north-east face. Shortly after, the rock mass is redirected towards the west (see East-West component of the force history) and together with the entrained mass of the ice glacier develops a second peak. This maximum corresponds to a prominent peak in the ice mass, related to the glacier impact and, together with the rock, related to the flow above the steep glacier forefield and the bed outcrop. As the initial momentum gets dissipated over time through friction, a large mass of sediments is entrained and increases its kinetic energy. This component is dominating the total signal during the final deposition phase.

The 2019 Flüela Wisshorn rock-and snow avalanche Geographical and geological context
The Flüela Wisshorn is a 3085 m high mountain peak located 2 km northwest of the Pass dal Flüela at the border between the two municipalities of Davos and Zernez (CH). The main summit is divided from the south-west peak by an almost flat ridge constituted by Augengneiss and characterised by permafrost conditions. On March 19 2019 shortly after midnight, a rock fall of more than 250×10 3 m 3 detached from the west flank of the mountain (38). The collapse and the downhill movement of the rock mass released a large snow slab, which developed in a snow avalanche and altogether entrained large amounts of snow from the slopes. Most of the rock mass deposited in the combe below the rock wall within a kilometer from the detachment, while the snow avalanche mixed with some debris almost reached the road of the Flüelapass and came to rest in the vicinity of the construction of the Wägerhütta, as visible in Fig. 4a. A considerable amount of rock deposited in correspondence of a rock glacier, which in the next years is expected to accelerate in its creep due to an increase in the driving stresses (12). Seismic measurements from the nearby stations of Jenatschalp in Dischma valley indicate that the rock avalanche lasted not more than one minute. At the time of the event, according to field measurements and observations in nearby snow stations, the snowpack on the mountain flank was more than two meters thick with maximum values above three meters.

Model set-up
We simulate the rock mass detachment and the subsequent process cascade, including snowpack mechanics and rock-and snow avalanche dynamics. The volume of the rock detachment is directly obtained by subtracting two digital elevation models (DEMs) obtained prior and after the collapse. DEMs were obtained by means of airborne laserscan and by uncrewed aerial vehicles surveys. The high resolution terrain model used for the flow area has been obtained in summer, thus it is processed with a Gaussian filter with a radius of 10 meters in order to account for the smoothing effect of the snow cover. The rock mass is modelled with a Drucker-Prager model, the snow mass with the Modified Cam Clay model (27). For details see the Methods and Tab. 1.  The mechanical behaviour of the snow cover is explicitly modelled for an area sufficient to completely cover the terrain interested by the mass flow. The height of the snow cover is set to H snow = 4 m, slightly exceeding the observed snowpack thickness. This value is overestimated on purpose, so that the dynamics of the rock and the snow avalanche are de-facto directly influenced by mechanical interaction with the snow cover and neither come in direct contact with the no-slip boundary at the base of the snowpack. The stability of the snowpack prior to the rock mass failure is guaranteed by a value of M = 1.2, corresponding to a realistic equivalent friction angle of ϕ = 30 • (76). Outside this area, a confining snowpack is imposed as boundary condition in order to recreate realistic stress conditions for the simulation (see Fig. 4). The resolution of the model is set to 1 m, resulting in 17.5×10 6 material points.

Results and discussion
The unstable rock mass is the trigger of the process cascade. Fractures propagate in the fragmenting rock mass, which begins its downward motion as soon as the simulation starts. As a consequence, snow is entrained mainly at the front and at the bottom of the rock avalanche. Complex erosion and deposition patterns can be observed in the modelling, including secondary slab detachments, roll waves and fingering patterns. Lateral levees form in the middle and in the lowermost part of the avalanche, where the slope angle is low. Although we observe secondary snow slab detachments in the model, the extent of the observed slab at the release area is not reproduced correctly for our set of parameters. A possible explanation of this mismatch is that, in our simulation, the mechanical properties of the snow cover are assumed in first approximation homogeneous over the entire catchment. A more detailed parametrization, with elevation and aspect dependency of the snow mechanical properties, would likely allow to realistically reproduce the spatially variable snow cover properties and its stability, thus the correct extent of the slab detachment. These processes emerge naturally from the modelling, due to the combination of the complete 3D description of mass and momentum conservation and realistic constitutive relations (see the rendering in Fig. 4 and the movies in the supplementary material).
Our approach allows the explicit simulation of the rock fragmentation of the unstable rock mass during the detachment from the rock wall and along the flow path. In the deposition area, several boulders are discernible in the simulation results, coherent with the observed deposit, where the largest boulder is up to 20 m length along the major axis. Most of the mobilised mass deposits below the first steep slope in the terrain depression between 2600 and 2800 m a.s.l. after about 40 seconds, in good agreement with the seismic observations. Following the process cascade, part of the snow and a smaller fraction of the rock mass have enough energy to overcome the overdeepening and re-gain speed in their motion towards the valley bottom and the pass road. The mixed mass flow finally comes at rest in correspondence of gentler slope, 110 seconds after its start. Unfortunately, no other data are available for discussion of the avalanche dynamics, but more can be analysed regarding the consequences of the event in the environment. In Figure 4a, we directly compare an image of the Flüela Wisshorn after the event with a rendering of our simulation. Visually, a general agreement is apparent in terms of run-out distance, shape of the avalanche and density of the rock mass. The deposition pattern can be analysed in detail with aid of digital elevation models prior and after the event and comparing our results with photographic material obtained a few days after the avalanche. Figure 4b shows the deposition pattern along the longitudinal profile of the avalanche, starting at the top of the detachment area and terminating at the furthest position in the deposition zone. The difference between the initial (dashed line) and the final height (continuous-dotted line) of the rock and snow masses gives insights into the erosion and deposition patterns. For a distance of almost 200 m from the origin, most of the snow cover has been eroded at the front of the rock avalanche, which in turn deposited several meters of rocks from the tail of the moving rock mass. This effect displays no clear signal in the observed elevation difference, possibly due to a net close-to-zero balance, but it is clearly confirmed from the optical imagery and is reproduced in the modelling. After about 400 m of distance, the deposited mass increases drastically up to more than 10 meters of height. Further down slope, at about 800 m along the profile, the deposit reaches a maximum value of about 15 meters. The deposition patterns are well reproduced in our modelling in terms of both their longitudinal distribution and magnitude.
The maps shown in Figure 4c allow extending the analysis from a linear to a planar perspective. In the upper part, it is apparent that the general shape of the deposit, the spatial location and extent of its maximum values, as well as secondary patterns, such as the lobe-shaped deposition in correspondence of the rock glacier front, are closely reproduced. In the lower part, below the rock glacier, no elevation signal can be recognised in the observations. Comparing this with our results might indicate an overestimation of the modelled deposit height. This effect can be explained by a possible overestimation of the entrained mass in this zone because an increase in snow cohesion and friction with decreasing elevation was not accounted for. However, no conclusion can be drawn on this point, because the modelled values are still in the uncertainty range of the measurements (∼ 2 m). Also, despite no clear signal in the elevation data, the optical imagery clearly shows the extent and provides a general indication about the density of the rock mass in the lower part, which is in accordance with our results.
The 2020 collapse of the Whymper hanging glacier

History of instabilities
Below the summit of the Grand Jorasses on the south side of the Mont Blanc massif, at an elevation of 4000 m a.s.l., is located the Whymper hanging glacier. This unbalanced cold ramp glacier has become famous in the scientific community as well as in the media due to its history of recurrent instabilities, periodically threatening the village of Planpincieux and its surroundings in Val Ferret (IT) (44). Several glacier collapses have been documented with volumes ranging from few tens of thousands to several hundreds of thousands cubic meters (45). The terrain below the hanging glacier is steep and strongly crevassed and divides in four possible tracks, also reaching the Planpincieux glacier.
Here, we investigate the small event of November 2020, which despite its limited volume , is well constrained. The volume, separated from the rest of the glacier by a deep crevasse, detached from the left orographic side of the cold based glacier (see Fig. 5a). The rupture took place above the rock/ice interface, similarly to previous instabilities (18). For the investigated event, no data are available about the motion of the avalanche (duration, velocity), but some optical images taken after the event allow to reconstruct the run-out distance in the two main flow tracks (Whymper and Planpincieux) and in general the deposition pattern of this event. In-situ surveys in the unstable area also provide the digital elevation models necessary for our simulation. Regarding the snow cover properties, according to observations at nearby monitoring sites and numerical simulations, at the time of the event the snowpack was around 1 meter thick in the lower part and not exceeding 2 meters in the upper part with dense snow with rounded grains (Aosta Avalanche Service, personal communication). It is important to note that in the month of October of the same year, a few weeks prior to our study case, another volume of ice of similar dimension but slightly different position broke off from the hanging glacier. As a consequence, the large crevasses on the Planpincieux track have been filled with snow and ice, effectively smoothing the glaciated terrain.

Model set-up
We simulate the detachment of the ice mass, the mechanics of the snowpack and the dynamics of the resulting mixed-avalanche. High resolution digital terrain models are available prior and after the event in a small region around the hanging glacier, which are used to reconstruct with precision the unstable volume. The terrain of the entire valley, used for the computation of the avalanche dynamics, has been obtained by the local public authorities by means of airborne laser scanning in summer 2008, twelve years before the event. This time shift might cause substantial differences in the glaciated terrain, with possible elevation changes up to several meters and variations in the crevasses patterns. The ice mass is modelled with a Drucker-Prager model, the snow mass with the Modified Cam Clay model as for in the Flüela Wisshorn case.
For details see the Methods and Tab. 1. In order to realistically initialise the snow cover in this extremely steep terrain, we impose 2 m of snow cover only on areas with a slope angle gentler than 50°. The mechanics of the snowpack is modelled over a large area (blue area in Fig. 5) around the flow path; initial and boundary conditions are set in the same manner as in the case of the Flüela Wisshorn. To account for the filling effect of the October collapse, the strongly crevassed glacier terrain on the Planpincieux track has been smoothed with a Gaussian filter with a radius of 20 meters. Simulations without this smoothing are also performed for the sake of comparison.

Results and discussion
The simulation of the process cascade begins with the break off of the ice volume at the Whymper hanging glacier (see Fig. 5 and Supplementary movies). In Figure 5b, the evolution of the front velocity along the two different paths can be analysed against the horizontal displace-  ments. The unstable mass quickly hits the underlying hanging glacier and reaches a maximum in velocity (similar values are obtained for all the lines in the graph) with values of more than 30 m s −1 . After this high-altitude plateau, half of the mass flows straight into the steep Whymper couloir, losing contact with the ground due to the high kinetic energy and the abrupt change in slope of the terrain at the border of the hanging glacier. After several tens of meters of free falling in the air, the mass re-gains contact with the terrain and flows until the end of the rock ridges where it peaks at about 50 m s −1 in velocity before the terrain widens on a moderately steep terrace (see Fig. 5). Here, the avalanche slowly decreases its velocity before coming to rest, at a position in good agreement with the run-out distance reconstructed on the basis of webcam imagery and expert opinion (blue bar in Fig. 5c). The rest of the mass flows towards the Planpincieux glacier, on relatively gentler and deeply crevassed glacier terrain. As explained in the model set-up, two scenarios have been simulated for this branch: one considering smoothed terrain (due to a previous event -dark green line), the other one simply using the available terrain data (light green line). Both scenarios reach another velocity maximum at the divergence between the Planpincieux and Jorasses glaciers, before rapidly decreasing their kinetic energy after about one kilometer of distance from the detachment, where they create multiple fingering patterns. Despite a larger uncertainty in the definition of the deposition and the run-out distance in this area, our modelling results are again in good agreement with the observations.

Discussion
The modelling approach we propose has the ability to accurately simulate complex mass movements and process cascades on the basis of solid physical and mathematical formulations. This important achievement is made possible through the coupling of a three-dimensional modelling framework based on finite-strain elastoplasticity with advanced constitutive models for different geo-materials. Interactions between different materials is naturally handled with MPM. The results of the four benchmark simulations presented in the previous section demonstrate the potential of this model and illustrate its practical applicability to real-world events. The major advantages and limitations of our approach can be discussed in the following points: • The adopted mathematical description of the constitutive laws of different materials is deeply grounded in more than a century of research in solid and geo-mechanics. This has important implications for the simulation of alpine mass movements and process cascades. In particular, in contrast to the parameters of the Voellmy model used in most depth-averaged modelling approaches (10,48), parameters of the constitutive laws used in this study have a true physical meaning and can be directly obtained on the basis of laboratory experiments (e.g. triaxial test). These aspects confer a strong predictive capability to our physically based model, bringing us one step forward in the simulation of unprecedented events and process cascades and in the management of natural hazards.
• Another important aspect is the explicit description of entire complex systems in three dimensions, including the mechanical behaviour of eroded material on entire slopes on complex terrain. The solution of the entire system of equations constituted by conservation of mass and momentum for compressible materials allows to resolve vertical stresses and velocity profiles under different loading conditions. This level of detail is necessary for the simulation of complex three dimensional effects such as rotational movements, complex trajectories non adherent to the ground, dilatancy, compaction and complex interaction processes. Additionally, the use of elastoplasticity enables the simultaneous simulation of meta-stable and highly dynamic systems (e.g. initially stable snowpack disturbed by an avalanche). The latter favours an unprecedented description of erosion, entrainment, and deposition, which in turn results in the natural emergence in the model of complex phenomena that control interaction processes between different materials and strongly influences the dynamics of process cascades. Some notable examples that have been achieved in our four benchmark simulations are breaking waves in water, lateral levees, fingering patterns and secondary avalanche releases in snow, fracture propagation and fragmentation in rock and ice, ballistic trajectories of individual rock boulders and erosion/deposition waves.
• The mathematical foundations of the proposed model are not new (36,39,72). Yet, its application to complex full-scale real-world events has been hampered for decades due to high computational costs. Here, by taking advantage of parallel computing and highend CPUs, we were able to simulate mass movements lasting up to 200 seconds using a number of material points ranging from 10 to 20 millions, in less than 70 hours. In details, we used a personal workstation, equipped with 36 Intel(R) Core(TM) i9-10980XE 3.00GHz CPUs and 128 Gigabytes of RAM. In order to reduce the computational cost of our simulations, which depends on the elastic wave propagation speed, we artificially reduced the elastic modulus of rock and ice by a factor of 100 (see Tab. 1). While this would influence simulations of the release mechanisms, we verified that this assumption does not affect the mass movement dynamics, which is dominated by friction and cohesion. In the future, the computational capabilities can be further enhanced by relying on distributed memory parallelism, including multi-GPU simulations (82).
• Concerning current limitations of the model, it should be noted that thermal effects were not taken into account. In reality, mechanical properties of geo-materials, especially cohesion and friction, strongly depend on temperature. This is especially true for large mass movements in alpine settings, where water and ice play an important role in their dynamics (20,65). Mechanical properties can not only change due to elevation and thus initial temperature differences but also due to the mechanical work created during flow through friction and finally transformed into heat. In our simulations, we have presented a highly simplified approach to simulate phase change for glacier ice, but this was only valid under the assumption of immediate transformation given the huge forces exerted by the rock avalanche. In the future, the energy balance equation should be included together with empirically-based equations relating materials' friction and cohesion to tem-perature (14). An additional limitation is the rate independence of the plasticity model. The next step will thus be to include viscoplasticity, e.g., through Perzyna's overstress model (55), combined will well-known rheological models for granular materials (49). This addition would allow simulations of hyperconcentrated suspensions and viscous creep processes such as debris flows, glacier dynamics, lava streams, and permafrost creep (11,24,34,58,67). Finally, we should emphasise that our approach is still computationally very expensive and may not (at present state) be suitable for processes involving slow motion over very long distances such as debris or lava flows. In that case, we would suggest to use our approach on smaller scales in order to better understand mechanical processes at play and then inform more efficient depth-averaged models for the simulation of very long run-outs.

Conclusions and future perspectives
The risk from alpine mass movements and process cascades is a relevant problem in all mountain ranges across our planet, currently increasing because of both varying environmental factors affecting slope stability and steadily growing anthropic pressure in vulnerable regions. The hazard related to mass movements of unprecedented magnitude and process cascades cannot be assessed neither effectively mitigated by means of classical numerical models, but requires more sophisticated approaches. We have presented here a novel three-dimensional multi-phase model for the detailed simulation of complex alpine mass movements and process cascades. The proposed approach is based on finite strain elastoplasticity and the Material Point Method, a hybrid Eulerian-Lagrangian numerical technique used to solve mass and momentum conservation equations. We show the potential of our approach by performing four benchmark simulations of well documented historical events. Complex phenomena such as the formation of leveed channels, erosion-deposition waves, fragmentation, ballistic trajectories, surges and fingering patterns naturally emerge from the model and allow a detailed description of the simulated dynamics. The materials are described with constitutive laws whose parameters have a clear physical meaning and can be directly obtained on the basis of experimental studies, facilitating the parametrization of the model. As a consequence, our approach possesses a strong predictive character, which makes it suited for detailed investigations and forecast of unprecedented and complex events. Given the current challenges in the simulation of gravitational mass movements, future research should focus on i) further development of advanced numerical models which include fully thermo-hydro-mechanical coupling and rate-dependency ii) laboratory and field investigations of the rheological properties of the materials and their dependency on environmental conditions iii) process-based investigations of the influence of climate change on the magnitude and dynamics of large events and process cascades.