Afiqah Nadzirah

• D Andrew

S Rees
email: rees@bath.ac.uk

Abdelkader Mojtabi

B D Andrew

The Effect of Conducting Boundaries on the Onset of Convection in a Porous Layer Which is Heated from Below by Internal Heating

Keywords: Natural convection, Porous medium, Linear stability, Dispersion relation, Conducting bounding layers

The onset of convection in a porous layer heated from below is considered, and we determine how the presence of two solid but heat-conducting bounding plates of finite thickness alters the manner in which convection ensues. Heat is generated by the lower plate (with an insulating lower boundary), but the upper one is passive with a fixed upper boundary temperature. It is shown that this composite layer may mimic in turn one of the three different types of classical single-layer onset problems which are well-known in the literature. The type which is selected (or indeed whether it corresponds to a transitional case) depends quite critically on the precise values of the relative thickness of the solid layers and their conductivity ratio. It is also shown that care needs to be taken over declaring that the solid plates are thin: extreme values of the conductivity ratio can yield a stability criterion which appears to be different from that suggested by the imposed boundary conditions.

Introduction

The modelling of convective instability in porous layers now has a very strong heritage. The foundational works are by Lapwood (1948) and [START_REF] Horton | Convection currents in porous media[END_REF], and these determined the criterion for the onset of convective instability in a horizontal porous layer which is heated from below. In addition, the upper and lower boundaries were considered to be impermeable and subject to constant fixed temperatures, and the porous medium was isotropic and homogeneous. Other types of boundary condition, such as a constant heat flux or constant pressure, may be applied in a variety of different ways; Table 6.1 in [START_REF] Nield | Convection in Porous Media[END_REF] gives a comprehensive account of the critical values of the Rayleigh number and wavenumber for these different combinations. The effects of confinement, such as convection in a cuboidal cavity and in a cylindrical cavity (with a circular horizontal planform), were considered by [START_REF] Beck | Convection in a box of porous material saturated with fluid[END_REF] and [START_REF] Bau | Onset of convection in a permeable medium between vertical coaxial cylinders[END_REF], respectively. Unlike the corresponding pure fluid configurations, there exist aspect ratios of the cuboid and radii of the cylinder for which the critical Darcy-Rayleigh number is identical with that of the infinite unconfined layer.

The literature is now replete with the study of how the onset and subsequent development of convection is altered by the presence of different effects, such as Brinkman effects, non-Boussinesq effects, the presence of imperfections such as sinusoidally varying boundary temperatures, throughflow effects, rotation, layering, anisotropy, internal heating, local thermal nonequilibrium, more than one diffusing species and unsteady basic states; note that this list, while long, is not comprehensive. Although some work has been reported on configurations where two porous sublayers are separated by a solid intermediate sublayer (see [START_REF] Postelnicu | Thermal stability of two fluid porous layers separated by a thermal barrier[END_REF][START_REF] Jang | Thermal instability of two horizontal porous layers with a conductive partition[END_REF][START_REF] Patil | The onset of convection in a porous layer with multiple horizontal solid partitions[END_REF][START_REF] Straughan | Nonlinear stability of convection in a porous layer with solid partitions[END_REF] for systems of this general type), we are presently interested three-layer systems where a porous sublayer is bounded by two heat-conducting but impermeable solid sublayers. [START_REF] Mojtabi | The onset the convection in Horton-Rogers-Lapwood experiments: of conducting bounding plates[END_REF] considered such a three-layer system subject to constant heat flux boundary conditions, and this was motivated by the need to model an experimental system where solid boundaries do not have zero thickness and idealised boundary conditions cannot be applied in practice. It was found that the critical parameters (namely the Darcy-Rayleigh number and wavenumber) and the shape of the neutral curve depend quite strongly on the thickness of the solid boundaries and their conductivity relative to that of the porous sublayer. For example, when the solid sublayers are sufficiently highly conducting, the full system may be made to resemble a single porous layer system with constant temperature boundary conditions. In addition, there is a locus in parameter space where the neutral curve has a quartic minimum at zero wavenumber, and this locus separates the regime within which the more usual parabolic minimum at nonzero wavenumbers exists from that where the parabolic minimum is at zero wavenumber. A second study was undertaken by [START_REF] Rees | The effect of conducting boundaries on weakly nonlinear Darcy-Bénard convection[END_REF] where constant temperature boundary conditions were applied. There are many qualitative similarities between this study and the previous one by the same authors except that there is no quartic minimum, and the zero critical wavenumber limit is obtained as the solid conductivities tend to zero. [START_REF] Rees | The effect of conducting boundaries on weakly nonlinear Darcy-Bénard convection[END_REF] also provided a weakly nonlinear analysis to determine regimes within which postcritical convection takes the form of rolls or patterns with square planform; this was an extension of the work of [START_REF] Riahi | Nonlinear convection in a porous layer with finite conducting boundaries[END_REF] who considered solid layers of infinite thickness. A paper by [START_REF] Rees | The effect of conducting boundaries on Lapwood-Prats convection[END_REF] also considered the effect on the onset of convection of a horizontal pressure gradient. If was found that there is a thermal drag caused by the stationary solid boundaries which reduces the velocity of the cells along the layer.

In the present paper, we consider the linear stability of a configuration which has some similarities to those of [START_REF] Mojtabi | The onset the convection in Horton-Rogers-Lapwood experiments: of conducting bounding plates[END_REF] and [START_REF] Rees | The effect of conducting boundaries on weakly nonlinear Darcy-Bénard convection[END_REF]. Here, the lower solid sublayer is insulated on its lower surface, but its interior is subject to a uniform rate of heat generation, while the upper sublayer is subject to a constant temperature. It is shown in the Appendix that this may be mapped onto a system where the lower solid sublayer is subject to a constant heat flux but without heat generation within the sublayer. While this may seem to be a simple intermediate case between those of [START_REF] Mojtabi | The onset the convection in Horton-Rogers-Lapwood experiments: of conducting bounding plates[END_REF] and [START_REF] Rees | The effect of conducting boundaries on weakly nonlinear Darcy-Bénard convection[END_REF], three different limiting cases for the critical parameters are found in different asymptotic regions of parameter space, and this is a novel finding.

Porous layer

Solid layer

Solid layer

ẑ = -h ẑ = 0 ẑ = H ẑ = H + h ∂T 1 ∂ ẑ = 0 T = Tupper
Fig. 1 Definition sketch of the configuration being studied. The porous layer is sandwiched between two impermeable but thermally conducting layers

Governing Equations and Basic Solution

We consider a three-layer system as shown in Fig. 1. It consists of a central uniform porous layer (layer 2) within which convective flow may occur, and this is embedded between two thermally conducting but impermeable layers. These solid layers have identical properties, namely thickness, conductivity and heat capacity. A uniformly cold temperature is imposed on the upper surface of the upper solid layer (layer 3), while the lower layer (layer 1) is subject to uniform Joule heating but its lower surface is insulated. Thus, all the heat is generated within the lowest sublayer and passes out through the upper surface, and a potentially unstable temperature gradient is formed. The origin of the coordinate system is located at the bottom of the porous layer with x being the horizontal coordinate and ẑ the vertical coordinate.

The full set of governing equations are as follows. For layer 1 (-h ≤ ẑ ≤ 0):

(ρC) s ∂ T 1 ∂ t = k s ∂ 2 T 1 ∂ x2 + ∂ 2 T 1 ∂ ẑ2 + q;
(1) for the porous layer (0 ≤ z ≤ H ):

∂ û ∂ x + ∂ ŵ ∂ ẑ = 0, (2) û = - K μ ∂ p ∂ x , ŵ = - K μ ∂ p ∂ ẑ + ρ f ḡβ K μ (T 2 -T ref ), ( 3
) (ρC) pm ∂ T 2 ∂ t + (ρC) f u ∂ T 2 ∂ x + w ∂ T 2 ∂ ẑ = k pm ∂ 2 T 2 ∂ x2 + ∂ 2 T 2 ∂ ẑ2 ; (4)
and for the upper layer (H ≤ z ≤ H + h):

(ρC) s ∂ T 3 ∂ t = k s ∂ 2 T 3 ∂ x2 + ∂ 2 T 3 ∂ ẑ2 . ( 5
)
In these equations, the numerical subscripts refer to the appropriate sublayer, and T ref is the basic temperature of the interface between the porous sublayer and the lower solid sublayer for the purely conducting situation. The boundary and interface conditions are given by, ẑ = -h :

∂ T 1 ∂z = 0, ẑ = 0 : ŵ = 0, T 1 = T 2 , k s ∂ T 1 ∂ ẑ = k pm ∂ T 2 ∂ ẑ , ẑ = H : ŵ = 0, T 2 = T 3 , k pm ∂ T 2 ∂ ẑ = k s ∂ T 3 ∂ ẑ , ẑ = H + h : T 3 = T upper . (6)
The basic state consists of no flow ( û = ŵ = 0) and a steady ẑ-dependent temperature. It is a straightforward if tedious matter to determine the conduction profiles in each layer:

T 1 = T upper + q 3h 2 2k s + H h k pm - q k s (ẑ + h) 2 2 , T 2 = T upper + q h 2 k s - qh k pm (ẑ -H ), T 3 = T upper - qh k s (ẑ -H -h). ( 7
)
The temperature drop across the porous layer will be used as the temperature scale for the purposes of nondimensionalisation. This value is

T = q Hh k pm , (8)
and we see that the temperature, T ref , of the lower interface temperature is given by,

T ref = T upper + q h 2 k s + H h k pm . ( 9
)
We may now nondimensionalise by introducing the following scalings,

( x, ẑ) = H (x, z), t = H 2 k pm /(ρC) pm t, (û, ŵ) = k pm (ρC) f H (u, w), p = μk pm K (ρC) f p, T = T ref + T θ. (10)
Thus, Eqs. ( 2) and (3) become,

∂u ∂ x + ∂w ∂z = 0, ( 11
) u = - ∂ p ∂ x , w = - ∂ p ∂z + Ra θ, (12)
Here, the Darcy-Rayleigh number, Ra, is defined as,

Ra = ρ f ḡβ K H T μ(k pm /(ρC) f )) , ( 13
)
where T is given in terms of q in Eq. ( 8). The three heat transport equations become,

(ρC) s k s k pm (ρC) pm ∂θ 1 ∂t = ∂ 2 θ 1 ∂ x 2 + ∂ 2 θ 1 ∂z 2 + 1 dδ , ( 14
)
∂θ 2 ∂t + u ∂θ 2 ∂ x + w ∂θ 2 ∂z = ∂ 2 θ 2 ∂ x 2 + ∂ 2 θ 2 ∂z 2 , ( 15
) (ρC) s k s k pm (ρC) pm ∂θ 3 ∂t = ∂ 2 θ 3 ∂ x 2 + ∂ 2 θ 3 ∂z 2 , ( 16
)
where

d = k s k pm and δ = h H (17)
are the respective conductivity and thickness ratios between the solid and porous sublayers.

The boundary and interface conditions given in Eq. ( 6) become,

z = -δ : ∂θ 1 ∂z = 0, z = 0 : ŵ = 0, θ 1 = θ 2 , d ∂θ 1 ∂z = ∂θ 2 ∂z , z = 1 : ŵ = 0, θ 2 = θ 3 , ∂θ 2 ∂z = d ∂θ 3 ∂z , z = 1 + δ : θ 3 = -1 - δ d . ( 18
)
The apparently unusual boundary condition for θ 3 at z = 1 + δ has been defined in such a way that the temperature field varies between 0 and -1 within the porous layer. Thus, the basic solution is given by

θ = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ - z(z + 2δ) 2dδ (-δ ≤ z ≤ 0), -z (0 ≤ z ≤ 1), -1 - z -1 d (1 ≤ z ≤ 1 + δ). ( 19
)
While there is a unit change in temperature across the porous layer, the overall temperature drop over the whole system is 1 +3δ/2d, and this increases as either the thickness of the solid sublayers increase or their conductivity decreases. Five representative profiles of the basic temperature field are given in Fig. 2 for the case δ = 0.5; profiles for other values of δ display the same qualititative features. The insulating condition at z = -1 -δ is evident because of the zero gradient. By design, the temperature profile in the porous layer is independent of the values of both d and δ.

Figure 2 shows that there is little variation in the temperature within the solid sublayers when d is large. This happens quite naturally because of the need to conserve the heat flux across the interface. And then, when d takes small values, there is a very large variation in temperature within the solid sublayers.

Perturbation Analysis

The basic state consists of the temperature profiles given in Eq. ( 19) with no flow. Since Squire's theorem may be shown easily to apply, i.e. that all three-dimensional disturbances may be decomposed into a sum or integral of two-dimensional ones, we need to consider only two-dimensional disturbances. Therefore, we adopt a streamfunction/temperature formulation using, u = -ψ z and w = ψ x , which guarantees that the equation of continuity, Eq. ( 11), is satisfied. Although the heat transport equations in the two solid sublayers remain unchanged, Eqs. ( 12) and ( 15) take the following forms,

∂ 2 ψ ∂ x 2 + ∂ 2 ψ ∂z 2 = Ra ∂θ ∂ x , ( 20
)
∂θ 2 ∂t + ∂ψ ∂ x ∂θ 2 ∂z - ∂ψ ∂z ∂θ 2 ∂ x = ∂ 2 θ 2 ∂ x 2 + ∂ 2 θ 2 ∂z 2 . ( 21
)
We may now introduce a disturbance of infinitesimally small amplitude using,

ψ = ψ basic + , θ j = θ j, basic + j ( j = 1, 2, 3), ( 22
)
where ψ basic = 0 and the basic temperature fields are given in Eq. ( 19). On neglecting powers of the disturbances, the linearised stability equations take the form,

(ρC) s k s k pm (ρC) pm ∂ 1 ∂t = ∂ 2 1 ∂ x 2 + ∂ 2 1 ∂z 2 , ( 23
)
∂ 2 ∂ x 2 + ∂ 2 ∂z 2 = Ra ∂ 2 ∂ x , ( 24
)
∂ 2 ∂t = ∂ 2 2 ∂ x 2 + ∂ 2 2 ∂z 2 + ∂ ∂ x , ( 25
) (ρC) s k s k pm (ρC) pm ∂ 3 ∂t = ∂ 2 3 ∂ x 2 + ∂ 2 3 ∂z 2 . ( 26
)
These equations are to be solved subject to the homogeneous form of the boundary conditions given in (18), with the addition of = 0 on both z = 0 and z = 1. ( 27)

It may also be shown easily that the principle of exchange of stabilities applies, and therefore, the onset of convection is stationary and the time derivatives in Eqs. ( 23)-( 26) may be neglected. If the domain is of infinite horizontal extent, then we may Fourier-decompose disturbances into independent monochromatic cells. Thus, we set,

(x, z, t) → i (z)e ikx + c.c., j (x, z, t) → j (z)e ikx + c.c., ( 28
)
where k is the horizontal wavenumber. Equations ( 23)-( 26) reduce to

1 -k 2 1 = 0, ( 29
) -k 2 = Ra k 2 , ( 30
) 2 -k 2 2 = k 2 , ( 31
) 3 -k 2 3 = 0, ( 32
)
where primes denote derivatives with respect to z.

The equations for 1 and 3 may be solved analytically and converted into suitable boundary conditions for 2 upon application of the interface conditions, as follows. For Layer 1, Eq. ( 29) has the solution,

1 = E cosh(kz) + F sinh(kz). ( 33
)
The application of the boundary condition, 1 (-δ) = 0, means that this solution may now be written in the form,

The linear stability problem has now been reduced to solving Eqs. ( 30) and ( 31) subject to the two thermal boundary conditions given in ( 35) and ( 36), together with those in Eq. ( 27).

Dispersion Relation

Equations ( 30) and ( 31) form a fourth-order homogeneous system, and there are four independent solutions. These may be written in the following form,

= √ Ra A cosh λz + B sinh λz + C cos σ z + D sin σ z , ( 37
)
where A, B, C and D are arbitary constants, and where

λ = k( √ Ra + k), σ = k( √ Ra -k). ( 39
)
The Darcy-Rayleigh number is found always to lie above the curve, Ra = k 2 , and therefore, the value of σ is always purely real.

If we now apply the boundary conditions, ( 27), ( 35) and ( 36), then we obtain the following matrix/vector system for the four arbitrary constants:

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 0 1 0 cosh λ sinh λ cos σ sin σ -2kd tanh(kδ) λ 0 -σ λ sinh λ+ λ cosh λ+ σ sin σ - -λ cos λ- kd coth(kδ) cosh λ kd coth(kδ) sinh λ kd cos σ kd coth(kδ) sin σ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ A B C D ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 40
)
This system has nonzero solutions only when the determinant of the matrix is zero. After some manipulation, we find that this determinantal condition is equivalent to

λσ 1 -cosh λ cos σ -(1 + 2d 2 )k 2 sinh λ sin σ -kd tanh kδ + coth kδ σ sinh λ cos σ + λ cosh λ sin σ = 0. ( 41
)
It is impossible to solve this equation analytically to give Ra in terms of k, d and δ, and therefore, we resort to numerical means. A simple Newton-Raphson may be used easily to determine neutral curves as k varies for fixed values of d and δ. A slightly more complicated version was also written in order to find critical values, i.e. those values of k and Ra at the lowest point of the neutral curve. The manner in which this was undertaken is precisely as described in [START_REF] Rees | The effect of conducting boundaries on weakly nonlinear Darcy-Bénard convection[END_REF] and [START_REF] Rees | Onset of convection in porous layers with multiple horizontal partitions[END_REF] and is therefore omitted here. The codes were written in quadruple precision in Fortran90 because it was found that this amount of precision was required to obtain accurate computations when d and δ take values which are as small as 10 -10 and as large as 10 10 .

Results and Discussion

In what follows there are three sets of critical data which are useful references. For the classical Darcy-Bénard layer, the critical values of the Darcy-Rayleigh number and wavenumber are, CT/CT:

Ra c = 4π 2 , k c = π, ( 42
)
where the appelation, CT/CT, refers to the fact that both the upper and lower boundary conditions correspond to a constant temperature.

Onset Profiles

We begin by showing, in Fig. 3, a selection of onset profiles ( , ) so that the effect of variations in the two nondimensional parameters, d and δ, on those profiles may be seen. In both frames the temperature profile, , takes positive values, while the streamfunction disturbance profile, , takes negative values. The profiles have been scaled so that the maximum magnitude is precisely unity for easy comparison. Given our definition of the streamfunction, the values of turn out conveniently to be negative when is positive. Although we have set k = 2 in all cases, the qualitative trends shown in Fig. 3 do not change for different wavenumbers.

Figure 3a shows the effect of variations in d when δ = 0.4. When d takes large values, the solid sublayers are highly conducting and we see therefore that the disturbance temperature is almost perfectly constant; this situation is close to the CT/CT single-layer case mentioned above. At the opposite extreme, d 1, the solid layers are poor conductors, and this gives rise to relatively large variations in within the solid regions, as was seen for the basic conduction profiles, and only small variations within the porous sublayer. This is approximately the CHF/CHF case mentioned above. It is also interesting to note that the streamfunction profile hardly varies at all over the whole range of values of d.

Figure 3b shows the effect of changing the depth, δ, of the solid sublayers. We have chosen a fairly representative value for the conductivity ratio, d = 5, and this case shows an interesting transition in the profiles from one which is close to being of CT/CT type when δ is large (i.e. thick solid sublayers) to one which is almost indistinguishable from a CT/CHF type when δ is small. In such cases, the very thin solid sublayers appear to be yield stability properties and profiles which are identical to cases where the layers are entirely absent. We will see later that this initial hypothesis turns out not to be that straightforward in general, and attention must be paid to the magnitude of δ in comparison with that of d.

Neutral Curves

Figure 4 shows a small selection of neutral stability curves for some typical values of the parameters, d and δ, where each lies in the range, 10 -3 -10 3 . These curves show the standard unimodal form with a clearly identifiable minimum value; numerical investigations indicate that unimodal curves are always obtained for the present configuration. The uppermost curve (d = 10 3 ) in all four cases is almost exactly identical to that for the classical Darcy-Bénard problem (CT/CT) for which the critical values are Ra c = 4π 2 and k c = π. This is because the solid layers are highly conducting which means that temperature variations are rapidly damped, and therefore the interfaces act as perfectly conducting boundary conditions. That said, the critical values for the d = 10 3 curve when δ = 10 -3 are very slightly smaller than the classical values, but this is difficult to see graphically. On the other hand, the lowest curves (d = 10 -3 ) for all four cases change quite substantially as δ increases, and the critical values decrease quite markedly. Indeed, when δ is large, the highly insulating bounding layers mimic the CHF/CHF porous layer for which Ra c = 12 and k = 0.

On the other hand, from these graphs, it is possible to find out how quickly the neutral curves change as d varies. It is clear that there is almost no change in the location of the neutral curve when d increases from 10 -2 to 10 2 when δ = 10 -3 . But when δ = 10 3 , there is a rapid change as d increases from 10 -1 to just below 10 1 ; outside of these ranges there is little variation. While we have gleaned this information from Fig. 4, it is clear that we need to consider the variations in the critical parameters themselves, rather than by observing further examples of neutral curves.

An extensive set of computations has therefore been undertaken to determine the critical values of Ra and k in the ranges 10 -6 ≤ δ ≤ 10 6 and 10 -10 ≤ d ≤ 10 10 , and these are shown in Fig. 5. Referring to both frames in this figure, we see that the general behaviour of the critical values is that they pass from the CHF/CHF case via the CT/CHF case to the CT/CT case as d increases from exceptionally small values to exceptionally large values. When the (lowest, dotted line), 10 -2 , 10 -1 , 10 0 , 10 1 , 10 2 and 10 3 (uppermost) solid sublayers are comparable to the thickness of the porous layer, or larger, then there is a direct transition from CHF/CHF to CT/CT which is quite quick near d = 10 0 when d is increasing. But when δ < 10 -3 , there is range centred on d = 10 0 where the composite layer resembles the CT/CHF layer with very high accuracy. Figure 6 shows a close-up image of the central region of the wavenumber curves in Fig. 5. This shows that, when δ = 10 -5 , the composite layer mimics a CT/CHF single layer when 10 -2 < δ < 10 2 , and when δ = 10 -6 , this happens when 10 -3 < δ < 10 3 . In the close-up view, we also see that there is an undershoot in the critical wavenumber as d decreases whenever δ is sufficiently small. This presages the passage towards the CT/CHF equivalence. It is unclear what the physical reason for this might be, but a similar behaviour does not occur for the critical Darcy-Rayleigh number.

Summary

An alternative view of the curves shown in Fig. 5 is provided by Fig. 7, which shows contour representations of both Ra c and k c in the ranges, 10 -10 ≤ d, δ ≤ 10 10 . A grid with 100 inter- Fig. 5 Critical values of Ra and k as a function of log 10 d for δ = 10 -6 , 10 -5 , 10 -4 , 10 -3 , 10 -2 , 10 -1 , 10 0 (dotted) and 10 1 vals in both directions was used for these figures. We now see clearly how the (d, δ) plane is divided into three very distinct regions with smooth transitions between these regions. Generally, the right-hand side corresponds to being equivalent to a CT/CT single layer, the left-hand Fig. 6 Close-up view of the critical values of k as a function of log 10 d for 10 -5 , 10 -4 , 10 -3 , 10 -2 , 10 -1 , 10 0 (dotted) and 10 1 side to a CHF/CHF layer, while the triangular wedge is a CT/CHF configuration. Therefore, our composite layer is clearly capable of supporting three different types of behaviour.

Earlier, a statement was made that what constitutes a narrow solid layer is not a straightforward concept. In Fig. 6, if we choose δ = 10 -5 to represent a very thin pair of solid sublayers, then we should expect this to represent the single-layer CT/CHF case simply because 10 -5 is exceptionally small. However, the critical values correspond to the single-layer CT/CHF critical values only when d lies in the range 10 -4 < d < 10 4 . Outside of this range, we obtain the CHF/CHF case when d < 10 -9 and the CT/CT case when d > 10 5 . These extreme values for d overrule our standard intuition of what 'thin' means in the context of the stability. Of course, if one takes the alternative viewpoint that d is fixed, but δ decreases, then one will always eventually obtain a CT/CHF configuration.

It is clear from the two subfigures in Fig. 7 that the transitions between the three main regions take place at triangular region at the bottom of that plot and which are = 10 0 are also within 10 -6 of the CT/CHF value in Eq. ( 44). This corresponds almost perfectly with curves given in Fig. 6.

Conclusions

We have conducted a linear stability analysis to determine the effect on stability criterion of having two solid but identical bounding sublayers above and below a saturated porous layer where the lower sublayer generates heat internally. The equations were nondimensionlised in such a way that the Darcy-Rayleigh number is based on the temperature drop across the porous layer. There are two parameters, namely the ratio between the solid and porous sublayers, d and the relative thickness of the bounding plates, δ. Therefore, we have been able to provide a comprehensive account of the stability properties over a very wide range of values of the parameters. The stability criteria were determined by finding the zeros of a dispersion relation and by minimisation over the wavenumber.

It has been found that the addition of conducting layers above and below the porous layer allows the stability criteria to vary quite substantially and, in particular, to allow this composite layer to mimic one of three standard single-layer configurations, namely the CT/CT, CT/CHF and the CHF/CHF cases. This type of behaviour is quite different from that of [START_REF] Mojtabi | The onset the convection in Horton-Rogers-Lapwood experiments: of conducting bounding plates[END_REF] where the composite layer which they studied only displayed two types of limiting behaviour.

From the practical point of view of constructing an experimental configuration, we have shown that the presence of a physically thin pair of solid bounding plates does not necessarily guarantee that the stability properties correspond to the type (namely CT/CHF) corresponding to the boundary conditions which are applied. Rather, when δ is very small, only a certain range of values of d will yield the CT/CHF onset criterion.

We do not know at this stage whether the rolls which we have considered will form the stable mode of convection at values of the Darcy-Rayleigh number which are just above the critical values. Given the conclusions of [START_REF] Riahi | Nonlinear convection in a porous layer with finite conducting boundaries[END_REF] and [START_REF] Rees | The effect of conducting boundaries on weakly nonlinear Darcy-Bénard convection[END_REF], we feel that it is quite likely that square cells are likely to be preferred when the conductivity of the bounding sublayers is sufficiently small, but this aspect would need to be considered quantitatively elsewhere. However, it is very likely indeed that both the roll and square cell patterns will arise from supercritical bifurcations, and therefore, all the results presented here will also be valid should square cells be the preferred pattern post-onset.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix

The aim of this Appendix is to show that the problem solved in the main body of this paper is equivalent to solving one where the lower solid sublayer does not have heat generation but is subject to a constant heat flux Specifically Eq. ( 1) is replaced by,

(ρC) s ∂ T 1 ∂ t = k s ∂ 2 T 1 ∂ x2 + ∂ T 1 ∂ ẑ2 , ( 46
)
and the ẑ = -h boundary condition, Eq. ( 6) by

-k s ∂ T 1 ∂z = q. ( 47
)
This means that the basic temperature profiles become,

T 1 = T upper + q h k s + H k pm - q k s ẑ, T 2 = T upper + q h k s - q k pm (ẑ -H ), T 3 = T upper - q k s (ẑ -H -h). ( 48
)
The temperature drop across the porous layer is, therefore,

T = q H k pm , ( 49
)
and the interface temperature at ẑ = 0, T ref , is given by,

T ref = T upper + q h k s + H k pm . ( 50
)
Nondimensionalisation now takes place as given in Eq. ( 10) where T is given by Eq. ( 49).

After nondimensionalisation, Eq. ( 14) is replaced by

(ρC) s k s k pm (ρC) pm ∂θ 1 ∂t = ∂ 2 θ 1 ∂ x 2 + ∂ 2 θ 1 ∂z 2 , ( 51
)
and the z = -δ boundary condition in Eq. ( 18) is rendered homogeneous. Finally, the basic temperature profile in Eq. ( 19) becomes,

θ = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ - z d (-δ ≤ z ≤ 0), -z (0 ≤ z ≤ 1), -1 - z -1 d (1 ≤ z ≤ 1 + δ), ( 52
)
where the only change is in sublayer 1. The perturbation analysis described in §3 and the dispersion relation in §4 now apply to this alternative situation without change.

Fig. 2

 2 Fig. 2 Basic state temperature profiles for the cases with δ = 0.5: d = 0.25, 0.5 (short dashes), d = 1 (continuous) and d = 2, 4 (long dashes)

Fig. 3

 3 Fig. 3 Onset profiles for the streamfunction (lower curves) and the temperature (upper curves) for the wavenumber, k = 2: a δ = 0.4 with d = 0.01 (dotted), 0.1, 1, 10 and 100; b d = 5 with δ = 0.01, 0.03, 0.1, 0.3 and 1

Fig. 4

 4 Fig. 4 Neutral curves as a function of k for the indicated values of values of δ. Curves correspond to d = 10 -3(lowest, dotted line), 10 -2 , 10 -1 , 10 0 , 10 1 , 10 2 and 10 3 (uppermost)

d

  ∼ 10 0 10 0 < δ,d > 10 0 : dδ ∼ 1, d < 10 0 : d ∼ δ.

Fig. 7

 7 Fig. 7 Contours of Ra c and k c as functions of log 10 d and log 10 δ. Integer values of Ra c are used, while 20 equal intervals between 0 and π are used for k c . The dotted lines correspond to the CT/CHF case, Ra c = 27.097628 and k c = 2.326215

  k f , k pm , k s Conductivities of fluid, porous medium and solid, respectively

	K	Permeability
	p	Pressure
	q	Rate of heat generation
	Ra	Darcy-Rayleigh number
	t	Time
	T	Dimensional temperature
	u, w	Horizontal and vertical velocities, respectively
	x, z	Horizontal and vertical coordinates, respectively
	Greek characters
	β	Thermal expansion coefficient
	δ	Depth ratio
	T	Temperature drop across the porous layer
	θ	Nondimensional fluid temperature
		Perturbation temperature
	μ	Dynamic viscosity
	ρ	Density
	ψ	Streamfunction
		Perturbation streamfunction
	Superscripts and subscripts
	basic	Basic state
	f	Fluid
	pm	Porous medium
	ref	Lower interface temperature (basic state)
	upper	Upper surface
	c	Critical value
	s	Solid
	1	Lower solid layer
	2	Porous layer
	3	Upper solid layer
		Derivative with respect to z
	ˆDimensional quantity

= E cosh(kz) + tanh(kδ) sinh(kz) . (34) Application of the interface conditions at z = 0, namely that 1 =

and d 1 = 2 allows us to eliminate the arbitrary constant, E, to yield2 (0) = 2 (0)kd tanh(kδ),(35)which is regarded as a boundary condition for 2 that models perfectly the variation of 1 in the lower sublayer. A similar analysis for

yields 2 (1) = -2 (1)kd coth(kδ).(36)

= A cosh λz + B sinh λz -C cos σ z -D sin σ z,(38)