theoryCom1DFA.rst

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 (Zw2000,ZwKlSa2003).

Mass balance:

$$\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^ent represents the snow entrainment rate.

Momentum balance:

$$\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(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^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. ρ_ent < ρ) and F_i^res represents the resistance force due to obstacles (for example trees). This leads to in momentum-balance1:

$$\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 mass-balance1, we get:

$$\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(x, t) = z − s(x, y, t) = 0

The bottom surface is defined by :

F_b(x) = z − b(x, y) = 0

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

$$\begin{aligned} \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{aligned}$$

σ_i^(b) = (σ_kln_ln_k)n_i represents the normal stress at the bottom and τ_i^(b) = σ_ijn_j − σ_i^(b) represents the shear stress at the bottom surface. f describes the chosen friction model and are described in theoryCom1DFA: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, v₁, v₂, v₃), further referenced as Natural Coordinate System (NCS). In this NCS, v₁ is aligned with the velocity vector at the bottom and v₃ with the normal to the slope, i.e.:

$$\mathbf{v_1} = \frac{\mathbf{u}}{\left\Vert \mathbf{u}\right\Vert},\quad \mathbf{v_2} = \mathbf{v_3}\wedge\mathbf{v_1}, \quad \mathbf{v_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$$

Thickness averaged equations:

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 thickness,

$$\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$$

Small Lagrangian prism-like Control volume

Entrainment:

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^plo can be expressed as a function of the properties of the entrained snow (density ρ_ent and snow thickness h_ent), the velocity of the avalanche at the front $\overline{\mathbf{u}}$ and length w_f of the front (measured perpendicularly to the flow velocity $\overline{\mathbf{u}}$). It obviously only happens on the front of the avalanche:

$$\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^ero can be expressed as a function of the bottom area A_b of the control volume, the velocity of the avalanche $\overline{\mathbf{u}}$, the bottom shear stress τ^(b) and the specific erosion energy e_b:

$$\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 mass-balance1 to :

$$\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 + \frac{A_b}{\rho_0}\,\frac{\tau^{(b)}}{e_b}\,\left\Vert \overline{\mathbf{u}}\right\Vert$$

The force F_i^ent required to break the entrained snow from the ground and to compress it is expressed as a 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 thickness (Sa2007,SaFeFr2008,FiFrGaSo2013):

F_i^ent =  − w_f (e_s+ q^ent e_d)

Resistance:

The force F_i^res due to obstacles is expressed as a function of the characteristic diameter $\overline{d}$ and height h_res of the obstacles, the spacing s_res between the obstacles and an empirical coefficient c_w (see f-res). The effective height h^eff is defined as $\min(\overline{h}, h_{res} )$:

$$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}\|}$$

Resistance force due to obstacles (from FiKo2013)

Surface integral forces:

The surface integral is split in three terms, an integral over A_b the bottom x₃ = b(x₁, x₂), A_s the top x₃ = s(x₁, x₂, t) and A_h the lateral surface. Introducing the boundary conditions boundary-conditions leads to:

$$\begin{aligned} \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{aligned}$$

Which simplifies the momentum balance momentum-balance3 to:

$$\begin{aligned} \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{aligned}$$

The momentum balance in direction x₃ (normal to the slope) is used to obtain a relation for the vertical distribution of the stress tensor (Sa2007). Due to the choice of coordinate system and because of the kinematic boundary condition at the bottom, the left side of momentum-balance5 can be expressed as a function of the velocity $\overline{u}_1$ in direction x₁ and the curvature of the terrain in this same direction $\frac{\partial^2{b}}{\partial{x_1^2}}$ (Zw2000):

$$\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:

$$\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

Characteristic size of the avalanche along its path (from Zw2000, modified)

The previous equations momentum-balance5 and sigma33 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:     ε = H/L  
• Curvature:   λ = L/R  

The different properties involved are then expressed in terms of characteristic quantities L, H, g, ρ₀ and R (see fig-characteristic_size):

$$\begin{aligned} \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{aligned}$$

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

σ_ii = ρ₀ g H σ_ii^*

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

$$\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 be small compared to its length, L. Meaning that the equations are examined in the limit ε ≪ 1. It is then possible to neglect the last term in sigma33star which leads to (after reinserting the dimensions):

$$\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₃ = 0, the normal stress can be expressed as:

$$\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 momentum-balance5 requires to express the other components of the stress tensor. Here again a magnitude consideration between the shear stresses σ₁₂ = σ₂₁ and σ₁₃. The shear stresses are based on a generalized Newtonian law of materials, which controls the influence of normal stress 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 ∂x₁ and ∂x₂ are of the order of L, whereas ∂x₃ 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 σ₁₂ = σ₂₁ is negligible compared to σ₁₃. σ₁₃ is expressed using the bottom friction law $\tau^{(b)}_i = f(\sigma^{(b)},\overline{u},\overline{h},\rho_0,t,\mathbf{x})$ introduced in boundary-conditions.

In addition, a relation linking the horizontal normal stresses, σ_ii, i = (1, 2), to the vertical pressure distribution given by sigmab is introduced. In complete analogy to the arguments used by Savage and Hutter (SaHu1989) the horizontal normal stresses are given as:

σ_ii = K_(i) σ₃₃

Where K_(i) are the earth pressure coefficients (cf. ZwKlSa2003,Sa2004):

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

With the above specifications, the integral of the stresses over the flow height is simplified in equation momentum-balance5 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:

$$\begin{aligned} \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{aligned}$$

with

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

The mass balance mass-balance2 remains unchanged:

$$\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 + \frac{A_b}{\rho_0}\,\frac{\tau^{(b)}}{e_b}\,\left\Vert \overline{\mathbf{u}}\right\Vert$$

The unknown $\overline{u}_1$, $\overline{u}_2$ and $\overline{h}$ satisfy sigmab, momentum-balance6 and mass-balance3. In equation momentum-balance6 the bottom shear stress τ^(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 τ^(b) as a function of the flow state of the avalanche.

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

With

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

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

Mohr-Coulomb friction model

The Mohr-Coulomb friction model describes the friction interaction between twos solids. The bottom shear stress simply reads:

τ^(b) = tan δ σ^(b)

tan δ = μ is the friction coefficient (and δ the friction angle). The bottom shear stress linearly increases with the normal stress component σ^(b) (Zw2000,BaSaGr1999,WaHuPu2004,Sa2007).

With this friction model, an avalanche starts to flow if the slope inclination is steeper than the friction angle δ. In the case of an infinite slope of constant inclination, the avalanche velocity would increase indefinitely. This is unrealistic to model snow avalanches because it leads to over prediction of the flow velocity. The Mohr-Coulomb friction model is on the other hand well suited to model granular flow. Because of its relative simplicity, this friction model is also very convenient to derive analytic solutions and validate the numerical implementation.

Chezy friction model

The Chezy friction model describes viscous friction interaction. The bottom shear stress then reads:

τ^(b) = c_dyn ρ₀ ū²

c_dyn is the viscous friction coefficient. The bottom shear stress is a quadratic function of the velocity. (Zw2000,BaSaGr1999,WaHuPu2004,Sa2007).

This model enables to reach more realistic velocities for avalanche simulations. The draw back is that the avalanche doesn't stop flowing before the slope inclination approaches zero. This implies that the avalanche flows to the lowest local point.

Voellmy friction model

Anton Voellmy was a Swiss engineer interested in avalanche dynamics Vo1955. He first had the idea to combine both the Mohr-Coulomb and the Chezy model by summing them up in order to take advantage of both. This leads to the following friction law:

τ^(b) = tan δ σ^(b) + c_dyn ρ₀ ū²

This model is described as Voellmy-Fluid Sa2004,Sa2007, and the turbulent friction term ξ is used instead of c_dyn.

SamosAT friction model

SamosAT friction model is a modification of some more classical models such as Voellmy model Voellmy friction model. The basal shear stress tensor τ^(b) is expressed as (Sa2007):

$$\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{aligned} \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{aligned}$$

The minimum shear stress τ₀ defines a lower limit below which no flow takes place with the condition $\rho_0\,\overline{h}\,g\,\sin{\alpha} &gt; \tau_0$. α being the slope. τ₀ is independent of the flow thickness, 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⁰ 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 α, than without this effect. This effect is intended to avoid lateral creep of the avalanche mass (SaGr2009).

Wet snow friction type

Note

This is an experimental option to account for wet snow conditions, still under development and not yet tested. Also the parameters are not yet calibrated.

In addition, com1DFA provides an optional friction model implementation to account for wet snow conditions. This approach is based on the Voellmy friction model but with an enthalpy dependent friction parameter.

τ^(b) = μ σ^(b) + c_dyn ρ₀ ū²

where,

μ = μ₀ exp ( − enthalpy/enthRef)

The total specific enthalpy of the particles is initialized based on their initial temperature, specific heat capacity, altitude and their velocity (which is zero for the initial time step). Throughout the computation, the particles specific enthalpy is then computed following:

enthalpy = totalEnthalpy − g z − 0.5 ū²

Dam

The dam is described by a crown line, that is to say a series of x, y, z points describing the crown of the dam (the dam wall is located on the left side of the line), by the slope of the dam wall (slope measured from the horizontal, β) and a restitution coefficient (describing if we consider more elastic or inelastic collisions between the particles and the dam wall, varying between 0 and 1).

The geometrical description of the dam is given on the figure fig-DamToolSide. The dam crown line (x_crown) is projected onto the topography, which provides us with the dam center line (x_center). We compute the tangent vector to the center line (t_f). From this tangent vector and the dam slope, it is possible to compute the wall tangent vector (t_w). Knowing the wall tangent vector and height, it is possible to determine normal vector to the wall (n_w) and the foot line which is the intersection between the dam wall and the topography (x_foot).

When the dam fills up (flow thickness increases), the foot line is modified ($\mathbf{x_\text{foot}^\text{filled}} = \mathbf{x_\text{foot}} + \frac{h_v}{2} \mathbf{e_z}$). The normal and tangent vectors to the dam wall are readjusted accordingly.

Side view of the dam (cut view). x_crown describes the crown of the dam, x_center is the vertical projection of the crown on the topography (here the light blue line represents the topography). The tangent vector to the center line (t_f) is computed from the center line points. The tangent vector to the center line with the dam slope angle enable to compute the tangent (t_w) and normal (n_w) vector to the dam wall. Finally, this normal vector is adjusted depending on the snow thickness at the dam location (filling of the dam , n_w^filled)

In the initialization of the simulation, the dam tangent vector to the center line (t_f), foot line (x_foot) and normal vector to the wall (n_w) are computed. The grid cells crossed by the dam as well as their neighbor cells are memorized (tagged as dam cells).