High-Rayleigh-number convection in porous-fluid layers

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\'enard 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.


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) . 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, 1992Bagchi & Kulacki 2014;Su et al. 2015), or where restricted to the onset of convective instabilities (Chen & Chen 1988;Hirata et al. 2007Hirata et al. , 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 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 twolayer 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 singledomain 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 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 ℎ = ℎ = ℎ.
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.

Governing equations and numerical methods
2.1. The Darcy-Brinkman model Consider a two-dimensional system comprising a fluid-saturated porous medium of depth ℎ underlying an unconfined fluid region of depth ℎ . We locate the centre of the cell at height = 0, such that the whole system lies between −ℎ ℎ , as depicted in figure 1, and we introduce the length scale ℎ = (ℎ + ℎ )/2. For the majority of this paper we focus on the case of equal layer depths, where ℎ = ℎ = ℎ.
The fluid has a dynamic viscosity and density , and the porous medium is characterised by a uniform permeability 0 and porosity 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 Given local fluid velocity ℓ , the flux = ℓ 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 In the fluid later, the flow is governed by the Navier-Stokes equation, where is the pressure, while in the porous layer, the flux instead obeys Darcy's law, 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), where 1/ ( ) is a step function that goes from 1/ 0 for < 0 to zero for > 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 √ 0 . As a consequence, any spatial variation must have a typical length larger than √ 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, = 1 and 1/ = 0, which trivially gives the Navier-Stokes equation, whereas in the bulk of the porous medium, the damping term − / dominates over the inertial and viscous forces (provided 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 Balancing the viscous resistance and Darcy drag indicates that local viscous forces play a role over a length ℓ = √ 0 / 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 ℓ .
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 ∼ √ 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. (2007Hirata et al. ( , 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.

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 with 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 = 0 in the inertial terms. The temperature evolves according to an energy transport equation (Nield & Bejan 2013) where and are the heat capacity per unit of mass of the fluid and the porous matrix respectively, is the density of the porous matrix, and and 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.

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 where Δ > 0 is a constant. The upper and lower boundaries are rigid and impermeable, so Note that Darcy's law would only permit one velocity boundary condition on the boundary of the porous region at = −ℎ, 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 ∼ √ 0 / 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. 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 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 (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 = Δ ℎ, (2.12) and the associated free-fall timescale is * = ℎ/ * . Scaling lengths with ℎ, flux with * , time with * , temperature with Δ and pressure with * , we arrive at dimensionless equations where is the dimensionless flux, = ( − 0 )/Δ is the dimensionless temperature and ( ) is a step function that jumps from 1 in < 0 to 0 in > 0. In these equations, we have introduced three dimensionless numbers: (2.14) The Darcy number Da is a dimensionless measure of the pore scale √ 0 relative to the domain size ℎ, 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 ≫ 1. The Prandtl number compares momentum and heat diffusivities. The dimensionless layer depthsĥ andĥ are also, in general, variables; as noted above, in the majority of computations shown here we set these to be equal so thatĥ =ĥ = 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 (1), compared with ( √ RaDa/ √ Pr) in the porous layer. Given these scales, we can also introduce a porous Rayleigh number Ra 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 (2.13b), gives (2.15) 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, cancels out of the equations (see for instance (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 √ Da/ (as shown by (2.6)). In the following, we thus take = 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 section 2.7.
In the heat transport equation (2.13b), we reduce the number of control parameters by taking = Λ = 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 ∼ . The parameters and Λ 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 Λ 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 Λ.
Finally, to reduce the complexity of this study and maintain a focus on the key features of varying the driving buoyancy forces (i.e. Ra) and the properties of the porous matrix (i.e. 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 Ra and 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 = Λ = 1. First, the inertial terms must vanish in the porous medium, which demands that that is, the velocity scale in the porous medium must be much less than the (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 √ 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 ∼ √ RaDa/ √ Pr, yields a lengthscale ∼ Ra −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 50Ra −1 . For these structures to remain larger than the pore scale √ Da, we thus require RaDa 3/2 50Λ. (2.17) Note that the effect of violating this constraint is to amplify the importance of viscous resistance ∇ 2 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 Ra and 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.

Numerical method
The one-domain Darcy-Brinkman equations (2.13) are solved using the pseudo-spectral code D (Burns et al. 2020;Hester et al. 2021). The flow is decomposed into 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 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 √ RaDa, 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 = 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 Ra and Da for various fixed values of the porous velocity scale √ RaDa. The majority of simulations were carried out in a set-up where the heights of both the porous and the fluid layer were equalĥ =ĥ = 1, with resolutions × [ , ] = 1024 × [128, 256] below Ra = 10 8 and 1024 × [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 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 √ Ra, to ensure that the flow had reached a statistically steady state.

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, Da = 10 −5.5 , and equal layer depthsĥ =ĥ = 1, but with varying Ra in the range 10 4 Ra 10 9 . We use these to illustrate the basic features of high-Ra convective flow in the two-layer system.  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, ( ), are shown in figure 4(a), while the mean interface temperature = (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 in the bulk of the fluid.

Two different regimes depending on the stability of the porous medium
The porous layer, however, exhibits two different behaviours depending on the size of Ra. At low Rayleigh numbers (Ra 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 < 0 (see figure 3a,b), while the horizontally averaged temperature profiles ( ) appear to be linear (figure 4a). The corresponding interface temperature monotonically decreases with Ra (figure 4b) As the Rayleigh number is increased beyond Ra ∼ 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 , which reverses from decreasing to increasing with Ra around Ra ∼ 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 Ra = RaDa exceeds a critical value Ra ≃ 27.1 (Nield & Bejan 2017). At Da = 10 −5.5 , the critical Rayleigh number Ra such that RaDa = Ra is Ra ≃ 8.6 × 10 6 , which is reported in figure 4(b) and agrees well with the inversion of trend in . 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  (4.7) using C = 0.85 and the marked values C , as detailed in section section 4.2.
Figure 5: a: Nusselt number Nu = 2 √ Ra characterising heat transport across the two-layer system, together with predictions from section 4 for the porous-stable case (4.8) (red dashed) and porous-unstable (4.7) (black dotted and dot-dashed, with C = 0.85 and C as marked). The horizontal line marks Nu = 2. b: Root-mean-squared vertical velocity amplitude rms in the porous and fluid layers. The two grey lines indicate a scaling of √ RaDa (the characteristic speed in the porous layer). c: Ratio of the diffusive ( d,p ) and advective ( a,p ) fluxes to the total depth-averaged flux in the porous layer. In all three panels, Da = 10 −5.5 and the vertical line marks the threshold of porous convection estimated using (4.9).
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.

Characteristics of heat transport
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 ( ) ≡ − Ra −1/2 ′ , 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 The insets show zooms to make clearer the different levels of temporal variations, with horizontal segments indicating the porous turnover timescale and the free-fall timescale. Note that this simulation has been started from a diffusive state with a small noise.
purely diffusive system, Ra −1/2 /2, via the Nusselt number Nu ≡ 2 √ Ra , 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 (Ra 10 7 ), Nu appears to approach a horizontal asymptote Nu = 2, but once the porous layer is unstable 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 tends to zero. In this limit, → 1/ √ Ra and Nu → 2. The decreasing in figure 4(b) as 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 , 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 Nu(Ra, Da) until the following section.

Timescales, variability and statistically steady state
The governing equations (2.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 (1) in our free-fall normalisation. The second is given by diffusion, diff = √ Ra, and the third is the turnover timescale in the porous layer ∼ ( √ RaDa) −1 , which scales with the inverse of the porous speed scale. Because and measure advection and diffusion in the porous medium, these timescales are in a ratio = diff /Ra . In addition, we recall that ≫ 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 timescales should scale with the inverse of rms in each layer, as can be observed in figure 5(b): in the fluid layer, rms ∼ (1) with no systematic variation with , while in the porous layer, rms ∝ √ Ra at constant 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 ∼ 0.2 diff . Although the equilibration time is largely reduced by the use of continuation, we run all simulations over times that are similar to diff to ensure they are converged.

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.

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 √ RaDa ≪ 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 (i.e. R = Ra in the fluid layer and R = RaDa in the porous layer). If the boundary layer has mean thickness , then we can also introduce a local boundary-layer Rayleigh number , where = 1 in the porous layer or = 3 in the fluid layer to account for the different dependence on the height scale in the corresponding Rayleigh number (see (2.14) and (2.15)). Let the temperature difference across the boundary layer be Θ ( figure 7), and, for completeness, suppose that the region has an arbitrary dimensionless heightĥ.
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Θ and lengths byĥ, to give a new Rayleigh numberR = 2ΘĥR and boundary-layer depthˆ = /ĥ. Having rescaled thus, the Malkus-Howard approach would suggest that Rˆ = R , (4.1) for some critical value R , orˆ = (R /R) 1/ . The corresponding Nusselt number N for this rescaled system, given that the scaled temperature drop across the layer is 1/2, is Note that the actual, unscaled flux across the boundary layer is = Ra −1/2 Θ/ , which can thus be related to N from (4.1) and (4.2) by Thus specification of the critical value R yields a prediction of the flux through the system in terms of Θ andĥ. We can extract values of R from previous works that have determined experimentally or numerically the relation between N andR in either a fluid or a porous Rayleigh-Bénard cell. , suggest that (2R 1/ ) −1 ∼ 6−8×10 −2 , although no definitive observation of a well-developed Ra 1/3 law has been made so far and the 'true' asymptotic form of N(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 Ra ≫ 1 and √ RaDa ≪ 1 that can be compared to our simulations, as discussed in the next section.

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 N (R) of the rescaled Rayleigh number, rather than the asymptotic scaling. To achieve this, we simply replace (4.2) by the relationship N = N (ˆ ). (4.5) 10 5 10 6 10 7 10 8 10 9 R f = Ra × 2Θ fĥ 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 Θ ∼ 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 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 numberR in both the fluid layer (index ) and in the porous layer (index ), with the temperature drop across boundary layers Θ 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ĥ andĥ 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 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.  porous Rayleigh number is large so that C = 0.85 whereas C varies from 1 in the weakly confined limit ( √ RaDa 1) to 0.5 in the strongly confined limit ( √ RaDa ≪ 1).

The interface temperature
We can use the laws governing the heat flux to determine the interface temperature , 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 Θ across the boundary layers to the interface temperature . 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 (1) coefficients C , such that Θ = C and Θ = C (1 − ) (see figure 7). In single classical Rayleigh-Bénard cells, C , = 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 C , may take any value ranging from 0.5 to 1.
Given these coefficients, flux conservation between the layers yields which determines the interface temperature . In general, the fractions C , depend on the control parameters Ra and Da, as shown in figure 9.

The porous-convective regime
We first address the case where both layers are unstable to convection. Although C and C 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 = √ RaDa 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 C = C = 0.5. In the opposite, weakly confined limit where √ RaDa 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 Ra = 10 9 temperature profile in figure 4a for which √ RaDa = 10 −1 ). In such a limit, C = C = 1. The values of C and C 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 C varies with Ra it is also larger for larger values of √ RaDa. 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 Ra → ∞, as sketched in figure 9(c), C = 0.5 in the strongly confined limit ( √ RaDa ≪ 1 ) and C = 1 in the weakly confined limit ( √ RaDa 1). In the porous layer, the numerical data shows difference with the qualitative expactations for the values of C (see figure 9a). C is mainly a function of the porous Rayleigh number Ra = RaDa (figure 9a), and appears to approach an asymptote C ≃ 0.85. Interestingly, the trend of C with Ra 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 (≃ 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 C ≃ 0.85 rather reflects the temperature decrease in the bulk of the porous layer (see the Ra = 10 9 temperature profile in figure 4a).
Before showing asymptotic predictions based on these limiting values of C , , 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 Ra ≫ 1 for the two limiting cases C = 0.5 and C = 1.0 with C = 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 √ RaDa show interface temperatures approaching the limit C = 1 (figure 10a), while the more strongly confined simulations with √ RaDa ≪ 1 draw closer to C = 0.5 as Ra is increased (figure 10d), which suggests that these limits will become increasingly accurate for increasingly large .

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 = Ra −1/2 (1 − )/ĥ . A balance of these expressions should then replace (4.7) to predicted with the asymptotic model (4.13) as a function of Ra 2/3 Da. The two lines materialise the two limiting cases of strongly confined porous medium ( √ RaDa → 0, C = 0.85 and C = 0.5) and the weakly confined case ( √ RaDa 1, C = 0.85 and C = 1.0). Note that Ra 2/3 Da must remain smaller thant ∼ 10 to ensure that the constraint (2.17) is satisfied.
determine . Moreover, we know a priori that C = 0.5 since the flux is entirely diffusive at the interface = 0. Flux conservation thus reduces to 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 , we can also extract the Nusselt number for the two-layer system, Nu = 2 √ Ra = 2(1 − )/ĥ (see section 3), which was also shown giving good agreement with numerical data in figure 5(a) for the particular caseĥ = 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/ĥ when the porous layer remains stable.
In fact, prediction of 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 RaDaĥ (1 − ) reaches a critical value Ra ≃ 27.1 (Nield & Bejan 2017). Because → 0 at large Rayleigh numbers while the porous layer remains stable, this condition may be recast as at the threshold of porous convection. To leading order, the threshold is given by Ra = Ra /(Daĥ ), as already noted in section 3.1. For the value of Da used in section 3 wherê ℎ = 1, including the first-order correction increases the predicted threshold value of Ra by about 10%.

Asymptotic predictions
These relationships for , 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 (RaDa ≫ 1 while √ RaDa ≪ 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: N = 6.9 × 10 −3R ≡ R and N ≃ 7 × 10 −2R1/3 ≡ R 1/3 . (4.10) In the case of a stable porous medium, flux balance (4.8) reduces to with an asymptotic upper bound ∼ 7.3 ℎ Ra 1/4 (4.12) for Ra ≫ 1. In the case where both layers are convecting, (4.7) instead reduces to (4.13) which shows that the interface temperature is a function of the grouping Ra 2/3 Da alone in this limit. The height of the layersĥ andĥ 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 C , 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 Ra and Ra in the simulations. The figure demonstrates that decreases with decreasing Ra 2/3 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 . In the asymptotic limit of large Ra but very small Da, such that Ra 2/3 Da remains small, the interface temperature from (4.12) satisfies ∼ RaDa 3/2 , (4.14) for Ra 2/3 Da ≪ 1. Conversely, increases with increasing Ra 2/3 Da but evidently cannot increase without bound, which reflects the fact that the grouping Ra 2/3 Da cannot take arbitrarily large values. For Ra 2/3 Da (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)).

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, ( = −ĥ ) = 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 . Any penetrative flow must be driven by these horizontal variations, and so will have an amplitude ∼ √ RaDa . The heat transported on average by the penetrative flow scales with ∼ √ RaDa 2 , and so its contribution to the Nusselt number Nu is ∼ RaDa 2 . The contribution of the penetrative flow thus becomes significant, relative to the (1) diffusive flux through the layer, when RaDa 2 ∼ (1). But according to the upper bound in (4.12), RaDa 2 50 √ RaDa for Ra ≫ 1. The assumption that penetrative flows do not transport appreciable heat is therefore accurate for as long as 50 √ RaDa 1, or Ra 4 × 10 −4 Da −1 . Provided Da < (10 −5 ), this value of the porous Rayleigh number is always larger than the critical value Ra 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 Ra .
Given that in most geophysical systems we expect Da ≪ 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 Ra ≫ 1 penetrative convection can only occur if the Darcy number is such that either the constraint √ RaDa ≪ 1 (2.16) -which enforces that the flow in the porous medium remains confined with negligible inertial effects -or the constraint RaDa 3/2 50 (2.17) -which enforces that the flow lengthscales are larger than the pore scale -are violated.

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.

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 ( , ) = − Ra −1/2 ′ 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 ∼ ( √ RaDa) −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 (1) free-fall timescale. However, comparison between the time series at = 0 and = 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.

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 |ˆ ( , )| 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 ( = −0.5 in figure 13(a)) is almost entirely contained in harmonics smaller than ∼ 10 −2 . The energy of lower harmonics in the fluid layer ( 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 ( = 0 & = 1) being effectively identical above = 5 × 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 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 and with a rate that seems to increase with . 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.

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(|ˆ ( , )| 2 ) with ∈ [0, 0.6], we extract the decay rate for all frequencies below = 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 ≃ 5 × 10 −3 , the energy decay rate roughly increases like 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 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, ∼ 1/2 for oscillation frequency . 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 . As → 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 needed to ensure that the decay of the variability in the fluid layer is negligible. The power spectrum of the flux at = 0 and =ĥ are in a ratio ∼ exp(− ( )ĥ ). The cut-off frequency of such a low-pass filter, reached when exp(− ĥ ) ∼ 1/2 is roughly located at = = 10 −3 in the particular caseĥ = 1 according to figure 13(c). Therefore, using ∼ 3/4 , we predict in general that when 1/ = √ RaDa is smaller than ∼ĥ −4/3 , the temporal variability imposed by porous convection is sufficiently slow that it will be entirely transmitted across the fluid layer.

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, Ra, and small Darcy number, 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 Ra 27/(ĥ 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 Ra, we find that while flow in the porous layer remains stable (RaDaĥ 27) the interface temperature satisfies ∼ Ra −1/4 , whereas if the porous layer becomes unstable to convection, ∼ RaDa 3/2 , provided Da remains sufficiently small that RaDa 3/2 ≪ 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 √ RaDa < 10 −3ĥ −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).