1 Introduction

Natural convection arising by buoyancy-driven instability in a porous medium heated from below is usually called Darcy–Bénard convection. It is a relatively simple type of hydrodynamic stability. One simplifying element is the finite-amplitude undisturbed state with a linear temperature profile. Another simplifying element is the linear Darcy law with no other spatial differential operators than the gradient. The only nonlinearity is the convective term in the heat equation, which is linearized in the eigenvalue problem. The elementary onset problem is referred to as the Horton–Rogers–Lapwood (HRL) problem (Horton and Rogers 1945; Lapwood 1948).

The Darcy–Bénard eigenvalue problem is of fourth order. This coupled thermomechanical problem is physically more complicated than the mathematically related fourth-order problem of freely vibrating elastic plates. The theory for fourth-order problems is not highly developed because the biharmonic operator does not separate in orthogonal coordinate systems. Analytical methods for the Darcy–Bénard problem are available when the homogeneous boundary conditions are compatible with normal modes, inducing separation of variables. The price to pay for obtaining separability is some kind of degeneracy for the eigenvalue problem.

Solvable degenerate problems dominate the literature and give a background for understanding non-separable problems. A subtle degeneracy was discovered by Lyubimov (1975), finding two sets of eigenfunctions with coinciding critical Rayleigh numbers for two-dimensional (2D) porous cavities with conducting and impermeable walls. Nilsen and Storesletten (1990) calculated these thermomechanical eigenfunctions for a rectangular geometry, confirmed and extended by Rees and Lage (1996). Rees and Tyvand (2004) generalized the linear theory for 2D porous cavities with impermeable and conducting wall.

Mathematical difficulties arising from non-separability of the eigenvalue problem have a low priority in the literature. Tyvand and Storesletten (2018) showed that degeneracy is needed to obtain separable problems for vertical porous cylinders. Tyvand et al. (2019) solved the onset problem for a 2D porous rectangle with boundary conditions that defy analytical solutions.

Our onset problem has an additional length parameter compared with the HRL problem. The eigenvalue problem defies analytical solutions because its eigenfunctions do not separate. The uniform horizontal porous layer of infinite horizontal extent is designed with abrupt switches between Dirichlet and Neumann conditions along the upper and lower boundaries. A spatially concentrated instability is surrounded by penetrative convection laterally, beyond the two central borderlines where the boundary conditions switch.

Convection onset penetrating vertically has been studied by various authors. Convection in a porous layer with unstable thermal stratification may penetrate into a neighboring stable layer. This phenomenon may occur when the saturating fluid (water) has a density maximum, see Straughan (2004). Internal heat sources/sinks may also induce vertically penetrative convection. The book by Nield and Bejan (1998) describes thermally driven natural convection with penetration into passive domains.

McKibbin (1986) started the investigation of convection instability which is penetrative horizontally. He added physical realism for geothermal models by considering vertical porous layers with different properties. He maintained the confinement between two horizontal planes of the HRL problem and found locally triggered instability, penetrating laterally into locally stable domains. Rees and Tyvand (2009) performed a similar study with periodically varying permeability. Their 2D analysis was extended to 3D by Rees and Barletta (2014).

Ahmad and Rees (2016) considered the Laplacian thermal fields in solid blocks surrounding a porous box with convection onset. The two conducting blocks in their problem are similar to two surrounding isothermal domains where Darcy flow without buoyancy takes place. The Laplacian fields in the domains surrounding the unstable domain have passive spatial damping.

These previous papers considered heterogeneity in the horizontal direction, triggering local convection penetrating into more stable surroundings. The present model is different, with a homogeneous porous layer of infinite horizontal extent. Only the variable Dirichlet or Neumann conditions along the upper and lower boundaries are responsible for the penetrative convection. We will let a central domain of the layer have boundary conditions that trigger a local instability.

We will investigate how the onset criterion depends on the width of the unstable domain, and we will investigate the spatial decay of the thermomechanical field penetrating into the surrounding stable domains.

2 Mathematical Formulation of the Designed Physical Model

The present model is designed for the purpose of developing a simple spatially localized marginal state of convection in a porous layer. We choose the simplest possible geometry for a Darcy–Bénard onset problem, which is a uniform horizontal porous layer of infinite horizontal extent. The convective instability is concentrated in a small spatial domain without changing the simple geometry of a horizontal porous layer of unlimited extent. We will let the homogeneous thermomechanical boundary conditions change from Dirichlet to Neumann conditions at appropriate locations along the lower and upper boundaries of the porous layer. The marginal state of convection will then be a state of penetrative convection, with spatial decay as \(\vert x \vert \rightarrow \infty\).

The lower and upper boundaries of the horizontal porous layer are \(z=-\,h/2\) and \(z=h/2\). The z axis is directed vertically upward, where g denotes the gravitational acceleration. The velocity vector \({\mathbf{v}}\) has Cartesian components (u, v, w). The temperature field is T(x, y, z, t), with t denoting time. There is an undisturbed motionless state with a uniform vertical temperature gradient. The gravitational acceleration is written in vector form as \({\mathbf{g}}\). The porous medium is homogeneous and isotropic, with permeability K. The standard Darcy–Boussinesq equations for free thermal convection can be written

$$\begin{aligned}&\nabla p +\frac{\mu }{K} {{\mathbf{v}}} + \rho _0 \beta ~ (T - T_0) {\mathbf{g}}=0, \end{aligned}$$
(1)
$$\begin{aligned}&\nabla \cdot {\mathbf{v}}=0, \end{aligned}$$
(2)
$$\begin{aligned}&(\rho c_p)_m\frac{\partial T}{\partial t}+(\rho c_p)_f~ {\mathbf{v}}\cdot \nabla T = \lambda _m \nabla ^2 T. \end{aligned}$$
(3)

In these equations, p is the dynamic pressure, \(\beta\) is the coefficient of thermal expansion, \(\rho =\rho _0\) is the fluid density at the reference temperature \(T_0\), \(\mu\) is the dynamic viscosity of the saturating fluid, \(c_p\) is the specific heat at constant pressure, and \(\lambda _m\) is the thermal conductivity of the saturated porous medium. The subscript m refers to an average over the solid/fluid mixture, while the subscript f refers to the saturating fluid alone.

A semi-infinite left-hand portion of the lower boundary is taken to be perfectly heat-conducting and impermeable

$$\begin{aligned} T = T_0,\quad w = 0,\quad {\text{ at }}\quad z=-\frac{h}{2},\quad x<-\frac{l}{2}. \end{aligned}$$
(4)

where l is a length parameter that is defined below. The semi-infinite right-hand portion of the lower boundary is subject to a given heat flux, and with constant pressure (open boundary)

$$\begin{aligned} \frac{\partial T}{\partial z} = - \frac{\Delta T}{h},\quad {\mathbf{k}} \times {\mathbf{v}} = 0,\quad {\text{ at }}\quad z=-\frac{h}{2},\quad x>-\frac{l}{2}. \end{aligned}$$
(5)

This dynamic constant-pressure condition implies that \({\mathbf{k}} \times \nabla p = 0\) along a semi-infinite portion of the lower boundary \(z=-\,h/2\), and the cross-product of Darcy’s law (1) implies this flow condition of zero tangential velocity. \(\Delta T\) is the temperature difference across the layer in its undisturbed state. \(T_0\) is a reference temperature.

In our designed boundary conditions, we introduce a horizontal displacement length l for applying the same types of boundary conditions at the upper boundary. The length l is the horizontal distance between two straight lines: (1) the horizontal line \(x=-\,l/2\) where the boundary conditions switch from Dirichlet to Neumann type along the lower boundary. (2) The horizontal line \(x=l/2\) where the boundary conditions switch from Neumann to Dirichlet type along the upper boundary.

The semi-infinite right-hand portion of the displaced upper boundary is taken to be perfectly heat-conducting and impermeable

$$\begin{aligned} T = T_0 - \Delta T,\quad w = 0,\quad {\text{ at }}\quad z=\frac{h}{2},\quad x>\frac{l}{2}. \end{aligned}$$
(6)

The semi-infinite right-hand portion of the displaced upper boundary is subject to a given heat flux, and with constant pressure (open boundary)

$$\begin{aligned} \frac{\partial T}{\partial z} = - \frac{\Delta T}{h},\quad {\mathbf{k}} \times {\mathbf{v}} = 0,\quad {\text{ at }}\quad z=\frac{h}{2},\quad x<\frac{l}{2}. \end{aligned}$$
(7)

Thus, the boundary conditions for the portion \(x<-l/2\) of the lower plane \(z=-\,h/2\) are handpicked to be the same as the boundary conditions for the portion \(x>l/2\) of the upper plane \(z=h/2\). The same is true for the semi-infinite portion \(x>-l/2\) of the lower plane and the portion \(x>l/2\) of the upper plane. This means that the whole configuration of the porous layer has an antisymmetry around the origin, which is a reason that we have placed the origin in the middle of the porous layer.

Our open-boundary condition of constant pressure along semi-infinite horizontal planes must agree with establishing an undisturbed state of rest with a linear temperature profile. A nonzero buoyancy force at these open boundaries may induce a basic vertical flow. Therefore, we assume a uniform fluid layer of small thickness \(\delta\) compared with the layer thickness h, below and below each open boundary. Each layer of saturating fluid may itself be a porous layer of permeability much greater than K, next to an impermeable bottom. Mass balance forbids a vertical throughflow in the basic state. The fluid layers \(-\,\delta< z + h/2 < 0\) below the porous medium and \(0< z-h/2 < \delta\) above the porous medium serve one single purpose. This purpose is to prevent horizontal pressure gradients from arising along the open portions of these boundaries. The rigid exterior horizontal planes confining these open boundaries (\(z + h/2=-\,\delta ,\, x<0\) and \(z-h/2 = \delta , x>0\)) serve to block the basic vertical flow that might otherwise be present in a basic state of pure conduction with hydrostatic pressure. A perturbation velocity is allowed in order to maintain the requirement of constant pressure along an open boundary. The thin layers above and below an open (penetrative) boundary must be thick enough to serve as hydrostatic pressure reservoirs where pressure fluctuations are absorbed to prevent horizontal pressure gradients from arising.

In Fig. 1, we show two sketches of the physical model, where we omit the thin neighboring layers that guarantee for the realism of the open-boundary condition. The horizontal x, y plane is placed in the middle of the porous layer, with a vertical z axis. The mathematical symbols in Fig. 1 will be introduced and explained below.

Fig. 1
figure 1

Definition sketches of a porous layer of unit thickness and of infinite horizontal extent. These are perspective sketches of perpendicular cuts through the physical model, showing the boundary conditions for the 3D eigenfunctions \(\varTheta\) and \(\varPsi\) at the lower and upper plane. The upper sketch represents the perturbation temperature \(\varTheta\). The lower sketch represents the vector potential \(\varPsi\). Mixed boundary conditions apply at the lower boundary \(z=1/2\) and the upper boundary \(z=1/2\). Neumann conditions are marked in red, while Dirichlet conditions are marked in blue. The vertical z axis is directed in the direction opposite of gravity (g). The thermomechanical boundary conditions switch between the Dirichlet type and Neumann type in the two lines \((x=-\,L/2,z=-\,1/2)\) and \((x=L/2,z=1/2)\), in terms of dimensionless length

Full size image

2.1 Dimensionless Equations

From now on, we work with dimensionless variables. We reformulate the mathematical problem in dimensionless form by means of the transformations

$$\begin{aligned}&\frac{1}{h} (x, y, z) \rightarrow (x, y, z),\quad \frac{h}{\kappa _m }(u, v, w)\rightarrow (u, v, w),\quad h \nabla \rightarrow \nabla ,\\&\frac{1}{\Delta T} (T - T_0) \rightarrow T,\quad \frac{K }{\mu \kappa _m} (p - p_0) \rightarrow p,\quad \frac{(\rho c_p)_f \kappa _m }{(\rho c_p)_m H^2} t \rightarrow t, \end{aligned}$$

where \(\kappa _m = \lambda _m/(\rho _0 c_p)_f\) is the thermal diffusivity of the saturated porous medium. We denote the vertical unit vector by \({\mathbf{k}}\), directed upward.

The dimensionless governing equations can then be written

$$\begin{aligned}&{\mathbf{v}}+\nabla p - R~ T {\mathbf{k}}=0. \end{aligned}$$
(8)
$$\begin{aligned}&\nabla \cdot {\mathbf{v}}=0 \end{aligned}$$
(9)
$$\begin{aligned}&\frac{\partial T}{\partial t} +{\mathbf{v}} \cdot \nabla T = \nabla ^2 T. \end{aligned}$$
(10)

The dimensionless boundary conditions at the lower boundary are

$$\begin{aligned}&T = w = 0,\quad {\text{ at }}\quad z=-\frac{1}{2},\quad x<-\frac{L}{2}, \end{aligned}$$
(11)
$$\begin{aligned}&\frac{\partial T}{\partial z} = - \,1,\quad u = v = 0,\quad {\text{ at }}\quad z=-\frac{1}{2},\quad x>-\frac{L}{2}. \end{aligned}$$
(12)

We give the corresponding conditions at the upper boundary

$$\begin{aligned}&T = -1,\quad w = 0,\quad {\text{ at }}\quad z=\frac{1}{2},\quad x>\frac{L}{2}, \end{aligned}$$
(13)
$$\begin{aligned}&\frac{\partial T}{\partial z} = - \,1,\quad u = v = 0,\quad {\text{ at }}\quad z=\frac{1}{2},\quad x < \frac{L}{2}. \end{aligned}$$
(14)

where we have introduced the dimensionless displacement length \(L = l/h\).

Here, the Rayleigh number R is defined as

$$\begin{aligned} R =\frac{\rho _0 g \beta K \Delta T h}{\mu \kappa _m }. \end{aligned}$$
(15)

2.2 Basic Solution

The stationary basic solution of Eqs. (8)–(14) is given subscript “b”

$$\begin{aligned} {\mathbf{v}}_b = 0, \quad T_b = -\, z,\quad p_b = - \frac{1}{2} R~ z^2. \end{aligned}$$
(16)

This basic state of hydrostatic fluid has a linear temperature profile.

2.3 Linearized Perturbation Equations

In our stability analysis, we disturb the basic state (16) with perturbed fields

$$\begin{aligned} {\mathbf{v}} = {\mathbf{v}}_b + {\mathbf{v}},\quad T = T_b(z) + \varTheta ,\quad P = P_b(z) + p'. \end{aligned}$$
(17)

where the perturbations \({\mathbf{v}}, \varTheta , p'\) are functions of x, y, z and t. Linearizing Eqs. (8)–(10) with respect to perturbations and eliminating the pressure gives

$$\begin{aligned}&\nabla ^2 w = R \nabla _1^2 \varTheta , \end{aligned}$$
(18)
$$\begin{aligned}&\frac{\partial \varTheta }{\partial t} - w = \nabla ^2 \varTheta . \end{aligned}$$
(19)

With vanishing vertical vorticity, one single scalar function \(\varPsi (x,y,z)\) represents the 3D thermomechanical vector field. The incompressible flow field is given by this poloidal vector potential as \({\mathbf{v}} = \nabla \times (\nabla \times {\mathbf{k}} \varPsi )\), where \({\mathbf{k}}\) is the vertical unit vector. The velocity components are

$$\begin{aligned} (u, v, w) = \left( \frac{\partial ^2 \varPsi }{\partial x \partial z}, \frac{\partial ^2 \varPsi }{\partial y \partial z}, - \nabla _1^2 \varPsi \right) , \end{aligned}$$
(20)

where the horizontal Laplacian operator \(\nabla _1^2 = \partial ^2/\partial x^2 + \partial ^2/\partial y^2\) has been introduced. From Eq. (18), it follows that the perturbation temperature is given by

$$\begin{aligned} \varTheta = - \,R^{-1} \nabla ^2 \varPsi . \end{aligned}$$
(21)

Assuming a non-oscillatory marginal state, the heat Eq. (19) becomes

$$\begin{aligned} \nabla _1^2 \varPsi = \nabla ^2 \varTheta . \end{aligned}$$
(22)

The boundary conditions at the lower boundary are

$$\begin{aligned}&\varTheta = \varPsi = 0,\quad {\text{ at }}\quad z=-\frac{1}{2},\quad x<\frac{L}{2}, \end{aligned}$$
(23)
$$\begin{aligned}&\frac{\partial \varTheta }{\partial z} = \frac{\partial \varPsi }{\partial z} = 0,\quad {\text{ at }}\quad z=0,\quad z=-\frac{1}{2},\quad x>\frac{L}{2}. \end{aligned}$$
(24)

The conditions at the upper boundary are

$$\begin{aligned}&\varTheta = \varPsi = 0,\quad {\text{ at }}\quad z=\frac{1}{2},\quad x>\frac{L}{2}, \end{aligned}$$
(25)
$$\begin{aligned}&\frac{\partial \varTheta }{\partial z} = \frac{\partial \varPsi }{\partial z} = 0,\quad {\text{ at }}\quad z=\frac{1}{2},\quad x<\frac{L}{2}. \end{aligned}$$
(26)

The coupled eigenvalue problem for the dimensionless eigenfunctions \(\varTheta\) and \(\varPsi\) is already illustrated in Fig. 1 above. The upper sketch in Fig. 1 shows the mixed thermal boundary conditions, while the lower sketch shows the mixed conditions for the flow. In our designed physical model, the eigenfunctions \(\varTheta\) and \(\varPsi\) obey the same homogeneous boundary conditions everywhere, either of Dirichlet type (with blue color markings) or of Neumann type (with red color markings). The lengths that we have introduced in Fig. 1 are dimensionless.

3 The 3D Eigenvalue Problem Reduced to 2D in x and z

With infinite horizontal extent, the 3D solution includes a Fourier component with wavenumber k in the y direction, prescribing the eigenfunctions

$$\begin{aligned} \varPsi = \chi (x, z) {\hbox {e}}^{i k y},\quad \varTheta = \theta (x, z) {\hbox {e}}^{i k y}, \end{aligned}$$
(27)

where i denotes the imaginary unit. The 2D perturbation equations for the eigenfunctions \(\theta (x,z)\) and \(\chi (x,z)\) are

$$\begin{aligned}&\left( \nabla ^2 - k^2\right) \chi = -\, R~ \theta , \end{aligned}$$
(28)
$$\begin{aligned}&\left( \nabla ^2 - k^2\right) \theta = \left( \frac{\partial ^2}{\partial x^2} - k^2\right) \chi . \end{aligned}$$
(29)

We will equip our 2D eigenvalue problem in the vertical xz-plane with boundary conditions. The conditions at the lower boundary are

$$\begin{aligned}&\theta = \chi = 0,\quad {\text{ at }}\quad z=-\frac{1}{2},\quad x<-\frac{L}{2}, \end{aligned}$$
(30)
$$\begin{aligned}&\frac{\partial \theta }{\partial z} = \frac{\partial \chi }{\partial z} = 0,\quad {\text{ at }}\quad z=-\frac{1}{2},\quad x>-\frac{L}{2} \end{aligned}$$
(31)

The conditions at the upper boundary are

$$\begin{aligned}&\theta = \chi = 0,\quad {\text{ at }}\quad z=\frac{1}{2},\quad x>\frac{L}{2}, \end{aligned}$$
(32)
$$\begin{aligned}&\frac{\partial \theta }{\partial z} = \frac{\partial \chi }{\partial z} = 0,\quad {\text{ at }}\quad z=\frac{1}{2},\quad x<\frac{L}{2}. \end{aligned}$$
(33)

3.1 The 2D Eigenvalue Problem Where We Introduce a Streamfunction

We will work with a modified eigenfunction \(\psi (x,z)\), which is the streamfunction in the x, z plane when \(k=0\). It is defined by

$$\begin{aligned} \psi = \frac{\partial \chi }{\partial x}. \end{aligned}$$
(34)

The perturbation equations for the eigenfunctions \(\theta (x,z)\) and \(\psi (x,z)\) are

$$\begin{aligned}&(\nabla ^2 - k^2) \psi = -\, R~ \frac{\partial \theta }{\partial x}, \end{aligned}$$
(35)
$$\begin{aligned}&(\nabla ^2 - k^2) \frac{\partial \theta }{\partial x} = \left( \frac{\partial ^2}{\partial x^2} - k^2\right) \psi . \end{aligned}$$
(36)

A 2D eigenvalue problem is formulated in the vertical xz-plane, to be equipped with boundary conditions. The conditions at the lower boundary are

$$\begin{aligned} \theta = \psi = 0,\quad {\text{ at }}\quad z=-\frac{1}{2},\quad x<-\frac{L}{2}, \end{aligned}$$
(37)
$$\begin{aligned} \frac{\partial \theta }{\partial z} = \frac{\partial \psi }{\partial z} = 0,\quad z=-\frac{1}{2},\quad x>-\frac{L}{2}, \end{aligned}$$
(38)

The conditions at the upper boundary are

$$\begin{aligned}&\theta = \psi = 0,\quad {\text{ at }}\quad z=\frac{1}{2},\quad x>\frac{L}{2}, \end{aligned}$$
(39)
$$\begin{aligned}&\frac{\partial \theta }{\partial z} = \frac{\partial \psi }{\partial z} = 0,\quad {\text{ at }}\quad z=1,\quad z=\frac{1}{2},\quad x<\frac{L}{2}. \end{aligned}$$
(40)

Figure 2 illustrates the 2D eigenvalue problem valid for the 3D physical model sketched in Fig. 1. The upper sketch introduces the mixed thermal conditions and the differential equation for \(\theta\). The lower sketch shows the corresponding differential equation and boundary conditions for \(\psi\). As in Fig. 1, the Dirichlet conditions are marked in blue, with the Neumann conditions in red. Figure 2 gives a 2D sketch showing the rectangular box defined by \(|x| < |L|/2\), designed for triggering local convection, intruding the two neighboring domains \(|x| > |L|/2\) as penetrative convection. The designed boundary conditions have antisymmetry around the origin.

Fig. 2
figure 2

Definition sketches of a vertical cross section for a porous layer of unit thickness and of infinite horizontal extent. The 3D eigenvalue problem is given a dimensionless 2D formulation in the x, z plane, where the perpendicular y direction is represented by the wavenumber k. Mixed boundary conditions apply at the lower boundary \(z=1/2\) and the upper boundary \(z=1/2\). The thermomechanical boundary conditions switch between the Dirichlet type and Neumann type in the two points with dimensionless positions \((-\,L/2,-\,1/2)\) and (L/2, 1/2)

Full size image

The 3D temperature field is

$$\begin{aligned} \varTheta = {\hbox {e}}^{i k y} \theta (x,z), \end{aligned}$$
(41)

with the 3D velocity components

$$\begin{aligned} u&= {\hbox {e}}^{i k y} \frac{\partial \psi }{\partial z} (x,z), \end{aligned}$$
(42)
$$\begin{aligned} v&= i k {\hbox {e}}^{i k y} \int \frac{\partial \psi }{\partial z} (x,z) {\mathrm{d}}x, \end{aligned}$$
(43)
$$\begin{aligned} w&= {\hbox {e}}^{i k y} \left( k^2 \int \psi (x,z) {\mathrm{d}}x - \frac{\partial \psi }{\partial x} (x,z) \right) , \end{aligned}$$
(44)

with real parts representing physical quantities.

4 Two-Dimensional Convection \((k=0)\)

We will first study 2D convection where \(k=0\), before looking for the possibility of an even lower Rayleigh number when k is nonzero. The 2D marginal states include streamline patterns, not available in 3D.

We will now present some numerical results for the marginal state of 2D convection. Figure 3 shows the critical Rayleigh number as a function of the displacement length L for \(k=0\).

Fig. 3
figure 3

Numerically calculated normalized Rayleigh numbers \(R/\pi ^2\) as a function of L for 2D convection with \(k=0\). The ten lowest modes of marginal stability are represented. The onset modes are numbered from 1 to 10 and classified in increasing order according to their Rayleigh numbers (ranked according to the numerical solution at each displacement length)

Full size image

In Fig. 4 we display the 2D thermomechanical eigenfunctions for gradually reduced displacement lengths: \(L=2\), \(L=1\), \(L=0\), \(L=-\,1\) and \(L=-\,2\). In all these cases, we have applied a numerical truncation. Dirichlet and/or Neumann conditions are applied at the lateral truncation boundaries \(x=\pm \,X\), since the numerical method does not handle the assumed infinite extent (\(X \rightarrow \infty\)). We have chosen \(X=5\) in Fig. 4, where the complete computational domain \(-\,X<x<X\) is displayed.

Fig. 4
figure 4

The lowest 2D eigenfunctions (with \(k=0\)) for the displacement lengths \(L=2\), \(L=1\), \(L=0\), \(L=0\), \(L=-\,1\), \(L=-\,2\). The thermal boundary condition at the numerical truncation walls \(x=\pm \, X\) is shown in the figure. When \(L \le 0\), we include both the Dirichlet condition and the Neumann condition. The truncation length is chosen as \(X=5\)

Full size image

We note the concentrated and closed streamlines near the two points \((x,z) = (-\,L/2,-\,1/2)\) and \((x,z) = (L/2,1/2)\) where the boundary conditions change abruptly, indicating singularities in the eigenfunctions. The thermal cell walls closest to the unstable domain have curvature, revealing that the temperature perturbation is non-separable in space. Only one thermal cell wall on each side has visible curvature, as the neighboring cell walls for the streamlines are almost vertical.

The situation where \(L \le 0\) is different, because the central domain around \(x=0\) will be stabilized. The instability will be triggered for \(|x| > |L|/2\), and the eigenfunctions will not vanish at infinity. The artificial sidewalls of numerical truncation \(x = \pm \,X\) will then overrule the eigenfunctions by setting their phase in the far field, and the critical Rayleigh number will be exactly \(R_c = \pi ^2\) for all values \(L \le 0\) when the horizontal porous layer extends to infinity. We include two different sets of boundary conditions for the plots with \(L \le 0\) in Fig. 4, denoting the applied thermal condition in each case. When we apply a Dirichlet condition for \(\theta\) at the truncation walls \(x = \pm \,X\), the conditions for \(\psi\) are of Neumann type. The opposite case of a thermal Neumann condition is accompanied by a Dirichlet condition for the streamfunction, corresponding to normal-mode-type separable eigenfunctions in the far field. The numerical values of \(R_c\) are slightly greater than \(\pi ^2\) when \(L \le 0\), for several reasons: (1) The stabilizing influence from the central domain around \(x=0\). (2) The applied truncation length \(X=5\) appears as small and quite restrictive for eigenfunctions that do not decay in space. (3) The options of Dirichlet- or Neumann-type thermal conditions in Fig. 4 cannot be expected to hit the lowest eigenfunctions. There is only one exception in Fig. 4, where the Rayleigh number is \(R=9.9182\), reasonably close to the exact value \(R_c = \pi ^2 = 9.8696\) for \(L\le 0\) with infinite horizontal extent. In this case, the truncation walls \(x=\pm \,X\) obey a Neumann condition, producing thermal cells of almost normal-mode shape with wavelength 4, even with an abrupt sign change through \(x=0\). The other case where the truncation walls obey a Dirichlet condition produces a broader and clearly non-separable thermal cell around \(x=0\).

4.1 The 2D Penetrative Eigenfunctions for \(|x| > |L|/2\)

The marginally stable 2D eigenfunctions (with \(k=0\)) will decay spatially into the domains of lateral penetration \(\vert x \vert > |L|/2\) (when \(L > 0\)). In our asymptotic analysis of this decay, we introduce a complex wavenumber \(\alpha (L) = \alpha _r + i \alpha _i\). It represents these locally stable domains intruded by penetrative flow from the marginally stable domain \(\vert x \vert < |L|/2\). A known Rayleigh number \(R=R_c\) is assumed for the marginal stability, where \(R_c < \pi ^2\) when \(L>0\). We know the form of the local 2D eigenfunction, when \(x \gg L/2\)

$$\begin{aligned} \theta = {\hbox {e}}^{i \alpha x} \sin \frac{\pi }{2} \left( z - \frac{1}{2}\right) . \end{aligned}$$
(45)

Inserting \(k=0\) in Eqs. (35)–(36) and eliminating \(\psi\) give the governing equation

$$\begin{aligned} \nabla ^4\theta + R \frac{\partial ^2 \theta }{\partial x^2} = 0, \end{aligned}$$
(46)

where the penetrative eigenfunction (45) produces the relationship

$$\begin{aligned} \left( \alpha ^2 + \frac{\pi ^2}{4}\right) ^2 - R \alpha ^2 = 0, \end{aligned}$$
(47)

which may be called the spatial dispersion relation of penetrative convection. The eigenvalue is the complex wavenumber \(\alpha = \alpha _r + i \alpha _i\), while R has the role of a fixed parameter. This spatial dispersion relation has two solution branches

$$\begin{aligned} \alpha ^2 + \frac{\pi ^2}{4} = \pm \, \sqrt{R} \alpha . \end{aligned}$$
(48)

The minus sign leads to the two eigenvalues relevant for \(x>L/2\)

$$\begin{aligned} \alpha _{1,2} = -\frac{\sqrt{R}}{2} \pm i \frac{\sqrt{\pi ^2 - R}}{2}. \end{aligned}$$
(49)

The penetrative temperature field is

$$\begin{aligned} \theta = A {\hbox {e}}^{i x \sqrt{R}/2} {\hbox {e}}^{-x\sqrt{\pi ^2 - R}/2} \sin \frac{\pi }{2} \left( z - \frac{1}{2}\right) , \end{aligned}$$
(50)

asymptotically valid as \(x \gg L/2\), and the real part represents physical temperature. The amplitude A is complex to account for a free phase angle. By Eq. (35) the streamfunction gets a similar form

$$\begin{aligned} \psi = B {\hbox {e}}^{i x \sqrt{R}/2} {\hbox {e}}^{-x\sqrt{\pi ^2 - R}/2} \sin \frac{\pi }{2} \left( z - \frac{1}{2}\right) , \end{aligned}$$
(51)

with a complex amplitude B. The vorticity Eq. (35) reduces to

$$\begin{aligned} \nabla ^2 \psi = -\, R~ \frac{\partial \theta }{\partial x}. \end{aligned}$$
(52)

with \(k=0\). Equations (50)–(51) combine with the spatial dispersion relation (47) to form the relationship between the penetrative eigenfunctions

$$\begin{aligned} \psi = - \,i \sqrt{R} \theta , \end{aligned}$$
(53)

in the asymptotic limit \(x \gg L/2\). We note the phase shift of a quarter of a wavelength between \(\psi\) and \(\theta\), well known from the HRL problem (Horton and Rogers 1945; Lapwood 1948).

The spatial decay rate is \(\sqrt{\pi ^2 - R}/2 = \sqrt{R_{{\mathrm{local}}} - R}/2\), where \(R_{{\mathrm{local}}}\) is the local Rayleigh critical number in the stable regions of penetration, while R is the smaller global Rayleigh number of externally triggered marginal stability. The spatial decay rate is weaker, the closer the globally marginal Rayleigh number is to the critical local Rayleigh number in the penetrated domain.

The penetrative eigenfunctions (50)–(51) have a dimensionless wavelength \(\lambda\) with asymptotic value \(\lambda = 4 \pi /\sqrt{R}\), in agreement with the wavelength 4 that would emerge if there had been local convection instead of penetrative convection. The wavelength starts at \(\lambda = 4\) with \(L=0\), where the critical Rayleigh number has its highest possible value \(R=\pi ^2\). Note that there is no physical wavelength in the domain \(\vert x \vert < L/2\) where the convection is triggered since the locally preferred mode is a uniform upwelling (or downwelling) flow, known from Nield (1968) as a limit case of zero Rayleigh number with zero wavenumber.

This limit \(R \rightarrow 0\) represents \(L \rightarrow \infty\) in the present model. The spatial decay of the penetrative eigenfunctions is stronger the wider the central domain L. The formulas (50)–(51) give the maximal decay factor \(\exp (- \,\pi \vert x \vert )\), in the limit of large L, where R is close to zero.

The asymptotic factor of decay for a marginal onset mode over its penetrative wavelength \(\lambda\) is given by \(\sqrt{\pi ^2 - R} \lambda /2 = 2 \pi \sqrt{R_{{\mathrm{local}}}/R -1}\), where the local critical Rayleigh number \(R_{{\mathrm{local}}}\) is \(\pi ^2\) for the present choice of boundary conditions.

In Tables 1 and 2, we check the validity of this asymptotic theory of penetrative convection, by comparison with our numerical computations. We consider only the case \(L=0.5\), with the most unstable mode of 2D convection (\(k=0)\). Admittedly, 3D convection is preferred in this case, but we had to choose a relatively small value of L for convergence in the tables, avoiding too strong spatial decay.

Table 1 shows successive cell wall positions \(x=x_m\) for the penetrative flow cells, with the numerically evaluated thermal eigenfunction at each flow cell wall. The decay factor over half a wavelength \(\lambda _m/2=x_{m+1} - x_m\) is computed and compared with the asymptotic formula. The numerical half-wavelength is compared with its asymptotic value. Table 2 shows successive cell wall positions \(x=x_n\) for the thermal cells, with the numerically evaluated flow eigenfunction at each flow cell wall. The agreement with the asymptotic theory is excellent, apart from the first half-wavelength near the central domain of marginal stability. We observe fluctuations in the numerical eigenfunctions, as their values get very small for large x.

Table 1 Cell wall positions \(x=x_m\) for successive penetrative flow cells, defined by \(\psi (x_m,y,0)=0\)
Full size table
Table 2 Cell wall positions \(x=x_n\) for successive penetrative thermal cells, defined by \(\theta (x_n,y,0)=0\)
Full size table

5 Three-Dimensional Convection

We will take into account the possibility that 3D disturbances can be more unstable than 2D disturbances. We will search for the preferred 3D mode of convection, where the wavenumber \(k=k_c(L)\) is selected to determine the minimum Rayleigh number \(R=R_c(L)\).

Our numerical results are obtained by the commercial software Comsol Multiphysics, where we work with two end walls \(x= \pm \,X\) where we specify boundary conditions either of Dirichlet or Neumann type. In the special case \(L=0\), the eigenfunctions do not decay in the far field, where they have sinusoidal variations with x. The choice of X will then set the phases for the eigenfunctions, which would be undetermined for an unlimited horizontal domain. Once L is greater than zero, the eigenfunctions will decay with increasing \(\vert x \vert\), as we have discussed analytically above, for the 2D case.

In Fig. 5 we investigate marginal stability in 3D with a given wavenumber \(k=\pi /2\) in the y direction. This value is equal to the preferred wavenumber for local convection in the stable domains of penetration. When \(L \gg 1\), we know from Nield (1968) that the preferred eigenfunction in the x, z plane represents a uniform vertical flow, and the 3D eigenfunctions will have the asymptotic form

$$\begin{aligned} \theta = A {\hbox {e}}^{i k y},\quad \psi = B {\hbox {e}}^{i k y},\quad L \rightarrow \infty , \end{aligned}$$
(54)

where A and B are complex amplitudes, taking care of phase shift between these eigenfunctions. Inserting these eigenfunctions into the governing Eqs. (35) and (36) gives the asymptotic onset criterion

$$\begin{aligned} R(k) = k^2,\quad L \rightarrow \infty , \end{aligned}$$
(55)

with the special case \(R(\pi /2) = \pi ^2/4\) shown in Fig. 5, where we have added the preferred mode for two cases with smaller wavenumbers, \(k=\pi /2\) (dashed curve) and \(k=0\) (the 2D case, dotted). There is a clear trend that 2D modes are preferred, with possible exceptions only when \(0<L<1\). These plots confirm the asymptotic limit formula (55). Comparing Figs. 3 and 5 indicates that 2D convection is usually preferred at the expense of 3D convection.

Fig. 5
figure 5

Numerically calculated normalized Rayleigh numbers \(R/\pi ^2\) as a function of L for 3D convection with \(k=\pi /2\). The truncation length is chosen as \(X=10\). The onset modes are numbered from 1 to 10 and classified in increasing order according to their Rayleigh numbers (ranked according to the numerical solution at each displacement length)

Full size image

We will search for possible exceptional cases where a 3D onset mode with \(k \ne 0\) may be preferred. Figure 6 is calculated for that purpose, with the small displacement length \(L=0.5\) giving a narrow porous domain with Neumann conditions above and below. The wavenumber k in the perpendicular y direction is given increasing values, where we observe how the preferred cell structures in the x, z broaden with increasing k. Figure 6 shows that the 2D case is not preferred, as there is a slightly smaller Rayleigh number for \(k=0.5\). Figure 6 is calculated with a truncation length \(X=5\), and we note how the highest displayed value \(k=4\) gives a misleading solution with perpendicular isotherms at the truncation boundaries \(x=\pm \, X\) where a thermal Neumann condition is applied.

Figure 7 shows the 3D dependency of the Rayleigh number for the case \(L=0\), with no locally unstable domain and no penetrative convection. The convection cells are spatially periodic in the far field, and the critical Rayleigh number is \(R_c=\pi ^2\), dictated by the far field with its critical wavenumber \(\pi /2\). The disturbance with \(k=\pi /2\) in the y direction seems to represent the global minimum for the Rayleigh number in Fig. 7, yet this is an artifact imposed by the truncation boundaries \(x=\pm \,X\) where the wall conditions set a phase for each eigenfunction which will not fit exactly with the global minimum when \(k=0\). We present two versions of Fig. 7 to show how these phase effects influence the onset problem when the Dirichlet/Neumann conditions for the two eigenfunctions are interchanged.

A finite domain \(L=0.5\) with local instability is shown in Fig. 8. This case is more interesting because it may expose genuine 3D effects for the onset of penetrative convection, but the disadvantage is that we do not have any analytical methods for proving that the small 3D effects that we find are genuine. Figure 8 shows Rayleigh numbers \(R(k)/\pi ^2\) for the onset of convection, where we calculate a critical wavenumber \(k_c = 0.58324\) with truncation length \(X=5\). Its corresponding Rayleigh number is \(R(k_c)/\pi ^2 = 0.80287\), which is about 0.25 per cent below the value \(R(0)/\pi ^2 = 0.80491\) according to 2D theory. Due to the short truncation length, these results are uncertain, so we have recalculated them with a ten times higher truncation length \(X=50\), where we found \(R(k_c)/\pi ^2 = 0.800795\) and \(R(0)/\pi ^2 = 0.802641\). This gives the same trend of slight preference for a 3D onset mode of convection. Figure 8 shows only one mode much more unstable than the higher modes, and all of the higher modes repeat the preference of a 3D convection mode where \(k=\pi /2\) and the Rayleigh number is \(R=\pi ^2\).

Fig. 6
figure 6

The lowest 3D eigenfunctions (with \(L=0.5\)) for different values of k. The chosen thermal boundary condition at the numerical truncation wall is of Neumann type. The truncation length is chosen as \(X=5\)

Full size image
Fig. 7
figure 7

Rayleigh number \((R/\pi ^2)\) as a function of \(k^2\) for 3D convection with \(L=0\). The ten lowest eigenfunctions are represented. In the upper plot, a thermal Neumann condition is applied at the truncation walls \(x=\pm \,X\), where \(X=10\). The lower plot applies a thermal Dirichlet condition. The dots represent the global minimum \(R=\pi ^2\) occurring for 3D convection with \(k=\pi /2\)

Full size image
Fig. 8
figure 8

Rayleigh number \((R/\pi ^2)\) as a function of \(k^2\) for 3D convection with \(L=0.5\). The ten lowest eigenfunctions are represented. Thermal Neumann condition is applied at the truncation walls \(x=\pm \,X\), where \(X=10\)

Full size image

6 Summary and Conclusions

In this paper, a simple physical model is introduced for the convection onset in a homogeneous and isotropic porous layer of infinite horizontal extent, heated uniformly from below. We have investigated the effects of switching the thermomechanical boundary conditions from Dirichlet to Neumann type at one location at the upper boundary, and oppositely at one location at the lower boundary, generating a transition zone with displacement length L between the two locations of switching boundary conditions. When L is positive, a domain of local convection emerges within this transition zone where the Neumann boundary condition allows an almost uniform vertical flow, with laterally penetrative convection in the more stable neighboring domains. The possible types of local convection have exact analytical solutions, with critical Rayleigh numbers either zero, \(\pi ^2\) or \(4 \pi ^2\).

The critical Rayleigh number as a function of L has been determined numerically, obtained with a finite numerical truncation width that is questionable when L is very small or zero. The general solution is 3D and with a wavenumber k in the y direction along the borderlines where the boundary conditions switch between Dirichlet type and Neumann type. This orthogonal wavenumber k transforms the 3D eigenvalue problem into a 2D problem in the xz-plane. The numerical results show that the marginal mode of convection onset is mostly 2D. A weak preference for 3D convection is found when L is smaller than one, but its critical Rayleigh number is only slightly lower than the Rayleigh number for 2D convection. The case \(L=0\) has no other length scale than the unit depth, with critical Rayleigh number \(R = \pi ^2\) for unlimited horizontal extent, known from Nield (1968).

We have investigated analytically the spatial dependency of the penetrative thermomechanical eigenfunctions, which decay exponentially as they intrude the locally stable domains. The dispersion relation for penetrative convection determines a complex wavenumber in such stable domains. The real part of the complex wavenumber settles the periodic cell structure, and these asymptotic predictions agree well with the relevant numerical results. The imaginary part of the complex wavenumber represents spatial decay of the intruding eigenfunctions, and it has also been confirmed by our numerical results for 2D convection.

The present theory for the penetrative onset of convection can be extended to the other sets of boundary conditions studied by Nield (1968), but the spatial dispersion relations must be evaluated numerically. A result for penetrative convection is the asymptotic spatial decay rate \(\sqrt{R_{{\mathrm{local}}} - R}/2\), according to Eq. (50), and we expect similar behavior for more complicated eigenfunctions. The penetrative eigenfunctions will generally be separable in the far field, while they are non-separable in the near field. Our numerical results indicate that weak local singularities exist at the boundary points where the boundary conditions switch between the Dirichlet type and the Neumann type. Such singularities seem not to influence the critical Rayleigh number.

The previous work on penetrative convection onset is mainly concerned with vertical penetration from unstable domains into stably stratified domains. Our model deals exclusively with horizontal sideways (lateral) penetration from marginally stable domains into locally stable domains. The phenomena of penetrative convection are rich, so we chose to design a theoretical model with as few physical parameters as possible. The model is simple physically, but the eigenvalue problem is complicated with no analytical solutions available since the separation of variables cannot be used. From Nield (1968) we know the simple limit case of local convection with zero Rayleigh number and zero wavenumber at marginal stability. Higher modes with finite wavenumber and finite Rayleigh number exist for local convection (Nield 1968), but our numerical results show that they will never be triggered. Thus, the finite Rayleigh number for finite gap width of the locally unstable domain will always be set by the neighboring domains of penetration. It also means that the spatial phase angle in the surrounding penetrative convection cells is indirectly set by these local domains themselves and not exported from the locally unstable domain which has uniform upwelling (or downwelling) with no sign change.

A physically simplifying element in the present theory is that the tendency of uniform upwelling/downwelling in the unstable zone prevents phase effects from originating there. Phase effects from recirculating cells in the unstable zone will give a more complicated interaction with the surroundings. We have not discussed the vortex structures with short length scale that appear in the numerical solutions near the points where the boundary conditions switch between Neumann type and Dirichlet type, since the solution is probably singular there, and the numerical results will then only represent the outer solution in a matched asymptotic expansion. Finite-amplitude convection with lateral penetration is a challenging topic for follow-up research.