# Two-fluid barotropic models for powder-snow avalanche flows

Yannick Meyapin, Denys Dutykh, Marguerite Gisclon

► To cite this version:

Yannick Meyapin, Denys Dutykh, Marguerite Gisclon. Two-fluid barotropic models for powder-snow avalanche flows. Computational Science and High Performance Computing IV, Springer Berlin Heidelberg, pp.219-232, 2011, 978-3-642-17770-5. 10.1007/978-3-642-17770-5_17. hal-00394437

HAL Id: hal-00394437

https://hal.science/hal-00394437v1

Submitted on 11 Jun 2009

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
**Chapter 1**

**Two-fluid barotropic models for powder-snow**
**avalanche flows**

Yannick Meyapin, Denys Dutykh and Marguerite Gisclon

**Abstract** In the present study we discuss several modeling issues of powder-snow avalanche flows. We take a two-fluid modeling paradigm. For the sake of simplicity, we will restrict our attention to barotropic equations. We begin the exposition by a compressible model with two velocities for each fluid. However, this model may become non-hyperbolic and thus, represents serious challenges for numerical methods. To overcome these issues, we derive a single velocity model as a result of a relaxation process. This model can be easily shown to be hyperbolic for any reasonable equation of state. Finally, an incompressible limit of this model is derived.

**1.1 Introduction**

Snow avalanches represent a serious problem for society in mountain regions. The avalanche winter of 1999 attracted a lot of attention to this hazardous natural phenomenon [1, 10]. Further development of mountain regions requires an adequate level of avalanche safety. Therefore, avalanche protective measures (deflecting and catching dams) become increasingly important [9]. During the same winter, several avalanches overran avalanche dams, underlining the need for further research in this field. Proper design of protecting structures necessitates profound understanding

Yannick Meyapin

LAMA, UMR 5127 CNRS, Université de Savoie, 73376 Le Bourget-du-Lac Cedex, France, e-mail: Yannick.Meyapin@etu.univ-savoie.fr

Denys Dutykh

LAMA, UMR 5127 CNRS, Université de Savoie, 73376 Le Bourget-du-Lac Cedex, France, e-mail: Denys.Dutykh@univ-savoie.fr

Marguerite Gisclon

LAMA, UMR 5127 CNRS, Université de Savoie, 73376 Le Bourget-du-Lac Cedex, France, e-mail: Marguerite.Gisclon@univ-savoie.fr
of the snow avalanches flow and of the interaction process with dams and other obstacles [6, 13].

Natural snow avalanches are believed to consist of three different layers: a dense core, a fluidised layer and a suspension cloud. Sometimes the surrounding powder cloud is absent and we speak about an avalanche in the flowing régime. Obviously, transition boundaries between these layers are not sharp and this classification is rather conventional.

The dense core consists of snow particles in persistent frictional contact [8]. The density is of the order of 300 kg/m³ and the depth of this layer does not exceed 3 m. The fluidised régime is characterized by particle's mean-free-paths up to several particle's diameters. This dynamics at microscopic level explains more fluid-like behaviour at large scales. The density of this layer is in the range of 50 - 100 kg/m³ and the height is about 3 - 5 m. To model successfully this kind of flows it is crucial to know the complex fluid rheology. Finally, these two interior layers can be covered by the powder cloud which is a turbulent suspension of snow particles in the air. The density ranges from 4 to 20 kg/m³ and an avalanche in aerosol régime can reach the height of 100 m or more [14]. This flow is driven essentially by turbulent advection and particles collisions are unimportant.

In the present study we are concerned with some questions of powder-snow avalanche modelling. Since the interface cannot be defined for this type of flows, we choose the modelling paradigm of two-phase flows. In this approach the governing equations of each phase are spatially averaged to come up with the description of the fluid mixture [7, 15].

It is known [14] that the front of such an avalanche can develop the speed<fnref n="1" /> $u_f \approx 100$ m/s. For comparison, the speed of sound $c_0$ in the air is about 300 m/s. It means that the local Mach number $Ma$ can reach the value of

$$Ma := \frac{u_f}{c_0} \approx 0.33.$$

Hence, compressible effects may become important. That is why, we begin our exposition with a compressible model. Then, we gradually simplify it to come up with an incompressible one at the end of the present article. The goal is achieved by taking the limit as the Mach number tends to zero.

The present article is organized as follows. In Section 1.2 we present a barotropic compressible two-phase model with two velocities. Then, this model is simplified in Section 1.3 using a velocity relaxation process. The incompressible limit of resulting system is derived in Section 1.4. Finally, several conclusions and perspectives are drawn out in Section 1.5.

<footnote n="1">¹ When we estimate the Mach number magnitude, the particle characteristic velocity should be taken. However, this information is not easily accessible and we took the maximum front velocity. It can lead to some overestimation of the Mach number.</footnote>
## 1.2 Two-phase flow modelling

Let us consider a domain $\Omega \subseteq \mathbb{R}^3$ where a simultaneous flow of two barotropic fluids occurs. All quantities related to the heavy and light fluids will be denoted by + and – correspondingly. In view of application to snow avalanches, one can consider the heavy fluid of being constituted of snow particles and the light fluid is the air. When the mixing process is extremely complicated and it is impossible to follow the interface between two fluids, the classical modelling procedure consists in applying a volume average operator [7, 15]. Thereby, we make appear two additional variables $\alpha^{\pm}(\mathbf{x}, t)$, $\mathbf{x} \in \Omega$ which are called the volume fractions and defined as:

$$ \alpha^{\pm}(\mathbf{x}, t) := \lim_{\substack{|d\Omega| \to 0 \\ \mathbf{x} \in d\Omega}} \frac{|d\Omega^{\pm}|}{|d\Omega|}, $$

the heavy fluid occupies volume $d\Omega^{+} \subseteq d\Omega$ and the light one the volume $d\Omega^{-} \subseteq d\Omega$ (see Figure 1.1) such that

$$ |d\Omega| = |d\Omega^{+}| + |d\Omega^{-}|. \qquad (1.1) $$

From the relation (1.1) it is obvious that $\alpha^{+}(\mathbf{x}, t) + \alpha^{-}(\mathbf{x}, t) \equiv 1, \forall \mathbf{x} \in \Omega$.

After performing the averaging process, one obtains two equations of mass and momentum conservation:

$$ \partial_t (\alpha^{\pm} \rho^{\pm}) + \nabla \cdot (\alpha^{\pm} \rho^{\pm} \mathbf{u}^{\pm}) = 0, \qquad (1.2) $$

$$ \partial_t (\alpha^{\pm} \rho^{\pm} \mathbf{u}^{\pm}) + \nabla \cdot (\alpha^{\pm} \rho^{\pm} \mathbf{u}^{\pm} \otimes \mathbf{u}^{\pm}) + \alpha^{\pm} \nabla p = \nabla \cdot (\alpha^{\pm} \tau^{\pm}) + \alpha^{\pm} \rho^{\pm} \mathbf{g}, \quad (1.3) $$

where $\rho^{\pm}(\mathbf{x}, t)$, $\mathbf{u}^{\pm}(\mathbf{x}, t)$, $\tau^{\pm}(\mathbf{x}, t)$ are densities, velocities and viscous stress tensors of each fluid respectively. Traditionally, the vector $\mathbf{g}$ denotes the gravity acceleration. We assume that both fluids share the same pressure<fnref n="2" /> $p = p^{\pm}(\rho^{\pm})$ and equations of state of each phase fulfill minimal thermodynamical requirements:

<figure><img src="image_2.png" /><figcaption>Fig. 1.1 An elementary fluid volume $d\Omega$ occupied by two phases.</figcaption></figure>

<footnote n="2">² In general, this kind of assumptions is reasonable, since relaxation processes will tend to equilibrate the system when time evolves.</footnote>
$$p^{\pm}(\rho^{\pm}) > 0, \quad \frac{\partial p^{\pm}(\rho^{\pm})}{\partial \rho^{\pm}} > 0, \quad \text{for } \rho^{\pm} > 0. \qquad (1.4)$$

In order to obtain a well-posed problem, governing equations (1.2), (1.3) should be completed by appropriate initial and boundary conditions.

If we assume both fluids to be Newtonian, the viscous stress tensor $\tau^{\pm}$ takes the following classical form:

$$\tau^{\pm} = \lambda^{\pm} \mathrm{tr} \mathbb{D}(\mathbf{u}^{\pm}) \mathbf{I} + 2\mu^{\pm} \mathbb{D}(\mathbf{u}^{\pm}), \quad \mathrm{tr} \mathbb{D}(\mathbf{u}^{\pm}) = \nabla \cdot \mathbf{u}^{\pm}, \qquad (1.5)$$

where $\mathbf{I} := (\delta_{ij})_{1\le i,j\le 3}$ is the identity tensor, $\mathbb{D}(\mathbf{u}) := \frac{1}{2}(\nabla\mathbf{u}^{\pm} + \mathbf{t}(\nabla\mathbf{u}^{\pm}))$ is the deformation rate and $\lambda^{\pm}, \mu^{\pm}$ are viscosity coefficients. For ideal gases, for example, these coefficients are related by Stokes relation $\lambda^{\pm} + \frac{2}{3}\mu^{\pm} = 0$. In application to powder-snow avalanches, viscosity coefficients $\lambda^{\pm}, \mu^{\pm}$ should be understood in the sense of eddy viscosity.

**Remark 1.** From physical point of view, presented here model (1.2), (1.3) is far from being complete. For example, one could supplement it by capillarity effects in the Korteweg form. Also we omitted all the terms which model mass, momentum and energy exchange between two phases. Generally, their form is strongly dependent on the physical situation under consideration.

**Remark 2.** Since we do not consider the total energy conservation equation, the fluids are implicitly assumed to be barotropic. In the absence of viscous stresses $\tau^{\pm}$, the flow is isentropic. This simplification can be adopted provided that important energy transfers do not occur. Non-isentropic flows are considered in [11].

**Remark 3.** While considering two-phase flows, it is useful to introduce several additional quantities which play an important rôle in the description of such flows. The mixture density $\rho$ and mass fractions $m^{\pm}$ are naturally defined as:

$$\rho(\mathbf{x},t) := \alpha^{+} \rho^{+} + \alpha^{-} \rho^{-} > 0, \quad \forall (\mathbf{x},t) \in \Omega \times [0, T],$$

$$m^{\pm} := \frac{\alpha^{\pm} \rho^{\pm}}{\rho}, \quad m^{+} + m^{-} = 1.$$

The total density $\rho$ is assumed to be strictly positive everywhere in the domain $\Omega$. Hence, the void creation is forbidden in our modeling.

Important quantities $\rho, m^{\pm}$ will appear several times below.

In principle, one could use equations (1.2), (1.3) to model various two-phase flows. However, this system remains quite expensive for large scale simulations required by real-life applications. The major difficulty comes from the advection operator associated to model (1.2), (1.3) which can be non-hyperbolic [3, 15]. In the next section we will derive a simplified two-fluid model which is proposed as a candidate for powder-snow avalanche compressible simulations.
## 1.3 Velocity relaxation

We would like to reduce the number of variables in the system (1.2), (1.3). The main idea is to introduce the common velocity field for both phases. For this purpose, we will introduce a relaxation term to the momentum conservation equation (1.3):

$$ \partial_t (\alpha^{\pm} \rho^{\pm} \mathbf{u}^{\pm}) + \nabla \cdot (\alpha^{\pm} \rho^{\pm} \mathbf{u}^{\pm} \otimes \mathbf{u}^{\pm}) + \alpha^{\pm} \nabla p = \nabla \cdot (\alpha^{\pm} \tau^{\pm}) + \alpha^{\pm} \rho^{\pm} \mathbf{g}^{\pm} \pm \frac{\kappa}{\epsilon} (\mathbf{u}^{\pm} - \mathbf{u}^{-}), \quad (1.6) $$

where $\kappa = O(1)$ is a constant and $\epsilon$ is a small parameter which controls the magnitude of the relaxation term. Physically this additional term represents the friction between two phases. In the following, we are going to take the singular limit as the relaxation parameter $\epsilon \to 0$. This is achieved with Chapman-Enskog type expansion. In this way, we constrain velocities $\mathbf{u}^{\pm}(\mathbf{x}, t)$ to tend to the common value $\mathbf{u}(\mathbf{x}, t)$. This technique has been already successfully applied to the Baer-Nunziato model [2] in [12].

The first step consists in rewriting the governing equations (1.2), (1.6) in the quasilinear form. To shorten notations, we will also use the material time derivative which is classically defined for any smooth scalar function $\varphi(\mathbf{x}, t)$ as

$$ \frac{d^{\pm} \varphi}{dt} := \frac{\partial \varphi}{\partial t} + \mathbf{u}^{\pm} \cdot \nabla \varphi. $$

**Lemma 1.** Smooth solutions to equations (1.2), (1.6) satisfy the following system:

$$ \alpha^{\pm} \frac{d^{\pm} p}{dt} + \rho^{\pm} (c_s^{\pm})^2 \frac{d^{\pm} \alpha^{\pm}}{dt} + \alpha^{\pm} \rho^{\pm} (c_s^{\pm})^2 \nabla \cdot \mathbf{u}^{\pm} = 0, \quad (1.7) $$

$$ \alpha^{\pm} \rho^{\pm} \frac{d^{\pm} \mathbf{u}^{\pm}}{dt} + \alpha^{\pm} \nabla p = \nabla \cdot (\alpha^{\pm} \tau^{\pm}) + \alpha^{\pm} \rho^{\pm} \mathbf{g}^{\pm} \pm \frac{\kappa}{\epsilon} (\mathbf{u}^{\pm} - \mathbf{u}^{-}), \quad (1.8) $$

where $(c_s^{\pm})^2 := \frac{\partial p^{\pm}}{\partial \rho^{\pm}}|_{s^{\pm}}$ represents the sound speed in each phase $\pm$.

*Proof.* This result follows from direct calculations. First of all, we remark that the mass conservation equation (1.2) can be rewritten using the material derivative as follows:

$$ \frac{d^{\pm}(\alpha^{\pm} \rho^{\pm})}{dt} + \alpha^{\pm} \rho^{\pm} \nabla \cdot \mathbf{u}^{\pm} = 0. \quad (1.9) $$

Using equations of state $p = p^{\pm}(\rho^{\pm})$, we can express the density material derivative in terms of the pressure and the sound speed:

$$ \frac{d^{\pm} \rho^{\pm}}{dt} = \frac{1}{(c_s^{\pm})^2} \frac{d^{\pm} p}{dt}. $$

Now, it is straightforward to derive equation (1.7) from (1.9).

Finally, if we multiply equation (1.9) by $\mathbf{u}^{\pm}$ and subtract it from the momentum conservation equation (1.6), we will get desired result (1.8).
Equations (1.7), (1.8) can be also recast in the matrix form which is particularly useful for further developments:

$$A(V_{\epsilon}) \frac{\partial V_{\epsilon}}{\partial t} + B(V_{\epsilon}) \nabla V_{\epsilon} = \nabla \cdot T(V_{\epsilon}) + S(V_{\epsilon}) + \frac{R(V_{\epsilon})}{\epsilon}, \quad (1.10)$$

where we introduced several notations. The vector $V_{\epsilon}$ represents four unknown physical variables $V_{\epsilon} := {}^t(p, \alpha^{+}, \mathbf{u}^{+}, \mathbf{u}^{-})$ and $\frac{\partial V_{\epsilon}}{\partial t} := {}^t(\partial_t p, \partial_t \alpha^{+}, \partial_t \mathbf{u}^{+}, \partial_t \mathbf{u}^{-})$ and $\nabla V_{\epsilon} := {}^t(\nabla p, \nabla \alpha^{+}, (\cdot \nabla)\mathbf{u}^{+}, (\cdot \nabla)\mathbf{u}^{-})$. Matrices $A(V_{\epsilon})$ and $B(V_{\epsilon})$ are defined as

$$A(V_{\epsilon}) := \begin{pmatrix} \alpha^{+} & \rho^{+}(c_{s}^{+})^{2} & 0 & 0 \\ \alpha^{-} & -\rho^{-}(c_{s}^{-})^{2} & 0 & 0 \\ 0 & 0 & \alpha^{+}\rho^{+}\mathbf{I} & 0 \\ 0 & 0 & 0 & \alpha^{-}\rho^{-}\mathbf{I} \end{pmatrix},$$

$$B(V_{\epsilon}) := \begin{pmatrix} \alpha^{+}\mathbf{u}^{+} & \rho^{+}(c_{s}^{+})^{2}\mathbf{u}^{+} & \alpha^{+}\rho^{+}(c_{s}^{+})^{2}\mathbf{I} & 0 \\ \alpha^{-}\mathbf{u}^{-} & -\rho^{-}(c_{s}^{-})^{2}\mathbf{u}^{-} & 0 & \alpha^{-}\rho^{-}(c_{s}^{-})^{2}\mathbf{I} \\ \alpha^{+}\mathbf{I} & 0 & \alpha^{+}\rho^{+}\mathbf{u}^{+} & 0 \\ 0 & \alpha^{-}\mathbf{I} & 0 & \alpha^{-}\rho^{-}\mathbf{u}^{-} \end{pmatrix}.$$

In these matrix notations the size of zero entries must be chosen to make the multiplication operation possible.

On the right hand side of (1.10), the work of viscous forces is denoted by symbol $\nabla \cdot T(V_{\epsilon}) := {}^t(0, 0, \nabla \cdot \tau^{+}, \nabla \cdot \tau^{-})$. The source term $S(V_{\epsilon}) := {}^t(0, 0, \alpha^{+}\rho^{+}\mathbf{g}, \alpha^{-}\rho^{-}\mathbf{g})$ incorporates the gravity force and $R(V_{\epsilon}) := {}^t(0, 0, \kappa(\mathbf{u}^{+} - \mathbf{u}^{-}), -\kappa(\mathbf{u}^{+} - \mathbf{u}^{-}))$ contains the relaxation terms.

Since we expect the limit $V_{\epsilon} \to V$ to be finite as $\epsilon \to 0$, necessary the limiting vector $V$ lies in the hypersurface $R(V) = 0$. In terms of physical variables, it implies $\mathbf{u}^{+} \equiv \mathbf{u}^{-}$. Consequently, we find our solution in the form of the following Chapman-Enskog type expansion:

$$V_{\epsilon} = V + \epsilon W + O(\epsilon^2).$$

After substituting this expansion into (1.10) and taking into account that $R(V) \equiv 0$, at the leading order in $\epsilon$ one obtains:

$$A(V) \frac{\partial V}{\partial t} + B(V) \nabla V = \nabla \cdot T(V) + S(V) + R'(V)W, \quad (1.11)$$

where

$$R'(V) := \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & \kappa \mathbf{I} & -\kappa \mathbf{I} \\ 0 & 0 & -\kappa \mathbf{I} & \kappa \mathbf{I} \end{pmatrix}$$

Henceforth, we make a technical assumption of the presence of both phases in any point $\mathbf{x} \in \Omega$ of the flow domain. Mathematically it means that $0 < \alpha^{+} < 1$. Since $\alpha^{+} + \alpha^{-} = 1$, the same inequality holds for $\alpha^{-}$. Otherwise, the relaxation
process physically does not make sense and we will have some mathematical technical difficulties.

Under the aforementioned assumption, the matrix $A(V)$ is invertible. Hence, we can multiply on the left both sides of (1.11) by $PA^{-1}(V)$ where the projection matrix $P$ is to be specified below:

$$ \mathbf{P} \frac{\partial V}{\partial t} + \mathbf{PA}^{-1}(V)\mathbf{B}(V)\nabla V = \mathbf{PA}^{-1}(V)\nabla \cdot \mathbf{T}(V) + \mathbf{P}\tilde{\mathbf{R}}'(V)\mathbf{W} + \mathbf{PA}^{-1}(V)\mathbf{S}(V), \quad (1.12) $$

where $\tilde{\mathbf{R}}'(V) := \mathbf{A}^{-1}(V)\mathbf{R}'(V)$ and has the following components

$$ \tilde{\mathbf{R}}'(V) = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{\mathbf{\kappa}}{\alpha^+ \rho^+} \mathbf{I} & -\frac{\mathbf{\kappa}}{\alpha^+ \rho^+} \mathbf{I} \\ 0 & 0 & -\frac{\mathbf{\kappa}}{\alpha^- \rho^-} \mathbf{I} & \frac{\mathbf{\kappa}}{\alpha^- \rho^-} \mathbf{I} \end{pmatrix}. $$

The vector of physical variables $V$ has four (in 1D) components $^t(p, \alpha^+, \mathbf{u}, \mathbf{u})$ and only three are different. In order to remove the redundant information, we will introduce the new vector $U$ defined as $U := {}^t(p, \alpha^+, \mathbf{u})$. The Jacobian matrix of this transformation can be easily computed:

$$ \mathbf{J} := \frac{\partial V}{\partial U} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \mathbf{I} \\ 0 & 0 & \mathbf{I} \end{pmatrix}. $$

In new variables equation (1.12) becomes:

$$ \mathbf{P}\mathbf{J}\frac{\partial U}{\partial t} + \mathbf{PA}^{-1}(U)\mathbf{B}(U)\mathbf{J}\nabla U = \mathbf{PA}^{-1}(U)\nabla \cdot \mathbf{T}(U) + \mathbf{P}\tilde{\mathbf{R}}'(U)\mathbf{W} + \mathbf{PA}^{-1}(U)\mathbf{S}(U). \quad (1.13) $$

Now we can formulate two conditions to construct the matrix $\mathbf{P}$. First of all, the vector $\mathbf{W}$ is unknown and we need to remove it from equation (1.13). Hence, we require $\mathbf{P}\tilde{\mathbf{R}}'(V) = 0$. Then, we would like the governing equations to be explicitly resolved with respect to time derivatives. It gives us the second condition $\mathbf{P}\mathbf{J} = \mathbf{I}$. The existence and effective construction of the matrix $\mathbf{P}$ satisfying two aforementioned conditions

$$ \mathbf{P}\tilde{\mathbf{R}}'(V) = 0, \quad \mathbf{P}\mathbf{J} = \mathbf{I}, $$

are discussed below. Presented in this section results follow in great lines [12].

We will consider a slightly more general setting. Let vector $V \in \mathbb{R}^n$ and its reduced counterpart $U \in \mathbb{R}^{n-k}$, $k < n$. In such geometry, $\tilde{\mathbf{R}}'(V) \in \mathrm{Mat}_{n,n}(\mathbb{R})$, $\mathbf{J} \in \mathrm{Mat}_{n,n-k}(\mathbb{R})$ and, consequently, $\mathbf{P} \in \mathrm{Mat}_{n-k,n}(\mathbb{R})$. Here, the notation $\mathrm{Mat}_{m,n}(\mathbb{R})$ denotes the set of $m \times n$ matrices with coefficients in $\mathbb{R}$. We have to say also that from algebraic point of view, matrices $\tilde{\mathbf{R}}'(V)$ and $\mathbf{R}'(V)$ are completely equivalent. Thus, for simplicity, in the following propositions we will reason in terms of $\mathbf{R}'(V)$.
**Lemma 2.** *The columns of the Jacobian matrix J form a basis of ker($R'(V)$).*

*Proof.* If we differentiate the relation $R(V) = 0$ with respect to $U$, we will get the identity $R'(V)J = 0$. It implies that $\text{range}(J) \subseteq \text{ker}(R'(V))$. By direct computation one verifies that $\text{dim range}(R'(V)) = k$. From the well-known identity $\text{range}(R'(V)) \oplus \text{ker}(R'(V)) = \mathbb{R}^n$, one concludes that $\text{dim ker}(R'(V)) = n-k$. But in the same time, the rank of J is equal to $n-k$ as well. It proves the result.

**Theorem 1.** *We suppose that for all V, range($R'(V)$) $\cap$ ker($R'(V)$) = {0} then it exists a matrix P $\in$ Mat$_{n-k,n}$(R) such that PR'(V) = 0 and PJ = I$_{n-k}$.*

*Proof.* Hypothesis range($R'(V)$) $\cap$ ker($R'(V)$) = {0} implies that range($R'(V)$) $\oplus$ ker($R'(V)$) = $\mathbb{R}^n$. From Lemma 2 it follows that range($J$) = $\text{ker}(R'(V))$. Thus, the space $\mathbb{R}^n$ can be also represented as a direct sum range($R'(V)$) $\oplus$ range($J$). We will define $\mathbf{P}$ to be the projection on $\text{ker}(R'(V)) \equiv \text{range}(J)$. Since obviously $R'(V) \in \text{range}(R'(V))$ and $J \in \text{range}(J)$, we have two required identities: $\mathbf{P}J = I_{n-k}$ and $PR'(V) = 0$.

Now, in order to compute effectively the projection matrix $\mathbf{P}$, we will construct an auxiliary matrix $\mathbf{D}(V) = [J^1, \dots, J^{n-k}, I^1, \dots, I^k]$, where $J^i$ is the column $i$ of the matrix $\mathbf{J}$ and $\{I^1, \dots, I^k\}$ are vectors which form a basis of range($R'(V)$). We remark that $\mathbf{PD}(V) = [I_{n-k}, 0]$. Lemma 2 implies that the matrix $\mathbf{D}(V)$ is invertible. Thus, the projection $\mathbf{P}$ can be computed by inverting $\mathbf{D}(V)$:

$$ \mathbf{P} = [I_{n-k}, 0] \cdot \mathbf{D}^{-1}(V). $$

Let us apply this general framework to our model (1.12), where $n = 4$ and $k = 1$. The matrix $\mathbf{D}(V)$ and its inverse $\mathbf{D}^{-1}(V)$ take this form:

$$ \mathbf{D}(V) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \mathbf{I} & \frac{\mathbf{I}}{\alpha^{+} \rho^{+}} \\ 0 & 0 & \mathbf{I} & -\frac{\mathbf{I}}{\alpha^{-} \rho^{-}} \end{pmatrix}, \quad \mathbf{D}^{-1}(V) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \mathbf{I}^{+} & \mathbf{I}^{-} \\ 0 & 0 & \frac{\mathbf{I}^{+} \mathbf{I}^{-} \rho}{\kappa} & \frac{\mathbf{I}^{+} \mathbf{I}^{-} \rho}{\kappa} \end{pmatrix}, $$

where $\mathbf{m}^{\pm}$ are mass fractions defined in Remark 3.

Now, the projection matrix $\mathbf{P}$ can be immediately computed:

$$ \mathbf{P} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \mathbf{m}^{+} \mathbf{I} & \mathbf{m}^{-} \mathbf{I} \end{pmatrix}. $$

Finally, after computing all matrix products $\mathbf{PA}^{-1}(U)\mathbf{B}(U)\mathbf{J}$, $\mathbf{PA}^{-1}(U)\nabla \cdot \mathbf{T}(U)$, $\mathbf{PA}^{-1}(U)\mathbf{S}(U)$ present in equation (1.13), we obtain the desired single velocity model:
$$ \frac{\partial p}{\partial t} + \mathbf{u} \cdot \nabla p + \rho c_s^2 \nabla \cdot \mathbf{u} = 0, \qquad (1.14) $$

$$ \frac{\partial \alpha^+}{\partial t} + \mathbf{u} \cdot \nabla \alpha^+ + \alpha^+ \alpha^- \delta \nabla \cdot \mathbf{u} = 0, \qquad (1.15) $$

$$ \rho \frac{\partial \mathbf{u}}{\partial t} + \rho(\mathbf{u} \cdot \nabla)\mathbf{u} = \rho \mathbf{g} + \nabla \cdot \boldsymbol{\tau}, \qquad (1.16) $$

where $\rho = \alpha^+\rho^+ + \alpha^-\rho^-$ is the mixture density and $c_s^2$ is the sound velocity in the mixture which is determined by this formula:

$$ \rho c_s^2 := \frac{\rho^+ \rho^- (c_s^+)^2 (c_s^-)^2}{\alpha^- \rho^+ (c_s^+)^2 + \alpha^+ \rho^- (c_s^-)^2}, $$

and $\delta$ is given by

$$ \delta := \frac{\rho^+(c_s^+)^2 - \rho^-(c_s^-)^2}{\alpha^- \rho^+(c_s^+)^2 + \alpha^+ \rho^-(c_s^-)^2}. $$

Finally, $\tau := \lambda \operatorname{tr}\mathbb{D}(\mathbf{u})\mathbf{I} + 2\mu\mathbb{D}(\mathbf{u})$ is the viscous stress tensor of the mixture. Viscosity coefficients $\lambda, \mu$ are naturally defined as

$$ \lambda := \alpha^{+} \lambda^{+} + \alpha^{-} \lambda^{-}, \quad \mu := \alpha^{+} \mu^{+} + \alpha^{-} \mu^{-}. $$

Equations (1.14) – (1.16) can be recast in the conservative form which is more convenient for numerical computations and theoretical analysis. To achieve this purpose, we replace the pressure $p$ in (1.14) by $\rho^\pm$ using the equation of state:

$$ \frac{\partial \rho^{\pm}}{\partial t} + \mathbf{u} \cdot \nabla \rho^{\pm} + \frac{\rho c_s^2}{(c_s^{\pm})^2} \nabla \cdot \mathbf{u} = 0. $$

The last equation is multiplied by $\alpha^\pm$, the second equation (1.15) is multiplied by $\rho^\pm$ and we sum them to come up with two mass conservation equations. Transformation of the momentum conservation equation (1.16) is straightforward. The resulting conservative system takes this form:

$$ \partial_t (\alpha^{\pm} \rho^{\pm}) + \nabla \cdot (\alpha^{\pm} \rho^{\pm} \mathbf{u}) = 0, \qquad (1.17) $$

$$ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla p = \nabla \cdot \boldsymbol{\tau} + \rho \mathbf{g}. \qquad (1.18) $$

These equations represent a barotropic version of the four-equations model proposed in [4, 5].

It can be shown that the advection operator of the model (1.17), (1.18) is hyperbolic for any reasonable equation of state (1.4). Moreover, this system contains fewer variables which allow more efficient computations required in practice.
## 1.4 Incompressible limit

The main scope of this paper is certainly around compressible two-fluid models. However, we decided to derive an incompressible limit of the single velocity model (1.17), (1.18) for the case when acoustic effects should be filtered out. The presence of acoustic waves represent, for example, a major restriction for the time step, if an explicit scheme is used.

For the sake of simplicity, we will neglect dissipative effects which do not affect the acoustic wave propagation. Thus, in this section we consider the following system of equations:

$$ \partial_t (\alpha^{\pm} \rho^{\pm}) + \nabla \cdot (\alpha^{\pm} \rho^{\pm} \mathbf{u}) = 0, \qquad (1.19) $$

$$ \rho \partial_t \mathbf{u} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p = \rho \mathbf{g}. \qquad (1.20) $$

$$ (1.21) $$

For convenience, we rewrite equation (1.18) in nonconservative form.

In order to estimate the relative importance of various terms, we introduce dimensionless variables. The characteristic length, time, and velocity scales are denoted by $l$, $t_0$ and $U_0$ respectively. For example, $l$ may be chosen as the diameter of the fluid domain $\Omega$, $t_0$ is the biggest vortex turnover time and $U_0$ is the typical flow velocity. The density and the sound velocity scales are chosen to be those of the heavy fluid, i.e. $\rho_0^+$ and $c_{0s}^+$ correspondingly. Since we are interested in acoustic effects, the natural pressure scale is given by $\rho_0^+ (c_{0s}^+)^2$. If we summarize these remarks, dependent and independent dimensionless variables (denoted with primes) are defined as:

$$ \mathbf{x}' := \frac{\mathbf{x}}{l}, \quad t' := \frac{t}{t_0}, \quad \mathbf{u}' := \frac{\mathbf{u}}{U_0}, \quad (\rho^{\pm})' := \frac{\rho^{\pm}}{\rho_0^{\pm}}, \quad p' := \frac{p}{\rho_0^{+} (c_{0s}^{+})^2}. $$

*Remark 4.* There is nothing to do for the volume fractions $\alpha^{\pm}$, since this quantity is dimensionless by definition.

After dropping the tildes, nondimensional system of equation becomes:

$$ \mathrm{St} \partial_t (\alpha^{\pm} \rho^{\pm}) + \nabla \cdot (\alpha^{\pm} \rho^{\pm} \mathbf{u}) = 0, \qquad (1.22) $$

$$ \mathrm{St} \rho \partial_t \mathbf{u} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} + \frac{1}{\mathrm{Ma}^2} \nabla p = \frac{1}{\mathrm{Fr}^2} \rho \mathbf{g}, \qquad (1.23) $$

where several scaling parameters have appeared:

* Strouhal number $\mathrm{St} := \frac{l}{U_0 t_0}$. In this study we will assume the Strouhal number to be equal to one, i.e. $t_0 = \frac{l}{U_0}$.

* Mach number $\mathrm{Ma} := \frac{U_0}{c_{0s}^+}$ which measures the relative importance of the flow speed and the sound speed in the medium.
* Froude number $Fr := \frac{U_0}{\sqrt{gL}}$ compares inertia and gravitational force. This parameter will not play an important rôle in the present study.

All physical variables $\alpha^{\pm}$, $\rho^{\pm}$, $p$ and $\mathbf{u}$ are expanded in formal series in powers of the Mach number:

$$ \varphi = \varphi_0 + \text{Ma}\varphi_1 + \text{Ma}^2\varphi_2 + \dots, \quad \varphi \in \{\alpha^{\pm}, \rho^{\pm}, p, \mathbf{u}\}. \quad (1.24) $$

Formal expansion (1.24) is then substituted into the system (1.22), (1.23). At the orders $\text{Ma}^{-2}$ and $\text{Ma}^{-1}$, we obtain

$$ \nabla p_0 = \nabla p_1 = 0. $$

In other words, $p_0 = p_0(t)$ and $p_1 = p_1(t)$ are only functions of time. At the order $\text{Ma}^0$ we get the following system of equations:

$$ \partial_t(\alpha_0^{\pm}\rho_0^{\pm}) + \nabla \cdot (\alpha_0^{\pm}\rho_0^{\pm}\mathbf{u}_0) = 0, \quad (1.25) $$

$$ \rho_0 \partial_t \mathbf{u}_0 + \rho_0 (\mathbf{u}_0 \cdot \nabla) \mathbf{u}_0 + \nabla \pi = \frac{1}{\text{Fr}^2} \rho_0 \mathbf{g}, \quad (1.26) $$

where by $\pi$ we denote $p_2$.

Using the same asymptotic expansion (1.24), one can show that at the leading order we keep usual relations between densities and volume fractions:

$$ \alpha_0^{\pm} + \alpha_0^{\mp} = 1, \quad \rho_0 = \alpha_0^{\pm}\rho_0^{\pm} + \alpha_0^{\mp}\rho_0^{\mp}. \quad (1.28) $$

In order to investigate the behaviour of $\rho_0^{\pm}$, we will invert the equation of state<fnref n="3" /> $\rho^{\pm} = \rho^{\pm}(p) = (p^{\pm})^{-1}(p)$ and expand it in powers of $\text{Ma}$:

$$ \rho^{\pm}(p) = \rho^{\pm}(p_0) + \text{Ma} \left. \frac{\partial \rho^{\pm}}{\partial p} \right|_{p_0} p_1 + \text{Ma}^2 \left( \left. \frac{\partial \rho^{\pm}}{\partial p} \right|_{p_0} p_2 + \left. \frac{\partial^2 \rho^{\pm}}{\partial p^2} \right|_{p_0} p_1^2 \right) + \text{O}(\text{Ma}^3) $$

On the other hand, from (1.24) we know that

$$ \rho^{\pm} = \rho_{0}^{\pm} + \text{Ma}\rho_{1}^{\pm} + \text{Ma}^{2}\rho_{2}^{\pm} + \dots $$

Matching these expansions at two lowest orders shows that $\rho_{0,1}^{\pm}$ are functions only of the time variable:

$$ \rho_0^{\pm} = \rho^{\pm}(p_0(t)) :=: r_0^{\pm}(t), \quad \rho_1^{\pm} = \left. \frac{\partial \rho^{\pm}}{\partial p} \right|_{p_0(t)} p_1(t) :=: r_1^{\pm}(t). $$

It is possible to show that $\rho_{0,1}^{\pm}$ are just constants. Consider the Gibbs relation which reads

<footnote n="3">³ The function $p = p^{\pm}(\rho^{\pm})$ is invertible since it is a strictly increasing function $\frac{\partial p}{\partial \rho^{\pm}} > 0$.</footnote>
$$T^{\pm}ds^{\pm} = de^{\pm} - \frac{p}{(\rho^{\pm})^2}d\rho^{\pm}.$$

Since we consider isentropic flows, $ds^{\pm} \equiv 0$ and, consequently, the Gibbs relation takes a much simpler form:

$$de^{\pm} = \frac{p}{(\rho^{\pm})^2}d\rho^{\pm}. \qquad (1.29)$$

It can be shown by considering the total energy conservation equation [11], that the internal energy $e^{\pm}$ naturally scales with $U_0^2$. After dividing (1.29) by $dt$ and switching to dimensionless variables, equation (1.29) takes the following form (after dropping the primes):

$$\frac{de^{\pm}}{dt} = \frac{p}{Ma^2(\rho^{\pm})^2} \frac{d\rho^{\pm}}{dt}.$$

Expanding $e^{\pm}$ in the series (1.24) and looking at two leading terms, leads to the desired result:

$$\frac{d\rho_{0,1}^{\pm}}{dt} = 0 \quad \Rightarrow \quad \rho_{0,1}^{\pm} = \text{const.}$$

The incompressibility condition $\nabla \cdot \mathbf{u}_0 = 0$ is obtained by summing up mass conservation equations (1.25) and taking into account relation (1.28).

If we summarize all developments made above and switch back to dimensional variables, the resulting incompressible system will become:

$$\partial_t \alpha^{\pm} + \nabla \alpha^{\pm} \cdot \mathbf{u} = 0, \qquad (1.30)$$

$$\nabla \cdot \mathbf{u} = 0, \qquad (1.31)$$

$$\rho \partial_t \mathbf{u} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla \pi = \rho \mathbf{g} + \nabla \cdot \tau, \qquad (1.32)$$

where we dropped the index 0 and added again dissipative effects. Viscous stress tensor $\tau$ is still defined by expression (1.5), as in compressible case. In this case, we can speak about two-fluid Navier-Stokes equations. This system of equations (1.30) – (1.32) is much easier to solve numerically than its compressible analogue (1.17), (1.18). In particular, this simplification is due to removed stiffness of acoustic waves.

## 1.5 Conclusions and perspectives

In this study we presented several barotropic two-fluid models which can be used for numerical simulation of powder-snow avalanche flows. One of the main objectives of this paper was to reveal the connection between barotropic models with single and two velocities. The extension to more general fluids is in progress [11].

Our exposition began with compressible two-phase model (1.2), (1.3) possessing two velocity variables. Then, using a relaxation process, we constrained the system to have a common velocity for both phases. Mathematically it was achieved with a Chapman-Enskog type expansion. Resulting model (1.17), (1.18) is hyperbolic for
any reasonable equation of state (1.4). Finally, two-fluid Navier-Stokes equations (1.30) – (1.32) were derived as an incompressible limit of the single velocity model (1.17), (1.18).

Hence, we presented three different two-fluid models which are related by formal derivation procedures. Simplifications made above, represent a good trade-off between accuracy and computational complexity. The final choice should be made after determining the flow régime and main goals of the simulation.

We did not incorporate yet any turbulence modeling. In this study we were focused essentially on the advection operators. However, it is obvious that the physical flow under consideration is fully turbulent in its aerosol part [14]. As the first physical approximation, turbulence effects can be taken into account by adding eddy viscosity terms and, thus, by modifying the viscous stress tensor $\tau$. It will be done in future studies.

## Acknowledgement

The authors would like to acknowledge the University of Savoie for the PPF grant linked to the project: "Mathématiques et avalanches de neige, une rencontre possible?". The support from the Research network VOR (Professors Jacky Mazars and Denis Jongmans) and Cluster Environnement through the program "Risques gravitaires, séismes" is also acknowledged.

We would like to thank Professor Carmen de Jong for interesting discussions around snow avalanches. Special thanks go to our colleagues and friends Didier Bresch and Céline Acary-Robert for their continuous help and support. Finally, the second author thanks Professors Jean-Michel Ghidaglia and Frédéric Dias for introducing him to the beautiful field of two-phase flows.

## References

1. Ancey, C., Bain, V., Bardou, E., Borrel, G., Burnet, R., Jarry, F., Kolbl, O., Meunier, M.: Dynamique des avalanches. Presses polytechniques et universitaires romandes (Lausanne, Suisse) (2006)

2. Baer, M., Nunziato, J.: A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials. International journal of multiphase flow **12**(6), 861–889 (1986)

3. Bresch, D., Desjardins, B., Ghidaglia, J.M., Grenier, E.: On global weak solutions to a generic two-fluid model. To appear in Archive for Rational Mechanics and Analysis (2009)

4. Dias, F., Dutykh, D., Ghidaglia, J.M.: A two-fluid model for violent aerated flows. Submitted to Comput. & Fluids (2008)

5. Dutykh, D.: Mathematical modelling of tsunami waves. Ph.D. thesis, École Normale Supérieure de Cachan (2007)

6. Dutykh, D., Acary-Robert, C., Bresch, D.: Numerical simulation of powder-snow avalanche interaction with an obstacle. Submitted to Applied Mathematical Modelling (2009)

7. Ishii, M.: Thermo-Fluid Dynamic Theory of Two-Phase Flow. Eyrolles, Paris (1975)
8. Issler, D.: Experimental information on the dynamics of dry-snow avalanches. In: K. Hutter, N. Kirchner (eds.) Dynamic Response of Granular and Porous Materials Under Large and Catastrophic Deformations, vol. 11. Springer, Berlin (2003)

9. Johannesson, T., Gauer, P., Issler, P., Lied, K.: The design of avalanche protection dams. Tech. rep., European Commission (2009)

10. Lied, K.: Satsie: Avalanche studies and model validation in europe. Tech. rep., European Commission (2006)

11. Meyapin, Y., Dutykh, D., Gisclon, M.: Velocity and energy relaxation in two-phase flows. Submitted to Eur. J. Mech. B/Fluids (2009)

12. Murrone, A., Guillard, H.: A five equation reduced model for compressible two phase flow problems. J. Comput. Phys. **202**, 664–698 (2005)

13. Naaim-Bouvet, F., Naaim, M., Bacher, M., Heiligenstein, L.: Physical modelling of the inter-action between powder avalanches and defence structures. Nat. Hazards Earth Syst. Sci. **2**, 193–202 (2002)

14. Rastello, M., Hopfinger, E.: Sediment-entraining suspension clouds: a model of powder-snow avalanches. J. Fluid. Mech. **509**, 181–206 (2004)

15. Rovarch, J.M.: Solveurs tridimensionnels pour les écoulements diphasiques avec transferts d'énergie. Ph.D. thesis, Ecole Normale Supérieure de Cachan (2006)