Europa's ocean translates interior tidal heating patterns to the ice-ocean boundary

The circulation in Europa's ocean determines the degree of thermal, mechanical and chemical coupling between the ice shell and the silicate mantle. Using global direct numerical simulations, we investigate the effect of heterogeneous tidal heating in the silicate mantle on rotating thermal convection in the ocean and its consequences on ice shell thickness. Under the assumption of no salinity or ocean-ice shell feedbacks, we show that convection largely transposes the latitudinal variations of tidal heating from the seafloor to the ice, leading to a higher oceanic heat flux in polar regions. Longitudinal variations are efficiently transferred when boundary-driven thermal winds develop, but are reduced in the presence of strong zonal flows and may vanish in planetary regimes. If spatially homogeneous radiogenic heating is dominant in the silicate mantle, the ocean's contribution to ice shell thickness variations is negligible compared to tidal heating within the ice. If tidal heating is instead dominant in the mantle, the situation is reversed and the ocean controls the pole-to-equator thickness contrast, as well as possible longitudinal variations.

thermally driven flows in Europa's ocean, neglecting salinity and feedback effects of the ice • Heterogeneous tidal heating in the mantle modifies the mean circulation in Europa's ocean and could drive large-scale thermal winds • The tidal heating anomaly in latitude is efficiently translated upwards, leading to a higher heat flux into the ice shell at the poles

Supporting Information:
Supporting Information may be found in the online version of this article.
our study neglect salinity effects as well as those due to phase changes at the ice-ocean boundary.Soderlund et al. (2014) suggested that the heat flux from Europa's ocean is higher at the equator, thereby melting and thinning the ice at low latitudes, a configuration called equatorial cooling.However, this work and following studies (Amit et al., 2020;Kvorka & Čadek, 2022;Soderlund, 2019) considered Europa's ocean to be uniformly heated by the silicate mantle, which is possibly unrealistic if tidal heating exceeds radiogenic heating in the mantle.
Unlike radiogenic heating, tidal heating is spatially heterogeneous (e.g., Beuthe, 2013;Sotin et al., 2009;Tobie et al., 2005).For Europa's silicate mantle, it is predicted to be significantly greater at the poles than at the equator, with longitudinal variations of order 2 at low latitudes as well (Figure 1 and Figure S2 in Supporting Information S1).To measure how strongly heterogeneous this pattern is, we introduce the parameter q* = Δq/q 0 , the ratio between the maximum heat flux difference (between the poles and the sub or anti-Jovian points) and the laterally averaged heat flux, q 0 .For pure tidal heating, q* ∼ 0.91, that is, heat flux variations are of the same order as the mean heat flux (Beuthe, 2013;Tobie et al., 2005).However, tidal heating is superposed on radiogenic heating, which is nominally homogeneous (q* = 0).The relative amplitude of tidal heating compared to radiogenic heating in Europa's silicate mantle is poorly constrained (Hussmann et al., 2016) and has likely varied during Europa's evolution (Běhounková et al., 2021).Therefore, q* can a priori lie anywhere between 0 and 1 (Figure 1) and we consider a sweep in q* to evaluate the ocean's response in various scenarios.
How rotating thermal convection responds to a heterogeneous heating has been investigated for Earth's outer core (e.g., Davies & Mound, 2019;Dietrich et al., 2016;Gubbins et al., 2011;Mound et al., 2019;Olson et al., 2015;Sahoo & Sreenivasan, 2020).It was found that the heat flux pattern at the core-mantle boundary could not reach the inner core (Davies & Mound, 2019).However, Earth's outer core is a much thicker shell compared to Europa's ocean, and we explore less rotationally controlled regimes compared to the geostrophic regime investigated for the Earth (Davies & Mound, 2019;Mound & Davies, 2017).In the context of icy moon oceans, motivated by convection patterns in high-pressure ices of Titan or Ganymede (Choblet et al., 2017), it was recently showed that small scale heating patterns do not influence convection at a global scale, and that the ocean cooling configuration (polar or equatorial) is unaffected (Terra-Nova et al., 2022).Narrow high-heating bands at the seafloor of Enceladus would also be mixed by the ocean and only have a diffuse imprint beneath the ice (Kang, Marshall, et al., 2022).Large-scale heating patterns could lead to different conclusions, because large-scale temperature gradients can drive a mean thermal-wind circulation (Dietrich et al., 2016).In the present work, we investigate (a) how the large-scale tidal heating of Europa's mantle impacts its ocean circulation, (b) whether large-scale thermal anomalies from the seafloor can be transposed up to the ice-ocean boundary, and (c) if resulting oceanic heat flux variations could affect the ice crust equilibrium.

Materials and Methods
To address these questions, we model the ocean by a fluid contained in a rotating spherical shell and driven by (a) vertical convection due to the heating of the ocean from below, and (b) horizontal convection (Gayen & Griffiths, 2022;Hughes & Griffiths, 2008) due to lateral variations of the basal heat flux.A competition between the vertically forced and horizontally forced convection will take place (Couston et al., 2022) and determine the properties of the large-scale circulation.In the presence of rotation, thermal winds develop due to horizontal temperature gradients (Dietrich et al., 2016, and references therein).To investigate this, we perform direct numerical simulations (DNS) of turbulent thermal convection with an imposed heat flux at the bottom boundary.The basal heat flux is either homogeneous (representing radiogenic heating) or increasingly heterogeneous following the tidal heating pattern in the mantle.These methods are described in more detail below.
From the ocean model, we extract relative heat flux variations obtained at the ice-ocean boundary.In Section 3.5, we use them as inputs in an ice thickness model (see Appendix B and Text S2 in Supporting Information S1) to evaluate the possible impact of a heterogeneous oceanic heat flux on the ice shell thermal equilibrium.Note that this is not a coupled model, that is, the ocean dynamics are solved separately from the ice shell equilibrium.Therefore, important feedback processes are missing (e.g., topography of the ice-ocean boundary, freezing point depression with pressure).The ice thickness maps obtained here are not meant to be definitive quantitative predictions for Europa, but are rather used as a tool to compare the relative effects of oceanic heat flux variations versus tidal heating variations within the ice. 10.1029/2023AV000994 3 of 21

Set of Equations
We numerically model the ocean circulation by solving for turbulent thermal convection in rotating spherical shells using the open-source code MagIC (Christensen et al., 2001;Schaeffer, 2013;Wicht, 2002) (URL: https:// magic-sph.github.io/).The ocean is modeled as a fluid of kinematic viscosity ν, density ρ, thermal diffusivity κ, heat capacity c p , and thermal expansivity α, contained in a spherical shell of outer and inner radii r o and r i , thickness D = r o − r i , rotating at a uniform rate Ω about the vertical axis z.The aspect ratio of the shell, defined as η = 1 − D/r o is fixed at 0.9, close to the aspect ratio of Europa's ocean (Table A1).In the following, all of the variables are made dimensionless using the shell thickness D for the length scale and the viscous time D 2 /ν for the timescale.The temperature scale is Dq 0 /k, where  0 = ⟨bot⟩ℎ = ⟨ |=  ⟩ℎ is the horizontal average of the imposed basal heat flux and k = ρc p κ is the thermal conductivity of the fluid.In the Boussinesq approximation, and in the rotating frame of reference, the set of dimensionless equations governing the system are: where Equation 1 accounts for conservation of momentum (Navier-Stokes equation), Equation 2 accounts for energy conservation and Equation 3 accounts for mass conservation.Here, e r and e z are the radial and vertical unit vectors, respectively, where the vertical is aligned with the rotation axis.T, p, and u are the dimensionless fluid temperature, pressure and velocity, and the gravity field is taken to be linear g = g o r/r o .At a fixed aspect ratio, the solution of the system of Equations 1-3 depends on the relative value of three dimensionless parameters: the flux Rayleigh number Ra F , the Ekman number E, and the Prandtl number Pr: (4) Since we considered a heterogeneous basal heat flux, a fourth dimensionless parameter is introduced to represent the relative basal heat flux anomaly,  (5) Here, Δq bot is the heat flux difference between the maximum and the minimum basal heat flux.Given the tidal pattern (Figure 1), the maximum heat flux point is always at the pole, and the minimum heat flux is located at the sub-Jovian and anti-Jovian points.

Boundary Conditions
The top and bottom boundaries are impenetrable and either both no-slip or both stress-free.We enforce angular momentum conservation in the fluid bulk.No-slip mechanical boundary conditions (BC) could seem the most natural ones to employ, but given the artificially high viscous forces in our simulations (and in all DNS of convection in a 3D shell), necessary to make the problem computationally tractable, no-slip conditions lead to the development of thick Ekman boundary layers which are not planetary relevant (e.g., Glatzmaier, 2002;Kuang & Bloxham, 1997).We therefore also consider stress-free BC.
We use a Neumann boundary condition at the bottom of the ocean-fixed flux q bot (θ, ϕ)-to represent the heat flux from the silicate mantle.We use a Dirichlet boundary condition at the top of the ocean-fixed temperature T top -to represent the fixed melting temperature.Note that we also investigated cases where the temperature is fixed at both boundaries in order to compare with and extend previous studies mostly done with Dirichlet BC (Amit et al., 2020;Kvorka & Čadek, 2022;Soderlund, 2019;Soderlund et al., 2014).These Dirichlet simulations are described in the Supplementary Information.
To impose a varying heat flux pattern on the bottom boundary, the map represented in Figure 1c is decomposed into spherical harmonics (Figure S2 in Supporting Information S1).The bottom heat flux can thus be expressed as where θ is the colatitude and ϕ the longitude.Here ℓ and m denote the spherical harmonic degree and order respectively, and the spherical harmonics     are completely normalized: with     the associated Legendre functions.We keep only the most significant modes, and therefore truncate the decomposition at degree ℓ = 4 (Figure S2 in Supporting Information S1).The corresponding spherical harmonic coefficients,    bot are provided in Table 1.For the ocean numerical simulations, a sweep in q* is performed by keeping the average heat flux fixed and progressively increasing the amplitude of all the other coefficients (l, m) ≠ (0, 0).The relative amplitude of the modes (l, m) ≠ (0, 0) between each other is always the same, so that the pattern is geometrically fixed.Coefficients for Dirichlet simulations are also provided in Table S2 in Supporting Information S1.

Dimensionless Parameters
In this work, the goal is to survey how different possible contributions of tidal heating in the mantle would impact ocean circulation and heat transport, not to quantify tidal heating in Europa's mantle.Instead, we consider a wide range of possible heat flux amplitudes at the seafloor, going from 6 to 46 mW m −2 .The lower bound corresponds to the case of pure radiogenic heating (Tobie et al., 2003) and negligible tidal heating, whereas the upper bound accounts for the possibility of a dominant, Io-like tidal dissipation (Howell, 2021) (see also Text S1 in Supporting Information S1).With the other physical parameters listed in Table A1, we obtain for Europa's ocean Pr ≈ 11, E = 7 × 10 −12 (∈ [5.0 × 10 −12 , 9.1 × 10 −12 ]) and Ra F = 5.4 × 10 27 (∈ [6.3 × 10 26 , 2.8 × 10 28 ]).
Using scaling laws of convection (Appendix A), a temperature-based Rayleigh (Ra T , Equation A1) can be estimated from the flux Rayleigh, allowing us to locate Europa on the regime diagram of Figure 2. Europa's ocean falls into the transitional regime of Gastine et al. (2016), and is less rotationally constrained than moons like Enceladus (Soderlund, 2019).Adapting the flux-based regime diagram of Long et al. (2020) leads to the same conclusion (see Figure S1 in Supporting Information S1), but their systematic study does not allow to delineate the boundary between transitional and non-rotating (NR) regimes.Therefore, we work with the regime diagram of Gastine et al. (2016).Because planetary parameters are out-of-reach with current computational capabilities, we work at more moderate parameters, chosen such that the system is in the same convective regime as Europa.Hence, we emphasize here a weakly rotating convective regime (Ra F = 1.26 × 10 9 , E = 3 × 10 −4 , Pr = 1, yellow star in Figure 2).Because the degree of rotational influence is debated (Ashkenazy & Tziperman, 2021;Bire et al., 2022), we provide results for a more rotationally constrained regime (blue stars in Figure 2) in Section 3.3 and in Text S5 in Supporting Information S1.
Note that DNS parameters are orders of magnitude wrong compared to realistic values (see also Table S4 in Supporting Information S1).One way to interpret the simulations' Rayleigh, Ekman and Prandtl numbers is to convert them into dimensional parameters.One possibility is to make Ω, D, g, α, ρ and c p equal to the planetary values (Table A1), leading to ν DNS = κ DNS = 77.6 m 2 s −1 and q DNS = 6.6 × 10 4 W m −2 .These values are unphysical, meaning that our simulations are both too strongly diffusive (viscous and thermal diffusion) and too strongly forced to overcome this extra-diffusion, similarly to all other numerical studies of icy moon ocean circulation using DNS (e.g., Amit et al., 2020;Kvorka & Čadek, 2022;Soderlund, 2019;Soderlund et al., 2014;Terra-Nova et al., 2022).Studying dependences with input parameters is hence required for any extrapolation to the real planetary system (see Section 3.4, where lower Ekman numbers are considered).
From negligible to dominant tidal heating in the silicate mantle, both the average heat flux, 〈q〉 h , and the relative heat flux variations, q*, should increase (Figure 1).The range of average heat fluxes implied is taken into account in Europa's vertical error bars in Figure 2, and it does not change the convective regime into which Europa falls.We hence neglect the change in the average heat flux and work at a fixed Rayleigh number for a given sweep in q*.However the amount of tidal heating relative to radiogenic heating will change the extent to which the heating of the ocean is heterogeneous.To illustrate this, we perform a sweep in q*, that is, gradually increase the bottom heat flux heterogeneity while keeping the pattern the same (top row in Figure 3).Hence, the imposed heat flux at the bottom boundary is either homogeneous (representing radiogenic heating, q* = 0) or increasingly heterogeneous following the tidal heating pattern (q* = [0.46, 0.87, 1.30, 1.72, 2.16]).Note that we chose to go beyond q* = 1 to have a better physical insight on the behavior of the system, and because an interesting transition was observed around q* ∼ 1 for Dirichlet simulations (see Figure S5 in Supporting Information S1).
Note.The coefficients for the tidal model correspond to the map of Figure 1c, with the coefficients normalized by the (ℓ = 0, m = 0) coefficient.For the ocean numerical simulations, the average heat flux is kept fixed (   =0,=0  is constant), but the amplitude of all the other coefficients (l, m) ≠ (0, 0) is progressively increased to increase q*.The relative amplitude of the modes (l, m) ≠ (0, 0) between each other is always the same.Note that the relative amplitudes between the coefficients (4,0), (4,2), and (4,4) compared to the mode (2,2) very slightly differ between the simulations and the tidal model because of an initial error in our tidal dissipation code.As these modes have the smallest amplitude, this does not significantly affect the tidal pattern nor the fluid's response.

Numerical Methods
MagIC solves for the system of Equations 1-3 in spherical coordinates, and employs a pseudo-spectral method, using Chebyshev polynomials in the radial direction and spherical harmonics in the longitudinal and latitudinal directions.The fast spherical harmonic transform library SHTns have been employed (Schaeffer, 2013) (URL: https://bitbucket.org/nschaeff/shtns/).The equations are integrated in time with a mixed implicit/explicit time stepping scheme.Here, we employ the commonly used CNAB2 second-order scheme which combines a Crank-Nicolson scheme for the implicit terms and a second-order Adams-Bashforth scheme for the explicit terms  2016) representing the supercriticality, Ra T E 3/4 , as a function of the Ekman number, for a fixed Prandtl number Pr = 1 (Direct Numerical Simulations) and Pr = 10 (icy moons).The gray shaded area represents a region where the flow is stable, Ra T < Ra c (stability threshold determined by Gastine et al. (2016) in the absence of horizontal forcing variations and for a shell aspect ratio of 0.6).Above that threshold, the flow is convectively unstable, and becomes fully non-linear once Ra T ≳ 6Ra c .Different regimes of convection are observed depending on the strength of the forcing compared to rotational effects.Each regime threshold is indicated by a black continuous line, and is taken from references GA16 (Gastine et al., 2016), KC22 (Kvorka & Čadek, 2022), EN14 (Ecke & Niemela, 2014) and JU12 (Julien et al., 2012).At a fixed E, as Ra T is increased, the convection is more vigorous, and the influence of rotation decreases, as indicated by the Geostrophic, Transitional, and Non-rotating regimes.The shaded blue region represents the region where Kvorka and Čadek (2022) identified a polar cooling behavior, that is, a higher heat flux at the pole.The dashed gray lines are lines of constant convective Rossby number   =  1∕2  Pr −1∕2 .The stars represent the simulations described in the present study (see also Table S5 in Supporting Information S1).(non-linear terms and Coriolis force).For full details on the numerical methods, we refer the reader to MagIC documentation available online (URL: https://magic-sph.github.io/).MagIC is parallelized using both OpenMP (URL: http://openmp.org/wp/) and MPI (URL: http://www.open-mpi.org/) and is designed to be used on supercomputers.The simulations were ran on the Lonestar6 system at the Texas Advanced Computing Center (URL: http://www.tacc.utexas.edu).To address differences between thermal and mechanical BC, rotation regimes, and value of q*, 79 numerical runs were performed in total, as listed in Table S5 in Supporting Information S1.The high Ekman number simulations are typically ran during 6 days (144 hr) on 128 CPUs.The lowest Ekman cases were typically ran during 14 days (336 hr) on 200 CPUs.

Mean Circulation
Figure 3 shows the mean oceanic circulation obtained for no-slip and stress-free BC for an increasing heterogeneity q*.In no-slip simulations with homogeneous heating (q* = 0), moderate zonal flows are present.By "zonal flows," we refer to the azimuthally averaged azimuthal velocity,  ⟨⟩ .The flow is prograde (eastward) at the equator, retrograde (westward) at intermediate latitudes, and then prograde again in the polar regions.Axial convection with columns aligned along the rotation axis develop in the equatorial region.
As q* is increased, the zonal flow amplitude decreases and the circulation becomes driven by the lateral variations of heating at the seafloor, a circulation called "thermal winds solution" hereafter, by analogy with thermal winds in atmospheres (Vallis, 2017).Thermal winds arise from the balance between buoyancy and Coriolis forces, and are described mechanistically in detail in Dietrich et al. (2016) in a configuration close to the present where the height of the spherical shell outer boundary measured from the equatorial plane outside the tangent cylinder (the vertical cylinder tangent to the seafloor), and s the cylindrical radius.The vertical velocity at z = ±H is found using non-penetration of the sloping boundaries,  (+) = −(−) = −⟨⟩∕ .Evaluating the previous equation on the equatorial plane (where s = r) leads to (Dietrich et al., 2016).This equation shows that, at the equator, in thermal wind balance, temperature gradients in longitude result in radial flows.At the equator, the tidal heating pattern is dominated by a m = 2 mode, with two zones of positive azimuthal heat flux gradient, and two zones of negative azimuthal gradient.The thermal wind balance hence predicts that two zones of downwelling and two regions of upwelling will develop, a feature which is indeed observed in our simulations (Figure 3a and Figure S4 in Supporting Information S1).
Nonlinearities due to temperature advection or inertia complicate this simple picture, and are responsible for the asymmetry between the narrow downwelling regions and the wide upwelling (see also Text S3 in Supporting Information S1).
With stress-free BC, the absence of viscous dissipation at the boundaries allows for the development of very strong zonal flows.For homogeneous heating (q* = 0), we retrieve the results of Soderlund et al. (2014) despite their use of a different thermal boundary condition: the zonal flows are retrograde at the equator, prograde at high latitudes, and there is a mean upwelling at the equator.When q* is increased, the mean circulation is modulated by the inhomogeneous heating, but its principal features remain unchanged.We note that when we subtract the zonal flow, the mean azimuthal circulation exhibits a pattern reminiscent of the thermal winds, but their amplitude is an order of magnitude smaller than the zonal flows amplitude.We call this solution the "zonal winds solution" hereafter.

Vertical Transfer of the Bottom Anomaly
We now investigate whether the ocean effectively transfers the tidal heating pattern from the seafloor to the ice-ocean boundary.Figures 4a and 4b shows time and zonally averaged heat flux profiles at the bottom and top boundaries.In the zonal winds solution, with homogeneous basal heating, the heat flux peaks at the equator due to the equatorial upwelling, consistently with Soderlund et al. (2014).As q* is increased, even if an equatorial peak persists, the top heat flux follows the basal heat flux and becomes greater at the pole even for the smallest heterogeneity considered (q* = 0.43).The thermal winds solution exhibits a very similar behavior.
To quantify the latitudinal cooling configuration, we define the dimensionless parameter where   is the time-averaged heat flux at the top or at the bottom of the ocean, 〈•〉 low is the average over latitudes lower than 10°, and 〈•〉 high over latitudes greater than 80°.When q h/l < 0 the heat flux is larger in the equatorial region, and the ocean is said to be in an equatorial cooling configuration.Conversely, the ocean is in polar cooling if q h/l > 0. Figure 4c confirms that the ocean is in an equatorial cooling configuration only for the homogeneous heating case.Furthermore, the ocean efficiently transfers the latitudinal variations of basal heating up to the ice-ocean boundary, with   ℎ∕ top ≈ 0.75  ℎ∕ bot .For dominant tidal heating, the heat flux at the poles would be systematically higher than at the equator, that is, the ocean would be in a polar cooling configuration.
We next measure the relative amplitude of longitudinal heat flux variations at the bottom and at the top of the ocean: Here, 〈•〉 lat denotes average in latitude and 〈•〉 h denotes horizontal average on the sphere.Figures 4d-4f show that in the thermal winds solution, longitudinal variations are almost entirely transferred upward thanks to the mean equatorial loops (Text S3 in Supporting Information S1).In the zonal winds solution, longitudinal variations are smeared, and only about 20% of the original anomaly persists at the ice-ocean boundary.

Sensitivity to Boundary Conditions and Rotation Regime
In each convective regime, we investigated the sensitivity to the thermal BC by imposing the temperature (Dirichlet BC) at the seafloor instead of the heat flux.We show in the Supplementary Information (Text S4 in Supporting Information S1) that the results are very similar.The only notable exception is that in the weakly rotating regime, with stress-free, temperature-imposed BC, an abrupt transition occurs between the zonal and thermal winds solutions (Figure 5c and Figure S5 in Supporting Information S1).The origin of this transition will be investigated further, but is beyond the scope of the present study as flux-imposed BC are the most planetary relevant.9) at the top versus at the bottom of the ocean.When q h/l < 0, the ocean is in an equatorial cooling configuration (higher heat flux at the equator).
If the points lie around the 1:1 line, it means that all of the basal heat flux anomaly in latitude is retrieved at the top, and that the ocean efficiently transfers heat heterogeneities upwards.The colors have the same meaning as in other panels.(f) Longitudinal heat flux anomalies (Equation 10) at the top versus at the bottom of the ocean.
Given the uncertainty in Europa's ocean convective regime (Appendix A), there has been a debate on its cooling configuration (polar vs. equatorial).We focused here on a single, weakly rotating regime, but from Figure 2, Europa falls within the polar cooling regime identified by Kvorka and Čadek (2022).To sample this polar cooling regime, we performed more rotationally constrained simulations, while remaining in the transitional regime of Gastine et al. (2016).The regime of these simulations is called "moderately rotating" hereafter and correspond to blue stars on Figure 2. We show in Supporting Information S1 (Text S5) that the results are qualitatively unaffected by the change in the convective regime: thermal winds develop with no-slip BC, whereas zonal flows are dominant with stress-free BC.Latitudinal heat flux variations from the seafloor are mirrored at the ice-ocean boundary, and longitudinal variations are reduced if zonal winds develop.

Heat Anomaly
To explore the effect of going toward more planetary relevant parameters (smaller Ekman and larger Rayleigh), we reduce E and increase Ra T while remaining in a given convective regime, that is, keeping the same degree of rotational influence (Figure 2).For both the weakly rotating (yellow-red stars) and moderately rotating (blue stars) regimes, we explored four different Ekman numbers E ∈ [3 × 10 −4 , 1.5 × 10 −4 , 6 × 10 −5 , 3 × 10 −5 ] computational costs reasonable, this systematic study was done for Dirichlet, stress-free BC only.We chose Dirichlet thermal BCs because the choice of Ra T to be in the correct convective regime can be directly picked from the regime diagram of GA16 (Figure 2), and stress-free BCs because they allow us to sample both the zonal winds and thermal winds solutions at the same time given the existence of the transition.The corresponding set of simulations is listed in Table S5 in Supporting Information S1 and represented in Figure 2. As an equivalent to q* for Neumann simulations, we define the bottom temperature anomaly ΔT* = ΔT h /ΔT v where ΔT h is the peak-to-peak temperature difference at the seafloor.For each pair (Ra T , E), we perform an entire sweep in ΔT* from 0 to 1.63 (see Text S4.1 in Supporting Information S1), except for the smallest Ekman number for which we select a few ΔT*.
Figures 5a-5f represents the relative heat flux anomalies at the top versus at the bottom of the ocean for all the simulations of the systematic study (circles), plus the simulations with Neumann BCs described previously (triangles).Figures 5g-5j shows the fraction of the anomaly retrieved at the top of the ocean, as a function of the Ekman number.Figures 5a, 5b, 5d, 5e, 5g, and 5i shows that in both the weakly rotating and moderately rotating regimes, the latitudinal anomaly is efficiently transferred across all Ekman numbers investigated.It is slightly reduced for the weakly rotating cases because of the equatorial cooling configuration (higher heat flux at the equator) of this regime when the bottom heating is homogeneous (Kvorka & Čadek, 2022).The heat flux latitudinal anomaly is on the contrary slightly amplified at the top of ocean for moderately rotating cases due to the polar cooling configuration of this regime when the bottom heating is homogeneous (Kvorka & Čadek, 2022).
The preservation of latitudinal heating variations up to the ice-ocean interface is qualitatively in line with Gastine and Aurnou (2023) who underline that different convective dynamics occur in the equatorial region compared to the polar regions, even with a homogeneous bottom forcing.This regionalization of convection is proposed to be due to the increased misalignment between the rotation vector and gravity from the pole to the equator, as well as the extent of the tangent cylinder, fixed by the shell aspect ratio.The more decoupled the polar and equatorial regions are, the better they should reflect their respective (local) bottom heating, and the better large-scale latitudinal variations should be preserved.To verify this hypothesis, our study calls for more detailed analysis of the heat transport (de)coupling between polar and equatorial regions in weakly rotating, moderately rotating and rapidly rotating (RR) regimes.Besides, to ascertain extrapolation of this behavior to real Europa parameters (E ∼ 10 −11 , Ra F ∼ 10 27 ), accurate asymptotic trends have to be established on a larger range than what we could cover here (E ∈ [3 × 10 −4 , 3 × 10 −5 ]), which is a computational challenge.
For longitudinal variations, in the weakly rotating regime, Figure 5c shows that it is hard to argue for a clear change in the efficiency of the longitudinal variations transfer as Ekman is reduced.This is especially true for the thermal winds solution (after the jump in panel (c)), for which the fraction of the longitudinal anomaly retrieved at the top of the ocean,   lon top/bot =  lon top ∕ lon bot , stays close to about 75%.For the zonal winds solution however,   lon top/bot decreases as we explore more planetary relevant regimes.In the weakly rotating regime and for the cases representing pure tidal heating (gray vertical line on Figure 5c),   lon top/bot decreases from 20% to 15.5% and 14.2% at E = [3 × 10 −4 , 1.5 × 10 −4 , 6 × 10 −5 ] respectively (Figure 5h).In the moderately rotating regime,   lon top/bot decreases from about 80% to 47% from E = 3 × 10 −4 to 6 × 10 −5 (Figure 5j).
Figure 6 shows heat flux profiles and maps for weakly rotating cases in the pure tidal heating scenario, for Dirichlet and Neumann BCs (ΔT* = 0.65 and q* = 0.87), at various Ekman numbers.As already said above, the latitudinal anomaly is transferred entirely up to the ice-ocean boundary for all Ekman numbers investigated (Figures 6a and 6d), and the ocean is in a polar cooling configuration.Figures 6b-6e, and 6f also shows that significant longitudinal variations are present at the top of the ocean, even if   lon top/bot decreases.For Dirichlet BCs,   lon top/bot decreases from 18% to 13.0% at E = [3 × 10 −4 , 6 × 10 −5 ] respectively.With Neumann BCs, it decreases from 18% at E = 3 × 10 −4 to 15% at E = 6 × 10 −5 (see also Figure 5h).Because of their computational cost, we did not perform additional cases at even smaller Ekman and higher Rayleigh.It remains to be investigated if the longitudinal anomaly reaches a plateau at even smaller Ekman numbers and higher Rayleigh, or if the heat flux at the top of the ocean becomes zonally symmetric.

Flow Speed
Typical convective flow speed in planetary regimes can be estimated using scaling laws for the convective Reynolds number, defined here as a dimensionless root-mean-squared (rms) velocity:   = √ 2∕ = rms∕ , where E k is the time-averaged, dimensionless kinetic energy from which the contribution of the axisymmetric mode is subtracted, and V s is the dimensionless fluid volume.In the transitional regime between RR and NR convection, there is no clear asymptotic scaling for the convective Reynolds number (Gastine et al., 2016) 2023) found, based on energetic constraints, that ocean currents in Europa should be of the order of a few centimeters per second, that is, intermediate values between our scaling and the RR scaling.

Ice Thickness
To illustrate the possible impact on ice shell equilibrium, two end-member scenarios are considered: (a) a radiogenic-dominated scenario (q* = 0) where tidal heating in the silicate mantle is negligible and (b) a tidally dominated scenario (q* = 0.87), where tidal heating significantly exceeds radiogenic heating.Oceanic heat flux maps used as inputs for each end-member scenario are represented in Figures 7a and 7d, where we focus on the zonal winds solution represented in Figure 3 (weakly rotating, Neumann, stress-free cases q* = 0 and q* = 0.87).We showed that relative latitudinal heat flux variations obtained from the simulations are robust for the range of parameters explored.Assuming that this pole-to-equator heat flux difference continues to hold for real Europa parameters (see Section 3.4), we compute oceanic heat flux maps by extracting the relative heat flux variations from the simulations, and multiplying them by realistic average heat fluxes for each end-member.We caution again that longitudinal variations could be erased once in planetary regimes and should hence be viewed as upper bounds.In scenario 1 (2), we assume 5 mW/m 2 (37 mW/m 2 ) as the average oceanic heat flux at the ice-ocean interface (see also Appendix B).
The resulting ice thickness is computed using the model of Nimmo et al. (2007), described in Text S2 in Supporting Information S1.In a nutshell, this model assumes a purely conductive ice shell, where the local ice thickness depends on the surface temperature (which varies with latitude), the temperature at the bottom of the shell, assumed to be at the melting temperature, the oceanic heat flux, and the volumetric tidal heat production in the ice.It neglects altogether heat diffusion in the horizontal (latitudinal and longitudinal directions), viscous relaxation, variations of the melting temperature with pressure, and latent heat due to phase change.The laterally varying fields used as inputs for the ice thickness model-oceanic heat flux, tidal strain rate and surface temperature-are represented in Figure B1 in Appendix B.
Both oceanic heating and internal ice heating vary laterally and could induce ice thickness variations.To estimate their relative contributions, we first focus on the ice thickness resulting from a heterogeneous oceanic heating, while the internal ice heating is taken homogeneous and equal to its average value (Figures 7b and 7e).In the scenario where radiogenic heating is dominant, the ice shell is thinner in the equatorial region due to higher oceanic heat flux and higher surface temperature.If tidal heating from the mantle is considered, this latitudinal difference is reduced or even reversed and longitudinal variations appear.Unlike latitudinal variations, longitudinal heat flux variations are sensitive to the choice of mechanical BC.If strong zonal winds develop, they could be swept out once planetary regimes are reached (Section 3.4).
Endogenic processes within the ice will add complexity to this picture.To illustrate this, we add lateral variations of internal ice heating, arising from the spatially variable tidal strain rate (Figure B1(a)).Figure 7c shows that the dominant radiogenic heating scenario leads to thick ice, between 21 and 37 km (average thickness of 25 km).The total internal ice heating is 5.6 × 10 11 W (equivalent to 19 mW/m 2 at the ice-ocean boundary), and is responsible for almost all of the ice thickness variations, which thus follow the tidal strain rate pattern within the ice.The fact that the oceanic heat flux is higher by a few mW m −2 at the equator leads, nevertheless, to slightly thinner ice at the equator.The dominant tidal heating scenario (Figure 7f) would lead to globally thinner ice, because of the higher oceanic heat flux, with an average ice shell thickness of 12 km and variations between 9 and 17 km.The Time-averaged heat flux at the top of the ocean for q* = 0 (pure radiogenic heating scenario) in the zonal winds solution (stress-free).The maps are averaged over 0.15 diffusive times, as in Figure 3.(b) Corresponding ice thickness maps assuming that the internal heating in the ice is laterally homogeneous (see Figure B1).(c) Ice thickness maps taking into account that the internal heating in the ice is laterally inhomogeneous.(d, e, f) Similar maps for the case q* = 0.87 (dominant tidal heating scenario).The black curves on panels (c) and (f) represent the zonally averaged ice thickness, to emphasize latitudinal variations.Given the idealized nature of our model, we warn the reader that neglected processes within the ice shell (e.g., viscous flow, convection) as well as feedback effects of phase changes on the ocean dynamics could act to smooth the thickness variations obtained here.
total internal ice heating is 1.6 × 10 11 W (5 mW/m 2 at the ice-ocean boundary).The ice thickness variations are almost entirely determined by the oceanic heat flux.Minimum thicknesses are found at high latitudes and at the equator and largest thicknesses at low latitudes (∼25°), at the longitudes of the sub and anti-Jovian points.The tidal heating within ice slightly reinforces the variations already induced by the oceanic heat flux.The thermal winds solution leads to similar conclusions (Figure S3 in Supporting Information S1).As discussed in Section 3.4, longitudinal variations obtained in Figure 7f represents upper bounds as they may not be maintained in planetary regimes.
The ice thickness maps of Figure 7 are used to illustrate which heating variations dominate the thermal equilibrium of a conductive ice shell (oceanic heating vs. tidal heating in the ice).Again, we warn the reader that they should not be considered as an accurate quantitative prediction of ice thickness variations on Europa.In the real system, many other processes within the ice shell as well as feedback effects on the ocean dynamics will likely act to smooth the thickness variations obtained here.Regarding endogenic ice processes, lateral thickness variations could be erased by viscous flow (Kamata & Nimmo, 2017;Nimmo, 2004), or by solid-state convection (e.g., Ashkenazy et al., 2018;McKinnon, 1999;Tobie et al., 2003).From Galileo limb profiles, Nimmo et al. (2007) argue that lateral shell thickness variations do not exceed 7 km on Europa.This is consistent with the idea that the thickness variations obtained here are upper bounds of what could be expected on Europa.
As discussed in Kamata and Nimmo (2017), melting and freezing of ice need to occur to actively maintain any ice thickness gradient against the lateral ice flow.Equation 6in Nimmo (2004) provides an estimate of the ice flow timescale assuming a Newtonian rheology: where η b is the basal viscosity of ice, Δρ is the density contrast between the ice shell and the ocean, R is the universal gas constant, T m is the melting temperature, T s is the surface temperature, D is the average ice thickness, E is the activation energy and λ is the wavelength of the topography.Using the same parameters as in the ice thickness model (Table S2 in Supporting Information S1), Δρ = 100 kg/m 3 , and λ = 2,450 km (pole-to-equator topography), we obtain τ flow = 115 Myr for a 12km-thick ice, and 13 Myr for a 25km-thick ice.For a steady-state to exist, the melting/freezing timescale needs to be equal to the ice flow timescale (Kamata & Nimmo, 2017;Shibley et al., 2023).The heat flux at the ice-ocean boundary that would provide energy to sustain this melt rate can be estimated as q melt ∼ ρ i LΔH/τ flow , where ρ i = 920 kg m −3 is density of ice, L = 330 kJ kg −1 is latent heat of ice and ΔH is the ice thickness contrast.Assuming a typical pole-to-equator change of ice thickness ΔH ∼ 5 km, we obtain q melt ∼ 0.4 mW m −2 for a 12km-thick ice, and q melt ∼ 3.8 mW m −2 for a 25 km-thick ice, which is not unrealistic.That being said, given the huge uncertainty in the basal ice viscosity, the ice flow could be an order of magnitude faster or slower (Nimmo, 2004).Hence, a 5 km pole-to-equator thickness contrast could either be impossible to sustain due to a fast viscous relaxation, or be effectively sustained by corresponding melting/freezing at the ice-ocean boundary.Note that the Newtonian rheology underestimates the flow timescale compared to more realistic rheologies, especially for long wavelength topographies (Nimmo, 2004), supporting further the possibility of a steady-state.
Regarding feedback effects on the ocean dynamics, an important missing process in our model is the phase change at the ice-ocean boundary.The freezing point depression with pressure, and therefore with ice thickness, could induce a feedback on the underlying convection.Ice thickness variations of 5 km on Europa would lead to lateral temperature variations of about 460 mK at the ice-ocean boundary (Equation 5 in Kang, Mittal, et al. (2022)).These lateral temperature variations could force horizontal convection in the ocean, in a similar fashion as lateral heat flux variations at the seafloor.As investigated by Kang (2023), the formation of warm water under thin ice regions, and cold water under thick ice regions, may build up a stratified layer underneath the ice and prevent the bottom heat flux pattern from reaching the ice-ocean boundary (see also Zhu et al. (2017)).As discussed above, any ice thickness variation would need to be actively maintained by freezing and melting to compensate the viscous flow, and the associated brine and fresh-water injections could then drive its own large-scale circulation (Kang, Mittal, et al., 2022;Lobo et al., 2021).Including salinity effects also adds the possibility of double-diffusive convection (Vance & Brown, 2005;Wong et al., 2022).The large-scale topography due to ice thickness variations could also induce its own feedback on the ocean dynamics.Finally, mechanically driven and electro-magnetically-driven flows could also arise, but how they would compete with buoyancy-driven flows remains an open problem (Soderlund et al., 2024, and references therein).

Conclusion
In this study, we investigate if rotating thermal convection in Europa's ocean can transfer the heterogeneous tidal heating pattern at the seafloor up to the ice-ocean boundary.We use an idealized model where we neglect salinity and the possible feedback of phase change and topography at the ice-ocean boundary.Under these assumptions, we found that latitudinal heat flux variations from the seafloor are essentially mirrored at the ice-ocean boundary.Hence, if tidal heating in the silicate mantle is dominant compared to radiogenic heating, it would impose a polar cooling configuration for Europa (higher oceanic heat flux at the poles), independently of the convective regime considered.Unlike radiogenic heating, tidal heating in the silicate mantle of Europa could also participate in producing longitudinal heat flux variations.The amplitude of these longitudinal variations is however strongly dependent on the ability of thermal winds to develop in planetary regimes, and they could be smoothed, particularly in the presence of zonal winds.Quantifying if longitudinal variations persist in more extreme regimes than those investigated here remains a computational challenge.
Using the conductive ice thickness model of Nimmo et al. (2007), we conclude with two end-member scenarios.
If radiogenic heating is dominant in the silicate mantle, the heat flux variations from the ocean would be too weak to influence the thermal equilibrium of the ice shell, which would be mostly controlled by internal tidal heating in the ice.If tidal heating in the silicate mantle is in excess of the radiogenic contribution, the oceanic heat flux would be high enough to control the ice thickness.The ocean heat flux would at least reflect latitudinal variations of heat flux at the seafloor, leading to relatively thin polar ice.
As detailed in the previous section, the present ocean model is missing important physical processes.Coupled models able to take into account simultaneously the phase change and the nonlinear ocean dynamics driven by various forcing mechanisms are sorely needed to investigate feedback effects.The freezing point depression, but also the large-scale topography at the ice-ocean boundary should be incorporated in future models of buoyancy-driven ocean circulation.Furthermore, tides, libration and precession can excite a variety of waves in the ocean (see the review by Soderlund et al. (2024), and references therein).The nonlinear interaction of these waves can trigger boundary or bulk-filling turbulence, and induce dissipation which could contribute to the heat budget of the ocean.The feedback between convection and mechanically driven flows deserves further attention.
Studying coupled ice and ocean processes constitute exciting challenges for planetary modelers, both from the fundamental perspective and to prepare for interpretation of Europa Clipper (Phillips & Pappalardo, 2014;Roberts et al., 2023) and JUICE (Grasset et al., 2013) observations.Future observations from Europa Clipper may be able to disentangle between the two scenarios proposed here, (a) radiogenic-dominated heating in the silicate mantle, resulting in a thick ice shell, with a higher equatorial oceanic heat flux and (b) tidally dominated heating in the silicate mantle, resulting in a thin ice shell, with a higher polar oceanic heat flux.This might be achieved by using motional magnetic induction (Vance et al., 2021), constraining the depth of the ice-ocean interface and thermal variations in the ice using the ice-penetrating radar (Blankenship et al., 2018;Kalousová et al., 2017), constraining the ice shell thickness using gravity measurements (Mazarico et al., 2023), or estimating the asymmetry of the ocean thickness using the induced magnetic field (Styczinski et al., 2022).

Appendix A: Icy Moon Parameters and Regimes of Rotating Convection
Table A1 contains the dimensional and dimensionless physical parameters used to plot the regime diagram of Figure 2. The range of Europa's plausible heat fluxes (6-46 mW m −2 ) is determined the following way: Assuming that the radiogenic heating in Europa's mantle is of 6-7 mW m −2 at the seafloor (Tobie et al., 2003), the lower bound corresponds to the case of a negligible tidal heating in the silicate mantle.The upper bound corresponds to the upper bound of the log-normal distribution considered by Howell (2021), that is, a tidal heat flux of about 10 −1.5 ≈ 32 mW m −2 at the base of the ice (see their Figure 3).At the seafloor, this would lead to a heat flux of 39 mW m −2 , to which a radiogenic heating of 7 mW m −2 is added.This upper bound accounts for the possibility of Io-like tidal dissipation in Europa's mantle (Howell, 2021), that is, a total dissipated power of 1000 GW.The other parameters are taken from Soderlund (2019).For Europa, with the average parameters and their bounds listed in Table A1, we obtain the dimensionless parameters Pr ≈ 11, E = 7 × 10 −12 (∈ [5.0 × 10 −12 , 9.1 × 10 −12 ]) and Ra F = 5.4 × 10 27 (∈ [6.3 × 10 26 , 2.8 × 10 28 ]).
Systematic studies of rotating thermal convection in spherical shells have shown that various regimes of convection can take place, depending on the degree of rotational influence.Figure 2  Note.The dimensional physical parameters are taken from Soderlund (2019) except for Europa's heat flux: g is the gravitational acceleration, Ω the rotation rate, ν the kinematic viscosity, κ the thermal diffusivity, R the moon's radius, D ice the ice thickness, D ocean the ocean thickness, q the heat flux per unit area, ρ the ocean density, c p the ocean's heat capacity, and α the ocean's thermal expansivity.For the dimensionless parameters, η is the aspect ratio, Pr the Prandlt number, E the Ekman number, Ra F the flux Rayleigh number.In global circulation models (GCM) where small-scale turbulence is unresolved, its effect is parametrized by using effective viscosities and diffusivities which are orders of magnitude larger than molecular values, and often different in horizontal and vertical directions (e.g., Kang, 2023).(viscous and thermal diffusion) and too strongly forced.Since the dimensional heat flux imposed in the DNS is unphysically high, using the DNS heat flux in combination with realistic ice parameters is inconsistent and would lead to unphysically thin ice.In the ice thickness model, we therefore use the DNS relative heat flux variations, multiplied with a realistic mean heat flux (see Figures B1c-B1f) to obtain realistic dimensional ice thicknesses.This procedure implicitly assumes that the relative variations will stay the same when reaching planetary regimes (smaller Ekman, higher Rayleigh).Section 3.4 shows that for the range of parameters studied here, this is true for latitudinal variations, but longitudinal variations could be The longitudinal variations presented here should hence be considered as upper bounds.

Table A1 Dimensional and Dimensionless Physical
In the radiogenic-dominated scenario, an average heat flux of 5 mW m −2 at the top of the ocean is used.The corresponding heat flux maps are represented on Figure B1c and B1e.At the top of the ocean, the heat flux would globally vary between 4 and 7 mW/m 2 , with a maximum at the equator, and no longitudinal variations.In the tidally dominated scenario, we use the upper bound of 46 mW m −2 at the bottom of the ocean (Table A1), that is, 37 mW/m 2 at the top of the ocean.Ojakangas and Stevenson (1989).Average value   2  = 3.3 × 10 −20 s −2 .This value is used for the cases of Figure 7 where homogeneous internal ice heating is assumed.(b) Temperature at the surface of the ice calculated from Ojakangas and Stevenson (1989).(c), (e) Oceanic heat flux map estimated from the relative heat flux variations of the direct numerical simulations with q* = 0 (pure radiogenic heating) and stress-free (c) and no-slip (e) boundary conditions (BC).(d), (f) Oceanic heat flux estimated from simulations with q* = 0.87 (dominant tidal heating) and stress-free (d) and no-slip (f) BC.

Figure 1 .
Figure1.Patterns of heat flux per unit area at the seafloor in various scenarios.q* = Δq/q 0 , the ratio between the maximum heat flux difference and the laterally averaged heat flux, measures the heterogeneity of the heat flux.(a) If radiogenic heating is dominant in the silicate mantle, heating of the ocean is homogeneous (q* ≈ 0).(b) As the relative amplitude of tidal heating is increased, heat flux from the mantle becomes heterogeneous.For dominant tidal heating (c) relative heat flux variations become of order one (q* → 0.91).The tidal heating pattern is computed using the tidal model ofRoberts and Nimmo (2008) described in Text S1 in Supporting Information S1.(d) Heat flux pattern obtained at the end of Europa's thermal evolution model by Běhounková et al. (2021): the tidal heating represents about 30% of the total heating in the silicate mantle.Despite the occurrence of melting and mantle convection, the large-scale tidal heating pattern is persistent.Figure modified from Figure 2d in Běhounková et al. (2021), using the corresponding data set (Běhounková et al., 2020) licensed under CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/).

Figure 2 .
Figure 2. Regime diagram adapted fromGastine et al. (2016) representing the supercriticality, Ra T E 3/4 , as a function of the Ekman number, for a fixed Prandtl number Pr = 1 (Direct Numerical Simulations) and Pr = 10 (icy moons).The gray shaded area represents a region where the flow is stable, Ra T < Ra c (stability threshold determined byGastine et al. (2016) in the absence of horizontal forcing variations and for a shell aspect ratio of 0.6).Above that threshold, the flow is convectively unstable, and becomes fully non-linear once Ra T ≳ 6Ra c .Different regimes of convection are observed depending on the strength of the forcing compared to rotational effects.Each regime threshold is indicated by a black continuous line, and is taken from references GA16(Gastine et al., 2016), KC22(Kvorka & Čadek, 2022), EN14(Ecke & Niemela, 2014) and JU12(Julien et al., 2012).At a fixed E, as Ra T is increased, the convection is more vigorous, and the influence of rotation decreases, as indicated by the Geostrophic, Transitional, and Non-rotating regimes.The shaded blue region represents the region whereKvorka and Čadek (2022) identified a polar cooling behavior, that is, a higher heat flux at the pole.The dashed gray lines are lines of constant convective Rossby number   = 2576604x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023AV000994by University Of St Andrews University, Wiley Online Library on [19/12/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Figure 3 .
Figure 3. Mean oceanic circulation in a weakly rotating convective regime, with no-slip and stress-free mechanical boundary conditions (BC).The bottom heat flux is imposed to reproduce the silicate mantle heating pattern for increasing relative amplitude of tidal heating (Figure 1).All velocity maps correspond to the ocean mid-depth, and are averaged over 0.15 diffusive times to reveal the mean circulation, after a quasi steady-state has been reached.(a), (d): radial velocity.(b), (e): azimuthal velocity.(c), (f): azimuthal velocity from which the zonal flow was subtracted:   ′ =  − ⟨⟩ .The velocity is provided in dimensionless form as global Rossby number (D ∼ 111 km is the ocean thickness and Ω = 2.1 × 10 −5 s −1 the rotation rate).Note the different color scales for the azimuthal and radial flow, and for the two BC.Extrapolation of the typical convective flow speed to Europa's parameters is discussed in Section 3.4.2.

Figure 4 .
Figure 4. Relative heat flux variations measured at the top and at the bottom of the ocean for the zonal winds (first column) and the thermal winds (second column) solutions.(a), (b) Azimuthally averaged heat flux as a function of latitude.The dashed lines show the basal heat flux whereas the full lines are the top heat flux.From black to yellow, q* is increased from 0 to 2.16 (see Figure3).The two black lines show the tangent cylinder latitude (≈26°).(d), (e) Same data but for the latitudinally averaged heat flux as a function of longitude.(c) Latitudinal cooling parameter (Equation9) at the top versus at the bottom of the ocean.When q h/l < 0, the ocean is in an equatorial cooling configuration (higher heat flux at the equator).If the points lie around the 1:1 line, it means that all of the basal heat flux anomaly in latitude is retrieved at the top, and that the ocean efficiently transfers heat heterogeneities upwards.The colors have the same meaning as in other panels.(f) Longitudinal heat flux anomalies (Equation10) at the top versus at the bottom of the ocean.

Figure 5 .
Figure 5. Systematic for stress-free boundary conditions (BC) (BCs) in weakly rotating (a)-(c) and moderately rotating (d)-(f) regimes.Circles: Dirichlet bottom BC, Triangles: Neumann bottom BC.Colors correspond to the colors of the direct numerical simulations stars in Figure 2. (a), (d) q h/l at the top and at the bottom of the ocean.(b), (e) Anomaly in latitude at the top and at the bottom of the ocean.K lat is defined as K lon (Equation 10), except with a longitudinal average instead of the latitudinal average.(c), (f) Anomaly in longitude at the top and at the bottom of the ocean.In panel (c), the jump in the longitudinal anomaly at the top of the ocean is due to a transition from the zonal winds to the thermal winds solution.The transition occurs for all Ekman numbers investigated.In all panels, the vertical gray line represent bottom anomalies for pure tidal heating in the silicate mantle.(g)-(j) Fraction of the anomaly retrieved at the top of the ocean for simulations with the zonal winds solution.Open symbols: Dirichlet bottom boundary condition.Filled symbols and dashed line: Neumann bottom boundary condition.(g) Weakly rotating regime, latitudinal anomaly.(h) Weakly rotating regime, longitudinal anomaly.(i) Moderately rotating regime, latitudinal anomaly.(j) Moderately rotating regime, longitudinal anomaly.

Figure 6 .
Figure 6.Heat flux at the top of the ocean for simulations with stress-free boundary conditions and a tidally dominant scenario (ΔT* = 0.65 and q* = 0.87).Top row: Dirichlet bottom boundary condition.Bottom row: Neumann bottom boundary condition.(a), (d) Latitudinal heat flux profiles at the bottom (dashed lines) and at the top (full lines) of the ocean.〈•〉 ϕ denotes average in longitude and 〈•〉 h denotes horizontal average on the sphere.(b), (e) Longitudinal heat flux profiles at the bottom (dashed lines) and at the top (full lines) of the ocean.〈•〉 θ denotes average in latitude.(c), (f) Normalized heat flux maps at the top of the ocean.

Figure 7 .
Figure 7. Heat flux at the ice-ocean boundary in the zonal winds solution and corresponding ice thickness assuming a purely radially conducting ice shell.(a) Time-averaged heat flux at the top of the ocean for q* = 0 (pure radiogenic heating scenario) in the zonal winds solution (stress-free).The maps are averaged over 0.15 diffusive times, as in Figure3.(b) Corresponding ice thickness maps assuming that the internal heating in the ice is laterally homogeneous (see FigureB1).(c) Ice thickness maps taking into account that the internal heating in the ice is laterally inhomogeneous.(d, e, f) Similar maps for the case q* = 0.87 (dominant tidal heating scenario).The black curves on panels (c) and (f) represent the zonally averaged ice thickness, to emphasize latitudinal variations.Given the idealized nature of our model, we warn the reader that neglected processes within the ice shell (e.g., viscous flow, convection) as well as feedback effects of phase changes on the ocean dynamics could act to smooth the thickness variations obtained here.
the values of b and β given in Kvorka and Čadek (2022): b ≈ 14.73 and β ≈ 0.525.a Dimensionless parameters are based on molecular values of viscosity and diffusivity.
Parameters Used to Plot the Regime Diagram of Figure 2 2576604x, 2023, 6, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023AV000994by University Of St Andrews University, Wiley Online Library on [19/12/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Figure B1 .
Figure B1.Inputs of the ice thickness model.(a) Ice tidal strain rate calculated fromOjakangas and Stevenson (1989).Average value   2  = 3.3 × 10 −20 s −2 .This value is used for the cases of Figure7where homogeneous internal ice heating is assumed.(b) Temperature at the surface of the ice calculated fromOjakangas and Stevenson (1989).(c), (e) Oceanic heat flux map estimated from the relative heat flux variations of the direct numerical simulations with q* = 0 (pure radiogenic heating) and stress-free (c) and no-slip (e) boundary conditions (BC).(d), (f) Oceanic heat flux estimated from simulations with q* = 0.87 (dominant tidal heating) and stress-free (d) and no-slip (f) BC.

.
Re c should nevertheless be bounded by the NR and RR predictions, Note that our fit does not include any dependence on Pr because we did not vary this parameter.Using the parameters for Europa listed in TableA1, the predicted convective Reynolds are temperature Rayleigh numbers from the flux Rayleigh number, as defined in Appendix A. For