com1DFA DFA-Kernel theory

Warning

This theory has not been fully reviewed yet. Read its content with a critical mind.

Governing Equations for the Dense Flow Avalanche

The governing equations of the dense flow avalanche are derived from the incompressible mass and momentum balance on a Lagrange control volume ([Zwi00, ZKS03]).

Mass balance:

(1)\[\frac{d}{dt} \int\limits_{V(t)} \rho_0 \,\mathrm{d}V = \rho_0 \frac{dV(t)}{dt} = \oint\limits_{\partial V(t)} q^{\text{ent}} \,\mathrm{d}A\]

Where \(q^{\text{ent}}\) represents the snow entrainment rate.

Momentum balance:

(2)\[\rho_0 \frac{d}{dt} \int\limits_{V(t)} u_i \,\mathrm{d}V = \oint\limits_{\partial V(t)} \sigma^{\text{tot}}_{ij}n_j \,\mathrm{d}A + \rho_0 \int\limits_{V(t)} g_i \,\mathrm{d}V, \quad i=(1,2,3)\]

We introduce the volume average of a quantity \(P(\mathbf{x},t)\):

\[\overline{P}(\mathbf{x},t) = \frac{1}{V(t)} \int\limits_{V(t)} P(\mathbf{x},t) \,\mathrm{d}V\]

and split the area integral into :

\[\oint\limits_{\partial V(t)} \sigma^{\text{tot}}_{ij}n_j \,\mathrm{d}A = \oint\limits_{\partial V(t)} \sigma_{ij}n_j \,\mathrm{d}A + F_i^{\text{ent}} + F_i^{\text{res}}, \quad i=(1,2,3)\]

\(F_i^{\text{ent}}\) represents the force required to break the entrained snow from the ground and to compress it (since the dense-flow bulk density is usually larger than the density of the entrained snow, i.e. \(\rho_{\text{ent}}<\rho\)) and \(F_i^{\text{res}}\) represents the resistance force due to obstacles (for example trees). Which leads to in Eq.2:

\[\rho_0 \frac{dV(t) \overline{u}_i}{dt} = \rho_0 V \frac{d\overline{u}_i}{dt} + \rho_0 \overline{u}_i \frac{dV}{dt} = \oint\limits_{\partial V(t)} \sigma_{ij}n_j \,\mathrm{d}A + \rho_0 V g_i + F_i^{\text{ent}} + F_i^{\text{res}}, \quad i=(1,2,3)\]

Using the mass balance equation Eq.1, we get:

(3)\[\rho_0 V \frac{d\overline{u}_i}{dt} = \oint\limits_{\partial V(t)} \sigma_{ij}n_j \,\mathrm{d}A + \rho_0 V g_i + F_i^{\text{ent}} + F_i^{\text{res}} - \overline{u}_i \oint\limits_{\partial V(t)} q^{\text{ent}} \,\mathrm{d}A, \quad i=(1,2,3)\]

Boundary conditions:

The free surface is defined by :

\[F_s(\mathbf{x},t) = z-s(x,y,t)=0\]

The bottom surface is defined by :

\[F_b(\mathbf{x}) = z-b(x,y)=0\]

The boundary conditions at the free surface and bottom of the flow read:

(4)\[\begin{split}\left\{\begin{aligned} &\frac{dF_s}{dt} = \frac{\partial F_s}{\partial t} + u_i\frac{\partial F_s}{\partial x_i} =0 \quad & \mbox{at }F_s(\mathbf{x},t) =0 \quad & \mbox{Kinematic BC (Material boundary)}\\ &\sigma_{ij}n_j = 0 \quad & \mbox{at }F_s(\mathbf{x},t) =0 \quad & \mbox{Dynamic BC (Traction free surface)}\\ &u_in_i = 0 \quad & \mbox{at }F_b(\mathbf{x},t) =0 \quad & \mbox{Kinematic BC (No detachment)}\\ &\tau^{(b)}_i = f(\sigma^{(b)},\overline{u},\overline{h},\rho_0,t,\mathbf{x})\quad & \mbox{at }F_b(\mathbf{x},t) =0\quad & \mbox{Dynamic BC (Chosen friction law)} \end{aligned} \right.\end{split}\]

\(\sigma^{(b)}_i = (\sigma_{kl}n_ln_k)n_i\) represents the normal stress at the bottom and \(\tau^{(b)}_i = \sigma_{ij}n_j - \sigma^{(b)}_i\) represents the shear stress at the bottom surface. \(f\) describes the chosen friction model and are described in [sec:friction model]. The normals at the free surface (\(n_i^{(s)}\)) and bottom surface (\(n_i^{(b)}\)) are:

\[n_i^{(s,b)} = \frac{\partial F_{s,b}}{\partial x_i}\left(\frac{\partial F_{s,b}}{\partial x_j} \frac{\partial F_{s,b}}{\partial x_j}\right)^{-1/2}\]

Choice of the coordinate system:

The previous equations will be developed in the orthonormal coordinate system \((B,\mathbf{e_1},\mathbf{e_2},\mathbf{e_3})\), further referenced as Natural Coordinate System (NCS). In this NCS, \(\mathbf{e_1}\) is aligned with the velocity vector at the bottom and \(\mathbf{e_3}\) with the normal to the slope, i.e.:

\[\mathbf{e_1} = \frac{\mathbf{u}}{\left\Vert \mathbf{u}\right\Vert},\quad \mathbf{e_2} = \mathbf{e_3}\wedge\mathbf{e_1}, \quad \mathbf{e_3} = \mathbf{n^{(b)}}\]

The origin \(B\) of the NCS is attached to the slope. This choice leads to:

\[n^{(b)}_i = \delta_{i3}, \quad \left.\frac{\partial b}{\partial x_i}\right\rvert_{\mathbf{0}} = 0\quad \mbox{for} \quad i=(1,2),\quad \mbox{and} \quad u^{(b)}_2 = u^{(b)}_3 = 0\]

In this NCS and considering a prism-like Control volume, the volume content \(V(t) = A_b(t)\overline{h}\) is obtained by multiplication of the basal area of the prism, \(A_b\), with the averaged value of the flow-depth,

(5)\[\overline{h} = \frac{1}{A_b(t)}\int\limits_{A_b(t)} [s(\mathbf{x})-b(\mathbf{x})]\,\mathrm{d}A = \frac{1}{A_b(t)}\int\limits_{A_b(t)} h(\mathbf{x})\,\mathrm{d}A,\qquad \overline{u}_i = \frac{1}{V(t)}\int\limits_{V(t)} u_i(\mathbf{x})\,\mathrm{d}V\]

The snow entrainment is either due to plowing at the front of the avalanche or to erosion at the bottom. The entrainment rate at the front \(q^{\text{plo}}\) can be expressed function of the properties of the entrained snow (density \(\rho_{\text{ent}}\) and snow depth \(h_{\text{ent}}\)), the velocity of the avalanche at the front \(\overline{\mathbf{u}}\) and length \(w_f\) of the front cell (measured perpendicularly to the flow velocity \(\overline{\mathbf{u}}\)). It obviously only happens on the front cells of the avalanche (meaning that \(w_f\) is zero for inner parts of the avalanche):

(6)\[\oint\limits_{\partial V(t)} q^{\text{plo}}\,\mathrm{d}A = \int\limits_{l_{\text{front}}}\int_b^s q^{\text{plo}}\, \mathrm{d}{l}\,\mathrm{d}{z} = \rho_{\text{ent}}\,w_f\,h_{\text{ent}}\,\left\Vert \overline{\mathbf{u}}\right\Vert\]

The entrainment rate at the bottom \(q^{\text{ero}}\) can be expressed function of the bottom area \(A_b\) of the control volume, the velocity of the avalanche \(\overline{\mathbf{u}}\), the bottom shear stress \(\tau^{(b)}\) and the specific erosion energy \(e_b\)):

(7)\[\oint\limits_{\partial V(t)} q^{\text{ero}}\,\mathrm{d}A = \int\limits_{A_b} q^{\text{ero}}\, \mathrm{d}A = A_b\,\frac{\tau^{(b)}}{e_b}\,\left\Vert \overline{\mathbf{u}}\right\Vert\]

This leads in the mass balance Eq.1 to :

(8)\[\frac{\mathrm{d}V(t)}{\mathrm{d}t} = \frac{\mathrm{d}(A_b\overline{h})}{\mathrm{d}t} = \frac{\rho_{\text{ent}}}{\rho_0}\,w_f\,h_{\text{ent}}\,\left\Vert \overline{\mathbf{u}}\right\Vert + A_b\,\frac{\tau^{(b)}}{e_b}\,\left\Vert \overline{\mathbf{u}}\right\Vert\]

The force \(F_i^{\text{ent}}\) required to break the entrained snow from the ground and to compress it is expressed function of the required breaking energy per fracture surface unit \(e_s\) (\(J.m^{-2}\)), the deformation energy per entrained mass element \(e_d\) (\(J.kg^{-1}\)) and the entrained snow depth ([FFGS13, SFF+08, Sam07]):

\[F_i^{\text{ent}} = -w_f\,(e_s+\,q^{\text{ent}}\,e_d)\]

The force \(F_i^{\text{res}}\) due to obstacles is expressed function of the characteristic diameter \(\overline{d}\) and height \(h_{\text{res}}\) of the obstacles, the spacing \(s_{\text{res}}\) between the obstacles and an empirical coefficient \(c_w\) (see Fig. 4). The effective height \(h^{\text{eff}}\) is defined as \(\min\left\lbrace\begin{array}{l} \overline{h}\\h_{res}\end{array}\right\rbrace\):

\[F_i^{\text{res}} = -(\frac{1}{2}\,\overline{d}\,c_w/s^2_{\text{res}})\,\rho_0\,A\, h^{\text{eff}}\,\overline{u}^2\, \frac{\overline{u}_i}{\|\overline{u}\|}\]

Fig. 4 Resistance force due to obstacles (from [FiKo2013])

The surface integral is split in three terms, an integral over \(A_b\) the bottom \(x_3 = b(x_1,x_2)\), \(A_s\) the top \(x_3 = s(x_1,x_2,t)\) and \(A_h\) the lateral surface. Introducing the boundary conditions Eq.4 leads to:

\[\begin{split}\begin{aligned} \oint\limits_{\partial{V(t)}}\sigma_{ij}n_j\,\mathrm{d}A & = \int\limits_{A_b}\underbrace{\sigma_{ij}\,n_j^{(b)}}_{-\sigma_{i3}}\,\mathrm{d}A + \int\limits_{A_s}\underbrace{\sigma_{ij}\,n_j^{(s)}}_{0}\,\mathrm{d}A + \int\limits_{A_h}\sigma_{ij}\,n_j\,\mathrm{d}A\\ &= -A_b\overline{\sigma}_{i3}^{(b)} + \oint\limits_{\partial A_b}\left(\int_b^s\sigma_{ij}\,n_j\,\mathrm{d}x_3\right)\,\mathrm{d}l \end{aligned}\end{split}\]

Which simplifies the momentum balance Eq.3 to:

(9)\[\begin{split}\begin{aligned} \rho_0 V \frac{d\overline{u}_i}{dt} = & \oint\limits_{\partial A_b}\left(\int_b^s\sigma_{ij}\,n_j\, \mathrm{d}x_3\right)\,\mathrm{d}l -A_b\overline{\sigma}_{i3}^{(b)} + \rho_0 V g_i + F_i^{\text{ent}} + F_i^{\text{res}} - \overline{u}_i \oint\limits_{\partial V(t)} q^{\text{ent}} \,\mathrm{d}A,\\ &\quad i=(1,2,3) \end{aligned}\end{split}\]

The momentum balance in direction \(x_3\) (normal to the slope) is used to obtain a relation for the vertical distribution of the stress tensor ([Sam07]). Due to the choice of coordinate system and because of the kinematic boundary condition at the bottom, the left side of Eq.9 can be expressed function of the velocity \(\overline{u}_1\) in direction \(x_1\) and the curvature of the terrain in this same direction \(\frac{\partial^2{b}}{\partial{x_1^2}}\) ([Zwi00]):

\[\rho\,A_b\,\overline{h}\,\frac{\,\mathrm{d}\overline{u}_3}{\,\mathrm{d}t} = \rho\,A_b\,\overline{h}\,\frac{\partial^2{b}}{\partial{x_1^2}}\,\overline{u}_1^2,\]

rearranging the terms in the momentum equation leads to:

(10)\[\overline{\sigma}_{33}(x_3) = \rho_0\,(s-x_3)\left(g_3-\frac{\partial^2{b}}{\partial{x_1^2}}\,\overline{u}_1^2\right)+ \frac{1}{A_b} \oint\limits_{\partial A_b}\left(\int_{x_3}^s\sigma_{3j}\,n_j\,\mathrm{d}x_3\right)\,\mathrm{d}l\]

Non-dimensional Equations

Fig. 5 Characteristic size of the avalanche along its path (from [Zw2000], modified)

The previous equations Eq.9 and Eq.10 are can be further simplified by introducing a scaling based on the characteristic values of the physical quantities describing the avalanche. The characteristic length L, the thickness H, the acceleration due to gravity g and the characteristic radius of curvature of the terrain R are the chosen quantities. From those values, it is possible to form two non dimensional parameters that describe the flow:

• Aspect ratio: \(\qquad\qquad\varepsilon = H / L\qquad\)

• Curvature: \(\qquad\lambda = L / R\qquad\)

The different properties involved are then expressed in terms of characteristic quantities \(L\), \(H\), \(g\), \(\rho_0\) and \(R\) (see Fig. 5):

\[\begin{split}\begin{aligned} x_i &= L\, x_i^*\\ (dx_3,h,\overline{h}) &= H\,(dx_3^*,h^*,\overline{h}^*)\\ A_b &= L^2\, A_b^*\\ t &= \sqrt{L/\text{g}}\, t^*\\ \overline{u_i} &= \sqrt{\text{g}L}\,\overline{u_i}^*\\ \text{g}_i &= \text{g} \, \text{g}_i^*\\ \frac{\partial^2{b}}{\partial{x_1}^2} &= \frac{1}{R}\,\frac{\partial^2{b^*}}{\partial{x_1}^{*2}}\end{aligned}\end{split}\]

The normal part of the stress tensor is directly related to the hydro-static pressure:

\[\sigma_{ii} = \rho_0\,\text{g}\,H\,\sigma_{ii}^*\]

The dimensionless properties are indicated by a superscripted asterisk. Introducing those properties in Eq.10, leads to :

(11)\[\overline{\sigma^*}_{33} = \left(g^*_3-\lambda\frac{\partial^2{b^*}}{\partial{x_1^{*2}}}\,\overline{u}_1^{*2}\right) (s^*-x^*_3) + \underbrace{\varepsilon\oint\limits_{\partial A_b^*}\left(\int\limits_{x^*_3}^{s^*}\sigma^*_{31}\,\mathrm{d}x^*_3\right)\,\mathrm{d}l^*}_{O(\varepsilon)}.\]

The height, H of dense flow avalanches is assumed to me small compared to its length, L. Meaning that the equations are examined in the limit \(\varepsilon \ll 1\). It is then possible to neglect the last term in Eq.11 which leads to (after reinserting the dimensions):

(12)\[\overline{\sigma}_{33}(x_3) = \rho_0\,\underbrace{\left(g_3-\overline{u_1}^2\,\frac{\partial^2{b}}{\partial{x_1^2}}\right)}_{g_\text{eff}} \left[\overline{h}-x_3\right]\]

And at the bottom of the avalanche, with \(x_3 = 0\), the normal stress can be expressed as:

(13)\[\overline{\sigma}^{(b)}_{33} = \rho_0\,\left(g_3-\overline{u_1}^2\,\frac{\partial^2{b}}{\partial{x_1^2}}\right)\,\overline{h}\]

Calculating the surface integral in equation Eq.9 requires to express the other components of the stress tensor. Here again a magnitude consideration between the shear stresses \(\sigma_{12} = \sigma_{21}\) and \(\sigma_{13}\). The shear stresses are based on a generalized Newtonian law of materials, which controls the influence of normal stresse and the rate of deformation through the viscosity.

\[\tau_{ij} = \eta\left(\frac{\partial{u_i}}{\partial{x_j}}+\frac{\partial{u_j}}{\partial{x_i}}\right), ~ i\neq j\]

Because \(\partial x_1\) and \(\partial x_2\) are of the order of \(L\), whereas \(\partial x_3\) is of the order of \(H\), it follows that:

\[O\left(\frac{\sigma_{12}}{\sigma_{13}}\right) = \frac{H}{L} = \varepsilon \ll 1\]

and thus \(\sigma_{12} = \sigma_{21}\) is negligible compared to \(\sigma_{13}\). \(\sigma_{13}\) is exressed using the bottom friction law \(\tau^{(b)}_i = f(\sigma^{(b)},\overline{u},\overline{h},\rho_0,t,\mathbf{x})\) introduced in Eq.4.

In addition, a relation linking the horizontal normal stresses, \(\sigma_{ii}\), \(i = (1,2)\), to the vertical pressure distribution given by Eq.13 is introduced. In complete analogy to the arguments used by Savage and Hutter ([SH89]) the horizontal normal stresses are given as:

\[\sigma_{ii} = K_{(i)}\,\sigma_{33}\]

Where \(K_{(i)}\) are the earth pressure coefficients (cf. [Sal04, ZKS03]):

\[\begin{split}\sigma_{11} &= K_{x~akt/pass}\,\sigma_{33}\\ \sigma_{22} &= K_{y~akt/pass}^{(x~akt/pass)}\,\sigma_{33}\end{split}\]

With the above specifications, the integral of the stresses over the flow height is simplified in equation Eq.9 to:

\[\int\limits_b^s\sigma_{ij}\,\mathrm{d}x_3 = \int\limits_b^s K_{(i)}\,\sigma_{33}\,\mathrm{d}x_3 = K_{(i)}\,\frac{\overline{h}\,\sigma^{(b)}}{2}\]

and the momentum balance can be written:

(14)\[\begin{split}\begin{aligned} \rho_0\,A\,\overline{h}\,\frac{\,\mathrm{d}\overline{u}_i}{\,\mathrm{d}t} = &\rho_0\,A\,\overline{h}\,g_i + \underbrace{K_{(i)}\,\oint\limits_{\partial{A}}\left(\frac{\overline{h}\,\sigma^{(b)}}{2}\right)n_i\,\mathrm{d}l}_{F_i^{\text{lat}}} \underbrace{-\delta_{i1}\,A\,\tau^{(b)}}_{F_i^{\text{bot}}} \underbrace{- \rho_0\,A\,h_{\text{eff}}\,C_{\text{res}}\,\overline{\mathbf{u}}^2\,\frac{\overline{u_i}}{\|\overline{\mathbf{u}}\|}}_{F_i^{\text{res}}}\\ &- \overline{u_i}\,\rho_0\,\frac{\mathrm{d}\left(A\,\overline{h}\right)}{\mathrm{d}t} + F_i^{\text{ent}} \end{aligned}\end{split}\]

with

\[C_{\text{res}} = \frac{1}{2}\,\overline{d}\,\frac{c_w}{s_{\text{res}}^2}.\]

The mass balance Eq.8 remains unchanged:

(15)\[\frac{\mathrm{d}V(t)}{\mathrm{d}t} = \frac{\mathrm{d}\left(A_b\overline{h}\right)}{\mathrm{d}t} = \frac{\rho_{\text{ent}}}{\rho_0}\,w_f\,h_{\text{ent}}\,\left\Vert \overline{\mathbf{u}}\right\Vert + A_b\,\frac{\tau^{(b)}}{e_b}\,\left\Vert \overline{\mathbf{u}}\right\Vert\]

The unknown \(\overline{u}_1\), \(\overline{u}_2\) and \(\overline{h}\) satisfy Eq.13, Eq.14 and Eq.15. In equation Eq.14 the ground shear stress \(\tau^{(b)}\) remains unknown, and and a constitutive equation has to be introduced in order to completely solve the equations.

Friction Model

The problem can be solved by introducing a constitutive equation which describes the basal shear stress tensor \(\tau^{(b)}\) as a function of the flow state of the avalanche.

(16)\[\tau^{(b)}_i = f(\sigma^{(b)},\overline{u},\overline{h},\rho_0,t,\mathbf{x})\]

With

\[\begin{split}\begin{aligned} &\sigma^{(b)} \qquad &\text{normal component of the stress tensor}\\ &\overline{u} \qquad &\text{average velocity}\\ &\overline{h} \qquad &\text{average flow depth}\\ &\rho_0 \qquad &\text{density}\\ &t \qquad &\text{time}\\ &\mathbf{x} \qquad &\text{position vector}\end{aligned}\end{split}\]

Several friction models already implemented in the simulation tool are described here.

SamosAT Model

SamosAT friction model is a modification of some more clasical models such as Voellmy model. The basal shear stress tensor \(\tau^{(b)}\) is expressed as ([Sam07]):

\[\tau^{(b)} = \tau_0 + \tan{\delta}\,\left(1+\frac{R_s^0}{R_s^0+R_s}\right)\,\sigma^{(b)} + \frac{\rho_0\,\overline{u}^2}{\left(\frac{1}{\kappa}\,\ln\frac{\overline{h}}{R} + B\right)^2}\]

With

\[\begin{split}\begin{aligned} &\tau_0 \qquad &\text{minimum shear stress}\\ &R_s \qquad &\text{relation between friction and normal pressure (fluidization factor)}\\ &R \qquad &\text{empirical constant}\\ &R_s^0 \qquad &\text{empirical constant}\\ &B \qquad &\text{empirical constant}\\ &\kappa \qquad &\text{empirical constant}\end{aligned}\end{split}\]

The minimum shear stress \(\tau_0\) defines a lower limit below which no flow takes place with the condition \(\rho_0\,\overline{h}\,g\,\sin{\alpha} > \tau_0\). \(\alpha\) being the slope. \(\tau_0\) is independent of the flow depth, which leeds to a strong avalanche deceleration, especially for avalanches with low flow heights. \(R_s\) is expressed as \(R_s = \frac{\rho_0\,\overline{u}^2}{\sigma^{(b)}}\). Together with the empirical parameter \(R_s^0\) the term \(\frac{R_s^0}{R_s^0+R_s}\) defines the Coulomb basal friction. Therefore lower avalanche speeds lead to a higher bed friction, making avalanche flow stop already at steeper slopes \(\alpha\), than without this effect. This effect is intended to avoid lateral creep of the avalanche mass ([SG09])

Numerics

Principle

Mass Eq.15 and momentum Eq.14 balance equations as well as basal normal stress Eq.13 are solved numerically using a SPH method (Smoothed Particle Hydrodynamis) ([Mon92]) for the three variables \(\overline{\mathbf{u}}=(\overline{u}_1, \overline{u}_2)\) and \(\overline{h}\) by discretization of the released avalanche volume in a large number of mass elements. SPH in general, is a mesh-less numerical method for solving partial differential equations. The SPH algorithm discretizes the numerical problem within a domain using particles ([Sam07, SG09]), which interact with each-other in a defined zone of influence. Some of the advantages of the SPH method are that free surface flows, material boundaries and moving boundary conditions are considered implicitly. In addition, large deformations can be modeled due to the fact that the method is not based on a mesh. From a numerical point of view, the SPH method itself is relatively robust.

Method and discretization

Space discretization

The domain is discretized in particles. Each particle \(p_j\) is affected with the following properties: a mass \(m_{p_j}\), a depth \(\overline{h}_{p_j}\), a density \(\rho_{p_j}=\rho_0\) and a velocity \(\mathbf{\overline{u}}_{p_j}=(\overline{u}_{p_j,1}, \overline{u}_{p_j,2})\). Those particles are also projected on a regular grid (raster) and the mass distributed on each node of the raster (see Fig. 6). This leads to the following expression for the mass \(m_{c_n}\) of each node on the raster cell \(m_{c_n} = \sum\limits_{p_j}{m_{p_j}}\) (cells are indexed with \({c_n}\)).

Fig. 6 Particles in raster grid (from [FiKo2013])

Each grid node is also affected with a velocity \(\overline{\mathbf{u}}_{c_n}\) expressed as the sum of the momentum of each raster cell divided by the mass of the same cell:

\[\overline{\mathbf{u}}_{c_n} = \frac{\sum\limits_{p_j}{m_{p_j}}\overline{\mathbf{u}}_{p_j}}{\sum\limits_{p_j}{m_{p_j}}}\]

The flow depth \(\overline{h}_{c_n}\) can be deduced from the mass and area of the raster cell:

\[\overline{h}_{c_n} = \frac{m_{c_n}}{\rho_0\,A_{c_n}}\]

Back to the particles properties, the bottom area paired to each particle is related to the mass and flow depth of this one:

\[A_{p_j} = \frac{m_{p_j}}{\rho_0\,\overline{h}_{p_j}}\]

Method

The SPH method is introduced when expressing the flow depth and its gradient for each particle as a weighted sum of its neighbors ([LL10, Sam07]). The \(p\) in \(p_i\) is dropped (same applies for \(p_j\)):

(17)\[\begin{split}\overline{h}_i &= \frac{1}{\rho_0}\,\sum\limits_j{m_j}\,W_{ij}\\ \mathbf{\nabla}\overline{h}_i &= -\frac{1}{\rho_0}\,\sum\limits_{j}{m_j}\,\mathbf{\nabla}W_{ij}\end{split}\]

Where \(W\) represents the SPH-Kernel function and reads:

(18)\[\begin{split}W_{ij} = W(\mathbf{r_{ij}},r_0) = \frac{10}{\pi r_0^5}\left\{ \begin{aligned} & (r_0 - \left\Vert \mathbf{r_{ij}}\right\Vert)^3, \quad &0\leq \left\Vert \mathbf{r_{ij}}\right\Vert \leq r_0\\ & 0 , & r_0 <\left\Vert \mathbf{r_{ij}}\right\Vert \end{aligned} \right.\end{split}\]

\(\left\Vert \mathbf{r_{ij}}\right\Vert= \left\Vert \mathbf{x_i}-\mathbf{x_j}\right\Vert\) represents the distance between particle \(i\) and \(j\) and \(r_0\) the smoothing length. And its first derivative reads:

(19)\[\begin{split}\mathbf{\nabla}W_{ij} = -3.\frac{10}{\pi r_0^5}\left\{ \begin{aligned} & (r_0 - \left\Vert \mathbf{r_{ij}}\right\Vert)^2\,\frac{\mathbf{r_{ij}}}{r_{ij}}, \quad &0\leq \left\Vert \mathbf{r_{jl}}\right\Vert \leq r_0\\ & 0 , & r_0 <\left\Vert \mathbf{r_{ij}}\right\Vert \end{aligned} \right.\end{split}\]

The lateral pressure forces on each particle are calculated from the compression forces on the boundary of the particle. The boundary is approximated as a square with the base side length \(\Delta s = \sqrt{A_p}\) and the respective flow height. This leads to (subscript \(|_{,d}\) stands for the component in the \(d^{th}\) direction, \(d = {1,2}\)):

\[F_{i,d}^{\text{lat}} = K_{(d)}\oint\limits_{\partial{A_{i}}}\left(\int\limits_{b}^{s}\sigma_{33}\,n_d\,\mathrm{d}x_3\right)\mathrm{d}l\]

From equation Eq.14

\[F_{i,d}^{\text{lat}} = K_{(d)}\,\frac{\Delta{s}}{2}\left((\overline{h}\,\overline{\sigma}^{(b)}_{33})_{x_{d}- \frac{\Delta{s}}{2}}-(\overline{h}\,\overline{\sigma}^{(b)}_{33})_{x_{d}+\frac{\Delta{s}}{2}}\right) = K_{(d)}\frac{\Delta{s}^2}{2}\,\left.\frac{d\,\overline{h}\,\overline{\sigma}^{(b)}}{d\,x_d}\right\rvert_{i}\]

The product of the average flow depth \(\overline{h}\) and the basal normal pressure \(\overline{\sigma}^{(b)}_{33}\) reads (using equation Eq.13 and dropping the curvature acceleration term):

\[\overline{h}\,\overline{\sigma}^{(b)} = \overline{h}^2\,\rho_0\,\left(g_3-\overline{u_1}^2\,\frac{\partial^2{b}}{\partial{x_1^2}}\right) \approx \overline{h}^2\,\rho_0\,g_3\]

Which leads to, using the relation Eq.17:

(20)\[F_{i,d}^{\text{lat}} = K_{(d)}\,\rho_0\,g_3\,A_{i}\,\overline{h}_{i}\,.\,\left.\frac{d\,\overline{h}}{d\,x_d}\right\rvert_{i} = -K_{(d)}\,m_{i}\,g_3\,.\,\frac{1}{\rho_0}\,\sum\limits_{j}{m_{j}}\,\left.\frac{d\,W_{ij}}{d\,x_d}\right\rvert_{j}\]

The bottom friction forces on each particle depend on the chose friction model and reads for the SamosAT friction model (using equation Eq.13 for the expression of \(\sigma^{(b)}_{i}\)):

(21)\[F_{i,d}^{\text{bot}} = -\delta_{d1}\,A_{i}\,\tau^{(b)}_{i} = -\delta_{d1}\,A_{i}\,\left(\tau_0 + \tan{\delta}\,\left(1+\frac{R_s^0}{R_s^0+R_s}\right)\,\sigma^{(b)}_{i} + \frac{\rho_0\,\mathbf{\overline{u}}_{i}^2}{\left(\frac{1}{\kappa}\,\ln\frac{\overline{h}}{R} + B\right)^2}\right)\]

The resistance force on each particle reads (where \(h^{\text{eff}}_{j}\) is a function of the average flow depth \(\overline{h}_{j}\)):

(22)\[F_{i,d}^{\text{res}} = - \rho_0\,A_{i}\,h^{\text{eff}}_{i}\,C_{\text{res}}\,\|\overline{\mathbf{u}}_{i}\|\,\overline{u}_{i,d}\]

The term related to the entrained mass and mass balance \(- \overline{u_d}\,\rho_0\,\frac{\mathrm{d}(A\,\overline{h})}{\mathrm{d}t}\) now reads:

\[- \overline{u}_{i,d}\,\rho_0\,\frac{\mathrm{d}}{\mathrm{d}t}\,\left(A_{i}\,\overline{h}_{i}\right) = - \overline{u}_{i,d}\,A^{\text{ent}}_{i}\,q^{\text{ent}}_{i}\]

The mass of entrained snow for each particle depends on the type of entrainment involved (ploughing or erosion) and reads:

\[\rho_0\,\frac{\mathrm{d}}{\mathrm{d}t}\,\left(A_{i}\,\overline{h}_{i}\right) = \frac{\mathrm{d}\,m_{i}}{\mathrm{d}t} = A_{i}^\text{ent}\,q_{i}^{\text{ent}}\]

with

\[\begin{split}\begin{aligned} A_{i}^{\text{plo}} &= w_f\,h_{i}^{\text{ent}}= \sqrt{\frac{m_{i}}{\rho_0\,\overline{h}_{i}}}\,h_{i}^{\text{ent}} \quad &\mbox{and} \quad &q_{i}^{\text{plo}} = \rho_{\text{ent}}\,\left\Vert \overline{\mathbf{u}}_{i}\right\Vert \quad &\mbox{for ploughing}\\ A_{i}^{\text{ero}} &= A_{i} = \frac{m_{i}}{\rho_0\,\overline{h}_{i}} \quad &\mbox{and} \quad &q_{i}^{\text{ero}} = \frac{\tau_{i}^{(b)}}{e_b}\,\left\Vert \overline{\mathbf{u}}_{i}\right\Vert \quad &\mbox{for erosion}\end{aligned}\end{split}\]

Finaly, the entrainment force reads:

\[F_{i,d}^{\text{ent}} = -w_f\,(e_s+\,q_{i}^{\text{ent}}\,e_d)\]

Time discretization

The different particle properties are also discretized in time (using a constant time step \(\Delta \,t\)). The supscript \(k\) decribes which time step is considered (for the quantity \(a\), \(a(t_k)=a(k\Delta \,t)\) is refered as \(a^k\)). This leads, for the velocity, flow depth and mass to :

\begin{align} &\overline{u}_{i,d}(t_k) = \overline{u}_{i,d}^k \quad &\text{current time step,}\\ &\overline{u}_{i,d}(t_{k+1}) = \overline{u}_{i,d}^{k+1} \quad &\text{next time step,} \\ &\overline{h}_{i}(t_k) = \overline{h}_{i}^k = \overline{h}_{i} \quad &\text{current time step (supscript dropped for simplicity),}\\ &m_{i}(t_k) = m_{i}^k = m_{i} \quad &\text{current time step (supscript dropped for simplicity),}\\ &m_{i}(t_{k+1}) = m_{i}^{k+1} = m_{i} + \Delta m_{i} \quad &\text{next time step.} \end{align}

Descretizing the momentum balance Eq.14 in time enables to write the the velocity of the particle at the next time step:

\[u_d^{k+1} = \frac{u_d^k + \Delta{t}\,\left(g_d+\frac{F_d+F_d^\text{ent}}{m_p}\right)} {1 + \Delta{t}\left(\frac{\tau^{(b)}}{\rho_0\,\overline{h}\,\|\overline{u}\|^k}+C_\text{res}\,\|\overline{u}\|^k\right)} -u_d^k\,\frac{m_p}{m_p+\Delta{m}_p}\]

With

\[\begin{split}\begin{aligned} k \qquad &\text{current time step,}\\ k+1 \qquad &\text{next time step,}\\ \Delta{t} \qquad &\text{time step size.} \end{aligned}\end{split}\]

The new position of the particle (in the next time step \(k+1\) using “central” velocity) reads:

\[X_d^{k+1} = X_d^k + \frac{\Delta{t}}{2}(u_d^k + u_d^{k+1})\]

Fig. 7 Infinitesimal volume element and acting forces on it (from [FiKo2013])