# Averaged models for compressible two-phase stratified flows on thin domains

Nicolas Seguin, Khaled Saleh, Pierrick Le Vourc'H

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

To cite this version:

Nicolas Seguin, Khaled Saleh, Pierrick Le Vourc'H. Averaged models for compressible two-phase stratified flows on thin domains. 2025. hal-05096837

HAL Id: hal-05096837

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

Preprint submitted on 6 Jun 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# AVERAGED MODELS FOR COMPRESSIBLE TWO-PHASE STRATIFIED FLOWS ON THIN DOMAINS*

PIERRICK LE VOURC'H†, KHALED SALEH‡, AND NICOLAS SEGUIN†

**Abstract.** This paper deals with the derivation of compressible two-phase flow models. We use a thin domain approximation of a two-layer configuration governed by the Navier-Stokes equations, following the works [H. B. Stewart and B. Wendroff, J. Comp. Phys., 56 (1984)] and [V. H. Ransom and D. L. Hicks, *J. Comput. Phys.*, 75 (1988)]. In order to recover source terms and two-velocity models directly from this asymptotic analysis, we remove the continuity of the tangential component of the velocity at the interface, and properly scale friction and viscosity coefficients. We are then able to derive two-velocity one-pressure models, first in the barotropic case and then in the full Navier-Stokes-Fourier case.

**Key words.** Compressible fluids, two-phase flows, asymptotic analysis, thin-film approximation, stratified flows

MSC codes. 76T10, 35Q30, 35C20, 76A20

**1. Introduction.** In this paper, we present a derivation of averaged models for compressible two-phase flows. The usual process consists in applying averaging operators on a complete local description of the flow, see for instance the classical textbooks [10] and [14]. The application of this process to bifluid formulations of stratified compressible two-phase flows led to a variety of different averaged models in the 1980s, see for instance [20, 17, 18]. The averaging process being incomplete, the previous papers provide additional heuristic arguments to close the models. They also discuss the influence of the closure arguments on the resulting models, in particular when focusing on interaction terms between phases. The *ad hoc* closure conditions introduced in these papers — especially at the interface — make the derivation of two-velocity models possible. Such approaches have been followed in [9] to include gravity effects as well as mass and momentum transfers.

Recent applications of mathematical homogenization to multiphase compressible Navier-Stokes equations provided the first complete, and rigorous, derivation of one-dimensional averaged models including the mechanical interaction source terms [7, 5, 6, 12, 3, 13]. Indeed, provided that the velocity field is smooth enough — especially through the interfaces between fluids —, this method enables the extraction of averaged variables and equations by letting the frequency of density oscillations tend to infinity. Since this procedure has yet only been carried in one space dimension, the inherent smoothness of the velocity field causes the averaged model to be described by a unique velocity field, making it difficult to obtain two-velocity models.

In this paper, we propose an asymptotic derivation of two-velocity, one-pressure averaged models. They are composed of two equations for phasic mass conservation, two nonconservative equations for phasic momentum, and two nonconservative

<footnote>*Submitted to the editors on June 6, 2025.</footnote>

<footnote>**Funding:** N. Seguin's work was co-funded by the Agence Nationale de la Recherche (ANR) under the ANR MSMPhi project (ANR-23-CE40-0013) and the ANR HEAD project (ANR-24-CE40-3260).</footnote>

<footnote>†IMAG, ANGUS, Univ Montpellier, Inria, CNRS, Montpellier, France (pier-rick.le-vourc-h@umontpellier.fr, nicolas.seguin@inria.fr).</footnote>

<footnote>‡Institut de Mathématiques de Marseille, Université d'Aix-Marseille, CNRS, Marseille, France (khaled.seleh@univ-amu.fr).</footnote>
equations for phasic energies:

$$
\begin{equation} \tag{1.1}
\begin{aligned}
& \partial_t (\alpha_k \rho_k) + \nabla \cdot (\alpha_k \rho_k \mathbf{u}_k) = 0, \\
& \partial_t (\alpha_k \rho_k \mathbf{u}_k) + \nabla \cdot (\alpha_k \rho_k \mathbf{u}_k \otimes \mathbf{u}_k) + \nabla p_k = S_k^m, \\
& \partial_t (\alpha_k \rho_k E_k) + \nabla \cdot (\alpha_k (\rho_k E_k + p_k) \mathbf{u}_k) = S_k^e,
\end{aligned}
\end{equation}
$$

where $\alpha_k$ is the volume fraction of phase $k \in \{1, 2\}$, $\rho_k$ is its density, $\mathbf{u}_k$ its velocity, $p_k$ its pressure and $E_k$ its energy. The source terms in the momentum and energy equations include the contribution of external forces, drag forces and heat exchanges. These models are closed by the definitions of total energies,

$$
(1.2) \qquad E_k = \frac{1}{2} |\mathbf{u}_k|^2 + e_k,
$$

where $e_k$ is the specific internal energy of phase $k$, the equations of state

$$
\begin{equation} \tag{1.3}
\begin{aligned}
p_k &= p_k(\rho_k, \rho_k e_k), \\
\theta_k &= \theta_k(\rho_k, \rho_k e_k)
\end{aligned}
\end{equation}
$$

where $\theta_k$ is the temperature of phase $k$ (which may appear in $S_k^e$), and the equality of pressures

$$
(1.4) \qquad p_1(\rho_1, \rho_1 e_1) = p_2(\rho_2, \rho_2 e_2)
$$

which permits to deduce the void fractions since $\alpha_1 + \alpha_2 = 1$.

It is well-known that this model is not well-posed in the sense of Hadamard,
since the convective part entails complex eigenvalues [16, 19]. Several directions have
been explored to stabilize this model, like the addition of viscous phenomena, surface
tension [16], virtual mass [2, 15] or internal capillary forces [4]. It is also possible to
consider two-velocity and two-pressure models [1, 17], which are well-posed for smooth
solutions [8].

<figure><img src="image_3.png" /><figcaption>FIGURE 1. Flow configuration</figcaption></figure>

Our work carries a derivation of model (1.1)-(1.4). We start from a thin three-
dimensional domain: the vertical length is assumed to be much smaller than the hor-
izontal lengths, the ratio being described by the thickness parameter $\epsilon = D/L \ll 1$
(see Figure 1). The standard averaging process is replaced by an integration with
respect to the vertical variable, resulting in a two-dimensional averaged model. In order to obtain a two-velocity model, our derivation circumvents the smoothness of the velocity field by making two complementary choices on the three-dimensional model. First, instead of imposing the continuity of the velocity field through the interfaces between phases, we allow a tangential friction at the interface by prescribing a Navier condition on the relative tangential velocity. Second, in order to concentrate the efforts on the interaction of this interface condition with the averaging process, we study a simplified configuration: we assume an initial and persistent stratified configuration, in the sense that one fluid remains below the other, and the position of the interface can be parametrized by a function which only depends on the horizontal space variables $(x, y)$. We finally set walls at the bottom and at the top. After rescaling the equations, we perform an asymptotic expansion of the microscopic model with respect to $\epsilon$ and average the equations to obtain the two-dimensional model (1.1)-(1.4) at the principal order. Thanks to this particular setting, it is possible to understand the influence of the different parameters of the initial model on the averaged equations, with respect to the thickness parameter. Compared to the works of Stewart and Wendroff and Ransom and Hicks [20, 17], the introduction of Navier conditions and the asymptotic analysis enable the deduction of a two-velocity averaged model without additional heuristic arguments. Moreover, our framework enables to deduce interaction source terms. Since this framework is intrinsically multidimensional, the extension of mathematical homogenization described above seems out of reach for the moment.

This paper is divided into three parts. Section 2 introduces the configuration of the flow and the three-dimensional microscopic model in the barotropic case. Section 3 presents the averaged two-velocity, one-pressure model derived from the equations of the first section and develops the proof of this result. Section 4 is dedicated to the extension to the non-barotropic case, therefore adding two energy equations in the three-dimensional model. The resulting averaged model is described by two velocities, one pressure and two temperatures.

**2. Barotropic two-phase stratified flows.** In this section, we detail the configuration of the flow, with a particular care for the description of the stratification. We then introduce the compressible barotropic Navier-Stokes systems which describe the flow. These systems are completed by Navier boundary conditions, which model the friction of the fluids on the top and bottom of the domain. We also prescribe Navier conditions on the relative velocity of the fluids at the interface, which depict the friction between both fluids tangentially to the interface. Finally, we develop the different scalings which take place in our configuration and deduce a dimensionless version of the three-dimensional viscous model.

**2.1. Stratified configuration and definition of the interface.** We consider a stratified two-phase flow in a thin domain

$$ (2.1) \qquad \Omega_D = \{(\mathbf{x}, z) = (x, y, z) \in \mathbb{T}_L^2 \times (0, D)\}, \quad L, D > 0 $$

where $\mathbb{T}_L^2$ denotes the two-dimensional torus of period $L$.

The model we study represents a compressible two-phase flow where the fluids are supposed immiscible. The domain filled by each fluid is defined thanks to the top and bottom boundary and the interface between each fluid. We assume that the interface $\Gamma_D$ can be described by the use of continuous function $\alpha_1 : \mathbb{R}_+ \times \mathbb{R}_L^2 \to (0, 1)$ in the
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Flow configuration</figcaption></figure>

following way:

$$ (2.2) \qquad \forall t \in \mathbb{R}_+, \quad \Gamma_D(t) = \{ (\mathbf{x}, D\alpha_1(t, \mathbf{x})), \mathbf{x} \in \mathbb{T}_L^2 \}. $$

The existence of $\alpha_1$ ensures that the stratified configuration is preserved through time. Let us define a velocity field $\mathbf{v}_i : \mathbb{R}_+ \times \mathbb{R}^2 \to \mathbb{R}^3$ which describes the velocity of particles located at the interface. The displacement of the latter is linked to the component of $\mathbf{v}_i$ which is normal to the interface thanks to the kinematic equation

$$ (2.3) \qquad \partial_t(D\alpha_1) - \mathbf{v}_i \cdot \left( \sqrt{\|\nabla_h(D\alpha_1)\|^2 + 1} \right) \mathbf{n}_{\Gamma_D} = 0 \quad \text{in } \mathbb{R}_+ \times \mathbb{T}_L^2, $$

$$ \text{where} \qquad \left(\sqrt{\|\nabla_h(D\alpha_1)\|^2 + 1}\right) \mathbf{n}_{\Gamma_D} = \left((-\nabla_h(D\alpha_1) \quad 1)\right)^\mathrm{T} $$

is the unitary normal vector to the interface pointing upwards (see Figure 2) and $\nabla_h$ is the horizontal gradient operator $(\partial_x \quad \partial_y)^T$.

The domain $\Omega_{D,1}(t)$ filled by the first phase is described by the characteristic function $\chi_{D,2} = 1 - \chi_{D,1}$.

$$ \forall \mathbf{x} \in \mathbb{T}_L^2, \quad \chi_{D,1}(t, \mathbf{x}, z) = \begin{cases} 1 & \text{for } 0 < z < D\alpha_1(t, \mathbf{x}), \\ 0 & \text{otherwise.} \end{cases} $$

Similarly, the domain $\Omega_{D,2}(t)$ filled by the second phase is described by the characteristic function $\chi_{D,2} = 1 - \chi_{D,1}$.

**2.2. The compressible barotropic stratified model.** Let $k \in \{1, 2\}$. Fluid $k$ is described by its density $\rho_k : \mathbb{R}_+^* \times \Omega_{D,k} \to \mathbb{R}_+$, its velocity $\mathbf{v}_k = (u_k, v_k, w_k) : \mathbb{R}_+^* \times \Omega_{D,k} \to \mathbb{R}^3$ and its pressure $p_k : \mathbb{R}_+^* \times \Omega_{D,k} \to \mathbb{R}$. Since we consider viscous flows, we denote $\mu_k$ and $\lambda_k$ the shear and bulk viscosity of the fluid. This flow is described by two compressible, barotropic and viscous Navier-Stokes systems — one for each fluid. Let $h_0 = 0$, $h_1 = \alpha_1$ and $h_2 = 1$. We consider the following system for

$$ k \in \{1, 2\}, \quad t > 0, \quad (\mathbf{x}, z) \in \Omega_{D,k}(t) = \mathbb{T}_L^2 \times (Dh_{k-1}(t, \mathbf{x}), Dh_k(t, \mathbf{x})) : $$
$$ (2.4) \qquad \partial_t \rho_k + \nabla \cdot (\rho_k \mathbf{v}_k) = 0, $$

$$ (2.5) \qquad \partial_t (\rho_k \mathbf{v}_k) + \nabla \cdot (\rho_k \mathbf{v}_k \otimes \mathbf{v}_k) + \nabla p_k = \nabla \cdot \mathcal{S}_k, $$

$$ (2.6) \qquad p_k = p_k(\rho_k) $$

where the viscous tensor $\mathcal{S}_k : \mathbb{R}_+^* \times \Omega_{D,k} \to \mathcal{M}_3(\mathbb{R})$ associated with fluid $k$ is of the form

$$ (2.7) \qquad \mathcal{S}_k = 2\mu_k \frac{(\nabla \otimes \mathbf{v}_k + (\nabla \otimes \mathbf{v}_k)^T)}{2} + \lambda_k (\nabla \cdot \mathbf{v}_k)\mathbb{I} $$

and the equations of state (2.6) are regular and increasing. Equation (2.4) illustrates the conservation of mass. The next one is the momentum equation, which is obtained thanks to Newton's second principle.

**Remark 2.1.** In the equations above, we did not include any effects of gravity or surface tension. Their effects may help to prove the preservation of the stratified configuration and generally influence the stability of the configuration [16]. Moreover, we assume in all this work that the systems we study admit strong solutions.

To complete the system description, boundary and interface conditions are necessary. Denoting $\mathbf{v}_{k,h}$ the horizontal component of $\mathbf{v}_k$ for $k \in \{1, 2\}$, and

$$ \mathcal{S}_k = \begin{pmatrix} \mathcal{S}_{k,h} \\ \mathcal{S}_{k,z} \end{pmatrix} = \left( \begin{array}{c|c} \mathcal{S}_{k,h} & \mathcal{S}_{k,hz} \\ \hline \mathcal{S}_{k,zh} & \mathcal{S}_{k,zz} \end{array} \right), \quad \text{where } \mathcal{S}_{k,hh} \in \mathcal{M}_{2,2}(\mathbb{R}), $$

we prescribe Navier boundary conditions on the top and bottom of the domain:

$$ (2.8) \qquad \begin{cases} w_1|_{z=0} = 0, \\ -\mathcal{S}_{1,hz} + \kappa_1 \mathbf{v}_{1,h}|_{z=0} = 0, \end{cases} $$

$$ (2.9) \qquad \begin{cases} w_2|_{z=D} = 0, \\ \mathcal{S}_{2,hz} + \kappa_2 \mathbf{v}_{2,h}|_{z=D} = 0. \end{cases} $$

The first equation of each boundary condition depicts the impermeability of the domain. The second one describes the horizontal friction of the fluid on the boundary. The quantities $\kappa_1$ and $\kappa_2$ are constant friction coefficients at the bottom and top boundary respectively.

There remains to prescribe interface conditions. Defining the total stress tensor

$$ (2.10) \qquad \sigma_k = (-p_k + \lambda_k(\nabla \cdot \mathbf{v}_k))\mathbb{I} + 2\mu_k \mathcal{D}_k = -p_k\mathbb{I} + \mathcal{S}_k, $$

we prescribe the continuity of its normal component to the interface:

$$ (2.11) \qquad \sigma_1 \mathbf{n}_{\Gamma_D}|_{z=D\alpha_1} = \sigma_2 \mathbf{n}_{\Gamma_D}|_{z=D\alpha_1}. $$

In addition to the previous condition, we prescribe Navier conditions at the interface. They read

$$ (2.12) \qquad \begin{cases} \forall k \in \{1, 2\}, & (\mathbf{v}_k - \mathbf{v}_i) \cdot \mathbf{n}_{\Gamma_D}|_{z=D\alpha_1} = 0 \\ (\mathcal{S}_1 \mathbf{n}_{\Gamma_D} + \kappa_i(\mathbf{v}_1 - \mathbf{v}_2))_{\tan}|_{z=D\alpha_1} = 0, \end{cases} $$
where $X_{\tan|_{z=D\alpha_1}}$ denotes the component of the vector $X$ which is tangential to the interface and $\kappa_i$ is the friction coefficient at the interface. The first equation expresses the continuity of the velocity in the normal direction to the interface. The second one describes the friction of the fluids tangentially to the interface.

The orthogonal component of the Navier interface condition (2.12) enables us to rewrite the fraction equation (2.3) as

$$ (2.13) \qquad \partial_t(D\alpha_1) + \mathbf{v}_{1,h}(D\alpha_1) \cdot \nabla_h(D\alpha_1) = w_1(D\alpha_1). $$

*Remark 2.2.* There is no need of any piece of information about the tangential component of $\mathbf{v}_i$, the normal component being the only useful one in this paper (see (2.3)). Looking at the tangential component of equation (2.11), using the fact that $((-p_k + \lambda_k(\nabla \cdot \mathbf{v}_k))\mathbf{n}_{\Gamma_D})_{\tan|_{z=D\alpha_1}} = 0$ for $k \in \{1, 2\}$, we have

$$ (2.14) \qquad (\mathbb{S}_1 \mathbf{n}_{\Gamma_D})_{\tan|_{z=D\alpha_1}} = (\mathbb{S}_2 \mathbf{n}_{\Gamma_D})_{\tan|_{z=D\alpha_1}}. $$

We can then rewrite the second equation of the interface condition (2.12) with $\mathbb{S}_2$ instead of $\mathbb{S}_1$.

*Remark 2.3.* We do not introduce different notations for the scalar products in two and three dimensions. The context should make the meaning of the scalar product clear.

*Remark 2.4.* The choice of the interface condition (2.12) has a lot of implications on the model that is derived thanks to dimension reduction. For example, an alternative possibility is the prescription of a usual Dirichlet condition with respect to the relative velocity, i.e the continuity of the velocity field at the interface:

$$ \mathbf{v}_1(D\alpha_1) = \mathbf{v}_2(D\alpha_1). $$

*Remark 2.5.* The domain occupied by each fluid being advected by the flow, we have the equation

$$ (2.15) \qquad \forall k \in \{1, 2\}, \qquad \partial_t \chi_{D,k} + \mathbf{v}_k \cdot \nabla \chi_{D,k} = 0. $$

Thanks to the Navier interfacial condition, the previous equation is compatible with the fraction equation (2.3). Indeed, integrating equation (2.15) for the first fluid over $\Omega_D$ yields

$$ \int_{\Omega_D} (\partial_t \chi_{D,1} + \mathbf{v}_1 \cdot \nabla \chi_{D,1}) = 0. $$

$$ \text{Firstly,} \qquad \int_{\Omega_D} \partial_t \chi_{D,1} = \int_{\mathbb{T}_L^2} \partial_t \int_{(0,D)} \chi_{D,1} = \int_{\mathbb{T}_L^2} \partial_t(D\alpha_1) $$

$$ \text{Secondly,} \qquad \int_{\Omega_D} \mathbf{v}_1 \cdot \nabla \chi_{D,1} = \int_{\Omega_D} \div(\mathbf{v}_1 \chi_{D,1}) - \int_{\Omega_{D,1}} \div(\mathbf{v}_1). $$

Thanks to the divergence theorem and the impermeability condition (2.8),

$$ \int_{\Omega_D} \mathbf{v}_1 \cdot \nabla \chi_{D,1} = - \int_{\Gamma_D} \mathbf{v}_1 \cdot \mathbf{n}_{\Gamma_D} dS. $$

Since the surface $\Gamma_D$ is parametrized by the function $\varphi(x, y) = (x, y, D\alpha_1(x, y))$, a change of variables gives

$$ \int_{\Gamma_D} \mathbf{v}_1 \cdot \mathbf{n}_{\Gamma_D} dS = \int_{\mathbb{T}_L^2} (\mathbf{v}_{1,h}(D\alpha_1) \cdot \nabla_h(D\alpha_1) - w_1(D\alpha_1))(x, y)dxdy = 0. $$
The normal component of the Navier interface condition (2.12) then gives the compatibility between equations (2.15) and (2.13). The reasoning is the same for the second fluid.

**2.3. Rescaling the equations in the thin domain assumption.** Let $\varepsilon = D/L$ be a small parameter corresponding to the thin domain assumption. We then introduce dimensionless variables $\tilde{x} = x/L$, $\tilde{y} = y/L$ and $\tilde{z} = z/D = z/(\varepsilon L)$. Denoting $T$ the characteristic time of the flow, we also introduce $\tilde{t}/T$. Denoting $\mathbb{T}^2 = \mathbb{T}_1^2$, switching the variables allows to study the flow on the rescaled domain

$$ (2.16) \qquad \Omega = \mathbb{T}^2 \times (0, 1). $$

In the same spirit, we introduce characteristic dimensions for every unknown in the system: $U = L/T$ for the horizontal velocity, $W = \varepsilon U$ for the vertical velocity, $R$ for the density, $P = RU^2$ for the pressure, $M = RUL$ for the viscosity, and $K = RU$ for the friction coefficients. Denote the dimensionless unknowns $\tilde{u} = u/U$, $\tilde{w} = w/W$, $\tilde{\rho} = \rho/R$, $\tilde{p} = p/P$. Denote $\tilde{\mu} = \mu/M$ the inverse of the Reynolds number, $\tilde{\lambda} = \lambda/M$ and $\tilde{\kappa} = \kappa/K$. We state that the dimensionless viscosity and friction coefficients are constant parameters which may depend on $\varepsilon$. Moreover, we assume that their scaling is independent of the fluid they characterize.

Denote $\chi_k$ the rescaled version of $\chi_{D,k}$ for $k \in \{1, 2\}$, $\Omega_k$ the rescaled version of $\Omega_{D,k}$ and $\Gamma$ the rescaled version of $\Gamma_D$. Removing the $\tilde{\cdot}$ notation for the sake of clarity, and denoting $\Delta_h = \nabla_h \cdot \nabla_h$, we rewrite the equations (2.4)-(2.15) as follows:

**LEMMA 2.6.** *The rescaling of model (2.4)-(2.15) results in the following set of equations on the rescaled domain $\Omega$.*

• Dimensionless advection of the domain

$$ (2.17) \qquad \partial_t \chi_k + \mathbf{v}_k \cdot \nabla \chi_k = 0. $$

• Dimensionless fraction equation

$$ (2.18) \qquad \partial_t \alpha_1 + \mathbf{v}_{1,h}(\alpha_1) \cdot \nabla_h \alpha_1 = w_1(\alpha_1). $$

• Dimensionless mass conservation equations

$$ (2.19) \qquad \forall k \in \{1, 2\}, \qquad \partial_t \rho_k + \nabla_h \cdot (\rho_k \mathbf{v}_{k,h}) + \partial_z (\rho_k w_k) = 0, $$

• Dimensionless horizontal momentum equations

$$ (2.20) \qquad \forall k \in \{1, 2\}, \quad \partial_t (\rho_k \mathbf{v}_{k,h}) + \nabla_h \cdot (\rho_k \mathbf{v}_{k,h} \otimes \mathbf{v}_{k,h}) + \partial_z (\rho_k w_k \mathbf{v}_{k,h}) + \nabla_h p_k \\ = (\lambda_k + \mu_k) \nabla_h (\nabla \cdot \mathbf{v}_k) + \mu_k \Delta_h \mathbf{v}_{k,h} + \frac{\mu_k}{\varepsilon^2} \partial_{zz} \mathbf{v}_{k,h}, $$

• Dimensionless vertical momentum equations

$$ (2.21) \qquad \varepsilon(\partial_t(\rho_k w_k) + \nabla_h \cdot (\rho_k \mathbf{v}_{k,h} w_k) + \partial_z(\rho_k (w_k)^2)) + \frac{1}{\varepsilon} \partial_z p_k \\ = \frac{1}{\varepsilon}((\lambda_k + \mu_k)\partial_z(\nabla \cdot \mathbf{v}_k) + \mu_k \partial_{zz} w_k) + \mu_k \varepsilon \Delta_h w_k. $$

• Dimensionless Navier boundary conditions

$$ (2.22) \qquad \begin{cases} w_1|_{z=0} = 0, \\ -\frac{\mu_1}{\varepsilon} \partial_z \mathbf{v}_{1,h} + \kappa_1 \mathbf{v}_{1,h}|_{z=0} = 0, \end{cases} $$
$$ (2.23) \qquad \left\{ \begin{array}{l} w_2|_{z=1} = 0, \\ \frac{\mu_2}{\varepsilon} \partial_z \mathbf{v}_{2,h} + \kappa_2 \mathbf{v}_{2,h}|_{z=1} = 0. \end{array} \right. $$

• Dimensionless continuity of the normal stress at the interface

$$ (2.24) \qquad \left( \partial_x \alpha_1 (p_1 - \lambda_1 \nabla \cdot \mathbf{v}_1 - 2\mu_1 \partial_x u_1) - \mu_1 \partial_y \alpha_1 (\partial_y u_1 + \partial_x v_1) + \mu_1 \partial_x w_1 + \frac{\mu_1}{\varepsilon^2} \partial_z u_1 \right) \Big|_{z=\alpha_1} = (\cdots)_2 \Big|_{z=\alpha_1}, $$

$$ (2.25) \qquad \left( \partial_y \alpha_1 (p_1 - \lambda_1 \nabla \cdot \mathbf{v}_1 - 2\mu_1 \partial_y v_1) - \mu_1 \partial_x \alpha_1 (\partial_y u_1 + \partial_x v_1) + \mu_1 \partial_y w_1 + \frac{\mu_1}{\varepsilon^2} \partial_z v_1 \right) \Big|_{z=\alpha_1} = (\cdots)_2 \Big|_{z=\alpha_1}, $$

$$ (2.26) \qquad \left( -p_1 + \lambda_1 \nabla \cdot \mathbf{v}_1 + 2\mu_1 \partial_z w_1 - \mu_1 \partial_x \alpha_1 (\partial_z u_1 + \varepsilon^2 \partial_x w_1) - \mu_1 \partial_y \alpha_1 (\partial_z v_1 + \varepsilon^2 \partial_y w_1) \right) \Big|_{z=\alpha_1} = (\cdots)_2 \Big|_{z=\alpha_1}, $$

where $(\cdots)_2$ denotes the same expression as the left-hand side of every condition with the subscript 1 replaced by 2, except for $\alpha_1$.

• Dimensionless normal Navier interface condition

$$ (2.27) \qquad \left( -u_1 \partial_x \alpha_1 - v_1 \partial_y \alpha_1 + w_1 \right) \Big|_{z=\alpha_1} = \left( -u_2 \partial_x \alpha_1 - v_2 \partial_y \alpha_1 + w_2 \right) \Big|_{z=\alpha_1} $$

• Dimensionless tangential Navier interface condition

$$ (2.28) \qquad \mu_1 \varepsilon \left( \frac{1}{\varepsilon^2} \partial_z u_1 + \partial_x w_1 - 2 \partial_x \alpha_1 (\partial_x u_1 - \partial_z w_1) - \partial_y \alpha_1 (\partial_y u_1 + \partial_x v_1) \right. \\ \left. - \partial_x \alpha_1 \partial_y \alpha_1 (\partial_z v_1 + \varepsilon^2 \partial_y w_1) - (\partial_x \alpha_1)^2 (\partial_z u_1 + \varepsilon^2 \partial_x w_1) \right) \\ + \kappa_i (u_1 - u_2 + \varepsilon^2 \partial_x \alpha_1 (w_1 - w_2)) \Big|_{z=\alpha_1} = 0, $$

$$ (2.29) \qquad \mu_1 \varepsilon \left( \frac{1}{\varepsilon^2} \partial_z v_1 + \partial_y w_1 - 2 \partial_y \alpha_1 (\partial_y v_1 - \partial_z w_1) - \partial_x \alpha_1 (\partial_x v_1 + \partial_y u_1) \right. \\ \left. - \partial_y \alpha_1 \partial_x \alpha_1 (\partial_z u_1 + \varepsilon^2 \partial_x w_1) - (\partial_y \alpha_1)^2 (\partial_z v_1 + \varepsilon^2 \partial_y w_1) \right) \\ + \kappa_i (v_1 - v_2 + \varepsilon^2 \partial_y \alpha_1 (w_1 - w_2)) \Big|_{z=\alpha_1} = 0. $$

*Proof.* The introduction of the dimensionless variables and unknowns in the various equations, boundary conditions and interface conditions is performed using the chain rule. The characteristic function $\chi_k$ is advected by the dimensionless velocity, hence equation (2.17). The rescaling procedure immediately yields the fraction equation (2.18), the mass conservation equation (2.19) and the momentum equation, which we split into two components — the horizontal component (2.20) and the vertical component (2.21).

Rescaling the boundary conditions yields equations (2.22) and (2.23). Since the normal vector to the interface is now

$$ (2.30) \qquad \mathbf{n}_{\Gamma} = \frac{1}{\sqrt{1 + \varepsilon^2 |\nabla_h \alpha_1|^2}} \begin{pmatrix} -\varepsilon \nabla_h \alpha_1 \\ 1 \end{pmatrix}, $$
the rescaling of equation (2.11) detailed into every component gives (2.24), (2.25) and (2.26).

Now, the Navier interface conditions must be rescaled. Using (2.30), we can rescale the continuity of the velocity orthogonally to the interface and thus obtain (2.27). The tangential condition requires more computations. First, let us introduce two tangential vectors to the interface:

$$ \mathbf{t}_{\Gamma_D}^x = \begin{pmatrix} 1 \\ 0 \\ \partial_x(D\alpha_1) \end{pmatrix} \quad \text{and} \quad \mathbf{t}_{\Gamma_D}^y = \begin{pmatrix} 0 \\ 1 \\ \partial_y(D\alpha_1) \end{pmatrix} $$

Those two are rescaled as

$$ \mathbf{t}_{\Gamma}^x = \begin{pmatrix} 1 \\ 0 \\ \varepsilon \partial_x \alpha_1 \end{pmatrix} \quad \text{and} \quad \mathbf{t}_{\Gamma}^y = \begin{pmatrix} 0 \\ 1 \\ \varepsilon \partial_y \alpha_1 \end{pmatrix}. $$

Then, the tangential component of condition (2.12) is equivalent to the system

$$ \begin{aligned} (S_1 \mathbf{n}_{\Gamma_D} + \kappa_i (\mathbf{v}_1 - \mathbf{v}_2)) \cdot \mathbf{t}_{\Gamma_D}^x (D\alpha_1) &= 0, \\ (S_1 \mathbf{n}_{\Gamma_D} + \kappa_i (\mathbf{v}_1 - \mathbf{v}_2)) \cdot \mathbf{t}_{\Gamma_D}^y (D\alpha_1) &= 0. \end{aligned} $$

which is rescaled as (2.28)-(2.29).

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

**3. The asymptotic system in the barotropic case.** In this section, we present the asymptotic model derived from the compressible barotropic double Navier-Stokes system (2.4)-(2.13). The first subsection introduces the fundamental notion of average that will be used later, and states the theorem presenting the asymptotic model. In a second part, we develop the derivation process.

**3.1. Presentation of the averaged model.** Before presenting our asymptotic model, we need to introduce some notations. As explained earlier, integrating the Navier-Stokes equations causes the appearance of averaged unknowns. Since the majority of the terms in the Navier-Stokes equations are nonlinear in the compressible setting and since the multiplication and the integration do not commute, it is necessary to define different notions of averages.

**DEFINITION 3.1 (Averages).** Let $f$ be a function defined on the domain filled by fluid $k \in \{1, 2\}$. The average of $f$, which we denote $\overline{f}$, is defined by

$$ (3.1) \qquad \forall (t, \mathbf{x}) \in \mathbb{R}_+ \times \mathbb{T}^2, \qquad \overline{f}(t, \mathbf{x}) = \frac{1}{\alpha_k} \int_{h_{k-1}}^{h_k} f(t, \mathbf{x}, z) \mathrm{d}z. $$

The Favre average of $f$ with respect to $\rho_k$, which we denote $\langle f \rangle$, is defined for all $(t, \mathbf{x}) \in \mathbb{R}_+ \times \mathbb{T}^2$ by

$$ (3.2) \qquad \langle f \rangle(t, \mathbf{x}) = \frac{\overline{\rho_k f}}{\overline{\rho_k}} = \frac{1}{\alpha_k \overline{\rho_k}} \int_{h_{k-1}}^{h_k} \rho_k f(t, \mathbf{x}, z) \mathrm{d}z. $$

**THEOREM 3.2.** Assume the following scalings for viscosity and friction coefficients:

$$ (3.3) \qquad \forall k \in \{1, 2\}, \qquad \begin{cases} \mu_k = \widehat{\mu_k} \varepsilon^\tau, & \lambda_k = \widehat{\lambda_k} \varepsilon^\tau, & 0 < \tau < 2; \\ \kappa_i = \widehat{\kappa_i} \varepsilon, & \kappa_k = \widehat{\kappa_k} \varepsilon^\xi, & \xi \ge 1, \end{cases} $$
where the coefficients $\hat{\mu}_k$, $\hat{\lambda}_k$, $\hat{\kappa}_i$ and $\hat{\kappa}_k$ are positive constants independent of $\varepsilon$.

Define

$$ \delta_{\xi} = \begin{cases} 1 & \text{if } \xi = 1, \\ 0 & \text{if } \xi > 1. \end{cases} $$

Let $(\alpha_1, \rho_1, \mathbf{v}_1, p_1, \alpha_2, \rho_2, \mathbf{v}_2, p_2)$ be a strong solution to the Navier-Stokes system with Navier boundary and interface conditions (2.4)-(2.13).

The functions $(\alpha_1, \overline{p_1}, \langle \mathbf{v}_{1,h} \rangle, \overline{p_1}, \alpha_2, \overline{p_2}, \langle \mathbf{v}_{2,h} \rangle, \overline{p_2})$ satisfy the following closed system:

• **Averaged mass conservation equations**

$$ (3.4) \qquad \forall k \in \{1, 2\}, \qquad \partial_t(\alpha_k \overline{p_k}) + \nabla_h \cdot (\alpha_k \overline{p_k} \langle \mathbf{v}_{k,h} \rangle) = 0. $$

• **Averaged momentum equations**

$$ (3.5) \quad \begin{aligned} \forall k \in \{1, 2\}, \quad & \partial_t(\alpha_k \overline{p_k} \langle \mathbf{v}_{k,h} \rangle) + \nabla_h \cdot (\alpha_k \overline{p_k} \langle \mathbf{v}_{k,h} \rangle \otimes \langle \mathbf{v}_{k,h} \rangle) + \nabla_h (\alpha_k \overline{p_k}) \\ & = p_i \nabla_h \alpha_k - \delta_{\xi} \hat{\kappa}_k \langle \mathbf{v}_{k,h} \rangle + (-1)^k \hat{\kappa}_i (\langle \mathbf{v}_{1,h} \rangle - \langle \mathbf{v}_{2,h} \rangle), \end{aligned} $$

where $p_i$ is any convex combination of the averaged pressures $\overline{p_k}$. This equation is satisfied up to an error of order $\varepsilon^\tau + \varepsilon^{2-\tau} + (1 - \delta_\xi)\varepsilon^{\xi-1}$.

• **Averaged equations of state**

$$ (3.6) \qquad \forall k \in \{1, 2\}, \qquad \overline{p_k} = p_k(\overline{p_k}). $$

This equation is satisfied up to an error of order $\varepsilon^\tau$.

• **Equality of the pressures**

$$ (3.7) \qquad \overline{p_1} = \overline{p_2}. $$

This equation is satisfied up to an error of order $\varepsilon^\tau$.

• **Saturation of the domain**

$$ (3.8) \qquad \alpha_1 + \alpha_2 = 1. $$

*Remark 3.3.* The approximation errors mentioned above show the necessity of prescribing these specific scalings of viscosity and friction coefficients. Our work provides the approximation errors and the necessary conditions for the model to be valid. The exponents $\xi$ and $\tau$ are independent.

*Remark 3.4.* In the tangential component of the Navier interface condition (2.12), we can replace the term $\kappa_i(\mathbf{v}_1 - \mathbf{v}_2)$ by a quadratic one $\kappa_i|\mathbf{v}_1 - \mathbf{v}_2|(\mathbf{v}_1 - \mathbf{v}_2)$. The resulting term in the averaged model becomes $\kappa_i|\langle \mathbf{v}_1 \rangle - \langle \mathbf{v}_2 \rangle|(\langle \mathbf{v}_1 \rangle - \langle \mathbf{v}_2 \rangle)$ and causes no change to the error approximation.

*Remark 3.5.* The scalings detailed in (3.3) result in an asymptotic model that is compressible and inviscid. In other words, our asymptotic limit corresponds to moderate Mach numbers and high Reynolds numbers.

**3.2. Proof of Theorem 3.2: derivation of the asymptotic model.** The proof of Theorem 3.2 is divided into several steps. Starting from the rescaled equations (2.17)-(2.29), we first simplify the vertical momentum equation and the interface conditions by neglecting terms of high order in $\varepsilon$. Then, we integrate the mass conservation equation and the horizontal momentum equation in the vertical direction. This procedure causes boundary terms to appear in the equations, which must be closed using the rescaled boundary and interface conditions. Finally, we prescribe specific profiles for the viscosity and friction coefficients for the asymptotic model to be valid.
**3.2.1. Simplification of the vertical momentum equation and interface conditions.** Some of the equations and conditions listed in Lemma 2.6 are very complex once rescaled. However, some terms are of high order in $\epsilon$, and therefore become very small in our setting ($\epsilon \ll 1$). Dropping these terms allows a first simplification of the system. To make sure that the derived system is really an approximation of the initial double Navier-Stokes system, we keep trace of the errors that we cause along the way using Landau notations. Since we have assumed that the scaling of the viscosity and friction coefficients is independent of the fluid, we drop the subscript $k$ in the Landau notations.

Performing the asymptotic analysis impacts the vertical component of the momentum equation (2.21) as well as interface conditions (2.26), (2.28) and (2.29) in the following way.

• Asymptotic vertical momentum equation

$$ (3.9) \quad \forall k \in \{1, 2\}, \quad \partial_z p_k = (\lambda_k + \mu_k) \partial_z (\nabla \cdot \mathbf{v}_k) + \mu_k \partial_{zz} w_k + O(\epsilon^2). $$

• Asymptotic continuity of the vertical component of the normal stress at the interface

$$ (3.10) \quad \left. \left( -p_1 + \lambda_1 \nabla \cdot \mathbf{v}_1 + 2\mu_1 \partial_z w_1 - \mu_1 \nabla_h \alpha_1 \cdot \partial_z \mathbf{v}_{1,h} \right) \right|_{z=\alpha_1} \\ = (\cdots)_2 \Big|_{z=\alpha_1} + O(\mu\epsilon^2). $$

• Asymptotic tangential Navier interface condition

$$ (3.11) \quad \left. \left( \frac{\mu_1}{\epsilon} \partial_z \mathbf{v}_{1,h} + \kappa_i (\mathbf{v}_{1,h} - \mathbf{v}_{2,h}) \right) \right|_{z=\alpha_1} = O(\mu\epsilon + \kappa_i\epsilon^2) $$

*Remark 3.6.* In the previous paragraph, we did not simplify the horizontal momentum equation (2.20) since the terms of order $\epsilon^2$ are not negligible. Indeed, multiplying (2.20) by $\epsilon^2$ and neglecting the quadratic terms in $\epsilon$ indeed yields

$$ (3.12) \qquad \mu_k \partial_{zz} \mathbf{v}_{k,h} = O(\epsilon^2). $$

The term $\mu_k \partial_{zz} \mathbf{v}_{k,h}/\epsilon^2$ is thus of order 1 and cannot be neglected in the equation.

**3.2.2. Averaging the equations.** Now, we perform the averaging of the mass conservation equation and momentum equations for each fluid in the vertical direction, therefore obtaining a two-dimensional system. We approximate the boundary terms which appear during the integration by terms which only depend on the averaged unknowns in order to close the system. We keep trace of the approximation errors during the whole averaging process. We perform the averaging for the first fluid, the second fluid being treated in the same way.

First, average the mass conservation equations (2.19). We multiply the mass conservation equation (2.19) by $\chi_1$ and we integrate it with respect to $z \in (0, 1)$ — which is the same as integrating the equation with respect to $z \in (0, \alpha_1)$:

$$ \int_0^1 \chi_1 (\partial_t \rho_1 + \nabla_h \cdot (\rho_1 \mathbf{v}_{1,h}) + \partial_z (\rho_1 w_1)) (z) dz = 0. $$

Since $\chi_1 \partial_t \rho_1 = \partial_t (\chi_1 \rho_1) - \rho_1 \partial_t \chi_1$ with $\partial_t \chi_1 = \partial_t \alpha_1 \delta_{z=\alpha_1}$,

$$ \int_0^1 \chi_1 \partial_t \rho_1(z) dz = \partial_t \int_0^{\alpha_1} \rho_1(z) dz - \rho_1(\alpha_1) \partial_t \alpha_1 = \partial_t (\alpha_1 \overline{\rho_1}) - \rho_1(\alpha_1) \partial_t \alpha_1. $$
Proceeding similarly for the second term of the equation, we get that

$$ \partial_t(\alpha_1 \overline{\rho_1}) - \rho_1(\alpha_1)\partial_t\alpha_1 + \nabla_h \cdot (\alpha_1 \overline{\rho_1 \mathbf{v}_{1,h}}) - \rho_1 \mathbf{v}_{1,h}(\alpha_1) \cdot \nabla_h(\alpha_1) \\ + \rho_1 w_1(\alpha_1) - \rho_1 w_1(0) = 0. $$

Using the fraction equation (2.18), the bottom Navier condition (2.22) and the relation $\overline{\rho_1 \mathbf{v}_{1,h}} = \overline{\rho_1} \langle \mathbf{v}_{1,h} \rangle$, we finally get the averaged mass conservation equation (3.4) for the first fluid.

The reasoning is the same for the second fluid, the equation being multiplied this time by $\chi_2$ and using the Navier condition (2.23). Note that the averaged mass conservation equations are exact: we did not make any approximations in order to get them.

The averaging of the momentum equations is more complex and is detailed in the following proposition.

**PROPOSITION 3.7.** *The closed horizontal momentum equations read*

$$ (3.13) \quad \forall k \in \{1, 2\}, \quad \partial_t(\alpha_k \overline{\rho_k} \langle \mathbf{v}_{k,h} \rangle) + \nabla_h \cdot (\alpha_k \overline{\rho_k} \langle \mathbf{v}_{k,h} \rangle \otimes \langle \mathbf{v}_{k,h} \rangle) + \nabla_h(\alpha_k \overline{\rho_k}) \\ = p_i \nabla_h \alpha_k - \frac{\kappa_k}{\varepsilon} \langle \mathbf{v}_{k,h} \rangle + (-1)^k \frac{\kappa_i}{\varepsilon} (\langle \mathbf{v}_{1,h} \rangle - \langle \mathbf{v}_{2,h} \rangle) + \mathcal{R}, $$

where

$$ \mathcal{R} = O\left(\mu + \kappa_i \varepsilon + \frac{\kappa_i(\kappa + \varepsilon) + \kappa(\kappa + \varepsilon)}{\mu} + \frac{(\varepsilon\kappa + \varepsilon^2)^2}{\mu^2}\right). $$

*Proof.* Once again, we develop the proof for $k = 1$. Multiply the horizontal momentum equation (2.20) by $\chi_1$ and integrate it with respect to $z \in (0, 1)$. Using the fraction equation (2.18), we obtain

$$ \begin{aligned} \partial_t(\alpha_1 \overline{\rho_1} \langle \mathbf{v}_{1,h} \rangle) + \nabla_h \cdot (\alpha_1 \overline{\rho_1} \langle \mathbf{v}_{1,h} \otimes \mathbf{v}_{1,h} \rangle) + \nabla_h(\alpha_1 \overline{\rho_1}) \\ = p_1(\alpha_1)\nabla_h(\alpha_1) + \frac{\mu_1}{\varepsilon^2}\partial_z \mathbf{v}_{1,h}(\alpha_1) - \frac{\mu_1}{\varepsilon^2}\partial_z \mathbf{v}_{1,h}(0) \\ + \int_0^{\alpha_1} \left( (\lambda_1 + \mu_1)\nabla_h(\nabla \cdot \mathbf{v}_1) + \mu_1 \Delta_h \mathbf{v}_{1,h} \right)(z) dz. \end{aligned} $$

Thanks to the Navier condition (2.22), the equation reads

$$ (3.14) \quad \partial_t(\alpha_1 \overline{\rho_1} \langle \mathbf{v}_{1,h} \rangle) + \nabla_h \cdot (\alpha_1 \overline{\rho_1} \langle \mathbf{v}_{1,h} \otimes \mathbf{v}_{1,h} \rangle) + \nabla_h(\alpha_1 \overline{\rho_1}) = p_1(\alpha_1)\nabla_h \alpha_1 \\ - \frac{\kappa_1}{\varepsilon} \mathbf{v}_{1,h}(0) + \frac{\mu_1}{\varepsilon^2} \partial_z \mathbf{v}_{1,h}(\alpha_1) + \int_0^{\alpha_1} \left( (\lambda_1 + \mu_1)\nabla_h(\nabla \cdot \mathbf{v}_1) + \mu_1 \Delta_h \mathbf{v}_{1,h} \right)(z) dz. $$

The asymptotic Navier condition (3.11) then gives

$$ \frac{\mu_1}{\varepsilon^2} \partial_z \mathbf{v}_{1,h}(\alpha_1) = \frac{\kappa_i}{\varepsilon} (\mathbf{v}_{2,h}(\alpha_1) - \mathbf{v}_{1,h}(\alpha_1)) + O(\mu + \kappa_i \varepsilon). $$

Injecting this new term into (3.14) causes an approximation of order $\mu$. The integral term of the equation being of the same order, we obtain

$$ (3.15) \quad \partial_t(\alpha_1 \overline{\rho_1} \langle \mathbf{v}_{1,h} \rangle) + \nabla_h \cdot (\alpha_1 \overline{\rho_1} \langle \mathbf{v}_{1,h} \otimes \mathbf{v}_{1,h} \rangle) + \nabla_h(\alpha_1 \overline{\rho_1}) = p_1(\alpha_1)\nabla_h \alpha_1 \\ - \frac{\kappa_1}{\varepsilon} \mathbf{v}_{1,h}(0) + \frac{\kappa_i}{\varepsilon} (\mathbf{v}_{2,h}(\alpha_1) - \mathbf{v}_{1,h}(\alpha_1)) + O(\mu + \kappa_i \varepsilon). $$

To close this equation, the pointwise values of the velocity and the pressure must be expressed thanks to their averaged values.
LEMMA 3.8. The pointwise values of the velocity $\mathbf{v}_{1,h}$ can be approximated as follows:

$$ (3.16) \qquad \forall z \in (0, \alpha_1), \quad \mathbf{v}_{1,h}(z) = \langle \mathbf{v}_{1,h} \rangle + O\left(\frac{\varepsilon\kappa + \varepsilon^2}{\mu}\right). $$

The pointwise values of the pressure can be approximated as follows:

$$ (3.17) \qquad p_1(\alpha_1) = p_i + O(\mu), $$

where $p_i$ is defined in Theorem 3.2.

The same approximations remain valid for $\mathbf{v}_{2,h}$ and $p_2$.

*Proof.* Recall that in Remark 3.6, we noticed that

$$ \forall z \in (0, \alpha_1), \quad \mu_1 \partial_{zz} \mathbf{v}_{1,h}(z) = O(\varepsilon^2). $$

Moreover, the bottom boundary Navier condition (2.22) gives us that

$$ \mu_1 \partial_z \mathbf{v}_{1,h}(0) = O(\kappa\varepsilon). $$

Thanks to a Taylor-Lagrange expansion,

$$ (3.18) \qquad \forall z \in (0, \alpha_1), \quad \mu_1 \partial_z \mathbf{v}_{1,h}(z) = O(\kappa\varepsilon + \varepsilon^2). $$

Another Taylor-Lagrange expansion gives

$$ \begin{aligned} \forall (z, z') \in (0, \alpha_1)^2, \quad \mathbf{v}_{1,h}(z) &= \mathbf{v}_{1,h}(z') + (z - z') \partial_z \mathbf{v}_{1,h}(z') + O(\varepsilon^2/\mu), \\ &= \mathbf{v}_{1,h}(z') + O\left(\frac{\varepsilon\kappa + \varepsilon^2}{\mu}\right). \end{aligned} $$

Taking the Favre average of the previous equation with respect to $z'$ yields (3.16). The reasoning is the same for $\mathbf{v}_{2,h}$ with $z \in (\alpha_1, 1)$. We use the top Navier condition (2.23) to get the estimate.

In the case of the pressure, the vertical component of the momentum equations (3.9) provide the estimates

$$ \left\{ \begin{array}{ll} \forall z \in (0, \alpha_1), & \partial_z p_1(z) = O(\mu), \\ \forall z \in (\alpha_1, 1), & \partial_z p_2(z) = O(\mu). \end{array} \right. $$

Therefore, the pointwise values of the pressures can be approximated at order $\mu$ by their average. However, interface condition (3.10) yields

$$ p_1(\alpha_1) = p_2(\alpha_1) + O(\mu). $$

Therefore, $\overline{p_1} = \overline{p_2} + O(\mu)$ and

$$ p_1(\alpha_1) = p_i + O(\mu), $$

where $p_i$ is defined in Theorem 3.2. Hence (3.16) and (3.17).

*Remark 3.9.* The computations above allow more freedom in the definition of $p_i$ than being a convex combination of the averaged pressures. Indeed, the proposition is still valid for any $p_i$ such that $\overline{p_1} - p_i = O(\mu)$.
Approximating these different pointwise values by their (Favre) average yields

$$
\begin{align*}
\partial_t (\alpha_1 \overline{\rho}_1 \langle \mathbf{v}_{1,h} \rangle) + \nabla_h \cdot (\alpha_1 \overline{\rho}_1 \langle \mathbf{v}_{1,h} \otimes \mathbf{v}_{1,h} \rangle) + \nabla_h (\alpha_1 \overline{p}_1) \\
&= p_i \nabla_h \alpha_1 - \frac{\kappa_1}{\varepsilon} \langle \mathbf{v}_{1,h} \rangle + \frac{\kappa_i}{\varepsilon} (\langle \mathbf{v}_{2,h} \rangle - \langle \mathbf{v}_{1,h} \rangle) \\
&\quad + O \left( \mu + \kappa_i \varepsilon + \frac{\kappa_i (\kappa + \varepsilon) + \kappa (\kappa + \varepsilon)}{\mu} \right).
\end{align*}
$$

Remark 3.10. Replacing the pointwise values of the velocities by their Favre average in the Navier interface condition (3.11) and using (3.18), we have

$$
(3.19) \qquad \langle \mathbf{v}_{2,h} \rangle - \langle \mathbf{v}_{1,h} \rangle = O \left( \frac{\kappa + \varepsilon + \mu\varepsilon}{\kappa_i} + \frac{\kappa\varepsilon + \varepsilon^2}{\mu} + \varepsilon^2 \right).
$$

Finally, we have to express the Favre average $\langle \mathbf{v}_{1,h} \otimes \mathbf{v}_{1,h} \rangle$ with respect to the tensor product of the Favre averages. This way, the resulting system will only depend on the averages of the density, the velocity and the pressure.

LEMMA 3.11. The Favre average of the tensor product of the velocity $\mathbf{v}_{1,h}$ with itself is approximated by the tensor product of its Favre averages.

$$
(3.20) \quad \langle \mathbf{v}_{1,h} \otimes \mathbf{v}_{1,h} \rangle = \langle \mathbf{v}_{1,h} \rangle \otimes \langle \mathbf{v}_{1,h} \rangle + O \left( \frac{(\varepsilon\kappa + \varepsilon^2)^2}{\mu^2} \right).
$$

*Proof.* As shown in the previous lemma,

$$
\forall z \in (0, \alpha_1), \quad \mathbf{v}_{1,h}(z) = \langle \mathbf{v}_{1,h} \rangle + f(z),
$$

$$
\text{where} \qquad f = O\left(\frac{\varepsilon\kappa + \varepsilon^2}{\mu}\right) \quad \text{and} \quad \langle f \rangle = 0.
$$

We can thus develop the tensor product $\mathbf{v}_{1,h} \otimes \mathbf{v}_{1,h}$ as

$$
\mathbf{v}_{1,h} \otimes \mathbf{v}_{1,h} = \langle \mathbf{v}_{1,h} \rangle \otimes \langle \mathbf{v}_{1,h} \rangle + 2 \langle \mathbf{v}_{1,h} \rangle \otimes f + f \otimes f.
$$

Since $f \otimes f = O((\kappa\varepsilon + \varepsilon^2)^2/\mu^2)$, taking the Favre average of the equation above yields (3.20). $\square$

Taking all of the above into account, we obtain the averaged momentum equation (3.13). For the second fluid, the reasoning is the same once we have rewritten the tangential Navier condition using Remark 2.2. $\square$

**3.2.3. Closing the system.** To get the set of equations that is presented in the main result, one must make sure that the momentum equations (3.13) remain valid when $\varepsilon$ becomes very small. Three conditions arise:

1. $\mathcal{R} \xrightarrow{\varepsilon \to 0} 0$, which means that the approximation errors vanish when $\varepsilon$ tends to 0;

2. $\langle \mathbf{v}_{1,h} \rangle - \langle \mathbf{v}_{2,h} \rangle$ must be well defined (i.e finite) when $\varepsilon$ tends to 0;

3. the friction terms in the momentum equations must have a limit when $\varepsilon$ tends to 0 that is independent of $\varepsilon$.

These conditions are the reason for the particular scaling requirements for the viscosity and friction coefficients presented in (3.3). Indeed, let us prescribe viscosity and friction scalings as follows:

$$
(3.21) \quad \forall k \in \{1, 2\}, \quad \mu_k = \widehat{\mu}_k \varepsilon^\tau, \quad \kappa_k = \widehat{\kappa}_k \varepsilon^\xi, \quad \kappa_i = \widehat{\kappa}_i \varepsilon^\zeta,
$$
where $\tau, \xi, \zeta \in \mathbb{R}_+$. Rewriting the approximation error $\mathcal{R}$ with the previous scalings yields

$$ \mathcal{R} = O(\varepsilon^{\tau} + \varepsilon^{1+\zeta} + \varepsilon^{\xi+\zeta-\tau} + \varepsilon^{1+\zeta-\tau} + \varepsilon^{2\xi-\tau} + \varepsilon^{1+\xi-\tau} + \varepsilon^{2+2\xi-2\tau} + \varepsilon^{3+\xi-2\tau} + \varepsilon^{4-2\tau}). $$

Satisfying the first condition adds the following constraints on the viscosity and friction scalings:

$$ \zeta + \xi > \tau, \quad 1 + \zeta > \tau, \quad 1 + \xi > \tau, \quad 2\xi > \tau \text{ and } 2 > \tau. $$

Next, in order for the averaged model to be valid, the difference between the averaged velocities must remain finite when $\varepsilon$ becomes very small. Looking back at the relation (3.19) and replacing the viscosities and friction coefficients by their scalings yields

$$ \langle \mathbf{v}_{1,h} \rangle - \langle \mathbf{v}_{2,h} \rangle = O(\varepsilon^{\xi-\zeta} + \varepsilon^{1-\zeta} + \varepsilon^{1+\tau-\zeta} + \varepsilon^{1+\xi-\tau} + \varepsilon^{2-\tau} + \varepsilon^2). $$

The second condition is thus satisfied if

$$ \xi \ge \zeta, \quad 1 \ge \zeta \text{ and } 1 + \tau \ge \zeta. $$

Finally, for the model to be valid, it must be independent of $\varepsilon$. Since $\langle \mathbf{v}_{1,h} \rangle - \langle \mathbf{v}_{2,h} \rangle$ is finite, we must ensure that the friction terms $\frac{\kappa_i}{\varepsilon}(\langle \mathbf{v}_{2,h} \rangle - \langle \mathbf{v}_{1,h} \rangle)$ and $\frac{\kappa_k}{\varepsilon}\langle \mathbf{v}_{k,h} \rangle$ for $k \in \{1, 2\}$ are well defined at the limit $\varepsilon \to 0$. This condition is satisfied if

$$ \xi \ge 1 \text{ and } \zeta \ge 1. $$

Taking into account all the different constraints results in the scalings (3.3):

$$ \forall k \in \{1, 2\}, \quad \mu_k = \hat{\mu}_k \varepsilon^{\tau}, \quad 0 < \tau < 2, \quad \kappa_k = \hat{\kappa}_k \varepsilon^{\xi}, \quad \xi \ge 1 \text{ and } \kappa_i = \hat{\kappa}_i \varepsilon. $$

Dropping the errors yields equations (3.4)-(3.5). Equation (3.8) is given by the definition of the volume fractions. Equation (3.7) has been proven in the demonstration of Lemma 3.8.

At this stage, this system is underdetermined. We indeed have eight unknowns and six equations. It remains to adapt the barotropy laws to the new set of equations to get equations (3.6). Let us prove this result in the case $k = 1$.

Using (3.9), we know that $\partial_z(p_1(\rho_1)) = O(\mu)$. Let $z^* \in (0, \alpha_1)$ such that $\rho_1(z^*) = \overline{\rho_1}$. Thanks to a Taylor-Lagrange expansion,

$$ \forall z \in (0, \alpha_1), \quad p_1(\rho_1(z)) = p_1(\rho_1(z^*)) + O(\mu). $$

Averaging the previous equation with respect to $z$ and dropping the error yields (3.6).

*Remark 3.12.* Compared to the Navier-Stokes system of the first section, the limit of the fraction equation (2.13) is not considered. Indeed, the volume fractions $\alpha_1$ and $\alpha_2$ are deduced from the equality of the averaged pressures (3.7) and the saturation condition (3.8).

**3.3. Recovering the one-velocity one-pressure model.** The choice of prescribing a Navier friction condition at the interface between the fluids is motivated by the desire to derive a two-velocity model. As mentioned in Remark 2.4, other choices for the interface conditions are possible. If we allow the friction coefficient $\kappa_i$ to decrease slower than $\varepsilon$ (i.e taking $\kappa_i = \hat{\kappa}_i \varepsilon^\zeta$ with $\zeta < 1$), then the difference between
the averaged velocities (3.19) tends to zero when $\epsilon$ becomes small. We then recover the model which would have resulted from the derivation with Dirichlet interface conditions with respect to the relative velocity:

$$ (3.22) \qquad (\mathbf{v}_1 - \mathbf{v}_2)|_{z=D\alpha_1} = 0. $$

This condition imposes the continuity of the velocity at the interface and leads to a one-velocity, one-pressure model, which reads:

**THEOREM 3.13.** Assume the same scalings for the viscosity and friction coefficients as in Theorem 3.2 (conditions (3.3)). Let $(\alpha_1, \rho_1, \mathbf{v}_1, p_1, \alpha_2, \rho_2, \mathbf{v}_2, p_2)$ be a solution to the Navier-Stokes system with Navier boundary and Dirichlet interface conditions (2.4)-(2.11), (2.13), (3.22).

Then the asymptotic system satisfies $\langle \mathbf{v}_{1,h} \rangle = \langle \mathbf{v}_{2,h} \rangle$ up to an error of order $\epsilon^{2-\tau}$. Denoting $\langle \mathbf{v}_h \rangle$ this velocity, the functions $(\alpha_1, \alpha_2, \overline{p_1}, \overline{p_2}, \langle \mathbf{v}_h \rangle, \overline{p_1}, \overline{p_2})$ satisfy the following closed system:

$$ (3.23) \qquad \partial_t(\alpha_k \overline{p_k}) + \nabla_h \cdot (\alpha_k \overline{p_k} \langle \mathbf{v}_h \rangle) = 0, $$

$$ (3.24) \qquad \partial_t(\overline{p}(\mathbf{v}_h)) + \nabla_h \cdot (\overline{p}(\mathbf{v}_h) \otimes \langle \mathbf{v}_h \rangle) + \nabla_h(\alpha_1 \overline{p_1} + \alpha_2 \overline{p_2}) = -\delta_{\xi}(\kappa_1 + \kappa_2)\langle \mathbf{v}_h \rangle, $$

$$ (3.25) \qquad \overline{p_k} = p_k(\overline{p_k}), $$

$$ (3.26) \qquad \overline{p_1} = \overline{p_2}, $$

$$ (3.27) \qquad \alpha_1 + \alpha_2 = 1 $$

where $k \in \{1, 2\}$ and $\overline{p} = \alpha_1 \overline{p_1} + \alpha_2 \overline{p_2}$. Same as before, the momentum equation is satisfied up to an error of order $\epsilon^\tau + \epsilon^{2-\tau} + (1-\delta_\xi)\epsilon^{\xi-1}$ and the equality of the pressures is satisfied up to an error of order $\epsilon^\tau$. This system is symmetrizable and therefore hyperbolic.

**Remark 3.14.** The equality of the velocities up to an error of order $\epsilon^{2-\tau}$ stems from the Dirichlet condition associated with the estimate (3.16). Indeed, since the velocities are equal at the interface, and since they don't vary much vertically, we have that

$$ \left\{ \begin{array}{ll} \mathbf{v}_{1,h}(\alpha_1) = \mathbf{v}_{2,h}(\alpha_1), & \forall k \in \{1, 2\}, \\ \mathbf{v}_{k,h}(\alpha_1) = \langle \mathbf{v}_{k,h} \rangle + O\left(\frac{\kappa\epsilon+\epsilon^2}{\mu}\right). & \forall k \in \{3, 4\}, \end{array} \right. $$

Once combined, and developing every coefficient with respect to its scaling,

$$ \langle \mathbf{v}_{1,h} \rangle = \langle \mathbf{v}_{2,h} \rangle + O(\epsilon^{1+\xi-\tau} + \epsilon^{2-\tau}). $$

Since $\xi \ge 1$, the equation is satisfied up to an error of order $\epsilon^{2-\tau}$.

**4. The Navier-Stokes-Fourier case.** This final part extends the results of the previous section to the non-barotropic case. The microscopic model must therefore be completed with two energy equations, and the equations of state must now include the specific internal energies. We reproduce the averaging process on these two energy equations and pay special attention to the heat transfer between the fluids. This will enable us to derive a two-velocity, one-pressure and two-temperature averaged model.

**4.1. The Navier-Stokes-Fourier system.** Let $k \in \{1, 2\}$. Let $\theta_k : \mathbb{R}_+^* \times \Omega_k \to \mathbb{R}$ be the temperature of fluid $k$, $e_k : \mathbb{R}_+^* \times \Omega \to \mathbb{R}$ its specific internal energy $k$ and $E_k = (1/2) |\mathbf{v}_k|^2 + e_k$ its specific total energy. Let $\beta_k$ be the thermal conductivity of
the fluid $k$. It is a constant parameter of the flow that might depend on $\epsilon$, similarly to the viscosities and friction coefficients. The flow is now described by two compressible Navier-Stokes-Fourier systems and two equations of state:

$$ (4.1) \qquad \partial_t \rho_k + \nabla \cdot (\rho_k \mathbf{v}_k) = 0, $$

$$ (4.2) \qquad \partial_t (\rho_k \mathbf{v}_k) + \nabla \cdot (\rho_k \mathbf{v}_k \otimes \mathbf{v}_k) + \nabla p_k = \nabla \cdot \mathbb{S}_k, $$

$$ (4.3) \qquad \partial_t (\rho_k E_k) + \nabla \cdot (\mathbf{v}_k (\rho_k E_k + p_k)) = \nabla \cdot (\mathbb{S}_k \mathbf{v}_k) + \nabla \cdot (\beta_k \nabla \theta_k), $$

$$ (4.4) \qquad p_k = p_k(\rho_k, \rho_k e_k), $$

$$ (4.5) \qquad \theta_k = \theta_k(\rho_k, \rho_k e_k). $$

The pressure $p_k$ and the temperature $\theta_k$ are defined thanks to a complete equation of state which assumes the existence of a specific entropy $s_k$ such that

$$ (4.6) \qquad \theta_k ds_k = p_k d\left(\frac{1}{\rho_k}\right) + de_k. $$

We keep Navier boundary conditions on the top and bottom of the domain (see (2.8) and (2.9)) and the interface conditions (2.11) and (2.12). Additional conditions are required to describe the heat flux at the boundary and interface. We prescribe the following Robin-type conditions, see for instance [11]:

$$ (4.7) \qquad \beta_1 \partial_z \theta_1 |_{z=0} = 0, $$

$$ (4.8) \qquad \beta_2 \partial_z \theta_2 |_{z=D} = 0, $$

$$ (4.9) \qquad -\beta_1 \nabla \theta_1 \cdot \mathbf{n}_{\Gamma_D} |_{z=D\alpha_1} = h_c(\theta_1 - \theta_2) |_{z=D\alpha_1}, $$

$$ (4.10) \qquad -\beta_2 \nabla \theta_2 \cdot \mathbf{n}_{\Gamma_D} |_{z=D\alpha_1} = h_c(\theta_1 - \theta_2) |_{z=D\alpha_1}, $$

where $h_c$ is the contact conductance for the interface. The first two conditions depict the adiabatic nature of the flow: the system is therefore thermically isolated. The last two conditions describe the heat exchange at the interface. More precisely, the heat flux at the interface is continuous and is expressed thanks to (4.9) and (4.10).

**4.2. The averaged model.** The averaging procedure of the previous section is still valid in the Navier-Stokes-Fourier case. The difference lies in the presence of energy equations for both fluids as well as new equations of state for the pressure and the temperature. In the following, we are therefore going to focus on averaging these elements, which works in a similar way than what has been done previously. The detailed averaged model is given below.

**THEOREM 4.1.** Assume the following scalings for the different coefficients in our system:

$$ (4.11) \qquad \forall k \in \{1, 2\}, \qquad \left\{ \begin{array}{lll} \mu_k = \widehat{\mu_k \epsilon^\tau}, & \lambda_k = \widehat{\lambda_k \epsilon^\tau}, & 0 < \tau < 2; \\ \kappa_i = \widehat{\kappa_i \epsilon}, & \kappa_k = \widehat{\kappa_k \epsilon^\xi}, & \xi \ge 1; \\ \beta_k = \widehat{\beta_k \epsilon^\gamma}, & & 0 \le \gamma < 2. \end{array} \right. $$

$$ \text{Define} \qquad \delta_\xi = \begin{cases} 1 & \text{if } \xi = 1, \\ 0 & \text{if } \xi > 1. \end{cases} \quad \text{and} \quad \delta_\gamma = \begin{cases} 1 & \text{if } \gamma = 0, \\ 0 & \text{if } \gamma > 0. \end{cases} $$
Let $(\alpha_1, \rho_1, \mathbf{v}_1, p_1, E_1, \theta_1, \alpha_2, \rho_2, \mathbf{v}_2, p_2, E_2, \theta_2)$ be a strong solution to the Navier-Stokes system with Navier boundary and interface conditions (2.13),(4.1)-(4.10). Assume that

$$ (4.12) \qquad \partial_z \nabla_h \theta_k = O(\epsilon^2/\beta). $$

The functions $(\alpha_1, \overline{p_1}, \langle \mathbf{v}_{1,h} \rangle, \overline{p_1}, \langle E_1 \rangle, \overline{\theta_1}, \alpha_2, \overline{p_2}, \langle \mathbf{v}_{2,h} \rangle, \overline{p_2}, \langle E_2 \rangle, \overline{\theta_2})$ satisfy the following closed system:

$$ (4.13) \qquad \partial_t(\alpha_k \overline{p_k}) + \nabla_h \cdot (\alpha_k \overline{p_k} \langle \mathbf{v}_{k,h} \rangle) = 0, $$

$$ (4.14) \qquad \partial_t(\alpha_k \overline{p_k} \langle \mathbf{v}_{k,h} \rangle) + \nabla_h \cdot (\alpha_k \overline{p_k} \langle \mathbf{v}_{k,h} \rangle \otimes \langle \mathbf{v}_{k,h} \rangle) + \nabla_h(\alpha_k \overline{p_k}) = p_i \nabla_h \alpha_k - \delta_\xi \widehat{\kappa_k} \langle \mathbf{v}_{k,h} \rangle + (-1)^k \widehat{\kappa}_i (\langle \mathbf{v}_{1,h} \rangle - \langle \mathbf{v}_{2,h} \rangle), $$

$$ (4.15) \qquad \partial_t(\alpha_k \overline{p_k} \langle E_{k,h} \rangle) + \nabla_h \cdot (\alpha_k (\overline{p_k} \langle E_{k,h} \rangle + \overline{p_k}) \langle \mathbf{v}_{k,h} \rangle) + p_i \partial_t \alpha_k = -\delta_\xi \widehat{\kappa_k} |\langle \mathbf{v}_{k,h} \rangle|^2 + (-1)^k \widehat{\kappa}_i (\langle \mathbf{v}_{1,h} \rangle - \langle \mathbf{v}_{2,h} \rangle) \cdot \langle \mathbf{v}_{k,h} \rangle + \delta_\gamma (\widehat{\beta_k} \nabla_h \cdot (\alpha_k \nabla_h \overline{\theta_k})) + (-1)^k \widehat{h_c} (\overline{\theta_1} - \overline{\theta_2}), $$

$$ (4.16) \qquad \overline{p_k} = p_k(\overline{p_k}, \overline{p_k} \langle e_k \rangle) $$

$$ (4.17) \qquad \overline{\theta_k} = \theta_k(\overline{p_k}, \overline{p_k} \langle e_k \rangle), $$

$$ (4.18) \qquad \overline{p_1} = \overline{p_2}, $$

$$ (4.19) \qquad \alpha_1 + \alpha_2 = 1, $$

where $k \in \{1, 2\}$ and $\widehat{h_c} = h_c/\epsilon$.

As in the barotropic case, equations (4.14) and (4.18) are satisfied up to an error of order $\epsilon^\tau + \epsilon^{2-\tau} + (1 - \delta_\xi)\epsilon^{\xi-1}$ and $\epsilon^\tau$ respectively. The averaged energy equation (4.15) is satisfied up to an error of order $\epsilon^\tau + \epsilon^{2-\tau} + \epsilon^{2-\gamma} + (1 - \delta_\gamma)\epsilon^\gamma$. The averaged equations of state (4.16)-(4.17) are satisfied up to an error of order $\epsilon^\tau$ and $\epsilon^{2-\gamma}$ respectively.

**Remark 4.2.** Replacing the equations of state (2.6) by the energy equations (4.3) and the equations of state (4.4)-(4.5) does not change the fact that the asymptotic averaged system has only one pressure.

**Remark 4.3.** The scaling of the thermal conductivity can be expressed in terms of the Péclet number $Pe = \text{RePr}$, where $\text{Pr} = c_p\mu/\beta$ is the Prandtl number, $c_p$ is the specific heat and $\text{Re}$ is the Reynolds number. In engineering applications, the Péclet number is often much larger than one, which is consistent with the fact that the Reynolds number is very large in our work. The prescribed scalings for the viscosity and the thermal conductivity enable us to browse different regimes for the Prandtl number.

**4.3. Proof of Theorem 4.1.** The proof of Theorem 4.1 is similar to the one in the barotropic case. In particular, the averaging of the mass conservation and momentum equations is the same. In this proof, we will only detail the specificities of the Navier-Stokes-Fourier case: rescaling the new unknowns and the energy equations, averaging of the equations of state and integrating the energy equations.

**4.3.1. Rescaling of the energy equations and the thermal flux conditions.** Following the procedure of the barotropic case, it is now necessary to rescale and simplify (by neglecting some terms, as in subsection 3.2.1) the energy equation (4.3). We introduce the following characteristic dimensions: $\mathcal{E} = U^2$ for the specific energy, $\Theta$ for the temperature, $B = RU^3L/\Theta$ for the thermal conductivity
and $H = RU^3/\Theta$ for the contact conductance. We then define the following rescaled variables: $\tilde{E} = E/\mathcal{E}$, $\tilde{\theta} = \theta/\Theta$, $\tilde{\beta} = \beta/B$ and $\tilde{h}_c = h_c/H$. Then, we remove of the $\tilde{\cdot}$ notation for the sake of clarity.

The rescaling procedure can be done term by term. Let us start with the left-hand side of the energy equation (4.3). Denoting $E_{k,h} = (1/2)|\mathbf{v}_{k,h}|^2 + e_k$ for $k \in \{1, 2\}$, the rescaled left-hand side of the energy equation reads

$$ \partial_t(\rho_k E_{k,h}) + \nabla \cdot (\mathbf{v}_k(\rho_k E_{k,h} + p_k)) + \varepsilon^2(\partial_t(\rho_k w_k^2/2) + \nabla \cdot (\mathbf{v}_k(\rho_k w_k^2/2))). $$

The terms of order $\varepsilon^2$ are neglected in the rest of the derivation.

Rescaling the thermal diffusion term is straightforward and yields

$$ \beta_k \Delta_h \theta_k + \frac{\beta_k}{\varepsilon^2} \partial_{zz} \theta_k. $$

The viscous term is the most technical one.

**LEMMA 4.4.** *After simplification, the rescaled viscous term in the energy equation (4.3) reads*

$$ (4.20) \qquad \frac{\mu_k}{\varepsilon^2} \partial_z (\mathbf{v}_{k,h} \cdot \partial_z \mathbf{v}_{k,h}) + O(\mu). $$

*Proof.* Developing the viscous term of the energy equation (4.3) yields

$$ \nabla \cdot (S_k \mathbf{v}_k) = \nabla \cdot (\lambda_k (\nabla \cdot \mathbf{v}_k) \mathbf{v}_k) + \mu_k \nabla \cdot (2D_k \mathbf{v}_k) $$

The bulk term does not need to be detailed as the rescaling does not affect its expression. Since

$$ 2D_k \mathbf{v}_k = \begin{pmatrix} 2u \partial_x u + v(\partial_y u + \partial_x v) + w(\partial_z u + \partial_x w) \\ u(\partial_x v + \partial_y u) + 2v \partial_y v + w(\partial_z v + \partial_y w) \\ u(\partial_x w + \partial_z u) + v(\partial_y w + \partial_z v) + 2w \partial_z w \end{pmatrix}, $$

the rescaled version of its divergence reads

$$ \begin{pmatrix} \partial_x \\ \partial_y \\ \frac{1}{\varepsilon} \partial_z \end{pmatrix} \cdot \begin{pmatrix} 2u \partial_x u + v(\partial_y u + \partial_x v) + w(\partial_z u + \varepsilon^2 \partial_x w) \\ u(\partial_x v + \partial_y u) + 2v \partial_y v + w(\partial_z v + \varepsilon^2 \partial_y w) \\ u(\varepsilon \partial_x w + \frac{1}{\varepsilon} \partial_z u) + v(\varepsilon \partial_y w + \frac{1}{\varepsilon} \partial_z v) + 2\varepsilon w \partial_z w \end{pmatrix} \\ = \nabla_h \cdot ((\nabla_h \otimes \mathbf{v}_{k,h} + (\nabla_h \otimes \mathbf{v}_{k,h})^T) \mathbf{v}_{k,h}) + \nabla_h \cdot (w_k \partial_z \mathbf{v}_{k,h} + \varepsilon^2 w_k \nabla_h w_k) \\ + \partial_z (\mathbf{v}_{k,h} \cdot \nabla_h w_k) + \frac{1}{\varepsilon^2} \partial_z (\mathbf{v}_{k,h} \cdot \partial_z \mathbf{v}_{k,h}) + 2w_k \partial_z w_k. $$

Multiplying the previous term by $\mu_k$ and adding the bulk term yields the following expression for the rescaled viscous term:

$$ \lambda_k \nabla \cdot ((\nabla \cdot \mathbf{v}_k)\mathbf{v}_k) + \mu_k \nabla_h \cdot ((\nabla_h \otimes \mathbf{v}_{k,h} + (\nabla_h \otimes \mathbf{v}_{k,h})^T)\mathbf{v}_{k,h}) \\ + \mu_k \nabla_h \cdot (w_k \partial_z \mathbf{v}_{k,h} + \varepsilon^2 w_k \nabla_h w_k) + \mu_k \partial_z (\mathbf{v}_{k,h} \cdot \nabla_h w_k) \\ + \mu_k \frac{1}{\varepsilon^2} \partial_z (\mathbf{v}_{k,h} \cdot \partial_z \mathbf{v}_{k,h}) + 2\mu_k w_k \partial_z w_k. $$

The prevailing term in the previous expression is of order $\mu/\varepsilon^2$. The rest of the terms are of order $\mu$. Knowing that the use of Navier interface conditions cause an approximation error of order $\mu$, we neglect those terms right away to ensure simpler equations. Hence (4.20). □
Finally, the rescaled version of the energy equation reads

$$ (4.21) \quad \begin{aligned} \partial_t (\rho_k E_{k,h}) + \nabla \cdot (\mathbf{v}_k (\rho_k E_{k,h} + p_k)) = & \frac{\mu_k}{\varepsilon^2} \partial_z (\mathbf{v}_{k,h} \cdot \partial_z \mathbf{v}_{k,h}) \\ & + \beta_k \Delta_h \theta_k + \frac{\beta_k}{\varepsilon^2} \partial_{zz} \theta_k + O(\mu + \varepsilon^2). \end{aligned} $$

The rescaled boundary and interface conditions (4.7)-(4.10) are

$$ (4.22) \qquad \beta_1 \partial_z \theta_1 \Big|_{z=0} = 0, $$

$$ (4.23) \qquad \beta_2 \partial_z \theta_2 \Big|_{z=1} = 0, $$

$$ (4.24) \qquad \left. \left( \varepsilon \beta_1 \nabla_h \alpha_1 \cdot \nabla_h \theta_1 - \frac{\beta_1}{\varepsilon} \partial_z \theta_1 \right) \right|_{z=\alpha_1} = h_c (\theta_1 - \theta_2) \Big|_{z=\alpha_1}, $$

$$ (4.25) \qquad \left. \left( \varepsilon \beta_2 \nabla_h \alpha_1 \cdot \nabla_h \theta_2 - \frac{\beta_2}{\varepsilon} \partial_z \theta_2 \right) \right|_{z=\alpha_1} = h_c (\theta_1 - \theta_2) \Big|_{z=\alpha_1}. $$

Remark 4.5. Multiplying equation (4.21) by $\varepsilon^2$, we get that for $k \in \{1, 2\}$,

$$ \beta_k \partial_{zz} \theta_k(z) = -\mu_k \partial_z (\mathbf{v}_{k,h} \cdot \partial_z \mathbf{v}_{k,h}) + O(\varepsilon^2 + \mu\varepsilon^2 + \varepsilon^4). $$

Since $\mu_k \partial_z (\mathbf{v}_{k,h} \cdot \partial_z \mathbf{v}_{k,h}) = \mu_k |\partial_z \mathbf{v}_{k,h}|^2 + \mu_k \mathbf{v}_{k,h} \cdot \partial_{zz} \mathbf{v}_{k,h}$, we can use the estimates (3.12) and (3.18) and remember the prescribed scalings for the friction and the viscosity coefficients to get that

$$ (4.26) \qquad \beta_k \partial_{zz} \theta_k(z) = O(\varepsilon^2) \quad \text{in } \Omega_k. $$

Using the boundary conditions (4.22) and (4.23), we get that for $k \in \{1, 2\}$,

$$ (4.27) \qquad \beta_k \partial_z \theta_k(z) = O(\varepsilon^2) \quad \text{in } \Omega_k. $$

Therefore, for $k \in \{1, 2\}$,

$$ (4.28) \qquad \theta_k(z) = \bar{\theta}_k + O(\varepsilon^2/\beta) \quad \text{in } \Omega_k. $$

These estimates along with the interface conditions (4.24) and (4.25) tell us that

$$ (4.29) \qquad h_c(\theta_1 - \theta_2) \Big|_{z=\alpha_1} = O(\varepsilon). $$

### 4.3.2. Integration of the energy equations.
Let $k = 1$. Since the integral is linear, we are going to detail the averaging of the different terms of the energy equation (4.21).

First, we integrate the left-hand side of the equation between $z = 0$ and $z = \alpha_1$. Similarly to section 2.3, we use the fraction equation (2.18) and the impermeability condition (2.22) to get

$$ \int_0^{\alpha_1} \left( \partial_t (\rho_1 E_{1,h}) + \nabla \cdot (\mathbf{v}_1 (\rho_1 E_{1,h} + p_1)) \right) dz = \partial_t (\alpha_1 \overline{\rho_1} \langle E_{1,h} \rangle) \\ + \nabla_h \cdot \left( \alpha_1 (\overline{\rho_1} \langle E_{1,h} \mathbf{v}_{1,h} \rangle + \overline{p_1} \overline{\mathbf{v}_{1,h}}) \right) - p_1(\alpha_1) \nabla_h \alpha_1 + p_1(\alpha_1) w_1(\alpha_1). $$

Using once again the fraction equation (2.18),

$$ \int_0^{\alpha_1} \left( \partial_t (\rho_1 E_{1,h}) + \nabla \cdot (\mathbf{v}_1 (\rho_1 E_{1,h} + p_1)) \right) dz = \partial_t (\alpha_1 \overline{\rho_1} \langle E_{1,h} \rangle) \\ + \nabla_h \cdot \left( \alpha_1 (\overline{\rho_1} \langle E_{1,h} \mathbf{v}_{1,h} \rangle + \overline{p_1} \overline{\mathbf{v}_{1,h}}) \right) + p_1(\alpha_1) \partial_t \alpha_1. $$
Replacing the velocities by their Favre average thanks to (3.16), we get

$$ \int_0^{\alpha_1} (\partial_t(\rho_1 E_{1,h}) + \nabla \cdot (\mathbf{v}_1(\rho_1 E_{1,h} + p_1))) dz = \partial_t(\alpha_1 \overline{p_1} \langle E_{1,h} \rangle) \\ + \nabla_h \cdot (\alpha_1 (\overline{p_1} \langle E_{1,h} \rangle + \overline{p_1}) \langle \mathbf{v}_{1,h} \rangle) + p_1(\alpha_1) \partial_t \alpha_1 + O\left(\frac{\varepsilon\kappa + \varepsilon^2}{\mu}\right). $$

Next, we average the viscous term. Using the Navier interface condition (3.11) and the bottom Navier condition (2.22), we get

$$ \int_0^{\alpha_1} \frac{\mu_1}{\varepsilon^2} \partial_z (\mathbf{v}_{1,h} \cdot \partial_z \mathbf{v}_{1,h}) dz = \frac{\mu_1}{\varepsilon^2} (\mathbf{v}_{1,h}(\alpha_1) \cdot \partial_z \mathbf{v}_{1,h}(\alpha_1) - \mathbf{v}_{1,h}(0) \cdot \partial_z \mathbf{v}_{1,h}(0)) \\ = -\frac{\kappa_i}{\varepsilon} (\mathbf{v}_{1,h}(\alpha_1) - \mathbf{v}_{2,h}(\alpha_1)) \cdot \mathbf{v}_{1,h}(\alpha_1) - \frac{\kappa_i}{\varepsilon} |\mathbf{v}_{1,h}(0)|^2 + O(\mu + \kappa_i \varepsilon). $$

Replacing the pointwise values of the velocities by their Favre average using (3.16), we get

$$ \int_0^{\alpha_1} \frac{\mu_1}{\varepsilon^2} \partial_z (\mathbf{v}_{1,h} \cdot \partial_z \mathbf{v}_{1,h}) dz = -\frac{\kappa_i}{\varepsilon} (\langle \mathbf{v}_{1,h} \rangle - \langle \mathbf{v}_{2,h} \rangle) \cdot \langle \mathbf{v}_{1,h} \rangle - \frac{\kappa_i}{\varepsilon} |\langle \mathbf{v}_{1,h} \rangle|^2 \\ + O\left(\mu + \kappa_i \varepsilon + \frac{(\kappa + \kappa_i)(\kappa + \varepsilon)}{\mu}\right). $$

Last, integrate the thermal diffusion term. Using the adiabatic boundary condition (4.7) and the interface condition (4.24), we get

$$ \int_0^{\alpha_1} \left( \beta_1 \Delta_h \theta_1 + \frac{\beta_1}{\varepsilon^2} \partial_{zz} \theta_1 \right) dz = \beta_1 \nabla_h \cdot (\alpha_1 \overline{\nabla_h \theta_1}) \\ - \beta_1 \nabla_h \alpha_1 \cdot \nabla_h \theta_1 (\alpha_1) + \frac{\beta_1}{\varepsilon^2} \partial_z \theta_1 (\alpha_1), \\ = \beta_1 \nabla_h \cdot (\alpha_1 \overline{\nabla_h \theta_1}) - \frac{h_c}{\varepsilon} (\theta_1(\alpha_1) - \theta_2(\alpha_1)). $$

Approximating the pointwise values of the temperatures by their averages thanks to (4.28), we get

$$ \int_0^{\alpha_1} \left( \beta_1 \Delta_h \theta_1 + \frac{\beta_1}{\varepsilon^2} \partial_{zz} \theta_1 \right) dz = \beta_1 \nabla_h \cdot (\alpha_1 \overline{\nabla_h \theta_1}) - \frac{h_c}{\varepsilon} (\overline{\theta_1} - \overline{\theta_2}) + O(\varepsilon^2/\beta). $$

The term $\overline{\nabla_h \theta_1}$ is unsatisfactory since it introduces a new unknown. To solve this problem, two strategies are possible. The first one is to express this term in terms of the averaged temperature $\theta_1$ using the lemma below.

LEMMA 4.6. Assume that for $k \in \{1, 2\}$

$$ \partial_z \nabla_h \theta_k = O(\varepsilon^2/\beta) \quad \text{in } \Omega_k. $$

Then, for $k \in \{1, 2\}$,

$$ (4.30) \qquad \nabla_h \cdot (\alpha_k \overline{\nabla_h \theta_k}) = \nabla_h \cdot (\alpha_k \overline{\nabla_h \theta_k}) + O(\varepsilon^2/\beta). $$

Proof. The proof is done with $k = 1$. Since

$$ (4.31) \qquad \alpha_1 \nabla_h \overline{\theta_1} = \alpha_1 \overline{\nabla_h \theta_1} + \nabla_h \alpha_1 (\theta_1(\alpha_1) - \overline{\theta_1}), $$
taking the horizontal divergence of the previous equation yields

$$ \nabla_h \cdot (\alpha_1 \nabla_h \overline{\theta_1}) = \nabla_h \cdot ((\alpha_1 \overline{\nabla_h \theta_1})) + \Delta_h \alpha_1 (\theta_1(\alpha_1) - \overline{\theta_1}) \\ + \nabla_h \alpha_1 \cdot ((\nabla_h \theta_1)(\alpha_1) + \partial_z \theta_1(\alpha_1) \nabla_h \alpha_1 - \nabla_h \overline{\theta_1}). $$

Using the estimate (4.28),

$$ \nabla_h \cdot (\alpha_1 \nabla_h \overline{\theta_1}) = \nabla_h \cdot ((\alpha_1 \overline{\nabla_h \theta_1})) + \nabla_h \alpha_1 \cdot ((\nabla_h \theta_1)(\alpha_1) - \overline{\nabla_h \theta_1}) + O(\varepsilon^2/\beta). $$

$$ (4.32) \quad \text{Reusing (4.31) and estimate (4.28),} \qquad \nabla_h \overline{\theta_1} = \overline{\nabla_h \theta_1} + O(\varepsilon^2/\beta). $$

Therefore,

$$ \nabla_h \cdot (\alpha_1 \nabla_h \overline{\theta_1}) = \nabla_h \cdot (\alpha_1 \overline{\nabla_h \theta_1}) + \nabla_h \alpha_1 \cdot ((\nabla_h \theta_1)(\alpha_1) - \overline{\nabla_h \theta_1}) + O(\varepsilon^2/\beta). $$

Finally, thanks to the assumption (4.12), we can obtain equation (4.30).

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

The second way to solve the introduction of the unknown $\overline{\nabla_h \theta_1}$ is to scale the thermal conductivities as $\beta_k = \widehat{\beta}_k \varepsilon^\gamma$ for $k \in \{1, 2\}$ with $0 < \gamma < 2$. This allows to get rid of the heat diffusion term up to an error of order $\varepsilon^\gamma$. Note that in this case, no assumption about the vertical variation of $\nabla_h \theta_k$ is required. To synthesize the results, we assume that the thermal conductivities are scaled with $0 \le \gamma < 2$ as in (4.11) and that the hypothesis (4.12) is satisfied.

Assembling all the terms, we get the averaged energy equation

$$ (4.33) \quad \begin{aligned} \partial_t(\alpha_1 \overline{p_1} \langle E_{1,h} \rangle) + \nabla_h \cdot (\alpha_1 (\overline{p_1} \langle E_{1,h} \rangle + \overline{p_1}) \langle \mathbf{v}_{1,h} \rangle) + p_1(\alpha_1) \partial_t \alpha_1 = & - \frac{\kappa_1}{\varepsilon} |\langle \mathbf{v}_{1,h} \rangle|^2 \\ & + \frac{\kappa_i}{\varepsilon} (\langle \mathbf{v}_{1,h} \rangle - \langle \mathbf{v}_{2,h} \rangle) \cdot \langle \mathbf{v}_{1,h} \rangle + \beta_1 \nabla_h (\alpha_1 \nabla_h \overline{\theta_1}) - \frac{h_c}{\varepsilon} (\langle \theta_1 \rangle - \langle \theta_2 \rangle) + \mathcal{R}, \end{aligned} $$

$$ \text{where} \qquad \mathcal{R} = O\left(\mu + \kappa_i \varepsilon + \varepsilon^2 + \frac{(\kappa + \kappa_i + \varepsilon)(\kappa + \varepsilon)}{\mu} + \frac{\varepsilon^2}{\beta}\right). $$

Since we have that the averaged pressures are equal up to an error of order $\mu$, they are still indistinguishable in the averaged energy equation. Moreover, similarly to the momentum equation, the pointwise value $p_1(\alpha_1)$ can be replaced by any convex combination of the averaged pressures, that we will denote $p_i$.

Using the scalings (4.11), the averaged energy equation (4.15) holds up to an error of order $\varepsilon^\tau + \varepsilon^{2-\tau} + \varepsilon^{2-\gamma} + (1 - \delta_\gamma)\varepsilon^\gamma$.

**Remark 4.7.** The scaling for the contact conductance $h_c$ has been prescribed in order to get a two-temperature model. Note that the only constraint that we have is given by (4.29). Choosing a contact conductance that decreases slower than $\varepsilon$ would therefore imply that the difference of temperatures at the interface decreases to zero when $\varepsilon$ becomes very small, resulting in a one-temperature averaged model (since the temperatures do not vary sharply in the vertical direction).
### 4.3.3. Obtaining the averaged equations of state.
The equations of state are averaged in a similar way than in the barotropic case. However, now the pressures and the temperatures depend on two variables, so the proof needs a little adaptation. Let us prove the result for the pressure $p_1$ (the proof for $p_2$ and the temperatures are similar).

Recall that thanks to the asymptotic vertical momentum equation, we have that $\partial_z p_1(\rho_1, \rho_1 e_1) = O(\mu)$. Let $z^* \in (0, \alpha_1)$ such that $\rho_1(z^*) = \bar{\rho}_1$. Thanks to a Taylor-Lagrange expansion, we have that

$$ \forall z \in (0, \alpha_1), \quad p_1(\rho_1(z), \rho_1 e_1(z)) = p_1(\rho_1(z^*), \rho_1 e_1(z)) + O(\mu). $$

Then, consider $z' \in (0, \alpha_1)$ such that $\rho_1 e_1(z') = \bar{\rho}_1 e_1 = \bar{\rho}_1 \langle e_1 \rangle$. A second Taylor-Lagrange expansion yields

$$ \forall z \in (0, \alpha_1), \quad p_1(\rho_1(z^*), \rho_1 e_1(z)) = p_1(\rho_1(z^*), \rho_1 e_1(z')) + O(\mu). $$

$$ \text{Therefore, } \forall z \in (0, \alpha_1), \quad p_1(\rho_1(z), \rho_1 e_1(z)) = p_1(\bar{\rho}_1, \bar{\rho}_1 \langle e_1 \rangle) + O(\mu). $$

Averaging the previous equation with respect to $z$ and dropping the error yields (4.16).

The procedure is the same for the temperature, but using the estimate (4.27), we get that the approximation error is of order $\epsilon^2/\beta$.

*Acknowledgements.* The two last authors would like to thank Meriem Bahhi who nurtured the preliminary reflexions of this work.

## REFERENCES

[1] M. R. BAER AND J. W. NUNZIATO, *A Two-Phase Mixture Theory for the Deflagration-to-Detonation Transition (DDT) in Reactive Granular Materials*, International Journal of Multiphase Flow, 12 (1986), pp. 861–889, https://doi.org/10.1016/0301-9322(86)90033-9.

[2] D. BESTION, *The physical closure laws in the CATHARE code*, Nuclear Engineering and Design, 124 (1990), pp. 229–245, https://doi.org/10.1016/0029-5493(90)90294-8.

[3] D. BRESCH, C. BURTEA, AND F. LAGOUTIÈRE, *Mathematical justification of a compressible bifluid system with different pressure laws: a continuous approach*, Applicable Analysis, 101 (2022), pp. 4235–4266, https://doi.org/10.1080/00036811.2022.2103679.

[4] D. BRESCH, B. DESJARDINS, J.-M. GHIDAGLIA, AND E. GRENIER, *Global Weak Solutions to a Generic Two-Fluid Model*, Archive for Rational Mechanics and Analysis, 196 (2010), pp. 599–629, https://doi.org/10.1007/s00205-009-0261-6.

[5] D. BRESCH AND M. HILLAIRET, *Note on the Derivation of Multi-Component Flow Systems*, Proceedings of the American Mathematical Society, 143 (2015), pp. 3429–3443, https://doi.org/10.1090/proc/12614.

[6] D. BRESCH AND M. HILLAIRET, *A compressible multifluid system with new physical relaxation terms*, Annales scientifiques de l'École normale supérieure, 52 (2019), pp. 255–295, https://doi.org/10.24033/asens.2387.

[7] D. BRESCH AND X.-D. HUANG, *A Multi-Fluid Compressible System as the Limit of Weak Solutions of the Isentropic Compressible Navier-Stokes Equations*, Archive for Rational Mechanics and Analysis, 201 (2011), pp. 647–680, https://doi.org/10.1007/s00205-011-0400-8.

[8] F. COQUEL, J.-M. HÉRARD, K. SALEH, AND N. SEGUIN, *Two properties of two-velocity two-pressure models for two-phase flows*, Communications in Mathematical Sciences, 12 (2014), https://hal.science/hal-00788902/file/CHSS-1.pdf (accessed 2025-05-27).

[9] C. DEMAY AND J.-M. HÉRARD, *A compressible two-layer model for transient gas-liquid flows in pipes*, Continuum Mechanics and Thermodynamics, 29 (2017), pp. 385–410, https://doi.org/10.1007/s00161-016-0531-0.

[10] D. A. DREW AND S. L. PASSMAN, *Theory of Multicomponent Fluids*, vol. 135 of Applied Mathematical Sciences, Springer New York, New York, NY, 1999, https://doi.org/10.1007/b97678.

[11] D. W. HAHN AND M. N. OZISIK, *Heat Conduction*, Wiley, Hoboken, N.J, 3rd ed ed., 2012, https://doi.org/10.1002/9781118411285.
[12] M. HILLAIRET, *On Baer-Nunziato Multiphase Flow Models*, ESAIM: Proceedings and Surveys, 66 (2019), pp. 61–83, https://doi.org/10.1051/proc/201966004.

[13] M. HILLAIRET, H. MATHIS, AND N. SEGUIN, *Analysis of compressible bubbly flows. Part II: Derivation of a macroscopic model*, ESAIM: Mathematical Modelling and Numerical Analysis, 57 (2023), pp. 2865–2906, https://doi.org/10.1051/m2an/2023046.

[14] M. ISHI AND T. HIBIKI, *Thermo-Fluid Dynamics of Two-Phase Flow*, Springer, New York, 2nd ed., 2011.

[15] M. NDJINGA, *Influence of interfacial pressure on the hyperbolicity of the two-fluid model*, Comptes Rendus. Mathématique, 344 (2007), pp. 407–412, https://doi.org/10.1016/j.crma.2007.02.006.

[16] J. D. RAMSHAW AND J. A. TRAPP, *Characteristics, stability, and short-wavelength phenomena in two-phase flow equation systems*, Nuclear Science and Engineering, 66 (1978), pp. 93–102, https://doi.org/10.13182/NSE78-A15191.

[17] V. H. RANSOM AND D. L. HICKS, *Hyperbolic two-pressure models for two-phase flow*, Journal of Computational Physics, 53 (1984), pp. 124–151, https://doi.org/10.1016/0021-9991(84)90056-1.

[18] V. H. RANSOM AND D. L. HICKS, *Hyperbolic two-pressure models for two-phase flow revisited*, Journal of Computational Physics, 75 (1988), pp. 498–504, https://doi.org/10.1016/0021-9991(88)90125-8.

[19] H. B. STEWART, *Stability of two-phase flow calculation using two-fluid models*, Journal of Computational Physics, 33 (1979), pp. 259–270, https://doi.org/10.1016/0021-9991(79)90020-2.

[20] H. B. STEWART AND B. WENDROFF, *Two-phase flow: Models and methods*, Journal of Computational Physics, 56 (1984), pp. 363–409, https://doi.org/10.1016/0021-9991(84)90103-7.