Internally heated porous convection: an idealised model for Enceladus' hydrothermal activity

Recent astronomical data and geophysical modelling suggest that hydrothermal activity is ongoing under the ice crust of Enceladus, one of Saturn's moons. According to these models, hydrothermal flow in the porous, rocky core of the satellite is driven by tidal deformation that induces dissipation and volumetric internal heating. Despite the effort in the modelling of Enceladus' interior, systematic understanding---and even basic scaling laws---of internally-heated porous convection and hydrothermal activity are still lacking. In this article, using an idealised model of an internally-heated porous medium, we explore numerically and theoretically the flows that develop close and far from the onset of convection. In particular, we quantify heat-transport efficiency by convective flows as well as the typical extent and intensity of heat-flux anomalies created at the top of the porous layer. With our idealised model, we derive simple and general laws governing the temperature and hydrothermal velocity that can be driven in the oceans of icy moons. In the future, these laws could help better constraining models of the interior of Enceladus and other icy satellites.


Introduction
Enceladus, a 500 km-diameter icy satellite orbiting Saturn, has drawn a lot of attention since the first flybys operated by the Cassini probe in 2005. Pictures and in situ astrochemical measurement have revealed the presence of a water-vapour and ice plume ejected into outer space. It emerges along fractures in the ice crust at the south pole of Enceladus and is associated with a large heat-flux anomaly of 12.5 GW (Porco et al., 2006). Subsequent analyses have revealed that the ejected material contains silicate particles of nanometric size whose chemistry indicates that the water contained in the plume has been previously hot, liquid, and in contact with silicate rocks (Hsu et al., 2015;Sekine et al., 2015). Enceladus' plumes have since then been interpreted as evidence for hydrothermal activity occurring below the ice crust of Enceladus. This is a surprising implication because, unlike the Earth, Enceladus has radiated away all its initial heat, and its small size makes internal heating by radiogenic elements insufficient to explain the abnormal heat flux (Nimmo & Pappalardo, 2016;Choblet et al., 2017).
Building on the recent study of Lainey et al. (2017), Choblet et al. (2017) have recently proposed a self-consistent model to explain the hydrothermal activity based on internal heating by tides in Enceladus' water-saturated porous core. This model relies on recent findings regarding the interior of Enceladus. Underneath its ice crust, this satellite comprise a global subsurface ocean, with thickness varying from 30 to 50 km (Thomas et al., 2016). Below lies a core made of rocky material that remains undifferentiated and uncompacted owing to to the weakness of Enceladus' gravity field (Roberts, 2015;Choblet et al., 2017). The core is thus permeated with the water of the ocean; Choblet et al. (2017) estimate that the porosity ranges from 20 to 30% for a water-filled rocky core. Lastly, intense internal heating (Lainey et al., 2017) due to tidally-induced deformation and fric-tion heats the water and creates a porous flow with hot and narrow upwelling zones, possibly leading to hot spots of water flowing into the ocean (Choblet et al., 2017).
Hydrothermal porous convection driven by volumetric heating is a standard feature of the models of the interiors of icy moons (Travis et al., 2012;Travis & Schubert, 2015), even though it is in general thought to be driven by radiogenic decay or rock serpentinisation rather than tidal deformation (Nimmo & Pappalardo, 2016). However, general scaling laws governing heat transport by tidally-driven hydrothermal systems are still lacking. As well as better constraining the sensitivity of the model of Choblet et al. (2017) with respect to the values of control parameters, deriving such laws could prove useful to constrain heat transport in other icy worlds and possibly improve its parametrisation in thermal evolution models of icy worlds (Travis et al., 2012;Travis & Schubert, 2015).
In the present article, we explore systematically internally-heated porous convection close to and far from the onset of motion with numerical simulations and mathematical analysis. We use an idealised two-dimensional Cartesian model of a water-saturated porous layer with internal heating in order to reduce the complexity of the system as much as possible while retaining the key physical ingredients. With this model, we derive scaling laws governing the width of upwelling zones and the associated buoyancy flux at the bottom of the ocean. Despite their simplicity, our laws are in agreement with the model of (Choblet et al., 2017) when applied to the case of Enceladus. This kind of approach has a long history of use in convection studies. The most canonical model is the Rayleigh-Bénard set-up (a confined porous layer heated from below and cooled from the top) which has received a significant amount of study (Otero et al., 2004;Hewitt et al., 2012Hewitt et al., , 2014Hewitt & Lister, 2017). The more closely related case of Earth-like hydrothermal systems with a bottom heat flux and open top boundary has also been widely studied (see for instance Fontaine and Wilcock (2007); Coumou et al. (2008Coumou et al. ( , 2009). However, the results of these investigations are unlikely to apply to tidally-driven hydrothermal circulation because of inadequate boundary condition and the nature of the heat source. Very few systematic experimental and numerical studies have been devoted to internallyheated porous convection. Those that have are focused mostly on the onset of motion and average heat transport (Buretta & Berman, 1976;Nield & Kuznetsov, 2013;Hardee & Nilson, 1977;Kulacki & Ramchandani, 1975). Hence, these studies do not allow the derivation of scalings governing, for instance, the typical extent of upwelling zones or the associated thermal anomalies and fluid velocities, in the case of tidally-driven hydrothermal activity. That is our aim here. This paper is organised as follows. A first part is devoted to introducing our idealised model for an internally-heated saturated porous layer and identifying the relevant dimensionless parameters. In common with numerous convection set-ups, the intensity of heat-transporting motion is characterised by a unique dimensionless number, the Rayleigh number. We then carry out a stability analysis to determine the conditions under which convection happens. Afterwards, we describe and analyse numerical simulations of internallyheated porous convection, focusing in particular on the structure of the flow and the associated thermal anomalies. Lastly, we apply the laws derived from our idealised model to Enceladus to quantify the temperature anomalies and the typical hydrothermal velocities that can be induced in its ocean.

The model and its governing equations
Rather than modelling Enceladus' full fluid system comprising the water-saturated core and the ocean, we choose to dwell on the core only and to parametrise the core-ocean interaction via a boundary condition that will be specified hereafter, as in (Choblet et  al., 2017). The core of Enceladus is modelled by a two-dimensional porous medium of permeability k, which is saturated with water of viscosity µ. It lies beneath an ocean that we assume to be well mixed with a global temperature T 0 (see figure 1). The flow U = (U, W ) inside the porous core is modelled by Darcy's law, where P is the pressure, ρ is the density of water and g is the gravity field, pointing in the z direction. For the sake of simplicity, gravity is assumed to be constant over the layer, even though in increases with height in the cores of moons. In addition to Darcy's law, the flow is assumed to be incompressible, so that it must also satisfy a continuity equation, Water motion inside the core is driven by buoyancy and temperature differences. We model the effects of temperature on density assuming linear expansion of the fluid with temperature under the Boussinesq approximation, such that ρ = ρ 0 (1 − α(T − T 0 )) where ρ 0 is a reference density and α the thermal expansion coefficient. Darcy's law may thus be written as, where P = P + ρ 0 gz and Θ ≡ T − T 0 . Since the flow is driven by thermal anomalies Θ, we must introduce an equation modelling the transport of heat inside the porous medium. This is achieved using thermal energy conservation, in which a source term accounting for volume heat production is included: with κ the volume-averaged heat diffusivity inside the porous medium (i.e. of both water and the porous matrix together) and q is the internal heat source term. The latter is related to the volume heat production by tidal heating Q V via q = Q V /(ρc p ) where c p is the heat capacity per unit of mass of the porous medium and ρ its density. In this paper, we consider two idealised limits: either Q v is constant or it is assumed to vary laterally (i.e. in x) to model tidal heating inhomogeneities.

Boundary conditions
The bottom of the porous layer roughly corresponds to the core centre. We thus assume that there is no heat and mass flux crossing the bottom boundary, that is: The top of the layer at z = h is in contact with the ocean and must allow mass exchange between the core and the ocean. This is achieved by imposing a purely vertical velocity at the top, i.e. : The two layers are also thermally coupled, and we consider two possible boundary conditions for θ on the upper boundary. One first natural choice is to impose the temperature (on the upper boundary) to be the temperature of the ocean, i.e. , However, in this case, the advective heat flux driving hydrothermal activity W Θ(z = h) across the interface would vanish, which seems at odds with the idea that the water coming out the porous layer may drive a buoyant plume rising in the ocean. We could alternatively use another boundary condition where the temperature of water is left unchanged as it leaves the porous layer, while water enters with the imposed temperature of the ocean, that is, Such a boundary condition is a standard parametrisation of core-ocean interactions (Rabinowicz et al., 1998;Monnereau & Dubuffet, 2002;Cserepes & Lenkey, 2004;Choblet et al., 2017).
The thermal boundary conditions (7) and (8) may be regarded as two end-members of the fully coupled problem of the core-ocean interaction. In the case of slow ascent in the porous medium, diffusion from the ocean inside the core causes the temperature inside the porous medium to drop in the top boundary vicinity. Conversely, if the upwelling is fast, diffusion is not able to affect the temperature inside the ascending plume. As a side note, intermediary situations where ∂ z Θ(z = h) = −β with β > 0 could also be considered. Nevertheless, choosing between the two boundary conditions or parametrisation of β would require a demanding study of the fully coupled system involving both the ocean and the porous core. We instead carry out simulations using either boundary conditions (7) or (8), and we will find that the choice of boundary condition does not significantly affect the flow in the interior of the core.

Scaling the problem: dimensionless equations
First, all considered lengths are normalised by the height of the porous layer h. We must also define velocity and temperature scales, respectively denoted as U * and Θ * . Darcy's law (3) gives a simple relation between these two scales, Unlike in, say, Rayleigh-Bénard set-up, the temperature scale Θ * is not naturally imposed in the internally heated problem. We predict that in the non-linear regime, heat production and advection will be the dominant balance, leading to the following relation between velocity and temperature in (4) Both scales then can be written as a function of physical parameters as follows: Given these scales, we find that the system is governed by only one dimensionless parameter, a Rayleigh number comparing the relative importance of advection and diffusion, Note that other definitions have been considered for the Rayleigh number, depending in particular on the expected balance at play. For instance, Buretta and Berman (1976) choose velocity and temperature scales based on an advection and diffusion balance, rather than a balance between advection and heat production consider in (10), leading to a Rayleigh number Ra bb = Ra 2 .
Introducing the dimensionless temperature θ = Θ/Θ * and velocity u = U /U * , the dimensionless governing equations for a porous layer with internal heating are: where the pressure is rescaled by hU * and time by h/U * . The flow being incompressible and two-dimensional, it is particularly convenient to introduce a stream function ψ such that u = ∇ × (−ψe y ). The governing equations (13) become Lastly the boundary conditions are and either BC 1: θ(z = 1) = 0 , or BC 2: for the temperature. In addition, we assume horizontal periodic boundary conditions in the x direction.

Numerical modelling
We study this problem numerically with the code developed by (Hewitt et al., 2012). It solves Darcy's law and the heat transport equation (14) with finite differences in time and the vertical direction, and a Fourier transform in the horizontal direction (where periodic boundary conditions are used). The use of two staggered grids for the stream function ψ and the temperature field θ allows the numerical scheme to be second order accurate in space. The finite differences are second-order accurate in time as well. Anticipating strong gradients, a vertical stretched grid is implemented to ensure the boundary layers are well resolved. The numerical implementation of equations (14) will be tested in section 3.

The onset of convection
In this section, we investigate both theoretically and numerically the critical value of the Rayleigh number Ra above which a convective instability develops. The steady, purely diffusive base (u = 0, θ b ) state on which the instability develops is (c) Temperature field (θ1) and streamlines (iso-contours of ψ1) of the unstable mode at the onset of convection.
regardless of the upper thermal boundary condition. We look for perturbations of the base state of the form: such that |ψ 1 |, |θ 1 | θ b . The exponential terms allow to account for the existence of convective instability characterised by Re(σ) > 0. Using the ansatz (19), equations (14) to leading order in ψ 1 , θ 1 yield the following single fourth-order differential equation on the stream function: The invariance under translation along the x-axis allows further simplification by assuming that ψ 1 is a plane wave in x, that is ψ 1 =ψ 1 (z) exp(ikx). Equation (20) with the plane wave assumption yields the following ordinary differential equation for the functionψ 1 :ψ where σ is an unknown eigenvalue. We solve (21) using BC 1 in (16). (In fact, for this onset problem, BC 2 (17) gives an ill-posed system.) From this numerical solving, we find the lowest value of the Rayleigh number for which σ = 0 to be Ra = Ra c 5.894 at k = k c 1.751. Such a value for the critical Rayleigh number is close to the value 5.72 found experimentally and theoretical by Buretta and Berman (1976) in a system with closed boundary conditions. The marginal mode and its vertical structure functions (ψ 1 andθ 1 ) are shown in figure 2. The mode comprises a half-roll structure, with strong horizontal flow at the lower boundary and strong vertical flow at the upper boundary. The temperature deviation is maximised roughly half-way up the roll.
We have also use this theoretical investigation of the onset of convection to benchmark the numerical code. Simulations were carried out at values of the Rayleigh number Ra very close to the onset (|Ra − Ra c | ≤ 10 −1 typically). The horizontal extent of the domain is chosen to match approximately twice the wave length of the marginal mode. Computations were initiated with a small perturbation to the diffusive temperature profile (18). We observed an exponential growth or decay of the amplitude of the perturbation to the diffusive base state and found accurate reproduction of both the critical Rayleigh number and the growth or decay rate of the most unstable mode for nearby values of Ra. Figure 2 shows the excellent agreement between the theoretical and the computed vertical structure functionsθ 1 andψ 1 .

Non-linear heat transport by convection
In the following section, we investigate heat transport by convection for larger values of Ra. We first describe qualitatively the organisation of the flow as Ra is increased. We then show quantitatively that non-linear heat transport is dominated by advection, which constrains the typical size of hot plumes and thermal anomalies. We use both thermal boundary conditions (16) and (17) to find that the difference between them is negligible for large enough values of Ra.

Numerical process
Prior to delving into the results of the simulations, we explain how a typical numerical is carried out. The simulations are initialised with a random noise at a certain Rayleigh number Ra. After the initial growth of the instability, the flow reaches a statistically steady state, which is assessed by the convergence of the cumulative mean of the maximum temperature. The simulation is terminated once the steady state has lasted for 300 time units. The Rayleigh number is then switched to a new value, and the simulations is initiated with the last state of the previous one plus a small noise disturbance. A summary of all the numerical simulations that have been carried out is given in table 1.

Flow structures and organisation
To introduce the non-linear behaviour of the instability driven by internal heating, we first to illustrate typical flow patterns observed at different Rayleigh numbers. Figures 3 displays typical snapshots of the temperature field. At low Rayleigh number, i.e. for Ra c ≤ Ra < 20, the convection reaches a steady state with few plumes, be it for boundary condition BC 1 or BC 2 (see figure 3a). Similarly to the unstable mode at threshold, these plumes consist of half-rolls, although with steeper vertical gradients at the top boundary in the case of BC 1. For larger Rayleigh numbers (see 3b), the flow  (16) and (17)  exhibits an unsteady chaotic behaviour where usually two modes with different number of plumes alternate, thus inducing chaotic merging and growth of plumes. This situation ceases for Ra 600, at least for an aspect ratio L = 4: higher values of the Rayleigh number give rise to steady solutions with a large number of narrow plumes (see figure  3c) The only noticeable difference between the two boundary conditions is the existence of a thin thermal boundary layer when the top temperature is imposed (BC 1). Its thickness, of order Ra −1 , is set by a balance between vertical advection and diffusion. In addition, the high degree of similarity between the simulations carried out with different boundary condition suggests that the mixed boundary condition (BC 2) is reliable. Note that this is not the case below the threshold of the instability where flows that are highly sensitive to initial condition are observed. We therefore choose to use both boundary conditions in the study detailed hereafter, as long as Ra > Ra c . Lastly, we notice on these snapshot that θ = O(1) which confirms that the balance between advection and heat production drives the dynamics, a balance that was foreseen in section 2.3.

Advective heat transport
The qualitative analysis of snapshots carried out in the preceding section indicates that advection dominates heat transport. We propose in the following a quantitative analysis of the flow to support this assertion, in particular of the vertical temperature and heat flux profiles. This analysis will allow us to compare internally heated porous convection with the more classical Rayleigh-Bénard problem via the introduction of a generalised Nusselt number.

The mean temperature scale
It has been noted in the preceding section that the typical values of the temperature field remain of O(1). To better quantify this observation, we introduce a dimensionless temperature scale ∆θ = θ(z = 1) − θ(z = 0), where the operation · denotes horizontal and temporal average in the statistically steady state. Typical profiles of the horizontally averaged temperature are shown in figure 4. We note again the strong similarity between the two boundary conditions, especially at large Rayleigh numbers where they only differ by the presence of the top thermal boundary layer. The temperature scale ∆θ is also plotted in figure 4 as a function of the Rayleigh number for all simulations. We note it is well below the diffusive scaling ∆θ ∝ Ra even very close to the thresh- old of the instability. ∆θ = O(1) is a signature of efficient transport and vertically mixing of the thermal energy by the convective flows. In addition, we note a marked decrease of ∆θ at Ra 600, which corresponds to the transition from the chaotic to the steady regime. It indicates that the steady flow is even more efficient at transporting heat out of the system. Anomalous points may however be noticed; they are due to the locking of the simulation on a particular mode (i.e. a flow with a certain number of plumes) that remains stable as the Rayleigh number is slightly increased.

The advective flux
To further quantify heat transport in the strongly non-linear regime, we consider here the vertical heat flux, defined as which comprises an advective and a diffusive contribution. Time-averaged thermal energy conservation (13) prescribes a balance between vertical heat transport and volumetric heat production such that: In the asymptotic regime of high Rayleigh number, we expect that the heat produced is carried away by advection only, apart from the thermal boundary layer when it exists. In the bulk of the porous medium, we thus expect wθ(z) = z .
As it can be noticed in figure 5, the advective heat flux is well described by the asymptotic law (24) even at Rayleigh numbers as low as Ra = 29. For BC 1, this agreement

Nusselt number
It is interesting to assess how efficient the convecting system is at transporting heat relative to purely diffusive transport. It is quantified by a Nusselt number N that provides a comparison between the total heat flux (including advective and diffusive contributions) and the diffusive heat flux (see e.g. Goluskin (2016)), · denoting volume average, and where we have used that wθ − Ra −1 ∂ z θ = z = 1/2 and ∂ z θ = ∆θ from (22) and (23). Then, we define the mean Nusselt number N u to be the long-time average of N (t). Note that we retrieve that the transport is purely diffusive at threshold, since, at Ra = Ra c , ∆θ = Ra/2 so that N u = 1. Because ∆θ = O(1), we predict that, in the high Rayleigh number regime, N u ∝ Ra. For both boundary conditions BC 1 & 2, our simulations confirm this scaling down to Ra ∼ 20 (see figure 5b). There is a slight enhancement of the efficiency of heat transport as steady states emerge in the non-linear regime of the instability around Ra ∼ 500. The same scaling between N u and Ra is also found in the classical Rayleigh-Bénard set-up in a porous medium (Otero et al., 2004;Hewitt et al., 2012Hewitt et al., , 2014.

Plume scales
In section 4.2, we observed that as the Rayleigh number Ra is increased, the typical width of the plumes and their typical spacing decreases. We obtain a quantitative measure of the mean plume size p and separation ∆x p as a function of the Rayleigh number from the heat flux at the upper boundary. As shown in figure 6a, plumes produce a series of heat-flux peaks. At each time step, we record the mean plume widthˆ p and plume separation distance∆x p (t) over all plumes, and we define p and ∆x to be their long-time averages. Typical variability is given by the standard deviation ofˆ p and∆x p over time. The result of this process is shown in figure 6(b,c): both the plume width and separation exhibit the same scaling with the Rayleigh number, that is p , ∆x p ∝ Ra −1/2 , even close to the threshold. This power law can't be explained by linear theory, even at low Ra, as the mean separation between plumes does not coincide with the most unstable mode predicted by the linear stability analysis (see figure 6c). Instead, the typical scale of the plume is controlled by a balance between vertical advection, horizontal diffusion and heat production in (13), that is, Given that temperature contrast remains O(1), this balance demands both that the vertical velocity of the plume is O(1) and that the typical lateral scale of the plume must be proportional to Ra −1/2 .

Asymptotic plume solution
Building on the scalings governing the typical plume size found numerically and theoretically, we derive here fully non-linear solutions of the equations (13) in the asymptotic limit Ra → ∞. As explained below, the derivation of these equations allows us to understand the balance at play in the plume formation.
In the bulk of the porous medium, since the gradients are O(Ra 1/2 ) in the x direction and O(1) in the z direction, the incompressibility condition ∂ x u+∂ z w imposes a scaling on the ratio between u and w, that is u/w = O(Ra −1/2 ). We thus introduce the rescaled variablesx andû such that: x = Ra 1/2x and u = Ra −1/2û .
With these rescaled variables, the incompressibility condition is Taking the curl of Darcy's law in (13) yields that is, θ(x, z) = w(x, z)+θ(z) to leading order in Ra. Thus, Darcy's law compels the temperature and the vertical velocity to have the same horizontal variance. Lastly, the advection-diffusion equation in (13) with rescaled variables is where all terms appear to be of the same order.
Building on our numerical results, we seek steady solutions that are periodic in the x direction. We introduce an ansatz for the flow that is the lowest order truncation of a Fourier series, that is, we assume the velocity field to have the following form, which has no mean mass flux in either vertical or horizontal direction. According to the rescaled Darcy's law (29), the temperature becomes For the flow (31) to satisfy the incompressibility condition, the following relation is required: To determine the functions w 0 and θ, we use the advection-diffusion equation (30), which contains a mean and two harmonic (k and 2k) terms. Balancing the mean terms and using the continuity relation (33) simply yields a balance between vertical heat advection and heat production, that is: The harmonick terms correspond to a balance between horizontal diffusion and the vertical advection of the average thermal energy (or temperature) profile, As noticed in the preceding section, such a balance is responsible for setting the O(Ra 1/2 ) horizontal gradients. Lastly, the harmonic 2k terms coming from the non-linear term leads to exactly the same balance as Darcy's law, which is automatically satisfied by our ansatz.
We have, therefore, constructed a fully non-linear solution which is exact in the asymptotic limit Ra → ∞, and is given by wherek and θ 0 are O(1) but a priori unknown. Note that this solution only satisfies one boundary condition: the absence of mass flux at the bottom of the porous layer. The remaining boundary conditions, be it the absence of bottom heat flux, the purely vertical velocity at the top, or either of the thermal boundary condition BC 1 or BC 2, are all unmatched with the solution.   To draw a more quantitative comparison between the non-linear solution and the flow in one plume, we plot in figure 8 several horizontal cuts at different heights of the vertical and horizontal velocities. We find that in the bulk, the theoretical solution adequately describes the amplitude of the velocity variations. However, the model becomes quite inaccurate near the upper boundary where, as noted above, it does not satisfy the correct boundary conditions. In fact, this issue seems to lead to other inaccuracies in the model: it predicts a linear decrease in the mean temperature θ, whereas the numerical simulations show a more complex dependence on z (figure 4).

Conclusion on non-linear heat transport
Throughout this section, we have detailed the properties of heat transport by convection in strongly non-linear regimes. Based on several arguments, including temperature scale, heat flux and Nusselt number measurement, we have confirmed that heat transport is dominated by advection in the bulk of the porous medium. We have carried out simulations with two different top boundary conditions that are thought to be relevant to geophysical context: one where the top boundary temperature is imposed, and another where advective heat flux is conserved in upwellings and temperature is imposed in downwellings. We have confirmed that the two boundary conditions produce the same bulk flows. Lastly, we have shown that the typical plume size follows a Ra −1/2 power law, which is due to a balance between horizontal diffusion and vertical advection of heat. It is interesting to note that internally-heated and Rayleigh-Bénard convection are different regarding the typical plume scale: in the asymptotic regimes of large Ra, Hewitt et al. (2012) found that p scaled like Ra −0.4 , which they later suggested was a result of the stability of the plumes (Hewitt & Lister, 2017).
5 Accounting for the large scale modulation of tidal heating

A simple model
In this section, we briefly explore how the large-scale variations of tidal heating affect heat transport in internally heated porous media. This is important for the case of icy satellites such as Enceladus, for which heterogeneity of tidal heating have been shown to induce focusing of the heat flux where heating is the most intense (Choblet et al., 2017). We consider here a domain with aspect ratio L = 4 for which the volume production of heat q takes the form which is such that the mean heat production is unchanged compared the homogeneous case and the maximum heat production is located at the centre of the domain. In the following, we only illustrate heat modulation with ∆q = 0.5 in the case of the boundary condition BC 2. ∆q = 0.5 is a good proxy for tidal heating which bears latitudinal and longitudinal variations by about a factor 2 between minima and maxima.

Large scale flow and pulsatility
The introduction of internal heating large-scale modulation leads to the emergence of several striking features. The first one is the attraction of plumes towards the centre where the heating is the most intense. Although plumes may exist in the whole interior of the domain, they merge towards the centre, which results in a higher temperature region with larger heat flux anomaly, as illustrated in figure 9. Plume merging towards the centre is associated with a large-scale mean flow that is also shown in figure 9. Note that, even at Rayleigh numbers as high as Ra = 3000, small-scale features persist in the timeaveraged flow: despite the strong time variability, small plumes remain locked in the areas where heat production is minimal.
Advection of the plumes towards the largest internal heating region and the subsequent plume merging leads to pulsatility in the advective heat flux, as shown in figure 10. At intermediate Rayleigh number (Ra = 360), the flux is intermittent for homogeneous heating but it exhibits a quasi-periodic behaviour for a modulated heating. The typical period is of order one, i.e. it takes place over a convective time scale, and corresponds to the time needed for plume formation, advection towards the centre and merging. The effects of heterogeneous heating are even more striking at high Ra: the steady state observed in the homogeneous case is replaced by quick oscillations of the heat flux (see figure 10). They are due to the many plumes observed in the centre of the domain reaching the top boundary non synchronously (see figure 9).

Similarities with the homogeneous-heating case
Despite the existence of a mean flow and the pulsatile behaviour detailed in the preceding section, convection with heterogeneous internal heating bears many similarities with the homogeneous case. As already noticed earlier, small scale plumes are still present in the flow, and their typical width remains proportional to Ra −1/2 (see figure  11) but with increased temporal and spatial variability. This means that the balance between horizontal diffusion, heat production and vertical advection is still at play to determine the single plume dynamics.
Moreover, even if lateral variations of the mean temperature are obvious in figure  9, the horizontally averaged temperature follows a trend that is very close to the homogeneous case, as shown in 11. This comforting observation suggests that the spatial form of heating is not particularly important for the mean dynamics and scaling laws governing heat transport in an internally heated porous medium. Core radius (h) 186 km Kinematic viscosity (ν) 1 × 10 −6 m 2 .s −2 Thermal diffusivity (κ) 8 × 10 −7 m 2 .s −2 Water thermal expansion (α) 1.2 × 10 −3 Heat capacity (c p ) 1.3 × 10 3 J.K −1 .kg −1 Gravity (g) 0.1 m.s −2 6 Application of our idealised study to the case of Enceladus

Physical properties of Enceladus' core
To characterise convection inside Enceladus' core, and to compare our results to existing literature, we use the same physical parameters as in (Choblet et al., 2017). A set of fixed physical constants that are relevant to characterise heat transport are given in table 2. We reproduce the process used in (Choblet et al., 2017) and do not precisely specify the permeability k and internal heat production Q V values. Instead, we consider that k may range from 10 −15 m 2 to 10 −12 m 2 and that the tidal heating is between 10 GW and 40 GW. (The lower bound is directly inferred from the heat flux measurement at the South Pole of Enceladus (Porco et al., 2006).) Therefore, we draw a map of the behaviour of the system keeping the parameters of table 2 constant and varying both k and Q V .

The Rayleigh number inside Enceladus
As explained in the second section of this article, the overall behaviour of the system depends only on one dimensionless parameter, the Rayleigh number, defined in (12) The 200 K isotherm gives an idea of the liquid-vapour transition that arises around 500 K at the pressure reached at the core-ocean boundary. The 100 K isotherm gives the temperature upper bound derived from geochemical measurements (Sekine et al., 2015;Hsu et al., 2015).
system is always unstable to convection, although Ra does not reach very high values and peaks around 600.

Typical velocity and temperature inside the core of Enceladus
We have shown in section 4.3 that heat transport is mostly advective, even at values of the Rayleigh number that are close to the onset of the instability. In such a regime, the dimensionless temperature is O(1), so that the temperature scale Θ * is a good proxy for the temperature difference between the ocean and the core of Enceladus. Using the scalings introduced in section 2.3 allows us to write the typical temperature scale Θ * as a function of Ra: This scale is shown in figure 12. We find that depending on the parameters, the expected temperature difference Θ * between the porous core and the ocean ranges from a few tens of Kelvin to a few hundreds. Given that the temperature of the ocean of Enceladus is at most of order 300 K, and since the boiling point of water at the pressure relevant to the bottom of Enceladus' ocean is ∼ 500 K, Θ * = 200 K sets the upper bound for the validity of the present single phase model. Note that Hsu et al. (2015) have shown via the ice plume composition that the water flowing inside Enceladus has been in contact with rocks at a temperature of about 90 • C. Our idealised model allows constraining the permeability of Enceladus' core which could range from 10 −14 to 10 −13 m 2 , in the range of tidal heating considered.
The typical velocity scale U * of the flow in the core is given by a diffusive velocity κ/h augmented by a factor Ra, that is: The diffusive velocity scale amounts to 0.1 mm.yr −1 , and because Ra does not exceed 600, the Darcy flux remains below 6 cm.yr −1 . The hydrothermal activity at the bottom of Enceladus' ocean is therefore very different from the what is commonly observed in the Earth's oceans, where typical Darcy fluxes are rather of the order of a few meters per year (10 3 times larger). This difference is largely due to the much weaker gravitational acceleration in Enceladus.
Consequently, the convective time scale τ is: The typical variability timescale, for instance for the flux at the top boundary (see figure 10) is thus at least 2 million years. It is a very slowly evolving system; as a consequence, the presently observed symmetry breaking between the south and north pole of Enceladus, where tidal heating is the same (Choblet et al., 2017), could very well be attributed to the pulsatile behaviour happening over a timescale that is too long to be appreciated by our observations.

Hydrothermal velocity at the bottom of Enceladus' ocean
To evaluate the typical velocity of the buoyant hot water coming out of the core at the bottom of the ocean we must first evaluate the buoyancy flux. Following the observation of Rabinowicz et al. (1998); Monnereau and Dubuffet (2002); Choblet et al. (2017), we assume that the plumes observed in our simulations take the form of lines when transposed into three dimensions. The typical velocity of water in the ocean produced by a line source of buoyancy flux B is U h = B 1/3 (Morton et al., 1956;Woods, 2010). The buoyancy flux is given by (Woods, 2010) where the 1d integral is computed across an upwelling zone of typical extent p ∝ Ra −1/2 . As Θ and W are proportional to Ra, B scales like Ra 3/2 , or more explicitly, Focusing of the heat flux in narrow upwelling zones leads to enhanced values of (wθ)| z=h , as shown in figure 13-a. In the case of heterogeneous heating, focusing increases by a factor 10 the heat flux at the bottom of the ocean. Finally, the typical hydrothermal velocity is found to be about 1 cm/s, no matter what the permeability or the tidal heating are (see figure 13-b). This value is in agreement with the typical velocity found by (Choblet et al., 2017) with different scaling arguments relying on the power anomaly advected to the ocean floor.

Conclusion and discussion
Throughout this article, we have explored heat transport in a fluid-saturated, internally heated porous medium with an idealised Cartesian model. Our set-up is based an idealisation of the model of Choblet et al. (2017) describing the tidally-driven hydrothermal activity in the interiors of Enceladus. The behaviour of the system is governed by a single dimensionless number, the Rayleigh number Ra, which is an increasing function of both the permeability and the internal heat production. With the combination of numerical simulations and mathematical analysis, we have derived general laws governing hydrothermal activity driven by volumetric heating.
We have shown that heat transport in the porous medium is governed by advection. In this regime, the temperature difference between the porous matrix and the pure fluid ocean scales like Ra. This scaling enables use to constrain the plausible range of values for the permeability of Enceladus' core. According to Hsu et al. (2015), the temperature scale should be at most 100 K. For values of tidal heating that are consistent with the heat flux measurement at the surface of Enceladus, our scaling indicates that the permeability should be around 10 −13 m 2 . In our simulations, we have reproduced the observation of Choblet et al. (2017) that the upwelling zones tend to narrow as the Rayleigh number is increased, and that they concentrate where internal heating is the most intense. Our simulations show that the typical plume size follows a Ra −1/2 power law, which is imposed by a balance between vertical advection and horizontal diffusion of heat. This law governing the size of heat flux anomalies at the bottom of the ocean of Enceladus compels the typical buoyancy flux injected into the ocean to be proportional to Ra 3/2 . Over the range of tidal heating and permeability that are consistent with observational data, we have found that the typical hydrothermal velocity in the ocean of Enceladus is about 1 cm/s. Despite the idealisation of our model, such an order of magnitude is consistent with that of (Choblet et al., 2017) derived from an estimate of a typical heat flux anomaly. The model used here has also helped us to highlight the underpinning of heat transport in an internally heated porous layers. In particular, we have shown that the heat-transport efficiency, which has been characterised via a generalised Nusselt number, has the same scaling as the classical Rayleigh-Bénard convection in porous media (Otero et al., 2004;Hewitt et al., 2012Hewitt et al., , 2014.
We acknowledge the highly idealised nature of our approach, but we believe our study paves the way for more refined models of the interiors of Enceladus. A significant simplification is the assumption of a uniform gravity field. However, accounting for a linear vertical increase in g seems unlikely to result in significant differences. In particular, the balance between heat production, vertical advection and horizontal diffusion would still hold with a space-dependent gravity field, so the qualitative behaviour of the system should remain similar. Another important simplification is assuming homogeneous and isotropic permeability k. Since the core of small icy satellites such as Enceladus' is an aggregate of heterogeneous material, their permeability is unlikely to be uniform. It is not even clear whether coarse-grained modelling based on the assumption of strong confinement, that is, Darcy's law, is entirely relevant for the core of icy satellites, although, as we've seen, Darcy's law with small permeability is consistent with observed data. However, we do not believe any significant progress can be achieved in these directions without further constraining the core's small-scale structure.
In addition, the model we consider completely discards flows that are directly driven by the periodic tidal distortion. Although the tidal deformation field is purely incom-pressible in continuous, mean flows analogous to Stokes drift may result from the periodic motion of the porous matrix. Whether deformation-driven flows are comparable to buoyancy-driven flows remains to be quantified.
Lastly, there is a need to clarify the behaviour of the system at the top of the porous core and the coupling between the porous layer and the above ocean. We have stated in the second section that the two possible thermal boundary conditions used here (imposed temperature of free temperature in the upwellings) are the two end-members of the behaviour of the fluid at the interface. The imposed temperature condition could represent a very slow porous layer lying underneath a very well mixed ocean. Such a discrepancy is suggested by our estimates of the Darcy flux (∼ 1 cm.yr −1 ) compared to the hydrothermal velocity (∼ 1 cm.s −1 ). In this configuration, the water coming out of the core is at the same temperature as the ocean and is neutrally buoyant; there is then no hydrothermal activity in the sense of what we know at the bottom of the Earth's ocean. Nevertheless, it is associated to a diffusive heat flux anomaly on the subsurface ocean's floor which is likely to drive convection and mixing in the ocean. The observed chemical signature of contact with silicate rocks at high temperature (Hsu et al., 2015) could very well happen below the thin thermal boundary layer at the top of the core. Moreover, current thermal evolution models of icy moons rely on parametrisation of hydrothermallydriven convection in the sub-surface ocean that are based on the classical Rayleigh-Bénard problem (Travis et al., 2012;Travis & Schubert, 2015). It is, however, not clear at all whether such parametrisation actually applies to the present system where ocean convection is driven by strong and localised heterogeneities of either the advective of the diffusive heat flux at the bottom boundary In short, it remains difficult to produce a definitive statement about the thermal structure of the subsurface ocean without a careful study of the coupled system with two very different typical evolution timescales for each medium.