Brief communication: Depth-averaging of 3D depth-resolved MPM simulation results of geophysical flows for GIS visualization

Significant advances in full-3D modeling of geophysical flows have provided deeper insights into complex processes and predictive potential. However, practical application in the natural hazard community remains limited due to inadequate GIS integration of simulation results. This study addresses the oxymoronic transformation of 3D depth-resolved MPM simulation outputs into simplified depth-averaged results, such as flow depth and thickness, and slope-parallel and slope-normal velocities. Specifically, we present an algorithm that rasterizes scattered MPM outputs into a 2D format, enhancing their utility for hazard mapping and mitigation. We demonstrate our approach by applying it to an ice avalanche event, which is simulated using MPM and visualized in GIS. Notably, the 3D MPM shows slope-normal flow velocity components over terrain jumps, and our algorithm enables the identification of flow detachment from the terrain, which depth-averaged models typically neglect.

1 Introduction

With the emergence of computers, numerical modelling of geophysical flows (e.g., debris flows, rock, snow, and ice avalanches) has become an integrated part of natural hazards mapping and mitigation. Starting from the 1960s, depth-averaged models that simulate gravitational mass movements along a prescribed 2D topography were developed (see, e.g., the review by Eglit et al., 2020 for snow avalanche models; Hungr, 1995 for debris flows). Later, depth-averaged models have been extended to 3D topography, also aided by higher spatial coverage and resolution of Digital Terrain Models (DTMs). Several depth-averaged models have been proposed to model rock avalanches (e.g., Mergili et al., 2012), debris flows (e.g., Denlinger and Iverson, 2004; McDougall and Hungr, 2004; Iverson and George, 2014; Tayyebi et al., 2021), and snow avalanches (e.g., Christen et al., 2010; Vicari and Issler, 2025). The extensive scientific advancement and adoption of depth-averaged models can largely be attributed to their computational efficiency, as these models modify the 3D conservation equations by integrating them in the vertical direction or perpendicular to the topography – thus reducing the number of conserved flow variables that need to be resolved – and assume simplified rheologies. Additionally, their straightforward integration into Geographic Information Systems (GIS) has also contributed to their popularity, particularly among natural hazard agencies and practitioners focused on hazard mapping.

Recently, advancements in computing power have enabled the resolution of 3D conservation equations within a feasible (yet extensive) timeframe, making it possible to simulate geophysical flows at the slope scale. 3D models have been applied for debris flows (e.g., Kwan et al., 2015; Koo et al., 2018), snow avalanches (e.g., Gaume et al., 2019; Li et al., 2021; Kyburz et al., 2024), rock and ice avalanches (e.g., Cicoira et al., 2022), and landslides (e.g., Franci et al., 2020; Tran et al., 2024). Nonetheless, although these models are scientifically popular, their use in practical applications for evaluating and mitigating natural hazards remains limited. We believe that this gap can be attributed to multiple factors: (i) inexperience within the natural hazards community and scarcity of guidelines; (ii) limited laboratory experiments to determine the constitutive parameters of geomaterials and challenges in determining these parameters in the field; (iii) non-user-friendly software; (iv) higher computational cost compared to depth-averaged models; (v) poor integration of simulation results with GIS tools. The latter issue stems from a mismatch between the output of these 3D depth-resolved simulations (i.e., where conservation equations are also resolved in the direction normal to the terrain) and the 2D raster-based formats used in GIS platforms. Indeed, in these simulations, physical quantities such as material positions and their velocities vary continuously along the terrain-normal direction and are expressed in absolute 3D coordinates. However, 2D raster data can represent only a single value per grid cell, requiring some form of averaging and reduction of the depth-resolved physical quantities to the local topography. This dimensional mismatch limits the straightforward use of 3D simulation results in hazard mapping. Within this work, we aim at addressing this last point by developing a tool to export the results of fully 3D Material Point Method (MPM) simulations into a rasterized format. Specifically, our aim is to convert this depth-resolved MPM particle information into depth-averaged variables that can be saved in a 2D raster format and visualized in GIS. Han et al. (2020) and Su et al. (2024) already proposed algorithms to extract physical quantities from 3D Smoothed-Particle Hydrodynamics (SPH) simulations. Nevertheless, considering the significance of MPM for modeling gravitational mass movements, there is a need for a dedicated algorithm; our method also introduces novel elements, especially concerning the exported depth-averaged variables. We illustrate the practical application of the depth-averaging algorithm by presenting an ice avalanche event simulated using MPM and subsequently exported to GIS. We then examine the depth-averaged variables derived from the 3D MPM simulation and discuss potential future enhancements of the exporting tool.

2 Methods

2.1 3D depth-resolved MPM model and material parameters

MPM is a hybrid Eulerian–Lagrangian method commonly used to simulate the behavior of granular and fluid materials within a continuum mechanics framework. In our case study, we use a Drucker–Prager yield criterion to model the constitutive behavior of ice. However, alternative material models may also be employed within the MPM model to simulate other materials, such as snow or soil (e.g., Gaume et al., 2019; Cicoira et al., 2022). The Drucker–Prager yield function is defined as:

$$\begin{array}{}\text{(1)}& q-\mathit{\gamma }(p+a)\le \mathrm{0},\end{array}$$

where q and p are the deviatoric and mean stress, respectively. The material parameters γ and a represent friction and attraction, respectively. These can be derived from the Mohr-Coulomb yield criterion parameters – i.e., the friction angle φ and the cohesion c – given that the two yield criteria are assumed to be equivalent under triaxial compression conditions: $\mathit{\gamma }=\left(\mathrm{6}\sin \mathit{\phi }\right)/(\mathrm{3}-\sin \mathit{\phi })$ and $a=c/\tan \mathit{\phi }$. We assume a non-associative flow rule with zero dilatancy. The model incorporates softening, similar to the model of Cicoira et al. (2022), wherein we assume that the friction angle decreases from an initial φ_i=40° to a residual φ_r=20°, and cohesion decreases from an initial c_i=1 kPa to 0, depending on the accumulated deviatoric plastic strain ${\mathit{\epsilon }}_{\mathrm{d}}^{\mathrm{p}}$. Unlike the simulations of Cicoira et al. (2022), we do not model explicitly erosion of glacier's ice and snow cover. Instead, we assume that the basal friction is higher when the flow is in contact with the glacier (μ=0.65) than with bedrock or soil (μ=0.55). This frictional contrast yielded a reasonable match to the observed inundation area. One possible explanation for this assumption is that the erosion of compact, wet snow on the glacier surface may lead to reduction in flow mobility, as also suggested by Li et al. (2022). At t=0, the ice bulk density is set equal to ${\mathit{\rho }}_{\mathrm{0}}={\mathit{\rho }}_{p\mathrm{0}}=\mathrm{850}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, where ρ_p denotes the material points' density. The simulation is carried out with a grid size of gs =1 m, with n_g=6 material points per cell. A complete list of the material parameters is reported in the Supplement (S2).

Figure 1 Discretization of the topography in ASCII format and calculation of flow-related physical variables: (a) global and local coordinate system; (b) calculation of the surface normals and mapping of the material points to their corresponding cell, with definition of flow velocity components, which are parallel (U_pt) and normal (U_pn) to the terrain; (c) Schematic of flow depth h, defined as the maximum distance from the topography among the material points within a cell; and flow thickness t, defined as the difference between the maximum and minimum distances from the topography among all material points within a cell. The two quantities may differ notably when the flow detaches from the terrain.

2.2 Topography representation

MPM simulations are performed on a generic topography Σ, whose elevation is defined in a global Euclidean system (ENH) as H(E,N) (Fig. 1). We denote the latitude and longitude coordinates of sampled points on the surface as (E_s,N_s), and H_s their respective elevations on Σ. Each sampled coordinate (E_s,N_s) specifies the center of a cell with index (i,j), which is separated by a distance cs (the cell size of the DTM) from the adjacent cell (Fig. 1a):

$$\begin{array}{}\text{(2a)}& i={\displaystyle \frac{E_{\mathrm{s}}-E_{\mathrm{llc}}}{\mathrm{cs}}},\end{array}$$

$$\begin{array}{}\text{(2b)}& j={\displaystyle \frac{-N_{\mathrm{s}}+N_{\mathrm{llc}}+\mathrm{cs}\cdot (r-\mathrm{1})}{\mathrm{cs}}},\end{array}$$

where (E_llc,N_llc) define the planar coordinates of the centre of the lower left cell of the DTM. The number of cells of the discretized topography in the east and north directions is indicated by c (i.e., columns) and r (i.e., rows), respectively. The parameters (E_llc,N_llc), $(\mathrm{cs}\cdot c,\mathrm{cs}\cdot r)$, cs, define the location, extension and resolution of the DTM, respectively.

To project the simulation results on the topography, it is necessary to define a terrain orientation metric. The unit vector normal to the topography (outward-pointing) is expressed as (see, Fischer et al., 2012):

$$\begin{array}{}\text{(3)}& {\mathit{n}}_s={\displaystyle \frac{(-{\partial }_{E}H,-{\partial }_{N}H,\mathrm{1}{)}^T}{\sqrt{\mathrm{1}+({\partial }_{E}H{)}^{\mathrm{2}}+({\partial }_{N}H{)}^{\mathrm{2}}}}}.\end{array}$$

In the exporting code, the spatial gradients of the surface height are approximated using central difference at each grid point:

$$\begin{array}{}\text{(4a)}& {\partial }_{E}H(i,j)\approx {\displaystyle \frac{H(i+\mathrm{1},j)-H(i-\mathrm{1},j)}{\mathrm{2}\cdot \mathrm{cs}}},\end{array}$$

$$\begin{array}{}\text{(4b)}& {\partial }_{N}H(i,j)\approx {\displaystyle \frac{H(i,j-\mathrm{1})-H(i,j+\mathrm{1})}{\mathrm{2}\cdot \mathrm{cs}}}.\end{array}$$

2.3 Depth-averaging of MPM results

The results of the MPM simulation at a certain time t_e include particles (p) with their position ${\mathit{X}}_p\left(t_e\right)=(X_p,Y_p,Z_p)$ and velocity ${\mathit{U}}_p\left(t_e\right)=(U_p,V_p,W_p)$, expressed in the global coordinate system ENH. The algorithm initially finds the indexes (i_m,j_m) correspondent to the DTM's cell 𝒞_m to which each particle p belongs (Fig. 1b):

$$\begin{array}{}\text{(5a)}& i_m\left(t_e\right)=\underset{i\in \mathbb{Z}\cap [\mathrm{0},c-\mathrm{1}]}{\mathrm{argmin}}\left|{\displaystyle \frac{X_p\left(t_e\right)-E_{\mathrm{llc}}-\mathrm{cs}\cdot i}{\mathrm{cs}}}\right|,\end{array}$$

$$\begin{array}{}\text{(5b)}& j_m\left(t_e\right)=\underset{j\in \mathbb{Z}\cap [\mathrm{0},r-\mathrm{1}]}{\mathrm{argmin}}\left|{\displaystyle \frac{-Y_p\left(t_e\right)+N_{\mathrm{llc}}+\mathrm{cs}\cdot (r-\mathrm{1})-\mathrm{cs}\cdot j}{\mathrm{cs}}}\right|.\end{array}$$

Unlike Han et al. (2020) and Su et al. (2024), which project particles along the terrain normal, our approach uses vertical projection to assign particles to raster cells (see discussion in Sect. 2.3.2). The code offers the possibility to export the depth-averaged results to a new digital model with a different cell size (cs_e) than the original DTM's cell size (cs). The indices (i_e,j_e) for each cell 𝒞_e in the new model are thus computed using equations similar to Eqs. (5a) and (5b), but with cs_e.

2.3.1 Flow heights

Once the DTM's cell 𝒞_m has been identified, the elevation H_p of the surface at the planar coordinate (X_p,Y_p) is calculated imposing the coplanarity condition

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathit{n}}_s^T(i_m,j_m)& \cdot \left[X_p\left(t_e\right)-E_{\mathrm{s}}\left(i_m\right),Y_p\left(t_e\right)-N_{\mathrm{s}}\left(j_m\right),\right.\\ & {\left.H_p\left(t_e\right)-H_{\mathrm{s}}(E_{\mathrm{s}},N_{\mathrm{s}})\right]}^T=\mathrm{0}.\end{array}\end{array}$$

For each non-empty cell 𝒞_e, the maximum flow depth (defined as the maximum distance from the topography Σ in the surface-normal direction among all the material points within 𝒞_e, see Fig. 1c) can be found as

$$\begin{array}{}\text{(7)}& h(i_e,j_e,t_e)=\underset{{\mathit{X}}_p\left(t_e\right)\in {\mathcal{C}}_e}{max}\left[\left(Z_p\left(t_e\right)-H_p\right)\cdot n_{sH}(i_m,j_m)\right],\end{array}$$

where n_sH is the component of the normal-to-topography unit vector in the vertical direction.

Similarly, the maximum flow thickness (defined as the difference between the maximum and minimum distances from the topography Σ in the surface-normal direction among all material points within 𝒞_e, see Fig. 1c) is calculated as

$$\begin{array}{}\text{(8)}& \begin{array}{rl}t(i_e,j_e,t_e)& =h(i_e,j_e,t_e)-\underset{{\mathit{X}}_p\left(t_e\right)\in {\mathcal{C}}_e}{min}\left[\left(Z_p\left(t_e\right)-H_p\right)\right.\\ & \left.\cdot n_{sH}(i_m,j_m)\right].\end{array}\end{array}$$

In both Eqs. (7) and (8), we provisionally neglect the increase in flow height associated to the volume of the material points. The code also allows setting a percentile to get the flow depth and thickness from all heights and thickness values of all particles within each cell, thus potentially excluding scattered particles.

2.3.2 Flow velocities

The velocity U_p of each particle p within a cell 𝒞_e is decomposed into the slope-normal and slope-parallel components (Fig. 1b):

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathit{U}}_{\mathrm{pn}}({\mathit{X}}_p,t_e)=\left[{\mathit{U}}_p({\mathit{X}}_p,t_e)\cdot {\mathit{n}}_{\mathrm{s}}(i_m,j_m)\right]{\mathit{n}}_{\mathrm{s}}(i_m,j_m),\text{(10)}& {\displaystyle }{\mathit{U}}_{\mathrm{pt}}({\mathit{X}}_p,t_e)={\mathit{U}}_p({\mathit{X}}_p,t_e)-{\mathit{U}}_{\mathrm{pn}}({\mathit{X}}_p,t_e).\end{array}$$

The slope-normal and slope-parallel velocities can therefore be depth-integrated using (density-weighted) Favre-averaging:

$$\begin{array}{}\text{(11)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathit{U}}}_n(E,N,t_e)& ={\displaystyle \frac{\mathrm{1}}{{\int }_{\mathrm{0}}^h\mathit{\rho }(E,N,z,t_e)\mathrm{d}z}}\cdot \\ & \underset{\mathrm{0}}{\overset{h}{\int }}\mathit{\rho }(E,N,z,t_e){\mathit{U}}_n(E,N,z,t_e)\mathrm{d}z,\end{array}\text{(12)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathit{U}}}_t(E,N,t_e)& ={\displaystyle \frac{\mathrm{1}}{{\int }_{\mathrm{0}}^h\mathit{\rho }(E,N,z,t_e)\mathrm{d}z}}\cdot \\ & \underset{\mathrm{0}}{\overset{h}{\int }}\mathit{\rho }(E,N,z,t_e){\mathit{U}}_t(E,N,z,t_e)\mathrm{d}z,\end{array}\end{array}$$

where z is the coordinate normal to the topography and ρ is the flow bulk density. However, in the code, the integration is not performed in the normal-to-topography direction, but vertically for all particles within a cell 𝒞_e. Although this is formally not consistent with the definitions of the depth-averaged velocities of Eqs. (11) and (12), this calculation avoids inconsistencies on concave-up topographies (i.e., duplication of particles to multiple cells is avoided; see also the discussion in Hergarten, 2024, p. 783) and carries substantial simplifications in coding. We therefore calculate the depth-averaged velocities by replacing dz with the discrete particle volume 𝒱_p, and ρ with the particle density ρ_p:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathit{U}}}_n(i_e,j_e,t_e)& \approx {\displaystyle \frac{\mathrm{1}}{{\sum }_{{\mathit{X}}_p\left(t_e\right)\in {\mathcal{C}}_e}{\mathit{\rho }}_p({\mathit{X}}_p,t_e){\mathcal{V}}_p({\mathit{X}}_p,t_e)}}\cdot \\ & \sum _{{\mathit{X}}_p\left(t_e\right)\in {\mathcal{C}}_e}{\mathit{\rho }}_p({\mathit{X}}_p,t_e){\mathit{U}}_{\mathrm{pn}}({\mathit{X}}_p,t_e){\mathcal{V}}_p({\mathit{X}}_p,t_e)\\ & ={\displaystyle \frac{\mathrm{1}}{n_p\left(t_e\right)}}\sum _{{\mathit{X}}_p\left(t_e\right)\in {\mathcal{C}}_e}{\mathit{U}}_{\mathrm{pn}}({\mathit{X}}_p,t_e),\end{array}\end{array}$$

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathit{U}}}_t(i_e,j_e,t_e)& \approx {\displaystyle \frac{\mathrm{1}}{{\sum }_{{\mathit{X}}_p\left(t_e\right)\in {\mathcal{C}}_e}{\mathit{\rho }}_p({\mathit{X}}_p,t_e){\mathcal{V}}_p({\mathit{X}}_p,t_e)}}\cdot \\ & \sum _{{\mathit{X}}_p\left(t_e\right)\in {\mathcal{C}}_e}{\mathit{\rho }}_p({\mathit{X}}_p,t_e){\mathit{U}}_{\mathrm{pt}}({\mathit{X}}_p,t_e){\mathcal{V}}_p({\mathit{X}}_p,t_e)\\ & ={\displaystyle \frac{\mathrm{1}}{n_p\left(t_e\right)}}\sum _{{\mathit{X}}_p\left(t_e\right)\in {\mathcal{C}}_e}{\mathit{U}}_{\mathrm{pt}}({\mathit{X}}_p,t_e).\end{array}\end{array}$$

Since the (solid) mass of each particle, ρ_p𝒱_p, is constant in MPM¹, the expressions on the far right of Eqs. (13) and (14) are derived, where n_p is the number of particles in each cell 𝒞_e. Notably, in this context, the Favre-averaged velocities correspond to the cell-averaged velocities. In certain circumstances, however, the flow may dilate, which would lead to a reduction in the flow bulk density ρ from its initial value ρ_p0 – this change in bulk density is not always purely due to plastic volumetric deformation from the constitutive model (i.e., resulting in the particles' density changing from ρ_p0 to ρ_p; see, also, Li et al., 2021), but can also result from the formation of cracks and voids between granules of particles, as well as particles moving apart while air gets mixed into the flow (even though air is not explicitly modeled). Therefore, it would become necessary to replace ρ_p𝒱_p in Eqs. (13) and (14) with $[\mathit{\rho }\mathcal{V}={\mathit{\rho }}_{p\mathrm{0}}{\mathcal{V}}_{p\mathrm{0}}+{\mathit{\rho }}_a(\mathcal{V}-{\mathcal{V}}_{p\mathrm{0}}\left)\right]$, where ρ_a is the air density and (𝒱−𝒱_p0) is the extra-air volume ingested in the flow. Given that ρ_a is small and assuming limited flow dilation, we neglect the term ρ_a(𝒱−𝒱_p0).

2.3.3 Choice of exporting parameters for GIS visualization

GIS results are exported with a specified cell size cs_e. Choosing an appropriate value for cs_e is crucial: If the cells are too small, it may result in too many empty cells, especially in areas of the flow where the flow is dilated and particles are dispersed; Conversely, if the cells are too large, it might overestimate the flow extent, especially in areas where the flow has compacted. For an initial estimate, one can use the reference volume 𝒱 around a particle as a surrogate for cell size. This volume is derived from mass conservation of a MPM particle, neglecting the extra-air mass (see Sect. 2.3.2), such that ρ𝒱≈ρ_p0𝒱_p0, and therefore:

$$\begin{array}{}\text{(15)}& {\mathrm{cs}}_e\approx \sqrt[\mathrm{3}]{\mathcal{V}}\approx \sqrt[\mathrm{3}]{{\displaystyle \frac{{\mathit{\rho }}_{p\mathrm{0}}{\mathcal{V}}_{p\mathrm{0}}}{\mathit{\rho }}}}.\end{array}$$

As the flow bulk density ρ may not be constant spatially and temporally, the choice of cs_e requires some compromise and may cause the total flow mass (as derived from the exported flow height results) to be not conserved at all time steps. For instance, at the initial time step, results may be exported with cs${}_e\left(t_{\mathrm{0}}\right)\approx \sqrt[\mathrm{3}]{{\mathcal{V}}_{p\mathrm{0}}}=\sqrt[\mathrm{3}]{g{s}^{\mathrm{3}}/n_g}$. If, say, during the flow, ρ becomes smaller than ρ_p0, results should ideally be exported with cs_e(t_e)>cs_e(t₀). This becomes particularly difficult when results must be exported not just at a specific time step, but also require the export of spatial distribution of the maximum value of a certain depth-averaged variable ℱ (any of Eqs. 7, 8, 13, 14) over all time steps (this forcedly requires selecting a constant cs_e):

$$\begin{array}{}\text{(16)}& {\mathcal{F}}_{\max }(i_e,j_e)=\underset{t_e=k\cdot \mathrm{\Delta }{t}_e,k=\mathrm{0},\mathrm{1},\mathrm{\dots },\frac{t_{\mathrm{sim}}}{\mathrm{\Delta }{t}_e}}{max}\mathcal{F}(i_e,j_e,t_e),\end{array}$$

where t_sim is the simulated flow duration and Δt_e is the exported time step. Furthermore, to avoid spatial-loss of information of maximum values of the depth-averaged flow variables, Δt_e should be chosen as

$$\begin{array}{}\text{(17)}& \mathrm{\Delta }{t}_e\approx {\displaystyle \frac{{\mathrm{cs}}_e}{\Vert \stackrel{\mathrm{̃}}{\mathit{U}}{\Vert }_{\max }}},\end{array}$$

where $\Vert \stackrel{\mathrm{̃}}{\mathit{U}}{\Vert }_{\max }$ is the maximum (over time and space) depth-averaged speed projected in the EN plane.

3 Results and discussion

We first validate the exporting tool on a simplified case of a frictionless block moving on a 2D inclined plane and then falling from a cliff (see Supplement S1). Subsequently, we simulate in MPM the ice avalanche that resulted from the collapse of a 150 000 m³ portion of the Whymper serac (Courmayeur, Italy) in 1998 (Fig. 2a). The results are exported using cs_e=1 m (estimated from Eq. 15) and Δt_e=0.2 s and visualized in Fig. 2b–g. Since the current version of the code is not parallelized, and consequently considerably slow, we provisionally used a large time step Δt_e, about 10 times greater than what suggested by Eq. (17).

The simulated deposit is shown in Fig. 2b (flow depth h, normal to the terrain) and the flow inundation area is shown in Fig. 2c (plotted in terms of maximum depth-averaged slope-parallel flow speed). The runout and flow extent reasonably match the mapped extent of the event. Nonetheless, the simulated left bank lobe of the avalanche stops approximately 300 m up-slope compared to the mapped runout. Similarly, the right lobe of the simulated avalanche exhibits a shorter runout compared to the actual event. This reduced simulated mobility may result from excessive material deposition in the crevasses (the real terrain DTM at the time of the ice avalanche is unknown and could have had smaller crevasses than the DTM used for simulations) and not explicitly modeling the erodible glacier cover. Furthermore, we only made a few attempts to alter the parameters of the ice material in MPM (in particular, the residual friction and the basal frictions) without further optimization of the simulation results.

We now briefly analyze the influence of complex terrain topography on flow behavior, which can be effectively captured by the depth-averaged results, shown in Fig. 2d–i at t_e=50 s. In regions with curved topographies and jumps (e.g., cliffs, crevasses), high values of the depth-averaged slope-normal speed (see Fig. 2d) are calculated, whereas on more gentle terrain the ice flows approximately parallel to the terrain (see Fig. 2e). Specifically, on convex-up topographies along the flow direction (such as overhangs, or cliff edges), the ice flows in the out-of-slope direction (i.e., positive, blue values of ${\mathit{n}}_s^T\cdot {\stackrel{\mathrm{̃}}{\mathit{U}}}_n$). In contrast, on concave-up sections (such as gullies or terrain depressions) the flow is directed toward the topography (i.e., negative, red values of ${\mathit{n}}_s^T\cdot {\stackrel{\mathrm{̃}}{\mathit{U}}}_n$). Figure 2h shows the calculated depth-averaged velocities extracted along a linear profile of the glacier (the profile line is shown in Fig. 2d–g). A cliff is located between ≈60 and ≈80 m of this profile, causing the flow to separate from the topography (${\mathit{n}}_s^T\cdot {\stackrel{\mathrm{̃}}{\mathit{U}}}_n>\mathrm{0}$) at the jump and subsequently collide with the topography (${\mathit{n}}_s^T\cdot {\stackrel{\mathrm{̃}}{\mathit{U}}}_n<\mathrm{0}$) upon landing. In other areas where the curvature of the topography is smaller and there are no terrain jumps, the calculated slope-normal component is smaller, and mostly negligible compared to the slope-parallel velocity component. The flow detachment from the topography is also evident by comparing the flow thickness (Fig. 2f) to the flow depth (Fig. 2g). Figure 2i presents the comparison of these two parameters along the profile line: across the jump and upon landing (splash), the flow depth significantly exceeds the flow thickness – signifying the flow's detachment from the surface (see also Fig. 1c) – whereas the values of t and h are not so far apart when the flow sticks to the terrain (see also the additional analysis in the Supplement S3). Unlike in depth-averaged models, the detachment from the terrain in our 3D simulations allows to correctly capture that, while in airborne motion, the material does not experience basal flow resistance. Similarly, when the flow impacts the terrain, the 3D MPM can effectively capture the predominant motion towards the terrain, potentially leading to increased basal flow resistance, as well as compaction and/or dilation of the material. This may lead to more precise dynamics compared to depth-averaged models, particularly on complex terrains that exhibit jumps and substantial curvature². These spatial insights highlight how our export tool reveals terrain-driven flow dynamics that are not easily interpreted directly from raw 3D MPM outputs.

Figure 2 Depth-resolved MPM simulations and GIS visualization of exported depth-averaged results. (a) 3D visualization at t=50 s of the MPM simulation of the Whymper ice avalanche in 1998. (b) Exported final deposition depth (h(t_e=150 s)). (c) Maximum depth-averaged slope-parallel speed over all time steps (${\stackrel{\mathrm{̃}}{\mathit{U}}}_{t,max}$). (d) Depth-averaged slope-normal velocity (${\stackrel{\mathrm{̃}}{\mathit{U}}}_n$) at t_e=50 s. (e) Depth-averaged slope-parallel speed (${\stackrel{\mathrm{̃}}{\mathit{U}}}_t$) at t_e=50 s. (f) Flow thickness (t) at t_e=50 s. (g) Flow depth (h) at t_e=50 s. (h) Slope-normal (${\stackrel{\mathrm{̃}}{\mathit{U}}}_n$) and slope-parallel (${\stackrel{\mathrm{̃}}{\mathit{U}}}_t$) speeds at t_e=50 s along a profile line. (i) Flow thickness (t) and depth (h) at t_e=50 s along a profile line (black line in d–g).

4 Conclusion and outlook

In this study, we presented a methodology to transform 3D MPM simulation outputs into 2D rasterized formats that can be visualized with GIS tools. This methodology entails the conversion of depth-resolved MPM variables into depth-averaged variables. More specifically, we described the computation of flow height metrics, like flow depth and thickness, and flow velocity metrics, like Favre-averaged slope-normal and slope-parallel velocities. Thus, this approach facilitates harvesting the benefits of 3D numerical methods (Li et al., 2021; Cicoira et al., 2022), including the simulation of geophysical flows over complex topographies, terrain jumps, cavities, and interacting with mitigation structures, while enabling practitioners to visualize the model results on GIS maps similar to traditional depth-averaged simulation outputs. Therefore, this tool represents a step forward in using 3D MPM for hazard assessment and mapping of geophysical flows.

The application of the depth-averaging tool to the simulation of an ice avalanche highlighted how complex terrain, like with cliffs and crevasses, can significantly alter flow velocities, causing them to deviate from being parallel to the ground. The computation of slope-normal velocities and the distinction of flow depth (distance between the terrain surface and the highest flow particle) and flow thickness (projected distance between the highest and lowest particles in each cell) revealed that the material might detach from the terrain during jumps and impact the terrain upon landing, potentially influencing flow dynamics through dilation, compaction, and variations in internal and basal stresses. Additional research will be needed to study in detail the effect of complex topographical features on avalanche characteristics. Furthermore, we observed that crevasses along the flow path could lead to significant deposition, further influencing subsequent flow dynamics. Thus, selecting an appropriate DTM appears critical for simulating ice avalanches over glaciers: this needs to be further investigated in the future. Alternative constitutive models for ice, such as incorporating rate dependency, will be required in future studies to achieve better agreement between the simulated ice avalanche and the actual observed runout.

The depth integration of 3D simulation results presented in this work will allow future comparison between 3D MPM simulations and depth-averaged simulations (Wirbel et al., 2024) to verify the assumptions made in depth-averaged numerical models, such as neglecting flow velocities normal to the terrain, assuming the flow sticking to the terrain, constant bulk flow density, and uniform velocity profiles. Additional work is required to broaden the utility of our depth-averaging tool. In the future, we plan to export density and velocity profiles along the flow depth direction. These data could be crucial for determining the pressure distribution of avalanches impacting obstacles. Such vertical information (e.g., Kyburz et al., 2024) may, for example, be exported onto raster, to replace traditional empirical analytical models of impact pressure, which (typically) rely solely on depth-averaged density and velocity. The tool could also be used to determine topography changes and erosion rates when an erodible bed layer is explicitly modeled in MPM (e.g., Li et al., 2022).

Finally, a similar framework could be applied in the future to depth-average other 3D particle-based models' results, like SPH and DEM (Discrete Element Method). Recently, similar algorithms were proposed by Han et al. (2020) and Su et al. (2024) to extract depth-averaged physical quantities from 3D SPH simulations. Specifically, the algorithm by Han et al. (2020) assigns each 3D SPH particle to a boundary particle by projecting particle positions along the terrain-normal direction. Physical quantities (e.g., flow depth and basal velocities) are then computed for each boundary particle based on all SPH particles falling within its associated cell. In contrast, Su et al. (2024) identifies contributing SPH particles as those located within a cylindrical region whose longitudinal axis is centered on each boundary particle and aligned with the terrain normal. Our approach differs in that MPM particles are assigned to 2D raster cells based on vertical projection. While the computation of flow depth and velocity components is conceptually similar across these methods, our implementation includes additional features: the ability to exclude scattered particles from the flow depth calculation, and the explicit computation of the flow thickness. While Han et al. (2020) and Su et al. (2024) assign results to boundary particles, they can also export the data to 2D grids visualized in ParaView. Our method similarly produces rasterized outputs and can directly map quantities to a georeferenced 2D grid, which can be seamlessly exported to a GIS environment. Moreover, in contrast to the mesh-free SPH method, MPM solves momentum conservation on a background Eulerian grid. Thus, an alternative approach could rely on the direct extraction of physical quantities from the MPM code, straightforwardly utilizing the transferred velocities and masses at the grid nodes – thus eliminating the need for a separate export tool.