Nonlinear Convection in a Partitioned Porous Layer

: Convection in a partitioned porous layer is considered where the thin partition causes a mechanical isolation of the two identical sublayers from one another, but heat may neveretheless conduct freely. An unsteady solver that employs the multigrid method is employed to determine steady-state strongly nonlinear for values of the Darcy–Rayleigh number up to eight times its critical value. The predictions of linear stability theory are conﬁrmed and the accuracy of the computations are carefully monitored and controlled. It is found that the wavenumber for which the maximum rate of heat transfer is attained at any chosen value of the Darcy–Rayleigh number, Ra increases quite strongly from roughly 2.33 at onset to 6.25 when Ra = 200. It is also found that convection generally cannot take place with wavenumbers which are close to the left-hand branch of the neutral stability curve because nonlinear interactions favour modes selected from higher harmonics.


Introduction
Many authors have studied the effect of different types of layering on the convective instability caused by heating a layer of fluid-saturated porous media from below. The earliest studies begin with those of Georgitzha [1], who considered two cases: a two-layer system with slightly different permeabilities, and one with a weakly spatially varying permeability. Masuoka et al. [2] computed convection patterns in a two-layer system, and Rana, Horne and Cheng [3] modelled convection in three-layer system that represented the Pahoa reservoir in Hawaii. A much more systematc approach was then taken by McKibbin and O'Sullivan [4,5] who considered both the linearized theory and a two-dimensional weakly nonlinear theory. Of interest is the fact that, with two sublayers, it is possible to determine parameter ranges (i.e., of relative permeabilities and sublayer thicknesses) for which the neutral stability curve has two minima, rather than having the more usual unimodal form. Rees and Riley [6] revisited the weakly nonlinear theory in order to determine conditions on those parameters for which two-dimensional rolls do not form the preferred convection pattern, but where cells with square planform appear instead.
Other types of layering involve the sandwiching of the porous layer with solid but heat-conducting layers both above and below. The analysis of Mojtabi and Rees [7], who considered identical solid bounding layers, found that the onset criterion then depends very strongly indeed on the nature of the solid layers. At one extreme, the composite layer mimics well the classical Darcy-Bénard problem by having the critical Darcy-Rayleigh number, Ra c 4π 2 and a critical wavenumber, k c π; in general, this arises when the solid layers are relatively highly conducting. At the other extreme, the layer corresponds to Ra c 12 and k c = 0-this tends to happen when the solid layers are relatively poor conductors of heat (see, for example, the work of Falsaperla et al. [8] who consider the combined effect of rotation and Robin boundary conditions). Other examples of this type of work include the weakly nonlinear analyses of Riahi [9] and Rees and Mojtabi [10], the combined

Governing Equations
We consider convection in a porous layer which is heated from below and which is split into two identical sublayers by the presence of an impermeable, infinitesimally thin and perfectly heat conducting interface; see Figure 1. The height of the layer is 2L and therefore the sublayers each have height, L. The porous medium is isotropic and otherwise homogeneous with pemeability, K, and the dynamic viscosity of the fluid is µ. The lower surface of the layer is held at the uniform temperature, T h , while the upper surface is held at T c where T h > T c . We define a reference temperature drop according to and this represents the temperature drop across one of the two sublayers.
2014), nonlinear computations (Saleh et al. 2011), and internal conducting layers (Jang and Tsai 1988, Postelnicu 1999, Patil and Rees 2014, Straughan 2014 Of more relevance to the present work are the analyses of Genç and Rees (2011) and Rees and Genç (2011). The former paper considers a porous medium with two identical sublayers which are mechanically isolated from one another by a thin impermeable interface, where independent pressure gradients are applied to the two sublayers. The latter considers a system of N identical sublayers again with thin impermeable interfaces, but now considered as a free convective linear stability analysis. It was found that the large-N limit yields a critical Darcy-Rayleigh number and wavenumber which are more usually associated with the constant-heat-flux variant of the single-layer Darcy-Bénard problem, namely that Ra c ∼ 12 and k c ∼ 0. A weakly nonlinear analysis by Rees et al. (2014) shows that convection is always two-dimensional immediately post-onset. In the present paper we consider the two-sublayer case but extend the analysis of Rees and Genç (2011) to strongly nonlinear convection using numerical methods. Our aim here, therefore, is to determine the flow and temperature fields and to determine those wavenumbers which maximise the rate of heat transfer as a function of the Darcy-Rayleigh number.

Governing Equations
We consider convection in a porous layer which is heated from below and which is split into two identical sublayers by the presence of an impermeable, infinitesimally thin and perfectly heat conducting interface; see Figure 1. The height of the layer is 2L and therefore the sublayers each have height, L. The porous medium is isotropic and otherwise homogeneous with pemeability, K, and the dynamic viscosity of the fluid is µ. The lower surface of the layer is held at the uniform temperature, T h , while the upper surface is held at T c where T h > T c . We define a reference temperature drop according to and this represents the temperature drop across one of the two sublayers. We shall confine our attention to two-dimensional flow only; this is not a restriction for moderate values of the Darcy-Rayleigh number because the weakly nonlinear analysis of Rees et al. (2014) shows that the preferred convection pattern in the postcritical range is expected to be two-dimensional at first. On taking the velocity field to beû = (û, 0,ŵ), the conservation of mass yields, ∂û ∂x and Darcy's law, subject to the Boussinesq approximation, takes the form Figure 1. Depicting a porous layer split into two identical sublayers by an infinitesimally thin and impermeable interface.
We shall confine our attention to two-dimensional flow only; this is not a restriction for moderate values of the Darcy-Rayleigh number because the weakly nonlinear analysis of Rees et al. [19] shows that the preferred convection pattern in the postcritical range is expected to be two-dimensional at first. On taking the velocity field to beû = (û, 0,ŵ), the conservation of mass yields: and Darcy's law, subject to the Boussinesq approximation, takes the form The equation governing the transport of heat is given by The upper and lower surfaces together with the interface are impermeable and, therefore, we need to applyŵ = 0 atẑ = 0, L, and 2 L.
At the interface,ẑ = L, the continuity of both the temperature and the heat flux (i.e., temperature gradient) needs to be imposed.
The equations and boundary conditions are rendered dimensionless using the following transformations: and then the governing equations become and Given that the flow is two-dimensional, we may define the streamfunction, ψ, using which means that Equation (8) is satisfied automatically. Thus the elimination of p between Equations (9) and (10) means that the nondimensional governing equations for the porous medium and the boundary conditions become In the above the Darcy-Rayleigh number is given by and it is therefore based upon the height and the mean temperature drop across one sublayer. This convention was introduced in Genç and Rees [17] and Rees and Genç [18] in order to be able to compare critical values of Ra easily between the single-layer classical Darcy-Bénard problem and cases where there are any number of sublayers stacked above one another. Should one require a Darcy-Rayleigh number which is based upon the thickness and the temperature drop across the whole of the layer, rather than one sublayer, then it will be precisely 4 Ra.

Numerical Method
Equations (13) and (14) were approximated using a second order accurate finite difference discretisation on a uniform grid. In general, only one convection cell in the horizontal direction was computed, and therefore the boundary conditions given in Equation (15) were supplemented by the following sidewall conditions: which are exactly equivalent to one cell in a layer of infinite width. The heat transport equation was solved using the DuFort-Frankel scheme for each timestep and the computation proceeded until it converged to a steady convecting state, which was always found to happen whenever Ra ≤ 200. Indeed, computations undertaken using the critical wavenumber (i.e., x max = π/k c = 1.350517) suggest that the flow becomes unsteady when the Darcy-Rayleigh number is in the vicinity of 300. The Neumann conditions for temperature on the sidewalls were applied using the fictitious point method. Finally, the continuity of both the temperature and the temperature gradient at the interface at z = 1 needs to be applied. It is clear that the only potentially problematic term is ψ z θ x . This represents the advection of heat along the interface. The ψ z term was discretised using a one-sided forward difference approximation to represent the flow along the upper surface of the interface, and a one-sided backward difference for the flow along the lower surface. These were then averaged, and the resulting formula is then identical to the standard central difference approximation which one would apply in the absence of the interface. We note that, although one-sided differences have first order accuracy, our application of this is on a one-dimensional subset of a two-dimensional domain, and therefore second order overall accuracy is maintained; this is demonstrated numerically in the Appendix.
The Poisson's equation for the streamfunction was then solved at each timestep using the multigrid correction scheme algorithm with line-relaxation Gauss-Seidel as the basic smoother. It was essential to use line-relaxation because cell-aspect ratios deviated from being close to unity when wavenumbers were small and large. In the initial development of the code, the multigrid algorithm was applied twice: once in the upper cavity and once in the lower cavity.
A detailed discussion of the choice of the chosen grid may be found in the Appendix and therefore it suffices to say that a 48 × 64 grid was used in almost all cases. Four slightly different versions of the numerical code were developed. Code A finds the steady-state solution for a single pair of values of Ra and k, and the solution is used to draw the contours of the streamfunction and temperature profiles. Code B solves for a monotonically varying set of values for Ra for a given wavenumber. Code C varies k while Ra is fixed. Finally, Code D finds that value of k which maximises the rate of heat transfer for the given value of Ra, but then uses that information as an initial iterate for the next value of Ra, and so on. In this last code, k is incremented forward by 0.01 until the maximum value of Nu is found, and then one further case is executed. A parabola is then fitted through the maximum value and the two values on either side, and the maximum found together with its corresponding wavenumber. During the testing of the code and the investigation of the qualitative behaviour of the flow, it became apparent that the following centro-symmetry relation is satisfied by all realised flows when Ra ≤ 200: Therefore, all four versions of the codes described above were then modified to ensure that this symmetry is obeyed. Thus, it was necessary only to compute ψ and θ in the lower sublayer, thereby doubling the speed of execution of the codes.
The chief measure of the strength of the convection, which is induced is the Nusselt number. This is defined to be the mean rate of heat transfer at the lower surface, and it is given by This integral is approximated numerically using a one-sided forward difference for the derivative and the Trapezium rule for the integral. The approximation to the derivative turns out to have second order accuracy because it is easily shown that ∂ 2 θ/∂z 2 = 0 on z = 0 using Equation (14).

Linearised Theory
The context of the present paper is the linearised theory contained in Rees and Genç [18], which deals with an arbitrary number of identical sublayers, each pair of which is separated by an infinitesimally thin impermeable interface. Dispersion relations, the solutions of which yield the neutral stability curves, were derived and it was found that, for N sublayers, the full dispersion relation could be factorised into N products that have an essentially identical form except for a numerical coefficient. For two layers, the present case, this factorisation yields the following: where The neutral curves corresponding to each of the expressions in Equation (20) are shown in Figure 2 where the first four modes are displayed. The first of the relations in Equation (20) may be rearranged into the form: which is the well-known neutral curve for the classical single-layer Darcy-Bénard problem. In Figure 2, the second and fourth modes correspond to n = 1 and n = 2 in Equation (22); this observation extends in the obvious way to higher modes. Rather strangely, the second dispersion relation in Equation (20) also applies to a single-layer case, but it is one where one of the boundaries is held at a constant temperature while the other is subject to a constant heat flux. In Figure 2, the first and third curves correspond to the first two solutions of this latter dispersion relation. Concentrating now on the first mode to appear as the Darcy-Rayleigh number increases, the disturbance temperature profile for the present two-layer configuration is, therefore, such that it has a zero vertical derivative at the interface. As will be seen below, this property does not carry forward into the nonlinear regime. Therefore, it is not possible to replace the present two-layer system with this particular one-layer system for the purposes of computing nonlinear convection. The neutral curves corresponding to each of the expressions in Eq. (20) are shown in Figure 2 where the first four modes are displayed. The first of the relations in Eq. (20) may be rearranged into the form, which is the well-known neutral curve for the classical single-layer Darcy-Bénard problem. In Fig. 2 the 2 nd and 4 th modes correspond to n = 1 and n = 2 in Eq. (22); this observation extends in the obvious way to higher modes. Rather strangely, the second dispersion relation in Eq. (20) also applies to a single-layer case but it is one where one of the boundaries is held at a constant temperature while the other is subject to a constant heat flux. In Fig. 2 the 1 st and 3 rd curves correspond to the first two solutions of this latter dispersion relation. Concentrating now on the first mode to appear as the Darcy-Rayleigh number increases, the disturbance temperature profile for the present two-layer configuration is therefore such that it has a zero vertical derivative at the interface. As will be seen below, this property does not carry forward into the nonlinear regime. Therefore it is not possible to replace the present two-layer system with this particular one-layer system for the purposes of computing nonlinear convection.   Figure 3 displays a selection of streamlines and isotherms, which were produced using code A in order to show how the flow evolves as the Darcy-Rayleigh number increases. We have chosen to use x max = 1 (i.e., k = π), so that each of the two sublayers occupy a computational domain with a unit aspect ratio. For this domain, convection begins once Darcy-Rayleigh number is larger than 29.209 977. When Ra = 30, we see that the isotherms deviate by only a small amount from the horizontal, and this indicates that a weak clockwise circulation is induced in each of the two sublayers. As Ra increases, the flow strength also increases, as one would expect for a supercritical bifurcation, and this is shown by the increasing deformation of the isotherms. An apparent up/down symmetry about z = 1 in the streamlines is lost when Ra is sufficiently large, and we now see the centro-symmetry referred to above in Equation (18).

Nonlinear Flows
An alternative way of viewing this flow is shown in the lower row of Figure 3 which displays the perturbation isotherms, i.e., When Ra = 30, which is just supercritical, we see that the perturbation isotherms are close to obeying the up/down symmetry, which is satisfied by the linearised disturbance at onset. At onset, these isolines are vertical as they pass through the interface. Thus, it is this symmetry that is referred to above as corresponding to the classical Darcy-Bénard layer with a constant temperature on one surface and a constant heat flux on the other. Now, as Ra increases, we see the hot spots in θ pert rising on the left of each sublayer and the cold spots descending. We also see the beginnings of boundary-layer formation on the horizontal surfaces.
When Ra = 30 we see that the isotherms deviate by only a small amount from the horizontal, and this indicates that a weak clockwise circulation is induced in each of the two sublayers. As Ra increases, the flow strength also increases, as one would expect for a supercritical bifurcation, and this is shown by the increasing deformation of the isotherms. An apparent up/down symmetry about z = 1 in the streamlines is lost when Ra is sufficiently large, and we now see the centro-symmetry referred to above in Eq. (18). Ra = 30 Ra = 50 Ra = 100 Ra = 200 An alternative way of viewing this flow is shown in the lower row of Fig. 3 which displays the perturbation isotherms, i.e.
When Ra = 30, which is just supercritical, we see that the perturbation isotherms are close to obeying the up/down symmetry which is satisfied by the linearised disturbance at onset. At onset, these In Figure 4, the Darcy-Rayleigh number has been set at Ra = 200, the largest of which is considered here, and the values of x max are 0.5, 1, 1.5 and 2, which are equivalent to k = 2π, π, 2 3 π and 1 2 π. While these streamline and isotherm plots appear superficially to be roughly the same, they hide some unusual behaviour. The largest value of ψ corresponds to x max = 2 out of these four cases, but the largest value of Nu corresponds to x max = 0.5. When x max = 2, the domain is larger and offers less obstruction to the flow, given that most of the time any randomly-chosen fluid particle is moving horizontally. However, this horizontal movement also allows for heat to conduct further into the interior of each sublayer, thereby decreasing the temperature gradient and the Nusselt number. That this is the case is confirmed by the isolines of θ pert . Close to the top left-hand side of the upper sublayer, it is clear that the θ pert isotherms are almost identical, and that the close spacing of these isotherms indicates that there is high local rate of heat transfer on the upper surface. However, when the two extreme cases are compared, we see that this high rate of heat transfer is maintained over almost all of the upper surface when x max = 0.5, but when x max = 2, the local rate of heat transfer has already changed sign before the top right corner has been reached. Figure 5 shows how the Nusselt number varies with Ra for a choice of values of the wavenumber. The continuous curves represent values of k which are above the critical value, and all of them show that there is a straightforward supercritical bifurcation to a strongly convecting flow as Ra increases, where Nu = 1 represents the conducting state. Confining ourselves to these curves, it is clear that the value of Ra at which the onset of convection occurs decreases as k decreases, precisely in line with the neutral curve shown in Figure 2. Numerically, our extrapolated onset criteria based on extrapolations of the value of Nu back to Nu = 1, which corresponds to conduction, agree with the exact dispersion relation values to within 0.15% at most. It is of interest to note that these curves cross one another as Ra increases, and it is clear that the wavenumber at which Nu achieves its largest value depends on the Darcy-Rayleigh number; this will be investigated in more detail in Figure 6, below.
In Fig. 4 the Darcy-Rayleigh number has been set at Ra = 200, the largest which is considered here, and the values of x max are 0.5, 1, 1.5 and 2, which are equivalent to k = 2π, π, 2 3 π and 1 2 π. While these streamline and isotherm plots appear superficially to be roughly the same, they hide some unusual behaviour. The largest value of ψ corresponds to x max = 2 out of these four cases, but the largest value of Nu corresponds to x max = 0.5. When x max = 2, the domain is larger and offers less obstruction to the flow, given that most of the time any randomly-chosen fluid particle is moving horizontally. But this horizontal movement also allows for heat to conduct further into the interior of each sublayer, thereby decreasing the temperature gradient and the Nusselt number. That this is so is confirmed by the isolines of θ pert . Close to the top left hand side of the upper sublayer it is clear that the θ pert isotherms are almost identical, and that the close spacing of these isotherms indicates that there is high local rate of heat transfer on the upper surface. However, when the two extreme cases are compared, we see that this high rate of heat transfer is maintained over almost all of the upper surface when x max = 0.5, but when x max = 2, the local rate of heat transfer has already changed sign before the top right corner has been reached.
x max = 0.5 x max = 1 x max = 1.5 x max = 2   where Nu = 1 represents the conducting state. Confining ourselves to these curves, it is clear that the value of Ra at which the onset of convection occurs decreases as k decreases, precisely in line with the neutral curve shown in Fig. 2. Numerically, our extrapolated onset criteria based on extrapolations of the value of Nu back to Nu = 1, which corresponds to conduction, agree with the exact dispersion relation values to within 0.15% at most. It is of interest to note that these curves cross one another as Ra increases, and it is clear that the wavenumber at which Nu achieves its largest value depends on the Darcy-Rayleigh number; this will be investigated in more detail in Fig. 6, below. The dotted line in Fig. 5 corresponds to k = 2, which is less than the value of k c . It has the same behaviour as those drawn with continuous lines, although the value of Ra at which convection starts has begun to increase once more. But the dotted-dashed line, for which k = 1.5, terminates at approximately Ra = 75; it is not possible to find steady solutions using our time-dependent code at smaller values of Ra. In our simulations for such cases we find that, while the one-cell perturbation grows at first, it eventually transforms into a three-cell flow where each of these new cells corresponds to the one cell which is obtained when solving using a wavenumber which is three times larger. For relatively small wavenumbers the width of the domain, x max is now quite wide and therefore we obtain a flow where the three new cells each occupy a region of horizontal width 1 3 x max . The reason for this will be discussed after Fig. 6 is introduced.   Figure 6 shows how the Nusselt number varies with wavenumber for a choice of values of Ra, and these were obtained using code C. In all cases Nu = 1 corresponds to conduction, and departure from this value corresponds to onset. The general behaviour of Nu as a function of k follows that of the classical Darcy-Bénard layer shown in de la Torres and Busse (1995), where the maximum rate of heat transfer is obtained at some value of the wavenumber which is intermediate between the two values corresponding to the neutral curve at the chosen value of Ra; this wouldn't necessarily be the case if parts of the neutral curve correspond to a subcritical bifurcation. The locus of maxima is shown in Fig. 6 as a dashed line, and these values were obtained using code D. As Ra increases, the wavenumber at which Nu attains its maximum value increases as Ra increases, which is the same qualitative behaviour as was found in de la Torres and Busse (1995). The detailed variation of the maximising value of k as Ra increases is shown in Fig. 7. This variation shows a surprising change in direction once Ra rises above 150, and therefore the accuracy of our computations were re-checked using a coarser grid (24 × 32 -double the step lengths in the cordinate directions), and a finer grid (96 × 128 -approximately 0.7 times the respective step lengths). The line corresponding to the coarser grid terminates at Ra = 180 because of resolution difficulties, but the difference between the coarse values and those of our chosen grid are quite small when Ra < 100. On the other hand, the finer grid is in close agreement with our chosen grid. At Ra = 200 we obtain k = 6.2728 using our chosen 48 × 64 grid and k = 6.2965 using the finer grid, a relative change of less than 0.4%. The dotted line in Figure 5 corresponds to k = 2, which is less than the value of k c . It has the same behaviour as those drawn with continuous lines, although the value of Ra at which convection starts has begun to increase once more. However, the dotted-dashed line, for which k = 1.5, terminates at approximately Ra = 75; it is not possible to find steady solutions using our time-dependent code at smaller values of Ra. In our simulations for such cases, we find that, while the one-cell perturbation grows at first, it eventually transforms into a three-cell flow where each of these new cells corresponds to the one cell that is obtained when solving using a wavenumber that is three times larger. For relatively small wavenumbers, the width of the domain, x max is now quite wide, and, therefore, we obtain a flow where the three new cells each occupy a region of horizontal width 1 3 x max . The reason for this will be discussed after Figure 6 is introduced. Figure 6 shows how the Nusselt number varies with wavenumber for a choice of values of Ra, and these were obtained using code C. In all cases Nu = 1 corresponds to conduction, and departure from this value corresponds to onset. The general behaviour of Nu as a function of k follows that of the classical Darcy-Bénard layer shown in de la Torres and Busse [20], where the maximum rate of heat transfer is obtained at some value of the wavenumber which is intermediate between the two values corresponding to the neutral curve at the chosen value of Ra; this would not necessarily be the case if parts of the neutral curve correspond to a subcritical bifurcation. The locus of maxima is shown in Figure 6 as a dashed line, and these values were obtained using code D. As Ra increases, the wavenumber at which Nu attains its maximum value increases as Ra increases, which is the same qualitative behaviour as was found in de la Torres and Busse [20]. The detailed variation of the maximising value of k as Ra increases is shown in Figure 7. This variation shows a surprising change in direction once Ra rises above 150, and therefore the accuracy of our computations were re-checked using a coarser grid (24 × 32-double the step lengths in the cordinate directions), and a finer grid (96 × 128-approximately 0.7 times the respective step lengths). The line corresponding to the coarser grid terminates at Ra = 180 because of resolution difficulties, but the difference between the coarse values and those of our chosen grid are quite small when Ra < 100. On the other hand, the finer grid is in close agreement with our chosen grid. At Ra = 200, we obtain k = 6.2728 using our chosen 48 × 64 grid and k = 6.2965 using the finer grid, a relative change of less than 0.4%. The second feature found in Fig. 6 is that, as the wavenumber increases, our unsteady code is unable to obtain solutions all the way back to the onset criterion. The various curves terminate before reaching the value of k representing the neutral curve, the sole exception being when Ra = 30. This is the same phenomenon as is mentioned above with regard to Fig. 5. We attempt to explain this using the case x max = π (i.e. k = 1). When k = 1, convection begins when Ra = 47.832420, whereas when k = 3, it begins when Ra = 28.600239. From the point of view of small-amplitude disturbances, a k = 1 disturbance at Ra = 200 and a k = 3 disturbance at the same value of Ra will have different exponential growth rates. On applying such a theory, details of which are omitted here, we find that the k = 1 disturbance grows with an amplitude that proportional to exp(12.15538t) while the k = 3 disturbance grows as exp(78.70295t); these growth rates remain accurate until nonlinear effects become significant and growth is attenuated. Now, if the k = 1 disturbance is initiated, then the nonlinear terms in the heat transport equation combine to yield terms which include the k = 3 disturbance as a component, and this grows more quickly than the original disturbance. It appears, therefore, that if k takes a value which is sufficiently close to the left hand branch of the neutral curve, then the 3k component, which has a substantially greater growth rate, will eventually overwhelm the original disturbance and form the ultimate steady state. At larger values of k, the difference between the two growth rates is much smaller, and the original disturbance continues to develop towards steady state. In general, a nonlinear interaction forming a 2k component has a different symmetry, and this term will generally decay and play only a passive role.
The dotted line in Fig. 6 demonstrates the above k to 3k transition as k decreases in the case Ra = 150. The steady-state solution which we obtain when Ra = 150 with k decreasing undertakes a sudden transition at k = 1.18 to the 3k solution, as described above. This dotted curve is, in effect, a three-times compressed version of the continuous curve. Subject to the fact that this '3k' solution The second feature found in Figure 6 is that, as the wavenumber increases, our unsteady code is unable to obtain solutions all the way back to the onset criterion. The various curves terminate before reaching the value of k representing the neutral curve, the sole exception being when Ra = 30. This is the same phenomenon as is mentioned above with regard to Figure 5. We attempt to explain this using the case x max = π (i.e., k = 1). When k = 1, convection begins when Ra = 47.832420, whereas when k = 3, it begins when Ra = 28.600239. From the point of view of small-amplitude disturbances, a k = 1 disturbance at Ra = 200 and a k = 3 disturbance at the same value of Ra will have different exponential growth rates. On applying such a theory, details of which are omitted here, we find that the k = 1 disturbance grows with an amplitude that proportional to exp(12.15538t) while the k = 3 disturbance grows as exp(78.70295t); these growth rates remain accurate until nonlinear effects become significant and growth is attenuated. Now, if the k = 1 disturbance is initiated, then the nonlinear terms in the heat transport equation combine to yield terms that include the k = 3 disturbance as a component, and this grows more quickly than the original disturbance. It appears, therefore, that if k takes a value which is sufficiently close to the left-hand branch of the neutral curve, then the 3k component, which has a substantially greater growth rate, will eventually overwhelm the original disturbance and form the ultimate steady state. At larger values of k, the difference between the two growth rates is much smaller, and the original disturbance continues to develop towards steady state. In general, a nonlinear interaction forming a 2k component has a different symmetry, and this term will generally decay and play only a passive role.
The dotted line in Figure 6 demonstrates the above k to 3k transition as k decreases in the case Ra = 150. The steady-state solution that we obtain when Ra = 150 with k decreasing undertakes a sudden transition at k = 1.18 to the 3k solution, as described above. This dotted curve is, in effect, a three-times compressed version of the continuous curve. Subject to the fact that this '3k' solution at k = 1.18 has relatively poor spatial resolution, its Nusselt number is essentially the same as was obtained for the 'k' solution when k = 3.54.

Conclusions
The literature is replete with linear stability analyses for convection in a variety of layered systems involving porous media, but apart from some very early work [1][2][3] on fairly restrictive systems composed of porous sublayers, no nonlinear computations have been reported. The present paper has considered a very simple modification to the classic Darcy-Bénard problem, one where a horizontal, infinitesimally thin but impermeable interface is placed centrally within the porous layer. This has the effect of modifying the linearised stablity problem to one which is exactly equivalent to that variant of the Darcy-Bénard problem where one surface is held at a constant temperature, but the other is subject to a constant heat flux. We have found that this analogy does not extend into the nonlinear regime, where, for Darcy-Rayleigh numbers that are below 200, the flow and temperature fields satisfy appropriated forms of centro-symmetry, rather than the up-down symmetry which applies to the linearised disturbances. We find that the most favoured wavenumber, at least in terms of that which maximises the mean rate of heat transfer, increases as the Darcy-Rayleigh number increases, showing the narrow cells transport more heat than those of roughly O(1) aspect ratio.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Accuracy Considerations
The aim of this Appendix is the justification of the choice of numerical grid that has been used for the bulk of the present paper.
The critical Darcy-Rayleigh number for this two-layer configuration is Ra c = 27.097628 and the critical wavenumber is k c = 2.326215. The value of this wavenumber is such that one single cell of the corresponding convection pattern has width, π/k c = 1.350517. Given that we are using a multigrid methodology to solve the Poisson's equation for the streamfunction, the number of intervals in each direction must have factors which are powers of two. Therefore, we solved the equations of motion on the three different uniform grids: 12 × 16, 24 × 32 and 48 × 64 for a single-cell solution at the critical wavenumber, while the Darcy-Rayleigh number was varied within the range, Ra c < Ra < 200.
The variation of the computed Nusselt number with Ra is given in Figure A1 for the three grids, and we see that there is very little difference in the value of Nu when Ra = 200 when considering the two finest grids. The coarsest grid clearly has a substantial error associated with it. In particular, Table A1 gives that data for Ra = 200, which is the most extreme case considered here, and it is clear that they are commensurate with a second order accurate numerical scheme.

Appendix. Accuracy Considerations
The aim of this Appendix is the justification of the choice of numerical grid which has been used for the bulk of the present paper.
The critical Darcy-Rayleigh number for this two-layer configuration is Ra c = 27.097628 and the critical wavenumber is k c = 2.326215. The value of this wavenumber is such that one single cell of the corresponding convection pattern has width, π/k c = 1.350517. Given that we are using a multigrid methodology to solve the Poisson's equation for the streamfunction, the number of intervals in each direction must have factors which are powers of 2. Therefore we solved the equations of motion on the three different uniform grids: 12 × 16, 24 × 32 and 48 × 64 for a single-cell solution at the critical wavenumber, while the Darcy-Rayleigh number was varied within the range, Ra c < Ra < 200.  Figure A1. The value of Nu as a function of Ra for k = k c using the three grids: 12 × 16 (dotted), 24 × 32 (dashed) and 48 × 64 (continuous).
The variation of the computed Nusselt number with Ra is given in Fig. A1 for the three grids and we see that there is very little difference in the value of Nu when Ra = 200 when considering the two finest grids. The coarsest grid clearly has a substantial error associated with it. In particular Table A1 gives that data for Ra = 200, which is the most extreme case considered here, and it is clear that they are commensurate with a second order accurate numerical scheme.  Figure A1. The value of Nu as a function of Ra for k = k c using the three grids: 12 × 16 (dotted), 24 × 32 (dashed) and 48 × 64 (continuous). A second way of determining the accuracy of our computations (which also provides some validation of the coding) is to use the Nusselt number for values of Ra, which are just above the critical value. The numerical code was run using k = k c and convergence was deemed to have taken place when the maximum change in θ between timesteps is less than 10 −8 for 50 successive timesteps.
Given that the Nusselt number is a linear function of (Ra − Ra c ) according to weakly nonlinear theory (which is valid when Ra − Ra c is sufficiently small), we extrapolated back to Nu = 1 from two different postcritical solutions near Ra c in order to estimate the value of Ra c . The result of this extrapolation process is given in Table A2, and it is clear these data are also commensurate with a numerical method which has second order accuracy in space. Given that the order of accuracy of the numerical method is reproduced faithfully, not only in its prediction of the value of Ra c , but also in the extreme case of Ra = 200, we decided that it was unnecessary to adopt a still finer grid for the present work. Therefore, we have used a grid with 64 intervals in the vertical direction and 48 horizontally. Computations at higher values of the Darcy-Rayleigh number will require more resolution.