Laëtitia Lebec
email: laetitia.lebec@ens-lyon.fr

Stéphane Labrosse

Adrien Morison

Paul J Tackley

Scaling of convection in high-pressure ice layers of large icy moons and implications for habitability

Keywords:

come

Introduction

The internal structure of large icy moons, as Ganymede or Titan, differs from that of smaller ones as Enceladus or Europa for which a direct contact exists between the core and the ocean under the icy surface [START_REF] Hussmann | 10.18 -Interiors and Evolution of Icy Satellites[END_REF]. Their composition has been studied from the data of Cassini-Huygens which explored Saturn's moons from 2004 to 2017 and Galileo missions which explored Jupiter's moons from 1995 to 2003. Currently, Juno is orbiting Jupiter and performed one close flyby of Ganymede [START_REF] Ravine | Ganymede Observations by JunoCam on Juno Perijove 34[END_REF], which will also be one of the targets of the JUICE mission [START_REF] Grasset | JUpiter ICy moons Explorer (JUICE): An ESA mission to orbit Ganymede and to characterise the Jupiter system[END_REF] around 2030 in order to better constrain several parameters of this moon. For large icy moons, as Ganymede, Callisto and Titan, owing to the phase diagram of water, a high-pressure (HP) ice layer of tetragonal ice VI and ice V, depending on the assumed thickness of the ice mantle [START_REF] Vance | Geophysical Investigations of Habitability in Ice-Covered Ocean Worlds[END_REF], may exist between the internal liquid ocean and the rocky core [START_REF] Hussmann | 10.18 -Interiors and Evolution of Icy Satellites[END_REF], in addition to the ice Ih at the surface of the moon. The consequences of this internal structure is that the core and the ocean are not in direct contact, which has long been considered as an obstacle to their habitability by strongly limiting or even removing all interactions between these two layers. Efficient exchanges of nutrients between the core and the ocean being one of the necessary conditions for life to appear at high-pressure in a deep water ocean [START_REF] Picard | Pressure as an environmental parameter for microbial life -A review[END_REF], these moons were considered less favorable to habitability than smaller moons like Europa [START_REF] Noack | Water-rich planets: How habitable is a water layer deeper than on Earth?[END_REF].

In recent years, several studies showed that, under certain conditions depending on the thickness of the HP ice layer, the viscosity of the ice and the heat flux from the core, periodically indirect contacts, or even direct contacts, could still be possible by convection through the HP ice layer. [START_REF] Choblet | Heat transport in the high-pressure ice mantle of large icy moons[END_REF] described a 3D spherical model of the HP ice shell including partial melting at the bottom of the ice layer, in contact with the core, the melt being extracted instantly to the ocean. Kalousová and Sotin (2018) and Kalousová et al. (2018) described a 2D cartesian two-phase convection model including partial melting and melt transport through the HP ice layer. These studies consider the pressure and temperature dependence of the viscosity. In addition, Kalousová and Sotin (2018) and Kalousová et al. (2018) also include viscosity variations owing to the presence of melt. All those studies concluded that for various values of ice viscosity, shell thickness and heat flux from the core, mass transfers is possible between the core and the ocean through the HP ice layer by convection in a two-phase medium. A key aspect common to these models is the possibility of melting of the HP ice at the rock interface and transfer trough the ice layer, as liquid channels or by solid-state convection.

While these previous papers consider the possibility of melting at the bottom of the HP ice layer and partial melting in the bulk, with various levels of approximation, they do not consider the dynamical effect of the solid-liquid phase equilibrium at the upper boundary of the HP ice layer. Indeed, several studies showed that a melting-freezing boundary drastically changes the condition for convection in a solid layer and the resulting convection characteristics at finite amplitude [START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF][START_REF] Labrosse | Rayleigh-Bénard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries[END_REF][START_REF] Agrusta | Mantle convection interacting with magma oceans[END_REF][START_REF] Morison | Convection in the primitive mantle in interaction with global magma oceans[END_REF]. Applied to our study, the interface between the liquid ocean and the HP ice layer can be set as a phase change boundary condition. Physically, convective stresses in the solid create a topography of the interface which can be erased by melting and freezing. It implies a non-zero radial velocity at the surface of the HP ice layer, which significantly increases the efficiency of the mass transfer through the HP ice layer. This effect can be modeled as a boundary condition for the ice layer controlled by the phase change coefficient Φ [START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF]Labrosse et al., 2018, and references therein). The aim of this paper is to evaluate the implications of the liquid-solid phase change boundary condition for convection in high pressure ice layer of large icy satellites and water planets.

To this end, we used a numerical model solving the convection in the HP ice layer subject to the phase change boundary condition. We explored systematically the numerical solution depending on the dimensionless control parameters (mainly the Rayleigh number and the phase change number) and obtained scaling laws relating output parameters (bottom temperature, flow velocity) as function of these input dimensionless numbers. We show in this paper how these scaling laws can be used to make prediction on specific planetary objects, like Titan and Ganymede.

We do not consider two-phase flow in the bulk of the ice layer, as done in some previous studies [START_REF] Choblet | Heat transport in the high-pressure ice mantle of large icy moons[END_REF]Kalousová et al., 2018;Kalousová and Sotin, 2018) to isolate the effect of the phase change boundary condition. We nevertheless compute, as a post-treatment, the amount of partial melting that we would expect if it was permitted, which allows us to compare our results to these previous studies.

In section 2 we present the physical and numerical models used to solve the convection problem. The results are described in section 3. These generic results are then applied to several bodies in subsection 3.3. In section 4 we discuss the limitations of the results and the possible future work. A summary and conclusion are provided in section 5. Model illustration for the interior of an ocean world with an HP ice layer in direct contact with a liquid ocean at top and a rocky core/mantle at the bottom (type IV or V of [START_REF] Lammer | Origin and Evolution of Planetary Atmospheres -Implications for Habitability[END_REF]type H3 of Noack et al., 2016). Our numerical models treat convection in the HP ice layer.

Model

Physical model

Governing equations

We want to quantify the heat and mass transfer by convection in HP ice layers of various planetary objects, from large icy satellites like Ganymede and Titan, to water exoplanets. The large number of control parameters with wide ranges of possible values among these objects can be reduced by use of dimensionless equations and numbers and the principle of physical similarity [START_REF] Barenblatt | Scaling, Self-Similarity, and Intermediate Asymptotics[END_REF]. We use the thickness 𝑑 of the layer as length scale, the diffusion time 𝑑 2 ∕𝜅 as time scale, with 𝜅 the thermal diffusivity. For an ice layer subjected to an imposed heat flux 𝑞 at the bottom, the relevant temperature scale is 𝑞𝑑∕𝑘, 𝑘 being the thermal conductivity. Because we consider an incompressible model using the Boussinesq approximation, the actual value of the temperature T (dimensional) is irrelevant to the dynamics and we set the surface temperature T𝑡𝑜𝑝 to be equal to the melting temperature. The dimensionless temperature is

𝑇 = 𝑘 T -T𝑡𝑜𝑝 𝑞𝑑 . ( 1
)
With the dynamic viscosity 𝜂, which is considered constant with pressure and temperature through the HP ice layer in our model, the thermal expansion coefficient 𝛼, the gravity acceleration 𝑔, the shell inner radius 𝑅 -(the upper radius being 𝑅 + = 𝑅 -+ 𝑑, see fig 1), the reference density 𝜌, we complete the set of input parameters except for that relevant to the phase change boundary condition that is introduced below ( § 2.1.2). Table 2 gives numerical values of these parameters for Ganymede and a water exoplanet. As is well known (e.g. [START_REF] Ricard | 02 -physics of mantle convection[END_REF], Rayleigh-Bénard convection in the Boussinesq approximation (and with the previous variables considered constant) is controlled by two dimensionless parameters (aside from the one associated with the phase change boundary), the Rayleigh number,

Ra 𝑞 = 𝛼𝑔𝑞𝜌𝑑 4 𝑘𝜅𝜂 , ( 2
)
and Prandtl number,

𝑃 𝑟 = 𝜂 𝜌𝜅 . ( 3
)
The definition of the Rayleigh number given in eq. 2 differs from the classical one obtained for an imposed temperature difference and is relevant to the case of an imposed heat flux [START_REF] Choblet | Heat transport in the high-pressure ice mantle of large icy moons[END_REF], as considered here.

For the solid HP ice layer, the Prandtl number is large enough (of order 1 × 10 20 , see table 2) to be considered infinite, meaning the inertia of the ice can be neglected in front of viscous forces. Then, the dimensionless conservation equations of mass, momentum and energy for thermal convection in the spherical HP ice shell under the Boussinesq approximation for a infinite Prandtl number are the following (e.g. [START_REF] Ricard | 02 -physics of mantle convection[END_REF]:

𝛁 ⋅ 𝒖 = 𝟎, (4)
0 = -𝛁𝑝 + 𝛁 2 𝒖 + Ra 𝑞 (𝑇 -T )r, (5)
𝜕𝑇 𝜕𝑡 + 𝒖 ⋅ 𝛁𝑇 = 𝛁 2 𝑇 , (6)
with 𝒖 = (𝑣, 𝑤) the velocity, 𝑝 the dynamic pressure, 𝑇 and T the temperature and steady-state conductive temperature (see section 2.2.1) and r the radial unit vector.

Solid/Liquid phase change at top boundary

The existence of a solid/liquid phase equilibrium at the top of the ice layer leads to a specific mechanical boundary condition in place of the classical non-penetrative one usually considered in convection models. This boundary condition and its implications for Rayleigh-Bénard convection has been the subject of a few previous papers [START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF][START_REF] Deguen | Thermal convection in Earth's inner core with phase change at its boundary[END_REF][START_REF] Labrosse | Rayleigh-Bénard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries[END_REF][START_REF] Agrusta | Mantle convection interacting with magma oceans[END_REF][START_REF] Bolrão | Timescales of chemical equilibrium between the convecting solid mantle and over-and underlying magma oceans[END_REF] where interested readers will find all the details of the derivation of this boundary condition. Here, we only recall the basic ideas leading to this boundary condition and its implications.

At the interface between the HP ice layer and the ocean, a phase change can occur in either direction, melting or freezing. The secular evolution of the planetary object leads to a net motion of the interface, which we do not consider here since it would require a full thermal evolution model that goes beyond the scope of this paper. On the other hand, the ice flowing vertically toward the interface can cross it by melting and, conversely, liquid water can solidify and accrete to the HP ice layer above down-welling currents. In practice, this happens because the convective stresses in the ice layer generate a topography of the ice-water interface. Assuming thermal equilibrium at the interface makes its temperature equal to that of melting everywhere, but at different depth. This is equivalent to creating lateral temperature differences at the same depth in the liquid layer. Convection in the liquid layer act against these lateral variations to erase them which in turn tends to erase the interface topography by melting and freezing. This requires transporting latent heat from regions of freezing that provide it to regions of melting that consume it. The behaviour of the boundary depends therefore on the competition between two processes, topography building by convection in the ice layer and its erasing by convective transfer in the liquid. Each process has its own timescale and their ratio decides which wins.

The timescale to build the topography at the interface is dictated by the equilibrium between the viscous stress in the ice layer and the weight of the topography and is

𝜏 𝜂 = 𝜂 Δ𝜌𝑔𝑑 , ( 7
)
with Δ𝜌 = 𝜌 𝑠 -𝜌 𝑙 the density difference across the interface, 𝜌 𝑠 and 𝜌 𝑙 being the density of the solid and the liquid, respectively. Using parameter values for Ganymede given in table 2 gives a timescale of order 1 to 125 yr depending on the value of 𝜂.

The timescale for erasing the topography is associated to convection in the liquid layer, with a typical flow velocity 𝑢 𝑙 , that transports latent heat 𝐿. A detailed analysis of this process [START_REF] Labrosse | Rayleigh-Bénard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries[END_REF][START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF][START_REF] Deguen | Thermal convection in Earth's inner core with phase change at its boundary[END_REF] gives

𝜏 𝜙 = 𝜌 𝑠 𝐿 𝜌 𝑙 𝑐 𝑝𝑙 𝑢 𝑙 | | | 𝜕𝑇 𝑚 𝜕𝑟 | | | , ( 8
)
with 𝑐 𝑝𝑙 the specific heat of the liquid. This timescale is difficult to estimate but assuming a flow velocity of about 3 × 10 -3 m s -1 (Computed using [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF]; Soderlund (2019)) and using parameter values from table 2 gives a timescale of order 15 yr. Note that this theory assumes a turbulent flow in the ocean, 𝑢 𝑙 being its typical RMS value. Considering 𝜏 𝜙 to be constant in time and uniform in space amounts to two assumptions. The typical timescales for the fluctuations in the liquid are assumed small compared to that for the dynamics in the solid. In addition, the turbulence is assumed isotropic. This second assumption could be relaxed by considering a laterally varying value of 𝜏 𝜙 to take into account the effect of rotation on the dynamics of the ocean. However, a proper theory is still lacking for this type of effect to be included.

The ratio of the two timescales,

Φ = 𝜏 𝜙 𝜏 𝜂 , ( 9
)
is called the phase change number and controls the behaviour of the boundary. For Ganymede, from the previous calculations, Φ would be in the range [0.1; 15]. This value is a very rough estimate and this range could be larger. The key point here is that the phase change at the top interface must be efficient. A detailed analysis combining the traction continuity across the boundary, the energy balance associated with latent heat and thermodynamic equilibrium at the boundary leads to a single boundary condition for the radial velocity which, in dimensionless form, is [START_REF] Labrosse | Rayleigh-Bénard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries[END_REF][START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF][START_REF] Deguen | Thermal convection in Earth's inner core with phase change at its boundary[END_REF][START_REF] Bolrão | Timescales of chemical equilibrium between the convecting solid mantle and over-and underlying magma oceans[END_REF])

Φ𝑤 + 2 𝜕𝑤 𝜕𝑟 -𝑝 = 0, (10)
with 𝑤 the radial velocity. The vertical velocity varies along the boundary, depending on the convection pattern, as can be seen on fig 2.c, but the mean radial velocity at the interface is null due to mass conservation. For this reason, in the following, we use the root-mean-square (RMS) of 𝑤 to evaluate the efficiency of the mass transfer at the top boundary.

Varying the value of Φ between 0 and ∞ makes the behaviour of the boundary evolve between end-members. In the limit case of Φ → ∞, the radial velocity at the boundary, 𝑤 𝑡𝑜𝑝 must tend to 0 (fig 3.b) and we recover the classical non-penetration boundary condition used in all previous studies of convection in icy satellites. This happens if heat transfer is inefficient in the ocean so that the weight of the topography limits its building and therefore the radial motion, as shown on fig. 2.a. On the other hand, for Φ → 0 fig. 2.b-c, the radial velocity is unconstrained and its radial gradient is set. This happens if the topography is erased by phase change faster than it is built so that its weight is never limiting radial flow. The boundary is then permeable, which is known to drastically change the convective regime (e.g. [START_REF] Ricard | Lifting the cover of the cauldron: Convection in hot planets[END_REF][START_REF] Monnereau | Is Io's mantle really molten[END_REF][START_REF] Agrusta | Mantle convection interacting with magma oceans[END_REF] compared to the other end-member that is usually considered. With the estimates discussed above, we expect a rather low value of the phase change parameter but we explore systematically the effects of changing its value, in the range [10 -2 , ∞), since it is rather ill-constrained.

Other boundary conditions

As can be seen on figure 1 and described by equation 2, our model is based on a fixed flux at the bottom boundary between the core and the HP ice layer and a fixed temperature at the upper boundary between the HP ice layer and the liquid ocean, which is set to 0 by rendering the temperature dimensionless. It leads to the following temperature boundary conditions:

𝑇 + = 0, (11)
( 𝜕𝑇 𝜕𝑟 ) - = -1, (12)
the + and -exponents referring to, respectively, the top and bottom boundaries. As the ice VI is in direct contact with a rocky core/mantle whose viscosity is orders of magnitudes larger than that of ice, we normally impose a non-penetrative no-slip (rigid) condition:

𝑢 -= 𝑤 -= 0. (13)
However, since the heat flux is fixed at the bottom, the temperature 𝑇 -varies along the bottom boundary, allowing, in some cases, pockets or even a global film of melt to form at the interface between the core and the HP ice layer, which could lubricate the boundary and lead to a free-slip BC, which would imply that:

𝜕 2 (𝑟𝑤) 𝜕𝑟 2 = 0. (14)
In practice, where melt is present at the bottom, the situation is likely intermediate between the free-slip and no-slip boundary condition, depending on the thickness of the melt layer and the roughness of the interface. For the sake of simplicity, we consider the two end-member situations of a free-slip or no-slip bottom boundary condition.

Numerical method

The first step of this study, for each choice of input parameters as listed in table 1, consists in finding the critical Rayleigh number Ra 𝑐 for the onset of convection and the associated mode. For a Rayleigh number lower than Ra 𝑐 , the heat transfer within the HP ice layer occurs only by conduction. The conditions for the onset of convection are computed using a linear stability analysis described in subsection 2.2.1. Having determined the critical Rayleigh number, we use a mantle convection code to compute finite amplitude solutions, as explained in §2.2.2.

Linear stability analysis

The convection equations and their boundary conditions always admit a steady motionless solution in which heat is transferred by conduction only. In the present case, the steady conduction temperature profile is:

𝑇 = 𝑟 -2 𝑟 -𝑟 -𝛾, ( 15
)
where 𝛾 = 𝑟 -∕𝑟 + , 𝑟 + (resp. 𝑟 -) the dimensionless radius of the upper (resp. lower) boundary of the HP ice layer, with 𝑟 + -𝑟 -= 1 by our choice of the thickness of the layer as length scale to render the equations dimensionless.

Using StabLinRB (https://github.com/amorison/stablinrb), a free computational tool for linear stability analysis developed by Stéphane Labrosse and Adrien Morison [START_REF] Labrosse | Rayleigh-Bénard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries[END_REF][START_REF] Morison | Timescale of overturn in a magma ocean cumulate[END_REF][START_REF] Morison | Convection in the primitive mantle in interaction with global magma oceans[END_REF], we looked for the stability of the conductive solution for each study case. Infinitesimally small perturbations are written as spherical harmonics with coefficient whose radial dependency is written as a sum of Chebyshev polynomials.

Using differentiation matrices for this Chebyshev collocation approach, the linearized problem for each mode is transformed into a generalized eigenvalue problem. The critical Rayleigh number is the one that makes the real part of the eigenvalue change sign and the mode which has the minimum critical Ra is the most unstable one. You can refer to the studies of [START_REF] Morison | Convection in the primitive mantle in interaction with global magma oceans[END_REF]; [START_REF] Labrosse | Rayleigh-Bénard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries[END_REF]; [START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF] and [START_REF] Deguen | Thermal convection in Earth's inner core with phase change at its boundary[END_REF] for further details on the method used in our study for linear stability analysis and associated spherical harmonics.

Finite amplitude models

The convection equations (4-6) subject to the boundary conditions described in § 2.1.2 and 2.1.3 are solved using the mantle convection code StagYY [START_REF] Tackley | Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid[END_REF] in the spherical annulus two-dimensional geometry [START_REF] Hernlund | Modeling mantle convection in the spherical annulus[END_REF]. The runs have been divided into four main categories, low super-criticality (Ra 𝑐 ≤ Ra 𝑞 ≤ 10Ra 𝑐 ) or high super-criticality (Ra 𝑞 > 10Ra 𝑐 ) for both free-slip and rigid mechanical boundary condition at the bottom. In each category, values of the phase change parameter Φ was varied systematically in the range 10 -2 to ∞ (no phase change), with two values of the aspect ratio, 𝛾 = 0.9 and 0.95 (see table 1). The runs are pursued until a statistically steady-state is reached, which is the case when all global diagnostics (mean temperature, heat fluxes at both boundaries, RMS velocity) only fluctuate around a time-independent mean. The chosen grid size depends on the simulation study-case and boundary conditions, in the form of 𝑛𝑦 𝑡𝑜𝑡 × 𝑛𝑧 𝑡𝑜𝑡 , being the number of points in the HP ice shell respectively in the horizontal and vertical directions. The resolution is deemed sufficient when the global energy balance is satisfied.

We performed calculations with increasing values of 𝑛𝑧 𝑡𝑜𝑡 until this criterion was satisfied, with 𝑛𝑦 𝑡𝑜𝑡 chosen such that the cells are nearly square (detailed values can be found in the table captions of the supplementary material (see appendix A)). A vertical grid refinement depending on the study case is also applied at top and bottom boundaries. The geometry of the problem is adapted to the need of each simulation in order to capture the relevant flow structure: the convection equations can be solved either on the entire ice shell or only a part of it (see § 3.2 for detailed explanations).

As shown by [START_REF] Agrusta | Mantle convection interacting with magma oceans[END_REF], the phase change boundary condition drastically changes the behaviour of the boundary layer (and therefore the convective dynamics) compared to the classical non-penetration BC. While in the classical situation, plumes originating from the opposite boundary have to turn when reaching a non-penetration boundary, which builds a boundary layer, they can here flow directly through the boundary by melting. This process still requires to adjust the temperature to that of the boundary but this happens on a distance much smaller than the thickness of a classical boundary layer. The thickness of a classical boundary layer is set by its stability, and typically scales as Ra -1∕3 for a fixed temperature free-slip situation. On the other hand, the thickness of a phase-change boundary layer is simply set by the rate of flow toward it 𝛿 ∼ 1∕𝑤 𝑡𝑜𝑝 and we will show below that it therefore scales as Ra

-1∕2 𝑞 .
This boundary layer is not driving any dynamics but needs to be resolved to balance heat correctly, which rapidly becomes prohibitive (in terms of grid spacing and therefore time-step) at high Ra 𝑞 .

For this reason, while the fixed temperature boundary condition (eq. 11) is applied at the top boundary for runs close to the critical value of the Rayleigh number, we adopted another thermal boundary condition for cases with a high Rayleigh number and small value of Φ, as initially introduced by [START_REF] Agrusta | Mantle convection interacting with magma oceans[END_REF]. The basic idea is to remove the very thin boundary layer from the calculation and apply a boundary condition relevant to what happens physically on the edge of that boundary layer. Mathematically, it is described as an intermediate condition between

Dirichlet and Neumann boundary conditions, which respectively corresponds to a fixed temperature or flux condition.

The intermediate Robin boundary condition is defined by the following equation [START_REF] Agrusta | Mantle convection interacting with magma oceans[END_REF]:

Γ𝜃 + (1 -Γ) 𝜕𝜃 𝜕𝑟 = 0, ( 16
)
where 𝜃 is the lateral deviation of temperature compared to the mean and Γ is an approximation of the Heaviside function:

Γ = 1 2 ⎡ ⎢ ⎢ ⎢ ⎣ 1 + tanh ⎛ ⎜ ⎜ ⎝ 𝜋 𝑤 0 2 -𝑤 𝑡𝑜𝑝 𝑤 0 2 ⎞ ⎟ ⎟ ⎠ ⎤ ⎥ ⎥ ⎥ ⎦ (17)
with 𝑤 0 a threshold velocity. In our study, we used 𝑤 0 = 5.10 -2 Ra 2 3

𝑞 , in order to be about half the typical scaling of the RMS velocity in convection [START_REF] Agrusta | Mantle convection interacting with magma oceans[END_REF]. This function allows the boundary condition to switch smoothly between the Neumann and Dirichlet options, depending on the flow velocity. For rapidly up-welling flows, 𝑤 𝑡𝑜𝑝 > 𝑤 0 , Γ ≃ 0 and we apply 𝜕𝜃 𝜕𝑟 = 0, which is expected for a hot plume rapidly melting in contact with the liquid. For downwelling or slowly upwelling ones, 𝑤 𝑡𝑜𝑝 ≤ 𝑤 0 , the normal fixed temperature condition can be imposed, 𝜃 = 0.

The validity of this approach has been checked by comparing the results to high resolution calculations using the Dirichlet BC (eq. 11).

Results

Onset of convection

Let us first discuss the results of the linear stability analysis that give the conditions for the onset of convection.

The figure 4 shows the convective solution at onset (linear) in the HP ice layer for various values of Φ and different bottom boundary conditions. The linear problem depends on the degree of spherical harmonics but is degenerate in terms of order. For ease of representation and for comparison with the spherical annulus geometry used for the finite amplitude calculations, we represent the sectoral harmonics (i.e. with order equal to degree) in the equatorial plane.

For a free-slip BC the critical Rayleigh number Ra 𝑐 and the number of convective rolls (spherical harmonics degree at onset) 𝑙 𝑐 of the first unstable mode are always lower than for a rigid BC. Also, both Ra 𝑐 and 𝑙 𝑐 increase with Φ (fig 5).

It is particularly interesting to note that the most unstable mode for a free-slip BC and small values of Φ (Φ ≲ 4) is a degree 1 mode. A similar result was obtained in previous studies with an imposed temperature at the bottom [START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF][START_REF] Morison | Timescale of overturn in a magma ocean cumulate[END_REF] and the present calculations extend that finding to the situation of an imposed bottom heat flux.

The impact of the boundary conditions can also be observed directly on the shape of flow lines. For small values of Φ, the phase change at the top is efficient and the flow lines are open at the top interface, while they are closed when Φ → ∞. Also, for a rigid BC at the bottom and a free-slip one at the top, the convective circulation is concentrated near the top interface since horizontal motion is not limited by that type of boundary.

Then, left and middle panels of figure 5 show respectively the evolution of the critical Rayleigh number Ra 𝑐 and critical wave-number 𝑙 𝑐 with the phase change number at the top boundary, Φ, depending on the bottom BC and the aspect ratio 𝛾 for the HP ice shell. In all cases the critical Rayleigh number for the onset of convection increases with Φ, with important variations in the range Φ = 1-100. For a free-slip BC the convection starts for lower values of Ra 𝑐 compared to a rigid BC, meaning that the convection is easier to start in that configuration. While the effect of the aspect ratio 𝛾 on Ra 𝑐 is limited, it is important for the wave-number, which is increasing with 𝛾, except for a free-slip BC with Φ ≲ 4. The effect of the aspect ratio on the wave-number for classical non-penetrating BC is well documented [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF] and is related to the linear stability in the case of planar layers. The large values of 𝛾 considered here makes the shell close to being planar but changing 𝛾 modifies the number of convective rolls of a given aspect ratio (the one obtained for plane layers) that can fit in. The degree 1 mode obtained for low values of Φ and a free-slip bottom boundary condition can be understood when considering the flow presented on figure 4. In this case, down-welling occurs on the hemisphere where freezing happens, while up-welling occurs on the other hemisphere, where melting occurs. The return flow entirely happens in the liquid state, which is not treated explicitly. Such a flow in the solid minimizes the amount of deformation and, therefore, viscous resistance. For this reason, its critical Rayleigh number is much smaller than the one obtained at large values of Φ. This mode requires free-slip around the core and is suppressed when using rigid boundary conditions.

Applying the linear stability analysis to a specific ocean world as Ganymede, the critical thickness 𝑑 𝑐 of the HP ice layer for the onset of convection can be computed. The right panel of figure 5 shows the evolution of this parameter with Φ, depending on the bottom BC and the aspect ratio 𝛾, using values listed in table 2, when considering reference values for the heat flux from the core 𝑞 = 10 mW m -2 and for the viscosity of the ice 𝜂 = 10 15 Pa s. Then, to start the convection on Ganymede considering these reference values, the minimum thickness has to be in 1.5-2.1 km range for a rigid BC and in 1-1.7 km range for a free-slip BC. For the reference HP ice thickness of 100 km, the Rayleigh is close to Ra 𝑞 = 5 × 10 8 (see eq. 28), which is far above the critical value to onset the convection. However, in the context of planetary evolution, the thickness of the HP ice layer is expected to evolve and it is interesting to consider a full range, from the critical value to the large nominal value just mentioned. This is done using the finite amplitude convection code StagYY ( §2.2.2).

Finite amplitude calculations and scaling laws

We ran calculations systematically exploring a wide range of values of the main dimensionless input parameters, the flow is dominant and the simulation has to be solved on the entire ice shell. The fact that the dominant degree is not 1 but rather a combination of 2 and 3, contrary to the prediction of the linear stability analysis, can be explained by the initial conditions chosen for this calculations, a random white noise added to the steady conduction solution. Even though the degree 1 mode has a larger linear growth rate, if other modes with only slightly lower growth rates have a larger initial amplitude, they can take over and dominate in the non-linear calculation. When Ra 𝑞 = 2Ra 𝑐 (right), the solution is already dominated by higher degree perturbations. In that case, as in the ones for Φ ≥ 10 or a rigid BC for which the critical harmonics 𝑙 𝑐 is at least around 10 (see fig 5), the simulations can be performed on a part of the shell to reduce the computational resources (detailed values can be found in the table captions of the supplementary material (see appendix A)).

The first targets of our study concern the mass flux between the ice layer and the overlying ocean and heat transfer efficiency by convection in the ice layer. The former is quantified by the root-mean-square (RMS, the mean being null by virtue of mass conservation) of the radial velocity at the top surface, 𝑤 𝑡𝑜𝑝 , while the latter is measured by the dimensionless mean temperature difference across the ice layer, Δ T , or its inverse, which, in our case of an imposed heat flux at the bottom, is precisely the Nusselt number, i.e. the ratio between the heat flux and what would be transported by steady-state conduction, Nu = 1∕Δ T . In the following, we consider the lateral mean temperature T , because as can be seen on fig. 3a, its variation along the shell can be important. Fully detailed results tables and figures corresponding to the exploration of the parameters space for all cases described in table 1 can be found in the supplementary material (see appendix A).

Starting with the low super-criticality cases (Ra 𝑞 ≤ 10Ra 𝑐 ), figure 7 shows the evolution of both diagnostics, 𝑤 𝑡𝑜𝑝 and Nu, as function of Ra 𝑞 for both choices of bottom BC and various values of Φ. As expected, the values of 𝑤 𝑡𝑜𝑝 and Nu rapidly increase with Ra 𝑞 from their values at onset (Ra 𝑞 = Ra 𝑐 ) of convection, 0 and 1, respectively. The increase rate with Ra 𝑞 gets larger when smaller values of Φ are considered. This behaviour is similar to the one obtained in a plane layer with different boundary conditions [START_REF] Agrusta | Mantle convection interacting with magma oceans[END_REF].

Figure 8 shows the global diagnostics 𝑤 𝑡𝑜𝑝 and Nu for both choices of bottom BC and various values of Φ for larger values of Ra 𝑞 (SC3 and SC4, see table 1). Since these diagnostics and Ra 𝑞 vary on several orders of magnitude, the plots use log-log scales. As expected, the values of 𝑤 𝑡𝑜𝑝 and Nu increase with Ra 𝑞 as power laws. We can also see that, for a given value of Ra 𝑞 , 𝑤 𝑡𝑜𝑝 and Nu increase when decreasing Φ, as the top boundary becomes progressively more permeable. The evolution of the Nusselt number saturates at some point which seems to depend on the value of Ra 𝑞 and the mechanical boundary condition at the bottom: for a rigid BC, the value for Φ = 100 is similar to that for lower values at large values of Ra 𝑞 but is intermediate for small values of Ra 𝑞 . For this reason, the exponent of the Nu = 𝑓 (Ra 𝑞 ) law is different for Φ = 100 than for other choices of that parameter. Each combination of choice of BC and Φ value leads to a power law relationship of the form:

𝑤 𝑡𝑜𝑝 = 𝑎 𝑤 Ra 𝛽 𝑤 𝑞 , ( 18
)
Nu = 1 Δ T = 𝑎 𝑇 Ra 𝛽 𝑇 𝑞 , (19)
(𝑎 𝑤 , 𝛽 𝑤 ) being the scaling coefficients for the top radial velocity and (𝑎 𝑇 , 𝛽 𝑇 ) being those for the Nusselt number. These scaling coefficient depend on the value of the phase change number Φ and the mechanical boundary condition at the bottom.

As can be seen on figure 8, the exponents 𝛽 𝑥 depend on the choice of boundary condition but weakly on Φ. On the other hand, the coefficients 𝑎 𝑥 depend strongly on Φ, as shown on fig. 9 on which 𝑎 𝑤 is plotted as function of Φ.

In the case of a rigid BC, we get 𝛽 𝑤 = 0.47-0.53, i.e. close to 1∕2, and 𝛽 𝑇 = 0.2-0.22, i.e. close to 1∕5. In the case of a free-slip BC, we obtain 𝛽 𝑤 = 0.45-0.55 which is also similar to 1∕2 but 𝛽 𝑇 = 0.22-0.24 is close to 1∕4. The differences in exponents between the free-slip and no-slip BCs have been explained with some theoretical models in a different context, notably using a fixed boundary condition on both boundaries (e.g. [START_REF] Roberts | Fast viscous bénard convection[END_REF]. Adapting that theory to the present situation falls beyond the scope of this paper but we can explain in a simpler way how 𝛽 𝑇 can be justified.

The 𝛽 𝑇 exponents can be obtained theoretically using the argument of the stability of boundary layers [START_REF] Malkus | The heat transport and spectrum of thermal turbulence[END_REF]: the thickness of the boundary layer is set by its stability, therefore its Rayleigh number, Ra 𝛿 . The existence of a boundary layer at the top interface depends on the value of Φ. As shown on the average temperature profiles on figure 3a, when Φ → 0, the boundary layer disappears at the top interface and only one boundary layer of thickness 𝛿 ′ has to be considered at the bottom interface, while there are two boundary layers of thickness 𝛿 ′ in the case Φ → ∞. For a free-slip BC at the bottom interface, when Φ → 0,

Ra 𝛿 = Ra 𝑞 Δ T 𝛿 3 = 𝐴 (20)
with Ra 𝛿 the Rayleigh number of the boundary layer, 𝛿 = 𝛿 ′ ∕𝑑 the dimensionless thickness of the boundary layer and 𝐴 a constant related but not exactly equal to the critical Rayleigh number for the instability of the boundary layer [START_REF] Howard | Convection at high Rayleigh number[END_REF][START_REF] Sotin | Three-dimensional thermal convection in an iso-viscous, infinite Prandtl number fluid heated from within and from below: applications to the transfer of heat through planetary mantles[END_REF]. As 𝑞 ′ = 𝑘 Δ T ′ 𝛿 ′ , primed variables being dimensional, is used to scale temperature, its With this simple argument, we obtain the exponent 𝛽 𝑇 = 1∕4. The coefficient 𝑎 𝑇 = 𝐴 -1 4 cannot easily be predicted from first principles [START_REF] Sotin | Three-dimensional thermal convection in an iso-viscous, infinite Prandtl number fluid heated from within and from below: applications to the transfer of heat through planetary mantles[END_REF].

On the other hand, when Φ → ∞, the coefficient 𝛽 𝑇 should be the same, but the variation of the 𝑎 𝑇 coefficient can be computed in order to better understand the figure 8. The derivation is similar but the temperature jump is now split between the two boundary layers. Neglecting the effect of curvature for simplicity, since the aspect ratio of the shell is large (see [START_REF] Sotin | Three-dimensional thermal convection in an iso-viscous, infinite Prandtl number fluid heated from within and from below: applications to the transfer of heat through planetary mantles[END_REF], for a more complete derivation), we can assume that each boundary layer takes an equal share of the total temperature jump, therefore

𝑅𝑎 𝛿 = Ra 𝑞 Δ T 2 𝛿 3 = 𝐴 (22)
and the dimensionless heat flux from the core is 1 = Δ T 2𝛿 . Then, 10 2 10 1 10 0 10 1 10 2 10 3 10 4 𝑞 over values of 𝑅𝑎 𝑞 .

As expected for such a reasoning, we get the same exponent but a different expression for the coefficient,

𝑎 𝑇 = 1 2 𝐴 -1 4 .
We can assume that the dynamics of the bottom boundary layer is entirely determined by its own stability, which implies that the value of 𝐴 is independent from Φ. In that case, we expect the Nusselt number to change by a factor of 2 between the Φ → 0 and Φ → ∞ limits, which is larger but similar to what can be seen on figure 8 for a free-slip BC.

Using a more sophisticated theory, [START_REF] Roberts | Fast viscous bénard convection[END_REF] obtained scaling laws for convection with both free-slip and noslip boundary conditions, both considering fixed temperature thermal boundary conditions. In both cases, he obtained a scaling relationship Nu = 𝑎Ra 𝛽 , with 𝛽 = 1∕3 for free-slip BCs and 𝛽 = 1∕5 for rigid BCs. In order to compare these scaling laws to our results, we need to account for the change of definition for the Rayleigh number:

Ra 𝑞 = Ra 𝑞 ′ 𝑑 𝑘Δ T ′ . ( 24
)
Then,

𝑞 ′ = Nu 𝑘Δ T ′ 𝑑 = 𝑎 𝑘Δ T ′ 𝑑 ( 𝛼𝑔𝜌Δ T ′ 𝑑 3 𝜅𝜂 ) 𝛽 (25)
which gives us 1∕Δ T ′ and

𝑁𝑢 = 1 Δ T = 1 Δ T ′ 𝑞 ′ 𝑑 𝑘 = 𝑎 1 𝛽+1 Ra 𝛽 𝛽+1 𝑞 . ( 26
)
Finally, for a rigid BC,

𝑁𝑢 = 𝑎 𝑇 Ra 1 6 𝑞 , ( 27
)
𝑎 𝑇 being different than previously (See figure 8.d for values).

Numerical applications to planetary objects

In this section, we show how the dimensionless results obtained in our study can be applied to various planetary objects using numerical values listed in table 2 to get dimensional results regarding the mass flux across the ice-ocean interface and the bottom temperature. Some parameters are assumed to be known, essentially because they are not expected to be too much in error with respect to the listed values. On the other hand, we express our results in a way that permits to see the effects of a change in the heat flux from the core, the thickness of the HP ice layer and the ice viscosity, over reference values, all having highly uncertain values. 𝑏 Valid for the reference thickness 𝑑 of the HP ice layer chosen for Ganymede. 𝑐 Bridgman (1912, 1937) (Valid for the water-rich exoplanet at the ice/ocean interface). 𝑑 Valid at ice/ocean interface [START_REF] Tchijov | Heat capacity of high-pressure ice polymorphs[END_REF].

Parameter

Application to Ganymede

Ganymede is the largest moon of the solar system, with a water-ice shell thickness in the 600-900 km range [START_REF] Hussmann | 10.18 -Interiors and Evolution of Icy Satellites[END_REF]. The HP ice layer thickness is estimated between 100 and 400 km (Kalousová et al., 2018).

For this numerical application we consider a HP ice layer reference thickness of 100 km, a reference heat flux from the core of 10 mW m -2 and a reference ice VI viscosity of 10 15 Pa s (Kalousová et al., 2018). With these parameters, the Rayleigh number can be computed as

Ra 𝑞 = 𝛼𝑔𝑞𝜌𝑑 4 𝑘𝜅𝜂 = 4.85 × 10 8 ( 𝑞 10 mW m -2 ) ( 𝑑 100 km ) 4 ( 𝜂 10 15 Pa s ) -1 (28)
It shows that the reference values give a large value of the Rayleigh number, of order 5 × 10 8 , and larger values can be obtained for a thicker layer, a larger heat flow from the core or a lower ice viscosity. The reference value is already quite large, in particular compared to the critical value for the onset of convection for 𝛾 = 0.9, which is in range

Ra 𝑐 =
) -0.47 , ( 29
) 1 Δ T =0.531Ra 0.2 𝑞 𝑘 𝑞𝑑 = 4.6 × 10 -2 K -1 ( 𝑞 10 mW m -2 ) -0.8 ( 𝑑 100 km ) -0.2 ( 𝜂 10 15 Pa s ) -0.2 . ( 30
)
The scaling relationship for the temperature difference across the ice layer can be used to compute the mean temperature at the ice-rock interface Tfrom the one at the ice-ocean interface 𝑇 + as: with 𝑇 + = 𝑇 + 𝑚 = 321 K (see table 2) the melting temperature of ice VI at the pressure relevant for the ice-ocean boundary (Kalousová et al., 2018;[START_REF] Bridgman | Water, in the liquid and five solid forms, under pressure[END_REF][START_REF] Bridgman | The phase diagram of water to 45,000 kg/cm 2[END_REF]. This temperature, obtained for our reference parameters (eq. 30), should be compared to the melting temperature of ice VI at the pressure relevant for the ice-rock boundary Kalousová et al., 2018;[START_REF] Bridgman | Water, in the liquid and five solid forms, under pressure[END_REF][START_REF] Bridgman | The phase diagram of water to 45,000 kg/cm 2[END_REF]. This means that taking the nominal values for all the parameters implies a large amount of melting at the bottom of the ice layer. All these parameters being quite uncertain, we can instead compute the value of any chosen parameter such that the bottom temperature is equal to the melting one, all other parameters being kept the same. The melting temperature values 𝑇 - 𝑚 and 𝑇 + 𝑚 respectively at the bottom and the top of the HP ice layer depend on the pressure at these corresponding depths and, then, for a fixed rocky core radius, on the thickness 𝑑 of the ice shell. For example, still in the case of a rigid bottom BC, for a fast phase change Φ = 10 -2 at the top boundary and taking the reference values for 𝑑 and 𝜂, it is possible to compute the minimum value 𝑞 𝑙 of the heat flux for which T -= 𝑇 - 𝑚 as follow:

T -= Δ T + 𝑇 + = 343K (31)
𝑇 - 𝑚 = 332 K (
𝑞 𝑙 = 𝑘Δ𝑇 𝑚 𝑑 𝑎 1∕(1-𝛽 𝑇 ) 𝑇 ( 𝛼𝑔𝜌𝑑 3 Δ𝑇 𝑚 𝜅𝜂 ) 𝛽 𝑇 ∕(1-𝛽 𝑇 ) = 4.3 mW m -2 (32)
with Δ𝑇 𝑚 = 𝑇 - 𝑚 -𝑇 + 𝑚 . Note that the expression in brackets is simply the Rayleigh number of the layer based on the temperature difference Δ𝑇 𝑚 . For a rigid BC, even a heat flux from the core as small as 4.3 mW m -2 could be sufficient to melt the bottom of the HP ice layer in Ganymede. This is further discussed in §4.

These calculations have also been performed for a free-slip BC and the results can be seen in table 3 for comparison.

As expected, the free-slip situation leads to a convection that is faster, more efficient to transfer heat and therefore to a cooler bottom temperature. Then, Tdoes not reach the melting temperature when using the reference values of table 2 for 𝑑, 𝑞 and 𝜂. Indeed, a minimum heat flux of 11.4 mW m -2 would be necessary in the case of a free-slip BC for the mean temperature to equal the melting temperature. However, the obtained temperature is the lateral mean one, which varies a lot along the shell (See fig. 3a), and is only one Kelvin under 𝑇 - 𝑚 , meaning that pockets of melt likely exist along the bottom boundary. In addition, the parameters chosen here for Ganymede are rather conservative and it is in fact likely that significant melting actually occurs.

Application to a water-rich exoplanet

To compare with a larger planetary object, a numerical application has been performed on an hypothetical ocean exoplanet modelled as described on figure 1, with a surface ocean instead of ice Ih. All the values considered for this application are coming from the study of [START_REF] Hernandez | Stability of high-temperature salty ice suggests electrolyte permeability in water-rich exoplanet icy mantles[END_REF] with some adjustments to match the conditions used here. We consider a rocky core/mantle of density 𝜌 𝑐 = 6 × 10 3 kg m -3 with a radius 𝑅 𝑐 = 5 000 km and a HP ice shell thickness of around 𝑑 = 260 km in order to have an aspect ratio 𝛾 = 0.95 as the scaling laws of part § 3.2 are computed for this value. The aspect ratio considered by [START_REF] Hernandez | Stability of high-temperature salty ice suggests electrolyte permeability in water-rich exoplanet icy mantles[END_REF] is ∼ 0.6 and the HP ice layer thickness considered in their study is much larger than the one used for the current numerical application. Therefore the melting temperature at the top boundary has to be adjusted and is computed as

𝑇 𝑚 = 𝑇 + 𝑚 + (𝑇 - 𝑚 -𝑇 + 𝑚 ) 𝑃 -𝑃 + 𝑃 --𝑃 + , ( 33
)
with

𝑃 = 𝑃 --𝜌𝑔𝑧. ( 34
)
The pressure considered in [START_REF] Hernandez | Stability of high-temperature salty ice suggests electrolyte permeability in water-rich exoplanet icy mantles[END_REF] at the bottom of the HP ice layer is in the range 30 -40 GPa for a melting temperature 𝑇 - 𝑚 = 800 K, which corresponds to ice VII. The melting temperature at the top of their HP Results to the numerical applications to an ocean exoplanet using parameters from table 2, for Φ = 𝟣𝟢 -𝟤 and 𝛾 = 0.95

Bottom BC Rigid Free-slip Ra 𝑞 𝟨.𝟩 × 𝟣𝟢 𝟫 𝟨.𝟩 × 𝟣𝟢 𝟫 𝑤 𝑡𝑜𝑝 (𝖼𝗆 𝗒𝗋 -𝟣 ) 𝟧𝟩 𝟨𝟣 1∕Δ T (𝖪 -𝟣 ) 𝟣 × 𝟣𝟢 -𝟤 𝟤.𝟦 × 𝟣𝟢 -𝟤 T -(𝖪)
𝟪𝟥𝟢 𝟩𝟩𝟤 ice layer is 𝑇 + 𝑚 = 367 K, which implies a pressure of about 2 GPa for ice VII. Then, we can compute the pressure 260 km above the ice/core boundary, which is the top of the HP ice layer in our case. This value should be in the range 25 -35 GPa and the melting temperature associated is around 730 K. The gravitational acceleration is computed as function of 𝑅 𝑐 and 𝜌 𝑐 as 𝑔 = 4𝜋𝐺𝜌 𝑐 𝑅 𝑐 ∕3 = 8.4 m s -2 with 𝐺 = 6.67 × 10 -11 m 3 kg -1 s -2 .

The heat flux from the core depends on the composition of the star, which will imply more or less radioactivity inside the planet. If we consider a star similar to ours, the heat flux is function of the core radius. For a given concentration in heat producing elements, the total radiogenic heating in the core scales as 𝑅 3 𝑐 while the surface scales as 𝑅 2 𝑐 and, therefore, its contribution to the heat flux density scales as 𝑅 𝑐 . As the value of 𝑅 𝑐 chosen for our application is ∼ 2.5-3

times higher than the one of Ganymede, the heat flux should be about 3 times higher. For this numerical application we consider a reference heat flux from the core of 30 mW m -2 and a reference ice VI viscosity of 10 17 Pa s. With these parameters, the Rayleigh number Ra 𝑞 , the radial velocity at ocean/ice interface 𝑤 𝑡𝑜𝑝 and the Nusselt number 1 Δ T for a rigid BC and a fast phase change Φ = 10 -2 at the top boundary have been computed as follows:

Ra 𝑞 = 𝛼𝑔𝑞𝜌𝑑 4 𝑘𝜅𝜂 = 6.7 × 10 9 ( 𝑅 𝑐 5 000 km ) ( 𝑞 30 mW m -2 ) ( 𝑑 260 km ) 4 ( 𝜂 10 17 Pa s ) -1 , ( 35
)
𝑤 𝑡𝑜𝑝 = 𝜅 𝑑 2.62 ×
) -0.47 , ( 36
)
1 Δ T = 𝑘 𝑞𝑑 5.31 × 10 -1 Ra 0.2 𝑞 = 1 × 10 -2 K -1 ( 𝑅 𝑐 5 000 km ) 0.2 ( 𝑞 30 mW m -2 ) -0.8 ( 𝑑 260 km ) -0.2 ( 𝜂 10 17 Pa s ) -0.2 . ( 37
)
The results can be found in table 4 for comparison between rigid and free-slip boundary conditions and with

Ganymede numerical application given in table 3. Obviously, the reference Rayleigh number is higher in the case of a large ocean exoplanet than for a Ganymede-like body, meaning a more efficient convection leading to a higher Nusselt number and a higher dimensionless Δ T , but a larger dimensional temperature gap between the two HP ice layer boundaries. But, the layer being ∼ 2.5 times thicker, the melting temperature gap between boundaries is also important and this should not have a huge impact on the melt production at the core/ice interface. For the same efficiency of the phase change at the top, the melting ability at the bottom of the HP ice layer, looking at the temperatures Tobtained, are comparable to the Ganymede case, but the difference of bottom temperature between rigid and free-slip BC seems to be more pronounced in the case of a large water-rich exoplanet. In fact, for both objects, the temperature difference across the ice layer is roughly twice larger for a rigid boundary condition than for a free-slip one. Scaling back to dimensional units for a larger object leads however to larger dimensional differences. In the case computed here, Tdoes not reach 𝑇 - 𝑚 for a free-slip BC but it does for a rigid boundary condition. First of all, as already discussed above for the application to Ganymede, the exact parameters for this application are subject to discussion and higher temperatures are possible. But the fact that a free-slip boundary condition leads to a lower value of the bottom temperature than a rigid one remains and this could lead to a conundrum: the free-slip BC that is justified by the presence of melt could predict its absence while the rigid BC that results from the absence of melt could predict its presence. In this case, we would expect a mixed situation with the presence of melt where the temperature is high and its absence elsewhere with a non-uniform boundary condition that depends on the local temperature. Such a situation could be the topic of future studies.

Finally, the vertical velocity at the top of the HP ice layer, 𝑤 𝑡𝑜𝑝 , is larger for this study case, meaning that the mass transfer between the core and the ocean should be more important.

Discussion

The model and results presented above are simple and concentrate on one effect not considered before, the phasechange boundary condition at the upper boundary. We showed that this effect alone is sufficient to considerably alter the solution, its structure and the efficiency of heat and mass transfer. Interestingly, the temperature profiles we obtain are similar to the ones Kalousová et al. (2018) and Kalousová and Sotin (2018) obtained with a mechanical conditions that is non-penetrative for the solid but considering the possibility of partial melting in the bulk of the ice layer and extraction of the liquid.

As mentioned before, in this study we considered the HP ice layer as an isoviscous material, for the sake of simplicity, while the viscosity is expected to depend on pressure, temperature and melt fraction [START_REF] Choblet | Heat transport in the high-pressure ice mantle of large icy moons[END_REF]Kalousová and Sotin, 2018;Kalousová et al., 2018). To get a hint on the potential effect of viscosity variations, we performed a single simulation with a depth-(pressure-) and temperature-dependent viscosity following the dimensionless Arrhenius law established so that the surface viscosity serves as reference in the Rayleigh number:

𝜂 = exp ( 𝐸 + (1 -𝑧)𝑉 𝑇 + 1 -𝐸 ) , (38)
with 𝐸 and 𝑉 the dimensionless activation energy volume. The values 𝐸 = 230 and 𝑉 = 4.6 are chosen for the simulation to be comparable with the ones of [START_REF] Choblet | Heat transport in the high-pressure ice mantle of large icy moons[END_REF]. This simulation has been ran for a rigid bottom boundary condition, Φ = 10 -2 and Ra 𝑞 = 10 6 . As can be seen on fig 10, the lowest viscosity is located at the bottom boundary and inside the up-welling hot plumes. This is consistent with the fact that the temperature is higher in these locations and could reach the melting temperature, implying the formation of melt and a lower viscosity. Beside that obvious difference in viscosity, the thermal structure of the solution is similar to that obtained with a constant the radial velocity at the top is 𝑅𝑎 𝑤 = 2.5 × 10 6 and the one corresponding to the Nusselt number is 𝑅𝑎 𝑁𝑢 = 8.1 × 10 6 , compared to the nominal value of 10 6 used for this calculation. It means that the calculation with this type of viscosity variation leads to results similar to those obtained with a constant viscosity but a slightly larger Rayleigh number. This is not surprising since the average viscosity is lower than one. It would however be worth exploring these effects more systematically in the future, especially if we take into account partial melting.

As already mentioned ( §3.2), for some values of the input parameters, we predict a temperature at the bottom of the ice layer that is larger than the melting temperature. This should result in the presence of melt, either everywhere or restricted to melt pockets, since the temperature is laterally variable. This melt, containing salts by interacting with the core, could be transported to the ocean by various processes, as modeled by Kalousová et al. (2018). This partial melting process has not been considered in this preliminary study, which only takes into account solid-state convection through the HP ice layer, but could have an important impact on the efficiency of the mass transfer. To justify the need of adding this property in future models, we computed the regions where melting would happen and the amount of melt produced in considering Ganymede's HP ice layer if the temperature field we obtained was not too affected by the presence of melt. This is obviously a very crude assumption but it can be useful to provide a first order estimate. For this calculation, we use the same reference values than previously used in § 3.3.1 and other parameters from table 2.

The variations of the melting temperature 𝑇 𝑚 of ice VI as function of pressure 𝑃 has been computed from eq. 33, with 𝑃 = 𝑃 + + 𝜌𝑔𝑧, using values given in table 2.

Comparing the actual temperature at each point to the melting value at the same position allows to find places that would be expected to partially melt. For regions where the temperature 𝑇 > 𝑇 𝑚 , the melting temperature at this depth, partial melting should occur and this should limit the temperature to the melting one, at least as long as the melt fraction is lower than 1. We compute the expected melt fraction by equating the energy associated with the temperature in excess to the melting temperature, Δ𝑇 = 𝑇 -𝑇 𝑚 , to latent heat of melting. The melt fraction is then

𝑓 𝑀 = 𝐶 𝑝 Δ𝑇 𝐿 ( 39
)
with 𝐶 𝑝 the heat capacity and 𝐿 the latent heat as defined in table 2. Since this calculation does not account for porous flow, the melt is assumed to stay in place, even for high fractions. In reality, differential motion of liquid and solid resulting from gravity would act to limit this melt fraction within the ice layer. The results presented in this part are therefore only indicative of the relative importance of melt as function of the input parameters.

Radius

Ra q = 1.39e + 05 Ra q = 2.78e + 05 Ra q = 1.39e + 06 Ra q = 2.78e + 06 Ra q = 2.78e + 07

Figure 12: This plot shows the horizontal average melt fraction 𝑓 𝑀 produced in the HP ice layer as function of the radius for a rigid bottom boundary with Φ = 10 -2 , 𝛾 = 0.95, for various values of Ra 𝑞 in the case of Ganymede, using reference parameters values given in table 2 for d & q, implying a viscosity in the range 𝟣𝟢 𝟣𝟨 -𝟣𝟢 𝟣𝟪 𝖯𝖺 𝗌 which corresponds more or less to the ranges defined for Ganymede. Note that these figures are an a posteriori prediction from calculations that do not consider partial melting in a self-consistent way (see text for details).

The melt fraction decreases with the Rayleigh number and reaches zero in the bulk of the layer for Ra 𝑞 ≳ 10 7 . In that case, the melt would refreeze during its ascent before remelting upon reaching the ice/ocean interface. As a mass transfer mechanism from the core to the ocean, this situation is likely less efficient than the cases with lower Rayleigh number which allow partial melt to be maintained in the whole layer. The amount of melt produced for low Ra 𝑞 cases in fig 12 reaches very high values which clearly falls out of the range where the assumption of no motion between phases fails. However, the Rayleigh number expected for Ganymede is at least Ra 𝑞 ∼ 5 × 10 8 as computed in § 3.3.1, for which the maximum melt mass fraction produced at the bottom of the HP ice layer would be smaller than for the ones shown in fig 12 .     Figure 13 shows spherical representations of 1∕5 𝑡ℎ of the HP ice shell, with a zoom on several plumes. The left panels depict the temperature field for various values of the heat flux coming from the core, the yellow line being the contour on which T = 𝑇 𝑚 . The expected melt fraction is shown on the right panels of the figure. Each panel is in fact based on the same snapshot of the results obtained for a rigid BC at the bottom of the HP ice layer when Ra 𝑞 = 10 3 × Ra 𝑐 ∼ 2.8 × 10 5 and Φ = 10 -2 and differ only by the choice of the heat flux to render the results dimensional. It shows that three options are conceivable depending on the heat flux.

In the first case, obtained for small values of the heat flux (𝑞 = 1 mW m -2 on fig 13), partial melting only occurs in hot plumes close to the upper boundary where it can reach a fraction of at most 0.5%. Actually taking into account partial melting in the model would make this melt rise faster to the ocean but would probably not change the overall dynamics. Salts can enter the ice layer by contact with the rocky core but is limited by their partitioning behaviour and by diffusion in both solids, a rather inefficient process. We expect therefore that the salt concentration is rather small in the ice layer and its effect on the dynamics should be limited [START_REF] Hernandez | Stability of high-temperature salty ice suggests electrolyte permeability in water-rich exoplanet icy mantles[END_REF]. However small in quantity, all the salts introduced at the bottom should be efficiently transported by convection in the solid.

For an intermediate value of the heat flux 𝑞 = 5 mW m -2 on fig. 13, we obtain localised melt pockets in the vicinity of hot plumes, both at their roots and close to the upper boundary, while these are not generally connected with each other. It means the liquid water would refreeze upon ascending in the ice layer and remelt again before reaching the overlying ocean. We expect that the liquid in contact with the rocky core would get enriched in salts, which would help keep it liquid while ascending and potentially making a continuous path toward the ocean, possibly as an independent porous flow (Kalousová and Sotin, 2018;Kalousová et al., 2018;[START_REF] Choblet | Heat transport in the high-pressure ice mantle of large icy moons[END_REF]. If it crystallises on the way up, the ice-salt mixture is likely to be denser than pure ice, which could slow down the up-welling flow, depending on the relative effects of salt and temperature on the density. A full treatment of that case requires a model that includes both a two-phase-flow approach (Kalousová and Sotin, 2018;Kalousová et al., 2018) and a proper treatment of composition variations. This will be the topic of our future studies.

Finally, for a large heat flux 𝑞 = 40 mW m -2 on fig 13, we get melt everywhere at the bottom that even reaches complete melting at places (melt fraction equal to 1). If the melt is so rich in salts that it is denser than pure ice, we expect a full layer of liquid below the ice layer, which would require applying a phase change boundary condition there as well as at the surface. Previous studies [START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF][START_REF] Morison | Timescale of overturn in a magma ocean cumulate[END_REF][START_REF] Morison | Convection in the primitive mantle in interaction with global magma oceans[END_REF] show that we should expect a totally different flow pattern dominated by a degree one mode of convection, a translation of the shell. This is a very efficient heat transport mechanism that would act to cool down the core very efficiently and freezing the dense water. Clearly, a more complete study of that scenario is needed to conclude. On the other hand, in the likely situation where the liquid water contains a small enough fraction of salt to stay less dense than the ice above, it is likely to transport it very efficiently toward the ocean as a porous flow. Even with our purely solid calculations, we predict a continuous connection between the melt layer at the bottom and the upper boundary, a situation that would be reinforced by the presence of salts in the water. We therefore expect this situation to be the most efficient one to transport salts from the core to the ocean.

Several choices have been made for the numerical applications and need to be discuss (See § 3.3). First, despite the range estimated for the value of Φ for Ganymede in § 2.1.2, we chose Φ = 10 -2 for the numerical application (See § 3.3.1). This range is a very rough estimate and according to figure 9, between Φ = 10 -2 and Φ ∼ 10, the coefficient 𝑎 𝑤 of the scaling law for the radial velocity hardly varies. Therefore, we chose the most extreme case among those we have studied for perform these calculations. Also, the scaling laws in § 3.2 are established for an aspect ratio in the range 0.9-0.95 (See table 1). For Ganymede, this ratio is not yet well constrained, with a value varying from 0.8 to 0.95 depending on the core radius, the thickness of the HP ice layer and the period considered in the different studies.

From the study of [START_REF] Bland | The orbital-thermal evolution and global expansion of Ganymede[END_REF], in the thermal history of Ganymede, the thickness of the HP ice layer could have vary between 100 km for a heat flux of about 25 mW m -2 and 500 km (∼700 km if the ocean is closed) for a heat flux of about 5 mW m -2 . The aim of this first numerical application was not to conclude about the actual efficiency of heat and mass transfers from the core to the ocean on Ganymede, but to give an idea of what kind of exchanges could happen through an HP ice layer for that type of planetary objects. Then, we chose to base this numerical application on the study of Kalousová et al. (2018), which considers a large range of combinations between the heat flux from the core and the thickness of the HP ice layer and we decided to represent the extreme case with the smallest values of 𝑞 and 𝑑, which means a minimum of melt at the bottom, in order to observe if melt could occur at the interface between the core and the HP ice shell even in this specific case. In addition, the main parameters currently not well constrained (𝑞, 𝑑 and 𝜂) are kept in clear in the application equations (See eq. 28, eq. 29 and eq. 30) and can be quickly modified to fit another planetary object, as Titan. Finally, future missions should one day allow to better constrain these parameters for various bodies and this study will allow to quickly conclude about the presence of melt at the bottom and the mass transfer efficiency through the HP ice layer.

Conclusions

This paper addresses the possibility of convection in a high pressure ice layer between a solid core and a liquid water ocean on icy/water-rich bodies and its heat and mass transfer efficiency. We include an effect not considered previously in this context: the solid-liquid phase change at the upper boundary and its implications for the mechanical boundary condition. As shown by previous papers in other contexts [START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF][START_REF] Deguen | Thermal convection in Earth's inner core with phase change at its boundary[END_REF][START_REF] Labrosse | Rayleigh-Bénard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries[END_REF][START_REF] Morison | Timescale of overturn in a magma ocean cumulate[END_REF][START_REF] Agrusta | Mantle convection interacting with magma oceans[END_REF], this leads to a non-zero vertical velocity at the ocean/ice interface, a decreased value for the critical Rayleigh number, a markedly different thermal structure and an enhanced heat and mass transfer efficiency. In order to conclude about the effect of this only aspect of the problem, convection of constant-viscosity pure solid water ice has been considered. Compared to previous papers that consider the effect of the phase change boundary condition, we consider here several aspects that are specific to the application of the HP ice layers in contact with a rocky core: we consider a thin spherical shell, with boundary conditions at the bottom not considered before, an imposed heat flux and either a rigid or free-slip mechanical condition. On the other hand, for the sake of simplicity, we have not included the effect of variable viscosity and partial melting that were considered by some previous studies [START_REF] Choblet | Heat transport in the high-pressure ice mantle of large icy moons[END_REF]Kalousová et al., 2018;Kalousová and Sotin, 2018). Future studies should be performed to combine all these effects. First of all, our results largely confirm those obtained on the effect of the phase change boundary condition in other contexts. The critical Rayleigh number decreases with a decrease of the phase change number, while the wavelength of the most unstable mode increases [START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF][START_REF] Labrosse | Rayleigh-Bénard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries[END_REF][START_REF] Morison | Convection in the primitive mantle in interaction with global magma oceans[END_REF]. The Rayleigh number expected for the HP ice layer of most icy satellites and planets is expected to be largely supercritical and, in that regime, the radial velocity at the upper boundary and heat transfer efficiency, as measured by the Nusselt number, scale as power laws of the Rayleigh number [START_REF] Agrusta | Mantle convection interacting with magma oceans[END_REF]. With the setup considered here, the exponents are close to 1∕2 and 1∕5 for the velocity and the Nusselt number, respectively, for a rigid BC at the bottom and 1∕2 and 1∕4 for a free-slip BC. The coefficients of the scaling laws increase when the phase change number is decreased, which eases the phase change, by roughly a factor 2 between end-members for the Nusselt number.

The scaling laws obtained in a dimensionless parameter space can be easily applied to any object by chosing the relevant values of all parameters. Section 3.3 presents applications to Ganymede and a hypothetical large water exoplanet, imposing values for the best constrained parameters and leaving the possibility of adjusting around reference values for the others, in particular the heat flux from the core, the ice layer thickness and its viscosity. In both cases, we find a typical RMS velocity across the upper boundary of the order of 50 cm yr -1 in the limit of a small phase change number. With such values, the mass exchange between the ice and the overlying ocean should be quite efficient and any salt added to the ice layer by interacting with the underlying core should be easily transported to the ocean.

Conversely, if the ocean starts already salted, the efficient mass transfer at the top of the HP ice layer should lead to a rapid chemical equilibrium [START_REF] Bolrão | Timescales of chemical equilibrium between the convecting solid mantle and over-and underlying magma oceans[END_REF].

Using the heat transfer scaling laws, we can predict the temperature at the bottom of the HP ice layer as function of the governing parameters. Using the reference values, we predict that the mean temperature at the bottom of the ice layer is close to the melting temperature, both for Ganymede and a large water planet. All parameters being equal, the temperature is lower for a free-slip boundary condition at the bottom than for a no-slip one but in both cases, it is close to the melting one. Considering the lateral variations of the temperature, it means that we should expect the presence of partial melt at least in some areas of the core-ice interface, and the mechanical boundary condition should be intermediate between the two end-members considered in this study. Chemical interaction with the rocky core by hydrothermal activity should enrich this melt in various salts, and its upward motion by porous flow should make it freeze and bring salts in the ice layer. The effect of this salt on convection remains to be studied in details.

Further studies focusing on salts and partial melting will complete the model in order to be more relevant to conclude about mass transfer efficiency through HP ice layer on icy/water-rich bodies. As several parameters are not well constraint yet, these studies will be useful when more information from future exploration missions will be available in order to conclude about the habitability of these bodies from this perspective.

Acknowledgments

This study benefited from discussion with Gaël Choblet and Thierry Alboussière and we gratefully acknowledge constructive comments from Gabriel Tobie and an anonymous reviewer. The support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources is also acknowledged.

  Figure1: Model illustration for the interior of an ocean world with an HP ice layer in direct contact with a liquid ocean at top and a rocky core/mantle at the bottom (type IV or V of[START_REF] Lammer | Origin and Evolution of Planetary Atmospheres -Implications for Habitability[END_REF] type H3 of Noack et al., 2016). Our numerical models treat convection in the HP ice layer.

Figure 2 :

 2 Figure 2: Snapshots of 1∕5 of the HP ice shell for a given Ra 𝑞 ∕Ra 𝑐 = 10 5 for Φ → ∞ (a) and Φ = 𝟣𝟢 -𝟤 (b). (c) is a zoom of (b) with the top radial velocity displayed in green, varying from ∼ -1300 to ∼ 3350 along the top boundary. The temperature scale on the right is common to all three panels.

Figure 3 :

 3 Figure 3: Figure of a) the minimum (dashed lines), maximum (dashdot lines) and horizontal average (solid lines) temperatures and b) the radial velocity through the HP ice layer for a given Ra 𝑞 ∕Ra 𝑐 = 10 5 for Φ = 10 -2 (blue lines) and Φ → ∞ (orange lines). The grey line in a) is for a null temperature.

Figure 4 :

 4 Figure 4: First unstable convective modes for Ra 𝑞 = Ra 𝑐 as function the bottom boundary condition (Free-slip or rigid) and the value of the phase change parameter Φ. The color gives the temperature anomaly with respect to the steady conduction profile (arbitrary units, amplitude not constrained by the linear stability analysis) and the contours give the stream function. The input parameters and output critical Ra and wavenumber 𝑙 𝑐 are provided inside each shell.

  Figure 5: The three panels are showing, from left to right respectively, the critical Rayleigh number, critical wave-number and the minimum thickness 𝑑 of the HP ice layer for the onset of convection as function of Φ, depending on the bottom boundary condition (BC) and the shell aspect ratio 𝛾.

Figure 6 :

 6 Figure 6: Snapshot of the perturbation temperature along the entire HP ice layer for Ra 𝑞 = 1.1Ra 𝑐 (left) and Ra 𝑞 = 2Ra 𝑐 (right) in the case of a free-slip BC for 𝛾 = 0.9 and Φ = 10 -2 .

Figure 7 :

 7 Figure 7: Cases Ra 𝑞 ∼ Ra 𝑐 for a rigid bottom BC (a)&(b) and a free-slip bottom BC (c)&(d). Left panels (a)&(c) show the dimensionless radial velocity at the top boundary of the HP ice layer, 𝑤 𝑡𝑜𝑝 . Right panels (b)&(d) show the dimensionless heat flux (Nusselt number). Both parameters are shown as functions of the Rayleigh Number Ra 𝑞 and for various values of Φ.

Figure 8 :

 8 Figure 8: Cases Ra 𝑞 ≥ 50Ra 𝑐 for a rigid bottom BC (a)&(b) and a free-slip bottom BC (c)&(d). Left panels (a)&(c) show the dimensionless radial velocity at the top boundary of the HP ice layer. Right panels (b)&(d) show the dimensionless heat flux (Nusselt number). Both parameters in function of the Rayleigh Number and for various values of Φ.

Figure 9 :

 9 Figure 9: Variation of the coefficient 𝑎 𝑤 of the scaling law for the radial velocity at the top boundary as function of the phase change number Φ (See eq. 18). For each value of Φ, it is computed as the mean of 𝑤 𝑡𝑜𝑝 ∕Ra 1 2

  65 -915, depending on the bottom BC and the value of Φ (see fig 4). All other parameters being kept identical, the critical value Ra 𝑐 = 292 for a rigid bottom BC when Φ = 10 -2 is reached for a thickness of the ice layer around 1.6 km (see fig 5).Considering first the situation for which the bottom boundary of the HP ice shell of Ganymede is a no-slip one (i.e. in the absence of liquid water to lubricate that boundary) and the phase change at the top boundary is fast, Φ = 10 -2 , with an aspect ratio 𝛾 = 0.95, from the scaling laws obtained in § 3.2 on figure 8:

  viscosity, with a boundary layer at the bottom from which hot plumes originate and a passive downward return flow whose temperature is that of the upper boundary. Comparing fig 11.a) and fig 11.b), we can see that our chosen viscosity law (eq.38) makes it increase by 10 2 with depth and vary by more than six orders of magnitude with temperature. Also, we can observe on fig 11.c) that the RMS radial velocity curve shape is slightly different than the one for an isoviscous calculation (as can be seen on the blue curve of fig 3.b). In both cases, the radial velocity increases with height in the layer, but the curvature are somewhat different. Of particular interest here are the local maximum reached in the bottom boundary layer, where the viscosity is lowest and the increase when getting close to the upper boundary, because of the decrease of the viscosity with height. When the statistically steady-state is reached, the time-averaged Nusselt number is Nu = 12.8 and the time-averaged radial velocity at the top is 𝑤 𝑡𝑜𝑝 = 264.9. From the scaling obtained in § 3.2 for a constant viscosity and the same boundary conditions (See fig 8.a and fig 8.b), the Rayleigh number corresponding to

Figure 10 :

 10 Figure 10: Snapshots of 1∕14 of the HP ice shell for Ra 𝑞 = 𝟣𝟢 𝟨 , Φ = 𝟣𝟢 -𝟤 and a rigid bottom boundary condition. The left (respectively right) panel shows the viscosity variations (respectively temperature) through the HP ice layer.

Figure 12

 12 Figure12shows the horizontal average of the melt mass fraction produced into the HP ice layer as function of the radius for various values of Ra 𝑞 , considering a rigid BC at the bottom of the ice shell, a phase-change number Φ = 10 -2 and the reference values of 𝑞 and 𝑑 identified in table 2, implying a viscosity in the range 10 16 -10 18 Pa s which corresponds more or less to the ranges defined for Ganymede. For all values of the Rayleigh number considered, the melt fraction profiles are similar to the mean temperature profile as shown on figure3. The melt fraction is largest at the bottom boundary and decays with radius in the boundary layer until it reaches a nearly constant value in the bulk. The cases with the lowest values of the Rayleigh number show a slight increase when getting close to the upper boundary.

Figure 11 :

 11 Figure 11: Radial profiles of the minimum (orange), maximum (green) and horizontal average (blue) viscosity (a) and temperatures (b). c) shows the RMS radial velocity profile through the HP ice layer. All the panels come from a varying viscosity simulation made with Ra 𝑞 = 𝟣𝟢 𝟨 , Φ = 𝟣𝟢 -𝟤 and a rigid bottom boundary condition.

Table 1

 1 Detailed parameters applied to all simulations for each study-case. When several numbers are written in cells, the model has been run for all possible combinations.

	Study case	Ra 𝑞 ∕Ra 𝑐	Bottom boundary Phase change at top	Φ	𝛾
	SC1 SC2	1.1, 2, 3, 4, 5, 6, 7, 8, 9, 10	Freeslip Rigid	True False True False	10 -2 , 10 -1 , 1, 10, 10 2 ∞ 10 -2 , 10 -1 , 1, 10, 10 2 ∞	0.9,0.95 0.9,0.95
	SC3 SC4	5.10 1 , 1.10 2 , 5.10 2 , 1.10 3 , 5.10 3 , 1.10 4 , 1.10 5	Freeslip Rigid	True False True False	10 -2 , 10 -1 , 1, 10, 10 2 ∞ 10 -2 , 10 -1 , 1, 10, 10 2 ∞	0.9,0.95 0.9,0.95

Table 2

 2 Detailed numerical application parameters for Ganymede fromKalousová et al. (2018) and a water-rich exoplanet from[START_REF] Hernandez | Stability of high-temperature salty ice suggests electrolyte permeability in water-rich exoplanet icy mantles[END_REF].

*

Are the reference values used in the numerical application in section 3.3.1

𝑎 Hussmann et al. (2015)

.

Table 3

 3 Results to the numerical applications to Ganymede using parameters from table 2, for Φ = 𝟣𝟢 -𝟤 and 𝛾 = 0.95 𝟣𝟢 𝟪 𝟦.𝟪𝟧 × 𝟣𝟢 𝟪 𝑤 𝑡𝑜𝑝 (𝖼𝗆 𝗒𝗋 -𝟣 ) 𝟦𝟤.𝟫 𝟦𝟩.𝟨 1∕Δ T (𝖪 -𝟣 ) 𝟦.𝟨 × 𝟣𝟢 -𝟤 𝟢.𝟣 T -(𝖪) 𝟥𝟦𝟥 𝟥𝟥𝟣

	Bottom BC	Rigid	Free-slip
	Ra 𝑞	𝟦.𝟪𝟧 ×

Table 4

 4

Convection in HP ice layer for large icy moons

). The yellow line on left panels shows the points where T = 𝑇 𝑚 , which delimits regions of partial melting. Each panel shows the whole computational domain (1∕5 𝑡ℎ of the shell) as well as a zoom on a small fraction of the domain. Note that these figures are an a posteriori prediction from calculations that do not consider partial melting in a self-consistent way (see text for details).

A. Supplementary material

The supplementary material containing the fully detailed results tables and figures corresponding to the exploration of the parameters space for all cases described in table 1 can be found at:


domain: phys sde

text:
Michel Yves Louge

Alexandre Valance

Nikuradse and Shields: turbulent roughness controls shear stress at incipient sediment transport

Keywords:

Plain Language Summary

The transport of sediments in rivers, oceans, planets and deserts begins when fluid imparts enough stress on the sediment surface, i.e. a sufficient force in a unit of area.

An inspection of historical data reveals that such incipient stress depends exclusively on the 'turbulent roughness' of the surface through a simple expression. This recognition implies that turbulent roughness, a quantity derived from measurements of velocity profiles above the surface, uniquely determines how much fluid friction must be exerted on a sediment bed to mobilize it.

For cohesionless, frictionless sediment buffeted by a turbulent boundary layer, [START_REF] Shields | Application of similarity principles and turbulence research to bed-load movement[END_REF] showed that particle detachment from a plane bed arises when the shear stress τ exerted by a flow of shear velocity u * ≡ τ /ρ over the erodible surface exceeds particle weight lightened by buoyancy. For material densities of solid and fluid ρ s and ρ, gravitational acceleration g and particle diameter d, his static balance of forces then produced a universal diagram predicting the minimum shear stress τ c for entrainment relative to the granular normal stress σ = (ρ s -ρ) g d that the surface layer of particles exerts on the underlying bed. [START_REF] Bagnold | The physics of blown sand and desert dunes[END_REF] later noted that, counterintuitively, the dimensionless threshold stress is smaller in air than water, and [START_REF] Andreotti | A lower-than-expected saltation threshold at Martian pressure and below[END_REF] showed that it is even lower for gases under reduced pressures.

As the upper limit of the ratio of tangential traction τ exerted on a grain held on the surface and the perpendicular buoyant weight σ in a unit of sediment area, Shields' dimensionless threshold µ ≡ τ c /σ has the structure of a Coulomb (1773) friction coefficient, i.e. the bed remains static for τ < σ × µ, and it mobilizes as τ exceeds this threshold.

Because Shields' work involved particles on a fluid-solid boundary, he inevitably referenced the first data set for turbulent flow on bumpy surfaces, which [START_REF] Nikuradse | Strömungsgesetze in rauhen Rohren (Laws of flow in rough pipes[END_REF] had staged in impermeable circular pipes of radius R p roughened by gluing sand grains to their interior wall. In these internal flows, Nikuradse reported the Darcy (1856) fluid friction factor f ≡ 2R p ∇p /(ρū 2 /2) that makes the pressure gradient ∇p along the pipe dimensionless with the mean kinetic energy density (ρ/2)ū 2 .

In this Letter, we re-examine these historical data to show how Shields friction µ is a simple function of the turbulent roughness z 0 made dimensionless with its invariant magnitude z 0∞ in the limit where shear velocity is large. This observation implies that z 0 /z 0∞ governs the Coulomb-like fluid friction on the sediment bed. We also confirm that the function can be corrected to capture a wide range of solid-to-gas density ratios [START_REF] Pähtz | Scaling laws for planetary sediment transport from DEM-RANS numerical simulations[END_REF].

First, we recall the structure of the turbulent boundary layer above a flat plate, which consists of a viscous sublayer, a buffer, and a turbulent 'core' with time-averaged velocity u satisfying u/u * a u + (1/κ) ln (zu * /ν) .

(

) 1
In this expression, κ 0.41 is von Kármán's constant, ν is the fluid kinematic viscosity, z is elevation above the surface, and a u is a parameter that depends on the nature of the solid wall. When basic anemometry musters insufficient spatial resolution to resolve the viscous and buffer layers, experiments focus instead on the core, writing its profile as

u/u * = (1/κ) ln (z/z 0 ) , (2)
thereby identifying the shear velocity u * and turbulent roughness z 0 as respective scales of speed and altitude. Eliminating u/u * between Eqs. ( 1) and (2) yields the general relation ln (u * z 0 /ν) = -a u κ .

(3)

In smooth pipes, the data of [START_REF] Saph | An experimental study of the resistances to the flow of water in pipes[END_REF] indicates that a u = a u0 is constant. In that case, because a smooth wall contains no discernible asperity beside microscopic features responsible for the vanishing slip velocity on its surface, calling z 0 a 'roughness' is a misnomer. Rather than being associated with geometry, z 0 on a smooth wall arises from the interplay of coherent structures generated by the turbulent boundary layer upstream [START_REF] Wu | Direct numerical simulation of turbulence in a nominally zero-pressure-gradient flat-plate boundary layer[END_REF]).

An instructive calculation hints at the fundamental meaning of z 0 in turbulent, fullydeveloped, steady pipe flow. Neglecting small contributions of the viscous and buffer layers to the overall flow rate, an exact integration of Eq. ( 1) across a pipe of radius R p yields 8) and ( 9). For Nikuradse's fully-turbulent flows in rough pipes (ln Re * > -4), lines in (B) are simpler, albeit equally acceptable models than those found in existing literature for cohesionless, frictionless sediment. For ln Re * < -4 (dotted line), Eq. ( 11) falls outside the range of Nikuradse's experiments. There, significant deviations can arise as the concept of z0 becomes moot, or as a gas is too rarefied to behave as a continuum.

the mean velocity

ū u * = 1 κ ln u * R p ν - 3 2κ + a u . (4)
Recalling the definition of the Fanning (1896) friction factor C f ≡ τ /(ρū 2 /2), which makes wall shear stress τ dimensionless with mean kinetic energy density, a force balance on a slice of pipe in steady, fully-developed, isothermal flow implies f = 4 C f . Substituting these relations in Eq. ( 3)-( 4) yields the general equation

ln z 0 R p = - 3 2 + κ 2 C f , (5)
in which u * or the pipe Reynolds number no longer appear [START_REF] Nikuradse | Strömungsgesetze in rauhen Rohren (Laws of flow in rough pipes[END_REF]. As a general expression for pipes, Eq. ( 5) implies that, given the outer scale R p of the system, ln(z 0 /R p ) is uniquely related to, -and therefore is a measure of -the Fanning friction factor.

If the system possesses no such outer dimension, a compelling question is whether turbulent roughness may yet be appropriately scaled to remain uniquely related to a friction coefficient. If so, does this coefficient arise from fluid mechanics, as Fanning's C f and Darcy's f , or is it akin to solid friction, like Shields' µ?

While u * rises with mean pipe Reynolds number, [START_REF] Nikuradse | Strömungsgesetze in rauhen Rohren (Laws of flow in rough pipes[END_REF] showed that ln z 0 does not decrease indefinitely in a rough pipe, as Eq. ( 3) would suggest, but that it reaches instead an asymptotic value set by the diameter of sand grains that he glued to the pipe's inside surface. Here, Nikuradse's rough walls produced values of z 0 directly associated with geometrical features of the boundary. By least-squares fitting the nearly invariant Darcy friction factor at the highest Reynolds number of his experiments, one obtains the classical result for a rough wall in the limit of large u * ,

z 0 = z 0∞ ≡ d/ , (6)
where = 34.6 ± 0.3 at 95% confidence. (Data in the Supporting Information).

Whereas the radius R p of Nikuradse's pipes is a conspicuous outer scale for his flows, one cannot identify a representative distance for fully-developed turbulent boundary layers of indefinite extent above a smooth plane. However, on a sediment surface, Eq. ( 6) implies that d/ has become the natural inner scale that makes elevation dimensionless. Similarly, the natural Reynolds number in an unbounded flow builds upon d/ and the only available velocity scale u * ,

Re * ≡ u * d/ ν . (7)
Staging several combinations of pipe radius and sand grain size, [START_REF] Nikuradse | Strömungsgesetze in rauhen Rohren (Laws of flow in rough pipes[END_REF] published extensive data on f vs pipe Reynolds number Re ≡ 2R p ū/ν. Remarkably, upon using Eq. ( 5) to convert f and Re to z † 0 vs Re * , his entire set collapses onto a single curve that is independent of pipe diameter (Fig. 1A) and is closely fitted by

z † 0 ≡ z 0 z 0∞ = exp (-a u0 κ) Re * + 1 -exp (-a z Re * ) . (8)
Because Eq. ( 8) has no outer scale, it has been universally applied to internal flows as well as external boundary layers with constant properties, unless geometrical roughness varies or the surface undulates on wavelengths d. [START_REF] Nikuradse | Strömungsgesetze in rauhen Rohren (Laws of flow in rough pipes[END_REF] hinted at the existence of this relation in his Fig. 11, while [START_REF] Jia | Hydrodynamic roughness induced by a multiscale topography[END_REF] recently derived its form after identifying d as the essential wavelength of microscopic flow over the boundary.

Despite its accuracy and simplicity, fluid mechanics textbooks do not mention Eq. ( 8).

Instead, they present the equivalent, albeit more complicated [START_REF] Moody | Friction factors for pipe flow[END_REF] chart, which offers practical correlations of f vs pipe Reynolds number and roughness d/(2R p ). Fitting Eq. ( 8) to the entire [START_REF] Nikuradse | Strömungsgesetze in rauhen Rohren (Laws of flow in rough pipes[END_REF] data set yields a u0 = 5.84±0.04 and a z = 1.35±0.03 at 95% confidence. The magnitude of a u0 , which extrapolates data in rough pipes to the small Re * 'smooth' limit, is consistent with the accepted value ∼ 5 that is derived from the seminal experiments of [START_REF] Saph | An experimental study of the resistances to the flow of water in pipes[END_REF] with rigorously smooth pipes [START_REF] Blasius | Das Aehnlichkeitsgesetz bei Reibungsvorgängen in Flüssigkeiten[END_REF][START_REF] Steen | Saph and Schoder and the friction law of Blasius[END_REF]. Figure 1A shows that straight segments also delineate three distinct regimes of [START_REF] Nikuradse | Strömungsgesetze in rauhen Rohren (Laws of flow in rough pipes[END_REF] data in a plot of ln z † 0 vs ln Re * ,

ln z † 0        -a u0 κ -ln Re * low Re * limit, ln Re * < (au 2 -κ au 0 ) (1+au 1 )
-a u2 + a u1 × ln Re * intermediate,

(au 2 -κ au 0 ) (1+au 1 )
< ln Re * < au 2 au 1 0 rough limit,

au 2 au 1 < ln Re * , (9)
where a u1 0.234 ± 0.015 and a u2 0.184 ± 0.010 at 95% confidence.

Meanwhile, [START_REF] Shields | Application of similarity principles and turbulence research to bed-load movement[END_REF] inferred the threshold shear stress by writing an overall force balance on the flow of water in an inclined flume. Unlike Nikuradse's measurements, the presence of side walls and erodible sediment complicated his analysis, and therefore produced greater scatter. However, as Fig. 1 illustrates, Shields' correlation for dimensionless basal shear stress µ at incipient sediment transport bears a striking resemblance to Eq. ( 8), including the transition Reynolds numbers Re * demarcating the regime transitions of Eq. ( 9).

This observation suggests the existence of a universal relation between z † 0 and µ.

To unveil it, we plot available data sets of threshold shear stress ln µ vs ln Re * as ln µ vs ln z † 0 by transforming ln Re * into ln z † 0 using universal Eq. ( 8). Carrying out this manipulation for fluids of densities spanning nearly six orders of magnitude, including water [START_REF] Guo | Empirical model for shields diagram and its applications[END_REF], air [START_REF] Pähtz | The cessation threshold of nonsuspended sediment transport across aeolian and fluvial environments[END_REF], and reduced-pressure gases with mean-free path smaller than the pore scale [START_REF] Swann | Experimentally-derived thresholds for windblown sand on Mars[END_REF][START_REF] Andreotti | A lower-than-expected saltation threshold at Martian pressure and below[END_REF], we confirm the density scaling

µ * ≡ µ × (ρ s /ρ) 1/3 (10)
that [START_REF] Pähtz | Scaling laws for planetary sediment transport from DEM-RANS numerical simulations[END_REF] identified, and we obtain the simple one-parameter correlation shown in Fig. 1C,

µ * = µ * ∞ × z 0 d/ , (11)
in which z 0∞ = d/ serves as reference for z 0 . Equivalently, τ c = τ * c∞ ×(ρ s -ρ) g z 0 (ρ/ρ s ) 1/3 , where τ * c∞ is the density-corrected threshold shear stress in the rough limit (Re * → ∞, ln z † 0 → 0), fitted to τ * c∞ = 2.35 ± 0.19 at 95% confidence within the fully-turbulent range Re * > -4 of Nikuradse's experiments. Here, the precision of Nikuradse's data implies that scatter is mainly inherited from the more challenging measurements of [START_REF] Shields | Application of similarity principles and turbulence research to bed-load movement[END_REF] and later authors.

Although one has not yet fully explained the exponent of (ρ/ρ s ) in Eq. ( 10), its value of 0.328 ± 0.016 fitted to Nikuradse's range is nearly indistinguishable from the (1/3) that [START_REF] Pähtz | Scaling laws for planetary sediment transport from DEM-RANS numerical simulations[END_REF] observed. Whereas Shields' static force balance only implied a sensitivity of τ c to the density difference (ρ s -ρ), the dependence of µ on the ratio (ρ/ρ s ) suggests that dynamical rates of momentum change in the fluid and solid phases are involved at the threshold [START_REF] Hjelmfelt | Motion of discrete particles in a turbulent fluid[END_REF], in addition to the mere weight and buoyancy that underpinned Shields' analysis.

Although some of Shields' original data involved ripples and bars, his correlation did not explicitly account for topography, in spite of the role altitude variations play in shaping the evolution of shear stress, static pressure and effective turbulent roughness along hills and bedforms [START_REF] Jackson | Turbulent wind flow over a low hill[END_REF][START_REF] Fourrière | Bedforms in a turbulent stream: formation of ripples by primary linear instability and of dunes by nonlinear pattern coarsening[END_REF][START_REF] Claudin | Basal pressure variations induced by a turbulent flow over a wavy surface[END_REF][START_REF] Jia | Hydrodynamic roughness induced by a multiscale topography[END_REF]. Neither did he consider that sediment porosity can produce a slip velocity at the surface [START_REF] Beavers | Boundary conditions at a naturally permeable wall[END_REF], thereby modifying the apparent turbulent roughness of a permeable bed or facilitate its mobilization through internal pressure gradients acting as a body force [START_REF] Musa | Pore pressure in a wind-swept rippled bed below the suspension threshold[END_REF]. In short, other characteristic distances, such as undulation wavelength or the square root of the bed permeability, can emerge. In their absence, d/ is the only available scale that makes z dimensionless in a fully-developed turbulent boundary layer on a uniform planar sediment surface prior to transport.

Meanwhile, Nikuradse's experiments with rough circular pipes may seem tenuously related to fully-developed turbulent boundary layers of indefinite thickness on flat ground, until one notices the independence of Eq. ( 8) on the outer scale of the flow. However, slip velocity and internal pressure gradients are features that Nikuradse's rough impenetrable pipes did not possess.

In summary, because Eq. ( 11) applies equally to the low Re * , intermediate and rough turbulent regimes of Eq. ( 9) without involving an outer scale, the uncanny resemblance of Eq. ( 9) to transport threshold regimes suggests that z 0 /(d/ ) governs Shields' Coulomblike friction µ exerted by the fluid on a uniform plane of sediment. However, in practical situations where the surface has long-range non-uniformity in topography or bumpiness, turbulent roughness may not be straightforward to predict. In this case, one must rise to the challenge of measuring it in the field to evaluate the stress threshold for particle entrainment.

Open Research Section

Historical data of [START_REF] Nikuradse | Strömungsgesetze in rauhen Rohren (Laws of flow in rough pipes[END_REF] and [START_REF] Shields | Application of similarity principles and turbulence research to bed-load movement[END_REF]

Figure 1 .

 1 Figure 1. (A) The historical data of Nikuradse (1933) for ln z † 0 vs ln Re * (bottom). Dashed and straight lines are Eqs. (8) and (9), respectively. (B) The compilation of Guo (2020) (triangles), Pähtz and Duràn (2018) (circles; simulations are darker), Swann et al. (2020) (diamonds) and Andreotti et al. (2021) (squares) for ln µ * vs ln Re * . (C) ln Re * has been converted to ln z † 0

  are available as tables and graphs in their original manuscripts. The workbook Nikuradse&ShieldsData.xlsx submitted as Supporting Information lists them all, together with other incipient stress data recently compiled by Pähtz and Duràn (2018), Guo (2020), Swann et al. (2020) and Andreotti et al. (2021).

Acknowledgments

We are grateful to Philippe Claudin, Bruno Andreotti, Lucie Guignier and Nicolas Le Dantec for illuminating discussions. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958, which allowed the authors to attend a program at the Kavli Institute for Theoretical Physics, and by the ISblue project, Interdisciplinary Graduate School for the Blue Planet (ANR-17-EURE-0015) that fostered related discussions on dune modeling.


Laura Frohn

Sandrine Skiba-Cassy

Frédéric Terrier

Pierre Aguirre

Cervin Guyomar

Sarah Maman Haddad

Julien Bobe

Christophe C Klopp

Benjamin Costas

Carla Teixeira Ferreira

A yeast product improves the efficiency of a processed animal protein diet in rainbow trout by modulating the gut response

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The total replacement of fishmeal (FM) in fish feeds is a major issue for the development of aquaculture. Research efforts have focused on alternative proteins such as plant proteins and processed animal proteins (PAP). However, fishmeal-free diets have not yet proven their efficiency. The use of new functional ingredients such as yeast products appears to be an interesting solution to improve the performance of these diets.

In this study, we compared a PAP-based diet (PAP, 17%) supplemented or not with a yeast extract (YE, 3%) to a fishmeal-based control diet (FM, 19%) in rainbow trout through a digestibility assay and a 12 weeks-growth trial. Growth performances, intestinal histology, immune and inflammatory parameters were analyzed, and mRNA expression on intestine and liver was determined by high throughput sequencing technology.

The diets differed only slightly in terms of nutrient digestibility coefficients which were overall high. Growth performances with PAP diets were lower compared to that obtained with a fishmeal-based diet. However, feeding trout a PAP diet with YE supplementation (PAP+YE) statistically allowed an improvement of animal growth. Moreover, feed intake appeared as the main driver of growth. In the intestine, dietary yeast extract supplementation increased villi size in the proximal part. Plasma immune parameters, including peroxidase and immunoglobulin levels, were significantly higher in fish fed the PAP+YE diet. We also observed that the diets affected the expression of a small number of genes in the gut and liver, mostly related to inflammation, iron transport, or oxidative stress.

In conclusion, the present study demonstrates the ability of yeast extract to improve the performance of fishmeal-free PAP diets in rainbow trout by potentially enhancing nutrient absorption capacity and modulating mechanisms related to gut health.

  by Lesaffre, 59700 Marcq-en-Baroeul, France Pinel Karine, Université de Pau & des Pays de l'Adour, E2S UPPA, NUMEA, 64310 Saint Pée-sur-Nivelle, France.