# High-Rayleigh-number convection in porous–fluid layers

Thomas Le Reun 1<EMAIL_ADDRESS>Duncan R. Hewitt 2  1 DAMTP, University of
Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK 2 Department of
Mathematics, University College London, UK

###### 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-timescale 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 et al.,
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 et al., 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 et al., 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 et al. (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 parameterised
boundary conditions rather than resolve both layers, and about how the
dynamics of flow in each layer communicate and influence 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, 2009a,
2009b).  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 timescales 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 timescales 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 setup (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 setup, 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 the fluid layer. 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 to 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 setup and governing equations are
introduced in section 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 (section 3), we show how previous works
on convection can be used to predict the thermal structure of the flow and the
heat it transports (section 4). In section 5 we discuss the temporal
variability of two-layer convection, before summarizing our findings and their
geophysical implications in section 6.

##  2 Governing equations and numerical methods

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$.

###  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
between $-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$.

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

$\phi(z)=\left\\{\begin{array}[]{cc}\phi_{0}&z<0\\\ 1&z\geq
0\end{array}\right.$ (1)

Given local fluid velocity $\boldsymbol{U}_{\ell}$, the flux
$\boldsymbol{U}=\phi\boldsymbol{U}_{\ell}$ reduces to the usual Darcy flux in
the porous medium, and to the fluid velocity in the unconfined fluid layer.

The flow is assumed to be incompressible everywhere, so

$\bnabla\cdot\boldsymbol{U}=0.$ (2)

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

$\rho\left[\displaystyle\p_{t}\boldsymbol{U}+\boldsymbol{U}\cdot\bnabla\boldsymbol{U}\right]=-\displaystyle\bnabla
P+\mu\bnabla^{2}\boldsymbol{U}+\rho\boldsymbol{g}$ (3)

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

$\boldsymbol{U}=\frac{K}{\mu}\left(-\displaystyle\bnabla
P+\rho\boldsymbol{g}\right).$ (4)

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

$\rho\left[\displaystyle\p_{t}\boldsymbol{U}+\boldsymbol{U}\cdot\bnabla\frac{\boldsymbol{U}}{\phi}\right]=-\phi\displaystyle\bnabla
P+\mu\bnabla^{2}\boldsymbol{U}+\phi\rho\boldsymbol{g}-\frac{\mu}{K}\phi\boldsymbol{U}.$
(5)

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 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 reduces to

$-\phi\displaystyle\bnabla
P+\mu\bnabla^{2}\boldsymbol{U}+\phi\rho\boldsymbol{g}-\frac{\mu}{K}\phi\boldsymbol{U}\sim
0.$ (6)

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 choose to model the two-layer system with a one-domain
formulation via the Darcy-Brinkman equation. Another classical alternative
formulation of the problem is the one introduced after 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 all 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). Yet,
there are known discrepancies between these two approaches that may not be
restricted to the interface, as shown by Le Bars & Worster (2006). In the
specific case of two-layer convection, Hirata et al. (2007, 2009a) have
demonstrated that the use of one-domain or two-domain formulation lead to
sightly different thresholds for the development of convective instabilities
in the two-layer system. As pointed out by Nield & Bejan (2017), there is
still debate on which of these formulations is the most adequate to model
flows in porous-fluid layers, definitive empirical evidence being still
awaiting. Although we choose to use the one-domain formulation of fluid-porous
convection, we will show numerical results that can be understood without
relying on the details of the interface flow were models might be inaccurate.

###  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

$\rho=\rho_{0}\left(1-\alpha(T-T_{0})\right)$ (7)

with $T_{0}$ being a reference temperature. The momentum equation under the
Boussinesq approximation follows from substituting 7 into the buoyancy term of
5, while letting $\rho=\rho_{0}$ in the inertial terms. The temperature
evolves according to an energy transport equation (Nield & Bejan, 2013)

$\overline{\phi}\p_{t}T+\boldsymbol{U}\cdot\bnabla T=\kappa\Lambda\nabla^{2}T$
(8)
$\displaystyle\mbox{with}~{}~{}~{}\overline{\phi}\equiv\frac{(1-\phi)c_{m}\rho_{m}+\phi
c\rho}{\rho c},~{}~{}~{}\kappa\equiv\frac{\lambda}{\rho
c}~{}~{}~{}\mbox{and}~{}~{}~{}\Lambda\equiv\phi+(1-\phi)\frac{\lambda_{m}}{\lambda}$
(9)

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 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–Benard cell (figure 1). Specifically,
for the temperature we set

$T(z=h)=T_{0},\qquad T(z=-h)=T_{0}+\Delta T,$ (10)

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

$U(z=\pm h)=V(z=\pm h)=0.$ (11)

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 5 allows for application of the no-slip condition in 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 6). It is not clear whether such a basal
boundary layer should arise in experimental situations. Nevertheless, 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 section 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

In order to extract the dimensionless equations that govern the two-layer
system we use a “free-fall” normalisation of 5 and 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,

${U^{*}}^{2}=\alpha\Delta Tgh,$ (12)

and the associated free-fall timescale 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

$\displaystyle\p_{t}\boldsymbol{u}+\boldsymbol{u}\cdot\bnabla\frac{\boldsymbol{u}}{\phi}=-\phi\bnabla
p+\sqrt{\frac{\Pr}{\mathrm{Ra}}}\bnabla^{2}\boldsymbol{u}+\phi\theta\boldsymbol{e}_{z}-\frac{\chi(z)}{\mathrm{Da}}\sqrt{\frac{\Pr}{\mathrm{Ra}}}\phi\boldsymbol{u}$
(13a)
$\displaystyle\overline{\phi}\p_{t}\theta+\boldsymbol{u}\cdot\bnabla\theta=\frac{\Lambda}{\sqrt{\mathrm{Ra}\Pr}}\nabla^{2}\theta$
(13b) $\displaystyle\bnabla\cdot\boldsymbol{u}=0,$ (13c)

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:

$\mathrm{Da}\equiv\frac{K_{0}}{h^{2}},~{}~{}~{}\Pr\equiv\frac{\nu}{\kappa}~{}~{}~{}\mbox{and}~{}~{}~{}\mathrm{Ra}\equiv\frac{\alpha
g\Delta Th^{3}}{\nu\kappa}~{}.$ (14)

The Darcy number $\mathrm{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 section 4.

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

$\mathrm{Ra}_{p}=\frac{\mathrm{Ra}\mathrm{Da}}{\Lambda}=\frac{\alpha g\Delta
TK_{0}h}{\nu\Lambda\kappa}~{}.$ (15)

###  2.5 Further simplifying assumptions

We simplify the complexity of 13 by noting that in the bulk of either the
fluid or the porous regions, $\phi$ cancels out of the equations (see for
instance 3 and 4). The porosity only affects 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{\mathrm{Da}}/\phi$ (as shown by 6).
In the following, we thus take $\phi=1$ in 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 section 2.7.

In the heat transport equation 13b, we reduce the number of control parameters
by taking $\overline{\phi}=\Lambda=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 $\overline{\phi}$ and $\Lambda$ 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 $\Lambda$ and thermal diffusion would be greatly
enhanced in the porous medium, reducing 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 $\Lambda$.

Finally, to reduce the complexity of this study and maintain a focus on the
key features of varying the driving buoyancy forces (i.e. $\mathrm{Ra}$) and
the properties of the porous matrix (i.e. $\mathrm{Da}$), we set the Prandtl
number to be $\Pr=1$ throughout this work.

###  2.6 Limits on the control parameters

Several constraints must be imposed on the control parameters $\mathrm{Ra}$
and $\mathrm{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 $\Pr=\Lambda=1$. First,
the inertial terms must vanish in the porous medium, which demands that

$\frac{\sqrt{\mathrm{Ra}}\mathrm{Da}}{\sqrt{\Pr}}\ll 1;$ (16)

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 lengthscale of the flow in the porous layer must be larger than
the pore scale $\sqrt{\mathrm{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 13b: such a balance, given typical velocity
$\sim\sqrt{\mathrm{Ra}}\mathrm{Da}/\sqrt{\Pr}$, yields a lengthscale
$\sim\mathrm{Ra}_{p}^{-1}$. In fact, simulations carried out in a porous
Rayleigh-Bénard cell (Hewitt et al., 2012) indicate that the narrowest
structures of the flow, which are thermal boundary layers, have a thickness of
at least $50\mathrm{Ra}_{p}^{-1}$. For these structures to remain larger than
the pore scale $\sqrt{\mathrm{Da}}$, we thus require

$\mathrm{Ra}\mathrm{Da}^{3/2}\lesssim 50\Lambda.$ (17)

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
13a, which would no longer reduce to Darcy’s law.

Figure 2 provides an overview of the space of control parameters $\mathrm{Ra}$
and $\mathrm{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 section 3.2 and which is
theoretically quantified in section 4.2.2.

Figure 2: Domain of existence and validity of the two-layer model with respect
to the control parameters $\mathrm{Ra}$ and $\mathrm{Da}$, given
$\Pr=\Lambda=1$. Each dot represent a simulation. The line
$\sqrt{\mathrm{Ra}}\mathrm{Da}=1$ marks the limit beyond which inertial terms
affect the flow in the porous medium, while the line
$\mathrm{Ra}\mathrm{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 $\mathrm{Ra}\mathrm{Da}=\mathrm{Ra}_{p}^{c}=27.1$ corresponds
to the threshold of thermal convection in a porous Rayleigh-Bénardcell with an
open-top boundary as discussed in section 3.1.

###  2.7 Numerical method

The one-domain Darcy-Brinkman equations 13 are solved using the pseudo-
spectral code Dedalus (Burns et al., 2020; Hester et al., 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 comprised of two layers whose interface must be accurately
resolved, each layer is discretised with its own Chebyshev grid, with
$[M_{p},M_{f}]$ nodes 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{\mathrm{Da}}$ regularisation length must be
resolved. Time evolution of the fields is computed with implicit-explicit
methods (Wang & Ruuth, 2008): the non-linear and Darcy terms in 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 timestep in the whoe two-layer
domain, with an upper limit given by $\sqrt{\mathrm{Ra}}\mathrm{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 $\mathrm{Ra}$
and $\mathrm{Da}$ for various fixed values of the porous velocity scale
$\sqrt{\mathrm{Ra}}\mathrm{Da}$.

The majority of simulations were carried out in a set-up where the heights of
both the porous and the fluid layer were equal $\hat{h}_{p}=\hat{h}_{f}=1$,
with resolutions $N\times[M_{p},M_{f}]=1024\times[128,256]$ below
$\mathrm{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 $\mathrm{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 timescale
$\sqrt{\mathrm{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, $\mathrm{Da}=10^{-5.5}$, and equal layer depths
$\hat{h}_{p}=\hat{h}_{f}=1$, but with varying $\mathrm{Ra}$ in the range
$10^{4}\leqslant\mathrm{Ra}\leqslant 10^{9}$. We use these to illustrate the
basic features of high-$\mathrm{Ra}$ convective flow in the two-layer system.

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

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 $\mathrm{Da}=10^{-5.5}$. The colour
scale is cut at $0.8$ to enhance the contrast in the fluid layer.

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,
$\overline{\theta}(z)$, are shown in figure 4(a), while the mean interface
temperature $\theta_{i}=\overline{\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 $\overline{\theta}$ in the bulk of the fluid.

The porous layer, however, exhibits two different behaviours depending on the
size of $\mathrm{Ra}$. At low Rayleigh numbers ($\mathrm{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 $\overline{\theta}(z)$
appear to be linear (figure 4a). The corresponding interface temperature
monotonically decreases with $\mathrm{Ra}$ (figure 4b) As the Rayleigh number
is increased beyond $\mathrm{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 $\mathrm{Ra}$ around $\mathrm{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 open-top boundary, convection occurs if the porous
Rayleigh number $\mathrm{Ra}_{p}=\mathrm{Ra}\mathrm{Da}$ exceeds a critical
value $\mathrm{Ra}_{p}^{c}\simeq 27.1$ (Nield & Bejan, 2017).  At
$\mathrm{Da}=10^{-5.5}$, the critical Rayleigh number $\mathrm{Ra}$ such that
$\mathrm{Ra}\mathrm{Da}=\mathrm{Ra}_{p}^{c}$ is $\mathrm{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 section 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 figure 4a) that has also been observed in
previous works on porous convection at large Rayleigh number Hewitt et al.
(2012). Lastly, the temperature contrast across the interface reduces as the
Rayleigh number is increased.

###  3.2 Characteristics of heat transport

Figure 4: Temporally averaged quantities for $\mathrm{Da}=10^{-5.5}$. (a)
Horizontally averaged temperature $\overline{\theta}(z)$ for values of
$\mathrm{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 25 and the black lines are asymptotic
predictions obtained by solving 24 using $\mathcal{C}_{p}=0.85$ and the marked
values $\mathcal{C}_{f}$, as detailed in section section 4.2.  Figure 5:  a:
Nusselt number $\mathrm{Nu}=2\sqrt{\mathrm{Ra}}\left<J\right>$ characterising
heat transport across the two-layer system, together with predictions from
section 4 for the porous-stable case 25 (red dashed) and porous-unstable 24
(black dotted and dot-dashed, with $\mathcal{C}_{p}=0.85$ and
$\mathcal{C}_{f}$ as marked). The horizontal line marks $\mathrm{Nu}=2$. b:
Root-mean-squared vertical velocity amplitude $w_{\mathrm{rms}}$ in the porous
and fluid layers. The two grey lines indicate a scaling of
$\sqrt{\mathrm{Ra}}\mathrm{Da}$ (the characteristic speed in the porous
layer). c: Ratio of the diffusive ($J_{\mathrm{d,p}}$) and advective
($J_{\mathrm{a,p}}$) fluxes to the total depth-averaged flux $\left\langle
J\right\rangle$ in the porous layer. In all three panels,
$\mathrm{Da}=10^{-5.5}$ and the vertical line marks the threshold of porous
convection estimated using 26.

The transition between the porous stable and the porous unstable case 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}-\mathrm{Ra}^{-1/2}\overline{\theta}^{\prime}$,
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 to what it would be in a purely
diffusive system, $\mathrm{Ra}^{-1/2}/2$, via the Nusselt number
$\mathrm{Nu}\equiv 2\sqrt{\mathrm{Ra}}\left<J\right>$, 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 ($\mathrm{Ra}\lesssim 10^{7}$), $\mathrm{Nu}$ appears
to approach a horizontal asymptote $\mathrm{Nu}=2$, but once the porous layer
is unstable $\mathrm{Nu}$ increases much more steeply beyond this value.

The behaviour below the threshold of convection is due to 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, $\left\langle J\right\rangle\to
1/\sqrt{\mathrm{Ra}}$ and $\mathrm{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 section 4.3, and defer more
detailed discussion and prediction of the behaviour of
$\mathrm{Nu}(\mathrm{Ra},\mathrm{Da})$ until the following section.

###  3.3 Timescales, variability and statistically steady state

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

The governing equations 13 reveal three different timescales 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_{\mathrm{diff}}=\sqrt{\mathrm{Ra}}$,
and the third is the turnover timescale in the porous layer
$\tau_{p}\sim(\sqrt{\mathrm{Ra}}\mathrm{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 timescales are in a ratio
$\tau_{p}=\tau_{\mathrm{diff}}/\mathrm{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 16).

The turnover timescales should scale with the inverse of $w_{\mathrm{rms}}$ in
each layer, as can be observed in figure 5(b): in the fluid layer,
$w_{\mathrm{rms}}\sim O(1)$ with no systematic variation with $Ra$, while in
the porous layer, $w_{\mathrm{rms}}\propto\sqrt{\mathrm{Ra}}$ at constant
$\mathrm{Da}$ in agreement with the scaling above. These two very different
timescales are also clearly visible in 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
due to flow in the porous layer. Such a time series also illustrate 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_{\mathrm{diff}}$. Although the equilibration time is largely
reduced by the use of continuation, we run all simulations over times that are
similar to $\tau_{\mathrm{diff}}$ to ensure they are converged.

##  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 in order to predict the behaviour here.

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

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.

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{\mathrm{Ra}}\mathrm{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. In order to apply this idea, let us consider
a general region with a Rayleigh number $R$ (i.e. $\mathrm{R}=\mathrm{Ra}$ in
the fluid layer and $\mathrm{R}=\mathrm{Ra}\mathrm{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 14
and 15). Let the temperature difference across the boundary layer be $\Theta$
(figure 7), and, for completeness, suppose that the region has an arbitrary
dimensionless height $\hat{h}$.

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\Theta$ and lengths by $\hat{h}$, to give
a new Rayleigh number $\hat{\mathrm{R}}=2\Theta\hat{h}\mathrm{R}$ and
boundary-layer depth $\hat{\delta}=\delta/\hat{h}$. Having rescaled thus, the
Malkus–Howard approach would suggest that

$\hat{\mathrm{R}}\hat{\delta}^{\beta}=\mathrm{R}_{c},$ (18)

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

$\mathrm{N}=\frac{1}{2\delta}=\frac{1}{2}\left(\frac{\hat{\mathrm{R}}}{\mathrm{R}_{c}}\right)^{1/\beta}.$
(19)

Note that the actual, unscaled flux $\left<J\right>$ across the boundary layer
is $\left<J\right>=\mathrm{Ra}^{-1/2}\Theta/\delta$, which can thus be related
to $\mathrm{N}$ from 18 and 19 by

$\left<J\right>=\frac{2\Theta}{\hat{h}\sqrt{\mathrm{Ra}}}\mathrm{N}.$ (20)

Thus specification of the critical value $\mathrm{R}_{c}$ yields a prediction
of the flux through the system in terms of $\Theta$ and $\hat{h}$. We can
extract values of $\mathrm{R}_{c}$ from previous works that have determined
experimentally or numerically the relation between $\mathrm{N}$ and
$\hat{\mathrm{R}}$ in either a fluid or a porous Rayleigh-Bénard cell. For
porous convection, Hewitt et al. (2012) found that
$(2\mathrm{R}_{c}^{1/\beta})^{-1}\simeq 6.9\times 10^{-3}$. For unconfined
fluid convection, the host of historical experiments reported in Ahlers et al.
(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\mathrm{R}_{c}^{1/\beta})^{-1}\sim 6-8\times 10^{-2}$, although no
definitive observation of a well-developed $\mathrm{Ra}^{1/3}$ law has been
made so far and the ‘true’ asymptotic form of $\mathrm{N}(\hat{\mathrm{R}})$
remains a hotly contested questions. Nevertheless, these values provide an
estimate for the heat flux in both the porous and the fluid layer in the
asymptotic limit $\mathrm{Ra}\gg 1$ and $\sqrt{\mathrm{Ra}}\mathrm{Da}\ll 1$
that can be compared to 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{\mathrm{R}})$ of the rescaled Rayleigh number, rather than
the asymptotic scaling. To achieve this, we simply replace 19 by the
relationship

$\mathrm{N}=\mathcal{N}(\hat{R}).$ (21)

(Equivalently, one could generalise the asymptotic results above to allow
$R_{c}$ to vary with $\hat{R}$ in the manner necessary to recover 21.) 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

$\mathcal{N}_{f}(\hat{\mathrm{R}}_{f})=0.162\,\hat{\mathrm{R}}_{f}^{0.284}.$
(22)

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{\mathrm{R}}_{p})=6.9\times
10^{-3}\,\hat{\mathrm{R}}_{p}+2.75$, which recovers the asymptotic
relationship in the limit $\hat{\mathrm{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

$\mathcal{N}_{p}(\hat{\mathrm{R}}_{p})=6.99\times
10^{-3}\,\hat{\mathrm{R}}_{p}+1.56~{}.$ (23)

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
about $\Theta\sim 0.8$. Note also that the linear coefficient of the fit in 23
is effectively the same as that found in the classical cell, which indicates
that for 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:  The Nusselt number $N$, extracted from the measured flux
$\left<J\right>$ and temperature contrast across the boundary layers via 20,
as a function of the rescaled Rayleigh number
$\hat{\mathrm{R}}=\mathrm{R}\,2\Theta\hat{h}^{\beta}$ in the fluid layer
(panel a and subscripts $f$) and in the porous layer (panel b and subscripts
$p$).  The values of $\Theta$, the temperature difference across the boundary
layers, is extracted from numerical data.  Numerical data are sorted by equal
values of the damping factor $\sqrt{\mathrm{Ra}}\mathrm{Da}$. In each panel,
fits from studies of the single-layer case 22 or 23 are included (dashed grey)
together with asymptotic predictions for $\hat{\mathrm{R}}\to\infty$ as
discussed in section 4.1.1 (dotted red). Panel b also includes data from
simulations of a single porous Rayleigh-Bénard cell with an open top boundary
on which the fit 23 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.

Figure 8 shows a comparison of the predictions of this theory and various flux
laws with our numerical data. We show the Nusselt number $\mathrm{N}$,
calculated in our simulations using 20, as a function of the rescaled Rayleigh
number $\hat{\mathrm{R}}$ in both the fluid layer (index $f$) and in the
porous layer (index $p$), with the temperature drop across boundary layers
$\Theta$ determined from numerical simulations. We have indicated for
reference the single-layer laws 22 and 23 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 22 found
by Cheng et al. (2015), but we note that our values of $\mathrm{N}_{f}$ lie
within the range of the scatter in the data upon which that fit is based in
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.

###  4.2 The interface temperature

Figure 9:  a & b: Values of $\mathcal{C}_{p}$ and $\mathcal{C}_{f}$ (see
figure 7) extracted from simulations for different porous velocity scale
$\sqrt{\mathrm{Ra}}\mathrm{Da}$. c: 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{\mathrm{Ra}}\mathrm{Da}\lesssim 1$) to $0.5$ in the strongly confined
limit ($\sqrt{\mathrm{Ra}}\mathrm{Da}\ll 1$).

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 20 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 $\Theta$ 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 $\Theta_{f}=\mathcal{C}_{f}\theta_{i}$ and
$\Theta_{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
advecion, and $\mathcal{C}_{f,p}$ may take any value ranging from $0.5$ to
$1$.

Given these coefficients, flux conservation between the layers yields

$\frac{2\mathcal{C}_{f}\theta_{i}}{\hat{h}_{f}}\mathcal{N}_{f}(2\mathcal{C}_{f}\theta_{i}\hat{h}_{f}^{3}\mathrm{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}\mathrm{Ra}\mathrm{Da}),$
(24)

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

####  4.2.1 The porous-convective regime

Figure 10:  Plot of the interface temperature $\theta_{i}$ as a function of
the Rayleigh number $\mathrm{Ra}$ for different values of
$\sqrt{\mathrm{Ra}}\mathrm{Da}$. The filled symbols are the numerical data and
the empty symbol represent is the predicted $\theta_{i}$ using equation 24
with the values of $\mathcal{C}_{p,f}$ extracted from the simulation and shown
in figure 9a & b. The black lines show the law expected for the interface
temperature in the two limiting cases of strongly confined porous medium
($\sqrt{\mathrm{Ra}}\mathrm{Da}\rightarrow 0$, $\mathcal{C}_{p}=0.85$ and
$\mathcal{C}_{f}=0.5$) and the weakly confined case
($\sqrt{\mathrm{Ra}}\mathrm{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 highlighting the
value at threshold of porous convection based on the prediction 26.

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 timescale $\tau_{p}=\sqrt{\mathrm{Ra}}\mathrm{Da}$ is very
small compared to the free-fall timescale ($=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{\mathrm{Ra}}\mathrm{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 $\mathrm{Ra}=10^{9}$ temperature profile in
figure 4a for which $\sqrt{\mathrm{Ra}}\mathrm{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 figure 9a and b, respectively. The qualitative
picture described above seems in agreement with the numerical observation in
the fluid coefficient: as shown in figure 9b, although $\mathcal{C}_{f}$
varies with $\mathrm{Ra}$ it is also larger for larger values of
$\sqrt{\mathrm{Ra}}\mathrm{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 $\mathrm{Ra}\to\infty$, as sketched in
figure 9(c), $\mathcal{C}_{f}=0.5$ in the strongly confined limit
($\sqrt{\mathrm{Ra}}\mathrm{Da}\ll 1$ ) and $\mathcal{C}_{f}=1$ in the weakly
confined limit ($\sqrt{\mathrm{Ra}}\mathrm{Da}\lesssim 1$).

In the porous layer, the numerical data shows difference with the qualitative
expactations for the values of $\mathcal{C}_{p}$ (see figure 9a).
$\mathcal{C}_{p}$ is mainly a function of the porous Rayleigh number
$\mathrm{Ra}_{p}=\mathrm{Ra}\mathrm{Da}$ (figure 9a), and appears to approach
an asymptote $\mathcal{C}_{p}\simeq 0.85$. Interestingly, the trend of
$\mathcal{C}_{p}$ with $\mathrm{Ra}_{p}$ is the same as in a single porous
layer with an open top (red stars in figure 9a), but with smaller values and
lower asymptote in that case ($\simeq 0.78$). In the confined porous limit,
this is due to the temperature drop at the interface being shared by both the
fluid and the porous layer. In the weakly confined limit where the interface
temperature drip tends to vanish, the value of $\mathcal{C}_{p}\simeq 0.85$
rather reflects the temperature decrease in the bulk of the porous layer (see
the $\mathrm{Ra}=10^{9}$ temperature profile in figure 4a).

Before showing asymptotic predictions based on these limiting values of
$\mathcal{C}_{p,f}$, we first use the actual extracted values to verify the
accuracy of 24 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 $\mathrm{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{\mathrm{Ra}}\mathrm{Da}$ show interface temperatures approaching the
limit $\mathcal{C}_{f}=1$ (figure 10a), while the more strongly confined
simulations with $\sqrt{\mathrm{Ra}}\mathrm{Da}\ll 1$ draw closer to
$\mathcal{C}_{f}=0.5$ as $\mathrm{Ra}$ is increased (figure 10d), which
suggests that these limits will become increasingly accurate for increasingly
large $Ra$.

####  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 20, but in
the porous layer it is simply
$\left<J\right>=\mathrm{Ra}^{-1/2}(1-\theta_{i})/\hat{h}_{p}$. A balance of
these expressions should then replace 24 to determine $\theta_{i}$. Moreover,
we know a priori that $\mathcal{C}_{f}=0.5$ since the flux is entirely
diffusive at the interface $z=0$. Flux conservation thus reduces to

$\frac{\theta_{i}}{\hat{h}_{f}}\mathcal{N}_{f}(\theta_{i}\hat{h}_{f}^{3}\mathrm{Ra})=\frac{1-\theta_{i}}{\hat{h}_{p}},$
(25)

which we solve numerically. The predicted interface temperature from this
model was 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,
$\mathrm{Nu}=2\sqrt{\mathrm{Ra}}\left<J\right>=2(1-\theta_{i})/\hat{h}_{p}$
(see section 3), which was also shown giving 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 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 $\mathrm{Ra}\mathrm{Da}\,\hat{h}_{p}(1-\theta_{i})$ reaches a
critical value $\mathrm{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

$\mathrm{Ra}=\frac{\mathrm{Ra}_{p}^{c}}{\mathrm{Da}\,\hat{h}_{p}}\left[1+\theta_{i}+O(\theta_{i}^{2})\right]$
(26)

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

####  4.2.3 Asymptotic predictions

Figure 11:  Plot of the interface temperature $\theta_{i}$ found in
simulations (symbols) and predicted with the asymptotic model 30 as a function
of $\mathrm{Ra}^{2/3}\mathrm{Da}$. The two lines materialise the two limiting
cases of strongly confined porous medium
($\sqrt{\mathrm{Ra}}\mathrm{Da}\rightarrow 0$, $\mathcal{C}_{p}=0.85$ and
$\mathcal{C}_{f}=0.5$) and the weakly confined case
($\sqrt{\mathrm{Ra}}\mathrm{Da}\lesssim 1$, $\mathcal{C}_{p}=0.85$ and
$\mathcal{C}_{f}=1.0$). Note that $\mathrm{Ra}^{2/3}\mathrm{Da}$ must remain
smaller thant $\sim 10$ to ensure that the constraint 17 is satisfied.

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
($\mathrm{Ra}\mathrm{Da}\gg 1$ while $\sqrt{\mathrm{Ra}}\mathrm{Da}\ll 1$). In
such a regime, the Rayleigh-Nusselt relations in convecting sub-layers follow
the asymptotic scalings discussed at the end of section 4.1.1:

$\mathcal{N}_{p}=6.9\times
10^{-3}\hat{\mathrm{R}}_{p}\equiv\alpha_{p}\hat{\mathrm{R}}_{p}~{}~{}~{}~{}\mbox{and}~{}~{}~{}~{}\mathcal{N}_{f}\simeq
7\times
10^{-2}\hat{\mathrm{R}}_{f}^{1/3}\equiv\alpha_{f}\hat{\mathrm{R}}_{f}^{1/3}.$
(27)

In the case of a stable porous medium, flux balance 25 reduces to

$\alpha_{f}\mathrm{Ra}^{1/3}\theta_{i}^{4/3}=\frac{1-\theta_{i}}{\hat{h}_{p}}$
(28)

with an asymptotic upper bound

$\theta_{i}\sim\frac{7.3}{\hat{h}_{p}\mathrm{Ra}^{1/4}}$ (29)

for $\mathrm{Ra}\gg 1$. In the case where both layers are convecting, 24
instead reduces to

$\theta_{i}^{4/3}=\frac{(2\mathcal{C}_{p})^{2}}{(2\mathcal{C}_{f})^{4/3}}\frac{\alpha_{p}}{\alpha_{f}}\,\mathrm{Ra}^{2/3}\mathrm{Da}\,(1-\theta_{i})^{2},$
(30)

which shows that the interface temperature is a function of the grouping
$\mathrm{Ra}^{2/3}\mathrm{Da}$ alone in this limit. The height of the layers
$\hat{h}_{p}$ and $\hat{h}_{f}$ do not appear in 30 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 30 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
$\mathrm{Ra}$ and $\mathrm{Ra}_{p}$ in the simulations. The figure
demonstrates that $\theta_{i}$ decreases with decreasing
$\mathrm{Ra}^{2/3}\mathrm{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 $\mathrm{Ra}$ but very small $\mathrm{Da}$, such that
$\mathrm{Ra}^{2/3}\mathrm{Da}$ remains small, the interface temperature from
29 satisfies

$\theta_{i}\sim\sqrt{\mathrm{Ra}\mathrm{Da}^{3/2}},$ (31)

for $\mathrm{Ra}^{2/3}\mathrm{Da}\ll 1$. Conversely, $\theta_{i}$ increases
with increasing $\mathrm{Ra}^{2/3}\mathrm{Da}$ but evidently cannot increase
without bound, which reflects the fact that the grouping
$\mathrm{Ra}^{2/3}\mathrm{Da}$ cannot take arbitrarily large values. For
$\mathrm{Ra}^{2/3}\mathrm{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 17).

###  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 section 3, we found that while weak flows do
exist in the porous-stable case, they do not enhance heat transport compared
to 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 section 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 due 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{\mathrm{Ra}}\mathrm{Da}\theta_{i}$. The
heat transported on average by the penetrative flow scales with
$w\theta\sim\sqrt{\mathrm{Ra}}\mathrm{Da}\theta_{i}^{2}$, and so its
contribution to the Nusselt number $\mathrm{Nu}$ is
$\sim\mathrm{Ra}\mathrm{Da}\theta_{i}^{2}$. The contribution of the
penetrative flow thus becomes significant, relative to the $O(1)$ diffusive
flux through the layer, when $\mathrm{Ra}\mathrm{Da}\theta_{i}^{2}\sim O(1)$.
But according to the upper bound in 29,
$\mathrm{Ra}\mathrm{Da}\theta_{i}^{2}\lesssim 50\sqrt{\mathrm{Ra}}\mathrm{Da}$
for $\mathrm{Ra}\gg 1$. The assumption that penetrative flows do not transport
appreciable heat is therefore accurate for as long as
$50\sqrt{\mathrm{Ra}}\mathrm{Da}\lesssim 1$, or $\mathrm{Ra}_{p}\lesssim
4\times 10^{-4}\mathrm{Da}^{-1}$. Provided $\mathrm{Da}<O(10^{-5})$, this
value of the porous Rayleigh number is always larger than the critical value
$\mathrm{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 $\mathrm{Ra}_{p}$.

Given that in most geophysical systems we expect $\mathrm{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
$\mathrm{Ra}\gg 1$ penetrative convection can only occur if the Darcy number
is such that either the constraint $\sqrt{\mathrm{Ra}}\mathrm{Da}\ll 1$ 16 \-
which enforces that the flow in the porous medium remains confined with
negligible inertial effects - or the constraint
$\mathrm{Ra}\mathrm{Da}^{3/2}\lesssim 50$ 17 \- which enforces that the flow
lengthscales 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 timescales while carrying the same average heat flux. The
dynamics in the unconfined fluid layer must, therefore, exhibit variability on
both a slow timescale imposed by the porous layer below and on a rapid
timescale 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}-\mathrm{Ra}^{-1/2}\overline{\theta}^{\prime}$ 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 to its average value, with a duration that is
controlled by the porous turnover timescale
$\tau_{p}\sim(\sqrt{\mathrm{Ra}\mathrm{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 timescale. 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.

###  5.2 Spectral content of the heat flux at different heights

Figure 12: Time series of the $x$-averaged flux
$J(z,t)=\overline{w\theta}-\mathrm{Ra}^{-1/2}\overline{\theta}^{\prime}$ at
different heights in a simulation carried out at $\mathrm{Ra}=10^{8}$ and
$\mathrm{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 timescale. The present simulation spans over more than one diffusive
timescale since $\tau_{\mathrm{diff}}=10^{4}$.

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
$|\hat{J}(\omega,z)|^{2}$ of the signals displayed in figure 12. First, the
spectral content that is inherent to fluid convection is easily distinguished
from the one that is imposed by porous convection. The variability of the flux
in the fluid layer ($z=-0.5$ in figure 13(a)) 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.

###  5.3 Spatial decay rate of low-frequency variability

Figure 13:  a: Power spectrum $|\hat{J}(\omega,z)|^{2}$ of the heat flux
timeseries shown in figure figure 12 and normalised by $\left\langle
J\right\rangle^{2}$. b: Spatial variations of the power spectrum
$|\hat{J}(\omega,z)|^{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 $\mathrm{Ra}=10^{8}$
with $\mathrm{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.

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(|\hat{J}(\omega,z)|^{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
13c. 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, since 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 timescale $\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}$ are 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 $\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{\mathrm{Ra}}\mathrm{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
comprised 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,
$\mathrm{Ra}$, and small Darcy number, $\mathrm{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 is strongly dependent on
whether the flow in the porous layer is unstable to convection, which we
demonstrated should occur if $\mathrm{Ra}\gtrsim 27/(\hat{h}_{p}\mathrm{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 $\mathrm{Ra}$, we find that while flow in the
porous layer remains stable ($\mathrm{Ra}\mathrm{Da}\hat{h}_{p}\lesssim 27$)
the interface temperature $\theta_{i}$ satisfies
$\theta_{i}\sim\mathrm{Ra}^{-1/4}$, whereas if the porous layer becomes
unstable to convection, $\theta_{i}\sim\sqrt{\mathrm{Ra}\mathrm{Da}^{3/2}}$,
provided $\mathrm{Da}$ remains sufficiently small that
$\mathrm{Ra}\mathrm{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 marginally depend on the choice of
the formulation of the two-layer convection problem (see the discussion in
section 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{\mathrm{Ra}}\mathrm{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 parameterised
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. In other words, 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 be due to 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 et al., 2019).

[Acknowledgements]The authors are grateful to Eric Hester for his help in
implementing the code Dedalus.

[Funding]This work was performed using resources provided by the Cambridge
Service for Data Driven Discovery (CSD3) operated by the University of
Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell
EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences
Research Council (capital grant EP/P020259/1), and DiRAC funding from the
Science and Technology Facilities Council (www.dirac.ac.uk). TLR is supported
by the Royal Society through a Newton International Fellowship (grant nr.
NIF\R1\192181).

[Declaration of interests]The authors report no conflict of interest.

## References

  * Ahlers et al. (2009) Ahlers, Guenter, Grossmann, Siegfried & Lohse, Detlef 2009 Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection. Reviews of Modern Physics 81 (2), 503–537.
  * Amit et al. (2020) Amit, Hagay, Choblet, Gaël, Tobie, Gabriel, Terra-Nova, Filipe, Čadek, Ondřej & Bouffard, Mathieu 2020 Cooling patterns in rotating thin spherical shells — Application to Titan’s subsurface ocean. Icarus 338, 113509.
  * Bagchi & Kulacki (2014) Bagchi, Aniruddha & Kulacki, Francis A. 2014 Natural Convection in Superposed Fluid-Porous Layers. New York, NY: Springer New York.
  * Beavers & Joseph (1967) Beavers, Gordon S. & Joseph, Daniel D. 1967 Boundary conditions at a naturally permeable wall. Journal of Fluid Mechanics 30 (1), 197–207.
  * Burns et al. (2020) Burns, Keaton J., Vasil, Geoffrey M., Oishi, Jeffrey S., Lecoanet, Daniel & Brown, Benjamin P. 2020 Dedalus: A flexible framework for numerical simulations with spectral methods. Physical Review Research 2 (2), 023068.
  * Chandesris et al. (2013) Chandesris, Marion, D’Hueppe, A., Mathieu, Benoit, Jamet, Didier & Goyeau, Benoit 2013 Direct numerical simulation of turbulent heat transfer in a fluid-porous domain. Physics of Fluids 25, 125110.
  * Chen & Chen (1988) Chen, F. & Chen, C. F. 1988 Onset of Finger Convection in a Horizontal Porous Layer Underlying a Fluid Layer. Journal of Heat Transfer 110 (2), 403–409.
  * Chen & Chen (1992) Chen, Falin & Chen, C. F. 1992 Convection in superposed fluid and porous layers. Journal of Fluid Mechanics 234 (-1), 97.
  * Cheng et al. (2015) Cheng, J. S., Stellmach, S., Ribeiro, A., Grannan, A., King, E. M. & Aurnou, J. M. 2015 Laboratory-numerical models of rapidly rotating convection in planetary cores. Geophysical Journal International 201 (1), 1–17.
  * Choblet et al. (2017) Choblet, Gaël, Tobie, Gabriel, Sotin, Christophe, Běhounková, Marie, Čadek, Ondřej, Postberg, Frank & Souček, Ondřej 2017 Powering prolonged hydrothermal activity inside Enceladus. Nature Astronomy p. 1.
  * Coumou et al. (2008) Coumou, D., Driesner, T. & Heinrich, C. A. 2008 The Structure and Dynamics of Mid-Ocean Ridge Hydrothermal Systems. Science 321 (5897), 1825–1828.
  * d’Hueppe et al. (2011) d’Hueppe, A., Chandesris, M., Jamet, D. & Goyeau, B. 2011 Boundary conditions at a fluid–porous interface for a convective heat transfer problem: Analysis of the jump relations. International Journal of Heat and Mass Transfer 54 (15), 3683–3693.
  * d’Hueppe et al. (2012) d’Hueppe, A., Chandesris, M., Jamet, D. & Goyeau, B. 2012 Coupling a two-temperature model and a one-temperature model at a fluid-porous interface. International Journal of Heat and Mass Transfer 55 (9), 2510–2523.
  * Favier et al. (2019) Favier, B., Purseed, J. & Duchemin, L. 2019 Rayleigh–Bénard convection with a melting boundary. Journal of Fluid Mechanics 858, 437–473.
  * Fontaine & Wilcock (2007) Fontaine, Fabrice J. & Wilcock, William S. D. 2007 Two-dimensional numerical models of open-top hydrothermal convection at high Rayleigh and Nusselt numbers: Implications for mid-ocean ridge hydrothermal circulation. Geochemistry, Geophysics, Geosystems 8 (7).
  * Goodman (2004) Goodman, Jason C. 2004 Hydrothermal plume dynamics on Europa: Implications for chaos formation. Journal of Geophysical Research 109 (E3).
  * Goodman & Lenferink (2012) Goodman, Jason C. & Lenferink, Erik 2012 Numerical simulations of marine hydrothermal plumes for Europa and other icy worlds. Icarus 221 (2), 970–983.
  * Hester et al. (2021) Hester, Eric W., Vasil, Geoffrey M. & Burns, Keaton J. 2021 Improving accuracy of volume penalised fluid-solid interactions. Journal of Computational Physics 430, 110043\.
  * Hewitt (2020) Hewitt, D. R. 2020 Vigorous convection in porous media. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 476 (2239), 20200111.
  * Hewitt et al. (2012) Hewitt, Duncan R., Neufeld, Jerome A. & Lister, John R. 2012 Ultimate Regime of High Rayleigh Number Convection in a Porous Medium. Physical Review Letters 108 (22).
  * Hirata et al. (2007) Hirata, S.C., Goyeau, B., Gobin, D., Carr, M. & Cotta, R.M. 2007 Linear stability of natural convection in superposed fluid and porous layers: Influence of the interfacial modelling. International Journal of Heat and Mass Transfer 50 (7-8), 1356–1367.
  * Hirata et al. (2009a) Hirata, S.C., Goyeau, B., Gobin, D., Chandesris, M. & Jamet, D. 2009a Stability of natural convection in superposed fluid and porous layers: Equivalence of the one- and two-domain approaches. International Journal of Heat and Mass Transfer 52 (1-2), 533–536.
  * Hirata et al. (2009b) Hirata, S. C., Goyeau, B. & Gobin, D. 2009b Stability of Thermosolutal Natural Convection in Superposed Fluid and Porous Layers. Transport in Porous Media 78 (3), 525–536.
  * Howard (1966) Howard, L. N. 1966 Convection at high Rayleigh number. In Applied Mechanics (ed. Henry Görtler), pp. 1109–1115. Berlin, Heidelberg: Springer.
  * Hsu et al. (2015) Hsu, Hsiang-Wen, Postberg, Frank, Sekine, Yasuhito, Shibuya, Takazo, Kempf, Sascha, Horányi, Mihály, Juhász, Antal, Altobelli, Nicolas, Suzuki, Katsuhiko, Masaki, Yuka, Kuwatani, Tatsu, Tachibana, Shogo, Sirono, Sin-iti, Moragas-Klostermeyer, Georg & Srama, Ralf 2015 Ongoing hydrothermal activities within Enceladus. Nature 519 (7542), 207–210.
  * Huguet et al. (2016) Huguet, Ludovic, Alboussière, Thierry, Bergman, Michael I., Deguen, Renaud, Labrosse, Stéphane & Lesœur, Germain 2016 Structure of a mushy layer under hypergravity with implications for Earth’s inner core. Geophysical Journal International 204 (3), 1729–1755.
  * Le Bars & Worster (2006) Le Bars, M. & Worster, M. Grae 2006 Interfacial conditions between a pure fluid and a porous medium: Implications for binary alloy solidification. Journal of Fluid Mechanics 550 (-1), 149\.
  * Le Reun & Hewitt (2020) Le Reun, Thomas & Hewitt, Duncan R. 2020 Internally Heated Porous Convection: An Idealized Model for Enceladus’ Hydrothermal Activity. Journal of Geophysical Research: Planets 125 (7), e2020JE006451.
  * Malkus (1954) Malkus, Willem V. R. 1954 The heat transport and spectrum of thermal turbulence. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 225 (1161), 196–212.
  * Nield & Bejan (2013) Nield, Donald A. & Bejan, Adrian 2013 Heat Transfer Through a Porous Medium. In Convection in Porous Media (ed. Donald A. Nield & Adrian Bejan), pp. 31–46. New York, NY: Springer.
  * Nield & Bejan (2017) Nield, Donald A. & Bejan, Adrian 2017 Convection in Porous Media. Springer.
  * Nimmo et al. (2018) Nimmo, F., Barr, A. C., Běhounková, M. & McKinnon, W. B. 2018 Geophysics and Tidal-Thermal Evolution of Enceladus. In Enceladus and the Icy Moons of Saturn. The University of Arizona Press.
  * Plumley & Julien (2019) Plumley, Meredith & Julien, Keith 2019 Scaling Laws in Rayleigh-Bénard Convection. Earth and Space Science 6 (9), 1580–1592.
  * Poulikakos et al. (1986) Poulikakos, D., Bejan, A., Selimos, B. & Blake, K. R. 1986 High Rayleigh number convection in a fluid overlaying a porous bed. International Journal of Heat and Fluid Flow 7 (2), 109–116.
  * Priestley (1954) Priestley, C. H. B. 1954 Convection from a Large Horizontal Surface. Australian Journal of Physics 7 (1), 176–201.
  * Soderlund (2019) Soderlund, Krista M. 2019 Ocean Dynamics of Outer Solar System Satellites. Geophysical Research Letters 46 (15), 8700–8710.
  * Soderlund et al. (2014) Soderlund, K. M., Schmidt, B. E., Wicht, J. & Blankenship, D. D. 2014 Ocean-driven heating of Europa’s icy shell at low latitudes. Nature Geoscience 7 (1), 16–19.
  * Spencer et al. (2018) Spencer, J. R., Nimmo, F., Ingersoll, Andrew P., Hurford, T. A., Kite, E. S., Rhoden, A. R., Schmidt, J. & Howett, C. J. A. 2018 Plume Origins and Plumbing: From Ocean to Surface. In Enceladus and the Icy Moons of Saturn (ed. Paul M. Schenk, Roger N. Clark, Carly J. A. Howett, Anne J. Verbiscer & J. Hunter Waite), pp. 163–174. Tucson, AZ: University of Arizona Press.
  * Spencer et al. (2006) Spencer, J. R., Pearl, J. C., Segura, M., Flasar, F. M., Mamoutkine, A., Romani, P., Buratti, B. J., Hendrix, A. R., Spilker, L. J. & Lopes, R. M. C. 2006 Cassini Encounters Enceladus: Background and the Discovery of a South Polar Hot Spot. Science 311 (5766), 1401–1405.
  * Su et al. (2015) Su, Yan, Wade, Aaron & Davidson, Jane H. 2015 Macroscopic correlation for natural convection in water saturated metal foam relative to the placement within an enclosure heated from below. International Journal of Heat and Mass Transfer 85, 890–896.
  * Travis et al. (2012) Travis, B. J., Palguta, J. & Schubert, G. 2012 A whole-moon thermal history model of Europa: Impact of hydrothermal circulation and salt transport. Icarus 218 (2), 1006–1019.
  * Travis & Schubert (2015) Travis, B. J. & Schubert, G. 2015 Keeping Enceladus warm. Icarus 250, 32–42.
  * Urban et al. (2014) Urban, Pavel, Hanzelka, Pavel, Musilová, Věra, Králík, Tomáš, Mantia, Marco La, Srnka, Aleš & Skrbek, Ladislav 2014 Heat transfer in cryogenic helium gas by turbulent Rayleigh–Bénard convection in a cylindrical cell of aspect ratio 1. New Journal of Physics 16 (5), 053042.
  * Wang & Ruuth (2008) Wang, Dong & Ruuth, Steven J. 2008 Variable step-size implicit-explicit linear multisept methods for time-dependent partial differential equations. Journal of Computational Mathematics 26 (6), 838–855.
  * Wells et al. (2019) Wells, Andrew J., Hitchen, Joseph R. & Parkinson, James R. G. 2019 Mushy-layer growth and convection, with application to sea ice. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 377 (2146), 20180165\.
  * Woods (2010) Woods, Andrew W. 2010 Turbulent Plumes in Nature. Annual Review of Fluid Mechanics 42 (1), 391–412.
  * Worster (1997) Worster, M. Grae 1997 Convection in mushy layers. Annual Review of Fluid Mechanics 29 (1), 91–122.