# Study of interface coupling in three-layer viscous fluid systems

Antoine Simon, Jorge César Brändle de Motta, Christophe Dumouchel, Marie-Charlotte Renoult

► To cite this version:

Antoine Simon, Jorge César Brändle de Motta, Christophe Dumouchel, Marie-Charlotte Renoult. Study of interface coupling in three-layer viscous fluid systems. Physica D: Nonlinear Phenomena, 2025, 477, pp.134670. 10.1016/j.physd.2025.134670. hal-05121586

HAL Id: hal-05121586

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

Submitted on 24 Nov 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Study of interface coupling in three-layer viscous fluid systems

A. Simon, J.C. Brändle de Motta<sup>1D</sup>, C. Dumouchel<sup>1D</sup>, M.-C. Renoult<sup>1D</sup>*

INSA Rouen Normandie, Univ Rouen Normandie, CNRS, Normandie Univ, CORIA UMR 6614, Rouen, F-76000, France

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

## ARTICLE INFO

Communicated by Dmitry Pelinovsky

**Keywords:**

Three-layer fluid system

Rayleigh–Taylor instability

Gravity wave

Linear stability analysis

Direct numerical simulations

## ABSTRACT

Linear stability analysis is performed to study the effect of viscosity and surface tension on a system composed of three superimposed immiscible incompressible Newtonian fluids under the gravity field. The eigenvalue problem is formulated in a general way to highlight an analytical solution and other solutions to be determined numerically for a given set of physical parameters. The behaviour of these solutions is analysed using a coupling parameter equal to the product between the thickness of the middle layer and the wave number of the disturbance. When this parameter is large enough, the study of this three-layer fluid system is reduced to the study of two-layer fluid subsystems. The solutions were determined for a three-layer fluid system of interest with a gravitationally unstable interface at the bottom and a gravitationally stable interface at the top to highlight respectively the Rayleigh–Taylor instability and a gravity wave as well as the coupling between these two phenomena. The temporal evolution of the physical quantities is obtained by solving the initial value problem. For this purpose, single-mode disturbances at the two interfaces are imposed as initial conditions. Direct numerical simulations performed with an in-house code initialized by the interface positions and the fluid velocities allow us to compare their temporal evolutions. Results show an excellent agreement with the theory until the amplitude of one interface becomes too large in the same way as the nonlinear effects of the flow. Three main regimes are observed depending on which solutions initially influence the flow the most.

## 1. Introduction

When two fluids are superimposed in the gravity field, two different types of flow may occur. If the top fluid is denser than the bottom one, Rayleigh–Taylor instability may develop at the interface between them. This instability gives rise to distinct digitations: spikes from the top to the bottom layer and bubbles from the bottom to the top layer. Whereas, if the bottom fluid is denser than the top one, a gravity wave whose amplitude will decrease with time could be generated.

Such two-layer fluid system was initially studied by Lord Rayleigh in the late 19th century, although he did not take into account the viscosity of the fluids and the surface tension between them [1]. He carried out linear stability analysis that he originally used to determine the law of falling away from unstable equilibrium in cases with a plane or cylindrical surface of separation [2]. This method was inspired by the work of Sir W. Thomson who determined the conditions under which a level surface of water is rendered unstable given the influence of wind on waves [3]. It consists of determining a state of equilibrium (here, an equilibrium position: the surface between two fluids considered flat) and then introducing a small disturbance. According to Fourier's theorem, this disturbance can be decomposed as a sum of waves that can be considered separately. In Cartesian coordinates each wave can be expressed under the shape $A \exp(jk \cdot x + \omega t)$, with

A being the wave amplitude, j the imaginary unit, k the wave vector, $x = (x, y, z)$ the position vector, t the time and $\omega$ the complex time coefficient comprising the growth rate of the wave, $\Re(\omega)$, and the pulsation of the wave, $\Im(\omega)$. $\Re(\omega)$ and $\Im(\omega)$ are the real and imaginary part of the time coefficient, respectively. Thus, a wave is said to be unstable if $\Re(\omega) > 0$; otherwise, it is stable. The dispersion relation that gives the relation between the time coefficient and the wave vector allows us to determine the stability condition of a system. Lord Rayleigh determined the following dispersion relation for his studied case:

$$ \omega^2 = \frac{\rho_2 - \rho_1}{\rho_2 + \rho_1} gk, \qquad (1) $$

where $\rho_1$ and $\rho_2$ are the density of the bottom and top layers, respectively, g is the norm of the gravity field and k is the wave vector norm. The solutions of this equation are the eigenvalues of this problem. Determining the spectrum of the system is called the resolution of the eigenvalue problem, with each eigenvalue corresponding to a normal mode. The system is said to be unstable if at least one of its normal modes is unstable. This means that for the case studied by Lord Rayleigh, the system is always unstable when the top layer is denser than the bottom layer.

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: renoultm@coria.fr (M.-C. Renoult).</footnote>
To go further in the mathematical description of the physical problem, it is important to consider the surface tension between fluids as well as the viscosity of the fluids.

Firstly, surface tension can be seen as an effect of the density ratio between the two fluids [4]. One consequence is a jump in pressure: on the concave side, the pressure exceeds that on the convex side by the Laplace pressure $\gamma_{12}(r_1^{-1} + r_2^{-1})$, where $\gamma_{12}$ is the surface tension and $r_1$ and $r_2$ are the radii of curvature of any two sections normal to the surface and to each other. Taking account of this pressure jump in the previous theory allows a new dispersion relation [5,6]:

$$ \omega^2 = \frac{\rho_2 - \rho_1}{\rho_2 + \rho_1} g k - \frac{\gamma_{12}}{\rho_2 + \rho_1} k^3. \quad (2) $$

A denser top layer compared with the bottom layer is no longer a sufficient condition for instability. Here, the wave number must also be less than a critical value called the cutoff wave number $k_c = ((\rho_2 - \rho_1)g/\gamma_{12})^{1/2}$.

Secondly, in the absence of viscosity, the tangential component of the velocity is bounded (note that the normal component of the velocity is continuous due to the conservation of mass). Taking viscosity into consideration removes this discontinuity [7]. While considering the viscosity of the fluids, the dispersion relation can no longer be solved analytically without some approximations. For example, disregarding the surface tension, Harrison [5] found that for $k\sqrt{v/\omega} \ll 1$, $v$ being the largest kinematic viscosity, the eigenvalue for the gravitationally unstable case ($\rho_2 > \rho_1$) can be approximated as the sum of different order terms with respect to $k\sqrt{v/\omega}$ (Fig. 1). The dominant term $\omega_1$ is the same as that found by Lord Rayleigh when $\rho_2 > \rho_1$:

$$ \omega_1 = \left( \frac{\rho_2 - \rho_1}{\rho_2 + \rho_1} g k \right)^{1/2}. \quad (3) $$

The second-order approximation term $\omega_2$ is given by:

$$ \omega_2 = \omega_1 - \omega_1^2 \frac{2k\rho_1\rho_2(v_1v_2)^{1/2}}{(\rho_2 + \rho_1)(\rho_2v_2^{1/2} + \rho_1v_1^{1/2})}, \quad (4) $$

where $v_1$ and $v_2$ are respectively the kinematic viscosity of the bottom and top layers.

Chandrasekhar [6] studied the limit of this case when $k\sqrt{v/\omega} \gg 1$ considering that the kinematic viscosities of the two fluids are the same. He found that the corresponding eigenvalue, $\omega_\infty$, becomes proportional to $k^{-1}$ (Fig. 1). He also showed that viscosity has no effect on the cutoff wave number. Thus, some analytical eigenvalues can be found for extreme values of $k\sqrt{v/\omega}$, but for intermediate values of $k\sqrt{v/\omega}$, they have to be determined numerically.

The study of this fluid system with two superimposed layers provides a better understanding of any fluid system in which one fluid is accelerated towards another fluid in contact with it [8]. This concerns natural phenomena such as interstellar medium and galaxy clusters, the inner edge of an astrophysical disk around a central black hole, accretion on the magnetospheres of neutron stars, the magnetospheres of Jupiter and Saturn and the sun's filamentary structure [9]. This also concerns artificial phenomena such as premixed combustion, reactive systems without premix, heat exchangers and sprays in internal combustion chambers, industrial coating with thin liquid films [9], inertial confinement fusion [10,11], Hall-Héroult cells, electric arc furnaces [12] and liquid metal batteries [13]. However, in the majority of cases, these systems have at least three fluid layers.

Three-layer fluid systems were studied by Mikaelian [14] without taking into account the viscosity of the fluids and the surface tension between them. For a given wavelength, he found that when the thickness of the middle layer tends to infinity, the eigenvalues of the system are identical to those obtained if each interface were treated independently; the interfaces are thus said to be decoupled. However, when the thickness of the middle layer tends to zero, the eigenvalues of the system are identical to those obtained if only one interface were treated between the two extreme layers plus a zero eigenvalue. As

<figure><img src="image_4.png" /><figcaption>Fig. 1. Dispersion relation of Chandrasekhar [6] in black for $\gamma_{12} = 0$ kg s$^{-2}$, $\rho_1 = 773$ kg m$^{-3}$, $\rho_2 = 1398$ kg m$^{-3}$ and $v_1 = v_2 = 10^{-5}$ m$^2$ s$^{-1}$ with analytical approximations for extreme values of $k\sqrt{v/\omega}$.</figcaption></figure>

expected, this work revealed that the thickness of the middle layer is an important parameter to study the coupling between interfaces for three-layer fluid systems. Sneyd [12] subsequently added an additional volume force: the Lorentz force due to an electric current flow in the fluids in order to determine the conditions under which this force can destabilize a three-layer fluid system with two gravitationally stable interfaces. Parhi et al. [15] studied the effect of both the surface tension and viscosity on three-layer fluid systems but without determining all the eigenvalues in the general case. The theoretical works cited previously are summarized in Table 1.

A few experimental works have been carried out in support of these more recent theories about three-layer fluid systems including the work of Jacobs et al. [16] for miscible fluids. Considering an unstable interface at the top and a stable interface at the bottom, their experiments show that the instability follows a self-similar development when the densities of the extreme layers are the same and the erosion of the bottom layer due to the turbulence generated by the upper interface is relatively small. This erosion decreases with respect to the density of the top layer and increases with respect to the density of the bottom layer. When the density of the bottom layer becomes smaller than the density of the middle layer the bottom interface becomes unstable and the erosion is similar to the one due to the Rayleigh-Taylor instability for the two-layer case. Adkins et al. [17] provided experimental data for immiscible fluids and compared their results with those of Mikaelian [14]. They studied two systems composed of a gravitationally unstable interface at the bottom and a gravitationally stable interface at the top. The growth rate of the bottom interface was measured by assuming that its flow was that of a Rayleigh-Taylor instability, whereas the flow of the upper interface was that of a gravity wave. These data then revealed some discrepancies with Mikaelian's theory. These discrepancies can be explained by the absence of viscosity and surface tension in Mikaelian's theory and by the fact that the coupling between the interfaces was not considered in the analysis of Adkins et al. An interesting extension of Mikaelian's theory is the current goal to consider the viscosity of the fluids and the surface tension between them. Moreover, in order to study the coupling between the two interfaces, in addition to solving the eigenvalue problem, the initial value problem needs to be solved. It consists of imposing initial conditions on the physical quantities; for instance, we can impose pure deformation (i.e., a deformed interface with a zero velocity) or pure
<table><caption>**Table 1**
Summary table of landmark theoretical works carried out on two- or three-layer fluid systems.</caption><thead><tr><th>Reference</th><th>Additional volume force</th><th>Viscosity</th><th>Surface tension</th><th>Layers</th><th>Resolution of the eigenvalue problem</th><th>Resolution of the initial value problem</th></tr></thead><tbody><tr><td>Rayleigh [1]</td><td>No</td><td>No</td><td>No</td><td>2</td><td>Yes</td><td>No</td></tr><tr><td>Harrison [5]</td><td>No</td><td>Yes</td><td>No</td><td>2</td><td>Yes</td><td>No</td></tr><tr><td>Chandrasekhar [6]</td><td>No</td><td>Yes</td><td>Yes</td><td>2</td><td>Yes</td><td>No</td></tr><tr><td>Proisperetti [19]</td><td>No</td><td>Yes</td><td>Yes</td><td>2</td><td></td><td>Yes</td></tr><tr><td>Mikaelian [14]</td><td>No</td><td>No</td><td>No</td><td>3</td><td>Yes</td><td>No</td></tr><tr><td>Sneyd [12]</td><td>Lorentz force</td><td>No</td><td>No</td><td>3</td><td>Yes</td><td>No</td></tr><tr><td>Parhi et al. [15]</td><td>No</td><td>Yes</td><td>Yes</td><td>3</td><td>No</td><td>No</td></tr><tr><td>Mikaelian [23]</td><td>No</td><td>Yes</td><td>Yes</td><td>2 (finite)</td><td>No</td><td>No</td></tr><tr><td>Vartdal et al. [22]</td><td>No</td><td>Yes</td><td>No</td><td>3</td><td>No</td><td>Yes (zero velocity)</td></tr><tr><td>Present work</td><td>No</td><td>Yes</td><td>Yes</td><td>3</td><td>Yes</td><td>Yes</td></tr></tbody></table>

impulse (i.e., a flat interface with a non-zero velocity) conditions. The number of conditions depends on the number of eigenvalues [18].

Solving the initial value problem has been used in the past and notably by Prosperetti [19]. He used it to study the damping of small-amplitude wave motion on the free surface between two fluids which has often been approximated from the least damped of the normal modes. This method allows to determine the temporal evolution of the interface initially in addition to its single-mode growth or decay determined from the normal mode approximation [20]. When the initial amplitude is too large, nonlinear effects do not allow solving the initial value problem to correctly determine the temporal evolution of the interface [21]. This method has also been used in the case of a three-layer system by Vartdal et al. [22] by neglecting the surface tension and considering initially perturbed interfaces but with a zero initial velocity within the fluids. In the present work we consider the initial velocity generated by the initial deformation of the interfaces.

The main goal of this paper is to study the coupling between the interfaces in three-layer viscous fluid systems solving successively the eigenvalue and initial value problems. In Section 2, the problem is exposed with all the considered assumptions. Then, the dispersion relation is determined. In Section 3, the eigenvalue problem is solved and compared with the experimental data of Adkins et al. [17]. A parametric study is also carried out to study the effects of the physical parameters on the coupling between the interfaces. In Section 4, the initial value problem is solved, which allows us to compare the temporal evolution of the physical quantities with direct numerical simulations.

## 2. Linear stability analysis

### 2.1. Formulation of the problem

We consider three superimposed fluids in Cartesian coordinates $(e_x, e_y, e_z)^T$ under isothermal conditions and under the gravity field $g = -ge_z$. The subscript $i \in \{1; 2; 3\}$ refers respectively to the bottom layer (1), middle layer (2) and top layer (3), while the subscript $ii+1 \in \{12; 23\}$ ($i \in \{1; 2\}$) refers respectively to the interfaces (1)–(2) and (2)–(3) between layers (1) and (2) and layers (2) and (3); the subscripts $x$, $y$ and $z$ refer respectively to the component in the direction $e_x$, $e_y$ and $e_z$. The fluids are considered immiscible, incompressible with densities: $\rho_1$, $\rho_2$ and $\rho_3$ as well as Newtonian with dynamic viscosities: $\mu_1$, $\mu_2$ and $\mu_3$ and kinematic viscosities: $\nu_1$, $\nu_2$ and $\nu_3$. The layers have infinite horizontal lengths. The extreme layers have infinite vertical depths while the middle one has a finite thickness $h$. For the present theory the surface tension at the interfaces $\gamma_{12}$ and $\gamma_{23}$ are considered (Fig. 2).

The physical quantities such as velocities $u_1$, $u_2$ and $u_3$, interface positions: $z_{12}$ and $z_{23}$ and pressures $p_1$, $p_2$ and $p_3$ can thus be determined from the conservation of mass and momentum equations:

$$ \nabla \cdot \mathbf{u}_i = 0 \quad \forall i \in \llbracket 1; 3 \rrbracket, \qquad (5) $$

$$ \rho_i (\mathbf{u}_{i,t} + (\mathbf{u}_i \cdot \nabla) \mathbf{u}_i) = -\nabla p_i + \nabla \cdot \bar{\mathbf{\tau}}_i + \rho_i g \quad \forall i \in \llbracket 1; 3 \rrbracket, \qquad (6) $$

<table><caption>**Table 2**
Fluid properties for the two cases studied by Adkins et al. [17]. Cases 1 and 2 refer to the three-layer systems (1)/(2)/(1) and (1)/(2)/(3) respectively.</caption><thead><tr><th>(i)</th><th>ρ<sub>i</sub><br>[kg m<sup>-3</sup>]</th><th>μ<sub>i</sub><br>[kg m<sup>-1</sup> s<sup>-1</sup>]</th><th>γ<sub>ii+1</sub><br>[kg s<sup>-2</sup>]</th></tr></thead><tbody><tr><td>(3)</td><td>1,3</td><td>1,8 × 10<sup>-5</sup></td><td></td></tr><tr><td>(2)</td><td>1398</td><td>6,2 × 10<sup>-3</sup></td><td>72 × 10<sup>-3</sup></td></tr><tr><td>(1)</td><td>773</td><td>3,3 × 10<sup>-3</sup></td><td>3,5 × 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Schema of the studied system.</figcaption></figure>

with $\bar{\mathbf{\tau}}_i$ being the viscous stress tensor in layer (i):

$$ \bar{\mathbf{\tau}}_i = \mu_i (\nabla \mathbf{u}_i + (\nabla \mathbf{u}_i)^T), \qquad (7) $$

and $u_{i,l}$ being the partial derivative of the velocity in the medium (i) with respect to time. Subsequently, the comma will be used to denote the partial derivatives.

To solve these volume equations it is necessary to know the behaviour of the physical quantities that compose them across the interfaces. For incompressible fluids without mass transfer, the normal component of the velocity should be continuous at the interfaces according to the mass balance equation. For viscous fluids, the tangential component of the velocity is also continuous at the interfaces if we consider no slip between adjacent fluids. Unlike its tangential component, the normal component of the viscous stress tensor should be bounded at the interfaces according to the momentum balance equation. These jump conditions can respectively be written as follows:

$$ \mathbf{u}_i \cdot \mathbf{n}_{ii+1} = \mathbf{u}_{i+1} \cdot \mathbf{n}_{ii+1} = z_{ii+1, i} \quad \text{at } z = z_{ii+1} \quad \forall i \in \llbracket 1; 2 \rrbracket, \qquad (8) $$

$$ \mathbf{u}_i \times \mathbf{n}_{ii+1} = \mathbf{u}_{i+1} \times \mathbf{n}_{ii+1} \quad \text{at } z = z_{ii+1} \quad i \in \llbracket 1; 2 \rrbracket, \qquad (9) $$

$$ \left( p_{i+1} - p_i - \gamma_{ii+1} \left( \frac{z_{ii+1,x,x}}{(1+z_{ii+1,x}^2)^{\frac{3}{2}}} + \frac{z_{ii+1,y,y}}{(1+z_{ii+1,y}^2)^{\frac{3}{2}}} \right) \right) \mathbf{n}_{ii+1} \\ = (\bar{\tau}_{i+1} - \bar{\bar{\tau}}_i) \cdot \mathbf{n}_{ii+1} \quad \text{at } z = z_{ii+1} \quad \forall i \in \llbracket 1; 2 \rrbracket, \qquad (10) $$
with $n_{il+1}$ the vector normal to interface (i)–(i + 1), directed from layer (i) to layer (i + 1):

$$ \mathbf{n}_{il+1} = \frac{1}{\sqrt{z_{il+1,x}^2 + z_{il+1,y}^2 + 1}} \begin{pmatrix} -z_{il+1,x} \\ -z_{il+1,y} \\ 1 \end{pmatrix}. \quad (11) $$

In the basic state, the fluids are considered at rest and the two interfaces between them are flat. The index * refers to this state. Thus, the velocities in each layer and the position of the interfaces can be written as follows:

$$ \mathbf{u}_i^* = \mathbf{0}, \quad (12) $$

$$ z_{il+1}^* = (i - 1)h. \quad (13) $$

Once the basic state is fully characterized, a small wave shape disturbance with a wave number: $k = (k_x, k_y, 0)^T$ of norm $k = (k_x^2 + k_y^2)^{1/2}$ and a complex number time coefficient $\omega$ is introduced around it. The index ' refers to the disturbed quantities. A physical quantity $Q$ can thus be written as follows:

$$ Q(x, y, z, t) = Q^*(x, y, z) + Q'(x, y, z, t). \quad (14) $$

Analysing the disturbance into normal modes, the dependence on $x, y$ and $z$ for space and $t$ for time of the disturbed part of a physical quantity $Q$ is seeked as follows:

$$ Q'(x, y, z, t) = f(z)e^{(jk_x x + jk_y y + \omega t)}, \quad (15) $$

with j beeing the imaginary unit and $f(z)$ a function of the variable $z$ to be determined by the volume equations.

## 2.2. Dimensional analysis and linearization

The twelve physical parameters of the problem are as follows: $\rho_1$, $\rho_2$, $\rho_3$, $v_1$, $v_2$, $v_3$, $\gamma_{12}$, $\gamma_{23}$, $g$, $h$, $k$ and $\omega$. The three fundamental units present are thus: $kg$, $m$ and $s$. According to the Vaschy-Buckingham theorem [24], we can define $12-3 = 9$ dimensionless numbers to characterize this problem. We propose the following set of dimensionless parameters:

$$ \begin{aligned} R_i &= \rho_i \rho_2^{-1} = \rho_i m_{ref}^{-1} l_{ref}^3 & \forall i \in [1; 3] & (R_2 = 1), \\ N_i &= v_i k^{3/2} g^{-1/2} = v_i l_{ref}^{-2} l_{ref} & \forall i \in [1; 3], \\ \Gamma_{il+1} &= \gamma_{il+1} k^2 \rho_2^{-1} g^{-1} = \gamma_{il+1} m_{ref}^{-1} l_{ref}^2 & \forall i \in [1; 2], \\ H &= kh = h l_{ref}^{-1}, \\ \Omega &= \omega g^{-1/2} k^{-1/2} = \omega t_{ref}, \end{aligned} \quad (16) $$

where we introduce a characteristic length $l_{ref} = k^{-1}$, a characteristic time $t_{ref} = (gk)^{-1/2}$ and a characteristic mass $m_{ref} = \rho_2 k^{-3}$. Subsequently, all the physical quantities $u_i$, $p_i$, $z_{il+1}$, the space unities $x$, $y$, $z$, the time unity $t$ and the operator $\nabla$ are replaced by $l_{ref} l_{ref}^{-1} \tilde{u}_i$, $m_{ref} l_{ref} l_{ref}^{-2} \tilde{p}_i$, $l_{ref} z_{il+1}$, $l_{ref} x$, $l_{ref} y$, $l_{ref} z$, $t_{ref} t$ and $l_{ref} l_{ref}^{-1} \tilde{\nabla}$ and are now dimensionless.

In the basic state, the dimensionless assumptions can be written as follows:

$$ \mathbf{u}_i^* = \mathbf{0}, \quad (17) $$

$$ z_{il+1}^* = (i - 1)H. \quad (18) $$

The pressure in layer (i) is defined up to a constant $p_{ci}$ by the momentum balance equation:

$$ p_i^* = -R_i z + p_{ci}, \quad (19) $$

with

$$ p_{ci} = p_i^*(z = (i - 1)H) + R_i(i - 1)H. \quad (20) $$

The jump condition for the momentum gives us the relation between these constants:

$$ p_{ci+1} = p_{ci} + (R_i - R_{i+1})(i - 1)H. \quad (21) $$

After disturbing the system in the way that we previously introduced, we obtain the dimensionless linearized volume equations below:

$$ j k_x k^{-1} u_{ix} + j k_y k^{-1} u_{iy} = -u_{iz,z}, \quad (22) $$

$$ j k_x k^{-1} p'_i = -R_i \Omega u_{ix} + R_i N_i (u_{ix,z,z} - u_{ix}), \quad (23) $$

$$ j k_y k^{-1} p'_i = -R_i \Omega u_{iy} + R_i N_i (u_{iy,z,z} - u_{iy}), \quad (24) $$

and

$$ p'_{iz} = -R_i \Omega u_{iz} + R_i N_i (u_{iz,z,z} - u_{iz}), \quad (25) $$

as $u'_i = u_i$. By summing the product of Eq. (23) by $k_x k^{-1}$ with the product of Eq. (24) by $k_y k^{-1}$ and using Eq. (22), we obtain a new expression for the pressure disturbance:

$$ p'_{i} = (-R_{i} \Omega u_{iz} + R_{i} N_{i} (u_{iz,z,z} - u_{iz})), z. \quad (26) $$

By comparing the derivative with respect to $z$ of this expression with Eq. (25), we obtain a differential equation on the vertical component of the velocity alone:

$$ (\Omega u_{iz} - N_i(u_{iz,z,z} - u_{iz}))_{,z,z} - (\Omega u_{iz} - N_i(u_{iz,z,z} - u_{iz})) = 0. \quad (27) $$

General solutions to this differential equation are linear combinations of $e^{\pm z}$ and $e^{\pm q_i z}$ with $q_i = (1 + \Omega N_i^{-1})^{1/2}$. Let $K$ be a $3 \times 4$ matrix defined such that the columns of the matrix $e^{zK}$ are the different functions of $z$ involved in $u_{iz}: e^{-z}$, $e^{-q_i z}$, $e^z$ and $e^{q_i z}$; its rows 1, 2 and 3 refer to the layers (1), (2) and (3):

$$ K = \begin{pmatrix} -1 & -q_1 & 1 & q_1 \\ -1 & -q_2 & 1 & q_2 \\ -1 & -q_3 & 1 & q_3 \end{pmatrix}. \quad (28) $$

The shapes of the vertical velocities in each layer can thus be defined as follows:

$$ u_{1z} = \sum_{n=3}^{4} A_n e^{K_{1n} z} e^{(jk_x x + jk_y y + \omega t)}, \quad (29) $$

$$ u_{2z} = \sum_{n=1}^{4} B_n e^{K_{2n} z} e^{(jk_x x + jk_y y + \omega t)}, \quad (30) $$

$$ u_{3z} = \sum_{n=1}^{2} C_n e^{K_{3n} z} e^{(jk_x x + jk_y y + \omega t)}. \quad (31) $$

As the velocity cannot be infinite, the coefficients $A_1, A_2, C_3$ and $C_4$ must be zero. The dimensionless linearized jump conditions can be written as follows:

$$ M \cdot (A_3 A_4 B_1 B_2 B_3 B_4 C_1 C_2)^T = \mathbf{0}, \quad (32) $$

where $M$ is an 8 x 8 matrix defined as follows:

$$ M = \begin{pmatrix} A_{12}^{+} & B_{12}^{-} & B_{12}^{+} & 0_{4,2} \\ 0_{4,2} & B_{23}^{-} & B_{23}^{+} & C_{23}^{-} \end{pmatrix}. \quad (33) $$

$A_{12}^{+}, B_{12}^{-}, B_{12}^{+}, B_{23}^{-}, B_{23}^{+}, C_{23}^{-}$ and $0_{4,2}$ are 4 x 2 matrices. The letters $A$, $B$ and $C$ refer to the coefficients of the matrix $M$ corresponding to the velocity in layers (1), (2) and (3), respectively. The subscripts 12 and 23 refer to the coefficients of the matrix $M$ corresponding to the interfaces (1)–(2) and (2)–(3), respectively. The exponents + and – refer to the coefficients of the matrix $M$ corresponding to the exponential functions of the positive values of $z$: $e^{+z}$ and $e^{-q_i z}$ and the negative values of $z$: $e^{-z}$ and $e^{+q_i z}$. $A_{12}^{+}, B_{12}^{-}, B_{12}^{+}, B_{23}^{-}, B_{23}^{+}$ and $C_{23}^{-}$ are expressed in Appendix A while $0_{4,2}$ is the zero matrix.

## 2.3. Dispersion relation

The dispersion relation is defined by the following equation:

$$ |M| = 0. \quad (34) $$
Solving this equation allows us to determine the eigenvalues of the problem posed, that is, the values of the complex time coefficient ($\Omega$ in its dimensionless form). Taking into account the viscosity only yields numerical solutions for the eigenvalues as highlighted in the introduction. We nevertheless found an analytical eigenvalue:

$$\Omega = -N_2. \quad (35)$$

This negative real eigenvalue corresponds to a stable mode that only appears when a three-layer fluid system and viscosity are considered.

To find all the solutions analytically, viscosity must be disregarded. Eqs. (36) and (41) are obtained from the matrix M by removing the columns corresponding to the terms in $e^{\pm q_1 z}$ as well as the rows corresponding to the continuity of the tangential velocity and the continuity of the tangential component of the stress tensor at each interface. Obviously, the viscosity is zero for the normal component of the viscous stress tensor at each interface. The matrix M is therefore a 4 × 4 matrix for the inviscid case instead of a 8 × 8 matrix for the viscous case. By also neglecting the surface tension (which does not change the analytical or numerical problem resolution), we find the dispersion relation of Mikaelian [14]:

$$-(1 - R_1)(1 - R_3) \sinh(H) + (\Omega)^2 (R_1 - R_1 R_3) [\cosh(H) + \sinh(H)] + (\Omega)^4 [(R_1 + R_3) \cosh(H) + (1 + R_1 R_3) \sinh(H)] = 0. \quad (36)$$

This is a second-degree polynomial with respect to $\Omega^2$, and its solutions can easily be determined analytically. For extreme values of H, Mikaelian [14] showed that these solutions tend towards those of two-layer fluid subsystems.

When $H \gg 1$, the hyperbolic functions $\sinh(H)$ and $\cosh(H)$ have $e^H/2$ as limit and the solutions are similar to those found when studying two subsystems independently: the two infinite layers (1) and (2) around the interface (1)–(2) and the two infinite layers (2) and (3) around the interface (2)–(3):

$$\Omega = \pm \left( (1 - R_1) (1 + R_1)^{-1} \right)^{1/2}, \quad (37)$$

and

$$\Omega = \pm \left( (R_3 - 1) (1 + R_3)^{-1} \right)^{1/2}. \quad (38)$$

The solutions represented by Eq. (37) are identical to those presented in the introduction (Eq. (1)) but formulated in a dimensionless manner.

When $H \ll 1$, the hyperbolic functions $\sinh(H)$ and $\cosh(H)$ have limits 0 and 1, respectively, and the solutions are similar to those found when studying the two infinite layers (1) and (3) around the interface (1)–(3) plus a zero eigenvalue:

$$\Omega = \pm \left( (R_3 - R_1) (R_1 + R_3)^{-1} \right)^{1/2}, \quad (39)$$

and

$$\Omega = 0. \quad (40)$$

Neglecting viscosity but no longer surface tension, the dispersion relation becomes:

$$(R_1 - 1 + \Gamma_{12})(1 - R_3 + \Gamma_{23}) \sinh(H) + \Omega^2 [(R_1 - R_3) \times (\cosh(H) + \sinh(H)) - \Gamma_{12}(\cosh(H) + R_3 \sinh(H)) - \Gamma_{23}(\cosh(H) + R_1 \sinh(H))] + \Omega^4 [(R_1 + R_3) \cosh(H) + (1 + R_1 R_3) \sinh(H)] = 0. \quad (41)$$

When $H \gg 1$, the eigenvalues of the previous dispersion relation become those of two subsystems with two separate layers. One corresponding to the interface (1)–(2) surrounded by fluids (1) and (2) of infinite depths. The other corresponding to the interface (2)–(3) surrounded by fluids (2) and (3) of infinite depths:

$$\Omega = \pm \left( (1 - R_1 - \Gamma_{12}) (1 + R_1)^{-1} \right)^{1/2}, \quad (42)$$

and

$$\Omega = \pm \left( (R_3 - 1 - \Gamma_{23}) (1 + R_3)^{-1} \right)^{1/2}. \quad (43)$$

These solutions are identical to those presented in the introduction (Eq. (2)) but formulated in a dimensionless manner. A stabilization term appears and the system is unstable if $1 - R_1 > \Gamma_{12}$ or $R_3 - 1 > \Gamma_{23}$.

When $H \ll 1$, the eigenvalues of the previous dispersion relation become those of a system with two superimposed layers (1) and (3) for which the surface tension between the layer (1) and (3) would be equal to the sum between the surface tensions between the layer (1) and (2) and (2) and (3) plus a zero eigenvalue:

$$\Omega = \pm \left( (R_3 - R_1 - (\Gamma_{12} + \Gamma_{23})) (R_1 + R_3)^{-1} \right)^{1/2}, \quad (44)$$

and

$$\Omega = 0. \quad (45)$$

A stabilization term also appears in this limit and here the system is unstable if $R_3 - R_1 > \Gamma_{12} + \Gamma_{23}$.

When only viscosity is not considered, when H tends to 0, the dispersion relation is thus similar to that of the two-layer case with an equivalent surface tension equal to the sum of the surface tension at each interface. The system then tends, as in the Mikaelian [14] case, towards a two-layer fluid system only if we consider the surface tension between the extreme layers as an equivalent surface tension depending on the properties of the middle layer.

Let us now return to our viscous problem whose solutions are not analytical (except $\Omega = -N_2$). The dispersion relation that accounts for both viscosity and surface tension will be very long to derive since M is a 8 × 8 matrix. Eventually, the very large dispersion relation will be more complicated to use than the matrix in a software for determining the eigenvalues. Interestingly, viscosity does not change the behaviour of the solutions for large values of H. To show this, $|M|$ can be written as the determinant of a block matrix as follows:

$$|M| = | - (B_{12}^{+} | 0_{4,2} )(B_{23}^{+} | C_{23}^{-})^{-1} (0_{4,2} | B_{23}^{-} )(B_{23}^{+} | C_{23}^{-}) + (A_{12}^{+} | B_{12}^{-})(B_{23}^{+} | C_{23}^{-}) |. \quad (46)$$

where $(M_1|M_2)$ denotes the concatenation between matrices $M_1$ and $M_2$. When the value of H is large enough, the concatenated matrix $(0_{4,2}|B_{23}^{-})$ tends to the zero matrix by 4 and $|M|$ can be simplified as follows:

$$|M| = |(A_{12}^{+} | B_{12}^{-})|(B_{23}^{+} | C_{23}^{-})|. \quad (47)$$

Eq. (46) is determined using the Gaussian elimination. It does not help to provide the dispersion relation for the general case but it allows to simplify the case when H becomes very large and thus to determine Eq. (47). The eigenvalues are then the solutions of these two equations:

$$|(B_{23}^{+} | C_{23}^{-})| = 0, \quad (48)$$

and

$$|(A_{12}^{+} | B_{12}^{-})| = 0. \quad (49)$$

$|(B_{23}^{+} | C_{23}^{-})| = 0$ is the same dispersion relation as in Chandrasekhar's book [6, p. 443], which takes into account the two infinite layers (1) and (2) around the interface (1)–(2) positioned at $z = 0$. $|(A_{12}^{+} | B_{12}^{-})| = 0$ is also the same relation but taking into account the two infinite layers (2) and (3) around the interface (2)–(3) positioned at $z = H$. As determined by Mikaelian [14] for the inviscid case without taking into account the surface tension, the two interfaces can be studied independently. Adding surface tension does not change this behaviour.

The assumption of ignoring the viscosity has the advantage of obtaining an analytical expression for the eigenvalues as a function of the other physical parameters. Within the scope of this paper, we impose a value for each physical parameter to determine numerically all the eigenvalues except for the analytical one for a given case. To obtain the dispersion relation, we vary H. A parametric study of this dispersion relation can therefore highlight the effect of a parameter on the eigenvalues and thus the temporal evolution of any single-mode disturbance. Although any hypothetical case could be analysed through this study, Adkins et al. [17] provided experimental data regarding the eigenvalues for two cases to which this theory applies.
<table><caption>**Table 3**
Variable parameters *k* and *h* corresponding to the two cases studied by Adkins et al. [17]. *ω*<sub>Ac</sub>, *ω*<sub>Mc</sub> and *ω*<sub>c</sub> refer respectively to the growth rate extracted by Adkins et al. [17], to the positive real eigenvalue predicted by Mikaelian theory and to the one predicted by this theory. *c* ∈ {1; 2} refers to case *c*.</caption><thead><tr><th>k<br>[m<sup>-1</sup>]</th><th>h<br>[mm]</th><th>ω<sub>A1</sub><br>[s<sup>-1</sup>]</th><th>ω<sub>A2</sub><br>[s<sup>-1</sup>]</th><th>ω<sub>M1</sub><br>[s<sup>-1</sup>]</th><th>ω<sub>M2</sub><br>[s<sup>-1</sup>]</th><th>ω<sub>1</sub><br>[s<sup>-1</sup>]</th><th>ω<sub>2</sub><br>[s<sup>-1</sup>]</th><th>|ω<sub>A1</sub> - ω<sub>1</sub>|<br>ω<sub>A1</sub></th><th>|ω<sub>A2</sub> - ω<sub>2</sub>|<br>ω<sub>A2</sub></th></tr></thead><tbody><tr><td rowspan="3">503</td><td>1</td><td>-</td><td>20,8</td><td>-</td><td>31,7</td><td>-</td><td>23,6</td><td>-</td><td>0,13</td></tr><tr><td>2</td><td>24,8</td><td>24,0</td><td>36,4</td><td>35,8</td><td>29,7</td><td>29,0</td><td>0,20</td><td>0,21</td></tr><tr><td>5</td><td>25,6</td><td>25,2</td><td>37,6</td><td>37,6</td><td>31,5</td><td>31,4</td><td>0,23</td><td>0,25</td></tr><tr><td rowspan="3">419</td><td>1</td><td>-</td><td>18,0</td><td>-</td><td>27,7</td><td>-</td><td>21,3</td><td>-</td><td>0,18</td></tr><tr><td>2</td><td>22,1</td><td>21,3</td><td>32,8</td><td>31,9</td><td>27,7</td><td>26,8</td><td>0,25</td><td>0,26</td></tr><tr><td>5</td><td>22,9</td><td>22,4</td><td>34,3</td><td>34,2</td><td>29,9</td><td>29,8</td><td>0,31</td><td>0,33</td></tr><tr><td rowspan="5">314</td><td>1</td><td>-</td><td>15,0</td><td>-</td><td>22,1</td><td>-</td><td>17,7</td><td>-</td><td>0,18</td></tr><tr><td>2</td><td>19,7</td><td>18,9</td><td>27,6</td><td>26,3</td><td>24,2</td><td>22,9</td><td>0,23</td><td>0,21</td></tr><tr><td>5</td><td>21,5</td><td>21,2</td><td>29,5</td><td>29,3</td><td>26,8</td><td>26,6</td><td>0,25</td><td>0,25</td></tr><tr><td>2</td><td>18,3</td><td>17,1</td><td>24,0</td><td>22,4</td><td>21,4</td><td>19,9</td><td>0,17</td><td>0,16</td></tr><tr><td>5</td><td>20,6</td><td>20,3</td><td>26,1</td><td>25,8</td><td>24,2</td><td>23,9</td><td>0,17</td><td>0,18</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Case 1 ≡ (1)/(2)/(1) (left) and case 2 ≡ (1)/(2)/(3) (right) studied by Adkins et al. [17].</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Representation of the decomposition of the three-layer fluid system that models case 1 of Adkins et al. [17] into two-layer fluid subsystems for extreme values of *H* according to the theory of Mikaelian [14].</figcaption></figure>

# 3. Resolution of the eigenvalue problem

## 3.1. *Studied case*

In this section, the experiments performed by Adkins et al. [17] are used to support our theoretical results. These experiments were performed with the magnetic levitation method to control the initial conditions by imposing a single-mode perturbation. The authors studied two cases where the bottom layer is oil, while the middle layer is water with paramagnetic salt in order to use magnetic levitation. For case 1, the upper fluid is the same as the bottom one (Fig. 3 on the left), while it is air for case 2 (Fig. 3 on the right). In both cases, the bottom interface is gravitationally unstable while the top one is gravitationally stable. As long as the wavenumber value is less than the cutoff wavenumber, solving the eigenvalue problem allows to determine two modes: an unstable one associated with a positive real eigenvalue and a stable one associated with two complex conjugate eigenvalues with negative real part.

The properties of the fluids are presented in Table 2. The variable parameters are the thickness of the intermediate layer and the wavenumber of the perturbation (the values are given in the first two columns of the Table 3). For a given pair of these parameters, the authors extracted the growth rate of the gravitationally unstable interface by assuming that the interfaces are decoupled: the gravitationally unstable interface is influenced only by the real eigenvalues, while the gravitationally stable interface is influenced only by the complex eigenvalues. The positive real eigenvalue is determined numerically by a zero-finding algorithm. The relative error between this theoretical prediction and the experimental value is between 13% and 33% depending on the cases (Table 3), the experimental values being

systematically lower than the theoretical values. Depending on the values of *k* and *h*, it is between 10% and more than 20% lower than that obtained by comparing to the theory of Mikaelian [14]. This relative error remains significant. This can probably be explained by the edge effects present in the experiments and also by the fact that the coupling of the interfaces was not taken into account while the gravitationally stable interface tends to become flat, and the closer the interfaces are, the more the amplitude of the deformation of the gravitationally unstable interface will decrease and the more the positive real eigenvalue will decrease. Before solving the initial value problem that allows to obtain the time evolution of the two interfaces as a function of all the eigenvalues, it is appropriate to analyse the influence of the physical parameters on the eigenvalues. Only case 1, i.e. the three-layer fluid system (1)/(2)/(1) (Fig. 3 on the left), is studied for simplification reasons. The number of parameters is reduced, six instead of nine, as $R_3 = R_1$, $N_3 = N_1$ and $\Gamma_{23} = \Gamma_{12}$. Note that this case does not depend on the direction of gravity. In the following analysis, the analytical eigenvalue ($\Omega = -N_2$) is neglected compared to the other eigenvalues, which are determined numerically.

## 3.2. *Comparison to the solutions of the two two-layer subsystems*

To better understand the evolution of the coupling between interfaces with respect to the thickness of the dimensionless medium layer, let us first decompose the system into two subsystems associated with each interface surrounded by fluids of infinite depth (Fig. 4). The first subsystem corresponds to an unstable interface surrounded by an infinite layer of water (above) and an infinite layer of oil (below). The second subsystem is the opposite, corresponding to a stable interface surrounded by an infinite layer of water (below) and an infinite layer of oil (above).
<figure><img src="image_8.png" /><figcaption>Fig. 5. Real parts (top) and imaginary parts (bottom) of the eigenvalues $\Omega$ plotted separately as a function of the thickness of the medium layer $H$. The dashed lines, round and diamond marks correspond, respectively, to the eigenvalues of the two subsystems, to the unstable and stable eigenvalues of the three-layer fluid system for $(R_1, N_1, N_2, \Gamma_{12}) = (5.53 \times 10^{-1}, 1.52 \times 10^{-2}, 1.61 \times 10^{-2}, 6.45 \times 10^{-2})$.</figcaption></figure>

of oil (above). For the first case, the spectrum of $M$ contains only one positive real eigenvalue, while for the second, there are two complex conjugate eigenvalues with a negative real part. The three-layer fluid system has one positive real eigenvalue and two complex conjugate eigenvalues with a negative real part:

$$S\rho(M) = \{\Omega_1, \Omega_2, \overline{\Omega_2}\}. \qquad (50)$$

This spectrum is solved as a function of $H$ for a set of fixed parameters $(R_1, N_1, N_2, \Gamma_{12}) = (5.53 \times 10^{-1}, 1.52 \times 10^{-2}, 1.61 \times 10^{-2}, 6.45 \times 10^{-2})$ obtained by considering the dimensional physical parameters of the Table 2 associated with the first value of $k$ of the Table 3. The solutions are plotted in Fig. 5. As determined previously, when the value of $H$ becomes very large, the eigenvalues of the system tend towards those of the two subsystems. These values are constant, because they do not depend on $H$. Remarkably, when $H$ tends to zero, the eigenvalues of the system under study do not tend to the values of a single-layer fluid system, with the upper layer superimposing the same fluid as in the lower one, as predicted by the theory of Mikaelian [14].

### 3.3. Parametric study

For this parametric study, the reference set of fixed parameters $(R_1, N_1, N_2, \Gamma_{12})$ corresponds to that of Fig. 5. Let us then vary the parameters $R_1, N_1$ and $N_2$ individually and analyse the evolution of the three eigenvalues $\Omega$ separated into real and imaginary parts with respect to the coupling parameter $H$. Note that the parameter $\Gamma_{12} = \gamma_{12}k^2\rho_2^{-1}g^{-1}$ is not analysed here but its value remains lower than $1-R_1$ in the parametric study to ensure that the gravitationally unstable interface is not stabilized by the surface tension effect. The asymptotic values of the previously defined two-layer fluid subsystems are added to the figures. This allows an additional point of analysis on the coupling between the interfaces, since all the eigenvalues of the three-layer fluid system tend to these asymptotes when $H$ becomes sufficiently large.

Concerning the parameter $R_1 = \rho_1\rho_2^{-1}$ (Fig. 6), the positive values of the real and imaginary asymptotes decrease when this parameter increases. It is consistent since when $R_1$ tends to unity, the system becomes stable. Second, the real parts of the eigenvalues corresponding to the three-layer fluid system tend more slowly towards their asymptotes when $R_1$ increases. This behaviour is not similar to the imaginary parts of $\Omega$, which tend more slowly to their asymptotes for median values of $R_1$. When $H=0$, the imaginary parts of $\Omega$ also depend on $R_1$: their absolute value decrease when $R_1$ increases.

Concerning the parameter $N_1 = \nu_1 k^{3/2} g^{-1/2}$ (Fig. 7), the values of the asymptotes vary in a similar way to the previous case. Generally, the behaviour differs from the previous case when $N_1 \gg 1$ where a cusp appears in the evolution of the negative real part of $\Omega$ corresponding

to the value of $H$ for which the imaginary parts of $\Omega$ differ from zero. From a certain value (not determined here), $N_1$ tends to make the oscillation of the relatively coupled interfaces disappear.

Concerning the parameter $N_2 = \nu_2 k^{3/2} g^{-1/2}$ (Fig. 8), as with the previous case, a cusp appears in the evolution of the negative real part of $\Omega$ when $N_2 \gg 1$. It corresponds to the value of $H$ for which all the imaginary parts of $\Omega$ tend towards zero. Unlike the previous case, from a certain value (not determined here), $N_2$ tends to make the oscillation of relatively decoupled interfaces disappear. Note that the cusps might be the signature of another branch of eigenvalue that would need to be explored further.

In the next part of this paper, the effects of all the modes on the temporal evolution of the interfaces will be studied using the resolution of the initial value problem.

## 4. Resolution of the initial value problem

### 4.1. Method

#### 4.1.1. Theoretical physical quantities

All the physical quantities $Q_i$ of the layer (i) can be written as follows:

$$Q_i = \Re \left( \sum_{n=1}^{N} Q_i^n(z) \exp \left( jk^{-1}(k_x x + k_y y) + \Omega^n t \right) \right). \qquad (51)$$

where $N$ is the number of modes and $Q_i^n$ the amplitude of $Q_i$ corresponding to the eigenvalue $\Omega^n$:

$$Q_i^n(z) = Q_{i1}^n e^{-z} + Q_{i2}^n e^{-q_i^2 z} + Q_{i3}^n e^{z} + Q_{i4}^n e^{q_i^2 z}. \qquad (52)$$

For example, $Q_i^n = A^n$, $Q_i^n = B^n$ and $Q_i^n = C^n$ for the velocity in the bottom, middle and top layers, respectively. For the following, we restrict ourselves to the two-dimensional case without loss of generality: $k_x = k$ and $k_y = 0$. It is interesting to consider the velocity in the middle layer, as this layer is in contact with the two interfaces. The initial conditions imposed at the interface could thus be written as a function of the amplitudes present in this velocity:

$$u_{2z} = \Re \left( \sum_{n=1}^{N} \sum_{m=1}^{4} B_m^n \exp \left( K_{2m}^n z \right) \exp \left( \Omega^n t \right) \exp \left( jx \right) \right). \qquad (53)$$

Let us now introduce the disturbance of the displacement ($\xi$ such as $\xi, t \equiv u$) expressed in the middle layer for the vertical component:

$$\xi_{2z} = \Re \left( \sum_{n=1}^{N} \sum_{m=1}^{4} B_m^n \exp \left( K_{2m}^n z \right) (\Omega^n)^{-1} \exp \left( \Omega^n t \right) \exp \left( jx \right) \right). \qquad (54)$$

At $z=0$, this quantity is equal to the disturbance of the bottom interface position. At $z=H$, it is equal to the disturbance of the top interface position:

$$\eta_{ii+1} = \Re \left( \sum_{n=1}^{N} H_{ii+1}^n \exp \left( \Omega^n t \right) \exp \left( jx \right) \right), \qquad (55)$$

with:

$$H_{ii+1}^n = \sum_{m=1}^{4} B_m^n \exp \left( K_{2m}^n (i-1)H \right) (\Omega^n)^{-1}. \qquad (56)$$

For the resolution of the eigenvalue problem, the eight equations (jump conditions) did not allow us to determine the eight velocity constants but rather two eigenvalues (plus their conjugates). The total number of velocity constants is thus sixteen. It is then necessary to impose two additional initial conditions ($N=2$) to be able to determine all the constants:

$$\eta_{ii+1}(t=0) = a_i \cos(x). \qquad (57)$$

where $i=1$ for the condition at the bottom interface and $i=2$ for the condition at the top interface. This means that:

$$\Re \left( \sum_{n=1}^{2} H_{ii+1}^{n} \right) = a_i. \qquad (58)$$
<figure><img src="image_9.png" /><figcaption>Fig. 6. Effect of $R_1$ for $N_1 = 1.52 \times 10^{-2}$, $N_2 = 1.61 \times 10^{-2}$ and $\Gamma_{12} = 6.45 \times 10^{-2}$ on the real parts (top) and imaginary parts (bottom) of the eigenvalues $\Omega$ plotted separately as a function of the thickness of the medium layer $H$. The dashed lines, round and diamond marks correspond to the eigenvalues of the two subsystems, the unstable and stable eigenvalues of the three-layer fluid system, respectively.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Effect of $N_1$ for $R_1 = 5.53 \times 10^{-1}$, $N_2 = 1.61 \times 10^{-2}$ and $\Gamma_{12} = 6.45 \times 10^{-2}$ on the real parts (top) and imaginary parts (bottom) of the eigenvalues $\Omega$ plotted separately as a function of the thickness of the medium layer $H$. The dashed lines, round and diamond marks correspond to the eigenvalues of the two subsystems, the unstable and stable eigenvalues of the three-layer fluid system, respectively.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Effect of $N_2$ for $R_1 = 5.53 \times 10^{-1}$, $N_2 = 1.52 \times 10^{-2}$ and $\Gamma_{12} = 6.45 \times 10^{-2}$ on the real parts (top) and imaginary parts (bottom) of the eigenvalues $\Omega$ plotted separately as a function of the thickness of the medium layer $H$. The dashed lines, round and diamond marks correspond to the eigenvalues of the two subsystems, the unstable and stable eigenvalues of the three-layer fluid system, respectively.</figcaption></figure>

Solving this new system allows us to determine all the velocity constants. The horizontal component of the velocity can be determined with the equation of mass conservation:

$$u_{ix} = u_{iz,z} \qquad (59)$$

For the two subsystems, the four jump conditions allowed us to find one eigenvalue (plus its conjugate). Thus only one initial condition ($N = 1$) is necessary to be able to determine all four constants. We consider $i = 1$ for the unstable subsystem and $i = 2$ for the stable one. This means that:
<figure><img src="image_12.png" /><figcaption>Fig. 9. Initial conditions on the interface position and velocity for the theory and simulations with the system corresponding to case 1 studied by Adkins et al. [17] for $H = 2.51 \times 10^{-1}$ (left), $H = 5.03 \times 10^{-1}$ (middle) and $H = 2.51$ (right). The coloured circles refer to Fig. 10.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Interface positions of interest at a given time. The colours magenta, red, blue and green refer to the lower interface at $x = 0$ and $x = \pi$ and the upper interface at $x = 0$ and $x = \pi$, respectively.</figcaption></figure>

$$ \Re(H_{it+1}^n) = a_i. \qquad (60) $$

For the unstable subsystem, the velocities in each layer can thus be written:

$$ u_{2z} = \Re \left( \sum_{m=1}^{2} B_m^n \exp(K_{2m}^n z) \exp(\Omega^n t) \right) \cos(x), \qquad (61) $$

and

$$ u_{1z} = \Re \left( \sum_{m=3}^{4} A_m^n \exp(K_{1m}^n z) \exp(\Omega^n t) \right) \cos(x). \qquad (62) $$

For the amplitudes of these velocities, we found the values:

$$ (A_3, A_4, B_1, B_2) = a_1 (3.16 \times 10^{-4}, -6.64 \times 10^{-5}, 2.95 \times 10^{-4}, -4.52 \times 10^{-5}). \qquad (63) $$

The temporal evolution of the interface position can be written as follows:

$$ \eta_{12}(x, t) = a_1 e^{\Omega^1 t} \cos(x). \qquad (64) $$

For the stable subsystem the velocities in each layer can be written as follows:

$$ u_{2z} = \Re \left( \sum_{m=3}^{4} B_m^n \exp(K_{2m}^n z) \exp(\Omega^n t) \right) \cos(x), \qquad (65) $$

and

$$ u_{3z} = \Re \left( \sum_{m=1}^{2} C_m^n \exp(K_{3m}^n z) \exp(\Omega^n t) \right) \cos(x). \qquad (66) $$

For the amplitudes of these velocities we found the values:

$$ (B_3, B_4, C_1, C_2) = a_2 (8.96 \times 10^{-1} + 4.44 \times 10^1 j, -6.06 + 6.42 j, -7.38 \times 10^{-1} + 4.16 \times 10^1 j, -4.43 + 3.62 j). \qquad (67) $$

The temporal evolution of the interface position can be written as follows:

$$ \eta_{23}(x, t) = \cos(kx) e^{\Re(\Omega^2)t} (\Re(\Omega^2)^2 + \Im(\Omega^2)^2)^{-1} [\cos(\Im(\Omega^2)t)a_2 + \sin(\Im(\Omega^2)t)((\Re(B_3) + \Re(B_4))\Im(\Omega^2) - (\Im(B_3) + \Im(B_4))\Re(\Omega^2))] $$

<caption>Table 4</caption>

$t_{10}$ as a function of the initial amplitude.

<table><tr><td>2πa<sub>1</sub><sup>-1</sup></td><td>50</td><td>100</td><td>200</td></tr><tr><td>t<sub>10</sub></td><td>4.6</td><td>6.3</td><td>7.6</td></tr></table>

$$ \times \Re(\Omega^2)] . \qquad (68) $$

The initial velocity in each layer with the initial position of the interfaces between them are the initial conditions of the direct numerical simulations presented in the next subsection (see Fig. 9).

### 4.1.2. Direct numerical simulations

To cross-validate the developed theory, direct numerical simulations are performed using the in-house code ARCHER [25,26]. This code has been used extensively in the past to study a wide range of cases such as atomization [25], microfluidics [27] and droplet fragmentation [28]. It has been also used as a reference for cross-validation for canonical cases such as Rayleigh–Taylor instability [29], Plateau–Rayleigh instability [30] and Lamb oscillations [31].

The code is based on a Cartesian staggered mesh. The continuity and momentum equations (Eqs. (5) and (6)) are coupled using a projection algorithm with an explicit Euler temporal scheme. The coupled level-set and volume-of-fluid method is applied to allow both second-order accuracy of the interface description and mass conservation. The viscous term is discretized using the Sussman scheme [32], while the convective flux terms are computed to ensure the consistency between the mass and the momentum equation [33]. The surface tension is taken into account by introducing the continuum source force into the Navier–Stokes equations [34].

The domain length in the horizontal direction is set to $2\pi/k$. The theoretical solution is developed for the infinite depth in the top and bottom layer. Nevertheless, this is not allowed for the direct numerical simulations. A sensibility study regarding the edge effects was performed, showing that the use of a domain height four times larger than the wavelength is sufficient, which is consistent with the deep-water approximation. In a similar way, a mesh convergence study was conducted, and the solutions for the two-layer case are converged for 64 mesh nodes in the wavelength. Nevertheless, in the case of three layers, the mesh is increased to ensure that the mesh size is large enough to describe the middle layer (at least 8 mesh nodes within the middle layer thickness). The top and bottom boundary conditions are set to symmetric conditions, whereas the horizontal periodic boundary conditions are set in the horizontal direction.

To initialize the level-set the exact distance to the interface is computed for every cell centre. The volume-of-fluid is then derived from the level-set. Since the velocity field given by Eqs. (29)–(31) are developed for small perturbations, the continuity at the interface is set at $z = 0$ and $z = h$. Thus, these equations cannot be directly used to initialize the numerical velocity. The exact position of the interface should be considered. Mapping is made for $z$ using the level-set field. The mapped values of $z$ ensure that the transitions between the velocities for each phase occur at $0$ and $h$.

For all the studied configurations, the extreme positions of each interface are extracted for comparison with the theory (Fig. 10).
<figure><img src="image_14.png" /><figcaption>Fig. 11. Comparison between the theoretical (dashed lines) and numerical (full lines) temporal evolution of the interface position at $x = 0$ (magenta) and $x = \pi$ (red) for the unstable subsystem and different values of $a_1$. A visualization of the flow evolution of the unstable two-layer subsystem for $2\pi a_1^{-1} = 100$ is given in the Supplementary material (Video 1).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Comparison between the theoretical (dashed lines) and numerical (full lines) temporal evolution of the interface position at $x = 0$ (blue) and $x = \pi$ (green) for the stable subsystem and different values of $a_2$. A visualization of the flow evolution of the stable two-layer subsystem for $2\pi a_2^{-1} = 10$ is given in the Supplementary material (Video 2).</figcaption></figure>

## 4.2. Two subsystems

### 4.2.1. Unstable case

Let us take the parameter values corresponding to the bottom two layers of the reference case: ($R_1, N_1, N_2, \Gamma_{12}$) = (5.53 × 10⁻¹, 1.52 × 10⁻², 1.61 × 10⁻², 6.45 × 10⁻²). They only highlight the Rayleigh-Taylor instability. Three values of the initial amplitude are studied: $2\pi a_1^{-1} \in \{50; 100; 200\}$. The previous result concerning this subsystem is compared in Fig. 11 with the numerical simulations for each value of the initial amplitude.

We firstly observe that up to a certain time, the theory coincides with the simulations, although the theory increasingly overestimates the amplitude of the interface position. This discrepancy is not surprising since when the amplitude of the disturbance becomes too large, the linear theory is no longer valid. Moreover, for the present linear theory the top and bottom of the interface are always symmetrical, they tend to move towards the bottom in the simulations. Let us introduce the quantity:

$$ \eta_{12}(x = 0, t) - \eta_{12}(x = \pi, t). \qquad (69) $$

This quantity is theoretically equal to zero. It allows us to determine the time from which non-linear effects are no longer negligible. Let us arbitrarily choose the time $t_{10}$ from which this quantity is equal to 10%. We observe as expected that $t_{10}$ increases when the initial amplitude decreases (Table 4). However, regardless of the value of the initial amplitude, $t_{10}$ corresponds to the time at which the amplitude reaches a value of unity.

### 4.2.2. Stable case

Let us take the parameter values corresponding to the bottom two layers of the reference case: ($R_3, N_3, N_2, \Gamma_{23}$) = (5.53 × 10⁻¹, 1.52 × 10⁻², 1.61 × 10⁻², 6.45 × 10⁻²). They allow us to only highlight a gravity wave. Three values of the initial amplitude are studied: $2\pi a_2^{-1} \in \{5, 10, 50\}$. The previous result concerning this subsystem is compared in Fig. 12 with numerical simulations for each value of the initial amplitude.

Firstly, we do not observe the same behaviour as with the previous case. Here, using a sufficiently small initial amplitude, the simulations support the theory. However, when the initial amplitude is too large,

the linear stability analysis cannot be used, because a shift and an error at the extreme amplitudes are observed. These errors increase with the initial amplitude. The cases where the initial amplitude is too important ($2\pi a_1^{-1} < 10$) will no longer be studied, since the linear theory is never valid.

## 4.3. Three-layer system

Let us take the parameter values corresponding to the reference case: ($R_1 = R_3, N_1 = N_3, N_2, \Gamma_{12} = \Gamma_{23}$) = (5.53 × 10⁻¹, 1.52 × 10⁻², 1.61 × 10⁻², 6.45 × 10⁻²). The initial amplitudes are set to the same value $a_1 = a_2 = a$. The three values of the initial amplitudes as with the unstable subsystem are studied: $2\pi a_1^{-1} = 2\pi a_2^{-1} = 2\pi a^{-1} \in \{50, 100, 200\}$. Five values of the middle layer thickness (i.e., the coupling parameter) are also studied: $H \in \{2.51 \times 10^{-1}, 3.77 \times 10^{-1}, 5.03 \times 10^{-1}, 1.01, 2.51\}$.

For all three values of the initial amplitude $a$ and five values of the middle layer thickness $H$ (total of 15 configurations), the amplitudes of the three velocities through $X^n$ are determined using the previously defined method. For each configuration, numerical simulations are performed for cross-validation purposes and to show the non-linear departure.

As expected with the results concerning the unstable subsystem, when the initial amplitude of the interface positions decreases the time during which the theory and the numerical simulations coincide increases (Fig. 16). As a general comment on these figures, we identify three main regimes concerning the interface behaviour depending on the thickness of the middle layer $H$ for our choice of parameters.

The first regime (see the first line of Fig. 16) occurs when $H$ is high enough ($H = 2.51$). In this regime, the behaviour of the two interfaces is independent: the amplitude of the lower interface increases exponentially, while the amplitude of the upper one oscillates as expected for decoupled interfaces. This result indicates that the interfaces can be considered decoupled for values $H$ lower than the deep water assumption one ($2\pi = 6.28$). This first regime is called the “decoupled regime”. For $H = 1.01$ we observe a departure from the “decoupled regime” (see the second line of Fig. 16), where the amplitude of the upper interface begins by oscillating, although it is quickly caught by the bottom interface whose amplitude continues to
<figure><img src="image_16.png" /><figcaption>Fig. 13. Temporal evolution of the position of the interfaces of the numerical simulations of the three-layer system for $H = 2.51$ and $2\pi a^{-1} = 100$. See video 3 in the additional material.</figcaption></figure>

increase exponentially. In the simulations, this is observed when the linear stability analysis is no longer valid.

The second regime (see the third line of Fig. 16) occurs when $H$ continues to decrease ($H = 5.03 \times 10^{-1}$). In this regime, the amplitude of the upper interface does not have time to oscillate. The mode associated with the positive real eigenvalue is the dominant mode, with the amplitude of the two interfaces increasing exponentially. This second regime is called the “unstable mode-dominant regime”. When $H$ continues to decrease, we observe a departure from the “unstable mode-dominant regime” (see the fourth line of Fig. 16) for $H = 3.77 \times 10^{-1}$, where the amplitude of the two interfaces begins by decreasing (or even oscillating) before increasing as for the previous regime.

The third regime (see the fifth line of Fig. 16) occurs when $H$ continues to decrease ($H = 2.51 \times 10^{-1}$). In this regime, the amplitude of the two interfaces clearly oscillates before the lower interface catches up with the upper one. The smaller $H$ is, the longer the amplitude of the interfaces oscillates. The modes associated with the imaginary eigenvalues are the dominant modes. This third regime is called the “stable mode-dominant regime”.

The temporal evolution of the interface positions for the cases $H = 2.51$, $H = 5.03 \times 10^{-1}$ and $H = 2.51 \times 10^{-1}$ is respectively shown in Figs. 13–15. These figures allow us to observe the three different regimes. In Fig. 13 we observe the decoupled regime where the upper interface

seems to behave like a gravity wave while the lower interface seems to behave like a Rayleigh–Taylor instability even adopting its well-known “mushroom shape” due to the growth of nonlinearities. In Fig. 14 we observe the unstable mode-dominant regime where the amplitude of the two interfaces increases exponentially without oscillating initially. In Fig. 15 we observe the stable mode-dominant regime where the amplitude of the two interfaces oscillates without increasing initially.

As for the unstable subsystem, we introduce the two quantities:

$$ \eta_{ii+1}(x = 0, t) - \eta_{ii+1}(x = \pi, t). \qquad (70) $$

The values of $t_{10}$ for which these quantities are equal to 10% are plotted for the two interfaces and the previous set of parameters $H$ and $a$ (Fig. 17). For both interfaces, the values of $t_{10}$ tend towards infinity as $H$ tends towards zero. This is not surprising since the three-layer fluid system tends towards a stable one-layer system. When $H$ increases, the values of $t_{10}$ for the bottom interface tend towards the values of the unstable subsystem. The values of $t_{10}$ for the top interface tend towards infinity (value of the stable subsystem). For sufficiently small initial amplitudes, $t_{10}$ strictly decreases as a function of $H$ for the bottom interface. For the top interface the minimum of $t_{10}$ is reached for $H$ of the order of unity.
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Temporal evolution of the position of the interfaces of the numerical simulations of the three-layer system for $H = 5.03 \times 10^{-1}$ and $2\pi a^{-1} = 100$. See video 5 in the additional material.</figcaption></figure>

**5. Conclusion**

To summarize, linear stability analysis was performed to study the effect of viscosity and surface tension in a system composed of three superimposed immiscible incompressible Newtonian fluids under the gravity field. The system of interest is composed of a gravitationally unstable interface at the bottom and a gravitationally stable interface at the top to highlight respectively the Rayleigh–Taylor instability and a gravity wave as well as the coupling between these two phenomena. The physical parameters correspond to a symmetrical system composed of a layer of water surrounded by two layers of oil studied experimentally in a previous work. The parameter that characterizes how these two flows are coupled together in the system is the product between the wave number of the perturbation and the middle layer thickness $H = kh$. When $H$ is large enough, this system tends towards two different two-layer fluid subsystems, which were also defined to maintain this asymptotic behaviour in all the studies. For the linear stability analysis, the study of the three-layer fluid system and the two two-layer fluid subsystems independently are almost the same when $H > 2\pi$. The eigenvalue problem was solved for this set of parameters to determine the dimensionless time coefficient of the perturbation $\Omega = \omega g^{-1/2}k^{-1/2}$ as a function of $H$. As $H$ increases the positive real eigenvalue increases and tends towards the asymptotic value of the growth rate of the

unstable two-layer fluid subsystem. As $H$ decreases the positive real eigenvalue tends to zero since the system tends to a single-layer fluid system. This system is not exactly the one formed by gathering the bottom and top layers since an equivalent surface tension depending on the properties of the middle layer causes the interface between the two identical layers to oscillate. The amplitude of this oscillation decreases with time. A parametric study was performed to show how this dispersion relation is influenced by the density ratio $R_1 = \rho_1\rho_2^{-1}$ and the dimensionless viscosity of each layer $N_i = \nu_i g^{-1/2}k^{3/2}$; $R_1$ tends to promote the coupling between the two interfaces, $N_1$ tends to make the oscillation of relatively coupled interfaces disappear, whereas $N_2$ tends to make the oscillation of relatively decoupled interfaces disappear. The dimensionless surface tension, $\Gamma_{12} = \gamma_{12}\rho_2^{-1}g^{-1}k^2$ was fixed to ensure that the bottom interface is not stabilized by surface tension ($\Gamma_{12} < 1 - R_1$). The temporal evolution of the physical quantities is obtained by solving the initial value problem. For this purpose, sinusoidal interfaces are initially imposed as initial conditions. Numerical simulations were made using the in-house code ARCHER initialized by both the initial position of the interfaces and the initial velocities in the entire domain. This allowed us to compare the temporal evolution of the interface positions, which supports these theoretical results until the amplitude of one interface becomes too large in the same way as the non linear effects of the flow. We identified three main regimes for the behaviour
<figure><img src="image_18.png" /><figcaption>Fig. 15. Temporal evolution of the position of the interfaces of the numerical simulations of the three-layer system for $H = 2.51 \times 10^{-1}$ and $2\pi a^{-1} = 100$. See video 7 in the additional material.</figcaption></figure>

of the interfaces when $H$ varies: two coupled regimes, the stable mode-dominant and the unstable mode-dominant regimes and a decoupled regime. When $H$ is close to zero, the amplitude of the complex eigenvalue mode of the two interfaces is much larger than that of the positive real eigenvalue mode. Then, when $H$ increases, the positive real eigenvalue mode is dominant for both interfaces. Finally, when $H$ continues to increase, the interfaces decouple. The amplitude of the positive real eigenvalue mode is dominant for the bottom interface (gravitationally unstable), while that of the complex eigenvalue mode is dominant for the top interface (gravitationally stable).

The perspectives of this work are to study the effects of including an additional volume force such as the Lorentz force for conductive fluids subjected to an electric current and then to explore the edge effects by considering both horizontal and vertical walls.

### CRediT authorship contribution statement

**A. Simon:** Writing – original draft, Validation, Software, Methodology, Formal analysis, Data curation, Conceptualization. **J.C. Brän-dle de Motta:** Writing – review & editing, Visualization, Validation, Supervision, Software. **C. Dumouchel:** Writing – review & editing, Supervision. **M.-C. Renoult:** Writing – review & editing, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

### Acknowledgements

The authors acknowledge the financial support provided by the French Agence Nationale de la Recherche (ANR) LaBEx EMC3 through the HILIMBA project (Grant No. 10-LABX-0009).

This work was granted access to the HPC resources of IDRIS, TGCC and CINES under the allocation A0172B10101 made by GENCI (Grand Equipement National de Calcul Intensif) and CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie) under the scientific project N. 201 704.

### Appendix A

In Section 2.2 the matrix $M$ was defined as a block matrix as follows:

$$M = \begin{pmatrix} A_{12}^{+} & B_{12}^{-} & B_{12}^{+} & 0_{4,2} \\ 0_{4,2} & B_{23}^{-} & B_{23}^{+} & C_{23}^{-} \end{pmatrix}.$$

Here are the expressions of these blocks:

$$A_{12}^{+} = \begin{pmatrix} 1 & 1 \\ 1 & q_1 \\ \frac{2R_1 N_1}{R_1 + \frac{R_1-1+\Gamma_{12}}{2\Omega^2}} - \frac{N_2-R_1 N_1}{\Omega} & \frac{(1+q_1^2)R_1 N_1}{R_1-1+\Gamma_{12}} - q_1 \frac{N_2-R_1 N_1}{\Omega} \end{pmatrix},$$
<figure><img src="image_19.png" /><figcaption>Fig. 16. Temporal evolution of the interface positions at $x = 0$ (red for $z_{12}$ and green for $z_{23}$) and $x = \pi$ (magenta for $z_{12}$ and blue for $z_{23}$) for different values of the initial amplitude for $H = 2.51$ (first line), $H = 1.01$ (second line), $H = 5.03 \times 10^{-1}$ (third line), $H = 3.77 \times 10^{-1}$ (fourth line) and $H = 2.51 \times 10^{-1}$ (fifth line) for numerical simulations (full lines) and theoretical predictions (dashed lines). The visualization of the middle column cases ($2\pi a^{-1} = 100$) are given in the Supplementary material (Videos 3 to 7, respectively). The visualization of the middle row cases ($H = 5.03 \times 10^{-1}$) are given in the Supplementary material (Videos 8, 5 and 9, respectively).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 17. $t_{10}$ as function of $H$ for each interface. The dashed, dotted, and dashdotted lines are the asymptotic values of the unstable subsystem for $2\pi a^{-1} = 200$, 100 and 50, respectively.</figcaption></figure>
$$B_{12}^{-} = \begin{pmatrix} -1 & -1 \\ 1 & q_2 \\ -2N_2 & -(1+q_2^2)N_2 \\ 1+\frac{R_1-1+\Gamma_{12}}{2\Omega^2}+\frac{N_2-R_1N_1}{\Omega} & \frac{R_1-1+\Gamma_{12}}{2\Omega^2}+q_2\frac{N_2-R_1N_1}{\Omega} \end{pmatrix},$$

$$B_{12}^{+} = \begin{pmatrix} -1 & -1 \\ -1 & -q_2 \\ -2N_2 & -(1+q_2^2)N_2 \\ -1+\frac{R_1-1+\Gamma_{12}}{2\Omega^2} - \frac{N_2-R_1N_1}{\Omega} & \frac{R_1-1+\Gamma_{12}}{2\Omega^2} - q_2\frac{N_2-R_1N_1}{\Omega} \end{pmatrix},$$

$$B_{23}^{-} = \begin{pmatrix} \exp(-H) & \exp(-q_2 H) \\ -\exp(-H) & -q_2 \exp(-q_2 H) \\ 2N_2 \exp(-H) & (1+q_2^2) N_2 \exp(-q_2 H) \\ \left(-1+\frac{1-R_1+\Gamma_{23}}{2\Omega^2}+\frac{R_3 N_3-N_2}{\Omega}\right) \exp(-H) & \left(\frac{1-R_1+\Gamma_{23}}{2\Omega^2}+q_2\frac{R_3 N_3-N_2}{\Omega}\right) \exp(-q_2 H) \end{pmatrix},$$

$$B_{23}^{+} = \begin{pmatrix} \exp(H) & \exp(q_2 H) \\ \exp(H) & q_2 \exp(q_2 H) \\ 2N_2 \exp(H) & (1+q_2^2) N_2 \exp(q_2 H) \\ \left(1+\frac{1-R_1+\Gamma_{23}}{2\Omega^2}-\frac{R_3 N_3-N_2}{\Omega}\right) \exp(H) & \left(\frac{1-R_1+\Gamma_{23}}{2\Omega^2}-q_2\frac{R_3 N_3-N_2}{\Omega}\right) \exp(q_2 H) \end{pmatrix},$$

$$C_{23}^{-} = \begin{pmatrix} -\exp(-H) & -\exp(-q_3 H) \\ \exp(-H) & q_3 \exp(-q_3 H) \\ -2R_3 N_3 \exp(-H) & -(1+q_3^2) R_3 N_3 \exp(-q_3 H) \\ \left(R_3+\frac{1-R_1+\Gamma_{23}}{2\Omega^2}+\frac{R_3 N_3-N_2}{\Omega}\right) \exp(-H) & \left(\frac{1-R_1+\Gamma_{23}}{2\Omega^2}+q_3\frac{R_3 N_3-N_2}{\Omega}\right) \exp(-q_3 H) \end{pmatrix}.$$

## Appendix B. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.physd.2025.134670.

### Data availability

Data will be made available on request.

### References

[1] Lord Rayleigh, Investigation on the character of the equilibrium of an incompressible heavy fluid of variable density, Proc. Lond. Math. Soc. s1-14 (1883) 170-177.

[2] Lord Rayleigh, On the stability, or instability, of certain fluid motion, Proc. Lond. Math. Soc. s1-11 (1880) 57-72.

[3] Sir J.J. Thomson, Hydrokinetic solutions and observations, Lond. Edinb. Dublin Philos. Mag. J. Sci. s4-42 (1871) 362-377.

[4] J.C. Maxwell, Capillary Action. Ency. Brit., ninth ed., 1876, pp. 56-71, 1875-89 5.

[5] W.J. Harrison, The influence of viscosity on the oscillations of superposed fluids, Proc. Lond. Math. Soc. s2-6 (1908) 396-405.

[6] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Instability, Oxford University Press, 1961.

[7] Sir H. Lamb, Hydrodynamics, Cambridge University Press, 1895.

[8] Sir G. Taylor, The instability of liquid surfaces when accelerated in a direction perpendicular to their planes, Proc. R. Soc. A: Math. Phys. Eng. Sci. 201 (1950) 192-196.

[9] Y. Zhou, Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence and mixing. I, Phys. Rep. 720-722 (2017) 1-136.

[10] Y. Zhou, R.J.R. Williams, P. Ramaprabhu, M. Groom, B. Thornber, A. Hillier, W. Mostert, B. Rollin, S. Balachandar, P.D. Powell, A. Mahalov, N. Attal, Rayleigh-Taylor and Richtmyer-Meshkov instabilities: A journey through scales, Phys. D: Nonlinear Phenom. 423 (2021) 132838.

[11] Y. Zhou, J.D. Sadler, O.A. Hurricane, Instabilities and mixing in inertial confinement fusion, Annu. Rev. Fluid. Mech. 57 (2025) 197-225.

[12] A.D. Sneyd, Stability of fluid layers carrying a normal electric current, J. Fluid Mech. 156 (1985) 223-236.

[13] D.H. Kelley, T. Weier, Fluid mechanics of liquid metal batteries, Appl. Mech. Rev. 70 (2018) 1-23.

[14] K.O. Mikaelian, Rayleigh-Taylor instabilities in stratified fluids, Phys. Rev. A 26 (1982) 2140-2158.

[15] S. Parhi, G. Nath, A sufficient criterion for Rayleigh-Taylor instability of incompressible viscous three-layer flow, g, Int. J. Eng. Sci. 29 (1991) 1439-1450.

[16] J.W. Jacobs, S.B. Dalziel, Rayleigh-Taylor instability in complex stratifications, J. Fluid Mech. 542 (2005) 251-279.

[17] R. Adkins, E.M. Shelton, M.C. Renoult, P. Carles, C. Rosenblatt, Interface coupling and growth rate measurement in multilayer Rayleigh-Taylor instabilities, Phys. Rev. Fluids 2 (2017) 1-8.

[18] L. Cottier, G. Brenn, M.-C. Renoult, Initial conditions to study the temporal behaviour of a viscoelastic liquid jet under perturbation, At. Sprays 31 (2021) 73-87.

[19] A. Prosperetti, Viscous effects on small-amplitude surface waves, Phys. Fluids 19 (1975) 195-203.

[20] A. Prosperetti, Motion of two superposed viscous fluids, Phys. Fluids 24 (1981) 1217-1223.

[21] P.K. Farsoiya, Y.S. Mayya, R. Dasgupta, Axisymmetric viscous interfacial oscillations - theory and simulations, J. Fluid Mech. 826 (2017) 797-818.

[22] M. Vartdal, A.N. Osnes, Linear motion of multiple superposed viscous fluids, Phys. Rev. E 99 (2019) 1-9.

[23] K.O. Mikaelian, Rayleigh-Taylor instabilities in finite-thickness fluids with viscosity and surface tension, Phys. Rev. E 54 (1996) 3676-3680.

[24] E. Buckingham, On physically similar systems; illustrations of the use of dimensional equations, Phys. Rev. 4 (1914) 345-375.

[25] T. Ménard, S. Tanguy, A. Berlemont, Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet, Int. J. Multiph. Flow 33 (5) (2007) 510-524.

[26] G. Vaudor, T. Ménard, W. Aniszewski, M. Doring, A. Berlemont, A consistent mass and momentum flux computation method for two phase flows. Application to atomization process, Comput. & Fluids 152 (2017) 204-216.

[27] J.B. Charpentier, J.C. Brändle de Motta, T. Ménard, Capillary phenomena in assemblies of parallel cylindrical fibers: From statics to dynamics, Int. J. Multiph. Flow 129 (2020) 103304.

[28] C. Deberne, V. Chéron, A. Poux, J.C. Brändle de Motta, Breakup prediction of oscillating droplets under turbulent flow, Int. J. Multiph. Flow 173 (2024) 104731.

[29] T. Chen, V. Chéron, G. Zhaoli, J.C. Brändle de Motta, T. Ménard, Lian-Ping Wang, Simulation of immiscible two-phase flows based on a kinetic diffuse interface approach, in: 10th International Conference on Multiphase Flow, 2019, p. 10.

[30] C. Dumouchel, W. Aniszewski, T.-T. Vu, T. Ménard, Multi-scale analysis of simulated capillary instability, Int. J. Multiph. Flow 92 (2017) 181-192.

[31] I. Roa, M.-C. Renoult, C. Dumouchel, J.C. Brändle de Motta, Droplet oscillations in a turbulent flow, Front. Phys. 11 (2023) 1173521.

[32] M. Sussman, K.M. Smith, M.Y. Hussaini, M. Ohta, R. Zhi-Wei, A sharp interface method for incompressible two-phase flows, J. Comput. Phys. 221 (2) (2007) 469-505.

[33] M. Rudman, A volume-tracking method for incompressible multifluid flows with large density variations, Internat. J. Numer. Methods Fluids 28 (2) (1998) 357-378.

[34] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (2) (1992) 335-354.