Filipe Terra-Nova
email: filipeterranova@gmail.com

Hagay Amit

Gaël Choblet

Gabriel Tobie

Mathieu Bouffard

Ondřej Čadek

The influence of heterogeneous seafloor heat flux on the cooling patterns of Ganymede's and Titan's subsurface oceans

Keywords:

Several icy moons of Jupiter and Saturn are known to possess deep water oceans. Heating in the rocky mantle underneath often produces heterogeneous heat flux patterns at the ocean's seafloor. How this internal ocean dynamically relates the seafloor to the surface ice shell is a crucial question to understand the long term evolution of icy moons. Here we investigate how a heterogeneous seafloor heat flux pattern affects the convection and heat transfer in the subsurface ocean of large icy worlds involving a high pressure ice layer beneath the seafloor such as Titan or Ganymede. We perform rotating convection simulations in a thin 3D spherical shell with a prescribed heterogeneous bottom heat flux inferred from 3D convection simulations of the underlying mantle (Choblet et al., 2017b). In our simulations, although the amplitude of imposed inner boundary heat flux heterogeneity is rather moderate, preferred longitudes of intense outer boundary heat flux are highly correlated with longitudes of intense inner boundary heat flux. In addition, a small imposed inner boundary large-scale order 2 pattern is amplified at the outer boundary heat flux by the convection in the thin shell. Lastly, deviations from axisymmetry and equatorial symmetry in the outer boundary heat flux increase with the main convection vigor and the amplitude of the inner boundary heterogeneity. In our models polar vs. equatorial cooling is mostly controlled by inertial effects, as was found by Amit et al. (2020) for homogeneous boundary conditions, with the latitudinally equilibrated inner boundary heterogeneity acting to reduce the amplitude of this effect. Our results support polar cooling for Titan's sub-surface ocean.

Introduction

The exploration of Jupiter's and Saturn's systems, by the Galileo and Cassini-Huygens missions respectively, has revealed that several moons harbor salty water oceans underneath their icy crusts (see [START_REF] Hussmann | Interiors and evolution of icy satellites[END_REF], for a review). Earlier predictions based on thermal evolution models (e.g. [START_REF] Kirk | Thermal evolution of a differentiated Ganymede and implications for surface features[END_REF][START_REF] Grasset | The cooling rate of a liquid shell in Titan's interior[END_REF][START_REF] Tobie | Titan's internal structure inferred from a coupled thermal-orbital model[END_REF] have thus been confirmed by various geophysical detection techniques. Magnetic induction caused by the moons' orbits in a tilted Jovian magnetosphere enabled the identification of global scale conductive layers in Europa, Ganymede and Callisto [START_REF] Khurana | Induced magnetic fields as evidence for subsurface oceans in Europa and Callisto[END_REF][START_REF] Khurana | Searching for liquid water in Europa by using surface observatories[END_REF][START_REF] Zimmer | Subsurface oceans on Europa and Callisto: Constraints from Galileo magnetometer observations[END_REF][START_REF] Kivelson | The permanent and inductive magnetic moments of Ganymede[END_REF]. Owing to the presence of an internal dynamo, the more ambiguous measurement in the case of Ganymede was later confirmed by Earth-based observation of auroral oscillations [START_REF] Saur | The search for a subsurface ocean in Ganymede with hubble space telescope observations of its auroral ovals[END_REF]. In the system of Saturn, low frequency electric waves observed in Titan's atmosphere by the Huygens probe first hinted at a deep conducting interface beneath the ice (e.g. [START_REF] Béghin | Analytic theory of Titan's Schumann resonance: Constraints on ionospheric conductivity and buried water ocean[END_REF] although this interpretation has been recently challenged [START_REF] Lorenz | Schumann resonance on titan : Huygens observations critically re-assessed and prospects for the dragonfly mission[END_REF]. The presence of an ocean was indicated by the observation of an abnormally high obliquity [START_REF] Baland | Titan's obliquity as evidence of a subsurface ocean?[END_REF] and further confirmed by the detection of internal gravitational tides [START_REF] Iess | The tides of Titan[END_REF], indicating the existence of a liquid layer at depth decoupling the outer shell from the deep solid interior.

No direct observations are available to constrain the dynamics of internal oceans, however long-wavelength topography and gravity data can be used indirectly to shed light on heat transfer conveyed by https://doi.org/10.1016/j.icarus.2022.115232 Received 19 March 2022; Received in revised form 22 July 2022; Accepted 16 August 2022 ocean dynamics. The combination of these two data sets provided constraints on the ice shell structure (average thickness and lateral variations, see Čadek et al. (2016, 2019) and [START_REF] Beuthe | Crustal control of dissipative ocean tides in Enceladus and other icy moons[END_REF] for Enceladus and [START_REF] Lefevre | Structure and dynamics of titan's outer icy shell constrained from cassini data[END_REF] and [START_REF] Kvorka | Does Titan's long-wavelength topography contain information about subsurface ocean dynamics?[END_REF] for Titan) which in turn can be interpreted in terms of melting/freezing at the ice/ocean interface caused by heterogeneity in the heat supplied by the buried oceans [START_REF] Kvorka | Does Titan's long-wavelength topography contain information about subsurface ocean dynamics?[END_REF][START_REF] Čadek | Long-term stability of Enceladus uneven ice shell[END_REF].

The global scale dynamics of convecting buried oceans is a major question for the long-term evolution, the surface properties and the habitability of icy moons, yet answers are only preliminary. In a pioneering study, [START_REF] Soderlund | Ocean-driven heating of Europa's icy shell at low latitudes[END_REF] found latitudinal dependence of heat transfer which, in the case of Europa, may explain the preferred occurrence of specific surface features (chaos terrains, salt deposits) at mid latitudes. This latitudinal modulation of the outer boundary heat flux has been described by [START_REF] Yadav | Effect of shear and magnetic field on the heat-transfer efficiency of convection in rotating spherical shells[END_REF] in the more general case of rotating convection. [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF] investigated thoroughly the various convective regimes simulated in a rotating spherical shell as functions of the vigor of convection and the rotation speed. In this framework, considering thinner shells appropriate for the geometry of ocean worlds, [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] further showed a transition in cooling patterns: from equatorial cooling close to the rapidly-rotating regime (a feature analyzed by [START_REF] Miquel | Equatorially trapped convection in a rapidly rotating shallow shell[END_REF], for such shallow shells) to polar cooling (and eventually, no latitudinal dependence) closer to the non-rotating regime.

In the case of the largest icy moons such as Ganymede and Titan, radioactive decay in the rocky core alone is sufficient to maintain the presence of a buried ocean, possibly favored by anti-freeze compounds. The phase diagram of water (e.g Choukroun andGrasset, 2007, 2010) indicates that the ocean is separated from the rocky core by a thick high-pressure ice mantle, initially envisioned as a barrier to chemical exchanges. However, the detection of a significant amount of radiogenic 40 Ar in Titan's atmosphere implies a relatively efficient chemical transfer throughout the thick hydrosphere to the surface [START_REF] Tobie | Episodic outgassing as the origin of atmospheric methane on Titan[END_REF][START_REF] Tobie | Titan's bulk composition constrained by Cassini-Huygens: implication for internal outgassing[END_REF]. Indeed, numerical simulations of solid-state convection in the high-pressure ice mantle showed that heat transfer is controlled by melt production at depth and extraction of meltwater to the ocean (Choblet et al., 2017b;Kalousová et al., 2018;Kalousová and Sotin, 2018), providing an efficient way to transport chemical species from the rocky core to the ocean. The heat flux pattern through the seafloor is therefore shaped by the outlets of heat pipes of hot, partially molten ice. Interestingly for its possible influence on the melting/freezing regime of the ocean ceiling (upper ocean-ice interface), this pattern is expected to be almost permanent over tens of million years: while velocities in the high-pressure ice mantle may typically amount to a few meters per year, hot convective instabilities do not exhibit significant amounts of lateral motion over the course of a simulation (i.e. with a typical duration of 10-50 Myr).

Here, we aim to investigate how this heterogeneous seafloor heat flux pattern affects the ocean's convective flow. In order to provide a wider understanding of the phenomena at play, we consider simulations that cover different cooling patterns which may correspond to different regimes of rotating convection [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF]. The prescribed heat flux at the bottom of the ocean reproduces the heterogeneous pattern exhibited at the outer boundary of 3D convection simulations in the high-pressure ice mantle (Choblet et al., 2017b). Our moderate amplitudes of heat flux heterogeneity are also chosen to reflect these results, i.e. much more modest than predicted for Enceladus, for example (Choblet et al., 2017a).

More specifically, we focus on the following questions: Is the latitudinal dependence of the heat flux at the top of the ocean affected by seafloor heterogeneity? Is there a correlation between the heat flux at the top of the ocean and the imposed seafloor heterogeneity?

In Section 2 we describe the numerical model for ocean convective flow. The results are presented in Section 3 and discussed in Section 4. We list our main conclusions in Section 5.

Methods

To model thermal convection in the internal oceans of Ganymede and Titan, we solve the set of hydrodynamics equations for the conservation of momentum, energy and mass for a rotating spherical shell in the Boussinesq approximation. These write, in non-dimensional form (e.g. [START_REF] Wicht | Inner-core conductivity in numerical dynamo simulations[END_REF],

𝐸 ( 𝜕𝐮 𝜕𝑡 + 𝐮 ⋅ ∇𝐮 -∇ 2 𝐮 ) + 2 ẑ × 𝐮 + ∇𝑃 = 𝑅𝑎 * 𝐫 𝑟 𝑜 𝑇 , ( 1
)
𝜕𝑇 𝜕𝑡 + 𝐮 ⋅ ∇𝑇 = 1 𝑃 𝑟 ∇ 2 𝑇 , ( 2
)
∇ ⋅ 𝐮 = 0, ( 3
)
where u is the fluid velocity, ẑ the direction of the rotation axis, 𝑃 the pressure, r the radial position vector, 𝑟 𝑜 the outer boundary radius and 𝑇 the temperature. Although the linear radial gravity profile in (1) may not be the most appropriate for buried oceans of icy moons, within a thin shell the gravity varies very little for any plausible radial profile (e.g. [START_REF] Vance | Geophysical investigations of habitability in ice-covered ocean worlds[END_REF]. Here we choose this gravity profile to ease the comparison with the results of [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF].

Eqs. ( 1)-( 2) are governed by three (internal) control parameters. The Ekman number represents the ratio of viscous to Coriolis forces:

𝐸 = 𝜈 𝛺𝐷 2 . ( 4
)
The Prandtl number is a ratio of diffusivities:

𝑃 𝑟 = 𝜈 𝜅 . ( 5
)
The modified Rayleigh number 𝑅𝑎 * represents the ratio of buoyancy vs. retarding forces. It is related to the conventional Rayleigh number 𝑅𝑎 by 𝑅𝑎 * = 𝑅𝑎 ⋅ 𝐸∕𝑃 𝑟, where the conventional heat flux based Rayleigh number is given by

𝑅𝑎 = 𝛼𝑔 𝑜 𝑞 0 𝐷 4 𝜅𝜈𝑘 . ( 6
)
In ( 4)-( 6), 𝛺 is the rotation rate, 𝜈 the kinematic viscosity, 𝐷 the shell thickness, 𝑔 𝑜 the gravitational acceleration at the outer boundary, 𝑞 0 the mean inner boundary heat flux, 𝑘 the thermal conductivity and 𝜅 the thermal diffusivity. Non-dimensional time is in units of viscous diffusion time 𝐷 2 ∕𝜈. In our models the heat flux is imposed from below, thus 𝑞 0 is the mean inner boundary heat flux. In the following, we explore control parameters (Table 1) that produce various cooling patterns.

In all simulations the inner to outer boundary radii ratio is set to 0.8, corresponding to a relatively thin shell compared e.g. to Earth's outer core aspect ratio of 0.35 (see Table 1 of [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] and only slightly thicker compared to the aspect ratio estimates of 0.84-0.96 for the subsurface oceans of Titan [START_REF] Vance | Geophysical investigations of habitability in ice-covered ocean worlds[END_REF] and Ganymede (Jara-Orué and Vermeersen, 2016). On the inner boundary, we impose a heat flux pattern (Fig. 1) based on output from simulations of high-pressure ice applied to Titan's mantle (Choblet et al., 2017b). As this interface corresponds to a phase change, prescribing a constant (melting) temperature may seem a more natural choice; in our simple framework, however, this boundary is impermeable while the most appropriate description is that of a permeable boundary (see Deguen et al., 2013, for example). As a consequence, prescribing a heterogeneous heat flux provides a way to describe the thermal influence of the high-pressure ice mantle. In our simulations, we ensure a posteriori that temperature variations along the bottom boundary are negligible. Note also that while we use the heat flux output from the solid high-pressure ice mantle as the value entering the bottom of the ocean, in reality the heat flux is expected to be discontinuous across the phase interface owing to the absorption or release of latent heat [START_REF] Kvorka | Does Titan's long-wavelength topography contain information about subsurface ocean dynamics?[END_REF][START_REF] Čadek | Long-term stability of Enceladus uneven ice shell[END_REF]. We posit that even accounting for this heat flux jump would not modify the pattern of heterogeneity which is precisely is the level of convective supercriticality for the corresponding cases with homogeneous boundary conditions, 𝑞 * 𝑖 is the amplitude of inner boundary heat flux heterogeneity, 𝓁 𝑚𝑎𝑥 is the spherical harmonic truncation degree and order and 𝑛 𝑟 is the number of grid points in the radial direction. The number of grid points in the longitudinal and latitudinal directions 𝑛 𝜙 and 𝑛 𝜃 in our setup are related to 𝓁 𝑚𝑎𝑥 by 𝑛 𝜙 = 3𝓁 𝑚𝑎𝑥 and 𝑛 𝜃 = 𝑛 𝜙 ∕2. 𝑅𝑒 is the Reynolds number (10), 𝑅𝑜 the Rossby number (11), 𝑅𝑜 𝑐 the convective Rossby number (12) and 𝑅𝑎 𝑅𝑎 𝑤 𝑇 a ratio that represents competition between thermal and rotational boundary layers ( 14).  is the duration of the simulation in units of advection times 𝐷∕𝑈 , where 𝐷 is the shell thickness and 𝑈 the RMS velocity amplitude in the shell. Fig. 1. Imposed pattern of inner boundary heat flux taken from the output of a high-pressure ice convection simulation applied to Titan's mantle (Choblet et al., 2017b). The heat flux is expanded up to spherical harmonic degree and order 30. what we consider here. The amplitude of the inner boundary heat flux lateral heterogeneity is quantified by [START_REF] Olson | The time averaged magnetic field in numerical dynamos with nonuniform boundary heat flow[END_REF]:

𝑞 * 𝑖 = 𝑞 𝑚𝑎𝑥 -𝑞 𝑚𝑖𝑛 2𝑞 0 , ( 7
)
where 𝑞 𝑚𝑎𝑥 and 𝑞 𝑚𝑖𝑛 are the maximum and minimum of the inner boundary heat flux, respectively. In this paper, 𝑞 * 𝑖 (where 𝑖 denotes inner boundary) is an additional (boundary) control parameter. The amplitude of this heat flux heterogeneity 𝑞 * 𝑖 is limited to moderate values (i.e. order unity) between 0.57 -2.68 (see Table 1), in agreement with model predictions for the heat flux lateral variability at Titan's seafloor (Choblet et al., 2017b). Fixed temperature is set on the outer boundary, so that the heat flux there is an output of the models. On both boundaries the mechanical boundary conditions are no-slip.

Note that the imposed inner boundary heat flux (Fig. 1) corresponds to one representative snapshot from a single model of high-pressure ice applied to Titan's mantle (case 4 𝑚 𝐷 of Choblet et al., 2017b). Obviously other high-pressure ice models with different combinations of control parameters may yield different patterns. However, our goal is to explore the extent to which the ocean dynamics alters a given bottom boundary heat flux. In that respect the pattern in Fig. 1 adequately serves our purposes and different patterns likely would not alter the main conclusions.

The convective supercriticality 𝑅𝑎 𝑇 ∕𝑅𝑎 𝑐 is calculated for corresponding cases with homogeneous inner boundary heat flux (𝑞 * 𝑖 = 0), where 𝑅𝑎 𝑐 is the critical Rayleigh number for the onset of convection and 𝑅𝑎 𝑇 is a temperature Rayleigh number given by [START_REF] Soderlund | Ocean dynamics of outer solar system satellites[END_REF])

𝑅𝑎 𝑇 = 𝑅𝑎 𝑁𝑢 , ( 8
)
with 𝑁𝑢 the Nusselt number

𝑁𝑢 = 𝑞 0 𝐷 𝑘𝛥𝑇 . ( 9
)
This definition of the Nusselt number ignores the effect of spherical geometry (e.g. Nu is not 1 for a conductive layer). The convective supercriticality values are given in Table 1.

Several diagnostic parameters are reported in Table 1. The ratio between inertial and viscous forces is the Reynolds number

𝑅𝑒 = 𝑈 𝐷 𝜈 , ( 10
)
where 𝑈 is the RMS velocity. The ratio between inertial and Coriolis forces is given by the Rossby number

𝑅𝑜 = 𝑈 𝛺𝐷 = 𝑅𝑒𝐸. ( 11
)
Following [START_REF] Soderlund | Ocean-driven heating of Europa's icy shell at low latitudes[END_REF] we also consider an alternative measure of inertia given by the convective Rossby number

𝑅𝑜 𝑐 = √ 𝑅𝑎𝐸 2 𝑃 𝑟 . ( 12
)
Additionally we define a horizontal convective Rossby number

𝑅𝑜 𝑐𝐻 = √ 𝑅𝑎 𝐻 𝐸 2 𝑃 𝑟 , ( 13
)
where 𝑅𝑎 𝐻 = 𝑅𝑎𝑞 * 𝑖 (where the subscript 'H' stands for 'horizontal') reflects the convective effect of the lateral heterogeneity. As in [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] we also monitor a ratio representing a competition between thermal and rotational boundary layers [START_REF] King | Boundary layer control of rotating convection systems[END_REF])

𝑅𝑎∕𝑅𝑎 𝑤 𝑇 = 0.1𝑅𝑎𝐸 3∕2 (14)
and focusing on the inner boundary heat flux anomaly we define

𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 = 0.1𝑅𝑎 𝐻 𝐸 3∕2 . ( 15
)
So far we presented parameters that characterize the global dynamics. Now we introduce parameters that characterize the pattern of outer boundary cooling. To assess whether polar or equatorial cooling governs, we follow [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] and calculate the ratio between the integrated heat flux inside and outside the tangent cylinder (TC) at the outer boundary 𝑞 ℎ∕𝑙 𝑜 (where 'h' and 'l' superscripts stand for high and low latitudes respectively) as follows. The co-latitude where the tangent cylinder [START_REF] Aurnou | Experiments on convection in Earth's core tangent cylinder[END_REF] intersects the outer boundary 𝜃 𝑡𝑐 is given by sin 𝜃 𝑡𝑐 = 𝑟 𝑖 ∕𝑟 𝑜 . For 𝑟 𝑖 ∕𝑟 𝑜 = 0.8, the TC co-latitude is 𝜃 𝑡𝑐 ≈ 53 • . The integrated heat flux inside the TC is given by

𝑞 ℎ 𝑜 = 1 𝑆 ℎ ( ∫ 2𝜋 0 ∫ 𝜃 𝑡𝑐 0 𝑞 𝑜 𝑑𝑆 + ∫ 2𝜋 0 ∫ 𝜋 𝜋-𝜃 𝑡𝑐 𝑞 𝑜 𝑑𝑆 ) (16)
and outside the TC by

𝑞 𝑙 𝑜 = 1 𝑆 𝑙 ∫ 2𝜋 0 ∫ 𝜋-𝜃 𝑡𝑐 𝜃 𝑡𝑐 𝑞 𝑜 𝑑𝑆, ( 17
)
where 𝑆 ℎ and 𝑆 𝑙 denote the outer boundary spherical surface areas inside and outside the TC, respectively. Polar vs. equatorial cooling is defined by

𝑞 ℎ∕𝑙 𝑜 = 𝑞 ℎ 𝑜 -𝑞 𝑙 𝑜 𝑞 ℎ 𝑜 + 𝑞 𝑙 𝑜 . ( 18
)
Note that here we use the full outer boundary heat flux 𝑞 𝑜 , not its anomaly 𝛿𝑞 𝑜 , to avoid division by zero. Also note that each integral is normalized by the appropriate area so we do not calculate total fluxes but mean fluxes (or flux density). The expression (18) has some desired properties. For polar cooling it is positive, for equatorial cooling it is negative. If one of the two cooling types is much stronger it will approach plus/minus 1, whereas if the two are comparable it will approach zero. We avoid quantifying the heterogeneous inner boundary impact on the cooling pattern by 𝑞 * 𝑜 , i.e. a quantity similar to (7) applied to the outer boundary heat flux, because unlike 𝑞 * 𝑖 , which relies on a fixed imposed pattern, 𝑞 * 𝑜 relies on local output values hence in practice it is biased by small-scale variability which might arise solely due to the finite simulation time of the runs. [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] overcame this problem by computing 𝑞 * 𝑜 from latitudinal averages. However, here we cannot adopt this approach because the heterogeneous inner boundary is expected to yield longitude dependence which we wish to quantify. Instead, we decompose the outer boundary heat flux into zonal equatorially symmetric (denoted by 'zs') and non zonal equatorially symmetric (denoted by 'nzs') parts. The zonal equatorially symmetric outer boundary heat flux is given by

𝑞 𝑧𝑠 𝑜 (𝜃) = ( 𝑞 𝑧 𝑜 (𝜃) + 𝑞 𝑧 𝑜 (𝜋 -𝜃) ) ∕2, (19)
where the 'z' superscript denotes zonal average, while the non zonal equatorially symmetric part is simply

𝑞 𝑛𝑧𝑠 𝑜 (𝜙, 𝜃) = 𝑞 𝑜 (𝜙, 𝜃) -𝑞 𝑧𝑠 𝑜 (𝜃), (20)
where 𝜙 and 𝜃 in ( 19) and ( 20) are longitude and co-latitude, respectively. We compute the RMS values of the non zonal equatorially 2018)). Note that the definition of the various numbers differ slightly for 𝐸 and 𝑅𝑎, explaining the different values. In particular, the heat flux value is used in the present study to compute Ra: it is assumed to result only from radiogenic decay in the rocky interior, assumed chondritic, at present -see values in [START_REF] Hussmann | Interiors and evolution of icy satellites[END_REF]. The plausible range for 𝑞 * 𝑖 is essentially unconstrained by observations: we use values obtained by models simulating convection in high-pressure ices (Choblet et al., 2017b). symmetric parts of the outer boundary heat flux anomaly 𝑅𝑀𝑆(𝛿𝑞 𝑛𝑧𝑠 𝑜 ) which is interpreted in terms of boundary control because without any boundary heterogeneity no persistent long-term deviations from axisymmetry or equatorial symmetry are expected. This is evident in e.g. Fig. 7 of [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] where the long-term time-average zonal outer boundary heat flux profiles exhibit high level of equatorial symmetry even inside the tangent cylinder.

The impact of the imposed inner boundary heterogeneity on the pattern of outer boundary heat flux is further evaluated by calculating the correlation and the maximal cross-correlation between the averages of the inner and outer boundary heat flux patterns along latitude lines  𝜙 and 𝑚𝑎𝑥( 𝑐𝑟 𝜙 ), respectively. Here we focus on longitudinal dependence because latitudinal variability in the outer boundary heat flux emerges even when the inner boundary heat flux pattern is homogeneous [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] whereas long-term longitudinal variability is expected only in the presence of boundary heterogeneity. While the correlation 𝐶 𝜙 quantifies the trend agreement between the inner and outer boundary heat flux patterns, the maximal cross-correlation 𝑚𝑎𝑥(𝐶 𝑐𝑟 𝜙 ) accounts for possible horizontal shifts between the imposed heat flux pattern at the bottom of the shell and the resulting heat flux at the top. Analogous shifts were reported in dynamo simulations between the imposed outer boundary heat flux pattern and the resulting time-average convection pattern [START_REF] Olson | The time averaged magnetic field in numerical dynamos with nonuniform boundary heat flow[END_REF][START_REF] Aubert | Detecting thermal boundary control in surface flows from numerical dynamos[END_REF][START_REF] Takahashi | Effects of thermally heterogeneous structure in the lowermost mantle on the geomagnetic field strength[END_REF].

The quantities that evaluate the boundary effect in Table 3 were computed based on the long-term time-average morphologies of the dynamical models. The averaging time of ∼77-10 000 advection times is over the full simulation run (see  in Table 1). For comparison, [START_REF] Mound | Regional stratification at the top of Earth's core due to core-mantle boundary heat flux variations[END_REF] reported an average over 37 advection times. The smallscale inner boundary heat flux pattern that we imposed (Fig. 1) likely converges slower than the commonly large-scale pattern imposed on the outer boundary in studies of Earth's core (e.g. [START_REF] Mound | Regional stratification at the top of Earth's core due to core-mantle boundary heat flux variations[END_REF]. Nevertheless, convergence tests using prolonged runs of cases with distinct 𝐸 values and corresponding spatial resolutions indicate that the main boundary heterogeneity scale features exhibit little variations.

Table 2 indicates plausible ranges for the dimensionless parameters that characterize the subsurface oceans of Titan and Ganymede. We note that considerable uncertainties remain on most physical properties of ocean water, not only because of an uncertain composition. In addition, the ocean thickness 𝐷, or the mean basal heat flux 𝑞 0 are unconstrained.

Results

In our hydrodynamic simulations we explored the ranges 𝐸 = 1 ⋅ 10 -3 -5 ⋅ 10 -5 and 𝑅𝑎 = 2 ⋅ 10 5 -8 ⋅ 10 7 , while 𝑃 𝑟 was fixed to 1. These 𝐸 values are much larger than expected for planetary conditions and likewise the 𝑅𝑎 values are too low, as in most numerical simulations of convection in deep liquid layers (e.g. [START_REF] Glatzmaier | Geodynamo simulations: how realistic are they?[END_REF]. Our analysis includes representative snapshots as well as time averages. To assess dynamical regimes of simulations with non-realistic parameters, Amit et al. ( 2020) localized their models on the regime diagram of Gastine is the same quantity applied to the inner boundary using the same latitude limits. 𝑅𝑀𝑆(𝑞 𝑛𝑧𝑠 𝑜 ) is the integrated non zonal symmetric heat flux (20) which represents the impact of the boundary heterogeneity.  𝜙 and 𝑚𝑎𝑥( 𝑐𝑟 𝜙 ) are the correlation and the maximum of the cross-correlation coefficients, respectively, between inner and outer boundary heat flux averaged along longitude lines.  2 ∕ 2 and  10 ∕ 10 are the outer/inner boundary ratio of non-dimensional heat flux power for order 2 and 10, respectively. (2016). However, such regime diagrams depend on 𝑅𝑎 and 𝐸 exclusively [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF][START_REF] Long | Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions[END_REF], while the little explored impact of 𝑞 * 𝑖 (which may vary depending on the particular imposed heat flux pattern) is absent in these diagrams. Since the main control parameters are roughly similar in our study and that of [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] with homogeneous boundary conditions, we compare our findings with those of [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] and assess the impact of the boundary heterogeneity in Section 5.

Polar vs. equatorial cooling

Competing effects of rotation and convection control the largescale pattern of outer boundary heat flux (see Amit et al., 2020, and references therein). Larger role of rotation effects leads to larger heat flux at low latitudes, whereas polar cooling emerges when convection is increased. Our outer boundary heat flux patterns encompass these two endmember dynamical regimes. Figs. 2 and3 show outer boundary heat flux anomaly distributions of our hydrodynamic models with equatorial and polar cooling, respectively.

The equatorial cooling patterns seen in Fig. 2 are characterized by low-and mid-latitude strips oriented in the north-south direction. Such invariance in the direction of the rotation axis indicates dominance of rotational effects (e.g. [START_REF] Busse | Thermal instabilities in rapidly rotating systems[END_REF]. These strips are not evenly distributed in longitude, especially in cases with stronger amplitude of imposed boundary heterogeneity or/and stronger main convection vigor. For the snapshots in Fig. 2, this longitudinal preference points to a possible consequence of either statistical preference due to boundary control or time-dependent chaotic convection. The latitude extent of the strips roughly reaches the edge of the TC (latitudes 37 • in our models' setup), though in some longitudes the meridional heat flux strips penetrate the TC. In models with increasing effects of rapid rotation (i.e. smaller Rossby numbers) these strips become increasingly thinner (compare Figs. 2a andb vs. 2c and d). Increasing amplitude of imposed inner boundary heterogeneity leads to destruction of some north-south strips (compare Figs. 2c andd vs. 2e and f). Some weaker heat flux anomalies inside the TC nevertheless appear in both snapshots and time averages of cases with equatorial cooling. On time averages, pronounced negative heat flux anomaly belts appear at the edges of the TC. As 𝑞 * 𝑖 increases from case 12 to 14, positive heat flux anomalies become more pronounced, especially inside the TC (compare Fig. 2d andf).

The polar cooling patterns seen in Fig. 3 are characterized by stronger high-latitude positive heat flux anomalies. The time-averages show more clearly that the high-latitude positive heat flux anomalies exceed the equatorial positive anomalies (also see positive 𝑞 ℎ∕𝑙 𝑜 for these cases in Table 3). Our cases 5-9 with the highest values of 𝑅𝑎 𝐻 are characterized by polar cooling and lack any coherent north-south strips in their outer boundary heat flux patterns. Instead, small-scale circular features characterize the instantaneous outer boundary heat flux. These high-latitude patches roughly extend until the edge of the TC.

To further illustrate the differences between polar and equatorial cooling, Fig. 4 compares polar and equatorial views of the outer boundary heat flux anomalies for case 12 (Fig. 4a andb) characterized by equatorial cooling vs. case 7 (Fig. 4c andd) characterized by polar cooling. As mentioned above, for models with stronger convection characterized by polar cooling the low-latitude north-south strips fade (Fig. 4c) and small-scale heat flux structures emerge at high latitudes inside the tangent cylinder (Fig. 4d). When rotation effects dominate, high-latitude positive anomalies are weaker and the heat flux is governed by north-south elongated strips outside the TC.

The latitude dependence of the imposed inner boundary heat flux is shown in Fig. 5 (blue curve). This latitudinal dependence is an outcome of a solid-state convection simulation in the high-pressure ice mantle, which includes no rotational (hence latitudinal) effect. In addition, the heat distribution is homogeneous in this model as tidal heating in the high-pressure ices and the rocky core is expected to be a very minor contribution and neglected in the simulation. Therefore there is a priori no reason for the imposed inner boundary heat flux pattern to contain a marked latitude dependence other than fortuitous. As its power spectrum peaks at relatively large spherical harmonic degrees (𝓁 ≃ 20; Choblet et al., 2017b), the average along longitude lines is expected to display modest largest-scale variability. Yet, somewhat incidentally it contains a rather substantial level of equatorial symmetry (red curve in Fig. 5). In contrast, as expected the integrated heat flux inside and outside the TC are rather comparable, albeit slightly more polar (see positive weak 𝑞 ℎ∕𝑙 𝑖 values in Table 3). Fig. 6 presents the latitude dependence of the time-averaged outer boundary heat flux patterns. The models are classified into two groups depending on the dominating cooling pattern (see Table 3). Cases dominated by equatorial cooling (negative 𝑞 ℎ∕𝑙 𝑜 in Table 3) are shown in Fig. 6a, whereas cases dominated by polar cooling (positive 𝑞 ℎ∕𝑙 𝑜 in Table 3) are shown in Fig. 6b. Fig. 6a highlights three characteristics of our equatorial cooling cases: positive anomalies at low latitudes with peaks at the equator, negative anomalies at high-latitudes with lows at the poles and pronounced negative heat flux anomalies at the edges of the TC. A somewhat simpler pattern of positive/negative anomalies inside/outside the TC respectively appears in the polar cooling cases (Fig. 6b).

The dynamical origin of the latitudinal distribution of outer boundary heat flux is presented in Fig. 7. In the equatorial cooling case, the flow consists of thin axial columns outside the TC which are invariant along the direction parallel to the rotation axis (Fig. 7a andc), typical to rapidly rotating flows (e.g. [START_REF] Busse | Thermal instabilities in rapidly rotating systems[END_REF]Olson et al., 1999;[START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF][START_REF] Long | Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions[END_REF]. In contrast, the radial velocity in a snapshot of the polar cooling case exhibits a 3D morphology (Fig. 7b andd) characteristic of inertial effects driven by vigorous convection. The zonal velocity profiles emphasize that indeed the convective activity in the polar cooling case is concentrated in the polar regions (Fig. 7d) whereas in the equatorial cooling case the axial columns efficiently transfer heat from the inner to the outer boundary at low and mid latitudes (Fig. 7c). • east of the central meridian inside the TC. This is best seen on long-term time averages (especially in Figs. 2f and3f). Indeed such a projection of the inner boundary heat flux pattern onto the outer boundary is expected for a thin shell, large amplitude of imposed inner boundary heterogeneity and larger role of convection vs. rotation in the dynamics. Stronger rotational effects are expected to smear latitudinal variability outside the TC but would still favor preferred longitudes where axial columns would expel more heat flux through the outer boundary (Fig. 2d andf). There is a hint of such a relation even when comparing the latitude dependence, for example the asymmetry between the positive/negative inner boundary heat flux anomaly in the North/South pole respectively (Fig. 5) appears also in most models (especially in the polar cooling cases) on the outer boundary (Fig. 6b).

Impact of boundary heterogeneity

We now focus on the longitudinal distributions of the time-averaged outer boundary heat flux which are expected to be uniform in the absence of boundary heterogeneity [START_REF] Olson | The time averaged magnetic field in numerical dynamos with nonuniform boundary heat flow[END_REF][START_REF] Amit | Towards more realistic core-mantle boundary heat flux patterns: a source of diversity in planetary dynamos[END_REF]. We compare these model outputs with the imposed inner boundary heat flux averaged along longitude lines (Fig. 8).

Fig. 9 presents the outer boundary heat flux averaged along longitude lines of the hydrodynamic models. Here we classified the models according to their correlation coefficient  𝜙 . Common to all cases is the dominance of the short wavelength that characterizes the inner boundary heat flux pattern (Fig. 8). The correlation coefficient seems to decrease with decreasing 𝐸 (from Fig. 9a to 9c), though this might be due to shorter averaging times for our lower 𝐸 cases (Table 1). In addition, outer boundary heat flux patterns of some models contain a large-scale order 2 signal (Fig. 9b). This order 2 pattern is particularly evident in the polar regions of polar cooling cases 5-9 (see e.g. Figs. 3d, f and4d) and also (to a lesser extent) at low latitudes of the equatorial cooling cases (e.g. Fig. 2f). Indeed, the outer boundary heat flux power of order 2 is amplified with respect to the inner boundary heat flux of the same order in almost all cases (Table 3), with the largest amplification in the polar cooling cases (Fig. 10). For comparison, the peak power of inner boundary heat flux at order 10 actually becomes weaker at the outer boundary for all cases (Table 3). We elaborate on the origin of the large-scale order 2 pattern in the Discussion.

For a more quantitative analysis we calculated the correlation and the maximum cross-correlation coefficients  𝜙 and 𝑚𝑎𝑥( 𝑐𝑟 𝜙 ), respectively, between the longitude-averaged outer vs. inner boundary heat flux (Table 3). The 𝑚𝑎𝑥( 𝑐𝑟 𝜙 ) coefficients account for possible azimuthal shifts between the inner and the outer boundary heat flux patterns. In all cases the correlation is positive and statistically significant. Based on the student's 𝑡-test [START_REF] Press | Numerical Recipes in FORTRAN: The Art of Scientific Computing[END_REF], for the spatial resolution of our models (Table 1), a level of 95% statistical significance corresponds to a global spatial correlation of about 0.2 [START_REF] Rau | Core flow inversion tested with numerical dynamo models[END_REF]. In all cases except case 12  𝜙 exceeds 0.2, while 𝑚𝑎𝑥( 𝑐𝑟 𝜙 ) exceeds 0.2 in all cases, indicating clear boundary control on the outer boundary heat flux. The largest values of  𝜙 and 𝑚𝑎𝑥( 𝑐𝑟 𝜙 ) are found for polar cooling case 5 and the smallest for equatorial cooling case 12. In general, models with polar cooling patterns show larger correlation improvement from  𝜙 to 𝑚𝑎𝑥( 𝑐𝑟 𝜙 ). We note that the low-latitude positive heat flux anomalies are somewhat larger in the northern hemisphere, as can be seen in the latitudinal profiles (Fig. 6). This is due to the slightly larger imposed inner boundary heat flux at low-latitudes of the northern hemisphere (Fig. 5). Obviously this north-south asymmetry is just pure chancethere is no physical reason for any hemispherical preference in the pattern of the inner boundary heat flux. We merely point here that this slight imposed hemispherical asymmetry is reflected in the resulting outer boundary heat flux.

Dependence on control parameters

It is important to note that the cooling type (equatorial or polar) is not affected by the amplitude of the inner boundary heterogeneity. Our models exhibit the same type of cooling patterns for different 𝑞 * 𝑖 . However, the amplitude of the cooling pattern is strongly affected by the boundary heterogeneity. In addition, the dependence of the cooling pattern on the internal control parameters (i.e. 𝐸 and 𝑅𝑎) is also modified by the boundary heterogeneity. We now turn to examine the dependence of our measure of polar vs. equatorial cooling 𝑞 ℎ∕𝑙 𝑜 (18) on the level of inertia in the system [START_REF] Soderlund | Ocean-driven heating of Europa's icy shell at low latitudes[END_REF][START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF]. We test the dependence of this quantity on the ratios representing competitions between thermal and rotational boundary layers 𝑅𝑎 𝑇 ∕𝑅𝑎 𝑤 𝑇 (as in ( 14) but based on the temperature Rayleigh number) and 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 (15), the latter capturing the convective effect of the inner boundary heat flux heterogeneity. The results are displayed in Fig. 11. In our models 𝑞 ℎ∕𝑙 𝑜 depends not only on the competition between rotation and convection [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] but also on the amplitude of the inner boundary heterogeneity 𝑞 * 𝑖 . Nevertheless, Fig. 11a shows that equatorial/polar cooling patterns characterized by negative/positive 𝑞 ℎ∕𝑙 𝑜 correspond to lower/higher 𝑅𝑎 𝑇 ∕𝑅𝑎 𝑤 𝑇 , respectively, in qualitative agreement with the main results of [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] for homogeneous boundary conditions (see red symbols in Fig. 11a), i.e. stronger inertial effects lead to polar cooling whereas rotational effects promote equatorial cooling. Also in agreement with [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF], our cases with 𝑅𝑎 𝑇 ∕𝑅𝑎 𝑤 𝑇 < 1 are classified as equatorial cooling, whereas our cases with 𝑅𝑎 𝑇 ∕𝑅𝑎 𝑤 𝑇 > 1 are characterized by polar cooling. However, our absolute 𝑞 ℎ∕𝑙 𝑜 values are significantly smaller than those obtained by [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF]. In addition, our 𝑞 ℎ∕𝑙 𝑜 values for the equatorial cooling cases slightly decrease with increasing 𝑅𝑎 𝑇 ∕𝑅𝑎 𝑤 𝑇 (Fig. 11a). These discrepancies are inevitably attributable to the imposed heterogeneous inner boundary heat flux in our simulations, but also to the lower Rossby numbers models that we studied compared to the models of [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] and [START_REF] Kvorka | A numerical model of convective heat transfer in Titan's subsurface ocean[END_REF]. However, accounting for 𝑞 * 𝑖 practically does not modify the qualitative dependence of 𝑞 ℎ∕𝑙 𝑜 on the measure of inertia (Fig. 11b).

Next in Fig. 12, we explore the dependence of the amplitude of the non zonal equatorially symmetric part of the outer boundary heat flux anomaly on the control parameters. This component of the heat flux is expected to arise solely due to the imposed inner boundary heterogeneity. Therefore, a natural candidate for the amplitude of the deviations of the outer boundary heat flux from axisymmetry and equatorial symmetry is the horizontal Rayleigh number 𝑅𝑎 𝐻 (e.g. [START_REF] Willis | Thermal core-mantle interaction: Exploring regimes for 'locked' dynamo action[END_REF][START_REF] Monteux | Giant impacts, heterogeneous mantle heating and a past hemispheric dynamo on mars[END_REF]. Fig. 12 indeed suggests that Fig. 10. Non-dimensional time-average azimuthal heat flux power of outer (colors) and inner (thick gray) boundary heat flux for four selected models, vs. order. The largest (blue and orange) and the smallest (red and green) amplifications of order 2 among the studied cases are shown. increasing 𝑅𝑎 𝐻 results in larger deviations from axisymmetry and equatorial symmetry caused by the imposed inner boundary heterogeneity. Attempts to explain 𝛿𝑞 𝑛𝑧𝑠 𝑜 with 𝑅𝑎 or 𝑞 * 𝑖 are much less satisfactory, and the same holds with 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 (not shown). For a single spherical harmonic imposed boundary condition, the order of the time-average convection pattern follows that of the imposed boundary condition, though a phase shift is possible [START_REF] Olson | The time averaged magnetic field in numerical dynamos with nonuniform boundary heat flow[END_REF][START_REF] Aubert | Detecting thermal boundary control in surface flows from numerical dynamos[END_REF][START_REF] Takahashi | Effects of thermally heterogeneous structure in the lowermost mantle on the geomagnetic field strength[END_REF]. When multi-harmonic boundary conditions are used, as in our case, convective mixing may lead to a more complex pattern distortion. Analogous Fig. 11, in Fig. 13 we explore the dependence of the correlations between the patterns of outer and inner boundary heat flux  𝜙 and 𝑚𝑎𝑥( 𝑐𝑟 𝜙 ) on 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 . Two branches are seen in Fig. 13. One branch associated with equatorial cooling cases and polar cooling cases with relatively low Rossby numbers exhibits large reduction in correlation for small decrease in 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 , whereas another branch associated entirely with polar cooling cases exhibits very weak reduction in correlation for changes in 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 . In addition, in models occupying the latter branch a more significant improvement in correlation appears when allowing for azimuthal shifts (compare Fig. 13a andb). Below we interpret these trends in terms of different dynamical regimes.

Discussion

We classified the latitudinal profiles of outer boundary heat flux into two types of cooling patterns: equatorial and polar (Fig. 6 and Table 3). [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] found that large inertia favors polar cooling whereas weak inertia leads to equatorial cooling. We find that 𝑞 ℎ∕𝑙 𝑜 , our measure of polar vs. equatorial cooling, depends on inertia in a similar way. Our polar cooling cases are systematically characterized by larger effective Rossby numbers than our equatorial cooling cases. However, the absolute 𝑞 ℎ∕𝑙 𝑜 values in our models are significantly lower than those obtained by [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] with homogeneous boundary conditions and similar effective Rossby numbers (Table 3 and Fig. 11a). Because the imposed inner boundary heterogeneity has comparable heat flux partitioning at the latitude ranges that correspond to inside and outside the TC intersection with the outer boundary (Fig. 5) as expected from solid convection simulations with no preferred latitude, it tends to reduce the latitudinal variability of the heat flux. Therefore, the resulting cooling pattern is less variable in latitude than that obtained with homogeneous boundary conditions (Fig. 11a). The extent of this reduction depends on the amount of zonal symmetric heat flux. The pattern imposed in this study contains a significant amount of zonal equatorial symmetry (Fig. 5). If projected to the outer boundary it would produce a close to balanced heat flux inside and outside the TC. Competing effects of boundary heterogeneity and internal control parameters hence lead to considerable reduction in the absolute values of 𝑞 ℎ∕𝑙 𝑜 (Fig. 11a). This trend persists even when incorporating 𝑞 * 𝑖 in the scaling law using 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 (Fig. 11b). The decrease in 𝑞 ℎ∕𝑙 𝑜 with increasing 𝑅𝑎 𝑇 ∕𝑅𝑎 𝑤 𝑇 in our equatorial cooling cases arises from the low Rossby numbers of these models which reside in the weakly non-linear parametric regime [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF] which was not explored by [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] and [START_REF] Kvorka | A numerical model of convective heat transfer in Titan's subsurface ocean[END_REF]. In the weakly non-linear regime, decreasing Rossby numbers approach the onset of convection and uniform latitudinal dependence of heat flux. In Fig. 11b there is a change in trend within the equatorial cooling cases at 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 ∼ 0.3. For larger 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 values 𝑞 ℎ∕𝑙 𝑜 increases with increasing 𝑅𝑎 𝑇 ∕𝑅𝑎 𝑤 𝑇 values, as in the homogeneous boundary condition models [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF][START_REF] Kvorka | A numerical model of convective heat transfer in Titan's subsurface ocean[END_REF]. This value of 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 ∼ 0.3 probably marks the limit between the weakly non-linear and transitional regimes for the pattern of inner boundary heat flux used in our study.

In all our models we set 𝑃 𝑟 = 1, which corresponds to a turbulent Prandtl number. This choice is motivated for comparison purposes with numerous previous studies (e.g. [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF][START_REF] Soderlund | Ocean dynamics of outer solar system satellites[END_REF][START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF]. However, in subsurface icy moons 𝑃 𝑟 = 10. [START_REF] Kvorka | A numerical model of convective heat transfer in Titan's subsurface ocean[END_REF] explored the dependence of the heat transfer in rotating convecting thin spherical shells on 𝑃 𝑟. They found for all 𝑃 𝑟 values that increasing effective Rossby numbers lead to a transition from equatorial to polar cooling and then back to equatorial cooling, in qualitative agreement with the results of [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF]. However, the peak level of polar cooling and the effective Rossby numbers at which it occurs depend on 𝑃 𝑟 [START_REF] Kvorka | A numerical model of convective heat transfer in Titan's subsurface ocean[END_REF]. It is therefore plausible that our results with a heterogeneous inner boundary heat flux may depend on 𝑃 𝑟.

Polar vs. equatorial cooling is linked to the role of strong zonal jets in the dynamics. Coriolis force deflects westward flow towards the rotation axis. Therefore at the equator zonal jets quench low-latitude radial flows and sweep thermal structures there hence leading to polar cooling [START_REF] Aurnou | Convective heat transfer and the pattern of thermal emission on the gas giants[END_REF][START_REF] Yadav | Effect of shear and magnetic field on the heat-transfer efficiency of convection in rotating spherical shells[END_REF], in contrast to the resulting spiraling motion at polar regions [START_REF] Olson | A polar vortex in the Earth's core[END_REF][START_REF] Cao | Geomagnetic polar minima do not arise from steady meridional circulation[END_REF]. Furthermore, [START_REF] Yadav | Effect of shear and magnetic field on the heat-transfer efficiency of convection in rotating spherical shells[END_REF] found that the increase in heat transfer with increasing Rayleigh number is much stronger in the polar regions than in the equatorial region. In addition, the amplitude of the zonal jets depends on the Rossby number in a nonmonotonic way [START_REF] Yadav | Effect of shear and magnetic field on the heat-transfer efficiency of convection in rotating spherical shells[END_REF]. These issues may pose problems in isolating the role of inner boundary heterogeneity in the latitudinal distribution of outer boundary heat flux.

The heterogeneous inner boundary heat flux leads to deviations from axisymmetry and equatorial symmetry. Above regions of high imposed inner boundary heat flux, plumes develop and carry more heat efficiently to the outer boundary. Consequently, these patches of high inner boundary heat flux project high heat flux at the same regions on the outer boundary, especially in models with large amplitude of inner boundary heterogeneity 𝑞 * 𝑖 (Figs. 123). Such a correlation is particularly expected when the shell is thin, 𝑞 * 𝑖 is large and convective effects dominate over rotational. More dominant rotational effects reduce the boundary effect outside the TC.

All models show a slight tendency of positive outer boundary heat flux towards the low latitudes of the northern hemisphere (Fig. 6), as in the imposed inner boundary heat flux pattern (Figs. 1 and5). Far more importantly, preferred longitudes of positive outer boundary anomalies are evident in the maps (Figs. 23) and in particular in the longitude averages (Fig. 9) of all models. The outer boundary heat flux longitude dependence is highly correlated with the inner boundary pattern (Fig. 9 and Table 3). The correlations degrade when more rapid rotation and weaker boundary heterogeneity are considered (Fig. 9 and Table 3). This correlation reduction follows two distinctive branches, both incorporated in the horizontal effective Rossby number 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 (Fig. 13). In equatorial and polar cooling cases with relatively low Rossby numbers, both of which appear in the dynamical transitional regime [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF] but close to the weakly nonlinear regime [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF], according to the scaling laws in this regime, the convective length scale (i.e. the thickness of the axial columns) decreases with decreasing Ekman number 𝐸 [START_REF] Aubert | A systematic experimental study of rapidly rotating spherical convection in water and liquid gallium[END_REF][START_REF] King | Flow speeds and length scales in geodynamo models: the role of viscosity[END_REF][START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF]. These thinner axial convective columns lead to stronger pattern distortion by convective mixing and thus reduction in correlation between the inner and outer boundary heat flux spatial distributions. In contrast, in polar cooling cases with relatively large Rossby numbers which also appear in the dynamical transitional regime [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF] but close to the non-rotating regime [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF], obviously rotation does not play a major role; here weaker boundary heterogeneity 𝑞 * 𝑖 results in weaker boundary control and hence lower correlations. We find that the correlation reduction due to this latter effect is much weaker. In addition, improved cross-correlations are more prominent in the polar cooling cases.

An intriguing feature of our results which is visible in the longitudedependent profiles of the outer boundary heat flux of some models (Fig. 9) is the emergence of a large-scale order 2 pattern, most notably in the polar cooling cases. The signature of this order 2 part in the imposed boundary condition is weak, though it exists; note that the order 2 peaks on the outer boundary (Fig. 9) coincide with inner boundary peaks (see Fig. 8). This weak imposed large-scale heterogeneity is amplified by convection across the thin shell. Convection may mix and hence filter small scales, i.e. enhance (relatively) large scales. [START_REF] Aubert | Thermo-chemical wind flows couple Earth's inner core growth to mantle heterogeneity[END_REF] obtained an order 1 pattern of co-density flux on the inner boundary from an imposed tomographic (i.e. dominantly order 2) pattern on the outer boundary of their geodynamo models. This filtering of small scales by mixing appears when convection is sufficiently vigorous; when convection is closer to onset, the imposed boundary heterogeneity is maintained on the other boundary [START_REF] Gubbins | Melting of the Earth's inner core[END_REF].

The boundary control effect is further quantified by the part of the outer boundary heat flux that is non zonal symmetric. This part increases with increasing horizontal Rayleigh number 𝑅𝑎 𝐻 (Fig. 12), i.e. stronger mean convection and larger amplitude of imposed heat flux heterogeneity lead to more longitude dependence and possibly more deviations from equatorial symmetry in the resulting outer boundary heat flux.

Our models span different values of 𝐸, 𝑅𝑎 and 𝑞 * 𝑖 , but we use a fixed 𝑃 𝑟 = 1 value in all cases. In addition, we use no-slip mechanical boundary conditions in all simulations, while under free-slip conditions rapid rotation may yield strong zonal flows which in turn may sweep the equatorial heat flux and favor polar cooling [START_REF] Aurnou | Convective heat transfer and the pattern of thermal emission on the gas giants[END_REF][START_REF] Yadav | Effect of shear and magnetic field on the heat-transfer efficiency of convection in rotating spherical shells[END_REF]. We also neglect all convective effects related to salts (buoyancy due to their crystallization or melting). Caution should thus be taken in the interpretation of the results and the projections to the icy moons since any result may additionally depend on the 𝑃 𝑟 value, the type of mechanical boundary condition or possibly the impact of salts on the convection.

Our results may have implications to the relations between different layers at the interior of icy moons. We hypothesize two main largescale characteristics of the heat flux at the top of the subsurface oceans. First, polar cooling [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF]. Second, a dominant largescale (order 2 in the case of the particular imposed pattern used here) longitude dependence. Moreover, can the heat flux at the bottom of the subsurface ocean be inferred from surface observations? Our results show that the outer boundary heat flux pattern is related to the imposed inner boundary heat flux pattern but some distortion may prevail. In addition, the top icy layer may further modulate the signal from the top of the ocean to the surface.

Conclusions

There are two differences between our models and those of [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF], both concern the inner boundary thermal condition. First, [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] set a fixed temperature whereas we set a fixed heat flux. Second, their inner boundary temperature is homogeneous whereas our heat flux is laterally varying. Although the amplitude of imposed inner boundary heat flux heterogeneity in our models is rather moderate (𝑞 * 𝑖 < 3, see Table 1), several important results arise from this heterogeneity:

• Persistent preferred longitudes of intense outer boundary heat flux are highly correlated with longitudes of intense inner boundary heat flux. For the equatorial cooling cases, this correlation decreases strongly with increasing rotation rate. In contrast, for the polar cooling cases, the correlation decreases weakly with decreasing amplitude of inner boundary heterogeneity. • In addition, a non-trivial large-scale order 2 pattern is amplified by the convection in the thin shell, which may have a positive retroaction on the dynamics of the high-pressure ice mantle. • Deviations from axisymmetry and equatorial symmetry in the outer boundary heat flux increase with the horizontal Rayleigh number which contains the amplitude of the inner boundary heterogeneity.

In contrast, polar vs. equatorial cooling is mostly controlled by inertial effects, as was found by [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] for homogeneous boundary conditions, with the latitudinally equilibrated inner boundary heterogeneity (Fig. 5 and Table 3) acting to reduce the amplitude of this effect (Fig. 11). [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF] applied their results to Titan's ocean by considering the parametric diagram of [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF]. They found that their models are in the transitional regime between rapidly rotating and non-rotating convection, while Titan's ocean (considering uncertainties in estimated parameters) is either in the same transitional regime or in the non-rotating regime. Because Titan's sub-surface ocean is expected to have stronger inertial effects than our models (see Fig. 10 of [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF], our results seem to suggest that the impact of the heat flux heterogeneity at the bottom of Titan's ocean on the latitudinal cooling pattern at the top of Titan's ocean is secondary. The thinning of Titan's ice shell and corresponding heat flux enhancement in polar regions [START_REF] Kvorka | Does Titan's long-wavelength topography contain information about subsurface ocean dynamics?[END_REF] can be explained by the polar cooling regime [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF], and does not require additional heat flux anomalies from the interior.

While the main degree-2 pattern in ice shell thickness and associated heat flux interpreted by [START_REF] Kvorka | Does Titan's long-wavelength topography contain information about subsurface ocean dynamics?[END_REF] can be explained by the oceanic dynamical regime [START_REF] Amit | Cooling patterns in rotating thin spherical shells -Application to titan's subsurface ocean[END_REF], the longitudinal variations observed at highest degree (see Fig. 11 of [START_REF] Kvorka | Does Titan's long-wavelength topography contain information about subsurface ocean dynamics?[END_REF] may reflect bottom driven heat flux anomalies. Even though the longitudinal variations are attenuated by ocean dynamics, especially for cases with high Rossby numbers, analysis of these longitudinal fluctuations may provide key constraints on the heat flux coming out of the high-pressure ice mantle, thus constraining its convective dynamics and the efficiency of chemical transfer. However, in order to consistently link the heat flux anomalies to both ocean and highpressure ice mantle dynamics, a more elaborated description of the ice/ocean interface taking into account melting/freezing processes (see [START_REF] Kvorka | Does Titan's long-wavelength topography contain information about subsurface ocean dynamics?[END_REF], for the upper ocean interface) would ultimately be required. As mentioned above, retroaction between the ocean and ice mantle may impact the heat flux pattern and possibly amplify some specific degrees. Future modeling efforts are needed to take into account a more appropriate boundary condition and its feedback on ocean/ice dynamics. Despite these limitations, the present work provides the first proof of concept on how ocean dynamics may be influenced by bottom heat flux anomalies.

  Fig. 1. Imposed pattern of inner boundary heat flux taken from the output of a high-pressure ice convection simulation applied to Titan's mantle (Choblet et al., 2017b). The heat flux is expanded up to spherical harmonic degree and order 30. (a) and (b) show different projections of the imposed heat flux pattern.

Fig. 2 .

 2 Fig. 2. Outer boundary heat flux anomalies in snapshots (left) and time-averages (right) for cases 1 (a,b), 12 (c,d) and 14 (e,f) which are characterized by equatorial cooling. Dashed latitude lines denote the tangent cylinder. Note the different color scales.

Fig. 3 .

 3 Fig. 3. As in Fig. 2 for cases 3 (a,b), 6 (c,d) and 8 (e,f) which are characterized by polar cooling.

Fig. 4 .

 4 Fig. 4. Equatorial (left) and North polar (right) views of the time-averaged outer boundary heat flux anomalies for cases 12 (a,b) characterized by equatorial cooling and 7 (c,d) characterized by polar cooling. Dashed latitude lines denote the tangent cylinder.

Fig. 5 .

 5 Fig. 5. Imposed inner boundary heat flux anomaly averaged along latitude lines (blue) and its equatorially symmetric part (red) vs. co-latitude. Vertical lines denote the tangent cylinder.

  Next we investigate the relation between the patterns of inner and outer boundary heat flux anomalies. Comparison of Fig. 1 with Figs. 2 and 3 shows traces of local boundary control in the form of regions of strongest heat flux on the outer boundary corresponding to most intense imposed heat flux on the inner boundary, see e.g. the features 60

Fig. 6 .

 6 Fig. 6. Outer boundary heat flux anomaly averaged along latitude lines vs. co-latitude. Vertical lines denote the tangent cylinder. Equatorial and polar cooling cases are shown in (a) and (b) respectively.

Fig. 7 .

 7 Fig. 7. Non-dimensional radial velocity in snapshots of cases 2 (left; equatorial cooling) and 7 (right; polar cooling) just below the Ekman boundary layer (top) and at zonal averages (bottom).

Fig. 8 .

 8 Fig. 8. Inner boundary heat flux averaged along longitude lines vs. longitude.

Fig. 9 .

 9 Fig. 9. Outer boundary heat flux anomaly averaged along longitude lines vs. longitude. Thick gray lines represent inner boundary heat flux.

Fig. 11 .

 11 Fig. 11. 𝑞 ℎ∕𝑙 𝑜 vs. 𝑅𝑎 𝑇 ∕𝑅𝑎 𝑤 𝑇 (a) and 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 (b). Diamonds/Circles represents equatorial/polar cooling patterns, respectively. Red symbols in (a) are the results from Amit et al. (2020). Note that in (a) the temperature Rayleigh number 𝑅𝑎 𝑇 is used for comparison with the models of Amit et al. (2020), whereas in (b) the flux Rayleigh number 𝑅𝑎 is used.

Fig. 12 .

 12 Fig. 12. 𝑅𝑀𝑆(𝛿𝑞 𝑛𝑧𝑠 𝑜 ) vs. 𝑅𝑎 𝐻 . Diamonds/Circles represent equatorial/polar cooling patterns, respectively.

Fig. 13 .

 13 Fig. 13. Correlations between outer and imposed inner boundary heat flux patterns vs. 𝑅𝑎 𝐻 ∕𝑅𝑎 𝑤 𝑇 (15). In (b) possible azimuthal shifts are accounted for. Diamonds/Circles represents equatorial/polar cooling patterns, respectively. Note slightly different y-axis scales in (a) vs. (b).

Table 1

 1 Models inputs and global outputs.

	Model	𝐸	𝑅𝑎	𝑅𝑎 𝑇 𝑅𝑎 𝑐	𝑞 * 𝑖	𝓁 𝑚𝑎𝑥	𝑛 𝑟	𝑅𝑒	𝑅𝑜	𝑅𝑜 𝑐	𝑅𝑎 𝑇 𝑅𝑎 𝑤	
	1	1e-03	2e+05	5.05	0.57	96		49	20.63	0.021	0.45	0.63	9693
	2	1e-03	2e+05	4.98	1.14	96		49	20.20	0.020	0.45	0.63	9290
	3	1e-03	1e+06	11.80	0.57	96		49	53.64	0.054	1.00	3.16	5099
	4	1e-03	1e+06	11.79	1.14	96		49	53.26	0.053	1.00	3.16	4823
	5	3e-04	1e+07	17.70	2.00	128	73	125.92	0.038	0.95	5.20	1850
	6	3e-04	5e+07	49.59	1.04	128	65	275.09	0.083	2.12	25.98	2132
	7	3e-04	5e+07	47.56	2.01	128	73	275.26	0.083	2.12	25.98	1479
	8	3e-04	5e+07	44.89	2.68	128	73	274.78	0.082	2.12	25.98	1313
	9	3e-04	8e+07	65.33	2.00	128	73	342.59	0.103	2.68	41.57	1537
	10	1e-04	1e+06	2.41	0.82	170	81	12.85	0.001	0.10	0.10	77
	11	1e-04	3e+06	5.83	0.82	170	81	36.04	0.004	0.17	0.30	188
	12	5e-05	1e+07	7.50	0.60	192	97	124.65	0.006	0.16	0.35	138
	13	5e-05	1e+07	7.48	1.20	192	97	120.83	0.006	0.16	0.35	134
	14	5e-05	1e+07	7.42	1.81	192	97	133.50	0.007	0.16	0.35	152
	𝐸 is the Ekman number (4), 𝑅𝑎 is the Rayleigh number (6),	𝑅𝑎 𝑇 𝑅𝑎 𝑐

Table 2

 2 Plausible ranges for the dimensionless numbers for Titan and Ganymede.

		𝑃 𝑟	𝐸	𝑅𝑎	𝑅𝑎 𝑇	𝑞 * 𝑖
	Titan	10	2e-12 -5e-11	1e+25 -1e+31	3e+20 -5e+22	0.2 -2
	Ganymede	10	7e-13 -3e-10	3e+27 -4e+30	5e+18 -1e+23	0.2 -2

Values showed for 𝑃 𝑟, 𝐸, 𝑅𝑎 and 𝑅𝑎 𝑇 are similar to those in Soderlund (

2019

) (e.g. interior models values from

Vance et al. (

Table 3

 3 Results for the outer boundary heat flux.

		Model	𝑞 𝑖 ℎ∕𝑙	𝑞 𝑜 ℎ∕𝑙	Cooling	RMS (𝛿𝑞 𝑛𝑧𝑠 𝑜 )	 𝜙	𝑚𝑎𝑥 ( 𝑐𝑟 𝜙 )	 2 ∕ 2	 10 ∕ 10
		1	0.004	-0.055	Equatorial	0.059	0.499	0.499	2.023	0.511
		2	0.008	-0.051	Equatorial	0.092	0.709	0.709	2.606	0.793
		3	0.004	0.011	Polar	0.094	0.571	0.585	1.545	0.520
		4	0.008	0.016	Polar	0.126	0.663	0.663	1.844	0.563
		5	0.016	0.024	Polar	0.210	0.717	0.734	1.907	0.640
		6	0.008	0.067	Polar	0.163	0.466	0.620	1.572	0.367
		7	0.016	0.072	Polar	0.184	0.520	0.675	3.964	0.941
		8	0.021	0.067	Polar	0.195	0.487	0.599	4.384	0.627
		9	0.015	0.078	Polar	0.174	0.480	0.630	3.964	0.644
		10	0.006	-0.009	Equatorial	0.042	0.419	0.520	1.167	0.253
		11	0.006	-0.035	Equatorial	0.098	0.399	0.416	2.446	0.338
		12	0.010	-0.033	Equatorial	0.155	0.164	0.203	0.398	0.230
		13	0.015	-0.030	Equatorial	0.152	0.417	0.421	1.129	0.325
		14	0.020	-0.039	Equatorial	0.159	0.413	0.432	1.126	0.178
	𝑞	ℎ∕𝑙 𝑜

represents equatorial vs. polar cooling (18). 𝑞 ℎ∕𝑙 𝑖

Acknowledgments

We thank Guy Moebs for the valuable technical support with the numerical dynamo simulations. Computations were performed on the CCIPL (Le Centre de Calcul Intensif des Pays de la Loire). The source of MagIC (numerical dynamo code) is available at https://magicsph.github.io/. F. T-N. acknowledges São Paulo Research Foundation (FAPESP), Brazil for grant 2018/07410-3. This study acknowledges the financial support from the French Agence Nationale de Recherche, project ANR COLOSSE (grant N • ANR-2020-CE49-0010).

Data availability

Numerical dynamo simulations code: https://magic-sph.github.io/.