Energetic Constraints on Ocean Circulations of Icy Ocean Worlds

The upward buoyancy flux is related to the upward flux of heat or compositional buoyancy. In this manuscript we assume a water ocean with dissolved salts, such that compositional buoyancy effects are encapsulated by the salinity. In general we thus allow buoyancy to be some function of potential temperature, salinity, and geopotential height (e.g., Young 2010):

$$\begin{eqnarray}&&b=\tilde{b}({\rm{\Theta }},S,z),\end{eqnarray} \tag{ 5 }$$

where Θ is the potential temperature, ⁷ S is the salinity, and z is the depth relative to some reference geopotential height level. The temperature and salinity evolve according to

$$\begin{eqnarray}&&\displaystyle \frac{D{\rm{\Theta }}}{{Dt}}={\kappa }_{T}{{\rm{\nabla }}}^{2}{\rm{\Theta }},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{DS}}{{Dt}}={\kappa }_{S}{{\rm{\nabla }}}^{2}S,\end{eqnarray} \tag{ 7 }$$

with the boundary conditions

$$\begin{eqnarray}&&-{\kappa }_{T}{\partial }_{z}{\rm{\Theta }}({z}_{\mathrm{bot}})=\displaystyle \frac{1}{\rho {c}_{p}}{{ \mathcal Q }}_{\mathrm{bot}},-{\kappa }_{T}{\partial }_{z}{\rm{\Theta }}({z}_{\mathrm{top}})=\displaystyle \frac{1}{\rho {c}_{p}}{{ \mathcal Q }}_{\mathrm{top}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&-{\kappa }_{S}{\partial }_{z}S({z}_{\mathrm{bot}})={{ \mathcal S }}_{\mathrm{bot}},-{\kappa }_{S}{\partial }_{z}S({z}_{\mathrm{top}})={{ \mathcal S }}_{\mathrm{top}},\end{eqnarray} \tag{ 9 }$$

where κ_T and κ_S are the molecular diffusivities of heat and salt, respectively, and ${{ \mathcal Q }}_{\mathrm{bot}/\mathrm{top}}$ and ${{ \mathcal S }}_{\mathrm{bot}/\mathrm{top}}$ denote the heat and salt fluxes through the bottom and top boundaries, respectively. For simplicity we do not include significant heat sources in the interior, although such sources could be added.

If the equation of state is nonlinear, the conservation equations for Θ and S cannot be readily translated into a conservation equation for buoyancy. The most general way to account for nonlinearities in the equation of state, which is discussed in the Appendix, is to introduce the dynamic enthalpy, which is related to the potential energy. A significantly simpler (and more intuitive) result can be obtained if we assume that horizontal variations in the thermal and haline expansion coefficients are small, such that we can approximate

$$\begin{eqnarray}&&b\approx \langle b\rangle +\langle \alpha \rangle ({\rm{\Theta }}-\langle {\rm{\Theta }}\rangle )+\langle \beta \rangle (S-\langle S\rangle ),\end{eqnarray} \tag{ 10 }$$

where α ≡ g⁻¹∂ _Θb and β ≡ g⁻¹∂_S b are the thermal and haline expansion coefficients, respectively, and 〈·〉 denotes a horizontal average at any given depth. Notice that α is typically positive (except for relatively fresh water at cold temperatures and modest pressures), while β, as defined here, is negative. We can then directly relate the total globally integrated upward advective buoyancy flux (which represents the source of kinetic energy in a buoyancy-driven flow) to the upward advective heat and salt fluxes:

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal B } & \equiv & \displaystyle \int {wbdV}\approx \displaystyle \int \left[g\langle \alpha \rangle \iint w{\rm{\Theta }}{dA}\right]{dz}\\ & & +\,\displaystyle \int \left[g\langle \beta \rangle \iint {wSdA}\right]{dz}\equiv {{ \mathcal B }}^{{\rm{\Theta }}}+{{ \mathcal B }}^{S},\end{array}\end{eqnarray} \tag{ 11 }$$

where we used that $\iint$ wdA = 0 due to volume conservation, and we defined ${{ \mathcal B }}^{{\rm{\Theta }}}$ and ${{ \mathcal B }}^{S}$ as the globally integrated net upward buoyancy fluxes associated with heat and salt advection, respectively.

The total upward heat and salt fluxes can be related to the sources and sinks of heat and salt via their conservation equations (see also sketches in Figures 1 and 2). In particular, the area-integrated vertical fluxes of heat and salt are constrained by the conservation of heat and salt above or below a respective level. In a statistically steady state, where the mean rate of change in Θ and S at any given depth is approximately zero, Equation (6) can be integrated over the volume above or below some depth z to obtain an equation for the vertical buoyancy flux at this depth:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int (w{\rm{\Theta }}-{\kappa }_{T}{\partial }_{z}{\rm{\Theta }}){dA}\\ =\,\displaystyle \frac{1}{\rho {c}_{p}}\left[\displaystyle \int {{ \mathcal Q }}_{\mathrm{bot}}{ \mathcal H }(z-{z}_{\mathrm{bot}}){dA}-\displaystyle \int {{ \mathcal Q }}_{\mathrm{top}}{ \mathcal H }(z-{z}_{\mathrm{top}})\ {dA}\right]\\ =\,\displaystyle \frac{1}{\rho {c}_{p}}\left[\displaystyle \int {{ \mathcal Q }}_{\mathrm{top}}{ \mathcal H }({z}_{\mathrm{top}}-z){dA}-\displaystyle \int {{ \mathcal Q }}_{\mathrm{bot}}{ \mathcal H }({z}_{\mathrm{bot}}-z)\ {dA}\right]\\ \equiv \,\displaystyle \frac{1}{\rho {c}_{p}}{{ \mathcal F }}^{{ \mathcal Q }},\end{array}\end{eqnarray} \tag{ 12 }$$

where ${ \mathcal H }$ is the Heaviside function (${ \mathcal H }(x)=0\ {\rm{for}}\,x\lt 0$ and ${ \mathcal H }(x)=1\ {\rm{for}}\,x\gt 0$), and we defined ${{ \mathcal F }}^{{ \mathcal Q }}$ as the the vertical heat flux across any depth z needed to balance the net heat and salt fluxes through the boundaries above and below the respective level. Similarly, Equation ((7)) can be integrated to obtain

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int ({wS}-{\kappa }_{S}{\partial }_{z}S){dA}\\ =\,\displaystyle \int {{ \mathcal S }}_{\mathrm{bot}}{ \mathcal H }(z-{z}_{\mathrm{bot}}){dA}-\displaystyle \int {{ \mathcal S }}_{\mathrm{top}}{ \mathcal H }(z-{z}_{\mathrm{top}})\ {dA}\\ =\,\displaystyle \int {{ \mathcal S }}_{\mathrm{top}}{ \mathcal H }({z}_{\mathrm{top}}-z){dA}-\displaystyle \int {{ \mathcal S }}_{\mathrm{bot}}{ \mathcal H }({z}_{\mathrm{bot}}-z)\ {dA}\\ \equiv \,{{ \mathcal F }}^{{ \mathcal S }}.\end{array}\end{eqnarray} \tag{ 13 }$$

Notice that interior sources of heat or salt could be included by modifying the definition of ${{ \mathcal F }}^{{ \mathcal Q }}$ and ${{ \mathcal F }}^{{ \mathcal S }}$, which generally need to balance any net sources or sinks above or below the respective level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These definitions allow us to express the globally integrated KE generation associated with heat and salt flux forcing as:

$$\begin{eqnarray}&&{{ \mathcal B }}^{{\rm{\Theta }}}\approx \mathop{\underbrace{\int g\langle \alpha \rangle \displaystyle \frac{{{ \mathcal F }}^{{ \mathcal Q }}}{\rho {c}_{p}}{dz}}}\limits_{{{ \mathcal B }}_{{ \mathcal F }}^{{\rm{\Theta }}}}+\mathop{\underbrace{\int g\alpha {\kappa }_{T}{\partial }_{z}\theta {dV}}}\limits_{{{ \mathcal B }}_{\mathrm{diff}}^{{\rm{\Theta }}}},\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{{ \mathcal B }}^{S}=\mathop{\underbrace{\int g\langle \beta \rangle {{ \mathcal F }}^{{ \mathcal S }}{dz}}}\limits_{{{ \mathcal B }}_{{ \mathcal F }}^{S}}+\mathop{\underbrace{\int g\beta {\kappa }_{S}{\partial }_{z}{SdV}}}\limits_{{{ \mathcal B }}_{\mathrm{diff}}^{S}}.\end{eqnarray} \tag{ 15 }$$

The first term on the rhs of Equations (14) and (15) can be interpreted as the potential energy source associated with heat and salt forcing at the lower and upper boundaries, while the second term is the potential energy source or sink due to diffusion (which provides a source of mechanical energy whenever the stratification is statically stable, and a sink when it is unstable).

If vertical variations in g α and g β are also small, thermal forcing at the boundaries provides a source of energy that can drive a circulation if, and only if,

$$\begin{eqnarray}&&\alpha \int {{ \mathcal F }}^{Q}{dz}\gt 0,\end{eqnarray} \tag{ 16 }$$

i.e., for α > 0, heating, on average, needs to occur at a greater depth than cooling, while for α < 0 cooling needs to occur at a greater depth than heating. Similarly, salt fluxes can drive circulation only if

$$\begin{eqnarray}&&\beta \int {{ \mathcal F }}^{S}{dV}\gt 0.\end{eqnarray} \tag{ 17 }$$

Since generally β < 0, this requires that salt needs to be removed at greater depth and added at shallower depth. Large vertical variations in g α and g β can be accounted for using the generalized relation in Equations (14) and (15), which allow us to compute the mechanical energy input associated with any given heat and salt flux boundary conditions more generally. When horizontal variations in the thermal and haline expansion coefficient are nonnegligible, the amount of energy that can be converted to kinetic energy depends on the specifics of the circulation, as elaborated on in the Appendix.

Equations (14) and (15) also highlight the potential role of vertical diffusion. Indeed, buoyancy gain and loss at the same level combined with vertical diffusion can drive circulation, which is sometimes referred to as horizontal convection (e.g., Sandström 1908; Paparella & Young 2002; Wunsch & Ferrari 2004; Hughes & Griffiths 2008). As the energy source in such a circulation is derived from diffusion, the energy dissipation has to go to zero in the limit of vanishing diffusivity. As a result, Paparella & Young (2002) argue that no "turbulent" circulation can be maintained in the limit of vanishing diffusivity (and viscosity), where "turbulence" is defined to imply a forward energy cascade with a dissipation rate that becomes independent of the viscosity. ⁸ In reality, molecular diffusion is not vanishingly small, and we will discuss its potential role below. It is also worth noting at this point that many numerical simulations employ strongly enhanced "eddy diffusivities," which can act as a substantial energy source in numerical simulations. However, it is important to remember that these parameterizations are meant to represent the effect of turbulent advection—that is, they parameterize the unresolved contribution to the wb term in Equation (3). The apparent energy source associated with parameterized "eddy diffusion" therefore needs to be interpreted as a conversion of unresolved turbulent kinetic energy to potential energy, and accordingly can be justified only if there exists a corresponding source of unresolved turbulent kinetic energy.

In equilibrium, the total sources and sinks of mechanical energy have to be in balance. The energy sources thereby provide a constraint on the energy dissipation, which in turn provides some constraint on the flow. Relating the energy dissipation rate to the kinetic energy of the flow itself is not straightforward, as it depends on the characteristics of the flow field, but we can make some progress by considering specific flow regimes and estimating the range of parameters and scales over which the flow regimes are expected to hold. In the following, we first discuss the relationship between the kinetic energy and the dissipation rate for a turbulent flow that is largely unaffected by rotation, as well as the conditions under which the assumption that rotation is negligible breaks down. We then derive an alternative relationship between the kinetic energy and the dissipation rate in the opposite limit where rotation is a leading-order effect and the most energetic flows are geostrophically balanced (i.e., the leading-order momentum balance is between the pressure gradient and Coriolis acceleration).

2.2.1. Isotropic Turbulence and the Role of Rotation

We first assume a fully turbulent flow unconstrained by the influence of rotation in the interior (i.e., away from the boundaries), and we assume that, in this limit, the dissipation rate is independent of the value of molecular viscosity. In particular, we start by considering Kolmogorov's theory for the kinetic energy spectrum in the inertial cascade of isotropic turbulence, which suggests that

$$\begin{eqnarray}&&E(k)={ \mathcal K }{\epsilon }^{2/3}{k}^{-5/3},\end{eqnarray} \tag{ 18 }$$

where E(k) is the kinetic energy spectrum as a function of the wavenumber, k, is the turbulent spectral KE flux, which in turn is equal to the dissipation rate, and ${ \mathcal K }\approx 1.5$ is the Kolmogorov constant (e.g., Vallis 2006). Integrating over the energy inertial range yields

$$\begin{eqnarray}&&{E}_{t}\approx 2{k}_{0}^{-2/3}{\epsilon }^{2/3},\end{eqnarray} \tag{ 19 }$$

where k₀ is the "injection scale" (where the KE spectrum flattens out), and E_t is the total turbulent kinetic energy in the inertial range.

We can use Equation (19) to estimate the kinetic energy of turbulent flows up to the largest scales of an isotropic turbulent energy inertial range. Isotropy may be broken by the geometry (e.g., the vicinity to a boundary), or by the effects of stratification or rotation. We will here consider in particular the importance of rotation, which can break isotropy and fundamentally change the nature of the turbulent flow throughout the water column.

The effect of rotation on the turbulent flow in the interior can be characterized by the flow Rossby number, which we here define as the ratio of the inverse "eddy turnover timescale," ${\tau }_{\mathrm{eddy}}^{-1}\sim \sqrt{{E}_{t}}{k}_{0}$, to the rotation rate, Ω. This ratio measures the relative magnitude of the nonlinear advection term to the Coriolis term in the momentum equation, and thus provides a useful and widely used measure for the role of rotation in the dynamics of high-Reynolds-number flows (e.g., Vallis 2006). Using Equation (19), and L ≡ 2π/k₀, the turbulent flow Rossby number can be estimated as

$$\begin{eqnarray}&&{\mathrm{Ro}}_{t}\equiv \displaystyle \frac{\sqrt{{E}_{t}}}{{\rm{\Omega }}L}\approx \displaystyle \frac{{\epsilon }^{1/3}}{{\rm{\Omega }}{L}^{2/3}}.\end{eqnarray} \tag{ 20 }$$

Equation (20) suggests a maximum length scale for turbulent flows unaffected by rotation (see also Fernando et al. 1991; Jones & Marshall 1993; Maxworthy & Narimousa 1994; Bire et al. 2022):

$$\begin{eqnarray}&&{L}_{\mathrm{rot}}\approx \displaystyle \frac{{\epsilon }^{1/2}}{{{\rm{\Omega }}}^{3/2}}.\end{eqnarray} \tag{ 21 }$$

Maxworthy & Narimousa (1994) and others found that, once rotation becomes important, the convective Rossby number of a rotating plume in the interior scales as

$$\begin{eqnarray}&&{\mathrm{Ro}}_{\mathrm{rc}}=\displaystyle \frac{{\epsilon }^{1/2}}{{{\rm{\Omega }}}^{3/2}H},\end{eqnarray} \tag{ 22 }$$

where H is the depth of the convecting fluid. Equation (22) is likely to be a better predictor of the convective flow Rossby number when ocean depth convection is strongly affected by rotation (i.e., L_rot < H). If we set L = H, we find ${\mathrm{Ro}}_{\mathrm{rc}}={\mathrm{Ro}}_{t}^{3/2}$, and hence both definitions give the same prediction for the critical length/depth scale at which Ro ≈ 1 and rotation becomes important, ⁹ which indeed also follows directly from dimensional analysis if we postulate that this scale shall depend only on and Ω. If this length scale is significantly smaller than the scale of the most energetic flows, we expect these flows to be strongly affected by rotation. For convective flows unconstrained by rotation, the natural length scale for the largest convective motion is the depth of the ocean, such that we expect the interior convective dynamics to become strongly affected by rotation if L_rot is much smaller than the depth of the ocean (see also Fernando et al. 1991).

2.2.2. Boundary-layer Dissipation in Geostrophic Dynamics

When rotation becomes of dominant importance it is likely that much of the energy becomes trapped in geostrophically balanced vortices and large-scale mean flows, which result from the upscale kinetic energy cascade associated with quasi-balanced turbulent motions (e.g., Vallis 2006). The lack of a forward energy cascade means that dissipation is likely to be limited mostly to turbulent boundary layers near the seafloor and the ice–ocean interface.

We can estimate the energy dissipation per unit area in a turbulent boundary layer as (e.g., Jansen 2016, and references therein):

$$\begin{eqnarray}&&\int {\epsilon }_{\mathrm{BL}}{dz}={c}_{D}| {U}_{g}{| }^{3},\end{eqnarray} \tag{ 23 }$$

where the integral on the lhs is over the depth of the turbulent boundary layer, c_D is the turbulent drag coefficient, and U_g is the characteristic near-surface geostrophic velocity. The value of c_D depends on the surface roughness, with c_D ≈ 0.0025 an empirical average value that is commonly used for the drag coefficient for Earth's seafloor (e.g., Egbert et al. 2004; Sen et al. 2008). The skin drag coefficient under smooth ice can be estimated to be around c_D ≈ 0.002, although rough morphology can significantly increase this value (e.g., Brenner et al. 2021). Lacking information about the ice roughness, and noting the order-of-magnitude nature of our estimates, we will here use c_D ≈ 0.001–0.01 at both the seafloor and under the ice sheet at the top. Experience from oceanic and atmospheric modeling has shown that drag coefficients of this order produce reasonable results for a wide range of flows (e.g., Smagorinsky et al. 1965; Egbert et al. 2004; Sen et al. 2008; Chen et al. 2018; Adcroft et al. 2019), which instills confidence that a similar coefficient can be used to obtain useful order-of-magnitude estimates for icy-moon oceans.

Averaging the energy dissipation in the top and bottom boundary layers over the depth of the whole water column gives an average energy dissipation rate per unit volume

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{BL}}\approx \displaystyle \frac{{2}^{5/2}{c}_{D}}{H}{E}_{g}^{3/2},\end{eqnarray} \tag{ 24 }$$

where H is the depth of the ocean, and E_g = 1/2∣U_g ∣² is the characteristic KE of the geostrophic flow. We can solve Equation (24) to get a crude estimate for the characteristic KE in a flow that is dominated by large-scale balanced dynamics, where KE dissipation occurs primarily in turbulent boundary layers:

$$\begin{eqnarray}&&{E}_{g}\approx \displaystyle \frac{{H}^{2/3}}{{2}^{5/3}{c}_{D}^{2/3}}{{\epsilon }_{\mathrm{BL}}}^{2/3}.\end{eqnarray} \tag{ 25 }$$

In general, balanced dynamics may coexist with intermediate-Rossby-number and high-Rossby-number convective turbulence and/or tidal waves, in which case interior energy dissipation may be important and thus _BL = − _int < , where _int denotes energy dissipation away from the boundary layers. It is also important to note that numerical simulations of planetary circulation typically use very high artificial viscosities for numerical stability, which can lead to a large, albeit probably unphysical, dissipation of balanced KE in the interior (e.g., Jansen 2016).

Given the energy dissipation rate, , and the rotation rate of the planetary body, Ω, we can use Equation (21) to estimate a critical scale L_rot above which we may expect turbulent motions to become strongly affected by rotation (Figure 4). Using the energy dissipation rates estimated above and Ω ≈ 8 × 10⁻⁵ s⁻¹ for Snowball Earth, 5 × 10⁻⁵ s⁻¹ for Enceladus, and 2 × 10⁻⁵ s⁻¹ for Europa, we find the scale at which convection is expected to become strongly affected by rotation from Equation (21) to be L_rot ∼ O(10 m) for Snowball Earth and Europa and L_rot ≲ 1 m for Enceladus. In addition to parametric uncertainties (quantified in Figure 4 and Table 2), it is possible that assumptions leading to Equation (21) may not hold. In particular, we here assumed that the interior dissipation (which enters the scaling for L_rot) is of similar order as the total kinetic energy dissipation. If the dissipation is dominated by boundary layers, the length scale where the interior flow would become strongly affected by rotation would be further reduced. Regardless, L_rot is much smaller than the ocean depth in all three cases, and we therefore expect deep convection and any other potential large-scale flow to be strongly affected by rotation. It is also worth noting that, for Snowball Earth and Europa, L_rot is likely to be much larger than the Kolmogorov scale where viscous dissipation occurs, ${L}_{\nu }={({\nu }^{3}/\epsilon )}^{1/4}$, which is expected to be on the order of a few centimeters, thus indicating that an isotropic inertial range is expected to exist at scales smaller than L_rot. On Enceladus it is not clear whether L_rot is significantly larger than the Kolmogorov scale, raising the possibility that no isotropic cascade exists. In that case, the relationship in Equation (18) may not be expected to hold at any scale, but the conclusion that the large-scale flow would be strongly affect by rotation would remain true, as indeed the flow at all scales (down to the viscous dissipation scale) would be strongly affected by rotation. We also note that, for all icy oceans considered here, L_rot is significantly smaller than the smallest length scales that can be resolved in global-scale numerical simulations of these oceans, suggesting that large-eddy-simulations that can explicitly resolve part of the isotropic turbulent inertial range are not feasible with realistic parameters (see Bire et al. 2022).

Zoom In Zoom Out Reset image size

3.1.1. Comparison to Existing Simulation Results

The importance of rotation in the oceans of Snowball Earth, Europa, and Enceladus is qualitatively consistent with the Snowball Earth simulations of Ashkenazy et al. (2013) and Jansen (2016), the Europa simulations of Ashkenazy & Tziperman (2021), and the Enceladus simulations of Kang et al. (2020, 2022) and Zeng & Jansen (2021). For illustration, Figures 5 and 6 show simulation results from Jansen (2016) and Zeng & Jansen (2021) for heating-driven flows in the oceans of Snowball Earth and Enceladus, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The flow field in the Snowball Earth simulation shows a horizontal eddy field characteristic of geostrophic turbulence (Figure 5(b)), with the flow largely barotropized—i.e., varying little in the vertical (Figure 5(b)). The characteristic vertical-flow Rossby numbers are small (Figure 5(a)), qualitatively consistent with the predictions in Section 2.2.1. However, we notice that the vertical flows here are associated with the quasi-geostrophic eddies rather than convection, as geostrophic eddies establish a statically stable stratification. Neither of the convective Rossby-number scalings in Equations (20) and (22) are therefore expected to apply to this flow (although the conclusion that rotation is important still holds, as it is indeed a necessary condition for the development of quasi-geostrophic eddies in the first place).

The vertical-flow Rossby numbers in the high-salinity Enceladus ocean simulations of Zeng & Jansen (2021) are even smaller and show a pattern of grid-scale convection (Figure 6(a)). The grid-scale convection is qualitatively consistent with the prediction that the scale at which rotation affects convective plumes (O(1m)) is much smaller than the model resolution (∼1 km), but unfortunately also implies that the details of the simulated flow are expected to be sensitive to the model resolution and somewhat arbitrary modeling choices, such as parameterized "eddy" viscosities and diffusivities (see Zeng & Jansen 2021). Nevertheless, the importance of rotation is robust and can also be seen in the organization of the zonal mean flow into jets that are approximately aligned with the axes of rotation (Figure 6(c)).

The result that flows are very strongly affected by rotation may appear at odds with simulation results for Europa's ocean by Soderlund et al. (2014), which show strong turbulent convection with relatively moderate flow Rossby numbers. The apparent contradiction can readily be understood by noting that the dimensional vertical heat flux in the DNS of Soderlund et al. (2014) would amount to about 10⁵ Wm⁻² after redimensionalizing with Europa's planetary parameters—versus a heat flux of 0.02–0.1 Wm⁻² assumed here. Increasing the vertical heat flux on Europa to 10⁵ Wm⁻² would yield L_rot ≈ 25 km, which in turn would make significant turbulent convection that is only relatively weakly affected by rotation plausible. Indeed, using Equation (20) with L ≈ H ≈ 100 km (for the depth of Europa's ocean) and a heat flux of 10⁵ Wm⁻² (which gives ≈ 5 × 10⁻⁶ m² s⁻³) we can estimate a Rossby number for full-depth turbulent convection of Ro_t ≈ 0.4, which is broadly consistent with the Rossby numbers of convective motions found in the simulations of Soderlund et al. (2014). The large implied dimensional heat flux in the DNS of Soderlund et al. (2014) is meant to account for the large viscosities and diffusivities that are necessary in numerical simulations. However, the energetic argument discussed here clarifies that the greatly amplified heat flux is expected to lead to much larger dimensional flow speeds. Hence, the magnitude of the dimensionalized velocities from the DNS of Soderlund et al. (2014), and the associated large-scale-flow Rossby numbers, should not be taken at face value.

In summary, our results suggest that the turbulent flow characteristics in the ocean's interior are strongly constrained by the rate of energy input, which in the case of thermally driven flows is given by the vertical heat flux, and the nature of the energy dissipation. The dissipation mechanisms and in some cases the energy source are likely to be misrepresented in numerical modeling studies, which cannot resolve all relevant scales of motion. The quantitative results of icy-moon ocean simulations therefore need to be interpreted with care.

At scales much larger than L_rot, we expect the motion to be strongly affected by rotation, leading to quasi-balanced flows that do not undergo a forward energy cascade. Energy dissipation then may be limited mostly to turbulent boundary layers, such that we can assume _BL ∼ . In this case, Equation (25) may provide a useful estimate for the total KE and thus typical flow speeds as a function of the energy dissipation rate, , and the depth of the ocean, H (Figure 7). Assuming an ocean depth of H ≈ 2 − 3 km for Snowball Earth (Ashkenazy et al. 2013; Yang et al. 2017), H ≈ 50 − 150 km for Europa (Anderson et al. 1998; Vance et al. 2018), and H ≈ 10–50 km for Enceladus (Čadek et al. 2016; Hemingway et al. 2018; Vance et al. 2018), we find

$$\begin{eqnarray}&&{E}_{g}\approx {10}^{-4}\mbox{--}{10}^{-3}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-2}\quad {\rm{for}}\,{\rm{Snowball}}\,{\rm{Earth}},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{E}_{g}\approx 5\times {10}^{-5}\mbox{--}3\times {10}^{-3}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-2}\quad {\rm{for}}\,{\rm{Europa}},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{E}_{g}\lesssim {10}^{-4}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-2}\quad {\rm{for}}\,{\rm{Enceladus.}}\end{eqnarray} \tag{ 32 }$$

These results suggest potential balanced flow velocities of up to a few cm s⁻¹ for Snowball Earth and Europa and up to one cm s⁻¹ for Enceladus. These estimates should generally be viewed as upper bounds as we here assumed that all energy input goes into balanced motions that are dissipated only in the boundary layers, while we ignored any potential additional dissipation in the interior, which would reduce the expected kinetic energy level.

Zoom In Zoom Out Reset image size

3.2.1. Comparison to Existing Simulation Results

Flow velocities on the order of centimeters per second are broadly consistent with the Snowball Earth simulations of Jansen (2016; Figure 5(b)), as well as the Snowball Earth simulations of Ashkenazy et al. (2013) and Europa simulations of Ashkenazy & Tziperman (2021). Figure (8)) shows the energy budget terms for the Snowball Earth simulation of Jansen (2016; see Figure 5(b) for the corresponding flow fields). As assumed in our scaling argument, the primary energy source is associated with the vertical heat flux imposed by the seafloor heating, whose magnitude ${{ \mathcal B }}_{{ \mathcal F }}^{{\rm{\Theta }}}/V\sim \epsilon \approx 2.5\times {10}^{-11}$ m² s⁻³ is consistent with our estimate in Equation (28). Somewhat more than half of this kinetic energy is dissipated at the boundary (where the model uses a quadratic drag, consistent with the assumption in Equation (23)), with the remainder dissipated by viscous dissipation in the interior. It is likely that the interior dissipation is unrealistically large in the simulations due to the limited resolution and required large "eddy" viscosity. However, the interior dissipation does not affect the order of magnitude of the velocity as long as boundary friction is a first-order contributor to dissipation, and hence the simulated velocities are consistent with the prediction in Equation (30).

Zoom In Zoom Out Reset image size

The high-salinity Enceladus ocean simulation of Zeng & Jansen (2021) has flow velocities on the order of 10⁻⁴ m s⁻¹ (Figure 6), which is broadly similar to results reported in Kang et al. (2020) but 2 orders of magnitude weaker than our estimate for the maximum flow speeds (Equation (32)). The parameters used in the Enceladus simulation of Zeng & Jansen (2021) are expected to place the energy input about a factor of 10 below our upper bound estimate, which by itself would translate to flow speeds that are only about a factor of 2 below the predicted maximum (see Equation (23)). Some discrepancies may be expected due to the parameterized boundary-layer drag. Instead of the quadratic boundary-layer drag assumed in Equation (24), the simulations of Zeng & Jansen (2021; and also Kang et al. 2020) apply a relatively strong linear drag parameterization. However, the stronger boundary-layer drag (which can readily be incorporated into the scaling law) can only explain a relatively small part of the misfit. Instead, interior viscous dissipation is key to explaining the much weaker velocities in the simulations, as shown explicitly in Figure 8(b) for the Enceladus simulation of Zeng & Jansen (2021). Strong dissipation of kinetic energy associated with the (numerically necessary) large interior viscosity implies that the total dissipation rate is much larger than the boundary-layer dissipation (i.e., ≫ _BL), thus leading to much weaker velocities. The large viscosities in the simulations (which cannot resolve the inertial cascade range) cannot be justified on physical grounds. Whether a dominant contribution from interior dissipation is likely to exist on Enceladus, or is purely an artifact of insufficient resolution and artificially large viscosity in the simulations, remains unknown and depends on how efficiently energy is transferred into balanced motions. In either case, the estimate in Equation (32) is expected to remain valid as an upper bound, although the true velocities may be much smaller.

The Europa simulations of Soderlund et al. (2014) as well as the Enceladus-like small Ekman number simulations of Soderlund (2019) suggest much larger velocities than estimated here (if redimensionalized with the actual ocean depth and rotation rate of the respective icy moons). Qualitatively, this discrepancy may be expected as a result of the much larger dimensional vertical heat flux, which leads to a much larger energy source, as discussed above. However, the scaling argument leading to Equation (25) is not expected to apply to the simulations of Soderlund et al. (2014) and Soderlund (2019), due to different boundary conditions, thus not allowing for a quantitative prediction of their simulation results. In particular, Soderlund et al. (2014) and Soderlund (2019) use smooth boundaries with free-slip boundary conditions, which implies no boundary drag.

If salinity and pressure are relatively low, the thermal expansion coefficient near the freezing point is negative. This scenario was first suggested for Europa by Melosh et al. (2004), although due to the modest pressures that are required for a negative thermal expansion coefficient, it appears most likely to be relevant for Enceladus, and was considered in the low-salinity Enceladus ocean simulations of Zeng & Jansen (2021) and Kang et al. (2022). To estimate the energetics of such an ocean, we assume a stably stratified layer below the ice sheet and above some depth z_strat, as found in the low-salinity simulation of Zeng & Jansen (2021; Figure 9). In the stratified layer, the temperature decreases upward from the "critical temperature," Θ_c , where α(Θ_c , z_strat) = 0 to the freezing point, Θ_f , at the bottom of the ice sheet, z_ice (where α(Θ_f , z_ice) < 0). For illustrative purposes, we assume that the seafloor and ice sheet are flat (i.e., both boundaries follow a geopotential height surface), such that ${{ \mathcal F }}^{{ \mathcal Q }}$ is vertically constant and equal to the total bottom heat flux. Assuming again that horizontal variations in the thermal expansion coefficient are small, Equation (14) gives the KE generation as

$$\begin{eqnarray}&&{{ \mathcal B }}^{{\rm{\Theta }}}=\int g\langle \alpha \rangle \left[\displaystyle \frac{{{ \mathcal F }}^{{ \mathcal Q }}}{\rho {c}_{p}}+A{\kappa }_{T}{\partial }_{z}\langle \theta \rangle \right]{dz}.\end{eqnarray} \tag{ 33 }$$

In the stratified layer, 〈α〉 < 0, such that the first term in the integral in Equation (33) amounts to an energy sink. One plausible solution is then that the vertical heat flux through the stratified layer is balanced by molecular diffusion (i.e., the second term in the integral in Equation (33)), which amounts to a conversion from internal to potential energy. In this case the stratified layer would only be a few tens of meters thick (Zeng & Jansen 2021), and most of the ocean would be convective with Θ > Θ_c and hence α > 0. The dynamics in the convective layer would follow scaling laws analogous to those presented above, although the very small thermal expansion coefficient at temperatures only marginally warmer than Θ_c implies that the energy input and hence maximum flow speeds would be much smaller than the upper bound estimated above.

Zoom In Zoom Out Reset image size

The simulations of Zeng & Jansen (2021) used a relatively large "eddy diffusivity" (κ_z = 5 × 10⁻⁵m² s⁻¹), chosen such as to obtain a well-resolved stratified layer with negative thermal expansivity, above a convective deep ocean with weakly positive thermal expansivity (Figure 9(a),(b)). The vertical heat flux through the stratified layer is balanced by this "eddy diffusion" (Figure 9(c)), such that from the perspective of the model's energetics, the two terms in the integral in Equation (33) balance in the stratified layer. Below the stratified layer, the thermal expansivity is positive such that convective upward heat flux releases kinetic energy that can balance dissipation. Integrated over the entire ocean, however, diffusion is the only net energy source and balances the energy loss associated with net upward heat flux and the viscous dissipation (Figure 9(d)). The simulations of Kang et al. (2022) used an even larger eddy diffusivity (κ = 5 × 10⁻³ m² s⁻¹), such that the entire depth of the ocean becomes stratified with the upward heat flux and associated downward buoyancy flux accomplished by "eddy diffusion" (not shown). However, in both studies the model's diffusivities, which are much larger than the molecular value, need to be interpreted as representing turbulent mixing. If these turbulent mixing rates are real, the associated vertical heat flux represents turbulent advection, which amounts to a negative ${{ \mathcal B }}^{{\rm{\Theta }}}$ - that is a conversion of turbulent kinetic energy to large-scale potential energy. An additional implied energy source is therefore needed to generate the turbulent kinetic energy and justify the applied eddy diffusivities. A possible source are tides and/or librations, to which we will return below.

Tidal forcing generates kinetic energy in the ocean (via the second term on the lhs of Equation (4)). The tidal perturbation potential is typically approximately vertically constant throughout the depth of the ocean, thus generating primarily barotropic (i.e., vertically constant) tides. In Earth's ocean these barotropic tidal waves are largely linear until they encounter shallow shelves, where the tidal amplitude increases and most of the tidal energy dissipation is assumed to occur (Wunsch & Ferrari 2004). The icy moons do not have shelf seas similar to Earth's ocean, although tidal dissipation may nevertheless be significantly enhanced in specific regions (e.g., Hay & Matsuyama 2019). Linear tides in the deep ocean do not contribute significantly to the transport of heat, salt, and other properties in the ocean, such that global models of Earth's ocean circulation can reproduce realistic large-scale circulations and tracer distributions without explicitly accounting for tides. However, it is possible that tidal waves on some icy moons become sufficiently nonlinear to lead to substantial rectified mean flows (e.g., Huthnance 1981; Brink 2011). As the current manuscript is focused on buoyancy-driven circulation, we will not further consider the possible role of rectified tidal flows here, but note that this remains a potentially important topic for further research.

Despite weak direct interactions with the large-scale flow, ocean tides can have an important effect on the large-scale circulation in the presence of a statically stable stratification (e.g., driven by freshwater fluxes at the ice–ocean interface). Barotropic tides interacting with a rough seafloor topography in the presence of a statically stable stratification can transfer their energy into baroclinic tides, i.e., internal waves that propagate both horizontally and vertically and can be associated with substantial vertical shears (e.g., MacKinnon et al. 2017 and references therein). These shears in turn can lead to instabilities and eventually the generation of turbulence within the water column. This pathway is believed to be one of the main drivers of 3D turbulence in Earth's deep ocean (Vic et al. 2019). Although most of the the turbulent kinetic energy is dissipated into heat, some fraction, usually denoted by the "mixing efficiency," Γ, is converted to potential energy via vertical mixing of the stratified water column:

$$\begin{eqnarray}&&{{ \mathcal B }}^{t}=\int {(\overline{w^{\prime} b^{\prime} })}^{t}{dV}=-\int {\rm{\Gamma }}{\epsilon }^{t}{dV},\end{eqnarray} \tag{ 35 }$$

where $-{{ \mathcal B }}^{t}$ is the potential energy generation due to tidally driven turbulent mixing, ${(\overline{w^{\prime} b^{\prime} })}^{t}$ is the vertical buoyancy flux associated with tidally driven turbulence, ^t is the tidally generated turbulent kinetic energy dissipation rate per unit mass, and Γ ≲ 0.2 is the mixing efficiency (e.g., Peltier & Caulfield 2003).

The potential energy generated by turbulent mixing in a stratified ocean can then drive a large-scale circulation that converts the potential energy back to kinetic energy via a positive conversion $\overline{w}\overline{b}$ (where the overbars denote the large-scale circulation). If turbulent mixing is the dominant source of potential energy (i.e., when buoyancy gain and buoyancy loss from external forcing occur at approximately the same depth) the conversion $\overline{w}\overline{b}$ is, on average, approximately equal and opposite to the downward buoyancy flux associated with tidally driven turbulence:

$$\begin{eqnarray}&&{ \mathcal B }\equiv \int {wbdV}=\int \overline{w}\overline{b}{dV}+\int {(\overline{w^{\prime} b^{\prime} })}^{t}{dV}\equiv \overline{{ \mathcal B }}+{{ \mathcal B }}^{t}\approx 0,\end{eqnarray} \tag{ 36 }$$

where $\overline{{ \mathcal B }}$ is the kinetic energy generated (and ultimately dissipated) by the mean flow. The energy cycle can then be summarized as follows. Tides generate turbulent kinetic energy, a fraction of which is converted into large-scale potential energy via downward turbulent buoyancy flux (with the remainder being dissipated into heat). This potential energy can then be converted back to kinetic energy (and ultimately be dissipated) by the large-scale circulation. This mechanism requires a statically stable stratification such that turbulent mixing transports buoyancy downward. In the absence of a statically stable stratification, we still expect tidally driven turbulence to contribute to the mixing of properties, but we do not expect it to contribute to driving large-scale circulation.

In addition to tides, librations and injection of fluid through fissures in the ice or solid core may provide sources of kinetic energy. Librations are expected to play a similar role to tides, and as for tides, their effect on the large-scale circulation is expected to depend on their ability to generate turbulence in the ocean interior and on the presence of a statically stable stratification. Injection of water through fissures in the boundaries (which is not formally included in the global kinetic energy budget in Equation (4), where we assumed no-normal-flow boundary conditions) is likely to drive mechanical turbulence primarily in the vicinity of the boundaries, but strong jets emanating from the ice shell, as proposed by Kite & Rubin (2016), may play an important role in the dynamics of the upper ocean. ¹² A full investigation of the effects of such jets is beyond the scope of this study, but remains an interesting subject for future research.

We begin by estimating the potential role of tidally driven turbulence in the Snowball Earth ocean. Although little is known about tidal energy dissipation during the Snowball Earth periods, the average tidal dissipation between the Precambrian and the present was likely somewhat smaller but of similar order of magnitude as the present-day value (Williams 2000; Green et al. 2017). The direct effect of the Snowball Earth ice sheet on tides was likely small (Wunsch 2016). For scaling purposes we therefore here assume a tidal energy input on the order of 10¹² W, which is comparable to the present-day value of around 3.5 × 10¹² W (Wunsch & Ferrari 2004). Assuming a mean ocean depth of about 2 km, an ocean area of 3.5 × 10¹⁴ m², and density ρ₀ ≈ 1000 kg m⁻³, the average tidal energy dissipation per unit mass is around ^t ≈ 10⁻⁹ m² s⁻³, which is a little over an order of magnitude larger than the estimated energy input due to thermal forcing (see Table 2). If a significant fraction of the dissipated tidal energy is dissipated in the stratified ocean interior and is relatively efficiently converted to potential energy (i.e., Γ ≳ 10%) it thus may play a significant role in driving a large-scale ocean circulation. In particular, one may envision a scenario where freezing and melting at the ice–ocean interface combined with tidally driven turbulent mixing, which pushes the light fresh water into the ocean interior, can drive a significantly stronger large-scale ocean circulation than the thermal forcing alone. This scenario appears to be relevant for the interpretation of the simulation results by Ashkenazy et al. (2013) and Jansen (2016), although the turbulent vertical mixing in these simulations is parameterized with prescribed vertical diffusivities (i.e., ${(\overline{w^{\prime} b^{\prime} })}_{t}\to -\kappa {\partial }_{z}\overline{b}$)and is thus not constrained by tidal energy dissipation. The present analysis nevertheless suggests that the vertical mixing may be justifiable on energetic grounds given the expected tidally driven turbulence, although it would require that a significant fraction of the tidal energy is dissipated in the stratified ocean interior (as opposed to being highly localized in boundary layers where the stratification is vanishingly small) and is relatively efficiently converted to potential energy.

Table 2. Overview of Estimates for Mechanical Energy Input (Per Unit Mass) Due to Thermal Forcing, $\epsilon ={ \mathcal B }/H$, the Critical Length Scale Where Rotation is Expected to Become Important, L_rot, the Convective Rossby Number for Rotating Convection, Ro_rc, and an Estimate of the Maximum Characteristic KE Energy (Assuming Small Interior Dissipation) of the Large-scale Geostrophic flow, E_g , for Snowball Earth, Europa, and Enceladus

	(m2 s−3)	Lrot(m)	Rorc	Eg (m2 s−2)
Snowball Earth	2.5–5 × 10−11	7–10	2–5 × 10−3	10−4–10−3
Europa	5 × 10−13–8 × 10−12	8–30	1–8 × 10−4	5 × 10−5–3 × 10−3
Enceladus	≲2 × 10−13	≲1	≲1 × 10−4	≲10−4

Download table as:  ASCIITypeset image

Tidal energy dissipation on Europa has been estimated from modeling to be on the order of 10⁹–10¹¹ W (Chen et al. 2014; Matsuyama et al. 2018; Hay & Matsuyama 2019). Assuming, as before, an ocean depth on the order of 10⁵ m, an area of around 3 × 10¹³ m², and a density ρ₀ ≈ 1000 kg m⁻³, we estimate a tidal energy dissipation per unit volume of around ^t ≈ 3 × 10⁻¹³–3 × 10⁻¹¹ m² s⁻³, which is between about 1 order of magnitude smaller to 2 orders of magnitude larger than the energy input by buoyancy forcing (see Table 2). Depending on the estimate for tidal energy dissipation, and the fraction of the dissipated tidal energy that is converted to potential energy, tidally driven turbulent mixing thus may or may not play a significant role in driving a large-scale ocean circulation on Europa. In particular, freezing and melting at the ice–ocean interface combined with tidally driven turbulent mixing may drive a significantly stronger large-scale ocean circulation then the thermal forcing alone, if, and only if, ocean tidal dissipation is near the upper end of these estimates and a large fraction of that energy contributes to vertical mixing against a stable stratification. A significant salt-driven circulation was found in the simulations of Ashkenazy & Tziperman (2021), although turbulent vertical mixing is parameterized with a prescribed Earth-like vertical diffusivity in the simulations and is not constrained by tidal energy dissipation. Better estimates of tidal energy dissipation in Europa's ocean are needed to determine whether such a large turbulent vertical diffusivity can be justified for stratified regions of Europa's ocean.

For Enceladus' ocean, energy dissipation associated with the eccentricity and obliquity tides has been estimated to be relatively small, with a total dissipation between 10 and 10⁴ W suggested by Matsuyama et al. (2018) and Hay & Matsuyama (2019). More recently, Rovira-Navarro et al. (2023) argued that much larger tidal dissipation may occur in the presence of strong stratification and baroclinic tide resonances, although dissipation rates larger than 10⁶ W remain unlikely for the expected range of stratifications. Tyler (2020) has argued that much larger dissipation rates are possible, accounting for virtually all of the observed O(10¹⁰ W) heat flux, although the calculations generally predict the total dissipation in the ice–ocean system, and it is likely that most of the dissipation in the suggested scenarios would indeed occur in the ice sheet (see Beuthe 2016). Assuming that the most likely range for ocean tidal dissipation on Enceladus is ∼10–10⁶ W, and using an average ocean depth of around 3 × 10⁴ m and an area of around 5 × 10¹¹ m², we find an average tidal energy dissipation rate anywhere between around ^t ≈ 7 × 10⁻¹⁹–7 × 10⁻¹⁴ m² s⁻³. The upper end of this range remains more than an order of magnitude smaller than the estimated maximum mechanical energy input via thermal forcing (see Table 2).

However, Enceladus' ocean may be more strongly affected by librations of the ice shell, which is decoupled from the solid interior (Thomas et al. 2016). It is not clear whether librations lead to relatively modest dissipation confined to the ice–ocean boundary layer, thus not contributing significantly to interior mixing, or whether the excitation of internal waves or elliptical instability (an instability that can break up elliptical streamlines) can lead to relatively strong mixing throughout the depth of the ocean (Lemasquerier et al. 2017; Wilson & Kerswell 2018; Rekier et al. 2019; Soderlund et al. 2020). At the extreme end, Wilson & Kerswell (2018) argue that it is possible that all of the O(10 GW) of heating on Enceladus may be a result of librational dissipation, which, divided by the mass of Enceladus' ocean would give ^t ∼ 10⁻⁹ m² s⁻³—many orders of magnitude larger than the possible mechanical energy input from thermal forcing. In this case it is likely that libration-driven turbulence would play the dominant role in ocean mixing. Narrowing down the large uncertainty in the energy dissipation rate associated with libration-driven turbulence will be key to constraining the turbulent diffusivity, which is an essential parameter in numerical simulations of Enceladus' ocean (Zeng & Jansen 2021; Kang et al. 2022).