2.1. The Darcy–Brinkman model

Consider a two-dimensional system comprising a fluid-saturated porous medium of depth $h_p$ underlying an unconfined fluid region of depth $h_f$. We locate the centre of the cell at height $z=0$, such that the whole system lies in the range $-h_p \leqslant z \leqslant h_f$, as depicted in figure 1, and we introduce the length scale $h = (h_p+h_f)/2$. For the majority of this paper, we focus on the case of equal layer depths, where $h_p = h_f = h$.

Figure 1. Schematic cartoon of the set-up considered in this paper. In almost all cases considered here, we take the layer depths to be equal, so $h_p=h_f = h$.

The fluid has a dynamic viscosity $\eta$ and density $\rho$, and the porous medium is characterised by a uniform permeability $K_0$ and porosity $\phi _0< 1$. We extend the definition of the porosity – that is, the local volume fraction of fluid – to the whole domain by introducing the step function

(2.1)\begin{equation} \phi(z) = \left\{\begin{array}{@{}ll} \phi_0 & z<0, \\ 1 & z\geq 0. \end{array} \right. \end{equation}

The flow is described by the local fluid velocity $\boldsymbol{U}_{\ell} $ in the unconfined fluid layer and by the Darcy or volume-averaged flux $\boldsymbol{U}_d = \phi_0 \boldsymbol{U}_{\ell}$ in the porous medium. While the latter quantity is, by necessity, coarse-grained over a larger scale (multiple pore scales) than the former, for notational convenience we can nevertheless introduce a single quantity $\boldsymbol{U} = \phi \boldsymbol{U}_{\ell}$ that reduces to each of these limits in the relevant domain. We will work in terms of this mean flux $\boldsymbol{U} = (U,V)$ throughout the domain.

The flow is assumed to be incompressible everywhere, so

(2.2)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0. \end{equation}

In the fluid layer, the flow is governed by the Navier–Stokes equation,

(2.3)\begin{equation} \rho [\partial_t \boldsymbol{U} + \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U} ] ={-} \boldsymbol{\nabla} P + \mu \boldsymbol{\nabla}^2 \boldsymbol{U} + \rho \boldsymbol{g}, \end{equation}

where $P$ is the pressure, while in the porous layer, the flux instead obeys Darcy's law,

(2.4)\begin{equation} \boldsymbol{U} = \frac{K}{\mu}(- \boldsymbol{\nabla} P + \rho \boldsymbol{g}). \end{equation}

We simulate the two-layer system using a one-domain approach in which both porous and unconfined fluid regions are described by a single Darcy–Brinkman equation (Le Bars & Worster Reference Le Bars and Worster2006),

(2.5)\begin{equation} \rho \left[\partial_t \boldsymbol{U} + \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} \frac{\boldsymbol{U}}{\phi} \right] ={-} \phi \boldsymbol{\nabla} P + \mu \boldsymbol{\nabla}^2 \boldsymbol{U} + \phi \rho \boldsymbol{g} - \frac{\mu}{K} \phi \boldsymbol{U}, \end{equation}

where $1/K(z)$ is a step function that goes from $1/K_0$ for $z<0$ to zero for $z>0$. As shown by Le Bars & Worster (Reference Le Bars and Worster2006), the Darcy–Brinkman formulation of the two-layer problem can be retrieved by carrying out a coarse-grained average of the flow over a few typical pore scales $\sqrt {K_0}$. As a consequence, any spatial variation must have a typical length larger than $\sqrt {K}$ for the model to remain valid. The Navier–Stokes equation and Darcy's law are retrieved from the Darcy–Brinkman equation (2.5) in the unconfined fluid and porous layers, respectively. In the fluid layer, $\phi = 1$ and $1/K = 0$, which trivially gives the Navier–Stokes equation, whereas in the bulk of the porous medium, the damping term $-\mu \phi \boldsymbol {U} /K$ dominates over the inertial and viscous forces (provided $K_0$ is sufficiently small), leading to a balance between the damping, pressure and buoyancy terms that yields Darcy's law. Just below the interface, however, viscous effects become important in the porous layer as the flow matches to the unconfined region above. Acceleration remains negligible, and the equations reduce to

(2.6)\begin{equation} - \phi \boldsymbol{\nabla} P + \mu \boldsymbol{\nabla}^2 \boldsymbol{U} + \phi \rho \boldsymbol{g} - \frac{\mu}{K} \phi \boldsymbol{U} \sim 0. \end{equation}

Balancing the viscous resistance and Darcy drag indicates that local viscous forces play a role over a length $\ell _r = \sqrt {K_0}/\phi _0$ below the the interface – i.e. a few times the pore scale. These forces regularise the velocity profile between the unconfined fluid and the porous medium through a boundary layer of typical length $\ell _r$.

To conclude, we model the two-layer system with a one-domain formulation via the Darcy–Brinkman equation. We note that this is not the only option: another classical formulation of the problem, for example, is the one introduced by Beavers & Joseph (Reference Beavers and Joseph1967) where the fluid and the porous layer are treated separately and their coupling is accounted by a semi-empirical boundary condition linking vertical velocity gradients and the velocity difference between the fluid and porous layers. These different models both feature the regularisation boundary layer at the fluid-porous interface over a length $\sim \sqrt {K_0}$ which is corroborated by many experiments, in particular those of Beavers & Joseph (Reference Beavers and Joseph1967). There are, however, some known discrepancies between these two approaches that may not be restricted to the interface (e.g. Le Bars & Worster Reference Le Bars and Worster2006; Hirata et al. Reference Hirata, Goyeau, Gobin, Carr and Cotta2007, Reference Hirata, Goyeau, Gobin, Chandesris and Jamet2009b). As pointed out by Nield & Bejan (Reference Nield and Bejan2017), there is ongoing debate on which of these formulations is the most adequate to model flows in mixed porous-fluid layers, with definitive empirical evidence still lacking.

2.2. Heat transport

We use the Boussinesq approximation to account for the effect of temperature-dependent density in the momentum equations: variations in temperature affect the buoyancy force but do not affect the fluid volume via conservation of mass. We further assume that any changes in viscosity, diffusivity or permeability associated with temperature variation are negligible. Although some of these assumptions may be questionable in complex geophysical settings, they are made here to allow a focus on the basic physics of these two-layer convecting systems.

In particular, we restrict our study to linear variations of density with temperature according to

(2.7)\begin{equation} \rho = \rho_0 ( 1 - \alpha (T-T_0) ), \end{equation}

with $T_0$ being a reference temperature. The momentum equation under the Boussinesq approximation follows from substituting (2.7) into the buoyancy term of (2.5), while letting $\rho = \rho _0$ in the inertial terms. The temperature evolves according to an energy transport equation (Nield & Bejan Reference Nield and Bejan2013)

(2.8)\begin{equation} \bar{\phi}\partial_t T + \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \kappa \varLambda \nabla^2 T, \end{equation}

with

(2.9a–c)\begin{equation} \bar{\phi} \equiv \frac{(1-\phi) c_m \rho_m + \phi c \rho}{\rho c },\quad \kappa \equiv \frac{\lambda}{\rho c}\quad \text{and} \quad \varLambda \equiv \phi + ( 1-\phi)\frac{\lambda_m}{\lambda}, \end{equation}

where $c$ and $c_m$ are the heat capacity per unit of mass of the fluid and the porous matrix, respectively, $\rho _m$ is the density of the porous matrix, and $\lambda$ and $\lambda _m$ are the thermal conductivities of the fluid and the porous matrix, respectively. Equation (2.8) assumes local thermal equilibrium between the porous matrix and the fluid.

2.3. Boundary conditions

We consider a closed domain with imposed temperature on the upper and lower boundaries, as in a classical Rayleigh–Bénard cell (figure 1). Specifically, for the temperature we set

(2.10a,b)\begin{equation} T(z=h) = T_0,\quad T(z={-}h) = T_0 + \Delta T, \end{equation}

where $\Delta T >0$ is a constant. The upper and lower boundaries are rigid and impermeable, so

(2.11)\begin{equation} U(z={\pm} h) = V(z={\pm} h) = 0. \end{equation}

Note that Darcy's law would only permit one velocity boundary condition on the boundary of the porous region at $z=-h$, but the higher-order derivative in the viscous term in (2.5) allows for application of the no-slip condition in (2.11) as well. This extra condition will induce a boundary layer of thickness $\sim \sqrt {K_0}/\phi _0$ at the base of the domain, just like the boundary-layer region at the interface (see (2.6)). It is not clear whether such a basal boundary layer should arise in experimental situations. Irrespective of whether this boundary is a physically realisable feature, we find that it plays no dynamical role here provided it is thinner than any dynamical lengthscale of the flow (and, in particular, thinner than the thermal boundary layer that can form at the base of the domain, as discussed in § 2.6.)

The horizontal boundaries of the domain are assumed to be periodic with the width of the domain kept constant at $4h$.

2.4. Dimensionless equations and control parameters

To extract the dimensionless equations that govern the two-layer system, we use a ‘free-fall’ normalisation of (2.5) and (2.8), based on the idea that a balance between inertia and buoyancy governs the behaviour of the fluid layer. Such a balance yields the free-fall velocity in the unconfined layer,

(2.12)\begin{equation} {U^*}^2 = \alpha \Delta T g h \end{equation}

and the associated free-fall time scale is $T^* = h/U^*$. Scaling lengths with $h$, flux with $U^*$, time with $T^*$, temperature with $\Delta T$ and pressure with $U^*$, we arrive at dimensionless equations

(2.13a)\begin{gather} \partial_t \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \frac{ \boldsymbol{u}}{\phi} ={-} \phi \boldsymbol{\nabla} p + \sqrt{\frac{\mathit{Pr}}{\mathit{Ra}}} \boldsymbol{\nabla}^2 \boldsymbol{u} + \phi \theta \boldsymbol{e}_z - \frac{ \chi(z)}{\mathit{Da}} \sqrt{\frac{\mathit{Pr}}{\mathit{Ra}}} \phi \boldsymbol{u}, \end{gather}

(2.13b)\begin{gather}\bar{\phi} \partial_t \theta + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta = \frac{\varLambda}{\sqrt{\mathit{Ra} \mathit{Pr}}} \nabla^2 \theta, \end{gather}

(2.13c)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

where $\boldsymbol {u}$ is the dimensionless flux, $\theta = (T-T_0)/\Delta T$ is the dimensionless temperature and $\chi (z)$ is a step function that jumps from $1$ in $z<0$ to $0$ in $z>0$. In these equations, we have introduced three dimensionless numbers:

(2.14a–c)\begin{equation} \mathit{Da} \equiv \frac{K_0}{h^2},\quad \mathit{Pr} \equiv \frac{\nu}{\kappa} \quad \text{and}\quad \mathit{Ra} \equiv \frac{\alpha g\Delta T h^3 }{\nu \kappa}. \end{equation}

The Darcy number $\mathit {Da}$ is a dimensionless measure of the pore scale $\sqrt {K_0}$ relative to the domain size $h$, and is thus typically extremely small. The Rayleigh number quantifies the importance of the buoyancy forces relative to the viscous resistance in the unconfined fluid layer, and the focus of this work is on cases where $Ra \gg 1$. The Prandtl number compares momentum and heat diffusivities. The dimensionless layer depths $\hat {h}_p$ and $\hat {h}_f$ are also, in general, variables; as noted above, in the majority of computations shown here, we set these to be equal so that $\hat {h}_p = \hat {h}_f = 1$, but we include the general case in the theoretical discussion in § 4.

Note that with this choice of scalings, the dimensionless velocity scale in the fluid layer is $O(1)$, compared with $O(\sqrt {\mathit {Ra}} \,\mathit {Da} /\sqrt {\mathit {Pr}})$ in the porous layer. Given these scales, we can also introduce a porous Rayleigh number $\mathit {Ra}_p$ to describe the flow in the porous layer. The porous Rayleigh number is the ratio between the advective and diffusive time scales in the porous medium, which, from the advection–diffusion ratio in (2.13b), gives

(2.15)\begin{equation} \mathit{Ra}_p = \frac{\mathit{Ra} \,\mathit{Da}}{ \varLambda} = \frac{\alpha g\Delta T K_0 h }{\nu \varLambda \kappa}. \end{equation}

2.5. Further simplifying assumptions

We simplify the complexity of (2.13) by noting that in the bulk of either the fluid or the porous regions, $\phi$ cancels out of the equations (see for example (2.3) and (2.4)). The porosity only affects (2.13) in the narrow boundary-layer region immediately below the interface and at the base of the domain, where it controls the regularisation length $\sqrt {\mathit {Da}}/\phi$ (as shown by (2.6)). In the following, we thus take $\phi = 1$ in (2.13a); the only effect of this is to change the regularisation length at the interface, a regularisation that must anyway be smaller than any dynamical lengths for the model to remain valid, as discussed in more detail in § 2.7.

In the heat transport equation (2.13b), we reduce the number of control parameters by taking $\bar {\phi } = \varLambda = 1$. For hydrothermal systems, water flows through a silicate rock matrix. The thermal diffusivity is typically a factor of two larger in the matrix than in the fluid, while $\rho _m c_m \sim \rho c$. The parameters $\bar {\phi }$ and $\varLambda$ are thus order-one constants that do not vary substantially from one system to another. This is perhaps less true in some industrial applications like the transport of heat through the metallic foam (Su et al. Reference Su, Wade and Davidson2015), where the thermal conductivity can be at least a hundred times larger in the matrix than in the fluid. This would lead to a large value of $\varLambda$ and thermal diffusion would be greatly enhanced in the porous medium, which would reduce its ability to convect. We do not consider such cases here, although we will find that cases where the porous medium is dominated by diffusive transport can be easily treated theoretically, and the theory could be straightforwardly adapted to account for varying $\varLambda$. Finally, to reduce the complexity of this study and maintain a focus on the key features of varying the driving buoyancy forces (i.e. $\mathit {Ra}$) and the properties of the porous matrix (i.e. $\mathit {Da}$), we set the Prandtl number to be $\mathit {Pr} = 1$ throughout this work.

2.6. Limits on the control parameters

Several constraints must be imposed on the control parameters $\mathit {Ra}$ and $\mathit {Da}$ to ensure that the Darcy–Brinkman model remains valid. We give these constraints in their most general form here, but recall from the previous section that all solutions in this work take $\mathit {Pr} = \varLambda = 1$. First, the inertial terms must vanish in the porous medium, which demands that

(2.16)\begin{equation} \frac{\sqrt{\mathit{Ra}} \,\mathit{Da}}{\sqrt{\mathit{Pr}}} \ll 1; \end{equation}

that is, the velocity scale in the porous medium must be much less than the $O(1)$ velocity in the unconfined fluid layer. Second, the continuum approximation that underlies Darcy's law requires that any dynamical length scale of the flow in the porous layer must be larger than the pore scale $\sqrt {\mathit {Da}}$; equivalently, the Darcy drag term must always be larger than local viscous forces in the bulk of the medium. We expect the smallest dynamical scales to arise from a balance between advection and diffusion in (2.13b): such a balance, given typical velocity $\sim \sqrt {\mathit {Ra}}\mathit {Da}/\sqrt {\mathit {Pr}}$, yields a length scale $\sim \mathit {Ra}_p^{-1}$. In fact, simulations carried out in a porous Rayleigh–Bénard cell (Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2012) indicate that the narrowest structures of the flow, which are thermal boundary layers, have a thickness of at least $50 \mathit {Ra}_p^{-1}$. For these structures to remain larger than the pore scale $\sqrt {\mathit {Da}}$, we thus require

(2.17)\begin{equation} \mathit{Ra} \,\mathit{Da}^{3/2} \lesssim 50 \varLambda. \end{equation}

Note that the effect of violating this constraint is to amplify the importance of viscous resistance $\nabla ^2 \boldsymbol {u}$ within the porous medium in (2.13a), which would no longer reduce to Darcy's law.

Figure 2 provides an overview of the space of control parameters $\mathit {Ra}$ and $\mathit {Da}$, where these various limits are identified and the parameter values in our numerical simulations are indicated. This plot also shows a line that approximately marks the threshold of convective instability in the porous medium, whose importance is discussed in § 3.2 and which is theoretically quantified in § 4.2.2.

Figure 2. Domain of existence and validity of the two-layer model with respect to the control parameters $\mathit {Ra}$ and $\mathit {Da}$, given $\mathit {Pr} = \varLambda = 1$. Each dot represents a simulation. The line $\sqrt {\mathit {Ra}}\mathit {Da}= 1$ marks the limit beyond which inertial terms affect the flow in the porous medium, while the line $\mathit {Ra} \mathit {Da}^{3/2} = 50$ gives an estimate of the point at which the smallest flow structures in the porous medium become comparable to the pore scale. The line $\mathit {Ra} \mathit {Da} = \mathit {Ra}_p^c = 27.1$ corresponds to the threshold of thermal convection in a porous Rayleigh–Bénard cell with an open-top boundary as discussed in § 3.1.

2.7. Numerical method

The one-domain Darcy–Brinkman equations (2.13) are solved using the pseudo-spectral code Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020; Hester, Vasil & Burns Reference Hester, Vasil and Burns2021). The flow is decomposed into $N$ Fourier modes in the horizontal direction, while a Chebyshev polynomial decomposition is used in the vertical direction. Because the set-up is composed of two layers whose interface must be accurately resolved, each layer is discretised with its own Chebyshev grid, with $[M_p,M_f]$ nodes for the porous and fluid layers, respectively. This ensures enhanced resolution close to the top and bottom boundary as well as at the interface where the $\sqrt {\mathit {Da}}$ regularisation length must be resolved. Time evolution of the fields is computed with implicit-explicit methods (Wang & Ruuth Reference Wang and Ruuth2008): the nonlinear and Darcy terms in (2.13) are treated explicitly while viscosity and diffusion are treated implicitly. The numerical scheme for time evolution uses second-order backward differentiation for the implicit part and second-order extrapolation for the explicit part (Wang & Ruuth Reference Wang and Ruuth2008). The stability of temporal differentiation is ensured via a standard CFL criterion evaluating the limiting time step in the whole two-layer domain, with an upper limit given by $\sqrt {\mathit {Ra}} \,\mathit {Da}$, which is never reached in practice. Rather, with the control parameters and resolution considered here, the time step is limited by the non-zero vertical velocity at the $z=0$ interface where the vertical discretisation is refined. The range of Rayleigh and Darcy numbers in our simulations is shown in figure 2. Note, with reference to this figure, that we carried out a systematic investigation of parameter space where the porous layer is unstable by varying $\mathit {Ra}$ and $\mathit {Da}$ for various fixed values of the porous velocity scale $\sqrt {\mathit {Ra}}\,Da$.

We use a $\mathcal{C}^{\infty}$-smooth step function for $\chi{(z)}$ in (2.13). The smoothing of the step is performed over a length $ 0.75 \sqrt{Da}$ which is slightly smaller than the regularisation length to ensure that the smoothing does not play any dynamical role. We note that a sharp step function could also be used directly without changing the statistical properties of the simulated flows.

The majority of simulations were carried out in a set-up where the heights of both the porous and the fluid layers were equal $\hat {h}_p = \hat {h}_f = 1$, with resolutions $N \times [M_p,M_f] = 1024 \times [128, 256]$ below $\mathit {Ra} = 10^{8}$ and $1024 \times [256, 512]$ above. As discussed in the following sections, we find that the porous layer in general absorbs more than half of the temperature difference, and so the effective Rayleigh number in the fluid layer is typically somewhat smaller than $\mathit {Ra}$. We used two methods to initiate the simulations. In a few cases, the initial condition was simply taken as the diffusive equilibrium state throughout the domain, perturbed by a small noise. In most cases, however, we proceeded by continuation, using the final output of a previous simulation similar control parameters as an initial condition for a new simulation.

Comparison between the two methods showed that they yielded the same statistically steady state, but the continuation approach reached it in the shortest time. In all cases, we ran simulations over a time comparable to the diffusive time scale $\sqrt {\mathit {Ra}}$, to ensure that the flow had reached a statistically steady state.

3.1. Two different regimes depending on the stability of the porous medium

Figure 3 shows snapshots of the temperature field taken for different simulations that have reached a statistically steady state. The corresponding profiles of the horizontally and temporally averaged temperature, $\bar {\theta }(z)$, are shown in figure 4(a), while the mean interface temperature $\theta _i = \bar {\theta }(0)$ is shown in figure 4(b). The fluid layer is convecting in all simulations, as attested by the presence of plumes and by the mixing of the temperature field that tends to create well-mixed profiles of $\bar {\theta }$ in the bulk of the fluid. The porous layer, however, exhibits two different behaviours depending on the size of $\mathit {Ra}$. At low Rayleigh numbers ($\mathit {Ra} \leq 10^6$ in figure 3), the porous layer is dominated by diffusive heat transport: there are no hot or cold plumes in the temperature field in $z < 0$ (see figure 3a,b), while the horizontally averaged temperature profiles $\bar {\theta }(z)$ appear to be linear (figure 4a). The corresponding interface temperature monotonically decreases with $\mathit {Ra}$ (figure 4b). As the Rayleigh number is increased beyond $\mathit {Ra} \sim 10^7$, the behaviour of the flow in the porous layer changes as it also becomes unstable to convection. This is attested by the visible presence of plumes in figure 3(c,d) and by the flattening of the horizontally averaged temperature profiles in figure 4(a). The signature of this transition is also very clear in the evolution of the interface temperature $\theta _i$, which reverses from decreasing to increasing with $\mathit {Ra}$ around $\mathit {Ra} \sim 10^7$ (figure 4b). The value of the Rayleigh number at which porous convection emerges can be roughly estimated from the stability of a single porous layer. In a standard Rayleigh–Bénard cell with an open-top boundary, convection occurs if the porous Rayleigh number $\mathit {Ra}_p = \mathit {Ra} \mathit {Da}$ exceeds a critical value $\mathit {Ra}_p^c \simeq 27.1$ (Nield & Bejan Reference Nield and Bejan2017). At $\mathit {Da} = 10^{-5.5}$, the critical Rayleigh number $\mathit {Ra}$ such that $\mathit {Ra} \mathit {Da} = \mathit {Ra}_p^c$ is $\mathit {Ra} \simeq 8.6 \times 10^{6}$, which is reported in figure 4(b) and agrees well with the inversion of trend in $\theta _i$. We will return to a more nuanced form of this argument in § 4.2.2. Lastly, as the Rayleigh number is increased beyond the threshold of porous convection, the porous plumes become thinner and more numerous, a behaviour that is similar to standard Rayleigh–Bénard convection in porous media (Hewitt et al. Reference Hewitt, Neufeld and Lister2012). In addition, the porous plumes become increasingly narrower at their roots in the thermal boundary layer, which causes a local minimum in the temperature profiles (see figure 4a) that has also been observed in previous works on porous convection at large Rayleigh numbers Hewitt et al. (Reference Hewitt, Neufeld and Lister2012).

Figure 3. Snapshots of the temperature field at different Rayleigh numbers in two stable (top panels) and two convective (bottom panels) cases for the porous medium. The Darcy number is kept at $\mathit {Da} = 10^{-5.5}$. The colour scale is cut at $0.8$ to enhance the contrast in the fluid layer.

Figure 4. Temporally averaged quantities for $\mathit {Da} = 10^{-5.5}$. (a) Horizontally averaged temperature $\bar {\theta }(z)$ for values of $\mathit {Ra}$ below (dashed) and above (solid) the threshold of porous convection; and (b) mean interface temperature $\theta _i$. The red dashed line is the diffusive prediction of (4.8) and the black lines are asymptotic predictions obtained by solving (4.7) using $\mathcal {C}_{p} = 0.85$ and the marked values $\mathcal {C}_f$, as detailed in § 4.2.

Lastly, as the Rayleigh number is increased beyond the threshold of porous convection, the porous plumes become thinner and more numerous, a behaviour that is similar to standard Rayleigh–Bénard convection in porous media (Hewitt et al. Reference Hewitt, Neufeld and Lister2012). In addition, the porous plumes become increasingly narrow at their roots above the thermal boundary layer, which causes a local minimumin in the temperature profiles (see figure 4a) that has also been observed in previous works on porous convection at large Rayleigh number (Hewitt et al. Reference Hewitt, Neufeld and Lister2012). Finally, note that the temperature contrast across the interface decreases as the Rayleigh number is increased; this is because the contrast between velocities in the porous and unconfined layers decreases as $\sqrt{ Ra Da}$ increases. At fixed $Da$, the model assumption that there is a separation of scales between these velocities must break down if $Ra$ is made sufficiently large (see the discussion in § 2.6).

3.2. Characteristics of heat transport

The transition between the porous-stable and the porous-unstable cases can be further identified by the analysis of global heat transport across the system. Heat transport is characterised by the horizontally averaged heat flux $J(z) \equiv \overline {w \theta } - \mathit {Ra}^{-1/2} \bar {\theta }'$, which is constant with height in a statistically steady state. As is standard in statistically steady convection problems, we measure the time-averaged enhancement of the heat flux, compared with what it would be in a purely diffusive system, $\mathit {Ra}^{-1/2} /2$, via the Nusselt number $\mathit {Nu} \equiv 2 \sqrt {\mathit {Ra}} \langle J\rangle$, where the angle brackets indicate a long-time average. The Nusselt number (figure 5a) is strongly influenced by the transition to instability in the porous layer. In the porous-stable case ($\mathit {Ra} \lesssim 10^7$), $\mathit {Nu}$ appears to approach a horizontal asymptote $\mathit {Nu} = 2$, but once the porous layer is unstable, $\mathit {Nu}$ increases much more steeply beyond this value.

Figure 5. (a) Nusselt number $\mathit {Nu} = 2 \sqrt {\mathit {Ra}} \langle J\rangle$ characterising heat transport across the two-layer system, together with predictions from § 4 for the porous-stable case (4.8) (red dashed) and porous-unstable (4.7) (black dotted and dot–dashed, with $\mathcal {C}_{p} = 0.85$ and $\mathcal {C}_{f}$ as marked). The horizontal line marks $\mathit {Nu} = 2$. (b) Root-mean-squared vertical velocity amplitude $w_{rms}$ in the porous and fluid layers. The two grey lines indicate a scaling of $\sqrt {\mathit {Ra}}\mathit {Da}$ (the characteristic speed in the porous layer). (c) Ratio of the diffusive ($J_{d,p}$) and advective ($J_{a,p}$) fluxes to the total depth-averaged flux $\langle J \rangle$ in the porous layer. In all three panels, $\mathit {Da} = 10^{-5.5}$ and the vertical line marks the threshold of porous convection estimated using (4.9).

The behaviour below the threshold of convection arises from the flux being predominantly diffusive in the porous layer. The total flux through the system is thus bounded above by a state in which all of the temperature contrast is taken up across the porous layer and the interface temperature $\theta _i$ tends to zero. In this limit, $\langle J \rangle \to 1/\sqrt {\mathit {Ra}}$ and $\mathit {Nu} \to 2$. The decreasing $\theta _i$ in figure 4(b) as $Ra$ increases towards the threshold reflects the approach to this limit. In fact, we find that despite the porous medium appearing to be stable to convection below the threshold, small amplitude flows still exist in this regime, as can be seen from the non-zero root-mean-squared vertical velocity in the porous layer for all $Ra$, shown in figure 5(b). Nevertheless, by computing the relative diffusive and advective contributions to the flux through the porous medium, we confirm that these flows have a negligible impact on heat transport below the threshold (figure 5c). We interpret the weak secondary porous flows in the porous-stable regime as a consequence of the horizontal variations in the interface temperature imposed by fluid convection, which are clearly visible in figure 3(a,b). As shown in figure 5(b,c), the strength of the porous flow increases dramatically as the porous layer becomes unstable, and it is only then that the advective contribution to the flux in the porous medium becomes significant. We return to discuss this induced flow below onset in the porous layer in § 4.3, and defer more detailed discussion and prediction of the behaviour of $\mathit {Nu}(\mathit {Ra},\mathit {Da})$ until the following section.

3.3. Time scales, variability and statistically steady state

The governing equations (2.13) reveal three different time scales that govern the variability of the two-layer system considered here. The first is the turnover time scale in the fluid layer, which is $O(1)$ in our free-fall normalisation. The second is given by diffusion, $\tau _{\mathit {diff}} = \sqrt {\mathit {Ra}}$, and the third is the turnover time scale in the porous layer $\tau _p \sim (\sqrt {\mathit {Ra}}\mathit {Da})^{-1}$, which scales with the inverse of the porous speed scale. Because $\tau _p$ and $\tau _d$ measure advection and diffusion in the porous medium, these time scales are in a ratio $\tau _p = \tau _{\mathit {diff}}/\mathit {Ra}_p$. In addition, we recall that $\tau _p \gg 1$ is required for the porous layer to be in the confined limit and for the Darcy–Brinkman model to hold (see (2.16)). The turnover time scales should scale with the inverse of $w_{rms}$ in each layer, as can be observed in figure 5(b): in the fluid layer, $w_{rms} \sim O(1)$ with no systematic variation with $Ra$, while in the porous layer, $w_{rms} \propto \sqrt {\mathit {Ra}}$ at constant $\mathit {Da}$ in agreement with the scaling above. These two very different time scales are also clearly visible in the time series of the heat flux difference across the two-layer system, as shown in figure 6. Fast oscillations driven by the fluid convection variability are superimposed onto longer time variations owing to flow in the porous layer. Such a time series also illustrates how the two-layer set-up reaches a statistically steady state, the latter being characterised by the flux difference averaging to $0$ over long times. In the particular case of figure 6, the simulation is initiated from the diffusive equilibrium plus a small noise and we observe the equilibration to occur after ${\sim }0.2\tau _{\mathit {diff}}$. Although the equilibration time is largely reduced by the use of continuation, we run all simulations over times that are similar to $\tau _{\mathit {diff}}$ to ensure they are converged.

Figure 6. (a) Time series of the flux difference between top and bottom at $\mathit {Ra} = 10^{7.75}$, $\mathit {Da} = 10^{-5.5}$. The flux difference is normalised by the volume and temporal average of the flux $\langle J \rangle$ and time is normalised by $\tau _{\mathit {diff}} = \sqrt {\mathit {Ra}}$. The insets show zooms to make clearer the different levels of temporal variations, with horizontal segments indicating the porous turnover time scale $\tau _p$ and the free-fall time scale. Note that this simulation has been started from a diffusive state with a small noise.

4.1. Predicting heat transport from individual layer behaviour

As observed in figure 4(a), when both layers are convecting, the temperature profiles in each layer are characterised by thin boundary layers at either the upper or lower boundaries of the domain, through which heat must diffuse. This structure is a generic feature of convection problems, and suggests that we may be able to generalise previous results and approaches used in standard Rayleigh–Bénard convection problems to predict the behaviour here.

4.1.1. An asymptotic approach based on boundary-layer marginal stability

Following the seminal approach of Malkus (Reference Malkus1954) and Howard (Reference Howard1966), we posit that in the asymptotic limit of large Rayleigh and small Darcy numbers such that $\sqrt {\mathit {Ra}}\mathit {Da} \ll 1$, the thermal boundary layers at the upper and lower boundaries of the domain are held at a thickness that is marginally stable to convection. To apply this idea, let us consider a general region with a Rayleigh number $R$ (i.e. ${R} = \mathit {Ra}$ in the fluid layer and ${R} = \mathit {Ra} \mathit {Da}$ in the porous layer). If the boundary layer has mean thickness $\delta$, then we can also introduce a local boundary-layer Rayleigh number $R \delta ^\beta$, where $\beta =1$ in the porous layer or $\beta =3$ in the fluid layer to account for the different dependence on the height scale in the corresponding Rayleigh number (see (2.14c) and (2.15)). Let the temperature difference across the boundary layer be $\varTheta$ (figure 7), and, for completeness, suppose that the region has an arbitrary dimensionless height $\hat {h}$.

Figure 7. A schematic of the temporal and horizontal average of the temperature in the two-layer system, with a focus on the fluid and porous boundary layers.

Suppose further that we want to rescale lengths and temperatures to compare this general case more directly with the standard Rayleigh–Bénard cell of unit dimensionless height and unit temperature difference. In such a cell, the temperature contrast across the boundary layers is $1/2$, which suggests that we need to rescale temperatures by $2 \varTheta$ and lengths by $\hat {h}$, to give a new Rayleigh number $\hat {R} = 2\varTheta \hat {h} {R}$ and boundary-layer depth $\hat {\delta } = \delta /\hat {h}$. Having rescaled thus, the Malkus–Howard approach would suggest that

(4.1)\begin{equation} \hat{R} \hat{\delta}^\beta = {R}_c, \end{equation}

for some critical value ${R}_c$, or $\hat \delta = ({R}_c/\hat {R})^{1/\beta }$. The corresponding Nusselt number ${N}$ for this rescaled system, given that the scaled temperature drop across the layer is $1/2$, is

(4.2)\begin{equation} {N} = \frac{1}{2\delta} = \frac{1}{2} \left(\frac{\hat{R}}{{R}_c}\right)^{1/\beta}. \end{equation}

Note that the actual, unscaled flux $\langle J\rangle$ across the boundary layer is $\langle J \rangle = \mathit {Ra}^{-1/2} \varTheta /\delta$, which can thus be related to ${N}$ from (4.1) and (4.2) by

(4.3)\begin{equation} \langle J \rangle = \frac{ 2 \varTheta}{\hat{h} \sqrt{\mathit{Ra}}} {N}. \end{equation}

Thus, specification of the critical value ${R}_c$ yields a prediction of the flux through the system in terms of $\varTheta$ and $\hat {h}$. We can extract values of ${R}_c$ from previous works that have determined experimentally or numerically the relation between ${N}$ and $\hat {R}$ in either a fluid or a porous Rayleigh–Bénard cell. For porous convection, Hewitt et al. (Reference Hewitt, Neufeld and Lister2012) found that $(2 {R}_c^{1/\beta })^{-1} \simeq 6.9 \times 10^{-3}$. For unconfined fluid convection, the host of historical experiments reported by Ahlers, Grossmann & Lohse (Reference Ahlers, Grossmann and Lohse2009) and Plumley & Julien (Reference Plumley and Julien2019), as well as more recent studies for instance by Urban et al. (Reference Urban, Hanzelka, Musilová, Králík, Mantia, Srnka and Skrbek2014) or Cheng et al. (Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015), suggest that $(2 {R}_c^{1/\beta })^{-1} \sim 6\text {--}8 \times 10^{-2}$, although no definitive observation of a well-developed $\mathit {Ra}^{1/3}$ law has been made so far and the ‘true’ asymptotic form of ${N}(\hat {R})$ remains a hotly contested question. Nevertheless, these values provide an estimate for the heat flux in both the porous and the fluid layer in the asymptotic limit $\mathit {Ra} \gg 1$ and $\sqrt {\mathit {Ra}}\mathit {Da} \ll 1$ that can be compared with our simulations, as discussed in the next section.

4.1.2. Generalising the flux estimate using Rayleigh–Nusselt laws

Our simulations remain limited to moderate Rayleigh numbers, mainly because of the flows through the interface that need to be accurately resolved, and so the asymptotic arguments outlined in the previous section may not be accurate. However, it is straightforward to generalise the asymptotic approach of the previous section to a case where the Nusselt number is some general function $\mathcal {N}(\hat {R})$ of the rescaled Rayleigh number, rather than the asymptotic scaling. To achieve this, we simply replace (4.2) by the relationship

(4.4)\begin{equation} {N} = \mathcal{N}(\hat{R}). \end{equation}

(Equivalently, one could generalise the asymptotic results above to allow $R_c$ to vary with $\hat {R}$ in the manner necessary to recover (4.4).) Over the range of fluid Rayleigh numbers considered here, Cheng et al. (Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015) determined an approximate fit to the fluid Nusselt number function $\mathcal {N}_f$ of

(4.5)\begin{equation} \mathcal{N}_f (\hat{R}_f) = 0.162\, \hat{R}_f^{0.284}. \end{equation}

In the porous layer, Hewitt et al. (Reference Hewitt, Neufeld and Lister2012) and Hewitt (Reference Hewitt2020) considered the equivalent porous Nusselt number $\mathcal {N}_p$ and found a correction to the asymptotic relationship of the form $\mathcal {N}_p (\hat {R}_p) = 6.9 \times 10^ {-3}\, \hat {R}_p + 2.75$, which recovers the asymptotic relationship in the limit $\hat {R}_p \rightarrow \infty$. In fact, we also carried out additional simulations of porous convection in a layer of height $\hat {h} = 1$ submitted to a maintained temperature difference of $1$, but in which the top boundary is open as it allows the fluid to flow in and out. An affine fit of the Nusselt number against the porous Rayleigh number in these simulations provided us with the law

(4.6)\begin{equation} \mathcal{N}_p (\hat{R}_p) = 6.99 \times 10^ {{-}3}\, \hat{R}_p + 1.56. \end{equation}

Note that in such a cell, the temperature difference across the bottom boundary layer is not $1/2$ as in a classical Rayleigh–Bénard cell, but rather approximately $\varTheta \sim 0.8$. Note also that the linear coefficient of the fit in (4.6) is effectively the same as that found in the classical cell, which indicates that for a sufficiently large porous Rayleigh number, the mechanisms controlling the lower boundary layer are the same regardless of the nature of the interface condition.

Figure 8 shows a comparison of the predictions of this theory and various flux laws with our numerical data. We show the Nusselt number ${N}$, calculated in our simulations using (4.3), as a function of the rescaled Rayleigh number $\hat {R}$ in both the fluid layer (index $f$) and in the porous layer (index $p$), with the temperature drop across boundary layers $\varTheta$ determined from numerical simulations. We have indicated for reference the single-layer laws (4.5) and (4.6), which are in close agreement with our data. We also include two simulations with different layers depths $\hat {h}_p$ and $\hat {h}_f$ to demonstrate the generality of the theory. In figure 8(a), we observe that the points lie slightly below the law (4.5) found by Cheng et al. (Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015), but we note that our values of ${N}_f$ lie within the range of the scatter in the data upon which that fit is based in the paper by Cheng et al. (Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015). Overall, the figure indicates that flux laws extracted from individual layers can be used to predict the flux in the two-layer system, after careful rescaling.

Figure 8. The Nusselt number $N$, extracted from the measured flux $\langle J\rangle$ and temperature contrast across the boundary layers via (4.3), as a function of the rescaled Rayleigh number $\hat {R} = {R}\,2 \varTheta \hat {h}^\beta$ in the fluid layer (panel a with subscript $f$) and in the porous layer (panel b with subscript $p$). The values of $\varTheta$, the temperature difference across the boundary layers, are extracted from numerical data. Numerical data are sorted by equal values of the damping factor $\sqrt {\mathit {Ra}}\mathit {Da}$. In each panel, fits from studies of the single-layer case (4.5) or (4.6) are included (dashed grey) together with asymptotic predictions for $\hat {R} \to \infty$ as discussed in § 4.1.1 (red dots). Panel (b) also includes data from simulations of a single porous Rayleigh–Bénard cell with an open-top boundary on which the fit (4.6) is based (red stars). All but two simulations had $\hat {h}_p = \hat {h}_f = 1$; the values of $\hat {h}_p$ for these two simulations are given in the legend.

4.2. The interface temperature

We can use the laws governing the heat flux to determine the interface temperature $\theta _i$, which is important for applications as it gives an order of magnitude of the average fluid temperature. The interface temperature must be set by the constraint of flux conservation between the two layers, which implies that the flux in (4.3) must be the same in both porous and fluid layers, in a statistically steady state. The caveat is that it is not clear how to relate the temperature difference $\varTheta$ across the boundary layers to the interface temperature $\theta _i$. We overcome this issue by saying that the temperature difference across the boundary layer is a fraction of the difference across the whole layer, and so introduce $O(1)$ coefficients $\mathcal {C}_{f,p}$ such that $\varTheta _{f} = \mathcal {C}_f \theta _i$ and $\varTheta _p = \mathcal {C}_p (1-\theta _i)$ (see figure 7). In single classical Rayleigh–Bénard cells, $\mathcal {C}_{f,p} = 0.5$ because both boundary layers are symmetric and diffusive. This is not true in the two-layer system where the transport across the interface includes contributions from diffusion and advection, and $\mathcal {C}_{f,p}$ may take any value ranging from $0.5$ to $1$. Given these coefficients, flux conservation between the layers yields

(4.7)\begin{equation} \frac{2 \mathcal{C}_f \theta_i}{ \hat{h}_f} \mathcal{N}_f (2 \mathcal{C}_f \theta_i \hat{h}_f^3 \mathit{Ra}) = \frac{2 \mathcal{C}_p (1-\theta_i)}{\hat{h}_p} \mathcal{N}_p (2 \mathcal{C}_p (1-\theta_i) \hat{h}_p \mathit{Ra} \mathit{Da}), \end{equation}

which determines the interface temperature $\theta _i$. In general, the fractions $\mathcal {C}_{f,p}$ depend on the control parameters $\mathit {Ra}$ and $\mathit {Da}$, as shown in figure 9.

Figure 9. (a,b) Values of $\mathcal {C}_p$ and $\mathcal {C}_f$ (see figure 7) extracted from simulations for a different porous velocity scale $\sqrt {\mathit {Ra}}\mathit {Da}$. (f) Schematic representation of two limiting temperature profiles in the limit of large Rayleigh numbers based on panel (a) and (b). The porous Rayleigh number is large so that $\mathcal {C}_p = 0.85$ whereas $\mathcal {C}_f$ varies from $1$ in the weakly confined limit ($\sqrt {\mathit {Ra}}\mathit {Da} \lesssim 1$) to $0.5$ in the strongly confined limit ($\sqrt {\mathit {Ra}}\mathit {Da} \ll 1$).

4.2.1. The porous-convective regime

We first address the case where both layers are unstable to convection. Although $\mathcal {C}_p$ and $\mathcal {C}_f$ may take any value between $0.5$ and $1$, we can still make qualitative predictions on their values depending on the control parameters. Because these coefficients describe the temperature difference between the bulk and the boundary layers, they are also a proxy for the temperature drop at the interface. In the case where the porous turnover time scale $\tau _p = \sqrt {\mathit {Ra}}\mathit {Da}$ is very small compared with the free-fall time scale ($=$1 in dimensionless units), the porous flow is strongly confined and heat transfer is purely diffusive at the interface so that $\mathcal {C}_p = \mathcal {C}_f = 0.5$. In the opposite, weakly confined limit where $\sqrt {\mathit {Ra}}\mathit {Da}$ is brought up to 1, the porous and fluid velocities become similar and the interface heat transfer becomes advective, which forces the interface temperature drop to vanish (see, for instance, the $\mathit {Ra} = 10^9$ temperature profile in figure 4(a) for which $\sqrt {\mathit {Ra}}\mathit {Da}=10^{-1}$). In such a limit, $\mathcal {C}_p = \mathcal {C}_f = 1$. The values of $\mathcal {C}_p$ and $\mathcal {C}_f$ extracted from numerical simulations are shown in figures 9(a) and 9(b), respectively. The qualitative picture described above seems to be in agreement with the numerical observation in the fluid coefficient: as shown in figure 9(b), although $\mathcal {C}_f$ varies with $\mathit {Ra}$, it is also larger for larger values of $\sqrt {\mathit {Ra}}\mathit {Da}$. Given there is not a clear collapse of this data, we consider in the following discussion on the interface temperature two extreme limits in the asymptotic limit $\mathit {Ra} \to \infty$, as sketched in figure 9(c), $\mathcal {C}_f = 0.5$ in the strongly confined limit ($\sqrt {\mathit {Ra}}\mathit {Da} \ll 1$) and $\mathcal {C}_f = 1$ in the weakly confined limit ($\sqrt {\mathit {Ra}}\mathit {Da} \lesssim 1$).

In the porous layer, the coefficient $\mathcal{C}_p$ is found to be mainly a function of $Ra_p = RaDa$ and appears to be reaching an asymptote $\mathcal{C}_p\simeq 0.85$ (see figure 9b). $\mathcal{C}_p$ follows a trend that is similar to the case of a single porous layer with an open-top boundary (the red stars in figure 9b), and only weakly depends on $\sqrt{Ra} Da$. This similarity with the single-layer case suggests that the turnover and mixing in the fluid is sufficiently rapid (compared to the motion in the porous medium) so that, from the point of view of the porous layer, the upper layer behaves like a reservoir of fluid at an approximately constant temperature. Given this, some departure from the observed insensitivity to $\sqrt{Ra}Da$ is expected when $\sqrt{Ra} Da \rightarrow 1$; indeed, some suggestion of this might be be visible at $\sqrt{Ra}Da = 10^{-1}$ (see figure 9b). However, despite following similar trends the values of $\mathcal{C}_p$ appear to be larger in the two-layer system than they are in the single porous medium where the asymptotic value is somewhat lower ($\mathcal{C}_p \simeq 0.78$). This is presumably a consequence of the temperature being imposed as a constant at every point on the upper boundary in the case of a single porous layer. This is a stronger constraint than in the two-layer system where hot porous plumes can locally heat the bottom of the fluid layer, hence reducing the interface temperature drop and resulting into larger values of $\mathcal{C}_p$.

Before showing the asymptotic predictions based on these limiting values of $\mathcal {C}_{p,f}$, we first use the actual extracted values to verify the accuracy of (4.7) in predicting the interface temperature (figure 10). We find the agreement between both to be within $10$ % relative error. This figure also shows the predicted interface temperature for $\mathit {Ra} \gg 1$ for the two limiting cases $\mathcal {C}_f = 0.5$ and $\mathcal {C}_f = 1.0$ with $\mathcal {C}_p = 0.85$ as discussed above. Apart from values close to the threshold of porous convection, all the numerical data lies between them. Moreover, the weakly confined simulations with larger $\sqrt {\mathit {Ra}}\mathit {Da}$ show interface temperatures approaching the limit $\mathcal {C}_f =1$ (figure 10a), while the more strongly confined simulations with $\sqrt {\mathit {Ra}}\mathit {Da} \ll 1$ draw closer to $\mathcal {C}_f = 0.5$ as $\mathit {Ra}$ is increased (figure 10d), which suggests that these limits will become increasingly accurate for increasingly large $Ra$.

Figure 10. Plot of the interface temperature $\theta _i$ as a function of the Rayleigh number $\mathit {Ra}$ for different values of $\sqrt {\mathit {Ra}}\mathit {Da}$. The filled symbols are the numerical data and the empty symbol represent is the predicted $\theta _i$ using (4.7) with the values of $\mathcal {C}_{p,f}$ extracted from the simulation and shown in figure 9(a,b). The black lines show the law expected for the interface temperature in the two limiting cases of strongly confined porous medium ($\sqrt {\mathit {Ra}}\mathit {Da} \rightarrow 0$, $\mathcal {C}_p = 0.85$ and $\mathcal {C}_f = 0.5$) and the weakly confined case ($\sqrt {\mathit {Ra}}\mathit {Da} \lesssim 1$, $\mathcal {C}_p = 0.85$ and $\mathcal {C}_f = 1.0$). The red line is the interface temperature in the case where the porous medium remains diffusive, the plus sign highlights the value at threshold of porous convection based on the prediction (4.9).

4.2.2. The porous-diffusive regime

If the porous layer remains stable, dominated by diffusive heat transfer, then the expression for the flux in the fluid layer is still given by (4.3), but in the porous layer, it is simply $\langle J \rangle = \mathit {Ra}^{-1/2} (1-\theta _i)/\hat {h}_p$. A balance of these expressions should then replace (4.7) to determine $\theta _i$. Moreover, we know a priori that $\mathcal {C}_{f} = 0.5$ because the flux is entirely diffusive at the interface $z=0$. Flux conservation thus reduces to

(4.8)\begin{equation} \frac{\theta_i}{\hat{h}_f} \mathcal{N}_f(\theta_i \hat{h}_f^3\mathit{Ra} ) = \frac{1-\theta_i}{\hat{h}_p}, \end{equation}

which we solve numerically. The predicted interface temperature from this model is shown in figure 4(b) and gives reasonably good agreement with the numerical data. Given $\theta _i$, we can also extract the Nusselt number for the two-layer system, ${\mathit {Nu} = 2 \sqrt {\mathit {Ra}} \langle J\rangle = 2 (1-\theta _i)/ \hat {h}_p}$ (see § 3), which was also shown to give good agreement with numerical data in figure 5(a) for the particular case $\hat {h}_p = 1$. Note that this prediction also agrees with an earlier result of Chen & Chen (Reference Chen and Chen1992) that the two-layer Nusselt number is bounded above by $2/\hat {h}_p$ when the porous layer remains stable.

In fact, prediction of $\theta _i$ allows for a more accurate prediction of the threshold of convection in the porous medium. Neglecting any temperature variations induced by fluid convection above, we expect that porous convection sets in when $\mathit {Ra} \mathit {Da}\, \hat {h}_p (1-\theta _i)$ reaches a critical value $\mathit {Ra}_p^c \simeq 27.1$ (Nield & Bejan Reference Nield and Bejan2017). Because $\theta _i \rightarrow 0$ at large Rayleigh numbers while the porous layer remains stable, this condition may be recast as

(4.9)\begin{equation} \mathit{Ra} = \frac{\mathit{Ra}_p^c}{\mathit{Da} \,\hat{h}_p} [ 1+ \theta_i + O(\theta_i^2)] \end{equation}

at the threshold of porous convection. To leading order, the threshold is given by $\mathit {Ra} = \mathit {Ra}_p^c /(\mathit {Da} \hat {h}_p)$, as already noted in § 3.1. For the value of $\mathit {Da}$ used in § 3 where $\hat {h}_p = 1$, including the first-order correction increases the predicted threshold value of $\mathit {Ra}$ by approximately $10\,\%$.

4.2.3. Asymptotic predictions

These relationships for $\theta _i$, and thus for the flux, simplify in the asymptotic regime of very large Rayleigh numbers and extremely low Darcy numbers that is relevant for geophysical applications ($\mathit {Ra} \mathit {Da} \gg 1$ while $\sqrt {\mathit {Ra}}\mathit {Da} \ll 1$). In such a regime, the Rayleigh–Nusselt relations in convecting sub-layers follow the asymptotic scalings discussed at the end of § 4.1.1:

(4.10a,b)\begin{equation} \mathcal{N}_p = 6.9 \times 10^{{-}3} \hat{R}_p \equiv \alpha_p \hat{R}_p \quad \text{and}\quad \mathcal{N}_f \simeq 7 \times 10^{{-}2} \hat{R}_f^{1/3} \equiv \alpha_f \hat{R}_f^{1/3}. \end{equation}

In the case of a stable porous medium, the flux balance (4.8) reduces to

(4.11)\begin{equation} \alpha_f \mathit{Ra}^{1/3} \theta_i^{4/3} = \frac{1-\theta_i}{\hat{h}_p }, \end{equation}

with an asymptotic upper bound

(4.12)\begin{equation} \theta_i \sim \frac{7.3}{\hat{h}_p \,\mathit{Ra}^{1/4}}, \end{equation}

for $\mathit {Ra} \gg 1$. In the case where both layers are convecting, (4.7) instead reduces to

(4.13)\begin{equation} \theta_i^{4/3} = \frac{(2 \mathcal{C}_p)^{2}}{(2 \mathcal{C}_f)^{4/3}} \frac{\alpha_p}{\alpha_f} \, \mathit{Ra}^{2/3} \mathit{Da} \,(1-\theta_i)^2, \end{equation}

which shows that the interface temperature is a function of the grouping $\mathit {Ra}^{2/3} \mathit {Da}$ alone in this limit. The height of the layers $\hat {h}_p$ and $\hat {h}_f$ do not appear in (4.13) because the heat flux is controlled by boundary layers whose widths are asymptotically independent of the depth of the domain (Priestley Reference Priestley1954).

The predictions of (4.13) are shown in figure 11 for the two limiting values of $\mathcal {C}_f$, along with the simulation data for reference. There is rough agreement between data and prediction, although the prediction slightly underestimates the data, presumably owing to the finite values of $\mathit {Ra}$ and $\mathit {Ra}_p$ in the simulations. The figure demonstrates that $\theta _i$ decreases with decreasing $\mathit {Ra}^{2/3} \mathit {Da}$: at constant Rayleigh number, decreasing the Darcy number decreases the efficiency of heat transport in the porous medium and, as a consequence, the porous layer must absorb most of the temperature difference, forcing a decrease in $\theta _i$. In the asymptotic limit of large $\mathit {Ra}$ but very small $\mathit {Da}$, such that $\mathit {Ra}^{2/3} \mathit {Da}$ remains small, the interface temperature from (4.12) satisfies

(4.14)\begin{equation} \theta_i \sim \sqrt{\mathit{Ra} \mathit{Da}^{3/2}}, \end{equation}

for $\mathit {Ra}^{2/3} \mathit {Da} \ll 1$. Conversely, $\theta _i$ increases with increasing $\mathit {Ra}^{2/3} \mathit {Da}$ but evidently cannot increase without bound, which reflects the fact that the grouping $\mathit {Ra}^{2/3}\mathit {Da}$ cannot take arbitrarily large values. For $\mathit {Ra}^{2/3}\mathit {Da} \gtrsim O(10)$, the flow structures in the porous medium will become smaller than the pore scale, breaking the assumption of Darcy flow there (see (2.17)).

Figure 11. Plot of the interface temperature $\theta _i$ found in simulations (symbols) and predicted with the asymptotic model (4.13) as a function of $\mathit {Ra}^{2/3} \mathit {Da}$. The two lines materialise the two limiting cases of a strongly confined porous medium ($\sqrt {\mathit {Ra}}\mathit {Da} \rightarrow 0$, $\mathcal {C}_p = 0.85$ and $\mathcal {C}_f = 0.5$) and the weakly confined case ($\sqrt {\mathit {Ra}}\mathit {Da} \lesssim 1$, $\mathcal {C}_p = 0.85$ and $\mathcal {C}_f = 1.0$). Note that $\mathit {Ra}^{2/3} \mathit {Da}$ must remain smaller than $\sim$10 to ensure that the constraint (2.17) is satisfied.

4.3. Penetrative convection in the porous-diffusive regime

We end this discussion of fluxes with a brief consideration of the issue of so-called ‘penetrative convection’ – significant subcritical flow in the porous layer – which has been a contentious subject in some previous studies. For instance, Poulikakos et al. (Reference Poulikakos, Bejan, Selimos and Blake1986) and more recently Bagchi & Kulacki (Reference Bagchi and Kulacki2014) reported numerical simulations of porous flows forced by fluid convection despite the porous Rayleigh number being sub-critical, but Chen & Chen (Reference Chen and Chen1992) found that convection should be confined to the fluid layer only. In the simulations detailed in § 3, we found that while weak flows do exist in the porous-stable case, they do not enhance heat transport compared with diffusion, in opposition to the results of Bagchi & Kulacki (Reference Bagchi and Kulacki2014) (see, for example, their figure 3.3). We argue here that provided the Darcy number is not too large, this is a general result: it is not possible to induce subcritical flow in the porous layer that has a significant impact on the heat flux through the layer, without violating the limitations of the model outlined in § 2.6.

While the lower boundary condition on the porous layer is uniform, $\theta (z=-\hat {h}_p) = 0$, the temperature at the the interface will, in general, display horizontal variations owing to fluid convection above. Because the temperature cannot be smaller than $0$, the amplitude of these horizontal variations is, at most, of the order of the interface temperature $\theta _i$. Any penetrative flow must be driven by these horizontal variations, and so will have an amplitude $w \sim \sqrt {\mathit {Ra}}\mathit {Da} \theta _i$. The heat transported on average by the penetrative flow scales with $w \theta \sim \sqrt {\mathit {Ra}}\mathit {Da} \theta _i^2$, and so its contribution to the Nusselt number $\mathit {Nu}$ is $\sim \mathit {Ra} \mathit {Da} \theta _i^2$. The contribution of the penetrative flow thus becomes significant, relative to the $O(1)$ diffusive flux through the layer, when $\mathit {Ra} \mathit {Da} \theta _i^2 \sim O(1)$. However, according to the upper bound in (4.12), $\mathit {Ra} \mathit {Da} \theta _i^2 \lesssim 50 \sqrt {\mathit {Ra}} \mathit {Da}$ for $\mathit {Ra} \gg 1$. The assumption that penetrative flows do not transport appreciable heat is therefore accurate, as long as $50 \sqrt {\mathit {Ra}}\mathit {Da} \lesssim 1$ or $\mathit {Ra}_p \lesssim 4 \times 10^{-4} \mathit {Da}^{-1}$. Provided $\mathit {Da} < O(10^{-5})$, this value of the porous Rayleigh number is always larger than the critical value $\mathit {Ra}_p^c$ for the onset of convection in the porous layer, and there will be no enhanced penetrative convection in the porous layer at subcritical values of $\mathit {Ra}_p$. Given that in most geophysical systems we expect $\mathit {Da} \ll 10^{-5}$, we conclude that, in general, diffusion controls the heat flow through the medium and penetrative convection plays a negligible role. In general, for $\mathit {Ra} \gg 1$ penetrative convection can only occur if the Darcy number is such that either the constraint $\sqrt {\mathit {Ra}}\mathit {Da} \ll 1$ (2.16) – which enforces that the flow in the porous medium remains confined with negligible inertial effects – or the constraint $\mathit {Ra} \mathit {Da}^{3/2} \lesssim 50$ (2.17) – which enforces that the flow length scales are larger than the pore scale – are violated.

5.1. Heat-flux variations with height

The contrast between imposed and inherent variability of fluid convection is first illustrated by time series of the horizontally averaged heat flux $J(t,z) = \overline {w\theta } - \mathit {Ra}^{-1/2} \bar {\theta }'$ at different heights, shown in figure 12. The flux in the porous layer experiences long-lived bursts of activity that can amount to up to a 50 % increase of the flux compared with its average value, with a duration that is controlled by the porous turnover time scale $\tau _p \sim (\sqrt {\mathit {Ra} \mathit {Da}})^{-1}$. In figure 12, the signature of these long-time variations can be traced up to the top of the fluid layer, where they are superposed on much faster variations in heat flux associated with the turbulent convective dynamics, which evolve on the $O(1)$ free-fall time scale. However, comparison between the time series at $z=0$ and $z=1$ in figure 12 reveals that the typical intensity of the bursts is notably weaker at the top of the fluid layer than at the interface. Therefore, fluid convection is not a perfect conveyor of the long-time, imposed variability, which it partially filters out.

Figure 12. Time series of the $x$-averaged flux $J (z,t) = \overline {w\theta }- \mathit {Ra}^{-1/2} \bar {\theta }'$ at different heights in a simulation carried out at $\mathit {Ra} = 10^8$ and $\mathit {Da} = 10^{-5.5}$. The heat flux is normalised by its averaged value over the whole domain and over time. Note that the scale of the $y$ axis is larger at mid-height in the fluid layer ($z = 0.5$). Time is normalised by the free-fall time scale. The present simulation spans over more than one diffusive time scale because $\tau _{\mathit {diff}} = 10^4$.

5.2. Spectral content of the heat flux at different heights

We use spectral analysis of the heat flux time series to better quantify the inherent and imposed variability of fluid convection and the latter's filtering effect on imposed variability. Figure 13(a) shows the power spectra $\vert \hat {J}(\omega,z)\vert ^2$ of the signals displayed in figure 12. First, the spectral content that is inherent to fluid convection is easily distinguished from that imposed by porous convection. The variability of the flux in the fluid layer ($z=-0.5$ in figure 13a) is almost entirely contained in harmonics smaller than ${\sim }10^{-2}$. The energy of lower harmonics in the fluid layer ($z\geq 0$) closely follows the spectral content in the porous layer, which confirms that it is primarily imposed by the porous flow. Higher harmonics are therefore controlled by fast fluid convection. This is particularly well illustrated by the spectra at the fluid boundaries ($z=0$ & $z=1$) being effectively identical above $\omega = 5 \times 10^{-2}$: we retrieve the top-down symmetry of classical Rayleigh–Bénard convection with imposed uniform temperature at the boundaries for these larger harmonics. The filtering effect of the fluid layer is visible at lower frequencies where the energy decays as $z$ increases. It is further quantified in figure 13(b), where we show the energy at particular frequencies as a function of depth through the fluid layer. We find that at low frequency, the energy decays exponentially with $z$ and with a rate that seems to increase with $\omega$. The spatial decay is lost for higher harmonics, which are driven directly by fluid convection; instead, the energy is maximised in the bulk of the fluid layer and decays at the boundaries.

Figure 13. (a) Power spectrum $\vert \hat {J}(\omega,z) \vert ^2$ of the heat flux time series shown in figure 12 and normalised by $\langle J \rangle ^2$. (b) Spatial variations of the power spectrum $\vert \hat {J}(\omega,z) \vert ^2$ in (a) with depth through the fluid layer, shown for several frequencies, the lower two being typical of porous convection and the larger two typical of fluid convection. (c) Vertical decay rate $r$ of the energy of low frequency harmonics of the heat flux, for $\mathit {Ra} = 10^{8}$ with $\mathit {Da}$ varied by an order of magnitude from $10^{-6}$ to $10^{-5}$ as indicated in the legend. The dotted line gives the power law $\omega ^{0.75}$ for reference.

5.3. Spatial decay rate of low-frequency variability

We attempt to quantify the filtering effect of fluid convection on the long-time variations by systematically measuring the spatial decay rate of the low-frequency harmonics.

By linear fitting of $\ln (\vert \hat {J} (\omega,z)\vert ^2 )$ with $z \in [0,0.6]$, we extract the decay rate $r$ for all frequencies below $\omega = 10^{-2}$, above which temporal variations become partially imposed by fluid rather porous convection. The result of this process is shown in figure 13(c). Despite some spread in the extracted values, they all appear to follow the same trend: below $\omega \simeq 5 \times 10^{-3}$, the energy decay rate roughly increases like $\omega ^{0.75}$, and above this value, the decay rate saturates, presumably because the harmonics are increasingly driven by fluid convection. That is, sufficiently slow variations imposed by the porous layer decay very slowly through the fluid layer, but more rapid variations decay faster. The exponential decay of the flux harmonics with $z$ is reminiscent of the problem of diffusion in a solid submitted to an oscillating temperature boundary condition. In that classical problem, the decay rate has a diffusive scaling, $r\sim \omega ^{1/2}$, for oscillation frequency $\omega$. Here, our results instead suggest that fluid convection acts as a sub-diffusive process on the low-frequency flux variations imposed by the porous medium. A quantitative explanation for such behaviour remains elusive. It does, however, at least seem reasonable that the decay rate should increase with $\omega$. As $\omega \to 0$, the decay rate must vanish, because the average heat flux through the system is conserved with height, but higher frequency variations in the porous layer will lead to localised bursts of plumes that are quickly mixed into the large-scale circulation of the fluid convection; the extra flux momentarily increase the temperature of the fluid but is not transmitted up to the top of the layer. Additional work, possibly with more idealised models of fluid convection submitted to flux variations, is required to fully understand the phenomenology that we have outlined here. Although our results on the filtering effect remain preliminary, we can still predict the typical porous turnover time scale $\tau _p$ needed to ensure that the decay of the variability in the fluid layer is negligible. The power spectrum of the flux at $z=0$ and $z=\hat {h}_f$ is in a ratio $\sim \exp (-r(\omega ) \hat {h}_f)$. The cut-off frequency of such a low-pass filter, reached when $\exp (-r \hat {h}_f) \sim 1/2$, is roughly located at $\omega = \omega _c = 10^{-3}$ in the particular case of $\hat {h}_f = 1$, according to figure 13(c). Therefore, using $r \sim \omega ^{3/4}$, we predict in general that when $1/\tau _p = \sqrt {\mathit {Ra}}\mathit {Da}$ is smaller than ${\sim }\hat {h}_f^{-4/3} \omega _c$, the temporal variability imposed by porous convection is sufficiently slow that it will be entirely transmitted across the fluid layer.