Ocean dynamics of outer solar system satellites

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 and is possible for Titan, and may help drive geologic activity via thermocompositional diapirism in the ice shell. Weak rotational influence allows concentric zonal flows and overturning cells with no preferred orientation. Predictions for Ganymede's ocean span all of these regimes.


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 & Pappalardo, 2016;Lunine, 2017). 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 & Saur, 2017). 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 & Goodman, 2009;Soderlund et al., 2014) or has mechanically driven flows (e.g., Tyler, 2008;Lemasquerier et al., 2017;Wilson & Kerswell, 2018). 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 et al., 2013). 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 & Barr, 2004;Soderlund et al., 2014).
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 & Pappalardo, 2014) and JUICE (Grasset et al., 2013) 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 of these results are discussed in Section 4.

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 et al., 2015). Cheng et al. (2018) combines asymptotic predictions, laboratory experiments, and numerical simulations to characterize the behavior of rotating thermal convection as a function of the dimensionless Ekman, Rayleigh, and Prandtl numbers. The Ekman number, E = ν/2ΩD 2 , represents the ratio of rotational to viscous timescales; thus, low Ekman numbers signify rapid rotation rates in planetary interiors. The Rayleigh number, Ra = αg∆T D 3 /νκ, is the ratio of the thermal diffusion time to the viscous buoyant rise time; large Rayleigh numbers denote strong buoyancy forcing. The Prandtl number, P r = ν/κ, defines the ratio of thermal to viscous diffusion timescales. Here, ν is kinematic viscosity, Ω is rotation rate, D is ocean thickness, α is thermal expansivity, g is gravitational acceleration, ∆T is superadiabatic temperature contrast, and κ is thermal diffusivity. Cheng et al. (2018) 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 GT U , 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 U N R , the influence of rotation is lost globally ("nonrotating heat transfer" regime). Significant debate exists in the community on the scaling laws for Ra GT U and Ra U N R ; thus, I consider the bounding scalings highlighted by Cheng et al. (2018) for each.
Using this regime diagram, one can predict the convective regime of a system if the Ekman, Rayleigh, and Prandtl numbers can be estimated. The Prandtl number depends only on fluid properties and is estimated to be P r ≈ 13 assuming temperatures near the freezing point and terrestrial salinities (Melosh et al., 2004;Nayar et al., 2016). The Ekman number is also relatively easy to calculate since it only requires assumptions about the fluid viscosity, rotation rate, and ocean thickness (see Table 1). I use the internal structure models of Vance et al. (2018) to determine D and consider three possible ocean thickness values to represent the uncertainty at present as well as over the satellites' evolutions. The nominal thickness is taken to be the average value of those predicted by their models (see Tables 5-8 of Vance et al. (2018)). As an upper limit, I consider the maximum combined liquid water-ice layer thicknesses of their models. The lower limit is assumed to be the minimum ocean thickness in their models or half of the nominal thickness, whichever is less. Enceladus has the largest nominal ocean Ekman number of E ∼ 2 × 10 −11 , while the ocean of Ganymede has the lowest at E ∼ 1 × 10 −12 .
The Rayleigh number is more difficult to estimate because it requires knowledge of the superadiabatic temperature contrast ∆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, N u = qD/ρC p κ∆T . Following Soderlund et al. (2014), I leverage N u−Ra scalings to solve for ∆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  Prandtl number, P r = ν/κ 13 13 13 13 Ekman number, E = ν/ΩD 2 2 × 10 −11 2 × 10 −12 3 × 10 −12 1 × 10 −12 Rayleigh number, Ra = αg∆T D 3 /νκ 10 16 − 10 21 10 18 − 10 23 10 19 − 10 22 10 17 − 10 25 theoretical limit of N u = 0.07Ra 1/3 (e.g., Gastine et al., 2015). Conversely, in the rapidly rotating limit, heat transfer is predicted to follow N u = 0.15Ra 3/2 (2E) 2 (Gastine et al., 2016). As a result, the temperature contrast is given by in the non-rotating regime and by in the rapidly-rotating regime. Here, I consider three ocean thicknesses (maximum, mean, minimum) as well as a range of heat flow values given their large uncertainties. As above, the nominal heat flux is set to the average value predicted by the Vance et al. (2018) models; the upper and lower bounds are assumed to be an order of magnitude larger and smaller. At Enceladus, I decrease the lower bound by another factor of six since the observed heat flux (Howett et al., 2011) is approximately sixty times greater than the radiogenic value (Chen et al., 2014). Temperature anomalies are found in the range of tenths to hundreds of milliKelvin. Rayleigh numbers span from Ra ∼ 10 16 for the lower Enceladus limit to Ra ∼ 10 25 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 et al., 2014;Mitri et al., 2014).  Table 1). Our numerical simulations are denoted by stars: red (E = 3.0 × 10 −5 ), magenta (E = 7.5 × 10 −5 ), cyan (E = 1.5 × 10 −4 ), blue (E = 3.0 × 10 −4 ), and black (E = 7.5 × 10 −4 ). The solid black line denotes the scaling for the onset of convection (Chandrasekhar, 1961), the yellow lines bound the range of predicted transitions from the geostrophic turbulence regime to the unbalanced boundary layer regime (Ecke & Niemela, 2014;Julien et al., 2012), and the red lines bound the range of predicted transitions from the unbalanced boundary layer regime to the nonrotating heat transfer regime (Gastine et al., 2016;Gilman, 1977). The Ra RoC=1 U N R and Ra EN 14 GT U scalings depend on the Prandtl number, P r; bold lines assume P r = 13 following estimates for icy satellite oceans while light lines assume P r = 1 as used in the simulations. Figure 1 plots the resulting estimates of the Ekman and Rayleigh numbers on the convective regime diagram defined by Cheng et al. (2018). 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 (c.f. Miquel et al., 2018). 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 GT transition.

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, 2002;Gastine & Wicht, 2012) 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, χ = r i /r o = 0.9. The system is further defined by the Ekman, Rayleigh, and Prandtl numbers. Following Soderlund et al. (2014), the boundaries are impenetrable, stress-free, and isothermal, and fluid motions derived from compositional buoyancy and orbital dynamics are neglected for simplicity.
Five cases that span a convective regime space consistent with the icy satellite ocean predictions ( Fig. 1) are considered. The Rayleigh and Prandtl numbers are fixed to Ra = 3.4 × 10 7 and P r = 1, and the Ekman number is increased from E = 3.0 × 10 −5 to E = 7.5 × 10 −4 . Hyperdiffusivities are not employed (c.f. Zhang & Schubert, 2000). The numerical grids have 73 radial points, 320 latitudinal points, and 640 longitudinal points for cases with E ≥ 7.5 × 10 −5 and 65 radial points, 640 latitudinal points, and 1280 longitudinal points for the E = 3.0 × 10 −5 case. Each case was initiated with a random temperature perturbation or restarted from a lower Ekman number case. Figure 2 shows the mean velocity fields of each model. In the lowest Ekman number case, Coriolis forces organize the flow into narrow radial 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 & Olson, 2001;Christensen, 2001;Heimpel et al., 2005;Gastine et al., 2014). When the Ekman number is increased to E = 7.5 × 10 −5 , the system becomes increasingly turbulent and the convective flows lose some their axial coherence. The zonal flow pattern is relatively unchanged at high latitudes, but one less jet forms at large cylindrical radii, causing the equatorial jet to become retrograde. The mean radial flows become larger scale with a pronounced upwelling near the equator and downwellings at mid-latitudes, essentially forming Hadley-like meridional circulation cells in each hemisphere. The high latitude flows still retain some axialization, however. Further increasing the Ekman number (1.5 × 10 −4 ≤ E ≤ 3.0 × 10 −4 ) reduces the influence of the Coriolis force sufficiently to allow mixing of absolute angular momentum (e.g., Gilman, 1978;Aurnou et al., 2007;Gastine et al., 2013). As a result, zonal flows are retrograde at large cylindrical radii and prograde closer to the rotation axis. The width of the equatorial upwelling also broadens and loses alignment with the rotation axis. Both mean zonal and radial flow speeds increase by a factor of five compared to the E = 7.5 × 10 −5 case. In the highest Ekman number case, 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 become increasingly concentric rather than vertical due to the increased role of viscous transport of angular momentum (Brun & Palacios, 2009).
These velocity field changes lead to distinct changes in heat flow patterns. Figure 3 shows the mean radial temperature gradient along the outer boundary. In the lowest Ekman number case, heat flow peaks at high latitudes with minima near the equator due to the axialization of convective flow structures and development of the strong equatorial jet (e.g., Aurnou et al., 2008). At intermediate Ekman numbers, formation of the Hadley-like circulation cells enhance low latitude heat flux with minima at mid-latitudes. Secondary peaks are found at high latitudes due to turbulent heat transfer associated with vertically ascending plumes. In the highest Ekman number case, the heat flow pattern reflects the development of multiple circulation cells that drift over time. The degree of mixing within the interior is reflected by the ∂T /∂r amplitude, with larger values denoting a more isothermalized interior and thinner thermal boundary layers.

Discussion
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. Enceladus' ocean may then be represented by the E = [3.0×10 −5 , 7.5×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 U N R = E −2 P r (e.g., Gilman, 1977) is assumed, then all of these oceans are near the center of UBL regime such that the E = [3.0 × 10 −5 , 7.5 × 10 −5 ] models would again be most appropriate. If the transition instead follows Ra Ga16 U N R = 100(2E) −12/7 (Gastine et al., 2016), then the E = [1.5 × 10 −4 , 7.5 × 10 −4 ] models would be most pertinent, with nominal parameter estimates best represented by the E = 1.5 × 10 −4 model.
Below, I discuss the implications for each satellite and look ahead to future missions. This analysis focuses on the heat flux pattern along the ice-ocean interface (i.e. ∂T /∂r| r=ro ). 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 & Perkin, 1986) 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 et al., 1994;Eicken et al., 1984), 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 & Barr, 2004;Soderlund et al., 2014). Alternatively, if the ice pump mechanism is not efficient, the ice shell may be more unstable to convection where it is relatively thick (Goodman, 2014).
For Enceladus, I predict the zonal flows to be characterized by multiple jets that alternate in direction ( Fig. 2A a,b). Converting model velocities to dimensional units U = ΩDRo, I expect peak zonal speeds of approximately 0.5 m/s. Meridional circulations are predicted to either be strongly aligned with the rotation axis with speeds up to a cm/s (Fig. 2B a) or to be concentrated in a low latitude upwelling with speeds up to a few cm/s (Fig. 2B b). Instantaneous velocities can locally be an order of magnitude larger for all models. As a result, heat flow along the ice-ocean interface has distinct peaks at either the poles (Fig. 3, red line) or at the equator and the poles secondarily (Fig. 3, magenta  line).
Measurements of Enceladus' 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 (Čadek et al., 2016;Beuthe et al., 2016). Inverting these measurements to infer the ocean heat flux along the ice-ocean interface, Čadek et al. (2019) find peak flux near the poles with a minima at the equator. These observations are consistent with the heat flux distribution of the E = 3.0 × 10 −5 model (Fig. 3), which is the appropriate analog if a low internal heat flux is assumed (Fig. 1).
Given the similarities in regime predictions for Titan, Europa, and Ganymede, they are considered together here. Assuming the nominal model and the Ra Ga16 U N R scaling, these satellites are predicted to have three zonal jets with retrograde flow outside the tangent cylinder (imaginary right cylinder aligned with the rotation axis and tangent to the seafloor) and prograde flow at higher latitudes ( Fig. 2A c). Jet speeds reach ∼ [2, 4, 5] m/s for [Titan, Europa, Ganymede], respectively. The local convective motions are fast (order m/s) and not strongly constrained by rotation; on average, they organize to form an equatorial upwelling with peak speeds of ∼ [5, 9, 12] cm/s (Fig. 2B c). Heat flow along the ice-ocean interface is maximum near the equator, with secondary maxima near the poles (Fig. 3, cyan line). Alternatively, for the Ra RoC=1 U N R = E −2 P r scaling, these satellites are predicted to behave similarly to Enceladus as discussed above, except with respect to dimensionalized flow speeds that will be a factor of ∼ [0.8, 1.4, 1.7] times faster.
Looking to Titan, the satellite's surface topography shows polar depressions compared to relatively elevated low latitudes (Nimmo & Bills, 2010;Choukroun & Sotin, 2012;Kvorka et al., 2018). The three intermediate Titan-relevant models have peak heat fluxes near the equator (Fig. 3, magenta, cyan, and blue lines). These peaks are significantly stronger than the secondary enhancements at high latitudes such that melting would cause the equator to be relatively thin compared to the poles, opposite to the observations. If this region were infilled with relatively pure marine ice, however, an equatorial bulge may result through Pratt isostasy. Associated thermocompositional diapirism could potentially lead to cryovolcanism that appears to be concentrated at low latitudes (e.g., Sohl et al., 2014). Alternatively, the lowest Ekman number model may be appropriate if we assume the Ra RoC=1 U N R = E −2 P r scaling and/or if a stable salinity gradient reduces the effective buoyancy forcing of the ocean (Ra). High oceanic heat flux and associated melting at high latitudes would then be consistent with the polar depression through Airy isostasy (Kvorka et al., 2018).
Europa has a young, tortured surface with geologic features indicating recent activity and the potential for ocean-derived materials (e.g., Figueredo & Greeley, 2003;Fischer et al., 2015). These chaos terrains appear to be located preferentially at low latitudes with a secondary prevalence near the poles (Leonard et al., 2018), and formation models suggest that they may be associated with upwelling diapirs (e.g., Collins & Nimmo, 2009;Schmidt et al., 2011) and marine ice accretion (Soderlund et al., 2014). No large gradients in ice shell thickness have been detected (Nimmo et al., 2007), suggesting an efficient ice pump (c.f. Nimmo, 2004). Given that most Europa-relevant models predict high oceanic heat flux at low latitudes with relatively low flux at mid-latitudes (Fig. 3, magenta, cyan, and blue lines), our new calculations continue to support the thermocompositional diapirism hypothesis.
The 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, 1980) and detection of hydrated salts suggests a subsurface briny layer of fluid (McCord et al., 2001), but no clear patterns are present. Mass anomalies were measured in the northern hemisphere during a single Galileo flyby (Palguta et al., 2006), 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.
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 & Pappalardo, 2014;Grasset et al., 2013). Ice penetrating radar will provide information on ice shell thermophysical structure and constrain ice-ocean exchange processes (e.g., Kalousová et al., 2017), while magnetometer measurements may allow probing of ocean currents through their induction of magnetic fields (e.g., Tyler, 2011).