The thermo-mechanical evolution of Titan's ice shell is primarily controlled by the mode of the heat transfer in the ice shell and the amount of heat coming from the ocean. Numerical models of thermal convection in Titan's ocean have suggested that the heat flux from the core can be strongly affected by the Coriolis effect, leading to a latitudinal gradient in the heat flux at the base of the ice shell [1,2,3]. These heat flux variations are of the order of mW/m2 and may induce phase changes at the ice/ocean interface, affecting the thickness of the ice shell. Until now, the effect of lateral variations in the heat flux has been studied under the assumption that the ice shell is highly viscous and the heat transfer in the ice occurs by conduction. In such a case, the heat flux from the ocean is negatively correlated with the ice shell thickness . A higher-than-average heat flux in polar regions then causes melting of ice and thinning of the ice shell, possibly explaining the fact that the poles on Titan are about 300 m lower than the equator .
The issue of whether Titan's ice shell is convecting or not has not yet been resolved. The mode of heat transfer depends on the ice shell thickness and the viscosity of ice, which in turn depends on the grain size and the temperature of the ocean. Since these parameters are not accurately known, the convection model cannot be excluded from analysis. The response of the convecting ice shell to heat flux variations from the ocean is likely to be different from that described for the conduction model . When heat is transferred by convection, the ice shell is on average warmer and less viscous than in the conductive case. The loss/gain of ice mass at the ice/ocean interface caused by melting/freezing is compensated by flow of low-viscosity ice and, therefore, the shape of the ice/ocean interface remains unchanged. The shape of the surface topography is difficult to assess a priori because it may depend on the vigor of convection.
In this study, we investigate the effect of spatial variations in the heat flux from the ocean on the behavior of a convecting ice layer. The models discussed below have been obtained by simultaneously solving the momentum equation, the transport equation for temperature and the mass conservation equation for an incompressible fluid. The temperature is fixed to 90 K on the upper boundary and to 265 K on the bottom (ice/ocean) boundary. The computational domain is 100 km thick and has an aspect ratio of 16:1. We assume that the deformation of ice is controlled by diffusion creep and the viscosity of ice only depends on temperature and grain size.
The surface is treated as a material boundary while the ice/ocean interface is open and its shape can change due to phase changes. The heat flux variations from the ocean enter the problem through the energy conservation law which relates the jump in the heat flux across the boundary with the relative motion of the ice/ocean interface [3,6,7]. Note that the heat flux across the phase boundary is generally not continuous (as assumed, for example, by ), but depends on the velocity of ice flow, which cannot be neglected if the viscosity of ice is low (< 1016 Pa.s). The lateral variation of the heat flux imposed at the bottom boundary of the domain is scaled to correspond to 50 % of the average heat flux (Fig. 1a).
In Figs. 1 and 2, we compare the results obtained for a grain size of 2 mm and 3 mm. Although these values are close to each other, the models represent two end-member cases discussed above. While in the former case (2 mm), the heat transfer is controlled by convection, in the latter case (3 mm), convection does not develop and the heat is transferred by conduction.
Inspection of Figs. 1 and 2 shows that the convection model gives the opposite sign of the surface topography than the conduction model (Fig. 1b) and predicts a significant (8 K) increase in ice temperature above the positive heat flux anomaly (Fig. 2a). The conduction and convection models also give different orientations of tectonic stresses - extension in the case of convection and compression in the case of conduction (Fig. 2b). Both models predict only small amplitudes of the geoid (Fig. 2c), suggesting that the equilibrium state can be approximated by either Airy or Pratt isostatic model, depending on the mode of heat transfer.
Variations in the heat flux from the ocean can significantly influence Titan’s topography and the distribution of stress near the surface. The topography is negatively correlated with the heat flux when the heat in the ice shell is transferred by conduction, whereas a positive correlation and more complex topography response is found for the convective heat transfer. In contrast to the topography, the heat flux variations have only a small effect on the geoid. A careful analysis of the geoid and topography data collected by Cassini could thus provides an insight into the processes occurring in Titan's interior and helps answer the question of which heat transfer mode dominates in Titan's ice shell.
This research was supported by the Czech Science Foundation through project No. 19-10809S. M.K. acknowledges the support from the Charles University project SVV-2020-260581. K.K. was supported by the Charles University Research program No. UNCE/SCI/023. G.C. and G.T. acknowledge the support from the ANR COLOSSe project.
 Soderlund, 2019, Geophys. Res. Lett. 46, 8700–8710.
 Amit et al., 2020, Icarus 338, 113509.
 Kvorka et al., 2018, Icarus 310, 149–164.
 Corlies et al., 2017, Geophys. Res. Lett. 44, 11754–11761.
 Nimmo and Bills, 2010, Icarus 201, 896–904.
 Čadek et al., 2019a, Icarus 319, 476–484.
 Čadek et al., 2019b, Geophys. Res. Lett. 46, 14299–14306.
 Ojakangas and Stevenson, 1989, Icarus 81, 242–270.