1. Introduction
The occurrence of convection heat transfer in a vertical porous layer saturated by a fluid is a topic that has been widely studied over the last few decades due to its potential interest for several applications. If the main practical use of this knowledge is for the thermal insulation of buildings and for devices such as breathing walls, the topic of convection in vertical porous layers may be important also for geophysics and for the design of filtration systems. The pioneering paper on this topic was that by Gill (Reference Gill1969). This author provided a straightforward and rigorous proof that a vertical porous layer with a homogeneous structure and subjected to a side heating via isothermal boundaries cannot display a multicellular pattern of convection heat transfer however large the Rayleigh number. This cornerstone result was obtained by modelling the flow through Darcy's law, the Boussinesq approximation and by assuming impermeable isothermal boundaries kept at different temperatures. The analysis carried out by Gill (Reference Gill1969) is grounded on the linear dynamics of perturbations superposed to the steady vertical buoyant flow in a conduction regime induced by a horizontal temperature gradient in the basic state.
The paper by Gill (Reference Gill1969) is the starting point of several later studies including the nonlinear extension of the stability proof (Straughan Reference Straughan1988; Flavin & Rionero Reference Flavin and Rionero1999), as well as variants involving the Prandtl–Darcy flow model (Rees Reference Rees1988) and local thermal non-equilibrium within the saturated porous material (Rees Reference Rees2011; Scott & Straughan Reference Scott and Straughan2013).
The use of momentum balance models for the seepage flow in the porous layer that include the velocity Laplacian term shows the possibility of a thermal instability and the emergence of convection patterns in the porous layer (Chen Reference Chen2004). This feature is an expected consequence of a momentum balance model similar to that for a fluid clear of porous material. In fact, linear convective instability emerges when a vertical and infinitely tall slot of fluid is bounded by impermeable plane walls at different temperatures (Vest & Arpaci Reference Vest and Arpaci1969). Kwok & Chen (Reference Kwok and Chen1987) present experimental evidence that convective instability in a vertical porous layer may arise when a sufficiently large temperature difference is forced across the layer boundaries. Those authors also propose possible explanations of the observed phenomenon based either on Brinkman's momentum transfer model or on the inclusion of variable viscosity effects (Kwok & Chen Reference Kwok and Chen1987).
The papers by Barletta (Reference Barletta2015) and by Shankar & Shivakumara (Reference Shankar and Shivakumara2022) offer apparently different approaches showing up the emergence of convective instability in a vertical porous slab subjected to side heating. The first one (Barletta Reference Barletta2015) focuses on the same flow set-up assumed by Gill (Reference Gill1969), but turning the impermeability boundary conditions into conditions of permeability by imposing a hydrostatic pressure distribution at the boundaries. The second one (Shankar & Shivakumara Reference Shankar and Shivakumara2022) considers a variant of Gill's set-up where the boundaries are impermeable, but the porous medium is heterogeneous with a transverse continuous change of the permeability. Both Barletta (Reference Barletta2015) and Shankar & Shivakumara (Reference Shankar and Shivakumara2022) proved that their variants of Gill's problem could lead to a condition of linear instability for sufficiently large Rayleigh numbers. Further recent explorations into different aspects of the onset of convection in porous vertical layers have been carried out by Barletta & Celli (Reference Barletta and Celli2021), Shankar, Shivakumara & Naveen (Reference Shankar, Shivakumara and Naveen2021) and Shankar, Naveen & Shivakumara (Reference Shankar, Naveen and Shivakumara2022).
The aim of this paper is to show that a common physical mechanism underlies the apparently different instabilities found by Barletta (Reference Barletta2015) and by Shankar & Shivakumara (Reference Shankar and Shivakumara2022). This task is achieved by envisaging a three-layer porous slab with impermeable isothermal boundaries kept at different temperatures. In this sandwiched porous slab, the two external layers are identical, while the core layer has different thermophysical properties. In particular, the external layers are considered as much more thermally conductive than the core. The basic buoyant flow in the internal core is identical to that devised by Gill (Reference Gill1969). Although the basic state is stationary with a purely vertical velocity, the core layer has interfaces to the external layers which may allow for a horizontal flow contribution when the basic state is perturbed. This circumstance is a relaxation of Gill's impermeability condition at the boundaries, if just the core layer is considered. Such a set-up is the link between the two studies by Barletta (Reference Barletta2015) and Shankar & Shivakumara (Reference Shankar and Shivakumara2022). In fact, the three-layer structure is effectively a horizontally heterogeneous medium with a piecewise-constant permeability. Furthermore, the internal core is a homogeneous porous layer with permeability conditions at the bounding interfaces to the external layers.
In this study, it is shown that the neutral stability condition is influenced by the permeability ratio between the external layers and the core layer and by the ratio between the thickness of the three-layer slab and that of the core layer. The role of such parameters in determining the stabilisation or destabilisation of the basic buoyant flow is identified by employing a numerical solution of the linear stability eigenvalue problem. This eigenvalue problem, formulated for the core layer, turns out to coincide with that solved by Barletta (Reference Barletta2015) in the asymptotic case where the external layers are much more permeable than the core. Hence, the neutral stability condition and the critical values of the wavenumber and of the Rayleigh number found by Barletta (Reference Barletta2015) are confirmed in the asymptotic case. Such an asymptotic condition turns out to be the most unstable relative to cases where the permeability ratio between the external layers and the core layer is finite. It is shown that the stability eigenvalue problem studied by Gill (Reference Gill1969) is recovered in the opposite asymptotic condition where the permeability ratio between the external layers and the core layer tends to zero. This is an expected result as, in this asymptotic case, the interfaces between the core layer and the external layers are effectively impermeable, thus reproducing just the same Gill's boundary conditions for the core layer.
2. Governing equations
Let us consider a vertical porous slab with a sandwiched structure. The 
$x$ axis is horizontal and perpendicular to the slab, the 
$y$ axis is also horizontal and the 
$z$ axis is vertical, so that the gravitational acceleration is 
$\boldsymbol {g} = - g \, \hat {\boldsymbol {e}}_z$, where 
$g$ is the modulus of 
$\boldsymbol {g}$ and 
$\hat {\boldsymbol {e}}_z$ is the unit vector of the 
$z$ axis.
As shown in figure 1, the core region 
$-L/2 \le x \le L/2$ is a layer of a porous material 
${\rm M}_1$, while the regions 
$-D/2 \le x < -L/2$ and 
$L/2 < x \le D/2$ are layers of a different porous material 
${\rm M}_2$. For the sake of simplicity, we assume that the slab extends without bounds over the 
$y$ and 
$z$ directions. The external boundaries 
$x=\pm D/2$ are impermeable and isothermal at temperatures 
$T_0 \pm \Delta T/2$, where 
$T_0$ and 
$\Delta T$ are constants, with 
$\Delta T > 0$.
Figure 1. Two-dimensional cross-section of the sandwiched porous slab in the 
$xz$ plane. The 
$y$ axis is perpendicular to the plane of the figure.
The flow model is based on the local mass, momentum and energy balance equations according to the Boussinesq approximation and to Darcy's law:
$$\begin{gather} \boldsymbol{\nabla} \cdot {\boldsymbol{u}_m} = 0, \end{gather}$$
$$\begin{gather}\displaystyle\frac{\mu}{K_m} \boldsymbol{u}_m ={-} \boldsymbol{\nabla}P_m + \rho_0 g \beta ( T_m - T_0 ) \hat{\boldsymbol{e}}_z , \end{gather}$$
$$\begin{gather}\sigma_m \displaystyle{{\partial T_m}\over{\partial t}} + \boldsymbol{u}_m \cdot \boldsymbol{\nabla}T_m = \alpha_m {\nabla ^2}{T_m},\quad m=1,2 , \end{gather}$$
where the index 
$m=1,2$ identifies the physical properties and the fields either in the core 
${\rm M}_1$ layer or in the external 
${\rm M}_2$ layers. Here, 
$t$ is time, 
$\boldsymbol {u}_m$ is the velocity field, 
$P_m$ is the dynamic pressure field and 
$T_m$ is the temperature field. We emphasise that 
$P_m$, the dynamic pressure, is the local difference between the pressure and the hydrostatic pressure, where the latter is evaluated with the reference fluid density 
$\rho _0$. The fluid saturating the three-layer slab is Newtonian with a dynamic viscosity 
$\mu$, a reference density 
$\rho _0$ and a coefficient of thermal expansion 
$\beta$.
The average properties of the saturated porous media are, with 
$m=1,2$: the permeability 
$K_m$, the thermal diffusivity 
$\alpha _m$ and the heat capacity ratio 
$\sigma _m$. The latter ratio is evaluated by dividing the volumetric heat capacity of each saturated porous medium by the volumetric heat capacity of the fluid. We also recall that the average thermal diffusivities, 
$\alpha _1$ and 
$\alpha _2$, are defined through the ratio of the average thermal conductivity of the fluid-saturated porous medium and the volumetric heat capacity of the fluid. The volumetric heat capacity of the fluid is the same for both media 
${\rm M}_1$ and 
${\rm M}_2$. Hence, the ratio 
$\alpha _2/\alpha _1$ is both a thermal diffusivity ratio and a thermal conductivity ratio.
2.1. Pressure–temperature formulation
Equation (2.1) can be expressed in a pressure–temperature formulation:
$$\begin{gather} {\nabla ^2}{P_m} = \rho_0 g \beta \displaystyle{{\partial T_m}\over{\partial z}}, \end{gather}$$
$$\begin{gather}\sigma_m \displaystyle{{\partial T_m}\over{\partial t}} - \frac{K_m}{\mu} \boldsymbol{\nabla}P_m \cdot \boldsymbol{\nabla}T_m + \frac{\rho_0 g \beta K_m}{\mu} ( T_m - T_0 ) {\displaystyle{{\partial T_m}\over{\partial z}}} = \alpha_m {\nabla ^2}{T_m},\quad m=1,2 . \end{gather}$$2.2. Boundary and interface conditions
The boundary and interface conditions are to be set at 
$x=\pm D/2$ and at 
$x=\pm L/2$, respectively:
$$\begin{gather} x ={\pm} D/2 : \quad {\displaystyle{{\partial P_2}\over{\partial x}}} = 0 , \quad T_2 = T_0 \pm \dfrac{\Delta T}{2}, \end{gather}$$
$$\begin{gather}x ={\pm} L/2 : \quad P_1 = P_2 , \quad K_1 {\displaystyle{{\partial P_1}\over{\partial x}}} = K_2 {\displaystyle{{\partial P_2}\over{\partial x}}}, \quad T_1 = T_2, \quad \alpha_1 {\displaystyle{{\partial T_1}\over{\partial x}}} = \alpha_2 {\displaystyle{{\partial T_2}\over{\partial x}}}. \end{gather}$$
In particular, the interface conditions reflect the momentum and energy balance across the planes at 
$x = \pm L/2$. They express the continuity of the pressure, of the normal component of the velocity, of the temperature and of the normal component of the heat flux density. The last interface condition (2.3) should employ the thermal conductivities instead of the thermal diffusivities, but we have already pointed out that the ratio 
$\alpha _2/\alpha _1$ coincides with the thermal conductivity ratio.
2.3. Dimensionless formulation
The mathematical model can be reformulated in dimensionless terms by scaling time, coordinates and fields as
$$\begin{gather} t^* = \dfrac{t}{\sigma_1 L^2/\alpha_1},\quad (x^*, y^*, z^*) = \dfrac{(x,y,z)}{L}, \quad \boldsymbol{u}^*_m = \dfrac{\boldsymbol{u}_m}{\alpha_1/L}, \end{gather}$$
$$\begin{gather}P^*_m = \dfrac{P_m}{\mu \alpha_1/K_1}, \quad T^*_m = \dfrac{T_m - T_0}{\Delta T}, \end{gather}$$with the dimensionless parameters
\begin{equation} R = \frac{\rho_0 g \beta \Delta T K_1 L}{\mu \alpha_1}, \quad a = \frac{D}{L}, \quad \xi =\frac{K_2}{K_1} , \quad \gamma = \frac{\alpha_2}{\alpha_1} \quad \tau = \frac{\sigma_2}{\sigma_1} . \end{equation}
Here, 
$a > 1$ and 
$R$ is the Darcy–Rayleigh number, hereafter called the Rayleigh number for conciseness. Equations (2.4) and (2.5a–e) allow one to rewrite (2.2) in a dimensionless form for the medium 
${\rm M}_1$:
$$\begin{gather} {\nabla ^2}{P_1} = R {\displaystyle{{\partial T_1}\over{\partial z}}}, \end{gather}$$
$$\begin{gather}{\displaystyle{{\partial T_1}\over{\partial t}}} - \boldsymbol{\nabla}P_1 \cdot \boldsymbol{\nabla}T_1 + RT_1 {\displaystyle{{\partial T_1}\over{\partial z}}} = {\nabla ^2}{T_1}, \end{gather}$$
and for the medium 
${\rm M}_2$:
$$\begin{gather} {\nabla ^2}{P_2} = R {\displaystyle{{\partial T_2}\over{\partial z}}}, \end{gather}$$
$$\begin{gather}\tau {\displaystyle{{\partial T_2}\over{\partial t}}} - \xi \boldsymbol{\nabla}P_2 \cdot \boldsymbol{\nabla}T_2 + R \xi T_2 {\displaystyle{{\partial T_2}\over{\partial z}}} = \gamma {\nabla ^2}{T_2}, \end{gather}$$while (2.3) reads
$$\begin{gather} x ={\pm} a/2 : \quad {\displaystyle{{\partial P_2}\over{\partial x}}} = 0 , \quad T_2 ={\pm}\dfrac{1}{2}, \end{gather}$$
$$\begin{gather}x ={\pm} 1/2 : \quad P_1 = P_2 , \quad {\displaystyle{{\partial P_1}\over{\partial x}}} = \xi {\displaystyle{{\partial P_2}\over{\partial x}}}, \quad T_1 = T_2, \quad {\displaystyle{{\partial T_1}\over{\partial x}}} = \gamma {\displaystyle{{\partial T_2}\over{\partial x}}}. \end{gather}$$In (2.6)–(2.8) and in the forthcoming analysis, the dimensionless fields, coordinates and time are denoted without the asterisks for simplicity of notation. This will not cause any ambiguity as we will only deal with dimensionless expressions, except when explicitly declared.
It is also worth saying that the second of equations (2.1) is rewritten in a dimensionless form as
$$\begin{gather} \boldsymbol{u}_1 ={-} \boldsymbol{\nabla}P_1 + R T_1 \hat{\boldsymbol{e}}_z , \end{gather}$$
$$\begin{gather}\boldsymbol{u}_2 ={-} \xi \boldsymbol{\nabla}P_2 + R \xi T_2 \hat{\boldsymbol{e}}_z , \end{gather}$$
for media 
${\rm M}_1$ and 
${\rm M}_2$, respectively.
3. The basic stationary flow
A basic stationary flow solution of (2.6)–(2.8) can be found such that
\begin{equation} \bar{P}_1 = 0 = \bar{P}_2, \quad {\displaystyle{{\partial {\bar T}_1}\over{\partial y}}} = 0 = {\displaystyle{{\partial {\bar T}_2}\over{\partial y}}}, \end{equation}
where the bar over the fields stands for basic state. Equation (3.1a,b) describes a two-dimensional, 
$y$-independent, flow regime where the pressure locally coincides with the hydrostatic pressure or, equivalently, the dynamic pressure is everywhere zero. As a consequence of (3.1a,b), (2.6)–(2.8) allow one to infer that, in the three-layer slab, the temperature profile 
$\bar {T}$ is a piecewise linear function of 
$x$ independent of 
$y$ and 
$z$. In fact, in the basic state, we have
$$\begin{gather} \bar{T}_1 = \dfrac{\gamma}{\gamma + a - 1} x, \end{gather}$$
$$\begin{gather}\bar{T}_2 =\begin{cases} \dfrac{1}{\gamma + a -1} x - \dfrac{\gamma - 1}{2 (\gamma + a - 1)}, \quad -\dfrac{a}{2} \le x <{-} \dfrac{1}{2},\\ \dfrac{1}{\gamma + a -1} x + \dfrac{\gamma - 1}{2 (\gamma + a - 1)}, \quad \dfrac{1}{2} < x \le \dfrac{a}{2}. \end{cases} \end{gather}$$Another consequence of (3.1a,b) is that (2.9) yields
\begin{equation} \bar{\boldsymbol{u}}_1 = (0, 0, R \bar{T}_1), \quad \bar{\boldsymbol{u}}_2 = (0, 0, R\xi \bar{T}_2), \end{equation}
which describes a purely vertical flow driven only by the buoyancy force. If 
$\bar {T}$ varies continuously along the range 
$-a/2 \le x \le a/2$, this is not the case for the vertical velocity component. Such a feature is a consequence of the different permeabilities of the porous media 
${\rm M}_1$ and 
${\rm M}_2$, so that continuity of the velocity profile occurs only for the special case 
$\xi =1$.
Figure 2 illustrates the basic solution by showing some plots of 
$\bar {T}$ versus 
$x$ over the whole range 
$-a/2 \le x \le a/2$. The geometrical ratio 
$a=2$ is prescribed, while different values of the conductivity ratio 
$\gamma > 1$ are considered. A focus on the case 
$\gamma > 1$ has been made since, as is discussed in § 3.1, our aim is to explore a condition of extremely high values of 
$\gamma$. It can be stressed that the case 
$\gamma = 100$ yields a temperature distribution which is almost uniform in the external layers (
$-a/2 \le x < -1/2$ and 
$1/2 < x \le a/2$). This behaviour is characteristic of the asymptotic case 
$\gamma \to \infty$.
Figure 2. Plots of 
$\bar {T}$ versus 
$x$ for 
$a=2$ and different values of 
$\gamma$.
3.1. Infinite thermal conductivity ratio
There are several practical cases where the limit 
$\gamma \to \infty$ is a fairly appropriate condition. Porous media with a very high thermal conductivity are the metal foams often employed in the design of heat exchangers, while the low-conductivity inner core can be devised as a slab of any thermal insulation material employed in the building industry. In the forthcoming stability analysis of the basic buoyant flow, we focus on this asymptotic condition. It is worth emphasising that the limit 
$\gamma \to \infty$ of (3.2) yields
\begin{equation} \bar{T}_1 = x, \quad \bar{T}_2 =\begin{cases} - \dfrac{1}{2}, \quad -\dfrac{a}{2} \le x <{-} \dfrac{1}{2},\\ \dfrac{1}{2}, \quad \dfrac{1}{2} < x \le \dfrac{a}{2}. \end{cases} \end{equation}
The basic state (3.4a,b) for the core layer 
${\rm M}_1$ is exactly that considered in the papers by Gill (Reference Gill1969) and by Barletta (Reference Barletta2015). Finally, we reckon that also the governing equations, the boundary and the interface conditions (2.7) and (2.8) undergo a marked simplification:
\begin{equation} {\nabla ^2}{P_2} = R {\displaystyle{{\partial T_2}\over{\partial z}}}, \quad {\nabla ^2}{T_2} = 0, \end{equation}with
$$\begin{gather} x ={\pm} a/2 : \quad {\displaystyle{{\partial P_2}\over{\partial x}}} = 0 , \quad T_2 ={\pm} \frac{1}{2}, \end{gather}$$
$$\begin{gather}x ={\pm} 1/2 :\quad P_1 = P_2 , \quad {\displaystyle{{\partial P_1}\over{\partial x}}} = \xi {\displaystyle{{\partial P_2}\over{\partial x}}}, \quad T_1 = T_2, \quad {\displaystyle{{\partial T_2}\over{\partial x}}} = 0 . \end{gather}$$4. Linearised perturbation dynamics
The onset of instability is studied by perturbing the basic state:
\begin{equation} P_m = \bar{P}_m + \epsilon \hat{P}_m , \quad T_m = \bar{T}_m + \epsilon \hat{T}_m , \quad m=1,2, \end{equation}
where 
$\epsilon >0$ is a perturbation parameter and the hat identifies the perturbation contributions to the pressure and temperature fields. The linear stability analysis is carried out by substituting (4.1a,b) into (2.6), (3.5a,b) and (3.6) and by neglecting the terms 
$\textit{O}\big(\epsilon ^2\big)$. We obtain, for the inner layer 
${\rm M}_1$,
\begin{equation} {\nabla ^2}{\hat{P}_1} = R {\displaystyle{{\partial {\hat{T}_1}}\over{\partial z}}},\quad {\displaystyle{{\partial {\hat{T}_1}}\over{\partial t}}} - {\displaystyle{{\partial {\hat{P}_1}}\over{\partial x}}} + Rx {\displaystyle{{\partial {\hat{T}_1}}\over{\partial z}}} = {\nabla ^2}{\hat{T}_1}, \end{equation}
and, for the external layers 
${\rm M}_2$,
\begin{equation} {\nabla ^2}{\hat{P}_2} = R {\displaystyle{{\partial {\hat{T}_2}}\over{\partial z}}},\quad {\nabla ^2}{\hat{T}_2} = 0, \end{equation}with
$$\begin{gather} x ={\pm} a/2 : \quad {\displaystyle{{\partial {\hat{P}_2}}\over{\partial x}}} = 0 , \quad \hat{T}_2 = 0, \end{gather}$$
$$\begin{gather}x ={\pm} 1/2 :\quad \hat{P}_1 = \hat{P}_2 , \quad {\displaystyle{{\partial {\hat{P}_1}}\over{\partial x}}} = \xi {\displaystyle{{\partial {\hat{P}_2}}\over{\partial x}}}, \quad \hat{T}_1 = \hat{T}_2, \quad {\displaystyle{{\partial {\hat{T}_2}}\over{\partial x}}} = 0 . \end{gather}$$Here, the features of the basic buoyant flow, defined by (3.1a,b) and (3.4a,b), have been employed.
4.1. Normal mode analysis
A normal mode expression of the perturbations 
$(\hat {P}_m, \hat {T}_m)$ is given by
\begin{equation} \left[\begin{array}{@{}c@{}}
\hat{P}_m \\ \hat{T}_m
\end{array}\right]=\left[\begin{array}{@{}c@{}} {f}_m(x) \\
{h}_m(x) \end{array}\right] {\rm e}^{{\rm i} (k_y y + k_z z
- \omega t)}, \quad \text{with}\ k_y \in \mathbb{R},k_z \in
\mathbb{R}, \omega \in \mathbb{C}, \ m = 1, 2,
\end{equation}
where 
$\boldsymbol {k} = (0, k_y, k_z)$ is the wavevector and 
$\omega$ is a complex angular frequency. The wavenumber is a positive quantity defined as the modulus of the wavevector, 
$k = |\boldsymbol {k}|$. The angular frequency is the real part of 
$\omega$, while the temporal growth rate is the imaginary part of 
$\omega$. By substituting (4.5) into (4.2a,b)–(4.4), one obtains for 
${\rm M}_1$
$$\begin{gather} f''_1 - k^2 f_1 - {\rm i} k_z R h_1 = 0 , \end{gather}$$
$$\begin{gather}h''_1 - (k^2 - {\rm i} \omega + {\rm i} k_z R x) h_1 + f'_1 = 0, \end{gather}$$
and for the external 
${\rm M}_2$ layers
$$\begin{gather} f''_2 - k^2 f_2 - {\rm i} k_z R h_2 = 0 , \end{gather}$$
$$\begin{gather}h''_2 - k^2 h_2 = 0, \end{gather}$$with the boundary and interface conditions
$$\begin{gather} x ={\pm} a/2 : \quad f'_2 = 0 , \quad h_2 = 0, \end{gather}$$
$$\begin{gather}x ={\pm} 1/2 : \quad f_1 = f_2 , \quad f'_1 = \xi f'_2, \quad h_1 = h_2, \quad h'_2 = 0 . \end{gather}$$
We can eliminate the dependence on the orientation of the wavevector 
$\boldsymbol {k}$ in (4.6)–(4.8) by defining a rescaled Rayleigh number, 
$S$, such that
\begin{equation} k S = k_z R. \end{equation}
$$\begin{gather} f''_1 - k^2 f_1 - {\rm i} k S h_1 = 0 , \end{gather}$$
$$\begin{gather}h''_1 - (k^2 - {\rm i} \omega + {\rm i} k S x) h_1 + f'_1 = 0 \end{gather}$$and
$$\begin{gather} f''_2 - k^2 f_2 - {\rm i} k S h_2 = 0 , \end{gather}$$
$$\begin{gather}h''_2 - k^2 h_2 = 0, \end{gather}$$respectively.
4.2. Eigenvalue problem for the core porous layer
Function 
$h_2$ is a solution of the second of equations (4.11) which, on account of (4.8), must satisfy the boundary conditions 
$h_2(\pm a/2)=0$. Then, it can be expressed as
\begin{equation} h_2 =\begin{cases} C_h \sinh \left(k \dfrac{a + 2 x}{2}\right), \quad -\dfrac{a}{2} \le x <{-} \dfrac{1}{2}, \\ \tilde{C}_h \sinh \left(k \dfrac{a - 2 x}{2}\right), \quad \dfrac{1}{2} < x \le \dfrac{a}{2}, \end{cases} \end{equation}
where 
$(C_h, \tilde {C}_h)$ are integration constants. However, (4.8) prescribes also the interface conditions 
$h'_2(\pm 1/2) = 0$. On account of (4.12), such conditions can be satisfied for 
$k>0$ only with 
$C_h = 0$ and 
$\tilde {C}_h = 0$. Hence, one has
\begin{equation} h_2 = 0 . \end{equation}
Function 
$f_2$ must be a solution of the first of equations (4.11) with 
$h_2 = 0$. Furthermore, on account of (4.8), 
$f_2$ satisfies the boundary conditions 
$f'_2(\pm a/2)=0$, so that
\begin{equation} f_2 =\begin{cases} C_f \cosh \left(k \dfrac{a + 2 x}{2}\right), \quad -\dfrac{a}{2} \le x <{-} \dfrac{1}{2}, \\ \tilde{C}_f \cosh \left(k \dfrac{a - 2 x}{2}\right), \quad \dfrac{1}{2} < x \le \dfrac{a}{2}, \end{cases} \end{equation}
where 
$(C_f, \tilde {C}_f)$ are integration constants.
Thus, from (4.14), the interface conditions 
$f_1(- 1/2) = f_2(- 1/2)$ and 
$f'_1(- 1/2) = \xi \, f'_2(- 1/2)$ given by (4.8) yield
\begin{equation} f_1 \left(-\dfrac{1}{2}\right) = C_f \cosh \left(k \dfrac{a - 1}{2}\right), \quad f'_1 \left(-\dfrac{1}{2}\right) = C_f \xi k \sinh \left(k \dfrac{a - 1}{2}\right). \end{equation}
By eliminating 
$C_f$, (4.15a,b) yields
\begin{equation} f'_1 \left(-\frac{1}{2}\right) - \xi k \tanh \left(k \dfrac{a - 1}{2}\right) f_1\left(-\frac{1}{2}\right) = 0 . \end{equation}
By employing the same method, the interface conditions 
$f_1(1/2) = f_2(1/2)$ and 
$f'_1(1/2) = \xi f'_2(1/2)$ given by (4.8) yield
\begin{equation} f_1 \left(\frac{1}{2}\right) = \tilde{C}_f \cosh \left(k \dfrac{a - 1}{2}\right), \quad f'_1 \left(\frac{1}{2}\right) ={-} \tilde{C}_f \xi k \sinh \left(k \dfrac{a - 1}{2}\right) . \end{equation}
By eliminating 
$\tilde {C}_f$, (4.17a,b) yields
\begin{equation} f'_1 \left(\dfrac{1}{2}\right) + \xi k \tanh \left(k \dfrac{a - 1}{2}\right) f_1 \left(\dfrac{1}{2}\right) = 0 . \end{equation}
On account of (4.10), (4.16) and (4.18), one can conclude that the linear stability analysis can be now formulated by focusing just on the inner core region 
$-1/2 \le x \le 1/2$ and solving the eigenvalue problem
$$\begin{gather} f''_1 - k^2 f_1 - {\rm i} k S h_1 = 0 , \end{gather}$$
$$\begin{gather}h''_1 - (k^2 - {\rm i} \omega + {\rm i} k S x) h_1 + f'_1 = 0, \end{gather}$$
$$\begin{gather}x ={\pm} \dfrac{1}{2}: \quad f'_1 \pm \xi k \tanh \left(k \dfrac{a - 1}{2}\right) f_1 = 0, \quad h_1 = 0. \end{gather}$$
After having solved (4.19) and determined 
$f_1$, one can determine also 
$f_2$ as
\begin{equation} f_2 =\begin{cases} f_1 \left(-\dfrac{1}{2}\right) \left[\cosh \left(k \dfrac{a - 1}{2}\right)\right]^{{-}1} \cosh \left(k \dfrac{a + 2 x}{2}\right), \quad -\dfrac{a}{2} \le x <{-} \dfrac{1}{2}, \\ f_1\left(\dfrac{1}{2}\right) \left[\cosh \left(k \dfrac{a - 1}{2}\right)\right]^{{-}1} \cosh \left(k \dfrac{a - 2 x}{2}\right), \quad \dfrac{1}{2} < x \le \dfrac{a}{2}, \end{cases} \end{equation}4.3. Features of the stability eigenvalue problem
We note that (4.16) and (4.18) are boundary conditions of the third kind for 
$f_1$. However, they are of a special type as the coefficient of the 
$f_1$ term depends on the wavenumber 
$k$. A similar circumstance occurs in the linear stability eigenvalue problems solved by Rees & Mojtabi (Reference Rees and Mojtabi2013) and by Mohammad, Rees & Mojtabi (Reference Mohammad, Rees and Mojtabi2017), where third-kind boundary conditions with a 
$k$-dependent coefficient were found for the temperature disturbances. In particular, in Rees & Mojtabi (Reference Rees and Mojtabi2013), Prats’ problem in a horizontal porous layer (Prats Reference Prats1966) was reconsidered by studying the effect of the non-zero thermal resistance of the horizontal boundary walls. Thus, from a mathematical viewpoint, the circumstances devised by those authors are comparable with those considered here and leading to (4.19). In fact, we are studying a finite and non-zero hydraulic (instead of thermal) resistance of the external porous layers. We mention that the case of third-kind boundary conditions for the temperature was considered in the analysis of Gill's problem for a vertical plane slab with permeable boundaries by Barletta, Celli & Ouarzazi (Reference Barletta, Celli and Ouarzazi2017), while third-kind boundary conditions for the pressure were investigated by Barletta, Celli & Rees (Reference Barletta, Celli and Rees2020). However, in Barletta et al. (Reference Barletta, Celli and Ouarzazi2017, Reference Barletta, Celli and Rees2020), the prescribed third-kind conditions feature constant 
$k$-independent coefficients, unlike the case defined by (4.19). We mention that third-kind boundary conditions with constant coefficients were also predicted for the pressure field by Nygård & Tyvand (Reference Nygård and Tyvand2010).
5. Discussion of the results
The basis for developing the stability analysis is the eigenvalue problem (4.19). Some important characteristics of the onset of instability can be gathered by exploring three significant asymptotic cases.
5.1. The limit 
$\xi \to \infty$
The limiting condition 
$\xi \to \infty$ embodies the case where the external 
${\rm M}_2$ layers are much more permeable than the core 
${\rm M}_1$ layer. We note that the example of metal foams for the 
${\rm M}_2$ layers is quite close to this condition as such metal foams are generally endowed with a large permeability. In any case, we assume that the large permeability of the 
${\rm M}_2$ layers is not so large as to suppress the validity of Darcy's law. In this case, the only change in the stability eigenvalue problem is a marked simplification of the boundary conditions. Equation (4.19) thus reads
$$\begin{gather} f''_1 - k^2 f_1 - {\rm i} k S h_1 = 0 , \end{gather}$$
$$\begin{gather}h''_1 - (k^2 - {\rm i} \omega + {\rm i} k S x) h_1 + f'_1 = 0, \end{gather}$$
$$\begin{gather}x ={\pm} \dfrac{1}{2}: \quad f_1 =0, \quad h_1 = 0. \end{gather}$$
We emphasise that (5.1) coincides with the stability eigenvalue problem for a homogeneous vertical porous layer with permeable boundaries solved by Barletta (Reference Barletta2015). This is an important result as the transition to instability in the limiting case 
$\xi \to \infty$ is defined by the neutral stability data obtained and discussed by Barletta (Reference Barletta2015). On taking the limit 
$\xi \to \infty$, one loses also the dependence on 
$a$.
5.2. The limits $\xi \to 0$ or $a \to 1$
An asymptotic case completely different from that discussed in § 5.1 is defined by the limit 
$\xi \to 0$. This limit describes a situation where the 
${\rm M}_2$ layers are much less permeable than the core 
${\rm M}_1$ layer. Strictly speaking, the external layers are quite close to impermeability if compared with the core layer. In this case, (4.19) markedly simplifies to
$$\begin{gather} f''_1 - k^2 f_1 - {\rm i} kS h_1 = 0 , \end{gather}$$
$$\begin{gather}h''_1 - (k^2 - {\rm i} \omega + {\rm i} k S x) h_1 + f'_1 = 0, \end{gather}$$
$$\begin{gather}x ={\pm} \dfrac{1}{2}: \quad f'_1 = 0, \quad h_1 = 0. \end{gather}$$The stability eigenvalue problem (5.2) coincides with that analysed by Gill (Reference Gill1969). This feature is completely unsurprising as Gill (Reference Gill1969) investigated the possible onset of instability in a homogeneous porous layer with impermeable isothermal boundaries kept at different temperatures. Gill (Reference Gill1969) proved that no instability is possible in this case.
It is significant that (5.2) represents also the limiting case 
$a \to 1$ for any finite 
$\xi$. This result is expected as 
$a \to 1$ means that the thickness of the 
${\rm M}_2$ layers tends to 
$0$. Thus, one has again a homogeneous porous layer, namely the 
${\rm M}_1$ layer, with the impermeable boundaries devised in Gill (Reference Gill1969).
5.3. The limit $a \to \infty$
When the external 
${\rm M}_2$ layers have an extremely large thickness, so that 
$D \gg L$, we have in mind the core 
${\rm M}_1$ layer surrounded by infinite 
${\rm M}_2$ media. This condition yields the limit 
$a \to \infty$ and the stability eigenvalue problem (4.19) simplifies to
$$\begin{gather} f''_1 - k^2 f_1 - {\rm i} k S h_1 = 0 , \end{gather}$$
$$\begin{gather}h''_1 - (k^2 - {\rm i} \omega + {\rm i}k S x) h_1 + f'_1 = 0, \end{gather}$$
$$\begin{gather}x ={\pm} \dfrac{1}{2}: \quad f'_1 \pm \xi k f_1 = 0, \quad h_1 = 0. \end{gather}$$5.4. The neutral stability curves
The solution of (4.19) leads to the determination of the neutral stability curves, namely the curves drawn in the 
$(k,S)$ plane which describe the condition of zero growth rate. In other words, the neutral stability curves are isolines of the imaginary part of the complex parameter 
$\omega$ corresponding to a zero value. Among the many neutral stability curves existing in the 
$(k,S)$ plane, our attention is focused on that displaying the lowest values of 
$S$ when input values of 
$a$ and 
$\xi$ are prescribed. In fact, the lowest neutral stability curve in the 
$(k,S)$ plane captures the parametric condition for the initiation of the instability. As this curve usually features an absolute minimum of 
$S$ for a given 
$k$, this minimum yields the threshold for the linear convective instability. The values of 
$k$ and 
$S$ for such a minimum are called critical and denoted with 
$k_c$ and 
$S_c$. In order to achieve these results, (4.19) must be solved numerically. There are several techniques available for the solution of stability eigenvalue problems formulated through systems of ordinary differential equations. Comprehensive surveys of and comparisons among the methods are available in Dongarra, Straughan & Walker (Reference Dongarra, Straughan and Walker1996) and in Straughan & Walker (Reference Straughan and Walker1996).
Our analysis of the neutral stability curves and of the critical values is carried out by employing a shooting method solution of (4.19) for input values of 
$(a, \xi )$. The code employed for the implementation of such a method is just an adaptation to the diverse boundary conditions of that described in Barletta (Reference Barletta2015). In that paper, the test of the numerical accuracy for this code is also discussed. An alternative and more extensive presentation of the shooting method for linear stability eigenvalue problems is also available in Barletta (Reference Barletta2019).
It is desirable to see the determined critical value of 
$S$, for given 
$(a, \xi )$, formulated with the Rayleigh number 
$R$ in order to fully understand the mode selection at the onset of the linear instability. This information is easily gathered from (4.9), as this equation clearly shows that the least value of 
$R_c$ able to reproduce a given 
$S_c$ is obtained when 
$k_z = k$. This means that the preferred modes at the onset of instability are those with 
$k_y = 0$, the so-called transverse normal modes.
The numerical solution of (4.19) leads to a first important remark: the linear transition to instability occurs with non-travelling modes. This means that the neutral stability curves, defined by a zero imaginary part of 
$\omega$, are characterised also by a zero real part of 
$\omega$, so that the phase velocity of the neutrally stable modes is zero. This feature is justified by the numerical solution of the eigenvalue problem (4.19) for several input data 
$(a, \xi )$.
Figures 3–5 display different frames for different values of 
$\xi$. In each frame, the neutral stability curves are displayed for distinct values of 
$a$, with the asymptotic case 
$a \to \infty$ reported for comparison as a dotted black line. Figures 3 and 4 also show the asymptotic case 
$\xi \to \infty$ as a dashed grey line for a comparison with the situation examined in Barletta (Reference Barletta2015). Figure 3, relative to 
$\xi = 100$ and 
$10$, displays an evident cluttering of the curves close to the asymptote 
$a \to \infty$ for the largest values of 
$a$. For the case 
$\xi =100$, the largest values of 
$a$ mean 
$a > 1.5$. In this case, there is no visible distinction between the asymptotes 
$\xi \to \infty$ and 
$a \to \infty$. We can interpret these findings by saying that, with 
$\xi = 100$, the analysis carried out in Barletta (Reference Barletta2015) yields a fair description of the linear onset of instability for thickness ratios 
$a$ down to 
$2$ or even smaller. Things are different when 
$\xi =10$ as one finds a marked difference between the asymptotes 
$\xi \to \infty$ and 
$a \to \infty$. This phenomenon turns out to be even more evident by exploring figure 4 with the cases 
$\xi = 5$ and 
$\xi = 2$. In figure 5, relative to 
$\xi = 1$ and 
$0.6$, the dashed grey line for 
$\xi \to \infty$ is not even drawn, as such cases are too far from the condition of isobaric boundary conditions, 
$f_1 = 0$, devised in Barletta (Reference Barletta2015).
Figure 3. Neutral stability curves in the 
$(k,S)$ plane for (a) 
$\xi =100$ and (b) 
$\xi = 10$ with different values of 
$a$ (solid lines). The dashed grey line corresponds to the limit 
$\xi \to \infty$, while the dotted black line describes the limit 
$a \to \infty$.
Figure 4. Neutral stability curves in the 
$(k,S)$ plane for (a) 
$\xi =5$ and (b) 
$\xi = 2$ with different values of 
$a$ (solid lines). The dashed grey line corresponds to the limit 
$\xi \to \infty$, while the dotted black line describes the limit 
$a \to \infty$.
Figure 5. Neutral stability curves in the 
$(k,S)$ plane for (a) 
$\xi =1$ and (b) 
$\xi = 0.6$ with different values of 
$a$ (solid coloured curves). The dotted black curves describe the limit 
$a \to \infty$.
There is a systematic trend gathered from figures 3–5 which ought to be emphasised. The basic flow tends to be destabilised as 
$a$ or 
$\xi$ increases. This conclusion can be inferred from figures 3–5 as the neutral stability curve moves upward when the value of either 
$a$ or 
$\xi$ decreases. In fact, if the neutral stability curve moves upward, larger and larger values of 
$S$ are needed to trigger the instability. There is a straightforward interpretation for such a trend. When 
$a$ gradually decreases, the overall width of the multilayer structure, 
$D$, tends to approach the width of the core layer, 
$L$. As a consequence, the stabilising effect of the impermeable external boundaries at 
$x=\pm D/2$, proved by Gill (Reference Gill1969), tends to affect directly the interfaces at 
$x=\pm L/2$, so that the core layer boundaries become gradually close to impermeability and, hence, to stability for every value of 
$S$. Just the same condition happens when the permeability ratio, 
$\xi$, tends to zero as demonstrated mathematically in § 5.2. It is also evident from figures 3–5 that the instability tends to be caused by smaller and smaller wavenumbers, 
$k$, for decreasing values of either 
$a$ or 
$\xi$.
5.5. Critical conditions for the instability
The minimum of the neutral stability curve in the 
$(k,S)$ plane identifies the onset of the linear instability, with the least value of 
$S$ needed for the convection cells to display a growing amplitude in time. Such a minimum identifies the critical values 
$k=k_c$ and 
$S=S_c$, which depend on both 
$a$ and 
$\xi$. The critical data are presented graphically in figures 6 and 7 where the trends of 
$S_c$ and 
$k_c$ versus 
$\xi$ are displayed. These figures reveal that the cases 
$a=20$, 10, 5 yield hardly distinguishable curves almost overlapped with the asymptotic dotted black line corresponding to the limit 
$a \to \infty$. A departure from this behaviour is found only for small values of 
$\xi$ in a range where the singularity of 
$S_c$ is displayed in figure 6 with a steep increase of the critical value of 
$S$ as 
$\xi$ decreases. We also report that the asymptotic regime 
$\xi \to \infty$ is approached by 
$S_c$ more and more rapidly the larger is the value of 
$a$. Thus, one can say that the case examined by Barletta (Reference Barletta2015) is one naturally emerging when both 
$a$ and 
$\xi$ are large enough. For instance, a glance at figure 6 may suggest 
$a \ge 5$ and 
$\xi \ge 20$ as a possible indication where the model of permeable boundaries employed by Barletta (Reference Barletta2015) can be considered as a fair enough approximation for practical purposes. Figures 6 and 7 also reveal that, for every 
$a>0$, 
$S_c$ becomes infinite and 
$k_c$ becomes zero when 
$\xi \to \xi _a > 0$. The constant 
$\xi _a$ is a decreasing function of 
$a$. The evaluation of 
$\xi _a$ is achieved by the solution of the stability eigenvalue problem in a large-
$S$ range where the accuracy of the numerical solver significantly decreases. A few values of 
$\xi _a$ relative to the values of 
$a$ considered in figures 6 and 7 are shown in table 1. Only three significant figures could be reported in this table due to the mentioned decrease in numerical accuracy. An important result is that no linear instability is found when 
$\xi \le \xi _a$.
Figure 6. Critical value of 
$S$ versus 
$\xi$ for 
$a$ decreasing from 
$20$ to 
$1.2$ (solid coloured curves). The dotted black curve describes the limit 
$a \to \infty$, while the dashed grey line denotes the asymptotic value for 
$\xi \to \infty$, 
$S_c = 197.0812$.
Figure 7. Critical value of 
$k$ versus 
$\xi$ for 
$a$ decreasing from 
$20$ to 
$1.2$ (solid coloured curves). The dotted black curve describes the limit 
$a \to \infty$, while the dashed grey line denotes the asymptotic value for 
$\xi \to \infty$, 
$k_c = 1.059498$.
Table 1. Values of 
$\xi _a$ versus 
$a$.
Critical values 
$k=k_c$ and 
$S=S_c$ are given in tables 2 and 3 versus 
$\xi$ for different decreasing values of 
$a$, from the asymptotic condition 
$a \to \infty$ to 
$a=1.2$. Critical data are not reported for the cases where 
$\xi \le \xi _a$. In fact, as pointed out above, a threshold to instability is found only for a 
$\xi$ larger than 
$\xi _a$.
Table 2. Critical values of 
$k$ and 
$S$.
Table 3. Critical values of 
$k$ and 
$S$.
5.6. Streamlines and isotherms
The visualisation of the streamlines and isotherms associated with the perturbation modes is quite important, especially at the critical conditions 
$k=k_c$ and 
$S=S_c$ defining the initiation of the convective instability. In fact, even if the flow patterns of convection are evaluated according to the linear theory, it is generally retained that such patterns are the starting point for the development of the nonlinear analysis of convection. This is particularly evident in the formulation of the weakly nonlinear theory of convective instability (Drazin & Reid Reference Drazin and Reid2004). In fact, the weakly nonlinear theory studies the interaction between normal modes in the vicinity of the critical condition.
Starting from (2.9) and (4.1a,b), one can define the perturbation velocity:
$$\begin{gather} \hat{\boldsymbol{u}}_1 ={-} \boldsymbol{\nabla}\hat{P}_1 + R \hat{T}_1 \hat{\boldsymbol{e}}_z , \end{gather}$$
$$\begin{gather}\hat{\boldsymbol{u}}_2 ={-} \xi \boldsymbol{\nabla}\hat{P}_2 + R \xi \hat{T}_2 \hat{\boldsymbol{e}}_z . \end{gather}$$
We mentioned that the onset of instability occurs with transverse modes, i.e. the two-dimensional perturbations in the 
$(x,z)$ plane with 
$k_y=0$ and 
$k=k_z$. Accordingly, only the velocity components 
$(\hat {u}_m, \hat {w}_m)$, with 
$m=1,2$, are non-zero for such modes. The mass balance constrains 
$\hat {\boldsymbol {u}}_m$ to be solenoidal. Thus, one can introduce the streamfunction 
$\hat {\varPsi }_m$, so that the divergence of 
$\hat {\boldsymbol {u}}_m$ is identically zero:
\begin{equation} \hat{u}_m = {\displaystyle{{\partial {\hat{\varPsi}_m}}\over{\partial z}}},\quad \hat{w}_m ={-}{\displaystyle{{\partial {\hat{\varPsi}_m}}\over{\partial x}}}. \end{equation}Following (4.5), we express the streamfunction as a normal mode:
\begin{equation} \hat{\varPsi}_m = \psi_m(x) \, {\rm e}^{{\rm i} ( k z - \omega t)}. \end{equation}Then, (4.5) and (5.4)–(5.6) yield
$$\begin{gather} {\rm i} k \psi_1(x) ={-} f'_1(x), \quad - \psi'_1(x) ={-}{\rm i} kf_1(x) + R h_1(x), \end{gather}$$
$$\begin{gather}{\rm i} k \psi_2(x) ={-} \xi f'_2(x),\quad - \psi'_2(x) ={-}{\rm i} k \xi f_2(x) + R\xi h_2(x). \end{gather}$$
As a consequence, the interface conditions at 
$x=\pm 1/2$ expressed by (4.8) imply a continuity of 
$\psi _m$ across such interfaces with a discontinuity in its first derivative when 
$\xi \ne 1$:
\begin{equation} x ={\pm} 1/2 : \quad \psi_1 = \psi_2 , \quad \xi \psi'_1 = \psi'_2 . \end{equation}
By employing (5.6) and (5.8), one realises that 
$\boldsymbol {\nabla } \hat {\varPsi }_m$ undergoes a cusp discontinuity at each interface where its 
$x$ component changes discontinuously, if 
$\xi \ne 1$, while the 
$z$ component is always continuous. Such cusp discontinuities of the streamlines at the interfaces are clearly visible in figures 8 and 9, except for the case 
$a=5$ and 
$\xi =1$ reported in the upper part of figure 9. These figures illustrate just four sample cases in order to display how the convective instability induces a cellular flow which penetrates the 
${\rm M}_2$ layers. In such external layers, the saturated porous medium is otherwise isothermal, at equilibrium with the external impermeable boundaries, 
$x=\pm a/2$. The latter feature emerges quite clearly from figures 8 and 9 as the cellular patterns for the isotherms are entirely confined in the core 
${\rm M}_1$ layer, that is, within the region 
$-1/2 \le x \le 1/2$. The red/blue colour code for the different media introduced in figure 1 is used also in figures 8 and 9 for convenience. We finally point out that the antisymmetry of the convection cells with respect to the vertical midplane reflects the antisymmetric thermal boundary conditions and the antisymmetric basic temperature profile triggering the onset of the convection cells.
Figure 8. Perturbation streamlines and isotherms for transverse modes (
$k_y=0$, 
$k_z=k$ and 
$S=R$) with 
$a=2$ and either 
$\xi =5$ or 
$\xi =50$, at critical conditions 
$k=k_c$ and 
$R=S=S_c$. The dotted black line is at 
$x=0$, while the dashed lines denote the interfaces 
$x=\pm 1/2$. The vertical range is over a period 
$0 \le z \le 2{\rm \pi} /k_c$.
Figure 9. Perturbation streamlines and isotherms for transverse modes (
$k_y=0$, 
$k_z=k$ and 
$S=R$) with 
$a=5$ and either 
$\xi =1$ or 
$\xi =50$, at critical conditions 
$k=k_c$ and 
$R=S=S_c$. The dotted black line is at 
$x=0$, while the dashed lines denote the interfaces 
$x=\pm 1/2$. The vertical range is over a period 
$0 \le z \le 2{\rm \pi} /k_c$.
6. Conclusions
The onset of convection in a vertical porous slab has been analysed by assuming a three-layer structure. In fact, an internal layer is sandwiched between two identical external layers having properties different from those of the core. The multilayer slab is saturated by a Newtonian fluid and the external layers are much more thermally conductive than the inner core. The multilayer slab is bounded by impermeable isothermal walls kept at different temperatures.
The dynamics of convection is governed by Darcy's law and by the thermal buoyancy modelled through the Boussinesq approximation. Being extremely conductive, the external layers are thermally passive but they can be penetrated by the convection cells arising in the core layer. The basic state features a vertical flow, with a piecewise-linear velocity profile, in a thermal conduction regime. Such a basic state reproduces in the inner layer that envisaged by both Gill (Reference Gill1969) and Barletta (Reference Barletta2015). Those authors considered a single-layer homogeneous slab with either impermeable boundaries (Gill Reference Gill1969) or permeable boundaries (Barletta Reference Barletta2015).
The linear perturbations of the basic flow have been studied by employing a modal analysis. The neutral stability curves have been obtained in the 
$(k, S)$ plane, where 
$k$ is the wavenumber and 
$S$ is a suitably rescaled Rayleigh number, which coincides with the Rayleigh number 
$R$ in the special case of two-dimensional transverse modes. The neutral stability condition depends on two governing parameters: the permeability ratio of the outer layers to the inner layer, 
$\xi$; and the width ratio of the whole slab to the inner layer, 
$a$. The critical conditions 
$k=k_c$ and 
$S=S_c$, corresponding to the point of minimum 
$S$ along a neutral stability curve, have been also evaluated.
The main results obtained from the linear stability analysis are the following:
• The most unstable perturbation modes are transverse. The transverse modes are two-dimensional lying in the vertical
$(x,z)$ plane, where the $x$ axis is horizontal and perpendicular to the slab and the $z$ axis is vertical.• The basic flow is destabilised by both the increase of
$a$ and the increase of $\xi$. The most unstable parametric set-up is the limiting case $\xi \to \infty$, where the neutral stability condition is not influenced by the value of $a$ as the stability eigenvalue problem coincides with that laid out and solved by Barletta (Reference Barletta2015).• For every value of
$a$, there exists a positive constant $\xi _a$ such that no linear instability is found for $\xi \le \xi _a$. The constant $\xi _a$ is a decreasing function of $a$.• The limiting case
$\xi \to 0$ drives the basic flow to linear stability for every $a$ and for every Rayleigh number. Such an asymptotic condition is that devised in the rigorous proof of stability presented by Gill (Reference Gill1969) in his classic paper. The interfaces between the internal layer and the external layers become, in fact, perfectly impermeable in this case.• Since the external layers have a thermal conductivity much larger than that of the core layer, they are isothermal with the same temperature as the neighbouring boundary. This happens both for the basic state and for the perturbed flow. Thus, the external layers are thermally passive at onset of instability. When the convection cells emerge at supercritical conditions, the streamlines of such cells penetrate the external layers showing up cusp discontinuities at the interfaces. The origin of the cusp discontinuities relies on the difference between the permeability of the external layers relative to that of the core layer.
The analysis carried out in this paper points out that the multilayer structure of a vertical porous slab with impermeable isothermal boundaries may induce the onset of convective instability whereas a single-layer vertical slab with the same boundary conditions displays no instability. This result points out a new connection between two apparently different types of instability: the instability in a vertical porous slab with a heterogeneous permeability structure and that in a homogeneous vertical slab with permeable boundaries. The former case, relative to a continuous type of heterogeneity, has been examined by Shankar & Shivakumara (Reference Shankar and Shivakumara2022). On the other hand, the latter case is retrieved from our analysis when the asymptotic condition 
$\xi \to \infty$ is examined, perfectly matching the mathematical formulation and the numerical results presented in Barletta (Reference Barletta2015).
Funding
The authors acknowledge financial support from grant no. PRIN 2017F7KZWS provided by the Italian Ministero dell'Istruzione, dell'Università e della Ricerca.
Declaration of interests
The authors report no conflict of interest.
Open access
