High-Rayleigh-number convection in porous–fluid layers

Abstract

We present a numerical study of convection in a horizontal layer comprising a fluid-saturated porous bed overlain by an unconfined fluid layer. Convection is driven by a vertical, destabilising temperature difference applied across the whole system, as in the canonical Rayleigh–Bénard problem. Numerical simulations are carried out using a single-domain formulation of the two-layer problem based on the Darcy–Brinkman equations. We explore the dynamics and heat flux through the system in the limit of large Rayleigh number, but small Darcy number, such that the flow exhibits vigorous convection in both the porous and the unconfined fluid regions, while the porous flow still remains strongly confined and governed by Darcy's law. We demonstrate that the heat flux and average thermal structure of the system can be predicted using previous results of convection in individual fluid or porous layers. We revisit a controversy about the role of subcritical ‘penetrative convection’ in the porous medium, and confirm that such induced flow does not contribute to the heat flux through the system. Lastly, we briefly study the temporal coupling between the two layers and find that the turbulent fluid convection above acts as a low-pass filter on the longer time-scale variability of convection in the porous layer.

1. Introduction

Heat transfer driven by flow exchange between a fluid-saturated porous bed and an overlying unconfined fluid arises in a variety of systems in engineering and geophysics. This is the case, for example, in various industrial cooling systems found in nuclear power generation, microelectronics or chemical engineering that require the circulation of fluid from an open channel into a fragmented medium (d'Hueppe et al. 2012; Chandesris et al. 2013; Su, Wade & Davidson 2015). A similar situation occurs during the progressive solidification of multi-component fluids, which creates a mushy solid through which liquid flows to transport heat and solute (Worster 1997). In geophysical contexts, this phenomenon is encountered below sea ice (Wells, Hitchen & Parkinson 2019) and in the Earth's core (Huguet et al. 2016). This work is particularly motivated by the physics of hydrothermal circulation, where a water-saturated, porous bed that is heated from below exhibits thermal convection that, in turn, drives buoyant plumes and convective motion in the overlying ocean. As well as being a well-known feature of the Earth's ocean, there is evidence of on-going hydrothermal activity under the ice crust of icy satellites such as Jupiter's moon Europa (Goodman 2004) or Saturn's moon Enceladus (Hsu et al. 2015). Unlike on Earth, the entire heat budget of these bodies is believed to be controlled by hydrothermal convection and, in particular, by the manner in which heat is transported through the rocky core and into the overlying oceans beneath their icy crusts (Travis, Palguta & Schubert 2012; Travis & Schubert 2015; Nimmo et al. 2018). Most previous works in this area have focussed on either the flow in the porous medium alone or on that in unconfined fluid alone, with the coupling between them modelled by a parametrised boundary condition. This is particularly the case for hydrothermal activity, for which there are numerous studies focussed either on the structure of the flow in the porous layer (see for instance Fontaine & Wilcock 2007; Coumou, Driesner & Heinrich 2008; Choblet et al. 2017; Le Reun & Hewitt 2020, among others), or on the buoyant plumes created in the ocean (Goodman 2004; Woods 2010; Goodman & Lenferink 2012), or on the induced large-scale oceanic circulation (Soderlund et al. 2014; Soderlund 2019; Amit et al. 2020). Travis et al. (2012) included both layers, but resorted to an enhanced diffusivity to parametrise flows in the sub-surface ocean to make calculations tractable. In all these cases, questions remain about how reasonable it is to use these parametrised boundary conditions rather than resolve both layers, and about how the dynamics of flow in each layer communicates and influences the flow in the other layer. Addressing these questions is the focus of this paper. Works that explicitly focus on the coupled transport across a porous and fluid layer are more numerous in engineering settings. However, they tend to either be focused on situations where inertial effects in the interstices of the porous layer play an important role (d'Hueppe et al. 2011; Chandesris et al. 2013), or on regimes where heat is mainly transported by diffusion through the porous layer (Poulikakos et al. 1986; Chen & Chen 1988, 1992; Bagchi & Kulacki 2014; Su et al. 2015) or where restricted to the onset of convective instabilities (Chen & Chen 1988; Hirata et al. 2007, 2009b; Hirata, Goyeau & Gobin 2009a). In general, these studies are difficult to apply to geophysical contexts and particularly to hydrothermal circulation: here, the large spatial and temperature scales and the typically relatively low permeabilities are such that the porous region can be unstable to strong convection while the flow through it remains inertia-free and well described by Darcy's law. In such a situation, there can be vastly different time scales of motion between the unconfined fluid, which exhibits rapid turbulent convection, and the porous medium, where the convective flow through the narrow pores is much slower. This discrepancy in time scales presents a challenge for numerical modelling, which is perhaps why this limit has not been explored until now.

In this paper, we explore thermal convection in a two-layer system comprising a porous bed overlain by an unconfined fluid. In particular, we focus on situations in which the driving density difference, as described by the dimensionless Rayleigh number, is large and heat is transported through both layers by convective flow, although for completeness we include cases in which there is no convection in the porous layer. The permeability of the porous medium, as described by the dimensionless Darcy number, is small enough that the flow through the medium is always inertia free and controlled by Darcy's law. As in some previous studies of this set-up (Poulikakos et al. 1986; Chen & Chen 1988; Bagchi & Kulacki 2014), we consider the simplest idealised system in which natural convection occurs, that is, a two-layer Rayleigh–Bénard cell. In this set-up, the base of the porous medium is heated and the upper surface of the unconfined fluid layer is cooled to provide a fixed destabilising density difference across the domain. Flow in such a system attains a statistically steady state, which allows for investigation of the fluxes, temperature profiles and dynamics of the flow in each layer, and of the coupling between the layers.

We carry out numerical simulations of this problem in two dimensions using a single-domain formulation of the two-layer problem based on the Darcy–Brinkman equations (Le Bars & Worster 2006). Using efficient pseudo-spectral methods, we are able to reach regimes where thermal instabilities are fully developed in both the porous and fluid layers. We demonstrate how to use previous results on thermal convection in individual fluid or porous layers to infer predictions of the heat flux and the temperature at the interface between the layers in our system. In addition, we revisit a long-standing controversy on the role of ‘penetrative convection’, i.e. flow in the porous medium that is actively driven by fluid convection above, and confirm that it is negligible in the limit where the pore scale is small compared with the size of the system. Lastly, we briefly address the temporal coupling between both layers and explore how fluid convection mediates the variability of porous convection.

The paper is organised as follows. The set-up and governing equations are introduced in § 2, where we also outline the main approximations of our model and, importantly, the limits on its validity. After presenting the general behaviour of the two-layer system and how it changes when the porous layer becomes unstable to convection (§ 3), we show how previous works on convection can be used to predict the thermal structure of the flow and the heat it transports (§ 4). In § 5, we discuss the temporal variability of two-layer convection, before summarising our findings and their geophysical implications in § 6.

2. Governing equations and numerical methods

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 2006),

(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 (2006), 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 (1967) 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 (1967). There are, however, some known discrepancies between these two approaches that may not be restricted to the interface (e.g. Le Bars & Worster 2006; Hirata et al. 2007, 2009b). As pointed out by Nield & Bejan (2017), 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 2013)

(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. 2015), 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 2012) 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. 2020; Hester, Vasil & Burns 2021). 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 2008): 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 2008). 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. An overview of two-layer convection

In this section, we describe the results of a series of simulations carried out at a fixed Darcy number, $\mathit {Da} = 10^{-5.5}$, and equal layer depths $\hat {h}_p = \hat {h}_f = 1$, but with varying $\mathit {Ra}$ in the range $10^4 \leqslant \mathit {Ra} \leqslant 10^9$. We use these to illustrate the basic features of high-$\mathit {Ra}$ convective flow in the two-layer system.

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 2017). 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. 2012). 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. (2012).

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. 2012). 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. 2012). 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. Modelling heat transport and interface temperature

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 (1954) and Howard (1966), 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. (2012) 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 (2009) and Plumley & Julien (2019), as well as more recent studies for instance by Urban et al. (2014) or Cheng et al. (2015), 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. (2015) 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. (2012) and Hewitt (2020) 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. (2015), 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. (2015). 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 (1992) 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 2017). 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 1954).

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. (1986) and more recently Bagchi & Kulacki (2014) reported numerical simulations of porous flows forced by fluid convection despite the porous Rayleigh number being sub-critical, but Chen & Chen (1992) 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 (2014) (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. Temporal coupling between the layers

In this two-layer set-up, heat is transported through two systems with very different response time scales while carrying the same average heat flux. The dynamics in the unconfined fluid layer must, therefore, exhibit variability on both a slow time scale imposed by the porous layer below and on a rapid time scale inherent to turbulent fluid convection. In turn, as heat is transported to the top of the two-layer cell, fluid convection must mediate and possibly filter the long variations of the porous activity in a manner that remains to be quantified.

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.

6. Conclusion

In this article, we have explored heat transport in a Rayleigh–Bénard cell composed of a fluid-saturated porous bed overlain by an unconfined fluid layer, using numerical simulations and theoretical modelling. The focus of the work has been on the geologically relevant limit of large Rayleigh number, $\mathit {Ra}$, and small Darcy number, $\mathit {Da}$, such that heat is transported through the system by vigorous convection but the flow within the porous medium remains inertia-free and well described by Darcy's law. To the best of our knowledge, this is the first study of porous–fluid two-layer systems in this limit.

Having identified suitable limits on the parameters, we demonstrated that the dynamics and heat flux through the two-layer system are strongly dependent on whether the flow in the porous layer is unstable to convection, which we demonstrated should occur if $\mathit {Ra} \gtrsim 27/(\hat {h}_p\mathit {Da})$. By suitably rescaling, we showed that flux laws from individual convecting fluid or porous layers could be used to predict both the flux through the two-layer system and the mean temperature at the interface between the two layers. In the asymptotic limit of large $\mathit {Ra}$, we find that while flow in the porous layer remains stable ($\mathit {Ra} \mathit {Da}\hat {h}_p \lesssim 27$), the interface temperature $\theta _i$ satisfies $\theta _i \sim \mathit {Ra}^{-1/4}$, whereas if the porous layer becomes unstable to convection, $\theta _i \sim \sqrt {\mathit {Ra} \mathit {Da}^{3/2}}$, provided $\mathit {Da}$ remains sufficiently small that $\mathit {Ra} \mathit {Da}^{3/2} \ll 1$.

We briefly investigated the role played by ‘penetrative convection’ (Bagchi & Kulacki 2014), i.e. subcritical flows in the porous medium driven by fluid convection above. If the fluid above is convecting, then some weak flow will always be driven in the porous layer by horizontal temperature variations imposed at the interface; but we show that they are always too weak to contribute significantly to heat transport, unless some of the model assumptions about flow in the porous medium are violated. Interestingly, the laws governing or limiting heat transport and the interface temperature in porous–fluid convection have been derived from the behaviour of each layer considered separately. They do not rely on the details of the flow at the interface, in particular, on the pore-scale boundary layer at the transition between the porous and the pure fluid flows. Therefore, the laws that we have identified here should only weakly depend on the choice of the formulation of the two-layer convection problem (see the discussion in § 2.1).

Lastly, we also briefly explored the manner in which rapid fluid convection mediates the long-time variations of activity in the porous layer. The amplitude of low-frequency temporal variations of flux imposed by porous convection decay exponentially through the unconfined fluid layer in a way that is reminiscent of a diffusive process. As a result, fluid convection acts as a low-pass filter on the bursts of activity in the porous layer. We predict that for $\sqrt {\mathit {Ra}}\mathit {Da} < 10^{-3} \hat {h}_f^{-4/3}$, the decay of low-frequency variations in the fluid layer are negligible and the variability of porous convection is entirely transmitted to the top of the two-layer system.

Before ending, we return to consider the question of how important it is to resolve both porous and fluid layers when studying these coupled systems in astrophysical or geophysical settings, rather than using parametrised boundary conditions. Our results in the limit of strong convection suggest that, while the details of convection in each layer are important for controlling the interface temperature and heat flux through the system, there is very little coupling in the dynamical structure of the flow between each layer. As may be noticed in figure 3, for example, fluid convection remains organised in large-scale rolls even if it is forced by several hot plumes from the porous medium below. Therefore, a localised ‘hot spot’ at the interface associated with, say, a strong plume in the porous medium, does not, in general, lead to any associated hot spot at the surface of the fluid layer.

This observation helps to justify the approach of various studies which neglect the dynamics of flow in the porous medium altogether (e.g. Soderlund 2019 and Amit et al. 2020 in the context of convection in icy moons). It is also an interesting point in the context of Enceladus, which is well known for sustaining a strong heat-flux anomaly at its South Pole that is believed to arise from hydrothermal circulation driven by tidal heating in its rocky porous core (Spencer et al. 2006; Choblet et al. 2017; Spencer et al. 2018; Le Reun & Hewitt 2020). While models predict that tidal heating does drive a hot plume in the porous core at the poles (Choblet et al. 2017), our results suggest that this hot spot is unlikely to be transmitted through the overlying ocean to the surface without invoking other ingredients such as rotation (Soderlund et al. 2014; Soderlund 2019; Amit et al. 2020) or topography caused by melting at the base of the ice shell (Favier, Purseed & Duchemin 2019).