Laterally Penetrative Onset of Convection in a Horizontal Porous Layer

The onset of Darcy–Bénard convection in an unlimited horizontal porous layer is studied theoretically. The thermomechanical boundary conditions of Dirichlet or Neumann type at the lower and upper plane are switched from one type to another, at certain values of the horizontal x-coordinate. A semi-infinite portion of the lower boundary is defined as thermally conducting and impermeable, while the remaining boundary is open and with given heat flux. At the upper boundary, the same thermomechanical conditions are applied, but with a relative spatial displacement L and in the opposite spatial order. A domain of local destabilization around the origin is generated between the lines of discontinuity x=±L/2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x = \pm \,L/2$$\end{document}. The marginal state of convection is triggered centrally, while it is penetrative in the domains exterior to the central domain. The onset problem is solved numerically, with a general 3D mode of disturbance, but 2D disturbances are preferred in most cases. The critical Rayleigh number is given as a function of the dimensionless gap width L and the wavenumber k in the y direction along the lines of discontinuity in the boundary conditions. An asymptotic formula for 2D penetrative eigenfunctions is shown to be in agreement with the numerical results.


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.

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 |x| → ∞.
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 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 . The porous medium is homogeneous and isotropic, with permeability K. The standard Darcy-Boussinesq equations for free thermal convection can be written In these equations, p is the dynamic pressure, is the coefficient of thermal expansion, = 0 is the fluid density at the reference temperature T 0 , is the dynamic viscosity of the saturating fluid, c p is the specific heat at constant pressure, and 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 heatconducting and impermeable (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) This dynamic constant-pressure condition implies that × ∇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. Δ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 The semi-infinite right-hand portion of the displaced upper boundary is subject to a given heat flux, and with constant pressure (open boundary) 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 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 − < z + h∕2 < 0 below the porous medium and 0 < z − h∕2 < 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 = − , x < 0 and z − h∕2 = , 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.

Dimensionless Equations
From now on, we work with dimensionless variables. We reformulate the mathematical problem in dimensionless form by means of the transformations where m = m ∕( 0 c p ) f is the thermal diffusivity of the saturated porous medium. We denote the vertical unit vector by , directed upward. Here, the Rayleigh number R is defined as

Basic Solution
The stationary basic solution of Eqs. (8)-(14) is given subscript "b" This basic state of hydrostatic fluid has a linear temperature profile.

Linearized Perturbation Equations
In our stability analysis, we disturb the basic state (16) with perturbed fields where the perturbations , , p ′ are functions of x, y, z and t. Linearizing Eqs. (8)-(10) with respect to perturbations and eliminating the pressure gives With vanishing vertical vorticity, one single scalar function (x, y, z) represents the 3D thermomechanical vector field. The incompressible flow field is given by this poloidal vector potential as = ∇ × (∇ × ) , where is the vertical unit vector. The velocity components are where the horizontal Laplacian operator ∇ 2 1 = 2 ∕ x 2 + 2 ∕ y 2 has been introduced. From Eq. (18), it follows that the perturbation temperature is given by Assuming a non-oscillatory marginal state, the heat Eq. (19) becomes The boundary conditions at the lower boundary are The conditions at the upper boundary are The coupled eigenvalue problem for the dimensionless eigenfunctions and 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 and 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.

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 where i denotes the imaginary unit. The 2D perturbation equations for the eigenfunctions (x, z) and (x, z) are We will equip our 2D eigenvalue problem in the vertical xz-plane with boundary conditions. The conditions at the lower boundary are The conditions at the upper boundary are

The 2D Eigenvalue Problem Where We Introduce a Streamfunction
We will work with a modified eigenfunction (x, z) , which is the streamfunction in the x, z plane when k = 0 . It is defined by The perturbation equations for the eigenfunctions (x, z) and (x, z) are A 2D eigenvalue problem is formulated in the vertical xz-plane, to be equipped with boundary conditions. The conditions at the lower boundary are The conditions at the upper boundary are 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 . The lower sketch shows the corresponding differential equation and boundary conditions for . 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. The 3D temperature field is with the 3D velocity components 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.
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 = ± X , since the numerical method does not handle the assumed infinite extent ( X → ∞ ). We have chosen X = 5 in Fig. 4, where the complete computational domain − X < x < X is displayed.
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 nonseparable 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 ≤ 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 = ± X will then overrule the eigenfunctions by setting their phase in the far field, and the critical Rayleigh number will be exactly R c = 2 for all values L ≤ 0 when the horizontal porous layer extends to infinity. We include two different sets of boundary conditions for the plots with L ≤ 0 in Fig. 4, denoting the applied thermal condition in each case. When we apply a Dirichlet condition for at the truncation walls x = ± X , the conditions for 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 2 when L ≤ 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 = 2 = 9.8696 for L ≤ 0 with infinite horizontal extent. In this case, the truncation walls x = ± 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.

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 |x| > |L|∕2 (when L > 0 ). In our asymptotic analysis of this decay, we introduce a complex wavenumber (L) = r + i i . It represents these locally stable domains intruded by penetrative flow from the marginally stable domain |x| < |L|∕2 . A known Rayleigh number R = R c is assumed for the marginal stability, where R c < 2 when L > 0 . We know the form of the local 2D eigenfunction, when x ≫ L∕2 Inserting k = 0 in Eqs. (35)-(36) and eliminating give the governing equation where the penetrative eigenfunction (45) produces the relationship which may be called the spatial dispersion relation of penetrative convection. The eigenvalue is the complex wavenumber = r + i i , while R has the role of a fixed parameter. This spatial dispersion relation has two solution branches The minus sign leads to the two eigenvalues relevant for x > L∕2 The penetrative temperature field is asymptotically valid as x ≫ 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 with a complex amplitude B. The vorticity Eq. (35) reduces to with k = 0 . Equations (50)-(51) combine with the spatial dispersion relation (47) to form the relationship between the penetrative eigenfunctions in the asymptotic limit x ≫ L∕2 . We note the phase shift of a quarter of a wavelength between and , well known from the HRL problem (Horton and Rogers 1945;Lapwood 1948).
The spatial decay rate is where R 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 with asymptotic value = 4 ∕ √ R , in agreement with the wavelength 4 that would emerge if there had been local convection instead of penetrative convection. The wavelength starts at = 4 with L = 0 , where the critical Rayleigh number has its highest possible value R = 2 . Note that there is no physical wavelength in the domain |x| < 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 → 0 represents L → ∞ 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(− |x|) , 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 is given by where the local critical Rayleigh number R local is 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 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.  Numerical values of m = (x m , y, 0) are given, with their decay factors | m ∕ m+1 | and half-wavelengths x m+1 − x m . The asymptotic formulas for spatial decay and wavelength are evaluated. The chosen displacement length is L = 0.5 , and we consider the preferred mode of 2D convection ( k = 0 ), with the computational domain − X < x < X chosen as X = 20 . The Rayleigh number of marginal stability is R = 7.904 1 3

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 = ± 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 |x| , as we have discussed analytically above, for the 2D case.
In Fig. 5 we investigate marginal stability in 3D with a given wavenumber k = ∕2 in the y direction. This value is equal to the preferred wavenumber for local convection in the stable domains of penetration. When L ≫ 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 (54) = Ae iky , = Be iky , L → ∞, Table 2 Cell wall positions x = x n for successive penetrative thermal cells, defined by (x n , y, 0) = 0 Numerical values of m = (x m , y, 0) are given, with their decay factors | n ∕ n+1 | and half-wavelengths x n+1 − x n . The asymptotic formulas for spatial decay and wavelength are evaluated. The chosen displacement length is L = 0.5 , and we consider the preferred mode of 2D convection ( k = 0 ), with the computational domain − X < x < X chosen as X = 20 . The Rayleigh number of marginal stability is R = 7.904 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 with the special case R( ∕2) = 2 ∕4 shown in Fig. 5, where we have added the preferred mode for two cases with smaller wavenumbers, k = ∕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. We will search for possible exceptional cases where a 3D onset mode with k ≠ 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 = ± 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 = 2 , dictated by the far field with its critical wavenumber ∕2 . The disturbance with k = ∕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 = ± 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)∕ 2 for the onset of convection, where we calculate a critical wavenumber k c = 0.58324 with Fig. 7 Rayleigh number (R∕ 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 = ± X , where X = 10 . The lower plot applies a thermal Dirichlet condition. The dots represent the global minimum R = 2 occurring for 3D convection with k = ∕2   Fig. 8 Rayleigh number (R∕ 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 = ± X , where X = 10 truncation length X = 5 . Its corresponding Rayleigh number is R(k c )∕ 2 = 0.80287 , which is about 0.25 per cent below the value R(0)∕ 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 )∕ 2 = 0.800795 and R(0)∕ 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 = ∕2 and the Rayleigh number is R = 2 .

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, 2 or 4 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 = 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 √ R 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.