# Topographic response to ocean heat flux anomaly on the icy moons of Jupiter and Saturn

Martin Kihoulou, Ondřej Čadek, Jakub Kvorka, Klára Kalousová, Gael Choblet, Gabriel Tobie

► To cite this version:

Martin Kihoulou, Ondřej Čadek, Jakub Kvorka, Klára Kalousová, Gael Choblet, et al. Topographic response to ocean heat flux anomaly on the icy moons of Jupiter and Saturn. Icarus, 2023, 391, pp.115337. 10.1016/j.icarus.2022.115337. hal-04244311

HAL Id: hal-04244311

https://nantes-universite.hal.science/hal-04244311v1

Submitted on 23 Jul 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Topographic response to ocean heat flux anomaly on
the icy moons of Jupiter and Saturn

Martin Kihoulou<sup>a</sup> (corresponding author), Ondřej Čadek<sup>a</sup>, Jakub Kvorka<sup>a</sup>,
Klára Kalousová<sup>a</sup>, Gaël Choblet<sup>b</sup>, Gabriel Tobie<sup>b</sup>

<sup>a</sup> Charles University, Faculty of Mathematics and Physics, Department of Geo-
physics, V Holešovičkách 2, 180 00 Praha 8, Czech Republic

<sup>b</sup> Laboratoire de Planétologie et Géodynamique, Université de Nantes, 2 rue de la Houssinière, 44322 Nantes, France

1 Highlights:

2

3 The outer shell of an icy moon is an open system exchanging mass with the ocean.

4

5 Topography depends on ice viscosity and the permeability of the phase boundary.

6

7 Large topography is only obtained for viscosities ≳ 10<sup>17</sup> Pa s.

8

9 Heat flux-induced topography of a convecting ice layer is smaller than 10 m.

10

11 The standard formulation of Airy isostasy can lead to an error of 5–15%.

12 *E-mail addresses:* martin.kihoulou@mff.cuni.cz (Martin Kihoulou), ondrej.cadek@mff.cuni.cz (Ondřej

13 Čadek), jakub.kvorka@mff.cuni.cz (Jakub Kvorka), klara.kalousova@mff.cuni.cz (Klára Kalousová), gael.choblet@univ-

14 nantes.fr (Gaël Choblet), gabriel.tobie@univ-nantes.fr (Gabriel Tobie)
Abstract

Recent studies of ocean dynamics suggest that the long-wavelength topography of some icy moons may reflect the phase transitions (melting/freezing) at the interface between the ice shell and the water ocean. Despite the obvious importance of phase changes in the evolution of icy moons, very little is known about how these processes influence their shape, gravity and near-surface stress. Here we address this issue by performing a series of numerical experiments in which we explore the thermo-mechanical response of an ice layer to heat flux variations imposed at the bottom (phase) boundary. We assume that the heat flux from the ocean consists of two components: the heat flux originating in the deep ocean and associated with the global ocean circulation, and the heat flux due to the flow of water generated by variations in the melting temperature along the deformed ice-water interface. The effect of salinity on the heat flux from the ocean is neglected. We demonstrate that the mass exchange between the ocean and the ice layer is a natural consequence of the ice-water phase transition and it occurs in both convection and conduction modes, regardless of whether the system is in equilibrium or not. The magnitude of the heat-flux induced topography strongly depends on the viscosity of ice and the flow of water controlled by the melting temperature. When heat transfer in the ice shell occurs by convection, the surface topography is dominated by small-scale convective features varying in time and its large-scale component does not exceed 10 m. When the viscosity of the ice shell is high (≥ 10¹⁶ Pa s) and the heat is transferred by conduction, the topography is negatively correlated with the heat flux from the ocean and its amplitude increases with increasing viscosity. Topographic amplitudes comparable to those observed on Titan, Enceladus and Dione are obtained only if the water flow associated with lateral variations in the melting temperature is neglected. This suggests that this flow may be too weak
to reduce the variations in ice shell thickness and the motion of water along the phase boundary is more likely to be controlled by other factors, such as variations in salinity and the presence of non-ice material. In addition, we show that the standard formulation of Airy isostasy can lead to an error of 5-15% in determining the variations in ice thickness and we propose a new formulation that takes into account the effect of thermal expansion.

Keywords: icy moons; phase transition; thermal convection; heat flux variations; dynamic topography
51

# 1. Introduction

52

53

54

55

56

57

58

59

60

There is compelling evidence that several moons of Jupiter (Europa and Ganymede) and Saturn (Titan and Enceladus) have oceans of liquid water beneath their solid ice surfaces (see, e.g., Nimmo and Pappalardo, 2016, for a review). A significant portion of the heat required to keep the subsurface oceans in a liquid state is generated in the silicate interior, with much of this energy coming from radioactive decay and/or tidal dissipation. The heat produced in the deep interior of the moon is transmitted through the ocean to the base of the ice shell, potentially controlling ice surface geomorphology (e.g., Goodman and Lenferink, 2012; Askenazy et al., 2018; Kvorka et al., 2018).

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

Most of the previous studies (e.g., Ojakangas and Stevenson, 1989; Nimmo and Bills, 2010) assumed that heat flowing into the ice shell is homogeneously distributed as a consequence of efficient thermal mixing in the ocean. This assumption has been questioned by recent studies of subsurface ocean dynamics which have demonstrated that the heat transfer in the ocean can be affected by the Coriolis effect, leading to an uneven distribution of heat flux at the base of the ice shell (Soderlund et al., 2014; Soderlund, 2019; Amit et al., 2020; Kvorka and Čadek, 2022; Terra-Nova et al., 2022). The models show that depending on the relative importance of rotation, the heat flux from the ocean may peak either at the equator or at the poles, with heat flux variations ranging from a few percent to more than 50% of the mean heat flux. It has been suggested that the heat flux from the ocean can also be affected by heterogeneity at the base of the ocean, originating either in a high-pressure ice layer (Choblet et al., 2017a; Kalousová et al., 2018) or a silicate interior (Choblet et al., 2017b; Běhounková et al., 2021). Another mechanism that can influence the distribution of the heat flux from the ocean involves a meridional overturning circulation driven by variations in salinity, as recently proposed by Lobo et al. (2021). Large-scale variations in heat
flux are required to explain the lateral variations in ice shell thickness on Enceladus inferred from gravity and topography data (Čadek et al., 2019b) and provide a possible explanation for large-scale polar depressions on Titan (Kvorka et al., 2018). Small-scale hydrothermal plumes have been predicted to be responsible for the formation of Europa's chaos regions and other small-scale features (Thomson and Delaney, 2001; Goodman et al., 2004; Goodman and Lenferink, 2012), while large-scale variations in heat flux might control the global dynamics of Europa's ice shell (Soderlund et al., 2014; Askenazy et al., 2018).

Although it is highly likely that the heat flux from the ocean is laterally variable, relatively little is known about how these variations influence the thermomechanical state of the ice shell above the ocean. More attention has been paid to the role of variable surface temperature and tidal heating in the ice shell (e.g., Ojakangas and Stevenson, 1989; Nimmo et al., 2007; Nimmo and Bills, 2010; Askenazy et al., 2018). While these phenomena are likely to affect the thickness of the ice shell in a similar way as the heat flux from the ocean, most of the existing models have a number of limitations and their results cannot be generalized beyond the specific contexts for which they were designed.

Until now, the effect of lateral variations in the heat flux from the ocean has been studied under the assumption that the heat transfer in the ice shell is dominated by conduction. Kvorka et al. (2018) have shown that in this case, a higher-than-average heat flux from the ocean is accompanied by melting of ice and thinning of the ice shell, while freezing of water and thickening of the ice shell occurs in areas where the heat flux from the ocean is lower than average. As a consequence, the resulting ice shell thickness is negatively correlated with the heat flux and its shape follows the Airy compensation scheme. A limitation of the study by Kvorka et al. (2018) is that it does not include lateral variations in material parameters and neglects the changes in the melting temperature along the deformed ice-water
105 interface.

106 When heat is transferred by convection, the ice shell is on average warmer
107 and less viscous than in the conductive case. The loss/gain of ice mass due to
108 melting/freezing at the ice-water interface is compensated by flow in the ice shell
109 and the interface is likely to coincide with an equipotential surface because the
110 ice near the base of the ice shell is not strong enough to support stress (e.g.,
111 Stevenson, 2000; Nimmo and Bills, 2010). The Airy model of isostatic compen-
112 sation is not relevant in this case since the topography of the ice-water boundary
113 is close to zero and the surface topography is generated by the convective flow.
114 It was hypothesized that the convecting ice shell might exhibit Pratt-like thermal
115 isostasy (see, e.g., Figs 2 and 9 in Nimmo and Bills (2010) and Kvorka and Čadek
116 (2022), respectively), but to date, no studies have addressed this issue in a rigorous
117 manner.

118 The question of whether convection occurs in the ice shell or whether heat
119 is transferred only by conduction is difficult because the answer depends on pa-
120 rameters that are not well known, such as the thickness of the ice shell and the
121 viscosity of ice, which in turn depends on the grain size and the temperature of
122 the ocean. When considering the possible role of convection in the heat transfer,
123 it is important to recognize that the ice shell is an open thermodynamic system
124 in the sense that it can exchange both energy and matter with the ocean. This
125 is different from traditional convective systems (such as the Earth's mantle) that
126 are closed, i.e., they have a fixed amount of matter and can only exchange energy
127 with the surroundings. When ice and liquid water are in contact, particles can
128 move from water to ice (freezing) or from ice to water (melting). The direction
129 of this motion and the resulting gains and losses of mass depend on the heat flux
130 at the ice-water interface. The fact that the bottom boundary of the ice shell is
131 permeable may influence the onset of convection as well as the geometry of flow
in the ice shell (Deguen et al., 2013; Deguen, 2013; Labrosse et al., 2018; Morison et al., 2019; Agrusta et al., 2020). The mathematical description of the phase boundary is different from that used for a closed boundary because it requires the knowledge of both temperature and heat flux.

To date, most studies investigating thermal convection in the outer shells of icy moons have assumed for simplicity that the ice shell is a closed system and neglected the mass transport across the ice-water interface (e.g., Tobie et al., 2003; Showman and Han, 2004; Barr and Pappalardo, 2005; Mitri and Showman, 2005; Freeman et al., 2006; Běhounková et al., 2013; Peddinti and McNamara, 2019; Kalousová and Sotin, 2020; Harel et al., 2020; Green et al., 2021). However, as shown by recent studies of the Earth's core and early magma ocean, the presence of a phase boundary may significantly influence the style of solid-state convection, leading to effects that could not have been predicted using a standard modeling approach (Deguen et al., 2013; Deguen, 2013; Labrosse et al., 2018; Morison et al., 2019; Agrusta et al., 2020). These studies argue that since the melting temperature is pressure-dependent, the deformation of the phase boundary is always accompanied by a flow of the liquid phase generated by temperature changes along the phase boundary. Depending on the speed of the flow, the heat transported by the liquid can enhance mass transport through the boundary and change the dynamics of convective heat transfer in the solid phase.

In this study, we investigate the behavior of an ice layer that is unevenly heated from the ocean. In order to understand the role of heat flux variations in the evolution of the layer, we perform a series of numerical experiments in which we vary the viscosity of ice and the parameter characterizing the water flow driven by lateral variations in the melting temperature. For each model, we evaluate the topography, the degree of compensation and the near-surface stress, i.e., the characteristics that can be measured or estimated by remote sensing. The main
159 purpose of this study is to provide a diagnostic tool for understanding the link between the heat flux from the ocean and the topographic response of the ice shell 160 and to show how this response depends on the mode of heat transfer. 161

## 2. Method 162

### 2.1. Governing equations and boundary conditions 163

The results presented in this study are obtained by simultaneously solving the transport equations for mass, momentum and energy in the Boussinesq approximation and the constitutive equation for a viscous fluid: 164 165 166

$$ \nabla \cdot \boldsymbol{v} = 0, \qquad (1) $$

$$ \nabla \cdot \boldsymbol{\sigma} - \rho_0(1 - \alpha(T - T_0))g\boldsymbol{e}_y = \mathbf{0}, \qquad (2) $$

$$ \rho_0 c_p \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) - \rho_0 c_p \boldsymbol{v} \cdot \nabla T, \qquad (3) $$

$$ \boldsymbol{\sigma} = -p\mathbf{I} + \eta (\nabla \boldsymbol{v} + (\nabla \boldsymbol{v})^T). \qquad (4) $$

Here, $\boldsymbol{v}$ is the velocity, $\boldsymbol{\sigma}$ is the stress tensor, $\rho_0$ is the density of ice at a reference temperature $T_0$, $\alpha$ is the coefficient of thermal expansion, $T$ is the temperature, $g$ is the gravitational acceleration, $\boldsymbol{e}_y$ is the unit vector pointing upward, $c_p$ is the specific heat capacity, $k$ is the thermal conductivity, $p$ is the pressure, $\mathbf{I}$ is the identity tensor, $\eta$ is the viscosity and the superscript $T$ denotes transposition of a matrix. 174 175 176 177 178 179

At time $t = 0$, the computational domain is rectangular with an aspect ratio of 16 (Fig. 1). This aspect ratio is a compromise between the computational cost, which increases with the width of the domain, and the realistic representation of the ice shell. As shown in studies of ocean dynamics (Soderlund, 2019; Amit et al., 2020; Kvorka and Čadek, 2022), the time-averaged heat flux from the ocean is dominated by zonal harmonics of degree 2–6, corresponding to a wavelength that
a) Boundary conditions

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

b) Heat flux from the ocean

<figure><img src="image_3.png" /><figcaption>Figure 1: a) Sketch of the computational domain and b) the lateral dependence of the heat flux, $q_o$, prescribed at the surface of the ocean.</figcaption></figure>

is much larger than the thickness of ice. The vertical boundaries of the domain are impermeable to mass and heat transfer and their position is fixed:

$$ \sigma_{xy} = 0, v_x = 0, \partial T / \partial x = 0 \text{ for } x = 0 \text{ and } x = l, \qquad (5) $$

where $x$ and $y$ are the Cartesian coordinates and $l$ is the length of the computational domain. The shape of the upper and bottom boundaries can change over time. The surface is stress-free and isothermal while the hydrostatic pressure and the melting temperature are prescribed at the ice-water interface:

$$ T = 90 \, \text{K}, \quad \mathbf{n} \cdot \boldsymbol{\sigma} = \mathbf{0} \quad \text{for } y = D + h_s(x, t), \qquad (6) $$

$$ T = 261.6 \, \text{K} + \frac{dT_m}{dy} h_{i/w}, \quad \mathbf{n} \cdot \boldsymbol{\sigma} = -n g (D \rho_i - h_{i/w} \rho_w) \quad \text{for } y = h_{i/w}(x, t), \qquad (7) $$

where $D$ is the average thickness of the ice layer, $\mathbf{n}$ is the unit outer normal, $T_m$ is the melting temperature of ice, $\rho_i$ is the density of ice, $\rho_w$ is the density of water, and $h_s$ and $h_{i/w}$ are the topography of the surface and the ice-water boundary,
respectively. The parameters of the model are chosen to roughly correspond to the ice shell of Titan. Note that $dT_m/dy$ depends on the gravitational acceleration and pressure ($dT_m/dy = -\rho_w g dT_m/dp$). On Enceladus, where $g$ is twelve times lower than on Titan and the average thickness of ice is less than 25 km, the value of $dT_m/dy$ would be a factor of 15 lower than that used in this study.

The evolution of the upper (material) boundary is described by the equation (e.g., Dingemans, 1997)

$$ \frac{\partial h_s}{\partial t} = v_y - \frac{\partial h_s}{\partial x} v_x, \qquad (8) $$

where $v = (v_x, v_y)$ is the velocity of ice flow at the surface. The shape of the bottom (phase) boundary is determined by two factors: the velocity of ice flow, $v$, and the velocity due to the phase change, $u$. Eq. (8) must then be modified as follows (Slattery et al., 2007):

$$ \frac{\partial h_{i/w}}{\partial t} = v_y + u_y - \frac{\partial h_{i/w}}{\partial x}(v_x + u_x), \qquad (9) $$

where

$$ \mathbf{u} = (u_x, u_y) = \frac{\mathbf{n} \cdot (\mathbf{q}_i - \mathbf{q}_w)}{L\rho_i} \mathbf{n}. \qquad (10) $$

Here, $L$ is the latent heat of ice and $q_i$ and $q_w$ are the heat flux in ice and water, respectively. Inspection of Eq. (9) shows that a steady state ($\partial h_{i/w}/\partial t = 0$) is reached when $u = -v$, i.e., when the loss/gain of ice mass caused by melting/freezing is compensated by flow in the ice shell. Eq. (10) can then be rewritten as

$$ -L\rho_i v \cdot n = n \cdot (q_i - q_w). \qquad (11) $$

If the topography of the ice-water interface is non-zero, the solution of Eqs. (1)-(4) is non-trivial ($v \neq 0$), implying that in steady state, the heat flux at the ice-water boundary must be discontinuous ($n \cdot q_i \neq n \cdot q_w$). This fact has been overlooked in most previous studies on thermally induced ice thickness variations (see, e.g., Ojakangas and Stevenson, 1989).
We assume that the heat flux from the ocean can be expressed as the sum of two contributions:

$$ q_w = q_o + q_\gamma, \qquad (12) $$

where $q_o$ is the heat flux originating in the deep ocean and associated with the global ocean circulation (cf. Soderlund, 2019; Amit et al., 2020; Kvorka and Čadek, 2022), and $q_\gamma$ is the heat flux due to the flow of water generated by variations in the melting temperature along the deformed ice-water interface. Following previous studies (Deguen et al., 2013; Deguen, 2013; Labrosse et al., 2018), $q_\gamma$ is assumed to be proportional to the phase boundary topography,

$$ q_\gamma \propto -\rho_w c_{p,w} v_w \frac{dT_m}{dp} \frac{dp}{dy} h_{i/w} \mathbf{n} = \rho_w^2 g c_{p,w} v_w \frac{dT_m}{dp} h_{i/w} \mathbf{n}, \qquad (13) $$

where $c_{p,w}$ is the specific heat capacity of water and $v_w$ is the speed of the water flow. Therefore, the heat flux $q_\gamma$ can be parameterized by a single parameter, $\gamma$, characterizing the response of the ocean to temperature variations along the phase boundary:

$$ q_\gamma = \gamma h_{i/w} \mathbf{n}. \qquad (14) $$

As shown by Labrosse et al. (2018), an increase in $\gamma$ leads to an enhancement of mass transport through the boundary and to a reduction of the boundary topography. Note that the notation used here is different from that used in Labrosse et al. (2018) where the degree of permeability of the boundary is parameterized in terms of the phase-change dimensionless number $\Phi = \rho_i(\rho_w - \rho_i)DLg\eta^{-1}\gamma^{-1}$. The reason why we use a different parameterization is that unlike $\Phi$, the parameter $\gamma$ is independent of the viscosity of ice, which is a function of stress in our model, and can be readily interpreted in terms of volumetric heating (the unit of $\gamma$ is W/m$^3$).

The heat flux associated with the global ocean circulation (Fig. 1b) is prescribed as follows:

$$ q_o(x) = -\bar{q} \left( 1 - \frac{1}{2} \cos(2\pi x/l) \right) \mathbf{n}, \qquad (15) $$
where $\bar{q}$ is the average heat flux at the ice-water interface after reaching a steady state. The magnitude of $\bar{q}$ depends on the mode of heat transfer in the ice layer and ranges from 6 mW/m$^2$ for the conductive case to 25 mW/m$^2$ for the convection case (for comparison, the heat flux on Titan is estimated to be between 3 and 7 mW/m$^2$, see Kalousová and Sotin (2020)). The amplitude of the heat flux changes ($\bar{q}/2$) is chosen in agreement with previous studies of ocean dynamics which show that the time-averaged heat flux from the ocean can vary by tens of percent (see, e.g., Fig. 4 in Kvorka and Čadek (2022)).

The definition of $q_o$ in terms of a cosine function allows the resulting topography to be extrapolated to larger scales and compared with the spherical harmonic models of topography on Titan and other icy moons. The extrapolation can be performed based on the results of previous studies (e.g., Kvorka et al., 2018) or topographic kernels calculated in spherical geometry. However, the reader should be aware that this extrapolation only indicates the trends and should not be used to interpret the observed topography of icy moons.

It is worth noting that the melting and freezing at the base of the ice shell can also be caused by tidal heating within the shell (e.g., Tobie et al., 2003; Nimmo and Bills, 2010; Kvorka et al., 2018) and/or laterally varying temperature at the surface of the moon (e.g., Askenazy et al., 2018). These effects are not considered here and will be left for a future study.

When the ice layer is subject to variations in heat flux from the ocean, the velocity field obtained by solving Eqs. (1)–(4) with conditions (5)–(10) is always non-zero, implying that the advection term, $v \cdot \nabla T$, cannot, in general, be omitted in the heat equation. In sections 3 and 4, we use the term “convection” to describe the models where the vertical heat transfer occurs due to a spontaneous flow of ice driven by buoyancy forces, while the term “conduction” refers to the case where the effect of advection is negligible except for a narrow low viscosity zone
<table><caption>Table 1: Creep parameters used in this study (adapted from Goldsby and Kohlstedt, 2001).</caption><thead><tr><th></th><th>T (K)</th><th>B (Pa<sup>-n</sup> K<sup>r</sup> m<sup>p</sup> s<sup>-1</sup>)</th><th>r</th><th>p</th><th>n</th><th>Q (kJ mol<sup>-1</sup>)</th></tr></thead><tbody><tr><td>η<sub>Vdiff</sub></td><td>—</td><td>1.4 &middot; 10<sup>-7</sup></td><td>1</td><td>2</td><td>1</td><td>59.4</td></tr><tr><td>η<sub>GBdiff</sub></td><td>—</td><td>3.0 &middot; 10<sup>-16</sup></td><td>1</td><td>3</td><td>1</td><td>49</td></tr><tr><td rowspan="2">η<sub>disl</sub></td><td>&lt; 258</td><td>3.1 &middot; 10<sup>-18</sup></td><td>0</td><td>0</td><td>4</td><td>60</td></tr><tr><td>&gt; 258</td><td>4.7 &middot; 10<sup>5</sup></td><td>0</td><td>0</td><td>4</td><td>180</td></tr><tr><td rowspan="2">η<sub>GBS</sub></td><td>&lt; 255</td><td>1.4 &middot; 10<sup>-13</sup></td><td>0</td><td>1.4</td><td>1.8</td><td>49</td></tr><tr><td>&gt; 255</td><td>1.1 &middot; 10<sup>16</sup></td><td>0</td><td>1.4</td><td>1.8</td><td>192</td></tr><tr><td>η<sub>BS</sub></td><td>—</td><td>7.1 &middot; 10<sup>-7</sup></td><td>0</td><td>0</td><td>2.4</td><td>60</td></tr></tbody></table>

near the base of the ice layer. The relative importance of heat advection can be characterized by the Nusselt number, $Nu = Q_{tot}/Q_{diff}$, where $Q_{tot}$ is the total heat flow over the upper boundary obtained by solving Eqs. (1)–(4) and $Q_{diff}$ is the heat flow computed from the one-dimensional thermal diffusion equation ($k(T)dT/dy = \text{const}$). While the Nusselt number of the convection models presented in section 3 is always larger than 1.9, the heat flow obtained for models of conductive heat transfer (section 4) is almost identical to that computed from the diffusion equation ($1 < Nu < 1.01$).

## 2.2. Material parameters

As shown in previous studies (e.g., Kvorka et al., 2018; Akiba et al., 2022), viscosity of ice is the key parameter determining the magnitude of both the surface and ice-water topography, which underlines the importance of realistic constitutive relations for predicting the correct topographic response. Based on the results of Goldsby and Kohlstedt (2001), we assume that deformation of ice is accommodated by four creep mechanisms: volume and grain boundary diffusion, dislocation creep,
basal slip (BS) and grain boundary sliding (GBS). The effective viscosity η of ice is then given by

$$ \frac{1}{\eta} = \frac{1}{\eta_{\text{Vdiff}}} + \frac{1}{\eta_{\text{GBdiff}}} + \frac{1}{\eta_{\text{disl}}} + \frac{1}{\eta_{\text{GBS}} + \eta_{\text{BS}}} + \frac{1}{\eta_{\text{max}}}, \quad (16) $$

where ηmax is the cut-off viscosity representing the effective viscosity of ice at low temperatures. Since the effective viscosity of low-temperature ice depends on a number of factors, such as the presence of ice contaminants, deformation history or porosity, and is not well constrained, we treat it as a free parameter. As will be shown in section 4, the choice of ηmax affects the topographic response of the ice layer only if ηmax ≪ 10²² Pa s.

The viscosities of individual creep mechanisms can be expressed as

$$ \eta_c = \frac{d^p T^r}{2B\sigma_{II}^{n-1}} \exp\left(\frac{Q}{RT}\right). \quad (17) $$

Here, index c represents the specific creep mechanism, d is the grain size, R is the universal gas constant, B, Q, p, r and n are material parameters (see Table 1), and σII is the second invariant of the deviatoric stress tensor. For simplicity, we assume that the grain size is constant in time and space. Note that the values of parameter B in Table 1 are different from the values of parameter A given in Table 5 in Goldsby and Kohlstedt (2001). The reason for this difference is that Eq. (17) is expressed in terms of σII, while the constitutive relations in Goldsby and Kohlstedt (2001) are formulated in terms of the differential stress (B = 3(n+1)/2A/2, see also Durham and Stern (2001)).

The thermal expansion coefficient, the thermal conductivity and the specific heat capacity depend on temperature as follows (Hobbs, 1974; Röttger et al., 1994; Feistel and Wagner, 2006):

$$ \alpha(T) = \frac{1}{V} \frac{\partial V}{\partial T}, \quad V(T) = A_0 + A_3 T^3 + A_4 T^4 + A_5 T^5 + A_6 T^6 + A_7 T^7, \quad (18) $$
<table><caption>Table 2: Parameters of the model. The value of $dT_m/dy$ is calculated using the formula in Marcolli (2017). For the creep parameters, see Table 1.</caption><thead><tr><th>Symbol</th><th>Meaning</th><th>Value</th></tr></thead><tbody><tr><td>A<sub>0</sub></td><td>parameter in Eq. (18)</td><td>128.2147 Å<sup>3</sup></td></tr><tr><td>A<sub>3</sub></td><td>parameter in Eq. (18)</td><td>-1.3152 · 10<sup>-6</sup> Å<sup>3</sup>K<sup>-3</sup></td></tr><tr><td>A<sub>4</sub></td><td>parameter in Eq. (18)</td><td>2.4837 · 10<sup>-8</sup> Å<sup>3</sup>K<sup>-4</sup></td></tr><tr><td>A<sub>5</sub></td><td>parameter in Eq. (18)</td><td>-1.6064 · 10<sup>-10</sup> Å<sup>3</sup>K<sup>-5</sup></td></tr><tr><td>A<sub>6</sub></td><td>parameter in Eq. (18)</td><td>4.6097 · 10<sup>-13</sup> Å<sup>3</sup>K<sup>-6</sup></td></tr><tr><td>A<sub>7</sub></td><td>parameter in Eq. (18)</td><td>-4.9661 · 10<sup>-16</sup> Å<sup>3</sup>K<sup>-7</sup></td></tr><tr><td>d</td><td>grain size</td><td>0.3 – 1000 mm</td></tr><tr><td>D</td><td>thickness of the ice shell</td><td>100 km</td></tr><tr><td>c<sub>1</sub></td><td>parameter in Eq. (20)</td><td>185 J K<sup>-1</sup>kg<sup>-1</sup></td></tr><tr><td>c<sub>2</sub></td><td>parameter in Eq. (20)</td><td>7.037 J K<sup>-2</sup>kg<sup>-1</sup></td></tr><tr><td>c<sub>p</sub></td><td>specific heat capacity</td><td>see Eq. (20)</td></tr><tr><td>h<sub>s</sub></td><td>topography of the surface</td><td>see Eq. (8)</td></tr><tr><td>h<sub>i/w</sub></td><td>topography of ice-water boundary</td><td>see Eq. (9)</td></tr><tr><td>g</td><td>gravitational acceleration</td><td>1.35 m s<sup>-2</sup></td></tr><tr><td>k</td><td>thermal conductivity</td><td>see Eq. (19)</td></tr><tr><td>k<sub>1</sub></td><td>parameter in Eq. (19)</td><td>0.4685 W m<sup>-1</sup>K<sup>-1</sup></td></tr><tr><td>k<sub>2</sub></td><td>parameter in Eq. (19)</td><td>488.12 W m<sup>-1</sup></td></tr><tr><td>l</td><td>length of the domain</td><td>1600 km</td></tr><tr><td>L</td><td>latent heat of ice</td><td>334 kJ kg<sup>-1</sup></td></tr><tr><td>R</td><td>universal gas constant</td><td>8.314 J mol<sup>-1</sup>K<sup>-1</sup></td></tr><tr><td>T<sub>0</sub></td><td>reference temperature</td><td>273 K</td></tr><tr><td>T<sub>i/w</sub></td><td>temperature at ice-water interface</td><td>equal to melting temperature</td></tr><tr><td>T<sub>m</sub></td><td>melting temperature</td><td>261.6 K+(dT<sub>m</sub>/dy)h<sub>i/w</sub></td></tr><tr><td>dT<sub>m</sub>/dy</td><td>vertical gradient of T<sub>m</sub></td><td>1.4 · 10<sup>-4</sup> K m<sup>-1</sup></td></tr><tr><td>T<sub>s</sub></td><td>temperature at the surface</td><td>90 K</td></tr><tr><td>α</td><td>coefficient of thermal expansion</td><td>see Eq. (18)</td></tr><tr><td>γ</td><td>scaling parameter in Eq. (14)</td><td>0 – 10 mW m<sup>-3</sup></td></tr><tr><td>η<sub>0</sub></td><td>viscosity at temperature T<sub>i/w</sub></td><td>10<sup>13</sup> – 10<sup>20</sup> Pa s</td></tr><tr><td>η<sub>max</sub></td><td>cut-off viscosity</td><td>10<sup>21</sup> – 10<sup>24</sup> Pa s</td></tr><tr><td>ρ<sub>0</sub></td><td>density of ice at temperature T<sub>0</sub></td><td>920 kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>i</sub></td><td>density of ice</td><td>depending on T</td></tr><tr><td>ρ<sub>w</sub></td><td>density of water</td><td>1000 kg m<sup>-3</sup></td></tr></tbody></table>
$$k(T) = k_1 + \frac{k_2}{T}, \qquad (19)$$

$$c_p(T) = c_1 + c_2 T, \qquad (20)$$

where $A_0$, $A_3-A_7$, $k_1$, $k_2$, $c_1$ and $c_2$ are constants whose values are given in Table 2. Most of the previous studies investigating the heat transfer across the outer shells of icy moons have assumed that the thermal expansivity is independent of temperature. Since the value of $\alpha$ is about five times smaller at the surface than at the base of the ice shell, this may have led to an overestimation of topographic amplitudes. As it will be discussed in section 5, an accurate approximation of $\alpha(T)$ is needed to correctly predict the surface topography and to assess the degree of isostatic compensation.

The problem described above is solved using the FEniCS software package (Logg et al., 2012; Alnaes et al., 2015). As the shape of the computational domain changes, the finite element mesh is modified using the Arbitrary-Lagrangian-Eulerian (ALE) method (Donea et al., 2004).

To understand the effect of the heat flux on the deformation of the ice layer, we perform a parametric study in which we vary parameters $d$ (0.3 – 1000 mm), $\eta_{max}$ (10$^{21}$ – 10$^{24}$ Pa s) and $\gamma$ (0 – 10 mW/m$^3$). The resulting models are discussed in terms of quantities that can be inferred from remote sensing data. In particular, we focus on the dynamic topography, the orientation of stresses near the surface and the degree of compensation. Note that the grain size in our simulations is as large as 1 m, which is significantly more than the values obtained from models of grain size evolution (e.g., Barr and McKinnon, 2007) and analysis of ice samples from the Antarctic ice sheet (Durand et al., 2006). Depending on the presence of impurities in ice, the grain size is expected to be between 1 and 10 mm, corresponding to $\eta_{Vdiff}$ at the melting point between 10$^{15}$ and 10$^{17}$ Pa s. However, in absence of impurities, the grain size may be larger than 100 mm, as observed in the accreted
ice above Lake Vostok (Souchez et al., 2003). The viscosity $\eta_{Vdiff}$ corresponding
to $d = 100$ mm is about $10^{19}$ Pa·s, which indicates that diffusion creep may no
longer be the dominant deformation mechanism in the warm part of the ice shell.

The reader should be aware that our model is simplified in two respects before
applying it to the real world. First, we assume that the ice layer is made of pure
water ice and we do not take into account the effects of other chemical compounds
(salts, clathrates, organics, etc.) that are likely to be present in the interior of
icy moons and can strongly influence the dynamics of phase transition and grain
size evolution (e.g., Han and Showman, 2005; Tobie et al., 2012; Buffo et al., 2020,
2021). Second, the effect of ocean circulation on the thermo-mechanical evolution
of the ice shell is approximated by Eq. (12), which is a clear simplification. The
relationship between the heat flux and the topography ($q_\gamma \propto h_{i/w}$), proposed by
Deguen et al. (2013), is intuitively plausible but it has not yet been verified. In
addition, it is not guaranteed that the heat flux from the ocean can be expressed
as a superposition of two independent components, $q_\gamma$ and $q_o$. In spite of these
limitations, we believe that our results provide a useful insight into the dynamics
of icy moons and can serve as a basis for more complex studies of ice-water systems.

362

3. Convective heat transfer

Our numerical simulations show that in the case of a 100 km thick ice layer,
convection occurs only when $d \lesssim 2$ mm, i.e., when the viscosity at the base of the
ice layer is lower than $\sim 2 \cdot 10^{15}$ Pa·s. First, we consider a model in which the
heat flux from the deep ocean does not vary laterally ($q_o = -\bar{q}n$) and the mass
transport across the ice-water boundary is only controlled by the parameter $\gamma$ (see
Eqs. (12)–(15)).

364

365

366

367

368

369
Fig. 2 shows four typical snapshots of the temperature field and flow characteristics obtained for the models with $d = 1$ mm and $\gamma$ ranging between 0 and 1 W/m$^3$. As already discussed by Labrosse et al. (2018), the oceanic heat flow, $q_{\gamma}$, due to lateral variation in the melting temperature has a decisive effect on the permeability of the phase boundary. When $q_{\gamma} = 0$ (Fig. 2a), the normal velocity of the phase change, $u_n$, is correlated with the heat flux in ice and is about the same magnitude as the velocity of ice ($u_n \approx -v_n$). The heat flux in ice varies laterally by tens of mW/m$^2$ while the melting rate is of the order of 1 mm/yr. Similar values of $q_i$ and $u_n$ are obtained for $\gamma = 0.1$ mW/m$^3$ (Fig. 2b). As in the case of $\gamma = 0$, ice melting occurs beneath the regions of upward ice flow (plumes) while freezing of water occurs beneath the regions of downwelling, suggesting that the phase change at the ice-water boundary is mainly controlled by heat flux anomalies in the ice. However, the heat flux $q_{\gamma}$ is comparable in magnitude with $q_i$ which indicates that the effect of lateral variations in the melting temperature can no longer be neglected. When $\gamma$ is further increased (to 10 mW/m$^3$, Fig. 2c), the heat flux from the ocean becomes the dominant driver of the phase change (cf. Labrosse et al., 2018). The amplitude of $q_{\gamma}$ is about 1.5 W/m$^2$, which is about 100 times more than that of $q_i$. Unlike the case where $\gamma \lesssim 0.1$ mW/m$^3$, freezing of water occurs beneath the regions of upward flow while ice melting occurs beneath the regions of downwelling. The change of the flow pattern is accompanied by the enhancement of the mass transported through the phase boundary, as evidenced by the shape of streamlines and the velocities at the base of the ice shell. This behavior is even more pronounced when $\gamma = 1$ W/m$^3$ (Fig. 2d). Note, that in this case, the variations in $q_i$ are negligible in comparison with the variations in $q_w$, which is the opposite to what is found for $\gamma = 0$.

The conclusions drawn from the numerical simulations illustrated in Fig. 2 are also valid in the case when the heat flux from the deep ocean varies laterally (Eq.
<figure><img src="image_4.png" /><figcaption>Figure 2: Snapshots of the temperature field and flow characteristics obtained for models with $d = 1$ mm, $q_o = \text{const}$ and a) $\gamma = 0$, b) 0.1 mW/m$^3$, c) 10 mW/m$^3$ and d) 1000 mW/m$^3$. The flow direction is shown by white streamlines while color shading shows the temperature. The upper line graphs show the heat flux $q_i$ (blue) and $q_w$ (red) at the ice-water boundary. The lower line graphs show the phase change velocity, $u_n \equiv n \cdot u$ (blue), and the normal flow velocity at the base of the ice shell, $v_n \equiv n \cdot v$ (red). Note that each line graph has a different vertical scale and $u_n$ is plotted with the opposite sign to facilitate the comparison with the ice flow velocity. The numerical simulations were carried out for an aspect ratio of 4 and $\eta_{\text{max}} = 10^{24}$ Pa s. To facilitate a comparison with previous studies, we also provide the values of the phase-change dimensionless number $\Phi$ (see section 2).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: Snapshots of the temperature field and surface topography (black lines) obtained for models with $d = 1$ mm, different values of $\gamma$ and laterally varying $q_o$. The red curves in the line graphs represent the long-wavelength (cosine) component of the topographic signal. All models are computed for $\eta_{max} = 10^{24}$ Pa s.</figcaption></figure>

15). As shown in Figs 3 and 4, the surface topography of the convecting layer (black lines) is characterized by several small-scale features associated with the convective flow and the ice thickness remains roughly constant on a large scale. The magnitude of the long-wavelength surface topography (red curves) is less than 10 m. Such topography would be extremely difficult to detect as it is much smaller than the short-wavelength topography associated with tectonic features, which can reach up to a few hundreds of meters (e.g., Nimmo and Schenk, 2006; Radebaugh et al., 2007; Giese et al., 2008).
<figure><img src="image_6.png" /><figcaption>Figure 4: As in Fig. 3 but for $d = 0.3$ mm.</figcaption></figure>

4. Conductive heat transfer

When the heat is transferred by conduction ($d > 2$ mm), the increase in heat flux from the ocean leads to an uplift of the ice-water interface and a large-scale subsidence of the surface, resulting in compressive stresses in the central part of the domain (Fig. 5). The resulting ice shell thickness is negatively correlated with the heat flux (cf. Fig. 1b) and its shape can be approximated by Airy isostasy (see also section 5). As already discussed in Labrosse et al. (2018), the variations in ice shell thickness are reduced when topography-induced variations in the heat flux are included.

Fig. 6 shows the dependence of the magnitude of the surface and ice-water interface topography on the parameters $d$, $\eta_{max}$ and $\gamma$. The amplitude of the surface topography increases with increasing grain size and decreases with increasing $\gamma$.
<figure><img src="image_7.png" /><figcaption>Figure 5: Surface topography (a), ice-water boundary topography (b) and surface stress (c) obtained for models with $d = 10$ mm, $\eta_{max} = 10^{24}$ Pa s and different values of $\gamma$. The surface stress is calculated as $\sigma_{xx}(x) = [\sigma(x, D+h_s) \cdot e_x] \cdot e_x$ where $\sigma$ is the Cauchy stress tensor.</figcaption></figure>
When $\gamma = 0$ (Fig. 6a), the magnitude of the surface topography ranges from 20 m (d = 3 mm, $\eta_0 \approx 10^{16}$ Pa s) to 250 m (d = 1000 mm, $\eta_0 \approx 10^{19}$ Pa s). These values are roughly consistent with those obtained by Kvorka et al. (2018) for the spherical degree $l = 8$ (the wavelength of 1600 km used in the present study corresponds to $l \approx 10$).

The magnitude of the topography is strongly influenced by the non-Newtonian behavior of ice at large grain sizes. If non-linear creep mechanisms were not included, the surface topography would exceed 2 km for $d > 30$ mm (not shown in Fig. 6). The reason for this difference is that the stress-dependent viscosity is significantly lower than the viscosity of diffusion creep when the topographic amplitude (and hence the magnitude of the deviatoric stress) is large. In contrast, the topography is not much affected by the choice of parameter $\eta_{max}$, characterizing the mechanical behavior of the uppermost part of the ice layer. The effect of $\eta_{max}$ on the deformation is significant only if $\eta_{max} < 10^{22}$ Pa s and $d > 10$ mm. When $\eta_{max} \gtrsim 10^{23}$ Pa s, all models with the same grain size predict about the same topography. Detailed analysis of the results shows that the relative importance of parameter $\eta_{max}$ decreases as the viscosity contrast across the layer is increased and is negligible when $\eta_{max}/\eta_0 \gtrsim 10^6$ (see also discussion in Kvorka et al. (2018)). The slope of the curves in Fig. 6a decreases with the increasing grain size but it does not reach the zero value even for $d = 1000$ mm. This means that even for the largest grain size considered, the velocity of ice near the phase boundary is non-zero, implying that $q_i \neq q_w$.

The inclusion of the heat flux $q_\gamma$, Eqs (13) and (14), does not change the shape of the topographic signal (Fig. 5) but dramatically affects its amplitude (Fig. 6b). Since the heat flux $q_\gamma \propto -h_{i/w}$ (Eq. 14) and the phase boundary topography is correlated with the heat flux from the deep ocean ($h_{i/w} \propto q_o$, Fig. 5), the heat flux $q_\gamma$ tends to counterbalance the lateral variations of $q_o$ and, as $\gamma$ increases, the total
heat flux from the ocean, $q_w = q_o + q_\gamma$, goes to a constant ($q_w \to \bar{q}$). When $\gamma = 0.1$ mW/m$^3$, the amplitude of the surface topography is less than 2.5 m, suggesting that the heat flux from the deep ocean has almost no effect on the shape of the ice shell if $\gamma$ is large.

The near-surface stress (Fig. 7) is equally sensitive to changes in the grain size as the topography. The magnitude of the stress increases with increasing $\eta_{max}$ (Fig. 7a) and decreasing $\gamma$ (Fig. 7b), reaching a value of $\sim 100$ kPa for $\gamma = 0$ and $d \gtrsim 10$ mm. If $\gamma$ is larger than 0.1 mW/m$^3$, the stress drops below 3 kPa.

Fig. 8 shows the flow of ice and the rate of the phase change for selected values of $d$ and $\gamma$. The flow is confined to a narrow boundary layer adjacent to the ice-water interface where it reaches a horizontal velocity of $\sim 1 - 10$ mm/yr. Inside the ice shell the velocity of ice flow is significantly lower than near the boundary confirming that the role of advection in heat transfer is negligible. The rate of the phase change (line graphs) is correlated with the heat flux imposed at the base of the layer and its magnitude decreases with increasing $\gamma$.

So far, we have considered the case where the average thickness of the ice layer is 100 km. The dependence of the topographic response on the thickness of ice, $D$, is illustrated in Fig. 9. If $\gamma = 0$, a reduction of $D$ from 100 to 25 km leads to a factor of 10 increase in response amplitudes. The topographic amplitude for $D = 25$ km is about 750 m (full orange line in Fig. 9a) and thus comparable to those observed on Enceladus or Titan. An increase in $\gamma$ from 0 to 0.1 mW/m$^3$ causes a dramatic decrease in the magnitude of response. The topographic amplitude obtained in this case for $D = 25$ km is less than 10 m (dashed orange line in Fig. 9a), suggesting that the effect of $\gamma$ on topography depends only weakly on the average thickness of the layer (cf. Fig. 6a).
<figure><img src="image_8.png" /><figcaption>Figure 6: Magnitude of surface (a, b) and ice-water interface (c, d) topography as a function of the grain size. The curves in panels a and c correspond to $\gamma = 0$ and different values of cut-off viscosity $\eta_{max}$, while panels b and d show the results obtained for $\eta_{max} = 10^{24}$ Pa s and different values of $\gamma$.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

Figure 7: Magnitude of stress at the surface ($\sigma_{xx}$) as a function of the grain size. The curves in panel a are computed for $\gamma = 0$ and $\eta_{max}$ ranging from $10^{21}$ to $10^{24}$ Pa s while panel b shows the results obtained for $\eta_{max} = 10^{24}$ Pa s and five different values of $\gamma$.

## 5. Degree of isostatic compensation

The gravity fields of Enceladus and Titan have only a small non-hydrostatic component (Iess et al., 2014; Durante et al., 2019), suggesting that the gravity anomalies associated with the long-wavelength topographic signal must be compensated by internal density anomalies. Since the ice shell is floating on a liquid water ocean, lateral variations in ice shell thickness provide a natural compensation mechanism. It has been shown that the degree of compensation depends not only on the thickness of the ice shell but also on its viscosity and the wavelength of the load (Čadek et al., 2017, 2019a; Kvorka et al., 2018; Akiba et al., 2022).

The current models of ice shells on Enceladus (e.g., Beuthe et al., 2016; Hemingway and Mittal, 2019; Čadek et al., 2016, 2019b), Dione (Beuthe et al., 2016; Zannoni et al., 2020) and Titan (e.g., Nimmo and Bills, 2010; Kvorka et al., 2018) are based on the assumption that for the ice shell to be in equilibrium, it is nec-
<figure><img src="image_10.png" /><figcaption>Figure 8: Steady-state flow field (color plots) and the rate of phase change (line graphs) obtained for different values of parameters *d* and *γ*. All models are computed for $\eta_{max} = 10^{24}$ Pa s. As in Fig. 2, the phase change velocity is plotted with the opposite sign.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Surface topography (a), ice-water boundary topography (b) and surface stress (c) obtained for average ice thickness *D* (Fig. 1) equal to 25 km (orange), 50 km (red) and 100 km (black) and $\gamma = 0$ (full lines) and 0.1 mW/m$^3$ (dashed lines). The simulations are carried out for $d = 10$ mm and $\eta_{max} = 10^{24}$ Pa s.</figcaption></figure>
essary for the degree of compensation,

$$C_{Airy}(x) = -\frac{h_s(x)\rho_i}{h_{i/w}(x)(\rho_w - \rho_i)}, \qquad (21)$$

to be equal to 1 (note that Eq. (21) is usually corrected for sphericity and the variations in gravity, see discussion in Čadek et al. (2019a)). It has been suggested that the shape of ice shells on Titan and Dione may follow the Pratt compensation scheme (Choukroun and Sotin, 2012; Čadek et al., 2021; Zannoni et al., 2020), but this mechanism is unlikely to explain the topography on Enceladus where the shell is only 20–25 km thick (Thomas et al., 2016).

The Airy isostatic concept, Eq. (21), assumes that the density of ice does not vary laterally and the degree of compensation depends only on the magnitude of the topographic loads. However, inspection of the models discussed in section 4 shows that even if the heat is transferred by conduction, the temperature (and hence density) of ice varies in the horizontal direction. The effect of density anomalies on the degree of compensation can be included in Eq. (21) as follows:

$$C_{corr}(x) = -\frac{h_s(x)\rho_i}{h_{i/w}(x)(\rho_w - \rho_i) - \int_0^D \bar{\rho}(y)\alpha(T(x,y)-\bar{T}(y))dy}, \qquad (22)$$

where $\bar{\rho}$ and $\bar{T}$ are the horizontally averaged density and temperature, respectively.

In Fig. 10a, $C_{Airy}$ and $C_{corr}$ are plotted as a function of the grain size for several values of the parameter $\gamma$. For a given value of $\gamma$, the degree of compensation $C_{corr}$ is always larger than $C_{Airy}$, suggesting that the use of the equation $C_{Airy} = 1$ in inverse modeling may lead to a slight (~5%) underestimation of ice thickness variations. The magnitude of $C_{corr}$ increases with increasing $\gamma$, approaching 1 for large $d$. This trend may be due to the fact that in Eq. (22), the effect of the topography and density anomalies is approximated by variations in pressure and the accuracy of this approximation increases as the magnitude of the topography and the deviatoric stress decrease (cf. Figs 6 and 7).
<figure><img src="image_12.png" /><figcaption>Figure 10: Degree of compensation, $C_{\text{Airy}}$ (dashed lines) and $C_{\text{corr}}$ (full lines), as a function of the grain size computed (a) for $\alpha(T)$ given by Eq. (18) and different values of parameter $\gamma$ and (b) for $\gamma = 0$ and different models of thermal expansion. All simulations were performed for $\eta_{\text{max}} = 10^{24}$ Pa s.</figcaption></figure>
Most of the previous studies of dynamic isostasy (e.g., Kvorka et al., 2018; Čadek et al., 2019a) have assumed that the flow of ice is driven by variations in ice shell thickness and the effect of thermal expansion on the density can be neglected. The contribution of the density anomalies to isostatic compensation is examined in Fig. 10b where $C_{Airy}$ and $C_{corr}$ are plotted as a function of the grain size for three different coefficients of thermal expansion: the temperature dependent $\alpha$ (Eq. 18), $\alpha = 1.6 \cdot 10^{-4} \text{ K}^{-1}$, corresponding to thermal expansion at the melting temperature, and $\alpha = 0$. The curves are shown only for $\gamma = 0$ but similar results were obtained for other values of $\gamma$. Inspection of the figure shows that $C_{Airy}$ varies by up to 15% depending on the model of $\alpha$ while the impact of thermal expansion on $C_{corr}$ is negligible (note that $C_{corr} = C_{Airy}$ if $\alpha = 0$). If the surface topography is known from observation, the ice-water interface topography can be calculated from Eq. (22):

$$h_{i/w} = -\frac{h_s}{C_{corr}} \frac{\rho_i}{(\rho_w - \rho_i)} + \frac{1}{\rho_w - \rho_i} \int_0^D \bar{\rho}\alpha(T - \bar{T})dy. \quad (23)$$

The relative magnitude of the thermal expansion term in Eq. (23) can be estimated by comparing the values of $C_{corr}$ and $C_{Airy}$ in Fig. 10. Since the lateral variations in temperature are negatively correlated with ice thickness (i.e., the two terms on the right-hand side have the same sign), the omission of thermal expansion in Eq. (23) leads to an underestimation of the amplitude of the ice-water topography by about 5% (compare the black lines in Fig. 10). On the other hand, when $\alpha = 1.6 \cdot 10^{-4} \text{ K}^{-1}$ and the thermal expansion term is included in Eq. (23), the value of $h_{i/w}$ is overestimated by 10–15%, as obvious when comparing $C_{Airy}$ computed for this value of $\alpha$ (dashed red line in Fig. 10b) with $C_{corr}$ computed for $\alpha(T)$ (full black line).
531

# 6. Comparison with the observation

532

533

534

535

536

537

538

539

540

541

542

543

544

There are only three icy moons where the long-wavelength non-hydrostatic topography is known with a reasonably high accuracy: Titan, Enceladus and Dione. The topography of these bodies has an amplitude between 300 m (Dione) and 1 km (Enceladus) and, in the case of Titan and Enceladus, is characterized by a strong degree-2 zonal component, suggesting a link between the topography and large-scale oceanic or atmospheric (on Titan) processes. Direct comparison of our results with the topographic heights on these moons is difficult because our simulations are performed in a planar geometry and the heat flux anomaly prescribed at the base of the ice layer has a wavelength of only 1600 km, which is about one tenth of Titan's circumference. This wavelength roughly corresponds to a spherical harmonic degree of 10, while the non-hydrostatic topography on Titan is dominated by zonal degree 2 and 4 harmonics and is only determined to degree and order 5.

545

546

547

548

549

550

551

552

553

554

555

The dependence of the topographic amplitude on the harmonic degree has been studied by Kvorka et al. (2018), who showed that the topography generated by a heat flux anomaly of degree 2 is about ten times larger than that obtained for degree 8. This suggests that the results shown in Fig. 6 could explain the long-wavelength topographic variations on Titan (about 500 m, see Durante et al. (2019)), provided that $\gamma$ is close to zero and the results can be generalized to spherical geometry. The large topographic variations on Enceladus, the circumference of which is ten times smaller than that of Titan, may be associated with the fact that the average thickness of the ice shell on Enceladus is only 20 – 25 km (Thomas et al., 2016) and the heat flux anomalies from the ocean are significantly larger than considered in our simulation (Čadek et al., 2019b).

556

557

However, as shown in Fig. 6b, large topographic amplitudes can only be obtained if $\gamma \approx 0$. As $\gamma$ increases, the amplitude of the surface topography decreases
and it becomes negligible for $\gamma \gtrsim 10^{-2}$ mW/m$^3$. Considering the large magnitude of the topography observed on Titan and Enceladus, this means that either (i) the real $\gamma$ is significantly smaller than $10^{-2}$ mW/m$^3$, or (ii) the processes described by parameters $q_o$ and $\gamma$ cannot be treated as independent of each other (implying that the heat flux from the ocean cannot be expressed in the form of Eq. (12)), or (iii) the concept of heat flux-induced long-wavelength topography proposed by Kvorka et al. (2018) is incorrect.

Alternative explanations for the long-wavelength topography on Titan have been proposed by, e.g., Choukroun and Sotin (2012) and Čadek et al. (2021). They include ethane precipitation in the polar regions, lateral variations in near-surface porosity and a variable thickness of the methane clathrate crust. On Europa, topography can also be affected by shallow melt pockets within the ice shell (Schmidt et al., 2011). On the other hand, the ocean heat flux scenario appears to be the most likely explanation for the large variations in ice shell thickness on Enceladus inferred from the analysis of Enceladus’ gravity field (e.g., Thomas et al., 2016; Van Hoolst et al., 2016; Beuthe et al., 2016; Čadek et al., 2016, 2019b; Hemingway et al., 2018).

The value of $\gamma$ that is relevant to subsurface oceans of icy moons is not known. If $\gamma$ represents the efficiency of the convective heat transfer in the ocean driven by variations in temperature along the phase boundary (Deguen et al., 2013; Deguen, 2013; Labrosse et al., 2018) then $\gamma \approx \rho_w c_p v_w dT_m/dy$. On Titan, $dT_m/dy \approx 1.4 \cdot 10^{-4}$ K/m, giving $\gamma \approx 6 \cdot 10^5 v_w$ (in mW/m$^3$). The results of our simulations show that even at very small flow rates ($v_w \gtrsim 10^{-8}$ m/s $\approx 0.3$ m/yr), the parameter $\gamma$ is large enough to erase the topography of the phase boundary. On the contrary, the large thickness variations inferred from topography and gravity data on Enceladus, Titan and Dione indicate that such oceanic flow along the phase boundary should be inhibited. This may be explained by the fact that the effect of temperature
variations on ocean circulation is relatively small and the flow in the real ocean
is mainly controlled by variations in salinity (e.g., Zeng and Jansen, 2021; Kang
et al., 2021). As argued by Zeng and Jansen (2021) in the case of Enceladus,
the temperature variation at the ice-ocean interface is about 0.1 K for thickness
variations of 20 km. This is significantly less than the temperature variation within
the ocean due to variations in salinity, which can exceed 1 K. Moreover, the phase
changes that occur at the upper boundary of a salty ocean may strongly affect
the salinity of water, possibly enhancing the heat flux in regions of melting, which
may subsequently lead to a positive feedback between heat flux and topography
(e.g., Lobo et al., 2021). The physical processes at the contact between ice and
water can also be influenced by other effects, such as the presence of clathrates
and organic compounds (Tobie et al., 2012; Sotin et al., 2021). The numerical
modeling performed in this study is a first step towards a consistent model of ice-
ocean interface evolution. Further modeling work is necessary to determine the
effect of salinity and other impurities on the coupling between the ice shell and
the ocean.

7. Conclusions

Based on our results we can conclude that the heat flux variations from the
ocean can produce a significant deformation of the surface only in the case when
the heat transfer through the ice shell occurs by conduction, the viscosity of ice
is relatively high and the parameter γ, characterizing the heat flux driven by
lateral variations in the melting temperature, is smaller than 10⁻³ mW/m³. The
amplitude of the topography increases with increasing viscosity of ice and, in the
specific case discussed in this paper, it can range from a few meters when the
viscosity is low and/or γ is large, to several hundred meters and more when the
viscosity at the base of the ice shell is higher than 10¹⁸ Pa s and γ ≲ 10⁻⁴ mW/m³.
611 When the heat transfer in the ice shell is dominated by convection, the variations
612 of heat flux from the deep ocean produce negligible topographic signal (~ 10 m).

613 A limitation of some previous studies which have investigated the effect of the
614 phase transition on the variations in ice shell thickness is that they have assumed
615 that in the absence of convection the ice shell behaves as a rigid body (i.e. η₀ → ∞
616 and v → 0). This has led to the wrong conclusion that the heat flux across the
617 ice-water boundary is continuous (cf. Equation (2.2) in Ojakangas and Stevenson
618 (1989)) and the advection term in the heat transfer equation can be neglected.
619 Our results show that even in the conductive case, the flow of ice must be taken
620 into account and the assumption of a rigid ice shell may result in a significant
621 overestimation of the magnitude of ice thickness variations.

622 If heat is transferred by conduction, the magnitude of the surface topography
623 strongly depends on the grain size and the additional heat flux due to lateral vari-
624 ation in the melting temperature at the base of the ice shell. In our model, the
625 grain size is treated as a model parameter and its value is kept constant during
626 the simulation. In reality, the grain size can vary in time and space, depending on
627 the amount of impurities (e.g., Durand et al., 2006) and the relative importance
628 of dislocation creep (e.g., Barr and McKinnon, 2007). Therefore it is likely that
629 the effect of large-scale heat flux variations from the ocean is further modified by
630 variations in the grain size, possibly leading to the development of topography
631 on smaller spatial scales and varying in time. Spatial variations of other material
632 parameters, such as the thermal conductivity, heat capacity and thermal expansiv-
633 ity, caused by variations in porosity of the near-surface material and the presence
634 of clathrates and organics in the ice shell, may also influence the topographic
635 response.

636 As in previous studies (Deguen et al., 2013; Deguen, 2013; Labrosse et al.,
637 2018; Morison et al., 2019; Agrusta et al., 2020), the additional heat flux due to
lateral variations in the melting temperature is assumed to be proportional to the topography of the phase boundary and is characterized by a single parameter (slope), here denoted by $\gamma$, which is in turn proportional to the speed of water flow. We find that even small values of $\gamma$, corresponding to a flow speed of $\sim 0.1$ m/yr, lead to a drastic (order of magnitude) reduction of the topography compared to the case where the additional heat flux is neglected. Although the concept of the heat flux controlled by variations in the melting temperature seems to reflect the physical conditions at the deformed phase boundary, the models with $\gamma \gtrsim 10^{-2}$ mW/m$^3$ can neither account for the large variations in ice shell thickness on Enceladus confirmed by analysis of gravity and topography data, nor provide a plausible explanation of the degree 2 topography on Titan. Even if the heat flux variations from the deep ocean are large as proposed by recent models of ocean dynamics, their effect is almost completely compensated by the flow of water due to temperature variations at the phase boundary. This indicates that either the description of the ocean heat flux used in our model (Eqs. (12)–(14)) is too simplistic, or the effect of variations in the melting temperature is limited due to the presence of salts and other compounds.

Compared to the effect of parameter $\gamma$, which controls the variations in ice thickness in the conductive model, the role of thermal expansion in the large-scale deformation of the ice shell is minor and can be neglected in the first approximation. However, thermally induced variations in the density are likely to be important for the solution of the inverse problem when the shape of the ice-water boundary and the average thickness of the ice shell are inferred from the observed shape and gravity. The results presented in section 5 can contribute to the ongoing discussion of isostatic compensation on icy moons (e.g., Beuthe et al., 2016; Hemingway and Matsuyama, 2017; Čadek et al., 2017, 2019a; Akiba et al., 2022). Underlying this discussion is the fact that the outer shell of an icy moon differs
from the crust of a terrestrial planet in that a first-order phase transition occurs at the bottom boundary. The stability of the ice shell thus requires an exchange of mass and energy between the ice shell and the ocean, implying that the isostatic compensation has a dynamic nature. We demonstrate that the flow which occurs in the conductive ice shell is driven not only by the variations in pressure at the ice-water boundary but also by the large-scale variations in the density of ice.

The physical processes explored in the present study, involving dynamic cou-
pling between the ice shell and the subsurface ocean, will be tested by future mea-
surements to be performed by Europa Clipper and JUICE on Europa, Ganymede
and Callisto. Gravimetric and altimetric measurements will be performed with
unprecedented accuracy and spatial resolution by the gravity experiment 3GM
(Cappuccio et al., 2020) and the laser altimeter GALA (Hussmann et al., 2019)
onboard JUICE at Ganymede and to a lesser extent by the gravity experiment and
ice-penetrating radar in altimetric mode (REASON) onboard Clipper at Europa
(Mazarico et al., 2021). Similar measurements will be performed at Callisto but
the coverage will be reduced due to a limited number of flybys. Similarly to what
has been achieved by Cassini on Titan and Enceladus, the precise determination
of long-wavelength topography and gravity anomalies will provide key constraints
on ice-water interface deflection. Precise magnetic measurements may also reveal
large-scale oceanic flows (Vance et al., 2021), providing constraints on the dynam-
ical regime of the ocean and potentially highlighting lateral variations in salinity
and correlation with ice thickness variations.

**Acknowledgements** We thank Stéphane Labrosse and an anonymous reviewer for helpful comments that improved the manuscript. This research was supported by the Czech Science Foundation through project No. 22-20388S. M.K. acknowledges the support from the Charles University project SVV-2021-260581. K.K.
691 was supported by the Charles University Research program No. UNCE/SCI/023.

692 G.C. and G.T. acknowledge the support from the ANR COLOSSE project and

693 from CNES in the preparation of the ESA JUICE and NASA Europa Clipper

694 missions.
695 References

696 Agrusta, R., Morison, A., Labrosse, S., Deguen, R., Alboussière, T., Tackley, P. J., Dubuffet, F., 2020. Mantle convection interacting with magma oceans.
697
698 Geophys. J. Int. 220, 1878–1892.

699 Akiba, R., Ermakov, A., Militzer, B., 2022. Probing the icy shell structure of ocean
700 worlds with gravity-topography admittance. Planet. Sci. J. 3, in press.

701 Alnaes, M. S., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M. E., Wells, G. N., 2015. The FEniCS Project Version 1.5. Archive of Numerical Software 3 (100).
702
703

704 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.
705
706

707 Askenazy, Y., Sayag, R., Tziperman, E., 2018. Dynamics of the global meridional
708 ice flow of Europa's icy shell. Nat. Astr. 2, 43–49.

709 Barr, A. C., McKinnon, W. B., 2007. Convection in ice I shells and mantles with self-consistent grain size. J. Geophys. Res. 112, 2012.
710

711 Barr, A. C., Pappalardo, R. T., 2005. Onset of convection in the icy galilean
712 satellites: Influence of rheology. J. Geophys. Res.-Planets 110, E12005.

713 Beuthe, M., Rivoldini, A., Trinh, A., 2016. Enceladus's and Dione's floating ice
714 shells supported by minimum stress isostasy. Geophys. Res. Lett. 43, 10,088–
715 10,096.

716 Buffo, J., Meyer, C., Parkinson, J., 2021. Dynamics of a solidifying icy aatellite
717 shell. J. Geophys. Res.-Planets 126, e2020JE006741.
718 Buffo, J., Schmidt, B., Huber, C., Walker, C., 2020. Entrainment and dynamics
719 of ocean-derived impurities within Europa's ice shell. J. Geophys. Res.-Planets
720 125, e2020JE006394.

721 Běhounková, M., Tobie, G., Choblet, G., Čadek, O., 2013. Impact of tidal heating
722 on the onset of convection in Enceladus's ice shell. Icarus 226, 898-904.

723 Běhounková, M., Tobie, G., Choblet, G., Kervazo, M., Daswani, M. M., Dumoulin,
724 C., Vance, S. D., 2021. Tidally induced magmatic pulses on the oceanic floor of
725 Jupiter's moon Europa. Geophys. Res. Lett. 48, e2020GL090077.

726 Čadek, O., Běhounková, M., Tobie, G., Choblet, G., 2017. Viscoelastic relaxation
727 of Enceladus's ice shell. Icarus 291, 31-35.

728 Čadek, O., Kalousová, K., Kvorka, J., Sotin, C., 2021. The density structure of
729 Titan's outer ice shell. Icarus 364, 114466.

730 Čadek, O., Souček, O., Běhounková, M., 2019a. Is Airy isostasy applicable to icy
731 moons? Geophys. Res. Lett. 46, 14,299-14,306.

732 Čadek, O., Souček, O., Běhounková, M., Choblet, G., Tobie, G., Hron, J., 2019b.
733 Long-term stability of Enceladus' uneven ice shell. Icarus 319, 476-484.

734 Čadek, O., Tobie, G., Van Hoolst, T., Massé, M., Choblet, G., Lefèvre, A., Mitri,
735 G., Baland, R.-M., Běhounková, M., Bourgeois, O., Trinh, A., 2016. Enceladus's
736 internal ocean and ice shell constrained from Cassini gravity, shape and libration
737 data. Geophys. Res. Lett. 43, 5653-5660.

738 Cappuccio, P., Hickey, A., Durante, D., Di Benedetto, M., Iess, L., De Marchi,
739 F., Plainaki, C., Milillo, A., Mura, A., 2020. Ganymede's gravity, tides and
740 rotational state from JUICE's 3GM experiment simulation. Planet. Space Sci.
741 187, 104902.
742 Choblet, G., Tobie, G., Sotin, C., Kalousová, K., Grasset, O., 2017a. Heat transport in the high-pressure ice mantle of large icy moons. Icarus 285, 252–262.

744 Choblet, G., Tobie, G., Sotin, C., Běhounková, M., Čadek, O., Postberg, F., Souček, O., 2017b. Powering prolonged hydrothermal activity inside Enceladus. Nat. Astr. 1, 841–847.

747 Choukroun, M., Sotin, C., 2012. Is Titan's shape caused by its meteorology and carbon cycle? Geophys. Res. Lett. 39, L04201.

749 Deguen, R., 2013. Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries. J. Earth Sci. 24, 669–682.

751 Deguen, R., Alboussière, T., Cardin, P., 2013. Thermal convection in Earth's inner core with phase change at its boundary. Geophys. J. Int. 194, 1310–1334.

753 Dingemans, M. W., 1997. Water Wave Propagation over Uneven Bottom, Part 1 – Linear Wave Propagation. World Scientific, Singapore, etc.

755 Donea, J., Huerta, A., Ponthot, J.-P., Rodríguez-Ferran, A., 2004. Arbitrary Lagrangian-Eulerian methods. In: Encyclopedia of Computational Mechanics. Vol. 1. John Wiley & Sons, Ltd, Ch. 14, pp. 413–437.

758 Durand, G., Weiss, J., Lipenkov, V., Barnola, J. M., Krinner, G., Parrenin, F., Delmonte, B., Ritz, C., Duval, P., Rothlisberger, R., Bigler, M., 2006. Effect of impurities on grain growth in cold ice sheets. J. Geophys. Res.-Earth Surface 111, F01015.

762 Durante, D., Hemingway, D. J., Racioppa, P., Iess, L., Stevenson, D. J., 2019. Titan's gravity field and interior structure after Cassini. Icarus 326, 123–132.
764 Durham, W. B., Stern, L. A., 2001. Rheological properties of water ice – Applications to satellites of the outer planets. Annu. Rev. Earth Planet. Sci. 29, 295–330.

767 Feistel, R., Wagner, W., 2006. A new equation of state for H₂O ice Ih. J. Phys. Chem. Res. Data 35, 1021–1047.

769 Freeman, J., Moresi, L., May, D. A., 2006. Thermal convection with a water ice I rheology: Implications for icy satellite evolution. Icarus 180, 251–264.

771 Giese, B., Wagner, R., Hussmann, H., Neukum, G., Perry, J., Helfenstein, P., Thomas, P. C., 2008. Enceladus: An estimate of heat flux and lithospheric thickness from flexurally supported topography. Geophys. Res. Lett. 35, L24204.

774 Goldsby, D. L., Kohlstedt, D. L., 2001. Superplastic deformation of ice: Experimental observations. J. Geophys. Res.-Solid Earth 106, 11017–11030.

776 Goodman, J. C., Collins, G. C., Marshall, J., T., P. R., 2004. Hydrothermal plume dynamics on Europa: Implications for chaos formation. J. Geophys. Res.-Planets 109, E03008.

779 Goodman, J. C., Lenferink, E., 2012. Numerical simulations of marine hydrothermal plumes for Europa and other icy worlds. Icarus 221, 970–983.

781 Green, A. P., Montesi, L. G. J., Cooper, C. M., 2021. The growth of Europa’s icy shell: Convection and crystallization. J. Geophys. Res.-Planets 126, e2020JE006677.

784 Han, L., Showman, A., 2005. Thermo-compositional convection in Europa’s icy shell with salinity. Geophys. Res. Lett. 32, L20201.
786 Harel, L., Dumoulin, C., Choblet, G., Tobie, G., Besserer, J., 2020. Scaling of heat transfer in stagnant lid convection for the outer ice shells of icy moons: influence of rheology. Icarus 338, 113448.

787 transfer in stagnant lid convection for the outer ice shells of icy moons: influence
788 of rheology. Icarus 338, 113448.

789 Hemingway, D. J., Iess, L., Tajeddine, R., Tobie, G., 2018. The interior of Enceladus. In: Schenk, P. M. et al. (Ed.), Enceladus and the Icy Moons of Saturn. Univ. of Arizona, Tucson, pp. 57–77.

790
791
Univ. of Arizona, Tucson, pp. 57–77.

792 Hemingway, D. J., Matsuyama, I., 2017. Isostatic equilibrium in spherical coordinates and implications for crustal thickness on the Moon, Mars, Enceladus, and elsewhere. Geophys. Res. Lett. 44, 7695–7705.

793
794
elsewhere. Geophys. Res. Lett. 44, 7695–7705.

795 Hemingway, D. J., Mittal, T., 2019. Enceladus's ice shell structure as a window
796 on internal heat production. Icarus, in press.

796
on internal heat production. Icarus, in press.

797 Hobbs, P. V., 1974. Ice Physics. Oxford Univ. Press.

798
799
800
801
802
803
804
805
806
Hussmann, H., Lingenauber, K., Kallenbach, R., Enya, K., Thomas, N., Lara, L., Althaus, C., Araki, H., Behnke, T., Castro-Marin, J., Eisenmenger, H., Gerber, T., de la Revilla, M., Huttig, C., Ishibashi, K., Jiménez-Ortega, J., Kimura, J., Kobayashi, M., Lötzke, H., Lichopoj, A., Lüdicke, F., Martínez-Navajas, I., Michaelis, H., Namiki, N., Noda, H., Oberst, J., Oshigami, S., Garcia, J., Rodrigo, J., Rösner, K., Stark, A., Steinbrügge, G., Thabaut, P., del Togno, S., Touhara, K., Villamil, S., Wendler, B., Wickhusen, K., Willner, K., 2019. The Ganymede laser altimeter (GALA): key objectives, instrument design, and performance. CEAS Space J. 11, 381–390.

799
Althaus, C., Araki, H., Behnke, T., Castro-Marin, J., Eisenmenger, H., Gerber,

800
T., de la Revilla, M., Huttig, C., Ishibashi, K., Jiménez-Ortega, J., Kimura,
801
J., Kobayashi, M., Lötzke, H., Lichopoj, A., Lüdicke, F., Martínez-Navajas,

802
I., Michaelis, H., Namiki, N., Noda, H., Oberst, J., Oshigami, S., Garcia, J.,
803
Rodrigo, J., Rösner, K., Stark, A., Steinbrügge, G., Thabaut, P., del Togno,
804
S., Touhara, K., Villamil, S., Wendler, B., Wickhusen, K., Willner, K., 2019.

805
806
The Ganymede laser altimeter (GALA): key objectives, instrument design, and
performance. CEAS Space J. 11, 381–390.

807 Iess, L., Stevenson, D. J., Parisi, M., Hemingway, D., Jacobson, R. A., Lunine, J. I., Nimmo, F., Armstrong, J. W., Asmar, S. W. Ducci, M., Tortora, P., 2014. The gravity field and interior structure of Enceladus. Science 344, 78–80.

808
809
The gravity field and interior structure of Enceladus. Science 344, 78–80.
810 Kalousová, K., Sotin, C., 2020. The insulating effect of methane clathrate crust on Titan's thermal evolution. Geophys. Res. Lett. 47, e2020GL087481.

812 Kalousová, K., Sotin, C., Choblet, G., Tobie, G., Grasset, O., 2018. Two-phase convection in Ganymede's high-pressure ice layer – Implications for its geological evolution. Icarus 299, 133–147.

815 Kang, W., Mittal, T., Bire, S., Campin, J., Marshall, J., 2021. How does salinity shape ocean circulation and ice geometry on Enceladus and other icy satellites? Sci. Adv. 8, eabm4665.

818 Kvorka, J., Čadek, O., 2022. A numerical model of convective heat transfer in Titan's subsurface ocean. Icarus 376, 114853.

820 Kvorka, J., Čadek, O., Tobie, G., Choblet, G., 2018. Does Titan's long-wavelength topography contain information about subsurface ocean circulation dynamics? Icarus 310, 149–164.

823 Labrosse, S., Morison, A., Deguen, R., Alboussière, T., 2018. Rayleigh-Benard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries. J. Fluid Mech. 846, 5–36.

826 Lobo, A. H., Thompson, A. F., Vance, S. D., Tharimena, S., 2021. A pole-to-equator ocean overturning circulation on enceladus. Nat. Geosci. 14, 185–189.

828 Logg, A., Mardal, K. A., Wells, G. N., 2012. Automated Solution of Differential Equations by the Finite Element Method. Springer, Berlin, Heidelberg, New York.

831 Marcolli, C., 2017. Ice nucleation triggered by negative pressure. Sci. Rep. 7, 16634.
832 Mazarico, E., Buccino, D., Castillo-Rogez, J., Dombard, A., Genova, A., Hussmann, H., Kiefer, W., Lunine, J., McKinnon, W., Nimmo, F., Park, R., Tortora, P., Withers, P., Roberts, J., Korth, H., Senske, D., Pappalardo, R., 2021. The Europa Clipper gravity/radio science investigation. In: 52nd Lunar and Planetary Science Conference, held virtually, 15-19 March, 2021, Contribution No. 2548, id.1784. Pasadena, CA.

838 Mitri, G., Showman, A. P., 2005. Convective-conductive transitions and sensitivity of a convecting ice shell to perturbations in heat flux and tidal-heating rate: Implications for Europa. Icarus 177, 447-460.

841 Morison, A., Labrosse, S., Deguen, R., Alboussière, T., 2019. Timescale of overturn in a magma ocean cumulate. Earth Planet. Sci. Lett. 516, 25-36.

843 Nimmo, F., Bills, B. G., 2010. Shell thickness variations and the long-wavelength topography of Titan. Icarus 208, 896-904.

845 Nimmo, F., Pappalardo, R. T., 2016. Ocean worlds in the outer solar system. J. Geophys. Res.-Planets 121, 1378-1399.

847 Nimmo, F., Schenk, P., 2006. Normal faulting on Europa: implications for ice shell properties. J. Struct. Geol. 26, 2194-2203.

849 Nimmo, F., Thomas, P. C., Pappalardo, R. T., Moore, W. B., 2007. The global shape of Europa: Constraints on lateral shell thickness variations. Icarus 191, 183-192.

852 Ojakangas, G. W., Stevenson, D. J., 1989. Thermal state of an ice shell on Europa. Icarus 81, 220-241.

854 Peddinti, D. A., McNamara, A. K., 2019. Dynamical investigation of a thickening ice-shell: Implications for the icy moon Europa. Icarus 329, 251-269.
856 Radebaugh, J., Lorenz, R. D., Kirk, R., Lunine, J. I., Stofan, E., Lopes, R., Wall, 857 S., 2007. Mountains on Titan observed by Cassini Radar. Icarus 192, 77–91.

858 Röttger, K., Endriss, A., Ihringer, J., Doyle, S., Kuhs, W. F., 1994. Lattice constants and thermal expansion of H₂O and D₂O ice Ih between 10 and 265 K. 859 Acta Cryst. 50, 644–648.

861 Schmidt, B. E., Blankenship, D. D., Patterson, G. W., Schenk, P. M., 2011. Active 862 formation of 'chaos terrain' over shallow subsurface water on Europa. Nature 863 479, 502–505.

864 Showman, A. P., Han, L. J., 2004. Numerical simulations of convection in Eu- 865 ropa's ice shell: Implications for surface features. J. Geophys. Res.-Planets 109, 866 E01010.

867 Slattery, J. C., Sagis, L., Oh, E. S., 2007. Interfacial Transport Phenomena, 2nd 868 ed. Springer, New York.

869 Soderlund, K. M., 2019. Ocean dynamics of outer Solar System satellites. Geophys. 870 Res. Lett. 46, 8700–8710.

871 Soderlund, K. M., Schmidt, B. E., Wicht, J., Blankenship, D. D., 2014. Ocean-872 driven heating of Europa's icy shell at low latitudes. Nat. Geosci. 7, 16–19.

873 Sotin, C., Kalousová, K., Tobie, G., 2021. Titan's interior structure and dynamics 874 after the Cassini-Huygens mission. Ann. Rev. Earth Planet. Sci. 49, 579–607.

875 Souchez, R., Jean-Baptiste, P., Petit, J., Lipenkov, V., Jouzel, J., 2003. What is 876 the deepest part of the Vostok ice core telling us? Earth Sci. Rev. 60, 131–146.

877 Stevenson, D. J., 2000. Limits on the variations of thickness of Europa's ice shell 878 (abstract). Lunar Planet. Sci. Conf. 31, 1506.
879 Terra-Nova, F., Amit, H., Choblet, G., Tobie, G., Bouffard, M., Čadek, O., 2022. The influence of heterogeneous seafloor heat flux on the cooling patterns of Ganymede's and Titan's subsurface oceans. Icarus, in press.

882 Thomas, P. C., Tajeddine, R., Tiscareno, M. S., Burns, J. A., Joseph, J., Loredo, T. J., Helfenstein, P., Porco, C., 2016. Enceladus's measured physical libration requires a global subsurface ocean. Icarus 264, 37–47.

885 Thomson, R. E., Delaney, J. R., 2001. Evidence for a weakly stratified European ocean sustained by seafloor heat flux. J. Geophys. Res.-Planets 106, 12,355–12,365.

888 Tobie, G., Choblet, G., Sotin, C., 2003. Tidally heated convection: Constraints on Europa's ice shell thickness. J. Geophys. Res.-Planets 108, 5124.

890 Tobie, G., Gautier, D., Hersant, F., 2012. Titan's bulk composition constrained by Cassini-Huygens: Implication for internal outgassing. Astrophys. J. 752, 125.

892 Van Hoolst, T., Baland, R. M., Trinh, A., 2016. The diurnal libration and interior structure of Enceladus. Icarus 277, 311–318.

894 Vance, S., Styczinski, M., Bills, B., Cochrane, C., Soderlund, K., GómezPérez, N., Paty, C., 2021. Magnetic induction responses of Jupiter's ocean moons including effects from adiabatic convection. J. Geophys. Res.-Planets 126, e2020JE006418.

897 Zannoni, M., Hemingway, D., Casajus, L. G., Tortora, P., 2020. The gravity field and interior structure of Dione. Icarus 345, 113713.

899 Zeng, Y., Jansen, M. F., 2021. Ocean circulation on Enceladus with a high- versus low-salinity ocean. Planet. Sci. J. 2, 151.