<figure><img src="image_1.png" /></figure>

HAL
open science

# Rotating convection with a melting boundary: An application to the icy moons

T. Gastine, B. Favier

► To cite this version:

T. Gastine, B. Favier. Rotating convection with a melting boundary: An application to the icy moons. Icarus, 2025, 429, pp.116441. 10.1016/j.icarus.2024.116441. hal-04946706

HAL Id: hal-04946706

https://hal.science/hal-04946706v1

Submitted on 14 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques 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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Research Paper

# Rotating convection with a melting boundary: An application to the icy moons

T. Gastine<sup>a, ∗</sup>, B. Favier<sup>b, ∗</sup>

<figure><img src="image_4.png" /></figure>

ARTICLE INFO

Keywords:

Icy satellites

Convection

Phase field

Numerical simulations

A B S T R A C T

A better understanding of the ice-ocean couplings is required to better characterise the hydrosphere of the icy moons. Using global numerical simulations in spherical geometry, we have investigated here the interplay between rotating convection and a melting boundary. To do so, we have implemented and validated a phase field formulation in the open-source code MagIC. We have conducted a parameter study varying the influence of rotation, the vigour of the convective forcing and the melting temperature. We have evidenced different regimes akin to those already found in previous monophasic models in which the mean axisymmetric ice crust transits from pole-ward thinning to equator-ward thinning with the increase of the rotational constraint on the flow. The derivation of a perturbative model of heat conduction in the ice layer enabled us to relate those mean topographic changes to the underlying latitudinal heat flux variations at the top of the ocean. The phase change has also been found to yield the formation of sizeable non-axisymmetric topography at the solid-liquid interface with a typical size close to that of the convective columns. We have shown that the typical evolution timescale of the interface increases linearly with the crest-to-trough amplitude and quadratically with the mean melt radius. In the case of the largest topographic changes, the convective flows become quasi locked in the topography due to the constructive coupling between convection and ice melting. The tentative extrapolation to the planetary regimes yields $\Theta(10^2 - 10^3)$ meters for the amplitude of non-axisymmetric topography at the base of the ice layer of Enceladus and $\Theta(10^3 - 10^4)$ meters for Titan.

# 1. Introduction

The presence of liquid water oceans has been evidenced on at least four icy satellites of our solar system: Europa, Ganymede, Titan and Enceladus (e.g. Soderlund et al., 2020). The determination of the structural properties of their hydrosphere is one of the main objective of the ongoing ESA's JUICE and NASA's Europa Clipper missions (Van Hoolst et al., 2024; Roberts et al., 2023). In this context, it is necessary to better understand the dynamics of the subsurface oceans and their interplay with their overlying ice shelves. Oceans of the icy moons are subjected to buoyancy forcing of thermal and solutal origins as well as mechanical forcing such as libration or tides (Soderlund et al., 2023). The aim of this study is to focus on thermally-driven flows and their interaction with the overlying ice crust.

Using global numerical simulations, previous parameter studies have delineated several dynamical physical regimes of rotating convection in spherical geometry (Gastine et al., 2016; Long et al., 2020). The relative influence of rotation is usually categorised in terms of

boundary regimes defined as combinations of the dimensionless governing parameters relevant to the fluid problem, namely the Ekman number *E*, the Rayleigh number *Ra* and the Prandtl number *Pr* (to be defined below). In the rapidly-rotating regime, the convective flow is governed by a quasi-geostrophic balance which imposes a strong alignment of the convective flows with the rotation axis. Multiple alternating zonal jets can form (Ashkenazy and Tziperman, 2021; Bire et al., 2022), reaching a typical width close to the so-called Rhines scale (Rhines, 1975; Cabanes et al., 2024). On the other hand, when the rotational constraint drops, the flow becomes three-dimensional and the mean zonal flows occupies a decreasing fraction of the kinetic energy content (Gastine et al., 2013; Yadav et al., 2016).

In addition to these zonal flow changes, the heat flux pattern also evolves with the rotational constraint (Gastine and Aurnou, 2023): heat flux at the outer boundary of the ocean may either peak near the equator (the so-called “equatorial cooling” regime) or close to the poles (the “polar cooling” regime) (Soderlund, 2019; Amit et al., 2020). The exact parameter combination of *Ra*, *E* and *Pr* which governs the transition

<footnote><sup>a</sup> Université Paris Cité, Institut de physique du globe de Paris, UMR 7154 CNRS, 1 rue Jussieu, F-75005, Paris, France</footnote>

<footnote><sup>b</sup> Aix Marseille Univ, CNRS, Centrale Med, IRPHE, Marseille, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: gastine@ipgp.fr (T. Gastine).</footnote>
between these two regimes is still a matter of debates and depends on the nature of the mechanical boundary condition (Kvorka and Čadek, 2022), the radius ratio of the subsurface ocean (Bire et al., 2022) and the coupling with the mean zonal flows (Gastine and Aurnou, 2023). Using stress-free boundary conditions, Kvorka and Čadek (2022) report a transition between equatorial and polar cooling when $Ra E^{12/7} Pr^{-1} \approx 1$, while Hartmann et al. (2024) rather favour $Ra E^{3/2} \sim 1$ when rigid boundaries are employed (see also Bire et al., 2022). Although the uncertainties remain sizeable, current estimates of the dimensionless numbers of the sub-glacial oceans of the icy moons would most likely place them in a weakly-rotating convection regime rather prone to polar cooling (Soderlund, 2019; Lemasquerier et al., 2023b).

In the case of a conducting ice layer, an increase of the heat flux at the solid-liquid interface goes along with a thinning of the overlying ice, such that the ice shell thickness is expected to be anti-correlated with the heat flux pattern atop the ocean (Kvorka et al., 2018; Kihoulou et al., 2023). This hypothesis prompts several authors to interpret observations such as the chaos terrain in Europa's equatorial regions (Soderlund et al., 2014), the polar depressions on Titan (Kvorka et al., 2018) or the poleward thinning of the ice crust on Enceladus (Čadek et al., 2019) in terms of the relative changes of the underlying convective heat flux (see also Kvorka and Čadek, 2024).

Recent large eddy simulations which account for the local changes of the melting temperature along the solid-liquid interface however challenged this interpretation (Ashkenazy and Tziperman, 2021; Kang et al., 2022). Because of the melting point dependence on pressure (e.g. Labrosse et al., 2018), the 20 km ice thickness variation between the South pole and the equator of Enceladus (Hemingway and Mittal, 2019) would approximately yield a 0.2 K latitudinal difference at the base of the ice crust (Lawrence et al., 2024). This thermal gradient could drive large-scale baroclinic flows and hamper the convective transport in the polar regions (Kang, 2023), hence questioning the relation between the heating pattern and the ice geometry. Using 2-D Cartesian models in which the oceanic flows are driven by the latitudinal variations of the melting temperature at the solid-liquid interface, Kang and Jansen (2022) suggest that the heat transport is more efficient on larger icy moons because of their higher gravity. Corresponding equilibrated ice shells are then expected to become flatter for increasing body sizes. Variations in salinity due to melting or freezing of water (Ashkenazy and Tziperman, 2021; Wong et al., 2022) or heat flux heterogeneities at the ocean's base or at the moon's surface (Terra-Nova et al., 2023; Lemasquerier et al., 2023b) are additional physical ingredients likely to modify the interplay between oceanic flows and the overlying ice layer. These competing interpretations also stem from the challenge to draw reliable scaling laws able to bridge the gap between the control parameters accessible to current numerical models and the relevant planetary regime (Jansen et al., 2023; Cabanes et al., 2024).

In current numerical simulations dedicated to the influence of ice thickness variations on the oceanic flows, the ice crust is however assumed to be static (e.g. Kang et al., 2022; Kang, 2023). This prohibits the dynamical generation of topographic features associated with the turbulent flows. The purpose of the present study is precisely to focus on the interplay between rotating convection in spherical geometry and generation of topography by melting or freezing of the overlying ice crust. To do so, we consider a phase field formulation which allows to model a two-phase fluid problem on a single fixed-grid domain (e.g. Beckermann et al., 1999). This type of approach is a smoothed approximation of the phase change which greatly simplifies the numerical implementation and has been shown to converge towards the exact moving-boundary formulation (e.g. Hester et al., 2020). Among many examples, let us cite the successful application of this method to the study of Rayleigh-Bénard convection (RBC) with a melting boundary in Cartesian geometry (Favier et al., 2019; Yang et al., 2023c), the generation of topography in a turbulent shear flow (Couston, 2021; Perissutti et al., 2024), the formation of pinnacles or scallops in melting ice (Weady et al., 2022), or the influence of the aspect ratio of an

iceberg on its melting (Hester et al., 2021a). All of these examples showed complex dynamical interactions between the convective flow pattern and the morphology of the phase change interface which can only be captured with an explicit treatment of its dynamics. To date, there is a scarce number of studies dedicated to rotating convection with a dynamical phase change boundary. Using numerical simulations in Cartesian geometry, Ravichandran and Wettlaufer (2021) examine the melting of a solid layer above a convecting liquid domain rotating about its vertical axis. They report on quasi-steady states in which the convective columns are locked in the topographic troughs and crests of the solid-liquid interface.

To our best knowledge, our present study is the first one to consider a phase field formulation applied to convection in a global spherical geometry. We would like to stress that this work should be regarded as a first incremental step in the process of improving our understanding of the dynamical interplay between oceanic flows and the overlying ice layers. As such, several effects relevant to the sub-glacial oceans of the icy moons, such as salinity, creep of ice or the pressure dependence of melting temperature, have been ignored in the present study.

The paper is organised as follows. The physical model including the phase field formulation and its numerical implementation and validation is discussed in Section 2. Results of the numerical simulations are described in Section 3, where we have split the analysis in terms of a mean-field approach, discussing first the axisymmetric topography and then the non-axisymmetric roughness. Tentative geophysical estimates are discussed in Section 4 before concluding in Section 5.

## 2. Model and methods

### 2.1. A phase field model

We consider a spherical shell gap of inner radius $r_i$ and outer radius $r_o$ with $d = r_o - r_i$ and $\eta = r_i/r_o$ which rotates with a constant rotation rate $\Omega$ about the z-axis. Convection is enforced by maintaining a fixed temperature contrast $\Delta T$ between both boundaries. We explore the uneven generation of topography associated with the freezing and melting which occurs at the fluid-solid interface for a melting temperature comprised between the imposed temperatures at the top and bottom boundaries. The solid phase is motionless and located in the outer part of the spherical volume. In the following, we consider a dimensionless formulation of the Navier-Stokes equations under the Boussinesq approximation with a constant kinematic viscosity $\nu$ and thermal diffusivity $\kappa$. We employ the viscous diffusion time $d^2/\nu$ as the reference time scale and the imposed temperature contrast $\Delta T$ as the temperature scale. To model the phase changes, we adopt the phase field formulation by Beckermann et al. (1999) combined with a volume-penalisation technique (Angot et al., 1999; Hester et al., 2021b). Practically, this method involves the time integration of a continuous scalar quantity $\phi$ which continuously varies from 0 in the liquid phase to 1 in the solid phase. A small dimensionless parameter $\epsilon = \lambda/d$, usually termed the Cahn number, then defines the ratio between the microscopic thickness of the transition between the two phases $\lambda$ and the macroscopic domain size which is here the shell gap $d$. Phase field methods represent a smoothed formulation of phase changes which are easier to implement numerically, especially when using pseudo-spectral methods, and converge to the exact moving boundary formulation in the limit of vanishing $\epsilon$ (e.g. Caginalp, 1986). Using the model recently derived by Hester et al. (2020), and used in many subsequent studies (e.g. Weady et al., 2022; Yang et al., 2023c), the governing equations for the velocity field $u$, the temperature $T$ and the phase field $\phi$ then read

$$ \nabla \cdot u = 0, \qquad (1) $$

$$ \frac{\partial u}{\partial t} + u \cdot \nabla u + \frac{2}{E} e_z \times u = - \nabla p + \frac{Ra}{Pr} g e_r + \nabla^2 u - \frac{1}{\tau_p \epsilon^2} \phi u, \qquad (2) $$
$$ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T = \frac{1}{Pr} \nabla^2 T + St \frac{\partial \phi}{\partial t}, \qquad (3) $$

$$ \frac{5}{6} St Pr \frac{\partial \phi}{\partial t} = a \nabla^2 \phi - \frac{1}{\epsilon^2} \phi(1 - \phi) [a(1 - 2\phi) + T - T_M]. \qquad (4) $$

In the above expressions $e_r$ and $e_z$ denote the unit vectors in the radial and axial directions while $g = r/r_o$ is the dimensionless gravity profile. $T_M$ is the dimensionless melting temperature, with $0 < T_M < 1$ delimiting the regimes for which both liquid and solid phases are present within the spherical shell. In the limit of vanishing phase-change interface $\epsilon$, Eq. (4) effectively imposes the so-called Stefan conditions (Huppert, 1990; Worster, 2000), which relate the interface velocity to the heat flux differences at the solid-liquid interface and ensure that the phase change interface is isothermal according to

$$ \mathbf{n} \cdot [\nabla T^{(S)} - \nabla T^{(L)}] = St Pr \mathbf{v} \cdot \mathbf{n}, \qquad (5) $$

$$ T = T_M, \qquad (6) $$

where $\mathbf{n}$ is the unit vector normal to the interface, $\mathbf{v}$ is the interface velocity and the superscripts (S) and (L) respectively correspond to the solid and liquid phases. Within the phase field formulation, the interface is implicitly defined by the isosurface $\phi = 1/2$.

For simplicity, our model assumes that thermal diffusivity and density are the same in both phases and we neglect any compositional effects (such as the dependence of the melting temperature on salinity for example). Density could increase by about 20% in Ganymede's hydrosphere, most of the contrast being accommodated in the fluid phase with values ranging from 950 kg/m³ in the ice crust to about 1200 kg/m³ at the ocean's base (e.g. Journaux et al., 2020, their Table 2). Bearing in mind the possible unknowns in the transport properties in the interiors of the icy moons, thermal diffusivity of ice is expected to be larger than the one of liquid water by about a factor 5–10 covering the range $\kappa \sim 10^{-7} - 10^{-6}$ m<fnref n="2" />/s (e.g. Abramson et al., 2001; Vance et al., 2018). Additionally, the melting temperature does not depend on pressure, an effect that is potentially important for many deep fluid systems from sub-glacial lakes on Earth (e.g. Couston, 2021) to oceans in icy satellites (e.g. Labrosse et al., 2018; Kang, 2023). While such thermobaric and salinity effects can be incorporated in our approach (see Hester et al., 2020; Yang et al., 2023b, for recent examples involving salinity), we choose to ignore them in this preliminary study and defer the question of their contributions to future studies.

The set of Eqs. (1)–(4) is governed by four physical dimensionless numbers: the Ekman number $E$, the Rayleigh number $Ra$, the Prandtl number $Pr$ and the Stefan number $St$ expressed by

$$ E = \frac{v}{\Omega d^2}, \quad Pr = \frac{v}{\kappa}, \quad Ra = \frac{\alpha g_0 \Delta T d^3}{\nu \kappa}, \quad St = \frac{L}{c_p \Delta T}. \qquad (7) $$

In the above expressions, $g_0$ denotes the gravity at the outer boundary, $c_p$ corresponds to the heat capacity, $\alpha$ to the thermal expansion coefficient and $L$ to the latent heat per unit mass associated with the solid-liquid transition.

The phase field formulation and the penalty method involve several additional dimensionless parameters compared to more classical monophasic models. The phase-change interface thickness controlled by the Cahn number $\epsilon$ is chosen to be smaller than the thickness of the Ekman boundary layers, which is the smallest relevant physical scale close to the phase change interface in our problem. This choice, while very constraining numerically, was motivated by systematic convergence tests carried out by gradually reducing $\epsilon$. While larger values of the interface thickness might be relevant for modelling mushy layers (Le Bars and Worster, 2006), we nevertheless opt for this conservative approach in order to ensure a meaningful comparison with monophasic simulations. In addition, our model also involves two other dimensionless quantities: a volume-penalisation parameter $\tau_p$ and a parameter related to the phase field model, $a = \gamma/\lambda\Delta T$, where $\gamma$ expresses the curvature dependence of the melting temperature expressed in units of $Km$ (see Beckermann et al., 1999). This

curvature effect is however very small in practice and only relevant when considering dendritic growth or other microscopic phenomena. To satisfy the standard isothermal Stefan condition (6), $\gamma$ must remain as small as possible, which is the case provided that the dimensionless parameter $a$ remains of order one while $\epsilon \ll 1$ (Hester et al., 2020). In practice, following Couston (2021) and Yang et al. (2023c), we adopt $a = 1$ for all the simulations considered in this study. Concerning the penalisation term $-\phi u/(\tau_p \epsilon^2)$ in Eq. (2), its role is to exponentially attenuate the velocity inside the solid phase, effectively treating it as a porous medium using a Darcy–Brinkman model (Le Bars and Worster, 2006). Although an optimal value for the parameter $\tau_p$ can be derived in simple configurations (Hester et al., 2021a), we follow a different approach in our setup which involves rotation and buoyancy force. The value of $\tau_p$ is practically adjusted for each simulation to ensure that the kinetic energy density of the solid phase remains smaller than that of the fluid phase by a factor of at least $10^{-4}$.

Considering a simplified model without rotation and buoyancy, Hester et al. (2020) report on the second-order convergence of the phase-field equation (4) with respect to the interface thickness $\epsilon$ towards the actual Stefan problem (5)–(6). In Appendix B, we examine the convergence properties of the phase field formulation by defining a weakly-supercritical quasi-stationary benchmark for rotating convection in spherical shell with a melting boundary. The convergence of our phase field formulation is found to be more complex than anticipated by Hester et al. (2020), with a gradual change in slope between first and second orders, similar to previous findings by Favier et al. (2019). This difference is attributed to the influence of rotation and buoyancy, two physical effects overlooked by Hester et al. (2020) in their analysis.

We adopt the spherical coordinates $(r, \theta, \phi)$ and assume rigid mechanical boundary conditions held at constant temperatures with $u(r = r_i) = u(r = r_o) = 0$ and $T(r = r_i) = 1$ and $T(r = r_o) = 0$. We assume Dirichlet boundary conditions on the phase field with the outer (inner) boundary in the solid (liquid) phase, i.e. $\phi(r = r_o) = 1$ and $\phi(r = r_i) = 0$.

## 2.2. Numerical methods

The numerical models considered in this study have been computed using the open source code MagIC<fnref n="1" /> combined with the SHTns<fnref n="2" /> library to handle the spherical harmonic transforms (see Schaeffer, 2013). To enforce the divergence-free constraint on the velocity field, $\mathbf{u}$ is expressed in terms of poloidal and toroidal potentials following

$$ \mathbf{u} = \mathbf{u}_P + \mathbf{u}_T = \nabla \times (\nabla \times W \mathbf{e}_r) + \nabla \times Z \mathbf{e}_r. $$

The quantities $W$, $Z$, $p$, $T$ and $\phi$ are expanded in spherical harmonics up to a degree and order $\ell_{\max}$ in the horizontal directions and in Chebyshev polynomials up to a degree $N_r$ in the radial direction. We employ the pseudo-spectral formulation of the equations introduced by Glatzmaier (1984) in which the nonlinear terms are handled in the physical grid space and then transformed to the spectral space using spherical harmonic transforms (for details, see for instance Christensen and Wicht, 2015).

The equations are advanced in time using implicit-explicit (IMEX) time schemes which handle the non-linear terms and the Coriolis force explicitly and the linear terms implicitly. The explicit treatment of the additional volume penalisation term entering Eq. (2) yields an extra constraint on the time step size $\delta t$ (Kolomenskiy and Schneider, 2009) compared to classical convection problems, such that

$$ \max(\delta t) < C \tau_p \epsilon^2, \qquad (8) $$

where $C$ is a Courant factor that depends on the IMEX time integrator. In this study, we employ the third order backward difference scheme (SBDF3) for most of the numerical simulations and the third-order

<footnote n="1">¹ https://github.com/magic-sph/magic.</footnote>

<footnote n="2">² https://gricad-gitlab.univ-grenoble-alpes.fr/schaeffn/shtns.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 1. (a) Example of a steep function $f(r) = \frac{1}{2}(1 + \tanh \frac{r_{M}}{2\epsilon})$ represented using the Gauss–Lobatto grid (GL), the KTE mapping by Kosloff and Tal-Ezer (1993) with $\alpha_1 = 0.993$, the JVH mapping by Jafari-Varzaneh and Hosseini (2015) and the BT mapping by Bayliss and Turkel (1992) with $\alpha_1 = 40$ and $\alpha_2 = 0.4$ for $\epsilon = 3 \times 10^{-3}$ and $r_M = 4.7$. (b) Absolute value of the coefficients of the discrete cosine transform of $f$ as a function of the degree of the Chebyshev polynomial $N$. (c) Grid spacing $\delta r$ as a function of $r$ for the same mappings.</figcaption></figure>

diagonally-implicit Runge–Kutta schemes ARS343 from Ascher et al. (1997) for the most turbulent configurations (for a comparison for convection problems, see Gopinath et al., 2022).

The phase field formulation involves large gradients of $\phi$ localised around the melting radius $r_M$. Here we employ a Chebyshev collocation technique combined with mapping functions capable of handling such steep fronts. Practically, the radial grid points are defined by

$$r_k = \frac{1}{2}(r_o - r_i)F(x_k) + \frac{1}{2}(r_o + r_i) \quad \text{with } k \in [1, N_r]$$

where $x_k = \cos[(k-1)\pi/(N_r - 1)]$ are the Gauss–Lobatto (GL) nodes and $F$ a mapping function. Of particular interest are the mapping functions which allow to refine the resolution around localised regions of rapid variations, such as the ones introduced by Bayliss and Turkel (1992) and Jafari-Varzaneh and Hosseini (2015) (hereafter respectively BT and JVH). The latter one has in addition been designed to minimise Gibbs phenomena when solving Allen–Cahn type equations. For the definitions of the mapping functions, see Appendix A. The melt radius $r_M$ is however susceptible to substantially vary over the horizontal directions in the case of large topography, making the tuning of the mapping parameters difficult. In that case, it is also worthy considering the mapping by Kosloff and Tal-Ezer (1993) (hereafter KTE) which simply redistributes the native Gauss–Lobatto points more evenly in the bulk of the domain.

As an illustrative example, Fig. 1(a) shows the radial profile of a steep tanh function centred around the radius $r_M$ (which is typically the kernel used in many phase field formulations, see Beckermann et al., 1999) using the GL and KTE grids with $N_r = 193$ and using the JVH and BT mappings with $N_r = 97$, while Fig. 1(b) shows the corresponding Chebyshev spectra as a function of the degree of the Chebyshev polynomial $N$. To stress the differences in terms of grid resolution, Fig. 1(c) shows the grid spacing $\delta r_k = |r_{k+1} - r_k|$ for $k \in [1, N_r - 1]$ as a function of $r$. GL and KTE are coarse with only 3-4 collocation points distributed around the steep changes of the function. Using the BT or the JVH mappings allows to reduce the grid spacing by a factor five around the transition radius $r_M$ for half the number of radial grid points used in the GL or KTE mappings. Comparison of the spectra in Fig. 1(b) also shows that a much faster convergence of the collocation method towards machine precision is achieved for the former two mappings. Of practical interest for fluid problems with boundary layers, we also note that the JVH mapping exhibits a slightly slower convergence than the BT mapping but retains a better resolution (i.e. a smaller grid spacing) at the spherical shell boundaries. The chosen value of the interface thickness $\epsilon$ yields a constraint on the maximum eligible grid spacing (Favier et al., 2019)

$$\max(\delta r) < 2\epsilon. \qquad (9)$$

This constraint combined with the time-step restriction due to the volume-penalisation (8) makes the computation of phase field models much more numerically-demanding than the classical rotating convection problems. For $\epsilon \approx 10^{-3} - 3 \times 10^{-3}$ considered here, the radial resolution is about two to three times greater in phase field models than in their monophasic counterparts to satisfy Eq. (9), while the time step size $\delta t$ is about a factor ten smaller to fulfil Eq. (8). This limits the range of physical parameters $E$ and $Ra$ accessible to global numerical modelling.

## 2.3. Diagnostics

The interface between the solid and liquid phases depends on the angular directions $(\theta, \varphi)$ and evolves in time. Following Favier et al. (2019) and Yang et al. (2023c), we choose to define it by

$$\phi(r = r_M, \theta, \varphi, t) = \frac{1}{2}. \qquad (10)$$

In the following, we employ angular brackets for spatial averaging and overbars for time-averaging such that

$$\langle f \rangle_V = \frac{1}{V} \int_V f dV, \qquad \langle f \rangle_S = \frac{1}{4\pi} \int_0^\pi \int_0^{2\pi} f \sin \theta d\theta d\varphi, \\ \langle f \rangle_\varphi = \frac{1}{2\pi} \int_0^{2\pi} f d\varphi, \qquad \bar{f} = \frac{1}{\tau_{\text{avg}}} \int_{t_0}^{t_0+\tau_{\text{avg}}} f dt,$$

where $t_0$ and $\tau_{\text{avg}}$ respectively denote the beginning and the width of the time-average interval, and $V$ is the spherical shell volume. For clarity, we define in the following the time-averaged mean melt radius

$$\xi_M = \langle r_M \rangle_S, \qquad (11)$$

and the corresponding thicknesses of the solid and liquid phases by

$$h_S = r_o - \xi_M, \quad h_L = \xi_M - r_i. \qquad (12)$$

The corresponding volumes are accordingly expressed by

$$V_S = \int_0^{2\pi} \int_0^{\pi} \int_{r_M(\theta, \varphi)}^{r_o} r^2 \sin \theta dr d\theta d\varphi, \quad V_L = V - V_S. \qquad (13)$$

For an easier comparison with standard monophasic convective models, we define effective quantities based on the actual mean thickness and temperature contrast of the fluid layer:

$$\eta_{\text{eff}} = \frac{r_i}{\xi_M}, \quad Ra_{\text{eff}} = Ra g_m \Delta T_M h_L^3, \quad E_{\text{eff}} = \frac{E}{h_L^2}, \qquad (14)$$

where $\Delta T_M = 1 - T_M$ is the temperature contrast in the liquid and $g_m = \xi_M/r_o$ is the gravity at $\xi_M$.
The kinetic energy content per unit volume can be decomposed in poloidal and toroidal contribution following

$$E_K = \frac{1}{2} \left( \langle \mathbf{u}_P^2 \rangle_V + \langle \mathbf{u}_T^2 \rangle_V \right) = \sum_{\ell=1}^{\ell_{\max}} E_{\ell}^P + \sum_{\ell=1}^{\ell_{\max}} E_{\ell}^T, \qquad (15)$$

where $E_{\ell}^P$ ($E_{\ell}^T$) respectively denote the poloidal (toroidal) kinetic energy content at the spherical harmonic degree $\ell$. The convective flow velocity is accordingly measured by a Reynolds number defined on the averaged liquid thickness $h_L$

$$Re_L = Re \sqrt{\frac{V}{V_L}} h_L, \quad Re = \sqrt{2E_K}. \qquad (16)$$

Heat transfer is characterised by the Nusselt number here defined at the spherical shell boundaries

$$Nu = \eta \frac{d\langle\overline{T}\rangle_S}{dr} \Big|_{r=r_i} = \frac{1}{\eta} \frac{d\langle\overline{T}\rangle_S}{dr} \Big|_{r=r_o}, \qquad (17)$$

where the factors involving the radius ratio $\eta$ come from the diffusive temperature gradient. Again, for comparison purpose with standard models without phase change, it is insightful to also define an equivalent Nusselt number

$$Nu_{\text{eff}} = Nu \frac{h_L}{\Delta T_M} \frac{\eta_{\text{eff}}}{\eta}. \qquad (18)$$

Following Schwaiger et al. (2019), the typical convective flow length-scale is defined as the time-average of the peaks of the poloidal kinetic energy spectra

$$\ell_U = \overline{\argmax_{\ell} E_{\ell}^P(t)}. \qquad (19)$$

To define a typical lengthscale for the topography of the solid–liquid interface, we proceed with a truncated spherical harmonic expansion such that

$$r_M(\theta, \varphi, t) \approx \sum_{\ell=0}^{\ell_{\max}} \sum_{m=-\ell}^{\ell} \xi_{\ell m}(t) Y_{\ell m}(\theta, \varphi), \qquad (20)$$

where $Y_{\ell m}$ denotes the spherical harmonic function of degree $\ell$ and order $m$. Similarly to $\ell_U$, we then define $\ell_\xi$ as the time average of the peaks of the spectra defined by

$$\ell_{\xi} = \overline{\argmax_{\ell} E_{\ell}^{\xi}(t), E_{\ell}^{\xi}(t) = \sum_{m=-\ell, m \neq 0}^{\ell} |\xi_{\ell m}|^2}. \qquad (21)$$

Corresponding lengthscales at any radius read $L_{[U,\xi]}(r) \approx \pi r / \ell_{[U,\xi]}$ (e.g. Backus et al., 1996, §3.6.3).

## 2.4. Parameter coverage

We carried out 26 numerical simulations with $\eta = 0.8$, $Pr = 1$ and $St = 1$ divided in four groups with fixed parameter pairs ($E, Ra$) and increasing values of the melting temperature $T_M$. The Stefan number $St$ defined by Eq. (7) controls the time-scale separation between the temporal evolution of the phase-change interface and that of the flow. While sub-glacial oceans are probably characterised by $St \gg 1$ (using $L \sim 3 \times 10^5$ J/kg, $c_p = 4 \times 10^3$ J/kg/K from Journaux et al. (2020) and a temperature gradient of 0.04 K/km (Vance et al., 2018) yields $St \approx 40$ for Europa), we choose to consider the less numerically-demanding value $St = 1$ which speeds up the melting and freezing dynamics leading to faster transients. Recent studies on thermal convection interacting with a phase change boundary have shown that the value of the Stefan number has only a marginal impact on the final quasi-stationary equilibrium (Rabbanipour Esfahani et al., 2018; Purseed et al., 2020; Yang et al., 2023c). Similarly, we choose to fix the value of the Prandtl number to $Pr = 1$ for numerical convenience, knowing that oceans are more likely to be characterised by $Pr \sim \mathcal{O}(10)$.

Spatial resolutions, control parameters and main diagnostics of these simulations are listed in . A good indicator of the rotational

<figure><img src="image_6.png" /><figcaption>Fig. 2. Numerical simulations carried out in this study located in a parameter space constructed using the effective Ekman and Rayleigh numbers $E_{\text{eff}}$ and $Ra_{\text{eff}}$ defined in (14). Symbols with a black rim correspond to the configurations with no solid phase (i.e. $T_M = 0$). Each type of symbols corresponds to fixed values of $Ra$ and $E$ and increasing values of $T_M$. The errorbars correspond to the control parameters constructed using $\min_{\theta,\phi}(\overline{r}_M)$ and $\max_{\theta,\phi}(\overline{r}_M)$. They highlight the maximum topography of each model. For comparison purpose, the critical Rayleigh numbers for onset of spherical shell convection with $r_i/r_o = 0.8$ come from the study by Barik et al. (2023). The solid black line corresponds to $R_{o_c} = 1$ (Eq. (22)).</figcaption></figure>

constraint on the convective flow is provided by the convective Rossby number

$$R_{o_c} = \sqrt{\frac{Ra_{\text{eff}} E_{\text{eff}}^2}{Pr}}, \qquad (22)$$

introduced by Gilman (1977). Fig. 2 shows the location of the simulations in terms of their effective Ekman and Rayleigh numbers once the system has reached a quasi-stationary state. Estimates of the transport properties in the subsurface oceans of the icy satellites yield Ekman and Rayleigh numbers that range from $E \sim 10^{-10}$ and $Ra \sim 10^{16}$ for smaller moons to $E \sim 10^{-13}$ and $Ra \sim 10^{24}$ for the larger ones (for details, see e.g. Table 1 in Soderlund, 2019), which is way outside the ranges considered in Fig. 2. We nonetheless see that our parameter choice yields configurations where $R_{o_c}$ exceeds unity with $E = 10^{-3}$ and $Ra \in [3 \times 10^6, 10^7, 3 \times 10^7]$ (triangles) and configurations where $R_{o_c} < 1$ with $E = 3 \times 10^{-4}$ and $Ra = 1.2 \times 10^7$ (squares) and $E = 3 \times 10^{-5}$ and $Ra = 2.5 \times 10^8$ (circles), a physical regime expected to be relevant to the subsurface oceans of the icy moons (e.g. Soderlund, 2019, their Fig. 1). For each group of simulations, an increase of $T_M$ goes along with a decrease of the effective Rayleigh number $Ra_{\text{eff}}$ and an increase of the effective Ekman number $E_{\text{eff}}$, as expected from the gradual decrease of the fluid layer thickness as $T_M$ increases. The errorbars attached to each symbol correspond to the changes of the effective control parameters when considering the extrema of the interface $\min_{\theta,\phi}(\overline{r}_M)$ and $\max_{\theta,\phi}(\overline{r}_M)$ to evaluate them. Overall, increasing $T_M$ at fixed ($E, Ra$) yields a parameter path along which (i) $R_{o_c}$ decreases, hence strengthening the rotational constraint; (ii) the supercriticality of the convective flow drops; and (iii) the topographic changes of the interface increase.

Most of the numerical simulations have been initiated from a diffusive state and a random noise temperature perturbation combined with a spherically-symmetric phase field tanh profile centred around $T_M$. Simulations are run until a quasi-stationary state is reached at which point diagnostics are computed. For the most demanding configurations, the simulations have been computed by restarting from setups
<figure><img src="image_7.png" /><figcaption>Fig. 3. Three-dimensional renderings of snapshots of four selected simulations. For each simulation, the inner sphere shows the temperature at $r = r_i + 0.02$ atop the inner thermal boundary layer, while the outer surface corresponds to the melting radius $r_M$. Equatorial and meridional slices show the temperature in the liquid and solid phase with two different separated colourmaps.</figcaption></figure>

with neighbouring parameters. While bistability is known to occur in convective systems with phase change (Purseed et al., 2020; Yang et al., 2023a), we have not observed such behaviour in our simulations.

## 3. Results

We first focus on selected simulations to highlight the interplay between the convective flow and the shape of the solid–liquid interface. Fig. 3 shows three-dimensional renderings of the temperature field for four selected simulations with an increasing rotational constraint on the flow, while Fig. 4 shows the corresponding Hammer projections of the melt radius $r_M$.

The first case with $E = 10^{-3}$, $Ra = 3 \times 10^7$ and $T_M = 0.5$ (Figs. 3a and 4a) corresponds to weakly-rotating convection with $Ro_c \approx 4$. Similar to classical RBC, the thermal plumes are radially-oriented and the convective pattern at the edge of the inner thermal boundary layer forms a network of thin sheet-like upwellings surrounding broad downwellings. The solid–liquid interface is almost spherical, with little variations with $r_M \in [4.89, 4.94]$ (we recall that we use $\eta = 0.8$ so that $r \in [4, 5]$) and no preferred direction, indicating the weak rotational constraint. The second model with $E = 10^{-3}$, $Ra = 10^7$ and $T_M = 0.6$ (Figs. 3b and 4b) corresponds to $Ro_c \approx 2$. The convective pattern is more laminar than in the previous case due to a weaker effective Rayleigh

number but remains similar regarding the influence of rotation. The most striking difference comes from the solid–liquid interface which now features large-scale topographic changes aligned with the rotation axis outside the effective tangent cylinder, defined by the colatitudes $\theta = [\arcsin \eta_{\text{eff}}, \pi - \arcsin \eta_{\text{eff}}]$. Those features are reminiscent of the so-called “banana cells”, i.e the convective columns outside the tangent cylinder of rotating spherical shells (e.g. Busse, 1970; Simitev and Busse, 2003). For that particular case, the amplitude of topography from trough to crest in the equatorial plane reaches about 15% of the shell gap. Similar to previous findings in non-rotating convection (e.g. Rabbanipour Esfahani et al., 2018; Yang et al., 2023c), cold plumes detach from the thermal boundary layer in the cusps of the solid–liquid interface (equatorial cut in Fig. 3b). The influence of rotation is more pronounced in the third case with $Ro_c \approx 0.7$ (Figs. 3c and 4c). This manifests itself by a stronger alignment of the thermal plume with the rotation axis inside the tangent cylinder (meridional cut in Fig. 3c) and vortical structures at the connection points of the convective sheets (radial cut in Fig. 3c). The solid–liquid interface is made up of two main contributions: (i) an overall large-scale latitudinal contrast with thinner ice crust inside the effective tangent cylinder than outside; (ii) small-scale columnar corrugations aligned with the rotation axis outside the tangent cylinder. The decrease of the effective Ekman number compared to the previous case goes along with smaller scale columnar
<figure><img src="image_8.png" /><figcaption>Fig. 4. Hammer projections of snapshots of the melt radius $r_M$ for the four simulations shown in Fig. 3. In each panel, the dashed lines correspond to the location of the effective tangent cylinder, i.e. $\theta = [\arcsin \eta_{eff}, \pi - \arcsin \eta_{eff}]$.</figcaption></figure>

convection and weaker topographic changes in the equatorial region. This trend is confirmed in the last simulation with the smallest Ekman number and $Ro_c \approx 0.3$ (Figs. 3d and 4d). In this rotationally-constrained configuration, most of the topographic changes are axisymmetric with a thinner ice thickness at the equator than at the poles. The melt radius also shows secondary peaks right at the location of the tangent cylinder, a feature reminiscent of the local flux maxima observed there in rotating spherical shell convection (see Fig. 3 in Gastine and Aurnou, 2023).

Overall the increasing influence of rotation goes along with a topography which transits from columnar troughs and crests outside the effective tangent cylinder to an almost axisymmetric profile modulated by small-scale roughness. This prompts us to separately investigate the large-scale axisymmetric topography and the non-axisymmetric features. In the following, we hence disentangle the axisymmetric topography defined by

$$ \bar{\xi} = \langle \overline{r_M} \rangle_{\phi}, \qquad (23) $$

from the fluctuating non-axisymmetric patterns defined by the following standard deviation

$$ \xi' = \overline{\langle (r_M^2)_{\phi} - \langle r_M \rangle_{\phi}^2 \rangle}^{1/2}. \qquad (24) $$

### 3.1. Axisymmetric topography

Fig. 5(a) shows latitudinal profiles of the mean axisymmetric melt radius $\bar{\xi}$ for a series of simulations with $E = 3 \times 10^{-4}$ and $Ra = 1.2 \times 10^7$. An increase in $T_M$ goes along with thicker ice with increasing latitudinal contrasts. For an easier comparison, Fig. 5(b) shows the corresponding latitudinal variations of the ice thickness normalised by its mean value, i.e. $(r_o - \bar{\xi})/(r_o - \xi_M)$. Relative changes in the ice thickness are found to be quite similar from one simulation to another, ranging between -0.2 and 0.2. The axisymmetric topography evolves from configurations with thicker ice at the equator than at the poles for $T_M < 0.7$ to the opposite for $T_M > 0.75$. Postponing for now the question of the control parameters which govern this transition, we note a good correlation between the profiles shown in Fig. 5(b) and the outer boundary heat flux obtained in the monophasic rotating convection models by Amit et al. (2020) (their Fig. 7) and by Kvorka and Čadek (2022) (their Fig. 10).

This raises the question of the feedback of the topography onto the flow. To examine this issue, we have computed the equivalent 26 rotating spherical shell simulations without phase change adopting the effective control parameters $Ra_{eff}$, $E_{eff}$ and $\eta_{eff}$. In addition, we have also carried out 5 extra simulations without solid phase (i.e. $T_M = 0$) for each $(Ra, E)$ pair. Fig. 6 shows a comparison between the mean radial profiles of temperature (panel a) and kinetic energy (panel b) for numerical models with a phase change (solid lines) and their equivalent counterparts without (dashed lines). Typical of rotating convection, the temperature profile of the reference case with no solid phase ($T_M = 0$) is made up of three distinct parts: two thermal boundary layers which accommodate most of the temperature contrast and a quasi-linear temperature drop in the fluid bulk. The corresponding kinetic energy profile features two localised maxima which mark the location of the edges of the Ekman boundary layers. The increase of $T_M$ goes along with a gradual decrease of the kinetic energy content due to the shrinking of the fluid region. Temperature then follows a quasi-linear conducting profile in the solid part. Models with or without phase changes almost perfectly overlap for $T_M \le 0.75$. They depart more strongly for the largest melting temperature $T_M = 0.8$, due to the largest amplitude of topography in this configuration (recall Fig. 5a). In view of the good agreement between the radial profiles, the average rms properties of the convective flow are found to be little affected by the phase change (see).

Fig. 7 shows an additional comparison between models with and without phase change for two configurations with large melting temperature to enhance the topographic changes. Panels (a) and (b) correspond to a weakly-rotating configuration with $E = 10^{-3}$, $Ra = 3 \times 10^7$, $T_M = 0.8$ and $Ro_c \approx 3.1$. In this setup, the solid-liquid interface presents significant changes in the azimuthal direction. Cold plumes which detach from the upper boundary layer are mostly localised in the cusps of the interface, while the hot upwellings are clustered in fluid regions with a thinner ice crust. Despite those topographic changes, the mean azimuthal zonal flows $\langle \overline{u_{\phi}} \rangle_{\phi}$ appear strikingly similar between the phase field case and its monophasic counterpart. Equatorial zonal flows are retrograde close to the outer boundary, a typical feature of spherical-shell rotating convection when the convective Rossby number exceeds one (e.g Gilman, 1977; Aurnou et al., 2007; Gastine et al., 2013; Yadav et al., 2016). Zonal flow gradients are predominantly radial, with no marked alignment with the axis of rotation, indicating a weak rotational constraint. Their energetic content amounts to about
<figure><img src="image_9.png" /><figcaption>Fig. 5. (a) Time and azimuthal average of the melt radius $\bar{\xi}$ (see Eq. (23)) as a function of colatitude for a series of numerical simulations with $E = 3 \times 10^{-4}$ and $Ra = 1.2 \times 10^7$ and increasing $T_M$. (b) Normalised thickness of the solid phase $(r_o - \bar{\xi})/(r_o - \bar{\xi}_M)$ as a function of colatitude. The dotted lines in panel (a) correspond to the mean melt radius $\bar{\xi}_{M}$, while the vertical segments mark the location of the effective tangent cylinder.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Comparison between models with a phase change and their equivalent convective models without for simulations with $E = 3 \times 10^{-4}$ and $Ra = 1.2 \times 10^7$. (a) Time-averaged radial profiles of temperature for increasing values of the melting temperature $T_M$. (b) Time-averaged radial profiles of kinetic energy. For comparison purpose, the grey lines in panels (a) and (b) correspond to the setup with no solid phase (i.e. $T_M = 0$). The dotted vertical lines in panel (b) correspond to $\bar{\xi}_M$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7. Comparison between models with a phase change (upper half of each panel) with equivalent convective models without (lower halves). Panels (a) and (c) show snapshots of the temperature in the equatorial plane, while panels (b) and (d) show the time and azimuthal average of $u_\varphi$ in a meridional plane. Panels (a) and (b) correspond to a weakly-rotating configuration with $E = 10^{-3}$, $Ra = 3 \times 10^7$, $T_M = 0.8$, while (c) and (d) correspond to a rotationally-constrained setup with $E = 3 \times 10^{-5}$, $Ra = 3 \times 10^8$, $T_M = 0.8$. The dashed lines mark the location of the solid-liquid interface.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 8. (a) Time and azimuthal average of the relative heat flux (27) expressed at the top of the liquid phase for a numerical simulation with $E = 3 \times 10^{-4}$, $Ra = 1.2 \times 10^7$ and $T_M = 0.55$. (b) Same quantities for a numerical simulation with $E = 3 \times 10^{-4}$, $Ra = 1.2 \times 10^7$ and $T_M = 0.7$. In both panels, the dashed vertical lines correspond to the location of the tangent cylinder.</figcaption></figure>

4% of the total kinetic energy in that case. Conversely, Fig. 7(c) and (d) correspond to a rapidly-rotating configuration with $E = 3 \times 10^{-5}$, $Ra = 2.5 \times 10^8$, $T_M = 0.8$ which yields $Ro_c \approx 0.3$. In this setup, topographic changes happen mostly in the latitudinal direction with a thinner ice at the equator than at the poles. As such, the fluid domain in the equatorial plane is therefore slightly thicker than its equivalent without phase change (Fig. 7c). The typical size and shape of the convective plumes nevertheless show a good agreement between the two simulations. Zonal flow profiles (Fig. 7d) take the form of a pair of alternated retrograde and prograde jets mostly localised outside the tangent cylinder, very much alike those obtained in monophasic spherical shell convection with rigid boundaries and comparable convective Rossby numbers $Ro_c \approx 0.3$ (e.g. Cabanes et al., 2024, rightmost simulation in their Fig. 5). Their energetic content remains however quite weak, about 3% of the total kinetic energy, a value in line with those reported by Cabanes et al. (2024) for similar effective Ekman numbers $E_{eff} \approx 10^{-4}$ (see their Table 1). Despite the sizeable axisymmetric topographic variations, zonal flows appear little affected by the phase change. Assessing the feedback between topography and zonal jets would require reaching lower (larger) Ekman (Rayleigh) numbers such that they represent a greater share of the kinetic energy.

We finally examine the differences in terms of heat flux at the top of the fluid layer. To do so, we derive in Appendix D a perturbative model for the diffusion of heat in a solid phase with an upper spherical boundary held at a constant temperature $T = 0$ and a lower quasispherical boundary with topographic changes held at $T = T_M$. Provided the amplitude of topography remains small compared to the radius, i.e. the spherical harmonic expansion coefficients from Eq. (20) fulfil $|\xi_{\ell m}|/\xi_M \ll 1$, we show in Appendix D that the temperature gradient along the interface reads

$$ \partial Q(\theta, \varphi) \approx \sum_{\ell \neq 0, m} \frac{\xi_{\ell m}}{\xi_M} f_{\ell}(\eta_S) Y_{\ell m}(\theta, \varphi), \qquad (25) $$

with

$$ f_{\ell}(\eta_S) = \frac{\ell - 1 + (\ell + 2)\eta_S^{2\ell+1}}{1 - \eta_S^{2\ell+1}}, \qquad (26) $$

and $\eta_S = \xi_M/r_o$ is the mean radius ratio of the solid phase and $\partial Q$ characterises the relative heat flux changes

$$ \partial Q(r_M, \theta, \varphi) = \frac{\frac{\partial \bar{T}}{\partial r}(r_M, \theta, \varphi) - \frac{\partial T_0}{\partial r}(\xi_M)}{\frac{\partial T_0}{\partial r}(\xi_M)}, \qquad (27) $$

with

$$ \frac{dT_0}{dr} = -\frac{T_M}{h_S \eta_S}, \qquad (28) $$

the zeroth order diffusive temperature gradient in a spherical shell of mean inner radius $\xi_M$ and outer radius $r_o$.

To evaluate this first-order model, Fig. 8 shows a comparison of the relative variations of the axisymmetric heat flux $\langle \partial Q \rangle_{\varphi}$ for two numerical simulations with $E = 3 \times 10^{-4}$ and $Ra = 1.2 \times 10^7$ which only differ from their melting temperature with $T_M = 0.54$ and a thin ice thickness $h_S = 0.15$ (panel a) and $T_M = 0.8$ and a thick ice shell $h_S = 0.65$ (panel b). To also examine the influence of the topography on the flow, we include in Fig. 8 the heat flux profiles of the equivalent models without phase change. In that case, the relative heat flux variations are expressed by

$$ \langle \partial Q \rangle_{\varphi} = \frac{\frac{\partial \langle \bar{T} \rangle_{\varphi}}{\partial r}(r_o, \theta) - \frac{\partial \langle T \rangle_{S}}{\partial r}(r_o)}{\frac{\partial \langle \bar{T} \rangle_{S}}{\partial r}(r_o)}. \qquad (29) $$

The heat flux profiles shown in Fig. 8(a) look very much alike those obtained by Amit et al. (2020) in their so-called “polar-cooling” configurations (see the third panel in their Fig. 7). In this setup, the heat flux is almost constant outside the tangent cylinder and increases gradually inside to reach local maxima in the polar regions. We observe, on the one hand, an excellent agreement between the actual heat flux obtained in the phase field simulation and the theoretical model expressed in Eq. (27), and on the other hand, very similar profiles with or without phase change. This latter observation indicates a negligible influence of the axisymmetric topography on the heat flux, which comes to no surprise given the small variations of the ice shell thickness in this configuration ($\xi \in [4.84, 4.88]$, see Fig. 5a). Fig. 8(b) corresponds to a configuration with a smaller convective Rossby number $Ro_c \approx 0.67$. In that case, and in line with previous findings of monophasic rotating convection (see, e.g. the first panel of Fig. 7 in Amit et al., 2020), the heat flux peaks in the equatorial region, a hallmark of the so-called “equatorial cooling regime”. More pronounced differences between the monophasic and the phase field configurations are however observed: while the perturbative model still correctly accounts for the actual heat flux, the configuration without phase changes now yields larger latitudinal heat variations than observed in the case with a phase change. This clearly indicates that the purely spherical analogue cannot account for the heat flux changes for configurations with large topographic variations (here $\xi \in [4.17, 4.46]$, see Fig. 5a).

## 3.2. Non-axisymmetric roughness

We now turn our attention to the non-axisymmetric topography, characterised in terms of the standard deviation of $r_M$ along the longitudinal direction (Eq. (24)).

We have already seen in Fig. 4 significant regionalised differences in the amplitude and size of the interface roughness depending on the strength of the rotational constraint. To illustrate this phenomenon, Fig. 9 shows mean latitudinal contrasts of $\xi'$ for two series of simulations with increasing melting temperature and with either a weak rotational constraint ($E = 10^{-3}$, $Ra = 3 \times 10^7$, panel a) or a strong influence of rotation ($E = 3 \times 10^{-5}$, $Ra = 2.5 \times 10^8$, panel b). For the former series of simulations, an increase of $T_M$ goes along with a gradual increase of the interface roughness. The configurations which are the least influenced by rotation, i.e. $T_M \le 0.6$ and $Ro_c > 3.5$,
<figure><img src="image_13.png" /><figcaption>Fig. 9. Mean latitudinal profile of non-axisymmetric topography (as defined in Eq. (24)) for two series of simulations with $E = 10^{-3}$, $Ra = 3 \times 10^{7}$ (panel a), and with $E = 3 \times 10^{-5}$, $Ra = 2.5 \times 10^{8}$ (panel b). In each panel, the vertical segments mark the location of the tangent cylinder.</figcaption></figure>

feature almost no latitudinal variation and $\xi'$ does not exceed 1% of the shell gap. When $T_M \ge 0.75$, the amplitude of the non-axisymmetric topography raises by one order of magnitude in the low-latitude regions reaching $\max_{\theta} \xi' \approx 0.1$. This corresponds to the columnar topographic changes visible in Fig. 4(b). The second series of simulations features more undulating profiles of overall weaker amplitude (Fig. 9b). The amplitude of $\xi'$ follows a non-monotonic behaviour with the increase of $T_M$. Latitudinal changes appear to be correlated with the variations of the mean axisymmetric ice thickness with localised minima at the equator and at the location of the tangent cylinder (see Fig. 3d).

To explore the time variability of the solid-liquid interface, we show in Fig. 10 Hovmöller diagrams of the melt radius in the equatorial plane $r_M(\pi/2, \varphi, t)$ for four simulations with $E = 10^{-3}$, $Ra = 3 \times 10^7$ and increasing values of $T_M$. The first two cases ($T_M = 0.5$ and $T_M = 0.6$) present rapid temporal variability on time scales of less than one tenth of the viscous diffusion time. Topography in the equatorial plane mostly drifts westward and the crest-to-trough amplitude is about 2-3% of the spherical shell gap. In the second case with $T_M = 0.6$, larger longer-lived topographic features with a larger amplitude are occasionally observed (e.g. $t \approx 5.2$ and $\varphi \approx 3.5$). This trend becomes more pronounced in the third and fourth cases with $T_M = 0.69$ and $T_M = 0.8$. Topographic contrast increases and remains stable and coherent over timespans that exceed the viscous diffusion time. We only observe a very slow westward drift of the topographic features. These latter two configurations are reminiscent of rotating convection models in Cartesian geometry by Ravichandran and Wettlaufer (2021) in which convective features are locked within the interface topography. In contrast with the moderate feedback of the axisymmetric topography on the convective flow discussed before, non-axisymmetric roughness hence yields significant differences with classical rotating convection without phase change. We recall that the Stefan number $St$ is here fixed to unity across all simulations. While systematically varying $St$ is outside the scope of the current study, it would influence the characteristic timescale of topographic features with larger Stefan numbers leading to slower topography dynamics.

To quantify the typical time of variation of the topography, we define the following auto-correlation function for each location $(\theta, \varphi)$:

$$C(\theta, \varphi, \tau) = \frac{\overline{r_M(\theta, \varphi, t+\tau)r_M(\theta, \varphi, t)}}{\overline{r_M^2(\theta, \varphi, t)}}.$$

The correlation time of the topography at each point $\tau_{\xi}(\theta, \varphi)$ is then defined as the full width at half maximum of $C$, i.e.

$$C(\theta, \varphi, \tau_{\xi}) = \frac{1}{2}. \qquad (30)$$

Fig. 11 shows an illustration of the time-averaged non-axisymmetric topography (panel a) alongside the correlation time $\tau_{\xi}(\theta, \varphi)$ for a model

<figure><img src="image_14.png" /><figcaption>Fig. 10. Longitudinal Hovmöller diagrams for the melt radius in the equatorial plane $r_M(\theta = \pi/2, \varphi, t)$ for four numerical simulations with $E = 10^{-3}$ and $Ra = 3 \times 10^7$ with $T_M = 0.5$ (panel a), $T_M = 0.6$ (panel b), $T_M = 0.69$ (panel c) and $T_M = 0.8$ (panel d).</figcaption></figure>

with $E = 10^{-3}$, $Ra = 10^7$ and $T_M = 0.6$. Given that the topography slowly evolves in that case, the time-averaging involved in panel (a) has been conducted over the typical thermal diffusion time of the ice crust $Prh_S^2 \approx 3.4 \times 10^{-2}$, a value that is about one order of magnitude below the typical time of variation of the topography in that case. For both diagnostics, important regionalised differences are clearly visible
<figure><img src="image_15.png" /><figcaption>Fig. 11. (a) Hammer projection of the time-averaged non-axisymmetric topography defined by $|r_M - \bar{\xi}|$ for a numerical model with $E = 10^{-3}$, $Ra = 10^7$ and $T_M = 0.6$. Time averaging has been conducted over one thermal diffusion time of the ice layer, i.e. $Prh_L^2 \approx 3.4 \times 10^{-2}$. (b) Corresponding Hammer projection of the correlation time of the topography expressed by Eq. (30). In each panel, the dashed lines mark the location of the effective tangent cylinder.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. Average of the correlation time of the topography (Eq. (30)) as a function of the theoretical scaling (32). The solid line corresponds to a polynomial fit to the data. Symbols convey the same meaning as in Fig. 2.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Comparison between the convective flow wavenumber $\ell_U$ and the roughness wavenumber $\ell_\xi$. The solid line corresponds to the equality $\ell_\xi = \ell_U$. Symbols carry the same meaning as in Fig. 2.</figcaption></figure>

between flow regions inside and outside the tangent cylinder. Most of the topographic changes, associated with the large scale columnar troughs and crests already seen in Fig. 4(b), indeed reside outside the tangent cylinder. This large-scale non-axisymmetric topography is correlated with a corresponding increase of the typical time $\tau_\xi$ which locally exceeds a few tenths of the viscous diffusion time.

To further analyse how $\tau_\xi$ relates to the roughness amplitude, we define

$$ \tilde{Q} = -\frac{\partial\langle\overline{T}\rangle_\varphi}{\partial r}, \quad Q' = \left[ \left\langle \left(\frac{\partial T}{\partial r}\right)^2 \right\rangle_\varphi - \left\langle\frac{\partial T}{\partial r}\right\rangle_\varphi^2 \right]^{1/2}, $$

following the decomposition already adopted in Eqs. (23)–(24). An horizontal and time average of Stefan’s condition (5) yields

$$ StPr \frac{d\xi_M}{dt} = 0 = \langle \tilde{Q} \rangle_S - \frac{T_M r_o}{(r_o - \xi_M)^\xi}, \quad (31) $$

where $\langle \tilde{Q} \rangle_S = Nu r_o / \xi_M^2$. The typical time variability of the interface roughness hence relates to the heat flux fluctuations, such that

$$ StPr \frac{d\langle \xi' \rangle_S}{dt} \approx \langle Q' \rangle_S. $$

Following Yang et al. (2023c), we make the additional hypothesis that the amplitude of the heat flux fluctuations are proportional to the average heat flux, i.e. $\langle Q' \rangle_S \sim \langle \tilde{Q} \rangle_S$. This latter assumption yields

$$ StPr \frac{d\langle \xi' \rangle_S}{dt} \sim Nu \frac{r_o}{\xi_M^2}, $$

and allows us to derive the following scaling relation for the typical time of variation of the interface roughness

$$ \langle \tau_\xi \rangle_S \sim \frac{StPr \xi_M^2 \langle \xi' \rangle_S}{Nu \xi_M^2}. \quad (32) $$

The proportionality between the correlation time and the Stefan number is expected since large $St$ implies large latent heat and therefore slower melting dynamics. This is also consistent with previous observation of melting rates being inversely proportional to the Stefan number (Favier et al., 2019). Fig. 12 shows the validity of this scaling relation, with little scatter in the data and a slope close to the expected value of one. This a posteriori validates the assumption $\langle Q' \rangle_S / \langle \tilde{Q} \rangle_S \approx \text{const.}$ retained in the derivation.

### 3.3. Roughness wavelength and amplitude

To evaluate the horizontal size of topography, Fig. 13 shows a comparison between the typical sizes of the convective flow and topography, both quantities being evaluated by the spherical harmonic degree where the corresponding spectra reach their maximum values (Eqs. (19) and (21)). The two wavenumbers are found to be broadly similar, indicating that the horizontal scale of topography follows that of the underlying convective pattern, in agreement with previous findings by e.g. Rabbanipour Esfahani et al. (2018) or Favier et al. (2019). There is less of an agreement at larger degrees (smaller lengthscales), for which several outliers exhibit topographic changes of larger length-scale than that of the convective flow. We recall that the Stefan number has been fixed to unity in all our simulations, whereas larger values could yield a smoothing of the small-scale topography variations. Given the moderate changes of $\ell_U$ across the parameter space studied here, we can however merely speculate on the possible increasing disagreement between $\ell_\xi$ and $\ell_U$ for smaller-scale convective flows.

For the configurations where the correlation time of topography exceeds the mean thermal diffusion time in the solid phase, i.e. $\tau_\xi >$
<figure><img src="image_18.png" /><figcaption>Fig. 14. Comparison between the actual time-averaged heat flux fluctuations $\partial Q$ (panel a) and the theoretical model expressed in Eq. (25) (panel b) for a simulation with $E = 10^{-3}$, $Ra = 10^7$ and $T_M = 0.6$. This is the same configuration as the one shown in Figs. 4(b) and 11. In each panel, the dashed lines mark the location of the effective tangent cylinder.</figcaption></figure>

$Pr h_S^2$, the roughness wavelength and amplitude can be related to the heat flux fluctuations using the perturbative model expressed in Eq. (25). Fig. 14 shows a comparison between the actual time-averaged heat flux variations $\partial Q$ (panel a) and the model (panel b) for the same numerical simulation already discussed in Fig. 11. Most of the heat flux variations are localised outside the tangent cylinder and locked in the columnar topography (Fig. 4b). The model – which we recall only retains the first-order contributions in terms of amplitude of topography – manages to accurately reproduce the observed heat flux variations. In this configuration, the non-axisymmetric topography promotes heat flux heterogeneities that can reach up to 20% of the average heat flux. At this stage, it is however unclear whether the locking phenomenon exemplified here will persist at stronger convective forcing (see e.g. Yang et al., 2023c) or is promoted by the large-scale convective pattern which develops at the moderate supercriticalities considered here.

Given the predominantly columnar nature of the interface roughness in this configuration, one can tentatively approximate the topographic changes by one single sectoral spherical harmonic function of degree and order $\ell_\xi$. Within this limit, the heat flux variations in the equatorial plane are approximated by

$$ \begin{aligned} \max_{\theta=\pi/2,\phi} |\partial Q| &\approx 2 \frac{|\xi_{\ell_\xi} \ell_\xi \ell_\xi|}{\xi_M} f_{\ell_\xi}(\eta_S) \max_{\theta=\pi/2,\phi} |Y_{\ell_\xi} \ell_\xi \ell_\xi|, \\ &\approx 2 \frac{|\xi_{\ell_\xi} \ell_\xi \ell_\xi|}{\xi_M} \frac{f_{\ell_\xi}(\eta_S)}{2^{\ell_\xi} \ell_\xi \ell_\xi!} \sqrt{\frac{(2\ell_\xi + 1)!}{4\pi}}. \end{aligned} \quad (33) $$

For the model with $E = 10^{-3}$, $Ra = 10^7$ and $T_M = 0.6$, one gets $\ell_\xi = 11$ and $|\xi_{\ell_\xi \ell_\xi}|/\xi_M \approx 3.4 \times 10^{-3}$ which yields $\max_{\theta=\pi/2,\phi} |\partial Q| \approx 0.1$, a value slightly underestimated compared to the actual extrema of the heat flux variations in the equatorial plane (Fig. 14a). Despite the rather crude single mode approximation involved in Eq. (33), this equation allows to directly relate the amplitude and dominant wavelength of the interface topography to the corresponding heat flux heterogeneities, which can prove useful for further order of magnitude estimates in the relevant geophysical regime.

## 4. Geophysical implications

We already noticed earlier that the ice shell is either thicker in the equatorial or in the polar regions depending on the influence of rotation (Fig. 5). Earlier studies of monophasic rotating convection by e.g. Amit et al. (2020) and Kvorka and Čadek (2022) characterised two heat flux regimes termed “polar” and “equatorial cooling” in which the heat flux is respectively larger in the polar or equatorial regions. The two studies however disagree on the control parameter which governs the transition between these two regimes. Amit et al. (2020) advocate that polar cooling occurs when $Ra E^{8/5} > 10$, while Kvorka and Čadek (2022) obtain this regime whenever $1 < Ra E^{12/7} < 10$. This latter combination of exponents on $Ra$ and $E$ is obtained by assuming a transition between the heat transfer of rapidly-rotating convection $Nu \sim Ra E^{3/2}$ (e.g. Julien et al., 2012) and its non-rotating RBC counterpart $Nu_{NR} \sim Ra^{1/3}$ (see Gastine et al., 2016). The differences

<figure><img src="image_19.png" /><figcaption>Fig. 15. Relative contrasts of ice thickness ($h_S^{h/1}$, Eq. (35)) and heat fluxes ($Q^{h/1}$, Eq. (34)) inside and outside the effective tangent cylinder as a function of $Ra_{eff} E^{3/2}$. The relative changes of the size of the solidus $h_S^{h/1}$ (high-opacity symbols) have been measured in the simulations with a phase change while the heat flux contrasts $Q^{h/1}$ (low-opacity symbols) come from the equivalent purely convective models. The symbols carry the same meaning as in Fig. 2. For comparison purposes, the measures of $Q^{h/1}$ coming from the monophasic simulations from Kvorka and Čadek (2022) with rigid boundaries (here KC2022) and from Hartmann et al. (2024) with $r_i/r_\phi = 0.8$ (here HSLV24) have been included. To account for the difference in the definition of the gravity profile, the Rayleigh numbers provided by Kvorka and Čadek (2022) have been rescaled by the average of their gravity profile over ours. For Enceladus, estimates of $h_S^{h/1}$ come from Čadek et al. (2019), while the uncertainties on $E$ and $Ra$ come from Soderlund (2019).</figcaption></figure>

in the scaling exponents between the two studies could arise from the different adopted mechanical boundary conditions.

Given the different regionalised dynamics in spherical shell convection (Wang et al., 2021; Gastine and Aurnou, 2023), we find however more appropriate to define this transition in terms of the scaling behaviour of polar convection. In the limit of rapid rotation and moderate supercriticalities, the heat transfer in the polar regions of spherical shells with rigid boundaries closely follows the scaling behaviour obtained in planar simulations, i.e. $Nu \sim Ra^3 E^4$ (King et al., 2012; Stellmach et al., 2014; Gastine and Aurnou, 2023). Assuming that the transition between equatorial and polar cooling is defined by the crossing between this scaling law and $Nu_{NR}$ then yields a parameter combination $Ra E^{3/2}$. This is in line with the recent study by Hartmann et al. (2024) who report an enhancement of the polar heat transport in spherical shell convection whenever $Ra E^{3/2} > 1$ using a set of numerical simulations with a fixed Rayleigh number of $Ra = 10^6$ and various Ekman numbers and radius ratios.

In terms of diagnostics, both regimes are usually characterised by defining the relative differences of the axisymmetric heat flux inside
and outside the tangent cylinder (e.g. Amit et al., 2020; Kvorka and Čadek, 2022; Bire et al., 2022):

$$Q^{h/l} = \frac{[\tilde{Q}]_0^{\theta_{TC}} - [\tilde{Q}]_{\theta_{TC}}^{\pi^{-\theta_{TC}}}}{[\tilde{Q}]_0^{\theta_{TC}} + [\tilde{Q}]_{\theta_{TC}}^{\pi^{-\theta_{TC}}}} \quad (34)$$

where $\theta_{TC} = \arcsin \eta_{eff}$ in the phase field models and $\theta_{TC} = \arcsin \eta$ without phase change. In the above expression, the square brackets have been employed to define angular averages between two colatitudes

$$[f]_{\theta_1}^{\theta_2} = \frac{1}{S_{\theta_1}^{\theta_2}} \int_{\theta_1}^{\theta_2} f(\theta) \sin \theta \, d\theta, \quad S_{\theta_1}^{\theta_2} = \int_{\theta_1}^{\theta_2} \sin \theta \, d\theta.$$

Here we accordingly define the relative differences in the axisymmetric ice thickness:

$$h_S^{h/l} = - \frac{[\tilde{\xi}]_{\theta_{TC}}^{\pi^{-\theta_{TC}}} - [\tilde{\xi}]_0^{\theta_{TC}}}{2r_o - [\tilde{\xi}]_0^{\theta_{TC}} - [\tilde{\xi}]_{\theta_{TC}}^{\pi^{-\theta_{TC}}}} \quad (35)$$

where the minus sign has been introduced to ensure that the variations of $h_S^{h/l}$ carry the same sign as those in $Q^{h/l}$.

Fig. 15 shows $h_S^{h/l}$ for the simulations with a phase change (coloured high-opacity symbols) and $Q^{h/l}$ for the corresponding models without phase change (coloured low-opacity symbols) as a function of $Ra_{eff} E_{eff}^{3/2}$. For comparison purposes, the simulations from Kvorka and Čadek (2022) with rigid boundaries and from Hartmann et al. (2024) with $r_i/r_o = 0.8$ have been included. Both diagnostics follow similar trends. When the influence of rotation is the strongest, i.e. $Ra_{eff} E_{eff}^{3/2} \lesssim 10$, the heat flux is larger outside the tangent cylinder and the ice is thinner there. For a limited range of parameters, $10 \lesssim Ra_{eff} E_{eff}^{3/2} \lesssim 100$, the heat flux is larger in the polar regions and ice is accordingly thicker outside the tangent cylinder. For larger values, both diagnostics taper off as the influence of rotation diminishes and the flow gradually loses its preferred axis. The limited number of simulations as well as the remaining scatter however prevent us to ascertain that $Ra_{eff} E_{eff}^{3/2}$ is the right parameter combination to describe the transition. In addition, recent simulations by Song et al. (2024) show that the heat transfer scaling rather adheres to $Nu \sim Ra^{3/2} E^2 Pr^{-1/2}$ in Cartesian geometry with rigid boundaries for $E < 3 \times 10^8$ and $Ra E^{4/3} > 10$. This is in line with the findings by Stellmach et al. (2014) who also report the same diffusivity-free scaling of rotating convection whenever stress-free boundary conditions are employed. In the geophysical regime relevant to the icy satellites, this is hence plausible that the transition between polar and equatorial cooling is rather governed by a $Ra E^{12/7}$ parameter combination.

With all these possible caveats in mind, we nevertheless tentatively locate in Fig. 15 the expected values for Enceladus for which ice shell thickness models have been devised. Using topography and gravity data from Cassini, Čadek et al. (2019) for instance derived a model of Enceladus' ice shell thickness which is about 30 km thick at the equator and reaches 15 (5) km near the North (South) pole. This would place this subsurface ocean in the polar cooling regime with $h_S^{h/l} \in [0.13, 0.20]$. Due to the uncertainties on the ocean thickness, the estimated values for E and Ra respectively span the intervals $[10^{-10}, 10^{-11}]$ and $[10^{16}, 10^{19}]$ (e.g. Soderlund, 2019; Cabanes et al., 2024). These estimates yield $Ra E^{3/2} > 10$ a parameter combination indeed compatible with the polar cooling regime obtained in numerical models. The amplitude of $h_S^{h/l}$ for Enceladus however exceeds the largest contrast obtained in our simulations as well as – to a lesser extent – the largest heat flux contrasts $Q^{h/l}$ obtained by Kvorka and Čadek (2022) and Hartmann et al. (2024) for their $r_i/r_o = 0.8$ simulations. Simulations by Bire et al. (2022) conducted at lower Ekman numbers also suggest a gradual drop of $Q^{h/l}$ when the radius ratio of the ocean increases (see their Fig. 12). At this stage, this is hence unclear whether the intrinsic heat flux fluctuations of rotating convection are susceptible to reach a sufficient amplitude at the geophysical parameters to explain the topographic changes of Enceladus' ice. The strong North-South asymmetry could also originate from large scale basal

heat flux heterogeneities or from temperature contrasts at Enceladus' surface (e.g. Lemasquerier et al., 2023b). Ice shell models of Titan are more uncertain but also suggestive of a polar cooling regime, with for instance $h_S^{h/l} \sim 0.05$ in the model by Lefevre et al. (2014). Due to the least rotational constraint on Titan's subsurface oceanic flows, estimates of Ra and E from Soderlund (2019) yield $Ra E^{3/2} \sim \mathcal{O}(10^4)$. Because of the scatter of the numerical simulations beyond $Ra E^{3/2} \sim 100$ (Fig. 15), it is uncertain whether such large scale topographic changes could be attributed to the convective fluctuations or are rather promoted by large scale thermal heterogeneities at the base or at the top of the fluid layer. At this stage, it is also important to recall that several physical ingredients such as salinity and thermobaric effects have been neglected in the current model and are likely to change the interplay between heat flux changes and mean topography (e.g. Kang, 2023).

Despite all the shortcomings of our models, the observation of long-lived stable non-axisymmetric corrugations of the solid–liquid interface in several simulations prompts us to attempt a re-scaling to the relevant geophysical regime. To do so, we make the following bold assumptions: (i) the non-axisymmetric topography can be described by one single sectoral model of size $\ell_U$; (ii) the ice layer is in a quasi-equilibrated conducting state that allows the application of the model Eq. (25) which relates the heat flux variations to the topographic wavelength and amplitude. Given the estimated average radius ratio of the ice layers of the icy satellites, the single mode approximation (Eq. (33)) can be further simplified to

$$\partial \xi \sim \frac{\partial Q}{\ell_U} \sim \frac{\partial Q \mathcal{L}_U}{\pi \xi_M} \quad (36)$$

where $\mathcal{L}_U$ is the dominant convective flow lengthscale. The studies of the transfer functions of basal heterogeneous heat fluxes through a rotating convective layer carried out by Terra-Nova et al. (2023) and Lemasquerier et al. (2023b) suggest $\partial Q \sim 2-4$ for Enceladus (see also Čadek et al., 2019, for similar estimates) and $\partial Q \sim 0.1-0.5$ for Titan (Choblet et al., 2017). To provide an estimate of the lengthscale $\mathcal{L}_U$, we further assume that the oceanic flows are rotationally-constrained and that the kinetic energy spectra of the zonal jets and the non-axisymmetric residuals reach their maxima at similar scales, in practice found by previous studies to be close to the Rhines scale (e.g Guervilly et al., 2019; Lemasquerier et al., 2023a). We therefore use the jet sizes derived by Cabanes et al. (2024) to estimate $\mathcal{L}_U$, such that $\mathcal{L}_U \sim 0.1D \sim 2$ km for Enceladus and $\mathcal{L}_U \sim 0.5D \sim 100$ km for Titan, with D the ocean depth. For both icy satellites, this yields $\partial \xi \sim \mathcal{O}(10^{-3}-10^{-2})$, depending on the amplitude of the heat flux variations. This translates to non-axisymmetric topographic changes at the solid–liquid interface that could reach up to $\mathcal{O}(10^2-10^3)$ meters for Enceladus and $\mathcal{O}(10^3-10^4)$ meters for Titan. Considering 2-D thermo-mechanical models of the ice response to basal heat flux variations, Kihoulou et al. (2023) demonstrated that the topographic changes in a conducting ice shell strongly depend on the ice viscosity as well as on the variations of the melting temperature with pressure (see their Fig. 6). Despite these uncertainties, their estimates of topographic changes are lower than ours, typically covering the range of $\mathcal{O}(10^2-10^3)$ m at the ice-water interface (see also Kang, 2023). Again, accounting for double-diffusive effects, ice creeping and temperature variations due to pressure changes along the deformed boundary are susceptible to change these conclusions.

## 5. Conclusion

Improving understanding of the dynamical coupling between the ocean and the overlying ice layer is of primary importance to better characterise the hydrosphere of the icy moons. In this study we have developed a novel approach to model rotating convection with a phase change using a phase field formulation (e.g. Beckermann et al., 1999).

To examine the interplay between rotating convection and a melting boundary, we have conducted a series of numerical simulations
in spherical geometry varying the control parameters – the Rayleigh number Ra and the Ekman number E – as well as the melting temperature. We have split the analysis between the large-scale axisymmetric topography and the non-axisymmetric features. For the former and in line with previous monophasic convection models, we have evidenced two regimes in which the mean axisymmetric heat flux either peaks at the equator or at the poles with a corresponding regionalised ice thinning. Transition between the two happens when RaE³/² ≈ 10, a parameter combination which is obtained when the scaling behaviour of rotating convection in the polar regions of the spherical shell crosses the one of classical RBC (Hartmann et al., 2024).

We have conducted numerical simulations of the equivalent monophasic setups and we have derived a perturbative model of the heat equilibrium in a quasi-spherical shell with bottom topographic changes. This enabled us to show that the relative heat flux variations at the top of the monophasic simulations usually provide a good guess of the actual flux as long as the mean axisymmetric topographic changes are small. As expected, departures between the configuration with a phase change and their counterparts without become more pronounced for large topographic variations.

Non-axisymmetric topography, termed roughness in this study, yields significant differences with the usual rotating convective flow in spherical shells. In particular, we have evidenced the formation of long-lived large-scale columnar troughs and crests which develop for intermediate rotational constraints. Convective upwellings are then locked in the topographic changes of the solid–liquid interface. Although our parameter coverage does not allow to determine the parameter combination which governs the transition to locked-in convection (see also Yang et al., 2023c), a tentative rescaling of the amplitude of this non-axisymmetric topography to the planetary regime yields $O(10^2 - 10^3)$ meters for Enceladus and $O(10^3 - 10^4)$ meters for Titan. Those values are substantially larger than current estimates by e.g. Corlies et al. (2017) or Kang (2023). These differences could possibly arise from several simplifications of our model. We do not include the effect of salinity (Ashkenazy and Tziperman, 2021; Wong et al., 2022), the pressure dependence of the melting temperature (Labrosse et al., 2018; Lawrence et al., 2024), the creep properties of ice (e.g. Weertman, 1983; Shibley and Goodman, 2024) the effect of tidal heating in the ice crust (Beuthe, 2019; Běhounková et al., 2021), or the thermal heterogeneities at the base of the ocean (Terra-Nova et al., 2023; Lemasquerier et al., 2023b) or at the moon’s surface (Weller et al., 2019). All those endogenic physical processes are possible candidates susceptible to raise the complexity of the results described here. Accounting for those effects could be the subject of future studies using the phase change formalism discussed here.

On a longer term, improving our understanding of the interplay between oceanic flows and phase changes in the icy moons will also require to conduct numerical simulations at more extreme parameters. Though numerically challenging (e.g. Song et al., 2024), this would allow to consolidate the scaling relation which governs the transition between equatorial and polar cooling and to study the coupling between mean zonal flows and topography (e.g. Hay et al., 2023).

### CRediT authorship contribution statement

**T. Gastine:** Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Investigation, Conceptualization.
**B. Favier:** Writing – review & editing, Validation, Methodology, Investigation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this article.

### Acknowledgements

The authors thank two anonymous reviewers for their suggestions that helped improve the manuscript. The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), France, under grant ANR-20-CE49-0010 (project COLOSSe). Numerical computations have been carried out on the S-CAPAD/DANTE platform at IPGP. All the figures have been generated using matplotlib (Hunter, 2007) and paraview (https://www.paraview.org). Several colourmaps come from the cmocean package by Thyng et al. (2016).

## Appendix A. Mapping functions for Chebyshev collocation method

When the fluid domain features regions of rapid changes, considering mappings can significantly improve the convergence of the Chebyshev pseudo-spectral approximation. The mapping introduced by Kosloff and Tal-Ezer (1993) aims at reducing the grid points clustering near the boundaries inherent in the native Gauss–Lobatto collocation grid points. The amplitude of the grid stretching is governed by a control parameter $α_1 ∈ [0, 1]$. The mapping function is expressed by

$$F(x) = \frac{\arcsin(\alpha_1 x)}{\arcsin \alpha_1} \qquad (37)$$

The main purpose of the mapping defined by Bayliss and Turkel (1992) is to refine the grid spacing around a particular interior point $x = α_2$. It is governed by two input parameters: $α_1$ which controls the stiffness of the grid refinement and $α_2 ∈ [-1, 1]$ which defines the centre of the mapping function. It is defined by

$$F(x) = α_2 + \frac{1}{α_1} \tan[\lambda(x - x_0)], \qquad (38)$$

with

$$λ = \frac{\arctan[α_1(1 - α_2)]}{1 - x_0}, \quad x_0 = \frac{\arctan[α_1(1 + α_2)] - 1}{\arctan[α_1(1 - α_2)] + 1}$$

Similarly to the mapping by Bayliss and Turkel (1992), the mapping by Jafari-Varzaneh and Hosseini (2015) was also introduced to handle steep localised fronts. It is governed by three input parameters: $α_1$ and $α_2$ retaining the same meaning as for the previous mapping and $α_3 ∈ [0.2, 0.9]$ being a small parameter that we keep to a fixed value of $α_3 = 0.4$. It is defined by

$$F(x) = α_2 + \frac{1}{α_1} \sinh \left\{ A \left[ \frac{\tan(x \arctan C)}{C} - 1 \right] + B \right\} \qquad (39)$$

with

$$A = \frac{1}{2} \{ \operatorname{argsinh}[α_1(1 - α_2)] + \operatorname{argsinh}[α_1(1 + α_2)] \},$$

$$B = \operatorname{argsinh}[α_1(1 - α_2)],$$

$$C = \left[ \left\{ \Im \left\{ \frac{1}{A} \left( \frac{i\pi}{2} - B \right) + 1 \right\} \right\} + α_3 \right]^{-1}.$$

## Appendix B. A benchmark for rotating convection in a spherical shell with a phase change

To validate the numerical implementation of the phase field method in the pseudo-spectral code MagIC, we consider a weakly nonlinear configuration of rotating convection close to onset. Similarly to classical benchmarks of monophasic convection in spherical geometry (e.g. Christensen et al., 2001; Marti et al., 2014), this allows to reach a saturated state which takes the form of steadily drifting convective columns. This quasi-stationarity allows to benchmark well-defined integrated diagnostics such as the total kinetic energy $E_K$ or the Nusselt number $Nu$ defined by (15) and (17) respectively. Here we adopt $E = 10^{-3}$, $Ra = 1.8 \times 10^5$, $Pr = 1$, $St = 1$, $T_M = 0.25$ for a radius ratio $η = 0.35$. For this parameter combination, the solid–liquid interface equilibrates
<figure><img src="image_20.png" /><figcaption>Fig. 16. Snapshots of the benchmark configuration with $E = 10^{-3}$, $Ra = 1.8 \times 10^{5}$, $Pr = 1$, $St = 1$, $T_M = 0.25$ and $r_i/r_o = 0.35$ for a numerical simulation with $\epsilon = 10^{-3}$. (a) Radial velocity $u_r$ in the equatorial plane. (b) Hammer projection of the melt radius $r_M$. The dashed line in panel (a) marks the location of the solid-liquid interface $r_M(\theta = \pi/2, \phi)$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 17. (a) Comparison of the radial profiles of kinetic energy for the benchmark configuration with decreasing values of the Cahn number $\epsilon$. (b) Comparison of the azimuthal average of the melt radius $\tilde{\xi}$ for decreasing values of $\epsilon$. (c) $L^2$ relative errors (Eq. (40)) of several diagnostic quantities as a function of $\epsilon$.</figcaption></figure>

around the mean radius $\xi_M \approx 1.09$, which yields the following effective quantities: $E_{eff} \approx 3.4 \times 10^{-3}$, $Ra_{eff} \approx 1.5 \times 10^4$ and $\eta_{eff} \approx 0.5$. For this Ekman number and radius ratio, the first unstable mode that becomes linearly unstable in the equivalent monophasic convection problem features an azimuthal wavenumber $m = 5$ and a critical Rayleigh number $Ra_c = 1.322 \times 10^4$ (Barik et al., 2023). This is closely followed by the $m = 6$ mode which has $Ra_c = 1.341 \times 10^4$. In order to define a reproducible reference case, we hence initiate all the numerical models using a thermal perturbation with a fivefold azimuthal symmetry.

Fig. 16 shows a snapshot of the radial velocity (panel a) and the solid-liquid interface (panel b) for a numerical simulation with $\epsilon = 10^{-3}$. The solution takes the form of convective columns confined in the lower half of the fluid domain. The interface features a fivefold symmetry corrugation outside the effective tangent cylinder. Due to the localised convective columns, the ice thickness is thinner near the equator than the poles.

Fig. 17 illustrates the convergence of the phase field method for this benchmark case when decreasing the Cahn number $\epsilon$. We recall that the phase-field formulation converges towards the original Stefan problem as $\epsilon$ vanishes. Panel (a) shows radial profiles of the kinetic energy, while panel (b) shows the latitudinal profile of the axisymmetric melt radius $\tilde{\xi}$. On decreasing $\epsilon$, the solid-liquid transition steepens, which goes along with a gradual increase of the kinetic energy content and an outward shift of the interface. Both diagnostics are suggestive of a gradual convergence towards an asymptotic solution as $\epsilon \to 0$. To examine this convergence in a more quantitative way, we compute in Fig. 17(c) the $L^2$ relative error between a given phase field model and the reference solution, here defined by the simulation with the lowest

value of $\epsilon = 7 \times 10^{-4}$. We consider several diagnostics, either based on mean radial profiles or mean latitudinal profiles:

$$ \|f\|_{L^2, r} = \left( \frac{\int_{r_i}^{r_o} [f(r) - f_{\text{ref}}(r)]^2 r^2 dr}{\int_{r_i}^{r_o} f_{\text{ref}}^2 dr} \right)^{1/2}, \quad (40) $$

$$ \|f\|_{L^2, \theta} = \left( \frac{\int_{0}^{\pi} [f(\theta) - f_{\text{ref}}(\theta)]^2 \sin \theta d\theta}{\int_{0}^{\pi} f_{\text{ref}}^2 \sin \theta d\theta} \right)^{1/2}, $$

where $f_{\text{ref}}$ in the above equations correspond to the reference solution. All diagnostics exhibit a similar trend with a slow convergence close to first order for the largest values of $\epsilon$ and a gradual increase towards the second order for $\epsilon \approx 10^{-3}$. We note that a similar transition from first to second order convergence has been observed in Favier et al. (2019). The convergence of our phase field model is therefore more complex than the second order reported by Hester et al. (2020). In addition to the known influence of the penalty coefficient $\tau_p$ on the convergence rate (Hester et al., 2021b), it is very likely that the influence of rotation, which was not considered in Hester et al. (2020), modifies the convergence behaviour of the phase field formulation.

To ease the future comparison of phase field models in spherical geometry, Table 1 lists the values of several integrated diagnostics for the benchmark configuration ($\xi_M$, $E_K$ and $Nu$) along with the temporal and spatial resolution employed to ensure an appropriate convergence of the models.

## Appendix C. Table of results

Lists the input dimensionless numbers, the main global diagnostics as well as the numerical parameters of all the direct numerical simulations considered in this study. All simulations have been computed using $\eta = r_i/r_o = 0.8$, $Pr = 1$ and $St = 1$.
<caption>Table 1</caption>

Table of results of the benchmark configuration. All the simulations have been conducted using the BT mapping with $\alpha_1 = 10$, the third order SBDF3 time scheme and a penalty coefficient $\tau_p = 0.6$.

<table><caption>Table 2</caption><thead><tr><th>ε</th><th>ξ<sub>M</sub></th><th>E<sub>K</sub></th><th>Nu</th><th>δt</th><th>(N<sub>R</sub>, ε<sub>max</sub>)</th></tr></thead><tbody><tr><td>7 × 10<sup>-3</sup></td><td>1.063</td><td>11.93</td><td>1.112</td><td>8 × 10<sup>-6</sup></td><td>(161, 213)</td></tr><tr><td>5 × 10<sup>-3</sup></td><td>1.068</td><td>13.04</td><td>1.111</td><td>4 × 10<sup>-6</sup></td><td>(161, 213)</td></tr><tr><td>4 × 10<sup>-3</sup></td><td>1.072</td><td>14.25</td><td>1.116</td><td>10<sup>-6</sup></td><td>(193, 213)</td></tr><tr><td>3 × 10<sup>-3</sup></td><td>1.076</td><td>15.86</td><td>1.123</td><td>8 × 10<sup>-7</sup></td><td>(257, 213)</td></tr><tr><td>2.5 × 10<sup>-3</sup></td><td>1.078</td><td>16.81</td><td>1.127</td><td>4 × 10<sup>-7</sup></td><td>(257, 213)</td></tr><tr><td>2 × 10<sup>-3</sup></td><td>1.081</td><td>17.82</td><td>1.132</td><td>3 × 10<sup>-7</sup></td><td>(257, 213)</td></tr><tr><td>1.5 × 10<sup>-3</sup></td><td>1.083</td><td>18.94</td><td>1.138</td><td>2 × 10<sup>-7</sup></td><td>(321, 320)</td></tr><tr><td>10<sup>-3</sup></td><td>1.085</td><td>20.09</td><td>1.144</td><td>10<sup>-7</sup></td><td>(385, 426)</td></tr><tr><td>7 × 10<sup>-4</sup></td><td>1.087</td><td>20.77</td><td>1.147</td><td>10<sup>-7</sup></td><td>(513, 533)</td></tr></tbody></table>

Table of results. All the simulations have been computed with a phase field parameter $a = 1$.

<table><thead><tr><th>T<sub>M</sub></th><th>ξ<sub>M</sub></th><th>E<sub>eff</sub></th><th>Ra<sub>eff</sub></th><th>Re<sub>ε</sub></th><th>Nu<sub>eff</sub></th><th>ℓ<sub>U</sub></th><th>ℓ<sub>ε</sub></th><th>ε</th><th>τ<sub>p</sub></th><th>(N<sub>R</sub>, ε<sub>max</sub>)</th><th>(Map, α<sub>1</sub>, α<sub>2</sub>)</th></tr></thead><tbody><tr><td colspan="12" style="text-align: left;">E = 10<sup>-3</sup>, Ra = 3 × 10<sup>6</sup></td></tr><tr><td>0.00</td><td>–</td><td>10<sup>-3</sup></td><td>3 × 10<sup>6</sup></td><td>297.3</td><td>9.54</td><td>12</td><td>–</td><td>–</td><td>–</td><td>(65, 426)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.44</td><td>4.87</td><td>1.325 × 10<sup>-3</sup></td><td>1.073 × 10<sup>6</sup></td><td>170.9</td><td>6.66</td><td>14</td><td>13</td><td>2 × 10<sup>-3</sup></td><td>0.40</td><td>(193, 426)</td><td>(JVH, 20, 0.73)</td></tr><tr><td>–</td><td>–</td><td>1.324 × 10<sup>-3</sup></td><td>1.072 × 10<sup>6</sup></td><td>174.8</td><td>6.92</td><td>14</td><td>–</td><td>–</td><td>–</td><td>(65, 426)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.49</td><td>4.83</td><td>1.455 × 10<sup>-3</sup></td><td>8.416 × 10<sup>5</sup></td><td>147.0</td><td>5.94</td><td>16</td><td>15</td><td>3 × 10<sup>-3</sup></td><td>0.20</td><td>(193, 341)</td><td>(JVH, 5, 0.64)</td></tr><tr><td>–</td><td>–</td><td>1.455 × 10<sup>-3</sup></td><td>8.416 × 10<sup>5</sup></td><td>154.6</td><td>6.44</td><td>16</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.69</td><td>4.58</td><td>2.972 × 10<sup>-3</sup></td><td>1.663 × 10<sup>5</sup></td><td>66.9</td><td>3.88</td><td>18</td><td>18</td><td>3 × 10<sup>-3</sup></td><td>0.40</td><td>(257, 213)</td><td>(KTE, 0.99, –)</td></tr><tr><td>–</td><td>–</td><td>2.972 × 10<sup>-3</sup></td><td>1.663 × 10<sup>5</sup></td><td>68.4</td><td>4.03</td><td>21</td><td>–</td><td>–</td><td>–</td><td>(65, 213)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td colspan="12" style="text-align: left;">E = 10<sup>-3</sup>, Ra = 10<sup>7</sup></td></tr><tr><td>0.00</td><td>–</td><td>10<sup>-3</sup></td><td>10<sup>7</sup></td><td>559.9</td><td>14.39</td><td>9</td><td>–</td><td>–</td><td>–</td><td>(97, 426)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.46</td><td>4.90</td><td>1.241 × 10<sup>-3</sup></td><td>3.825 × 10<sup>6</sup></td><td>324.0</td><td>9.48</td><td>12</td><td>11</td><td>5 × 10<sup>-3</sup></td><td>0.20</td><td>(97, 213)</td><td>(JVH, 20, 0.79)</td></tr><tr><td>–</td><td>–</td><td>1.241 × 10<sup>-3</sup></td><td>3.810 × 10<sup>6</sup></td><td>350.9</td><td>10.80</td><td>12</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.60</td><td>4.82</td><td>1.502 × 10<sup>-3</sup></td><td>2.094 × 10<sup>6</sup></td><td>254.5</td><td>8.48</td><td>12</td><td>11</td><td>3 × 10<sup>-3</sup></td><td>0.40</td><td>(129, 394)</td><td>(JVH, 10, 0.62)</td></tr><tr><td>–</td><td>–</td><td>1.506 × 10<sup>-3</sup></td><td>2.084 × 10<sup>6</sup></td><td>263.2</td><td>9.14</td><td>13</td><td>–</td><td>–</td><td>–</td><td>(65, 256)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.69</td><td>4.71</td><td>1.993 × 10<sup>-3</sup></td><td>1.037 × 10<sup>6</sup></td><td>181.0</td><td>6.90</td><td>14</td><td>13</td><td>3 × 10<sup>-3</sup></td><td>0.40</td><td>(193, 341)</td><td>(JVH, 7, 0.43)</td></tr><tr><td>–</td><td>–</td><td>1.993 × 10<sup>-3</sup></td><td>1.037 × 10<sup>6</sup></td><td>187.9</td><td>7.44</td><td>13</td><td>–</td><td>–</td><td>–</td><td>(65, 256)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td colspan="12" style="text-align: left;">E = 10<sup>-3</sup>, Ra = 3 × 10<sup>7</sup></td></tr><tr><td>0.00</td><td>–</td><td>10<sup>-3</sup></td><td>3 × 10<sup>7</sup></td><td>998.4</td><td>20.59</td><td>7</td><td>–</td><td>–</td><td>–</td><td>(97, 512)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.40</td><td>4.95</td><td>1.117 × 10<sup>-3</sup></td><td>1.508 × 10<sup>7</sup></td><td>673.9</td><td>14.96</td><td>8</td><td>7</td><td>3 × 10<sup>-3</sup></td><td>0.25</td><td>(257, 341)</td><td>(JVH, 20, 0.894)</td></tr><tr><td>–</td><td>–</td><td>1.117 × 10<sup>-3</sup></td><td>1.508 × 10<sup>7</sup></td><td>725.9</td><td>17.46</td><td>8</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.45</td><td>4.93</td><td>1.150 × 10<sup>-3</sup></td><td>1.321 × 10<sup>7</sup></td><td>630.0</td><td>14.34</td><td>8</td><td>8</td><td>3 × 10<sup>-3</sup></td><td>0.25</td><td>(257, 341)</td><td>(JVH, 10, 0.79)</td></tr><tr><td>–</td><td>–</td><td>1.150 × 10<sup>-3</sup></td><td>1.321 × 10<sup>7</sup></td><td>678.6</td><td>16.70</td><td>8</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.50</td><td>4.92</td><td>1.194 × 10<sup>-3</sup></td><td>1.130 × 10<sup>7</sup></td><td>584.4</td><td>13.64</td><td>9</td><td>9</td><td>3 × 10<sup>-3</sup></td><td>0.25</td><td>(257, 512)</td><td>(JVH, 20, 0.78)</td></tr><tr><td>–</td><td>–</td><td>1.195 × 10<sup>-3</sup></td><td>1.129 × 10<sup>7</sup></td><td>626.9</td><td>15.84</td><td>8</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.55</td><td>4.89</td><td>1.255 × 10<sup>-3</sup></td><td>9.400 × 10<sup>6</sup></td><td>533.6</td><td>12.88</td><td>9</td><td>10</td><td>3 × 10<sup>-3</sup></td><td>0.25</td><td>(257, 341)</td><td>(JVH, 10, 0.7)</td></tr><tr><td>–</td><td>–</td><td>1.255 × 10<sup>-3</sup></td><td>9.399 × 10<sup>6</sup></td><td>571.9</td><td>14.89</td><td>8</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.60</td><td>4.86</td><td>1.338 × 10<sup>-3</sup></td><td>7.543 × 10<sup>6</sup></td><td>478.1</td><td>12.04</td><td>10</td><td>10</td><td>3 × 10<sup>-3</sup></td><td>0.25</td><td>(257, 341)</td><td>(JVH, 10, 0.62)</td></tr><tr><td>–</td><td>–</td><td>1.338 × 10<sup>-3</sup></td><td>7.543 × 10<sup>6</sup></td><td>513.2</td><td>13.88</td><td>9</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.65</td><td>4.83</td><td>1.469 × 10<sup>-3</sup></td><td>5.694 × 10<sup>6</sup></td><td>417.1</td><td>11.06</td><td>10</td><td>11</td><td>3 × 10<sup>-3</sup></td><td>0.25</td><td>(257, 341)</td><td>(JVH, 10, 0.55)</td></tr><tr><td>–</td><td>–</td><td>1.469 × 10<sup>-3</sup></td><td>5.694 × 10<sup>6</sup></td><td>447.3</td><td>12.73</td><td>10</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.69</td><td>4.79</td><td>1.604 × 10<sup>-3</sup></td><td>4.384 × 10<sup>6</sup></td><td>371.2</td><td>10.63</td><td>11</td><td>14</td><td>3 × 10<sup>-3</sup></td><td>0.40</td><td>(193, 426)</td><td>(JVH, 5, 0.46)</td></tr><tr><td>–</td><td>–</td><td>1.604 × 10<sup>-3</sup></td><td>4.384 × 10<sup>6</sup></td><td>391.0</td><td>11.54</td><td>11</td><td>–</td><td>–</td><td>–</td><td>(65, 512)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.75</td><td>4.70</td><td>2.056 × 10<sup>-3</sup></td><td>2.391 × 10<sup>6</sup></td><td>277.2</td><td>8.81</td><td>14</td><td>13</td><td>3 × 10<sup>-3</sup></td><td>0.30</td><td>(193, 341)</td><td>(JVH, 5, 0.23)</td></tr><tr><td>–</td><td>–</td><td>2.056 × 10<sup>-3</sup></td><td>2.391 × 10<sup>6</sup></td><td>294.1</td><td>9.82</td><td>12</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.80</td><td>4.57</td><td>3.041 × 10<sup>-3</sup></td><td>1.035 × 10<sup>6</sup></td><td>187.3</td><td>6.86</td><td>14</td><td>13</td><td>3 × 10<sup>-3</sup></td><td>0.25</td><td>(257, 341)</td><td>(JVH, 5, 0.04)</td></tr><tr><td>–</td><td>–</td><td>3.041 × 10<sup>-3</sup></td><td>1.035 × 10<sup>6</sup></td><td>197.8</td><td>7.76</td><td>13</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td colspan="12" style="text-align: left;">E = 3 × 10<sup>-4</sup>, Ra = 1.2 × 10<sup>7</sup></td></tr><tr><td>0.00</td><td>–</td><td>3 × 10<sup>-4</sup></td><td>1.2 × 10<sup>7</sup></td><td>517.6</td><td>13.39</td><td>20</td><td>–</td><td>–</td><td>–</td><td>(97, 426)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.44</td><td>4.90</td><td>3.681 × 10<sup>-4</sup></td><td>4.848 × 10<sup>6</sup></td><td>308.5</td><td>9.52</td><td>26</td><td>26</td><td>10<sup>-3</sup></td><td>1.51</td><td>(193, 512)</td><td>(JVH, 100, 0.837)</td></tr><tr><td>–</td><td>–</td><td>3.681 × 10<sup>-4</sup></td><td>4.848 × 10<sup>6</sup></td><td>322.2</td><td>10.08</td><td>25</td><td>–</td><td>–</td><td>–</td><td>(65, 293)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.49</td><td>4.88</td><td>3.908 × 10<sup>-4</sup></td><td>4.015 × 10<sup>6</sup></td><td>275.2</td><td>8.78</td><td>27</td><td>26</td><td>3 × 10<sup>-3</sup></td><td>0.30</td><td>(193, 341)</td><td>(JVH, 50, 0.796)</td></tr><tr><td>–</td><td>–</td><td>3.908 × 10<sup>-4</sup></td><td>4.015 × 10<sup>6</sup></td><td>292.8</td><td>9.50</td><td>26</td><td>–</td><td>–</td><td>–</td><td>(65, 293)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.54</td><td>4.85</td><td>4.187 × 10<sup>-4</sup></td><td>3.245 × 10<sup>6</sup></td><td>246.3</td><td>8.25</td><td>28</td><td>27</td><td>3 × 10<sup>-3</sup></td><td>0.30</td><td>(129, 341)</td><td>(JVH, 50, 0.72)</td></tr><tr><td>–</td><td>–</td><td>4.181 × 10<sup>-4</sup></td><td>3.243 × 10<sup>6</sup></td><td>259.9</td><td>8.82</td><td>28</td><td>–</td><td>–</td><td>–</td><td>(97, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.69</td><td>4.69</td><td>6.244 × 10<sup>-4</sup></td><td>1.163 × 10<sup>6</sup></td><td>150.6</td><td>6.34</td><td>32</td><td>28</td><td>2 × 10<sup>-3</sup></td><td>0.60</td><td>(257, 426)</td><td>(BT, 10, 0.41)</td></tr><tr><td>–</td><td>–</td><td>6.244 × 10<sup>-4</sup></td><td>1.163 × 10<sup>6</sup></td><td>154.4</td><td>6.41</td><td>32</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.75</td><td>4.56</td><td>9.602 × 10<sup>-4</sup></td><td>4.777 × 10<sup>5</sup></td><td>97.3</td><td>4.79</td><td>36</td><td>14</td><td>2 × 10<sup>-3</sup></td><td>0.60</td><td>(257, 426)</td><td>(BT, 10, 0.135)</td></tr><tr><td>–</td><td>–</td><td>9.602 × 10<sup>-4</sup></td><td>4.777 × 10<sup>5</sup></td><td>98.3</td><td>4.83</td><td>38</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr><tr><td>0.80</td><td>4.35</td><td>2.511 × 10<sup>-3</sup></td><td>8.613 × 10<sup>4</sup></td><td>45.0</td><td>2.71</td><td>41</td><td>35</td><td>2 × 10<sup>-3</sup></td><td>0.30</td><td>(385, 213)</td><td>(BT, 5, –0.36)</td></tr><tr><td>–</td><td>–</td><td>2.511 × 10<sup>-3</sup></td><td>8.613 × 10<sup>4</sup></td><td>42.4</td><td>2.94</td><td>42</td><td>–</td><td>–</td><td>–</td><td>(65, 341)</td><td>(GL<sub>U</sub>, –)</td></tr></tbody></table>

(continued on next page)
<table><caption>Table 2 (continued).</caption><thead><tr><th>T<sub>M</sub></th><th>ξ<sub>M</sub></th><th>E<sub>eff</sub></th><th>Ra<sub>eff</sub></th><th>Re<sub>ε</sub></th><th>Nu<sub>eff</sub></th><th>ℓ<sub>U</sub></th><th>ℓ<sub>ξ</sub></th><th>ε</th><th>τ<sub>p</sub></th><th>(N<sub>R</sub>, ℓ<sub>max</sub>)</th><th>(Map, α<sub>1</sub>, α<sub>2</sub>)</th></tr></thead><tbody><tr><td colspan="13">E = 3 × 10<sup>-5</sup>, Ra = 2.5 × 10<sup>8</sup></td></tr><tr><td>0.00</td><td>–</td><td>3 × 10<sup>-5</sup></td><td>2.5 × 10<sup>8</sup></td><td>1968.3</td><td>28.87</td><td>28</td><td>–</td><td>–</td><td>–</td><td>(129, 512)</td><td>(GL<sub>−</sub>, −)</td></tr><tr><td>0.50</td><td>4.93</td><td>3.432 × 10<sup>-5</sup></td><td>1.008 × 10<sup>8</sup></td><td>1123.3</td><td>19.85</td><td>35</td><td>28</td><td>5 × 10<sup>-4</sup></td><td>1.20</td><td>(385, 682)</td><td>(JVH, 15, 0.82)</td></tr><tr><td>–</td><td>–</td><td>3.432 × 10<sup>-5</sup></td><td>1.008 × 10<sup>8</sup></td><td>1143.0</td><td>18.98</td><td>33</td><td>–</td><td>–</td><td>–</td><td>(129, 512)</td><td>(GL<sub>−</sub>, −)</td></tr><tr><td>0.59</td><td>4.90</td><td>3.674 × 10<sup>-5</sup></td><td>7.416 × 10<sup>7</sup></td><td>905.4</td><td>17.42</td><td>38</td><td>38</td><td>10<sup>-3</sup></td><td>0.60</td><td>(257, 512)</td><td>(JVH, 20, 0.79)</td></tr><tr><td>–</td><td>–</td><td>3.674 × 10<sup>-5</sup></td><td>7.416 × 10<sup>7</sup></td><td>952.0</td><td>16.86</td><td>35</td><td>–</td><td>–</td><td>–</td><td>(97, 512)</td><td>(GL<sub>−</sub>, −)</td></tr><tr><td>0.69</td><td>4.83</td><td>4.344 × 10<sup>-5</sup></td><td>4.297 × 10<sup>7</sup></td><td>669.8</td><td>13.98</td><td>41</td><td>30</td><td>10<sup>-3</sup></td><td>0.40</td><td>(257, 512)</td><td>(JVH, 20, 0.633)</td></tr><tr><td>–</td><td>–</td><td>4.347 × 10<sup>-5</sup></td><td>4.294 × 10<sup>7</sup></td><td>702.3</td><td>13.87</td><td>40</td><td>–</td><td>–</td><td>–</td><td>(97, 512)</td><td>(GL<sub>−</sub>, −)</td></tr><tr><td>0.75</td><td>4.74</td><td>5.495 × 10<sup>-5</sup></td><td>2.389 × 10<sup>7</sup></td><td>465.9</td><td>10.76</td><td>47</td><td>20</td><td>2 × 10<sup>-3</sup></td><td>0.15</td><td>(385, 512)</td><td>(KTE, 0.994, −)</td></tr><tr><td>–</td><td>–</td><td>5.482 × 10<sup>-5</sup></td><td>2.390 × 10<sup>7</sup></td><td>516.0</td><td>11.41</td><td>45</td><td>–</td><td>–</td><td>–</td><td>(97, 512)</td><td>(GL<sub>−</sub>, −)</td></tr><tr><td>0.80</td><td>4.60</td><td>8.388 × 10<sup>-5</sup></td><td>9.834 × 10<sup>6</sup></td><td>302.4</td><td>7.66</td><td>53</td><td>14</td><td>2 × 10<sup>-3</sup></td><td>0.13</td><td>(257, 426)</td><td>(JVH, 5, 0.24)</td></tr><tr><td>–</td><td>–</td><td>8.394 × 10<sup>-5</sup></td><td>9.824 × 10<sup>6</sup></td><td>324.6</td><td>8.40</td><td>55</td><td>–</td><td>–</td><td>–</td><td>(97, 512)</td><td>(GL<sub>−</sub>, −)</td></tr></tbody></table>

## Appendix D. Temperature diffusion in a quasi-spherical shell with bottom topographic changes

We solve for the temperature diffusion $\nabla^2 T(r, \theta, \varphi) = 0$ in a solid with a spherical upper boundary located at $r = r_o$ and a bottom boundary which features topographic changes attributed to heat flux variations coming from the underlying fluid layer. In a time-averaged sense, the latter boundary is described by the mean melt radius $\overline{r_M}(\theta, \varphi)$. We assume fixed temperature at both boundaries with $T(r = r_o, \theta, \varphi) = 0$ and $T[r = r_M(\theta, \varphi), \theta, \varphi] = T_M$. Let us introduce the following spherical harmonics expansions for the temperature and melt radius

$$ \begin{aligned} \overline{T}(r, \theta, \varphi) &= T_0(r) + \sum_{\ell=1}^{+\infty} \sum_{m=-\ell}^{\ell} \overline{T}_{\ell m}(r) Y_{\ell m}(\theta, \varphi), \\ \overline{r_M}(\theta, \varphi) &= \xi_M + \sum_{\ell=1}^{+\infty} \sum_{m=-\ell}^{\ell} \overline{\xi}_{\ell m} Y_{\ell m}(\theta, \varphi), \end{aligned} \quad (41) $$

where $Y_{\ell m}$ is the spherical harmonic of degree $\ell$ and order $m$. Assuming that the topographic changes are small compared to the mean melt radius $\xi_M$ (an approximation also considered by Kvorka and Čadek, 2024), the boundary condition at the solid–liquid interface, held at a constant temperature $T_M$, can be approximated by

$$ \begin{aligned} T_M &= \overline{T}(r_M, \theta, \varphi) \\ &\approx \overline{T}(\xi_M, \theta, \varphi) + \sum_{\ell \neq 0, m} \overline{\xi}_{\ell m} Y_{\ell m}(\theta, \varphi) \frac{dT_0}{dr}(r = \xi_M). \end{aligned} \quad (42) $$

where the quadratic terms that would involve products of spherical harmonics expansions have been neglected (for a similar approach applied to magnetic diffusion, see Styczinski et al., 2022).

At the leading order, $T_0$ is the solution of

$$ \frac{1}{r^2} \frac{d}{dr} \left( r^2 \frac{dT_0}{dr} \right) = 0, \quad T_0(r=r_o) = 0, \quad T_0(r=\xi_M) = T_M, \quad (43) $$

which yields

$$ \frac{dT_0}{dr}(r = \xi_M) = - \frac{T_M r_o}{h_S \xi_M}, \quad h_S = r_o - \xi_M, \quad (44) $$

The first order terms are solution of

$$ \begin{aligned} \frac{1}{r^2} \frac{d}{dr} \left( r^2 \frac{dT_{\ell m}}{dr} \right) - \frac{\ell (\ell + 1)}{r^2} \overline{T}_{\ell m} = 0, \quad \forall (\ell, m) \\ T_{\ell m}(r = r_o) = 0, \\ T_{\ell m}(r = \xi_M) = - \overline{\xi}_{\ell m} \frac{dT_o}{dr}(r = \xi_M). \end{aligned} \quad (45) $$

We seek for solution of the form $\overline{T}_{\ell m}(r) = \alpha r^{\ell} + \beta r^{-\ell-1}$. From the above boundary conditions, one gets

$$ \beta = -\alpha r_o^{2\ell+1}, \quad \alpha = \frac{T_M r_o \overline{\xi}_{\ell m} \xi_M^{\ell}}{h_S (\xi_M^{2\ell+1} - r_o^{2\ell+1})}. \quad (46) $$

This yields the following spherical harmonic coefficients of the temperature field in the solid phase

$$ \overline{T}_{\ell m}(r) = \frac{T_M r_o}{h_S} \frac{\overline{\xi}_{\ell m} \xi_M^{\ell}}{\xi_M^{2\ell+1} - r_o^{2\ell+1}} \left( r^{\ell} - \frac{r_o^{2\ell+1}}{r^{\ell+1}} \right), \quad \forall (\ell, m). \quad (47) $$

From this expression, one can compute the temperature gradient along the solidus, again neglecting the quadratic terms. Introducing

$$ \partial Q(r_M, \theta, \varphi) = \frac{\frac{\partial \overline{T}}{\partial r}(r_M, \theta, \varphi) - \frac{dT_0}{dr}(\xi_M)}{\frac{dT_0}{dr}(\xi_M)}, \quad (48) $$

one gets

$$ \partial Q \approx \sum_{\ell \neq 0, m} \frac{\overline{\xi}_{\ell m} \ell - 1 + (\ell + 2) \eta_S^{2\ell+1}}{1 - \eta_S^{2\ell+1}} Y_{\ell m}(\theta, \varphi), \quad (49) $$

where $\eta_S \equiv \xi_M/r_o$ is the radius ratio of the solidus.

## Data availability

Data will be made available on request.

## References

* Abramson, E.H., Brown, J.M., Slutsky, L.J., 2001. The thermal diffusivity of water at high pressures and temperatures. J. Chem. Phys. 115 (22), 10461–10463. http://dx.doi.org/10.1063/1.1418244.

* Amit, H., Choblet, G., Tobie, G., Terra-Nova, F., Čadek, O., Bouffard, M., 2020. Cooling patterns in rotating thin spherical shells - Application to Titan's subsurface ocean. Icarus 338, 113509. http://dx.doi.org/10.1016/j.icarus.2019.113509.

* Angot, P., Bruneau, C.-H., Fabrie, P., 1999. A penalization method to take into account obstacles in incompressible viscous flows. Numer. Math. 81 (4), 497–520. http://dx.doi.org/10.1007/s002110050401.

* Ascher, U.M., Ruuth, S.J., Spiteri, R.J., 1997. Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Appl. Numer. Math. 25, 151–167. http://dx.doi.org/10.1016/S0168-9274(97)00056-1.

* Ashkenazy, Y., Tziperman, E., 2021. Dynamic europa ocean shows transient Taylor columns and convection driven by ice melting and salinity. Nature Commun. 12, 6376. http://dx.doi.org/10.1038/s41467-021-26710-0.

* Aurnou, J., Heimpel, M., Wicht, J., 2007. The effects of vigorous mixing in a convective model of zonal flow on the ice giants. Icarus 190 (1), 110–126. http://dx.doi.org/10.1016/j.icarus.2007.02.024.

* Backus, G., Parker, R., Constable, C., 1996. Foundations of Geomagnetism. Cambridge University Press, ISBN: 0-521-41006-1.

* Barik, A., Triana, S.A., Calkins, M., Stanley, S., Aurnou, J., 2023. Onset of convection in rotating spherical shells: Variations with radius ratio. Earth Space Sci. 10 (1), e2022EA002606. http://dx.doi.org/10.1029/2022EA002606.

* Bayliss, A., Turkel, E., 1992. Mappings and accuracy for Chebyshev Pseudo-spectral approximations. J. Comput. Phys. 101 (2), 349–359. http://dx.doi.org/10.1016/0021-9991(92)90012-N.

* Beckermann, C., Diepers, H.J., Steinbach, I., Karma, A., Tong, X., 1999. Modeling melt convection in phase-field simulations of solidification. J. Comput. Phys. 154 (2), 468–496. http://dx.doi.org/10.1006/japh.1999.6323.

* Beuthe, M., 2019. Enceladus's crust as a non-uniform thin shell: II tidal dissipation. Icarus 332, 66–91. http://dx.doi.org/10.1016/j.icarus.2019.05.035.

* Bire, S., Kang, W., Ramadhan, A., Campin, J.-M., Marshall, J., 2022. Exploring ocean circulation on Icy moons heated from below. J. Geophys. Res. (Planets) 127 (3), e07025. http://dx.doi.org/10.1029/2021JE007025.
Busse, F.H., 1970. Thermal instabilities in rapidly rotating systems. J. Fluid Mech. 44, 441-460. http://dx.doi.org/10.1017/S0022112070001921

Běhounková, M., Tobie, G., Choblet, G., Kervazo, M., Melwani Daswani, M., Du-
moulin, C., Vance, S.D., 2021. Tidally induced magmatic pulses on the oceanic
floor of Jupiter's Moon Europa. Geophys. Res. Lett. 48 (3), e90077. http://dx.doi.
org/10.1029/2020GL090077.

Cabanes, S., Gastine, T., Fournier, A., 2024. Zonostrophic turbulence in the subsurface oceans of the Jovian and Saturnian moons. Icarus 415, 116047. http://dx.doi.org/10.1016/j.icalarus.2024.116047.

Čadek, O., Souček, O., Běhounková, M., Choblet, G., Tobie, G., Hron, J., 2019. Long-term stability of enceladus' uneven ice shell. Icarus 319, 476-484. http://dx.doi.org/10.1016/j.icalarus.2018.10.003.

Caginalp, G., 1986. An analysis of a phase field model of a free boundary. Arch. Ration. Mech. Anal. 92 (3), 205-245. http://dx.doi.org/10.1007/BF00254827.

Choblet, G., Tobie, G., Sotin, C., Kalousová, K., Grasset, O., 2017. Heat transport in the high-pressure ice mantle of large icy moons. Icarus 285, 252-262. http://dx.doi.org/10.1016/j.icalarus.2016.12.002.

Christensen, U.R., Aubert, J., Cardin, P., Dormy, E., Gibbons, S., Glatzmaier, G.A., Grote, E., Honkura, Y., Jones, C., Kono, M., Matsushima, M., Sakuraba, A., Takahashi, F., Tilgner, A., Wicht, J., Zhang, K., 2001. A numerical dynamo benchmark. Phys. Earth Planet. Inter. 128 (1-4), 25-34. http://dx.doi.org/10.1016/S0031-9201(01)00275-8.

Christensen, U.R., Wicht, J., 2015. 8.10 - Numerical dynamo simulations. In: Schu-
bert, G. (Ed.), Treatise on Geophysics (Second Edition), second ed. Elsevier, Oxford,
ISBN: 978-0-444-53803-1, pp. 245-277. http://dx.doi.org/10.1016/B978-0-444-
53802-4.00145-7.

Corlies, P., Hayes, A.G., Birch, S.P.D., Lorenz, R., Stiles, B.W., Kirk, R., Poggiali, V., Zekber, H., Hess, L., 2017. Titan's topography and shape at the end of the Cassini Mission. Geophys. Res. Lett. 44 (23), 11,754-11,761. http://dx.doi.org/10.1002/2017GL075518.

Couston, L.-A., 2021. Turbulent convection in subglacial lakes. J. Fluid Mech. 915, A31. http://dx.doi.org/10.1017/jfm.2021.38.

Favier, B., Purseed, J., Duchemin, L., 2019. Rayleigh-Bénard convection with a melting boundary. J. Fluid Mech. 858, 437-473. http://dx.doi.org/10.1017/jfm.2018.773.

Gastine, T., Aurnou, J.M., 2023. Latitudinal regionalization of rotating spherical shell convection. J. Fluid Mech. 954, R1. http://dx.doi.org/10.1017/jfm.2022.1010.

Gastine, T., Wicht, J., Aubert, J., 2016. Scaling regimes in spherical shell rotating convection. J. Fluid Mech. 808, 690-732. http://dx.doi.org/10.1017/jfm.2016.659.

Gastine, T., Wicht, J., Aurnou, J.M., 2013. Zonal flow regimes in rotating anelastic spherical shells: An application to giant planets. Icarus 225 (1), 156-172. http://dx.doi.org/10.1016/j.icalarus.2013.02.031.

Gilman, P.A., 1977. Nonlinear dynamics of Boussinesq convection in a deep rotating spherical shell. I. Geophys. Astrophys. Fluid Dyn. 8, 93-135. http://dx.doi.org/10.1080/0309127708240373.

Glatzmaier, G.A., 1984. Numerical simulations of stellar convective dynamos. I. The model and method. J. Comput. Phys. 55, 461-484. http://dx.doi.org/10.1016/0021-9991(84)90033-0.

Gopinath, V., Fournier, A., Gastine, T., 2022. An assessment of implicit-explicit time integrators for the pseudo-spectral approximation of Boussinesq thermal convection in an annulus. J. Comput. Phys. 460, 110965. http://dx.doi.org/10.1016/j.jcp.2022.110965.

Guervilly, C., Cardin, P., Schaeffer, N., 2019. Turbulent convective length scale in planetary cores. Nature 570 (7761), 368-371. http://dx.doi.org/10.1038/s41586-019-1301-5.

Hartmann, R., Stevens, R.J.A.M., Lohse, D., Verzicco, R., 2024. Toward understanding polar heat transport enhancement in subglacial oceans on Icy Moons. Geophys. Res. Lett. 51 (3), e2023GL105401. http://dx.doi.org/10.1029/2023GL105401.

Hay, H.C.F.C., Fenty, I., Pappalardo, R.T., Nakayama, Y., 2023. Turbulent drag at the ice-ocean interface of Europa in simulations of rotating convection: Implications for nonsynchronous rotation of the ice shell. J. Geophys. Res. (Planets) 128 (3), e2022JE007648. http://dx.doi.org/10.1029/2022JE007648.

Hemingway, D.J., Mittal, T., 2019. Enceladus's ice shell structure as a window on internal heat production. Icarus 332, 111-131. http://dx.doi.org/10.1016/j.icalarus.2019.03.011.

Hester, E.W., Couston, L.-A., Favier, B., Burns, K.J., Vasil, G.M., 2020. Improved phase-field models of melting and dissolution in multi-component flows. Proc. R. Soc. Lond. Ser. A 476 (2242), 20200508. http://dx.doi.org/10.1098/rspa.2020.0508.

Hester, E.W., McConnochie, C.D., Cenedese, C., Couston, L.-A., Vasil, G., 2021a. Aspect ratio affects iceberg melting. Phys. Rev. Fluids 6 (2), 023802. http://dx.doi.org/10.1103/PhysRevFluids.6.023802.

Hester, E.W., Vasil, G.M., Burns, K.J., 2021b. Improving accuracy of volume penalised fluid-solid interactions. J. Comput. Phys. 430, 110043. http://dx.doi.org/10.1016/j.jcp.2020.110043.

Hunter, J.D., 2007. Matplotlib: A 2D graphics environment. Comput. Sci. Eng. 9 (3), 90-95. http://dx.doi.org/10.1109/MCSE.2007.55.

Huppert, H.E., 1990. The fluid mechanics of solidification. J. Fluid Mech. 212, 209-240. http://dx.doi.org/10.1017/S0022112090001938.

Jafari-Varzaneh, H.A., Hosseini, S.M., 2015. A new map for the Chebyshev pseudospec-
tral solution of differential equations with large gradients. Numer. Algorithms 69,
95-108. http://dx.doi.org/10.1007/s11075-014-9883-3.

Jansen, M.F., Kang, W., Kite, E.S., Zeng, Y., 2023. Energetic constraints on ocean circulations of Icy Ocean Worlds. psJ 4 (6), 117. http://dx.doi.org/10.3847/PSJ/acda95.

Journaux, B., Kalousová, K., Sotin, C., Tobie, G., Vance, S., Saur, J., Bollengier, O., Noack, L., Rückriemen-Bez, T., Van Hoolst, T., Soderlund, K.M., Brown, J.M., 2020. Large ocean worlds with high-pressure ices. Space Sci. Rev. 216 (1), 7. http://dx.doi.org/10.1007/s11214-019-0633-7.

Julien, K., Rubio, A.M., Grooms, I., Knobloch, E., 2012. Statistical and physical balances in low Rossby number Rayleigh-Bénard convection. Geophys. Astrophys. Fluid Dyn. 106, 392-428. http://dx.doi.org/10.1080/03091929.2012.696109.

Kang, W., 2023. The modulation effect of ice thickness variations on convection in icy ocean worlds. MNRAS 525 (4), 5251-5261. http://dx.doi.org/10.1093/mnras/stad2638.

Kang, W., Jansen, M., 2022. On icy ocean worlds, size controls ice shell geometry. ApJ 935 (2), 103. http://dx.doi.org/10.3847/1538-4357/ac7a32.

Kang, W., Mittal, T., Bire, S., Campin, J.-M., Marshall, J., 2022. How does salinity shape ocean circulation and ice geometry on Enceladus and other icy satellites? Sci. Adv. 8 (29), eabm4665. http://dx.doi.org/10.1126/sciadv.abm4665.

Kihoulou, M., Čadek, O., Kvorka, J., Kalousová, K., Choblet, G., Tobie, G., 2023. Topographic response to ocean heat flux anomaly on the icy moons of Jupiter and Saturn. Icarus 391, 115337. http://dx.doi.org/10.1016/j.icalarus.2022.115337.

King, E.M., Stellmach, S., Aurnou, J.M., 2012. Heat transfer by rapidly rotating Rayleigh-Bénard convection. J. Fluid Mech. 691, 568-582. http://dx.doi.org/10.1017/jfm.2011.493.

Kolomenskiy, D., Schneider, K., 2009. A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles. J. Comput. Phys. 228 (16), 5687-5709. http://dx.doi.org/10.1016/j.jcp.2009.04.026.

Kosloff, D., Tal-Ezer, H., 1993. A modified Chebyshev pseudospectral method with an O (N⁻¹) time step restriction. J. Comput. Phys. 104 (2), 457-469. http://dx.doi.org/10.1006/jcph.1993.1044.

Kvorka, J., Čadek, O., 2022. A numerical model of convective heat transfer in Titan's subsurface ocean. Icarus 376, 114853. http://dx.doi.org/10.1016/j.icalarus.2021.114853.

Kvorka, J., Čadek, O., 2024. The role of subsurface ocean dynamics and phase transitions in forming the topography of icy moons. Icarus 412, 115985. http://dx.doi.org/10.1016/j.icalarus.2024.115985.

Kvorka, J., Čadek, O., Tobie, G., Choblet, G., 2018. Does titan's long-wavelength topography contain information about subsurface ocean dynamics? Icarus 310, 149-164. http://dx.doi.org/10.1016/j.icalarus.2017.12.010.

Labrosse, S., Morison, A., Deguen, R., Alboussière, T., 2018. Rayleigh-Bénard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries. J. Fluid Mech. 846, 5-36. http://dx.doi.org/10.1017/jfm.2018.258.

Lawrence, J.D., Schmidt, B.E., Buffo, J.J., Washam, P.M., Chivers, C., Miller, S., 2024. Ice-ocean interactions on ocean worlds influence Ice Shell topography. J. Geophys. Res. (Planets) 129 (2), e2023JE008036. http://dx.doi.org/10.1029/2023JE008036.

Le Bars, M., Werster, M.G., 2006. Interfacial conditions between a pure fluid and a porous medium: implications for binary alloy solidification. J. Fluid Mech. 550, 149-173. http://dx.doi.org/10.1017/S0022112005007998.

Lefevre, A., Tobie, G., Choblet, G., Čadek, O., 2014. Structure and dynamics of titan's outer icy shell constrained from cassini data. Icarus 237, 16-28. http://dx.doi.org/10.1016/j.icalarus.2014.04.006.

Lemasquerier, D.G., Bierson, C.J., Soderlund, K.M., 2023b. Europa's ocean translates interior tidal heating patterns to the ice-ocean boundary. AGU Adv. 4 (6), e2023AV000994. http://dx.doi.org/10.1029/2023AV000994.

Lemasquerier, D., Favier, B., Le Bars, M., 2023a. Zonal jets experiments in the gas giants' zonostrophic regime. Icarus 390, 115292. http://dx.doi.org/10.1016/j.icalarus.2022.115292.

Long, R.S., Mound, J.E., Davies, C.J., Tobias, S.M., 2020. Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions. J. Fluid Mech. 889, A7. http://dx.doi.org/10.1017/jfm.2020.67.

Marti, P., Schaeffer, N., Hollerbach, R., Cébron, N., Nore, C., Luddens, F., Guer-
mond, J.L., Aubert, J., Takehiro, S., Sasaki, Y., Hayashi, Y.Y., Simitev, R.,
Busse, F., Vantieghem, S., Jackson, A., 2014. Full sphere hydrodynamic and
dynamo benchmarks. Geophys. J. Int. 197 (1), 119-134. http://dx.doi.org/10.1093/gji/
ggt518.

Perissutti, D., Marchioli, C., Soldati, A., 2024. Morphodynamics of melting ice over turbulent warm water streams. Int. J. Multiph. Flow 181, 105007. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2024.105007.

Purseed, J., Favier, B., Duchemin, L., Hester, E.W., 2020. Bistability in Rayleigh-Bénard convection with a melting boundary. Phys. Rev. Fluids 5 (2), 023501. http://dx.doi.org/10.1103/PhysRevFluids.5.023501.

Rabbanipour Esfahani, B., Hirata, S.C., Berti, S., Calzavarini, E., 2018. Basal melting driven by turbulent thermal convection. Phys. Rev. Fluids 3 (5), 053501. http://dx.doi.org/10.1103/PhysRevFluids.3.053501.

Ravichandran, S., Wettlaufer, J.S., 2021. Melting driven by rotating Rayleigh-Bénard convection. J. Fluid Mech. 916, A28. http://dx.doi.org/10.1017/jfm.2021.223.

Rhines, P.B., 1975. Waves and turbulence on a beta-plane. J. Fluid Mech. 69, 417-443. http://dx.doi.org/10.1017/S0022112075001504.
Roberts, J.H., McKinnon, W.B., Elder, C.M., Tobie, G., Biersteker, J.B., Young, D., Park, R.S., Steinbrügge, G., Nimmo, F., Howell, S.M., Castillo-Rogez, J.C., Cable, M.L., Abrahams, J.N., Bland, M.T., Chivers, C., Cochrane, C.J., Dombard, A.J., Ernst, C., Genova, A., Gerekos, C., Glein, C., Harris, C.D., Hay, H.C.F.C., Hayne, P.O., Hedman, M., Hussmann, H., Jia, X., Khurana, K., Kiefer, W.S., Kirk, R., Kivelson, M., Lawrence, J., Leonard, E.J., Lunine, J.I., Mazarico, E., McCord, T.B., McEwen, A., Paty, C., Quick, L.C., Raymond, C.A., Retherford, K.D., Roth, L., Rymer, A., Saur, J., Scanlan, K., Schroeder, D.M., Senske, D.A., Shao, W., Soderlund, K., Spiers, E., Styczinski, M.J., Tortora, P., Vance, S.D., Villarreal, M.N., Weiss, B.P., Westlake, J.H., Withers, P., Wolfenbarger, N., Buratti, B., Korth, H., Pappalardo, R.T., The Interior Thematic Working Group, 2023. Exploring the interior of Europa with the Europa clipper. Space Sci. Rev. 219 (6), 46. http://dx.doi.org/10.1007/s11214-023-00990-y.

Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations. Geochem. Geophys. Geosyst. 14 (3), 751–758. http://dx.doi.org/10.1002/ggge.20071.

Schwaiger, T., Gastine, T., Aubert, J., 2019. Force balance in numerical geodynamics simulations: a systematic study. Geophys. J. Int. 219, S101–S114. http://dx.doi.org/10.1093/gji/ggz192.

Shibley, N.C., Goodman, J., 2024. Europa's coupled ice-ocean system: Temporal evolution of a pure ice shell. Icarus 410, 115872. http://dx.doi.org/10.1016/j. icarus.2023.115872.

Simitev, R., Busse, F.H., 2003. Patterns of convection in rotating spherical shells. New J. Phys. 5 (1), 97. http://dx.doi.org/10.1088/1367-2630/5/1/397.

Soderlund, K.M., 2019. Ocean dynamics of outer solar system satellites. Geophys. Res. Lett. 46 (15), 8700–8710. http://dx.doi.org/10.1029/2018GL081880.

Soderlund, K.M., Kalousová, K., Buffo, J.J., Glein, C.R., Goodman, J.C., Mitri, G., Patterson, G.W., Postberg, F., Rovira-Navarro, M., Rückriemen, T., Saur, J., Schmidt, B.E., Sotin, C., Spohn, T., Tobie, G., Van Hoolst, T., Vance, S.D., Vermeersen, B., 2020. Ice-ocean exchange processes in the Jovian and Saturnian satellites. Space Sci. Rev. 216 (5), 80. http://dx.doi.org/10.1007/s11214-020-00706-6.

Soderlund, K.M., Rovira-Navarro, M., Le Bars, M., Schmidt, B.E., Gerkema, T., 2023. The physical oceanography of ice-covered moons. Annu. Rev. Mar. Sci. 16, 25–53. http://dx.doi.org/10.1146/annurev-marine-040323-101355.

Soderlund, K.M., Schmidt, B.E., Wicht, J., Blankenship, D.D., 2014. Ocean-driven heating of Europa's icy shell at low latitudes. Nat. Geosci. 7 (1), 16–19. http://dx.doi.org/10.1038/ngeo2021.

Song, J., Shishkina, O., Zhu, X., 2024. Scaling regimes in rapidly rotating thermal convection at extreme Rayleigh numbers. J. Fluid Mech. 984, A45. http://dx.doi.org/10.1017/jfm.2024.249.

Stellmach, S., Lischper, M., Julien, K., Vasil, G., Cheng, J.S., Ribeiro, A., King, E.M., Aurnou, J.M., 2014. Approaching the asymptotic regime of rapidly rotating convection: Boundary Layers versus Interior dynamics. Phys. Rev. Lett. 113 (25), 254501. http://dx.doi.org/10.1103/PhysRevLett.113.254501.

Styczinski, M., Vance, S., Harnett, E., Cochrane, C., 2022. A perturbation method for evaluating the magnetic field induced from an arbitrary, asymmetric ocean world analytically. Icarus 376, 114840. http://dx.doi.org/10.1016/j.icarus.2021.114840.

Terra-Nova, F., Amit, H., Choblet, G., Tobie, G., Bouffard, M., Čadek, O., 2023. The influence of heterogeneous seafloor heat flux on the cooling patterns of Ganymede's and Titan's subsurface oceans. Icarus 389, 115232. http://dx.doi.org/10.1016/j. icarus.2022.115232.

Thyng, K.M., Greene, C.A., Hetland, R.D., Zimmerle, H.M., DiMarco, S.F., 2016. True colors of oceanography: Guidelines for effective and accurate colormap selection. Oceanography 29, 9–13. http://dx.doi.org/10.5670/oceanog.2016.66.

Van Hoolst, T., Tobie, G., Vallat, C., Altobelli, N., Bruzzone, L., Cao, H., Dirkx, D., Genova, A., Hussmann, H., Iess, L., Kimura, J., Khurana, K., Lucchetti, A., Mitri, G., Moore, W., Saur, J., Stark, A., Vorburger, A., Wieczorek, M., Aboudan, A., Bergman, J., Bovolo, F., Breuer, D., Cappuccio, P., Carrer, L., Cecconi, B., Choblet, G., De Marchi, F., Fayolle, M., Fienga, A., Futaana, Y., Hauber, E., Kofman, W., Kumamoto, A., Lainey, V., Molyneux, P., Mousis, O., Plaut, J., Puccio, W., Retherford, K., Roth, L., Seignovert, B., Steinbrügge, G., Thakur, S., Tortora, P., Tosi, F., Zannoni, M., Barabash, S., Dougherty, M., Gladstone, R., Gurvits, L.I., Hartogh, P., Palumbo, P., Poulet, F., Wahlund, J.-E., Grasset, O., Witasse, O., 2024. Geophysical characterization of the interiors of ganymede, callisto and Europa by ESA's JUpiter ICy moons explorer. Space Sci. Rev. 220 (5), 54. http://dx.doi.org/10.1007/s11214-024-01085-y.

Vance, S.D., Panning, M.P., Stähler, S., Cammarano, F., Bills, B.G., Tobie, G., Kamata, S., Kedar, S., Sotin, C., Pike, W.T., Lorenz, R., Huang, H.-H., Jackson, J.M., Banerdt, B., 2018. Geophysical investigations of habitability in ice-covered ocean worlds. J. Geophys. Res. (Planets) 123 (1), 180–205. http://dx.doi.org/10.1002/2017JE005341.

Wang, G., Santelli, L., Lohse, D., Verzicco, R., Stevens, R., 2021. Diffusion-free scaling in rotating spherical Rayleigh-Bénard convection. Geophys. Res. Lett. 48 (20), e95017. http://dx.doi.org/10.1029/2021GL095017.

Weady, S., Tong, J., Zidovska, A., Ristroph, L., 2022. Anomalous convective flows carve pinnacles and scallops in melting ice. Phys. Rev. Lett. 128 (4), 044502. http://dx.doi.org/10.1103/PhysRevLett.128.044502.

Weertman, J., 1983. Creep deformation of ice. Annu. Rev. Earth Planet. Sci. 11, 215. http://dx.doi.org/10.1146/annurev.ea.11.050183.001243.

Weller, M.B., Fuchs, L., Becker, T.W., Soderlund, K.M., 2019. Convection in thin shells of icy satellites: Effects of latitudinal surface temperature variations. J. Geophys. Res. (Planets) 124 (8), 2029–2053. http://dx.doi.org/10.1029/2018JE005799.

Wong, T., Hansen, U., Wiesehöfer, T., McKinnon, W.B., 2022. Layering by double-diffusive convection in the subsurface oceans of europa and enceladus. J. Geophys. Res. (Planets) 127 (12), e2022JE007316. http://dx.doi.org/10.1029/2022JE007316.

Worster, M.G., 2000. Solidification of fluids. In: Batchelor, G.K., Moffatt, H.K., Worster, M.G. (Eds.), Perspectives in Fluid Dynamics: A Collective Introduction to Current Research. Vol. 742, Cambridge University Press, Cambridge, ISBN: 978-0-521-53169-6, pp. 393–446.

Yadav, R.K., Gastine, T., Christensen, U.R., Duarte, L.D.V., Reiners, A., 2016. Effect of shear and magnetic field on the heat-transfer efficiency of convection in rotating spherical shells. Geophys. J. Int. 204 (2), 1120–1133. http://dx.doi.org/10.1093/gji/ggv506.

Yang, R., Howland, C.J., Liu, H.-R., Verzicco, R., Lohse, D., 2023a. Bistability in radiatively heated melt ponds. Phys. Rev. Lett. 131 (23), 234002. http://dx.doi.org/10.1103/PhysRevLett.131.234002.

Yang, R., Howland, C.J., Liu, H.-R., Verzicco, R., Lohse, D., 2023b. Ice melting in salty water: layering and non-monotonic dependence on the mean salinity. J. Fluid Mech. 969, R2. http://dx.doi.org/10.1017/jfm.2023.582.

Yang, R., Howland, C.J., Liu, H.-R., Verzicco, R., Lohse, D., 2023c. Morphology evolution of a melting solid layer above its melt heated from below. J. Fluid Mech. 956, A23. http://dx.doi.org/10.1017/jfm.2023.15.