Ocean dynamics of outer solar system satellites

TL;DR

This paper uses theoretical and numerical methods to analyze ocean currents and heat transfer in outer solar system satellites, revealing how rotation influences ocean dynamics and potential geological activity.

Contribution

It provides a comprehensive prediction of ocean circulation regimes in satellites like Enceladus, Titan, Europa, and Ganymede based on rotational influence, linking ocean dynamics to surface geology.

Findings

Enceladus and Titan likely have multiple zonal jets and high-latitude heat transfer.

Europa's ocean may feature Hadley-like circulation cells driving geologic activity.

Ganymede's ocean exhibits multiple circulation regimes.

Abstract

Ocean worlds are prevalent in the solar system. Focusing on Enceladus, Titan, Europa, and Ganymede, I use rotating convection theory and numerical simulations to predict ocean currents and the potential for ice-ocean coupling. When the influence of rotation is relatively strong, the oceans have multiple zonal jets, axial convective motions, and most efficient heat transfer at high latitudes. This regime is most relevant to Enceladus and possibly to Titan, and may help explain their long-wavelength topography. For a more moderate rotational influence, fewer zonal jets form, Hadley-like circulation cells develop, and heat flux peaks near the equator. This regime is predicted for Europa, where it may help drive geologic activity via thermocompositional diapirism in the ice shell, and is possible for Titan. Weak rotational influence allows concentric zonal flows and overturning cells with no preferred orientation. Predictions for Ganymede's ocean span multiple regimes.

Tables

Table 1. Table 1: Properties of icy satellite oceans in dimensional units and non-dimensional parameters. Kinematic viscosity from Nayar  \BOthers . ( \APACyear 2016 ) , thermal diffusivity from Abramson  \BOthers . ( \APACyear 2001 ) , interior model properties ( R S , D I ​ h , D o ​ c ​ e ​ a ​ n , q subscript 𝑅 𝑆 subscript 𝐷 𝐼 ℎ subscript 𝐷 𝑜 𝑐 𝑒 𝑎 𝑛 𝑞 R_{S},D_{Ih},D_{ocean},q ) from Vance  \BOthers . ( \APACyear 2018 ) , and thermodynamic properties ( ρ , C p , α 𝜌 subscript 𝐶 𝑝 𝛼 \rho,C_{p},\alpha ) at the mean ocean temperatures and pressures of the respective interior models from Vance ( \APACyear 2017 ) for MgSO 4 (0 and 10 wt %) and B. Journaux (personal communication) for seawater. Two outer ice shell thicknesses and three ocean compositions are considered for Enceladus and Europa, while three outer ice shell thicknesses and two ocean compositions are considered for Titan and Ganymede.

	Enceladus	Titan	Europa	Ganymede
Gravitational acceleration, $g$ [m/s2]	$0.1$	$1.4$	$1.3$	$1.4$
Rotation rate, $\Omega$ [s-1]	$5.3\times {10}^{-5}$	$4.6\times {10}^{-6}$	$2.1\times {10}^{-5}$	$1.0\times {10}^{-5}$
Kinematic viscosity, $\nu$ [m2/s]	$1.8\times {10}^{-6}$	$1.8\times {10}^{-6}$	$1.8\times {10}^{-6}$	$1.8\times {10}^{-6}$
Thermal diffusivity, $\kappa$ [m2/s]	$1.4\times {10}^{-7}$	$1.8\times {10}^{-7}$	$1.6\times {10}^{-7}$	$1.8\times {10}^{-7}$
Satellite radius, $R_S$ [km]	$252$	$2575$	$1561$	$2631$
Ice Ih thickness, $D_{Ih}$ [km]
Water	$51,10$	$141,74,50$	$30,5$	$134,70,5$
MgSO4 10 wt%	$50,10$	$149,86,58$	$30,5$	$157,95,26$
Seawater	$50,10$	$-$	$30,5$	$-$
Ocean thickness, $D_{ocean}$ [km]
Water	$11,53$	$130,369,420$	$97,124$	$119,361,518$
MgSO4 10 wt%	$13,63$	$91,333,403$	$103,131$	$24,287,493$
Seawater	$12,55$	$-$	$99,126$	$-$
Heat flux, $q$ [mW/m2]
Water	$16,81$	$14,18,20$	$24,119$	$16,20,107$
MgSO4 10 wt%	$16,83$	$14,17,19$	$24,123$	$15,18,25$
Seawater	$16,82$	$-$	$23,121$	$-$
Density, $\rho$ [${10}^3$ kg/m3]
Water	$1.00,1.00$	$1.11,1.14,1.14$	$1.04,1.04$	$1.11,1.14,1.14$
MgSO4 10 wt%	$1.11,1.11$	$1.20,1.23,1.24$	$1.15,1.14$	$1.19,1.23,1.24$
Seawater	$1.02,1.02$	$-$	$1.07,1.07$	$-$
Specific heat capacity, $C_p$ [${10}^3$ J/kg/K]
Water	$4.2,4.2$	$3.0,3.5,3.6$	$3.9,3.9$	$3.0,3.5,3.7$
MgSO4 10 wt%	$3.6,3.7$	$2.1,2.5,2.8$	$3.3,3.5$	$2.1,2.4,3.0$
Seawater	$4.0,4.0$	$-$	$3.8,3.8$	$-$
Thermal expansivity, $\alpha$ [${10}^{-4}$ K-1]
Water	$-0.5,-0.5$	$2.3,4.0,4.2$	$1.9,1.9$	$2.2,4.0,4.4$
MgSO4 10 wt%	$1.2,1.3$	$0.4,2.1,2.7$	$2.1,2.3$	$-0.1,1.9,3.2$
Seawater	$0.1,0.1$	$-$	$2.5,2.5$	$-$
Prandtl number, $Pr=\nu /\kappa$	13	10	11	10
Ekman number, $E=\nu /\Omega D^2$	${10}^{-10}-{10}^{-11}$	${10}^{-11}-{10}^{-12}$	${10}^{-12}$	${10}^{-10}-{10}^{-13}$
Rayleigh number, $Ra=\alpha g\Delta T{D}^3/\nu \kappa$	${10}^{16}-{10}^{19}$	${10}^{19}-{10}^{23}$	${10}^{20}-{10}^{22}$	${10}^{20}-{10}^{24}$
Radius ratio, $\chi =\frac{\left(R_S-D_{Ih}\right)}{\left(R_S-D_{Ih}-D_{ocean}\right)}$	$0.74-0.95$	$0.83-0.96$	$0.92-0.94$	$0.80-0.99$

Equations

$$\Delta T=7.3{\left(\frac{\nu }{\alpha g\rho C_p}\right)}^{1/4}{q}^{3/4}$$

$$\Delta T=2.1{\left(\frac{{\Omega }^4\kappa }{{\rho }^2{C}_p^2\nu {\alpha }^3{g}^3}\right)}^{1/5}(q^{2}D{)}^{1/5}$$

Full text

Plain Language Summary: The outer solar system is host to a large number of diverse satellites, many of which likely have global oceans beneath their outer icy shells. I use theoretical arguments and numerical models to make predictions about ocean currents and heat transfer across such oceans. Our results suggest that strong ocean currents exist in Enceladus, Titan, Europa, and Ganymede, and cause the transfer of heat to vary with latitude that may modify the overlying ice shell.

\draftfalse\journalname

Geophysical Research Letters

Institute for Geophysics, Jackson School of Geosciences, The University of Texas at Austin, Austin, Texas, USA

\correspondingauthor

Krista Soderlundkrista@ig.utexas.edu

{keypoints}

Ocean dynamics are important for the habitability of icy ocean worlds.

Strong ocean currents likely exist in Enceladus, Titan, Europa, and Ganymede.

Convective heat transfer in the ocean is predicted to vary with latitude, which would modify the thermophysical structure of the ice shell.

1 Introduction

Exploration of the outer solar system has shown that subsurface oceans may be relatively common in the interiors of icy satellites and dwarf planets (Nimmo \BBA Pappalardo, \APACyear2016; Lunine, \APACyear2017). Strong evidence exists for oceans in the Saturnian satellites Enceladus and Titan, with oceans also potentially present in Mimas and Dione. In the Jovian system, there is compelling evidence for a Europan ocean and oceans are predicted within Ganymede and potentially Callisto as well (c.f. Hartkorn \BBA Saur, \APACyear2017). Kuiper belt objects, such as Pluto, Charon, and the Neptunian satellite Triton, may also have subsurface oceans.

The presence of liquid water makes these ocean worlds compelling astrobiological targets. However, the dynamics of these oceans also play a role in promoting habitable environments. For example, heat and material exchange between the seafloor and ice shell will be enhanced if the ocean is unstable to convection (e.g., Vance \BBA Goodman, \APACyear2009; Soderlund \BOthers., \APACyear2014), has mechanically driven flows (e.g., Tyler, \APACyear2008; Lemasquerier \BOthers., \APACyear2017; Wilson \BBA Kerswell, \APACyear2018; Rovira-Navarro \BOthers., \APACyear2019), or has fluid motions driven by magnetic pumping (Gissinger \BBA Petitdemange, \APACyear2019). Currents and turbulence tend to mix the ocean waters, which implies the presence of strong thermal and compositional gradients near the top and bottom of the ocean that life may take advantage of. Mixing efficiency may vary spatially, so these motions are also important for the distribution of potential bionutrients. In addition, heat transfer from the ocean will influence where the ice shell melts and freezes. Melting of the ice shell and freezing of the ocean will impact the salt budget, especially near the ice-ocean interface, a habitable environment in analogous terrestrial ice shelves (e.g., Daly \BOthers., \APACyear2013). Moreover, accreted ice depleted in salts may have positive buoyancy and lead to upwelling thermocompositional diapirs in the ice shell that would link the surface and subsurface (e.g., Pappalardo \BBA Barr, \APACyear2004; Soderlund \BOthers., \APACyear2014).

Here, I focus on the ocean dynamics of Enceladus and Titan given the abundance of data from the Cassini mission and of Europa and Ganymede in preparation for the upcoming Europa Clipper (Phillips \BBA Pappalardo, \APACyear2014) and JUICE (Grasset \BOthers., \APACyear2013) missions. Scaling laws for rotating convection systems are used to make predictions about their convective behaviors in Section 2, and numerical models of global ocean convection characterizing the predicted regimes are presented in Section 3. Implications for icy satellites are explored in Section 4, and the challenges of extrapolating to realistic ocean conditions are discussed in Section 5.

2 Rotating Convection Scaling Laws

Convection characteristics depend critically on the relative importance of rotation, which tends to organize the fluid into columns aligned with the rotation axis, increase the critical Rayleigh number, constrain heat transfer efficiency, and drive zonal flows (e.g., Aurnou \BOthers., \APACyear2015). Cheng \BOthers. (\APACyear2018) combines asymptotic predictions, laboratory experiments, and numerical simulations to review the behavior of rotating thermal convection as a function of the dimensionless Ekman, Rayleigh, and Prandtl numbers. The Ekman number, $E=\nu/2\Omega D^{2}$, represents the ratio of rotational to viscous timescales; thus, low $E$ signifies rapid rotation rates in planetary interiors. The Rayleigh number, $Ra=\alpha g\Delta TD^{3}/\nu\kappa$, is the ratio of the thermal diffusion time to the viscous buoyant rise time; large $Ra$ denotes strong buoyancy forcing. The Prandtl number, $Pr=\nu/\kappa$, defines the ratio of thermal to viscous diffusion timescales. Here, $\nu$ is kinematic viscosity, $\Omega$ is rotation rate, $D$ is ocean thickness, $\alpha$ is thermal expansivity, $g$ is gravitational acceleration, $\Delta T$ is superadiabatic temperature contrast, and $\kappa$ is thermal diffusivity.

Cheng \BOthers. (\APACyear2018) identify five rotating convection regimes: columnar, plumes, geostrophic turbulence (GT), unbalanced boundary layer (UBL), and nonrotating heat transfer (NR) (see Fig. 1). Near onset, convection in the bulk fluid manifests as Taylor columns aligned with the rotation axis (“columnar” regime). With increased buoyancy forcing, the columns begin to deteriorate such that they no longer extend fully across the fluid layer (“plumes” regime). Convection eventually becomes vigorous enough for strong mixing in the bulk fluid (“geostrophic turbulence” regime). Despite the disappearance of coherent vertical structures, the Coriolis force still imposes a vertical stiffness on the flow field. These regimes are shown collectively on Figure 1. The influence of rotation is lost locally at Rayleigh numbers exceeding $Ra_{GTU}$, which corresponds to the breakdown of geostrophy (balance between Coriolis and pressure gradient forces) in the thermal boundary layers (“unbalanced boundary layers” regime). For Rayleigh numbers greater than $Ra_{UNR}$, the influence of rotation is lost globally (“nonrotating heat transfer” regime). As reviewed by Cheng \BOthers. (\APACyear2018), significant debate exists in the community on the scaling laws for $Ra_{GTU}$ and $Ra_{UNR}$. Rather than assume a single scaling law for each transition, I consider upper and lower bound scaling laws for each regime and highlight the resulting range of parameter space for each regime transition in Figure 1.

Using this regime diagram, one can predict the convective regime of a system if the Ekman, Rayleigh, and Prandtl numbers can be estimated (see Table LABEL:tab:physparam). The Prandtl number depends only on fluid properties and is estimated to be $Pr\sim 10$ for the satellite oceans (Abramson \BOthers., \APACyear2001; Nayar \BOthers., \APACyear2016). The Ekman number is also relatively easy to calculate since it only requires assumptions about the fluid viscosity, rotation rate, and ocean thickness. I use the internal structure models of Vance \BOthers. (\APACyear2018) (see their Tables 5-8) to obtain ocean thicknesses $D_{ocean}$ for six combinations of possible outer ice shell thicknesses and ocean compositions for each satellite. Enceladus’ ocean has the largest Ekman number of $E\sim\mathcal{O}(10^{-10})$, while the ocean of Ganymede has the lowest at $E\sim\mathcal{O}(10^{-13})$.

The Rayleigh number is more difficult to estimate because it requires knowledge of the superadiabatic temperature contrast $\Delta T$. One can derive an estimate, however, using the relationship between the Rayleigh number and the convective heat transfer efficiency as measured by the Nusselt number, $Nu=qD/\rho C_{p}\kappa\Delta T$. Following Soderlund \BOthers. (\APACyear2014), I leverage $Nu-Ra$ scalings to solve for $\Delta T$ algebraically and consider both non-rotating and rapidly rotating scaling laws to give end-member estimates. More recent scaling laws for rotating spherical shells are used here, however. In the non-rotating regime, heat transfer is expected to be independent of the Ekman number and follow the theoretical limit of $Nu=0.07Ra^{1/3}$ (e.g., Gastine \BOthers., \APACyear2015). Conversely, in the rapidly rotating limit, heat transfer is predicted to follow $Nu=0.15Ra^{3/2}(2E)^{2}$ (Gastine \BOthers., \APACyear2016). As a result, the temperature contrast is given by

[TABLE]

in the non-rotating regime and by

[TABLE]

in the rapidly-rotating regime.

Here, I assume the heat flux $q$ from each of the six interior models per satellite, noting that the lower $q$ estimates are associated with thicker ice Ih shells (Vance \BOthers., \APACyear2018). Although these values are generally consistent with the literature, the minimum heat fluxes tend to exceed those predicted for radiogenic heating in the mantle at present day (e.g., Bland \BOthers., \APACyear2009) and the upper bound for Ganymede is appropriate for a past active period (e.g., Dombard \BBA McKinnon, \APACyear2001). If $q$ is decreased by an order of magnitude, the lower bounds for $\Delta T$, and therefore $Ra$, only decrease by a factor of 2.5 per eqn. (2). Our global estimates also neglect spatial variations in heat flow that may be locally strong at Enceladus (Choblet \BOthers., \APACyear2017), for example, where narrow mantle upwellings can reach $1-5$ W/m2 (the global average, however, is in line with Vance \BOthers. (\APACyear2018)). If $q$ is increased by an order of magnitude, the $Ra$ upper bounds increase by a factor 5.6 per eqn. (3). As shown in Table LABEL:tab:physparam, Rayleigh numbers span from $Ra\sim\mathcal{O}(10^{16})$ for the lower Enceladus limit to $Ra\sim\mathcal{O}(10^{24})$ for the upper Ganymede limit. An important caveat to note here, however, is that these estimates do not include compositional contributions due to salinity gradients. This simplification may be especially significant for Titan since the ocean is hypothesized to have a high concentration of dissolved salts (Baland \BOthers., \APACyear2014; Mitri \BOthers., \APACyear2014).

Figures 1 and S1 plot the resulting estimates of the Ekman and Rayleigh numbers on the convective regime diagram. The oceans of Titan, Europa, and Ganymede are predicted to behave similarly since their estimated parameter spaces have considerable overlap. Since these estimates fall near the lower boundary between the UBL and NR regimes, I hypothesize that rotational effects do not dominate the turbulent local-scale convective flows. Conversely, rotation likely has a stronger influence on the ocean of Enceladus, which is also predicted to be primarily in the UBL regime, although extending into the GT transition.

3 Numerical Convection Models

Numerical models of global ocean convection are next used to characterize the currents and heat flow patterns. I utilized the pseudospectral code MagIC, version 5.6 (e.g., Wicht, \APACyear2002; Gastine \BBA Wicht, \APACyear2012) to simulate 3D, time-dependent, thermal convection of a Boussinesq fluid in a rotating spherical shell with geometry characterized by the ratio of inner to outer shell radii, $\chi=r_{i}/r_{o}=0.9$. The system is further defined by the Ekman, Rayleigh, and Prandtl numbers. Following Soderlund \BOthers. (\APACyear2014), the boundaries are impenetrable, stress-free, and isothermal. Compositional buoyancy, spatial variations in mantle heat flow, and mechanically driven flows are neglected for simplicity.

Seven models that span a convective regime space consistent with the icy satellite ocean predictions are considered (Fig. 1). In the first series, the Rayleigh and Prandtl numbers are fixed to $Ra=3.4\times 10^{7}$ and $Pr=1$, and the Ekman number is increased from $E=3.0\times 10^{-5}$ to $E=7.5\times 10^{-4}$. The second series of models increase the Rayleigh number from $Ra=2.4\times 10^{6}$ to $Ra=3.4\times 10^{7}$ for fixed $E=1.5\times 10^{-4}$ and $Pr=1$; higher $Ra$ values were not pursued to due computational limitations. Hyperdiffusivities are not employed (c.f. Zhang \BBA Schubert, \APACyear2000). The numerical grids have 73 radial points, 320 latitudinal points, and 640 longitudinal points for cases with $E\geq 7.5\times 10^{-5}$ and 65 radial points, 640 latitudinal points, and 1280 longitudinal points for the $E=3.0\times 10^{-5}$ case. Each model was initiated with a random temperature perturbation or restarted from a lower $E$ or $Ra$ case.

Figure 2 shows the mean velocity and temperature fields of each model across the $E$ parameter sweep, while Figure 3 shows the normalized heat flux along the outer boundary. In the highest Ekman number case (Fig. 2a), the zonal and radial flows have comparable magnitudes reminiscent of non-rotating convection. The radial flows have no preferred spatial orientation, while the zonal flows are concentric due to viscous transport of angular momentum (Brun \BBA Palacios, \APACyear2009). Ocean temperatures are nearly isothermal away from the boundaries, leading to localized heat flux perturbations along the ice-ocean interface (Fig. 3a).

When the Ekman number is decreased (Fig. 2b-c), homogenization of absolute angular momentum leads to zonal flows that are retrograde (westward) at large cylindrical radii and prograde (eastward) closer to the rotation axis (e.g., Gilman, \APACyear1978; Aurnou \BOthers., \APACyear2007; Gastine \BOthers., \APACyear2013). The mean radial flows become more organized with a pronounced upwelling near the equator and downwellings at mid-latitudes, essentially forming Hadley-like meridional circulation cells in each hemisphere. Upon further decreasing of the Ekman number (Fig. 2d), multiple zonal jets that alternate in direction develop, and the mean radial flows retain an equatorial upwelling that becomes more aligned with the rotation axis. Both mean zonal and radial flow speeds decrease by a factor of five compared to the $E=[3.0,1.5]\times 10^{-4}$ cases. In all three of these models (Fig. 2b-d), ocean temperatures are characterized by thin thermal boundary layers and warmer equatorial waters. Heat flux peaks at low latitudes (with minima at mid-latitudes) due to the mean overturning circulations with secondary peaks forming at high latitudes due to turbulent heat transfer associated with vertically ascending plumes (Fig. 3b-d).

In the lowest Ekman number case (Fig. 2e), Coriolis forces organize the flow into narrow structures that are aligned with the rotation axis. Reynolds stresses associated with these columns drive prograde equatorial flow with jets that alternate in direction at higher latitudes due to correlation locally between the azimuthal and cylindrically radial flow components (e.g., Aurnou \BBA Olson, \APACyear2001; Christensen, \APACyear2001; Heimpel \BOthers., \APACyear2005; Gastine \BOthers., \APACyear2014). Ocean temperatures are not well-mixed, especially at low latitudes, due to the axialized convective flows and strong equatorial jet (e.g., Aurnou \BOthers., \APACyear2008). Consequently, heat flow along the ice-ocean interface peaks at high latitudes with minima near the equator (Fig. 3e).

A similar trend from three-jet zonal flows, equatorial upwelling, and peak low latitude heat flux to multiple zonal jets, axialized convective flows, and peak high latitude heat flux is found as $Ra$ is decreased (Figs. S2 and S3).

4 Implications for Icy Satellites

In order to apply these models to icy satellite oceans, I assume that the velocity and temperature patterns extrapolate to more extreme parameters following the relative distance between regime boundaries (Fig. S1). Enceladus’ ocean may then be represented by the $E=[3.0\times 10^{-5},7.5\times 10^{-5}]$ models since both fall approximately between the GT-UBL regime transition and the lower bound of the UBL-NR transition. In contrast, the oceans of Europa, Ganymede, and Titan depend on the UBL-NR transition scaling used. If $Ra^{RoC=1}_{UNR}=E^{-2}Pr$ (e.g., Gilman, \APACyear1977) is assumed, then all of these oceans are near the center of UBL regime such that the $E=[7.5\times 10^{-5},1.5\times 10^{-4}]$ models would be most appropriate for these satellites. If the transition instead follows $Ra^{Ga16}_{UNR}=100(2E)^{-12/7}$ (Gastine \BOthers., \APACyear2016), then the $E=3.0\times 10^{-4}$ model would best characterize Europa and the $E=[3.0\times 10^{-4},7.5\times 10^{-4}]$ models would be most pertinent to Titan and Ganymede.

Below, I discuss the implications for each satellite. Regions with high heat flow are presumed to undergo enhanced melting, leading to ice shell thickness variations. However, large thickness disparities can set up a phenomena known as an ice pump (e.g., Lewis \BBA Perkin, \APACyear1986) where pressure-induced melting occurs where the ice shell is thick and re-accretes where the ice shell is thin, effectively reducing topography along the base of the ice shell. Since the accretion process is very efficient at excluding impurities in low temperature environments (e.g., Moore \BOthers., \APACyear1994; Eicken \BOthers., \APACyear1984), this marine ice may be salt-depleted compared to the overlying ice. The ice may, therefore, have positive buoyancy due to the associated thermal and compositional density anomalies and rise toward the surface in the form of convective diapirs (e.g., Pappalardo \BBA Barr, \APACyear2004; Soderlund \BOthers., \APACyear2014). Alternatively, if the ice pump mechanism is not efficient, the ice shell may be more unstable to convection where it is relatively thick (Travis \BOthers., \APACyear2012; Goodman, \APACyear2014).

For Enceladus, I predict the zonal flows to be characterized by multiple jets that alternate in direction (Fig. 2A d-e). Converting model velocities to dimensional units $U=\Omega DRo$, I expect peak zonal speeds of nearly a m/s depending on the ocean thickness assumed. Meridional circulations are predicted to either be strongly aligned with the rotation axis with speeds up to a few mm/s (Fig. 2B e) or be concentrated in a low latitude upwelling with speeds up to a few cm/s (Fig. 2B d). As a result, heat flow along the ice-ocean interface has distinct peaks at either the poles (Figs. 2C e, 3e) or at the equator and the poles secondarily (Figs. 2C d, 3d).

Measurements of Enceladus’ shape, gravitational field, and librational motions show that the ice shell is thin below the south pole and thick at the equator, with an intermediate thickness at the north pole (e.g., Čadek \BOthers., \APACyear2016; Beuthe \BOthers., \APACyear2016). Inverting these measurements to infer the oceanic heat flux along the ice-ocean interface, Čadek \BOthers. (\APACyear2019) find peak flux near the poles with a minima at the equator. This pattern implies upwelling of warm water at the poles and downwelling of cool water at low latitudes, which may be caused by ocean convection (Fig. 2e) and/or be a consequence of the pattern of tidal heating in the mantle (Choblet \BOthers., \APACyear2017). Considering the former, the $E=3.0\times 10^{-5}$ model is appropriate if a low internal heat flux is assumed (in contrast to Choblet \BOthers., \APACyear2017) or if the thermal expansion coefficient in our calculations is overestimated since $\alpha$ trends towards zero with decreasing salinity and becomes negative for freshwater (e.g., Nayar \BOthers., \APACyear2016; Feistel, \APACyear2010, see also Table LABEL:tab:physparam), which both effectively reduce the Rayleigh number and make rotational effects more important.

Europa’s ocean is predicted to have three zonal jets with retrograde equatorial flow that can reach m/s speeds (Fig. 2A b-c) or multiple zonal jets with retrograde equatorial flow and reduced speeds (Fig. 2A d) depending on the scaling law assumed. All Europa-relevant models, however, have an equatorial upwelling of warm water with peak speeds of roughly a few cm/s and enhanced heat transfer at low latitudes (Figs. 2B-C b-d, 3b-d).

The surface of Europa is riddled with geologic features indicating recent activity and the potential for ocean-derived materials (e.g., Figueredo \BBA Greeley, \APACyear2003; Fischer \BOthers., \APACyear2015). Chaos terrains, for example, appear to be located preferentially at low latitudes with a secondary prevalence near the poles (Leonard \BOthers., \APACyear2018), and formation models suggest that they may be associated with upwelling diapirs (e.g., Sotin \BOthers., \APACyear2002; Collins \BBA Nimmo, \APACyear2009; Schmidt \BOthers., \APACyear2011) and marine ice accretion (Soderlund \BOthers., \APACyear2014). No large gradients in ice shell thickness have been detected (Nimmo \BOthers., \APACyear2007), suggesting an efficient ice pump (c.f. Nimmo, \APACyear2004). Given our robust model predictions of high oceanic heat flux at low latitudes with relatively low flux at mid-latitudes, our new calculations continue to support the thermocompositional diapirism hypothesis.

Given the similarities in regime predictions for Titan and Ganymede, they are considered together here. Assuming the $Ra^{Ga16}_{UNR}$ scaling and upper $Ra$ estimates, the oceans would behave akin to a non-rotating system with no coherent heat transfer patterns (Figs. 2 C a, 3a). If the lower $Ra$ estimates are instead assumed, these satellites are predicted to have three-jet zonal flows with peak speeds up to a few m/s (depending on ocean thickness), Hadley-like circulation cells with peak speeds up to tens of cm/s (depending on ocean thickness), and maximum heat transfer near the equator (Figs. 2 A-C b, 3b). Alternatively, for the $Ra^{RoC=1}_{UNR}$ scaling, these satellites are predicted to behave similarly to Enceladus and Europa as discussed above (Figs. 2 A-C c-d, 3c-d), except with respect to dimensionalized flow speeds that could be considerably faster due to the larger ocean thicknesses.

Looking to Titan, the satellite’s surface topography shows polar depressions compared to relatively elevated low latitudes (e.g., Durante \BOthers., \APACyear2019) that are likely explained by ice shell thickness variations (Nimmo \BBA Bills, \APACyear2010; Hemingway \BOthers., \APACyear2013; Lefevre \BOthers., \APACyear2014) or ice shell density variations (Choukroun \BBA Sotin, \APACyear2012). As for Enceladus, geophysical measurements by Cassini have been used to infer the oceanic heat flux along the ice-ocean interface (Kvorka \BOthers., \APACyear2018). The pattern is spatially complex, but simplifies to peaks near the poles when only axisymmetric components are considered. In contrast, the ocean convection models predicted to be relevant for Titan have either no coherent heat flux pattern (Fig. 3a), peak flux and enhanced melting near the equator (Fig. 3b-c), or peak flux near both the equator and the poles (Fig. 3d), none of which are consistent with the observed long-wavelength topography assuming Airy isostasy. If ocean dynamics alone are responsible, this difference implies that either (1) the melted equatorial region in the intermediate scenario was infilled with less dense marine ice to form the equatorial bulge through Pratt isostasy or (2) the ocean has a stably-stratified salinity gradient that reduces the effective buoyancy forcing of the ocean ($Ra$) such that rotational effects become sufficient to maximize heat flow and melting at the poles (Fig. 3e).

Observational constraints for Ganymede are much more limited. The satellite’s ancient grooved terrains indicate a likely period of geologic activity in its early history (Lucchita, \APACyear1980) and detection of hydrated salts suggests a subsurface briny layer of fluid (McCord \BOthers., \APACyear2001), but no clear patterns are present. Mass anomalies were measured in the northern hemisphere during a single Galileo flyby (Palguta \BOthers., \APACyear2006), but the sparsity of data prohibit both characterization on a global scale and unique determination of their depth of origin. Consequently, there is no clear link at present between observations and the underlying ocean dynamics.

5 Discussion

Our results are broadly consistent with the literature. Moreover, by comparing our numerical models against those with different input parameters, we are able to assess their sensitivity to these aspects. For example, the satellite oceans are predicted to have geometries characterized by $\chi$ values ranging from $0.74$ to $0.99$ (Table LABEL:tab:physparam), compared to our models with fixed $\chi=0.90$. Gastine \BOthers. (\APACyear2013) found that anelastic columnar convection in thicker spherical shell geometries ($\chi=0.6$) is also characterized by a prograde equatorial jet with multiple, small-scale meridional circulations aligned with the rotation axis, which transitions to a regime with a retrograde equatorial jet and Hadley cell-like meridional circulations and ultimately a weakening of zonal flow speeds as the influence of rotation is decreased. Similarly, Aurnou \BOthers. (\APACyear2008) showed that heat transfer is inhibited at low latitudes and generally increases towards the poles for columnar convection in spherical shells with $\chi=[0.85,0.9]$; this result is different from Miquel \BOthers. (\APACyear2018), who obtained peak equatorial heat transfer in their asymptotic models of rapidly rotating convection near onset, where the polar regions are subcritical (e.g., Dormy \BOthers., \APACyear2004). Conversely, Brun \BBA Palacios (\APACyear2009) and Soderlund \BOthers. (\APACyear2013) showed that the equatorial heat transfer enhancement for less vertically stiff convection is robust for thicker shells ($\chi\leq 0.75$) and different thermal boundary conditions. Simulations with thinner spherical shells are computationally demanding and uncommon (c.f. DeRosa \BOthers., \APACyear2002). Furthermore, studies with a thin layer of stable stratification below the outer boundary, which may be expected in regions where thermal expansivity is negative (Melosh \BOthers., \APACyear2004, see also Table LABEL:tab:physparam), generally show similar trends (e.g., Heimpel \BOthers., \APACyear2015). Thus, these studies suggest that our results are not strongly sensitive to variations in ocean thickness or fluid properties with depth. Large spatial variations in ice shell thickness may enhance mechanically driven flows though (e.g., Lemasquerier \BOthers., \APACyear2017), which are not considered here.

The effects of different boundary conditions should also be considered. For example, we assumed stress-free mechanical boundaries in order to reduce the effects of viscosity (e.g., artificially large Ekman boundary layers) due to the large $E$ values of our models compared to the satellites (Table LABEL:tab:physparam; Kuang \BBA Bloxham, \APACyear1997). In models with no-slip boundary conditions, inertial effects tend to be reduced substantially and strong zonal flows can be inhibited, which can disrupt convection (e.g., Aurnou \BBA Heimpel, \APACyear2004; Jones, \APACyear2015). For sufficiently driven and rapidly rotating convection, however, no-slip boundaries do not necessarily have this inhibiting effect (Manneville \BBA Olson, \APACyear1996; Aubert \BOthers., \APACyear2001). Uniform fixed temperature boundary conditions were assumed because they enable a broader comparison with the literature and across the satellites, although fixed heat flux boundary conditions may be more appropriate along the seafloor. At moderate parameters, fixed flux conditions tend to promote larger convective scales (Sakuraba \BBA Roberts, \APACyear2009; Hori \BOthers., \APACyear2012) and spatial variations along the boundary can influence the flow and efficiency of heat transfer, especially near the interface (e.g., Dietrich \BOthers., \APACyear2016; Mound \BBA Davies, \APACyear2017, for anomalies along the outer boundary). At extreme (i.e. realistic) parameters, however, the solutions for both thermal conditions appear to converge for rapidly rotating convection (Calkins \BOthers., \APACyear2015) as well as Rayleigh-Bénard convection (Johnston \BBA Doering, \APACyear2009). This convergence implies that convective-scale spatial variations in boundary heat flow have a secondary influence on the interior convection (Calkins \BOthers., \APACyear2015). While significant effects may occur if the spatial scale of the thermal anomaly is comparable to the vertical scale of convection (e.g., Davies \BOthers., \APACyear2009), it is unclear whether these effects will persist across the entire fluid depth (Davies \BBA Mound, \APACyear2019).

Future numerical work should (1) strive for more realistic Ekman and Rayleigh numbers, (2) tackle the effect of boundary conditions, especially the Stephan-type condition at the top boundary due to melting/freezing of water along the interface and fixed heat flux along the bottom boundary, (3) consider both temperature and salinity buoyancy sources (e.g., Vance \BBA Brown, \APACyear2005; Jansen, \APACyear2016), and (4) couple convectively and mechanically driven flows (e.g., Le Bars \BOthers., \APACyear2015).

Future missions to the outer solar system may be able to better constrain the ocean flows and test the predictions of our calculations and convection models. Looking specifically to the Jovian system, the Europa Clipper and JUICE missions will determine the ocean thickness and salinity and may be able to place constraints on spatial variations of ice shell thickness (e.g., Phillips \BBA Pappalardo, \APACyear2014; Grasset \BOthers., \APACyear2013). Ice penetrating radar will provide information on ice shell thermophysical structure and constrain ice-ocean exchange processes (e.g., Kalousová \BOthers., \APACyear2017), while magnetometer measurements may allow probing of ocean currents through their induction of magnetic fields (e.g., Tyler, \APACyear2011).

Acknowledgements.

I thank Jonathan Aurnou, Baptiste Journaux, and Steve Vance for their helpful comments as well as Gabriel Tobie and Christophe Sotin for their constructive reviews. This work was supported by NASA Grant NNX14AR28G. Computational resources were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. The MagIC code is publicly available at https://magic-sph.github.io/contents.html. All data is provided within the publication pages.