Internal energy dissipation in Enceladus's ocean from tides and libration and the role of inertial waves

TL;DR

This study models energy dissipation in Enceladus's subsurface ocean due to tides and libration, finding libration as the main contributor but still insufficient to explain observed heat, with inertial wave resonances potentially amplifying dissipation.

Contribution

It extends previous models by including libration effects and demonstrates how inertial wave resonances can enhance tidal dissipation in icy moons.

Findings

Libration dominates linear energy dissipation in Enceladus's ocean.

Libration-induced dissipation is about 0.001 GW, much less than observed 10 GW.

Resonances with inertial modes can significantly increase dissipation.

Abstract

Enceladus is characterised by a south polar hot spot associated with a large outflow of heat, the source of which remains unclear. We compute the viscous dissipation resulting from tidal and libration forcing in the moon's subsurface ocean using the linearised Navier-Stokes equation in a 3-dimensional spherical model. We conclude that libration is the dominant cause of dissipation at the linear order, providing up to about 0.001 GW of heat to the ocean, which remains insufficient to explain the (about) 10 GW observed by Cassini. We also illustrate how resonances with inertial modes can significantly augment the dissipation. Our work is an extension to Rovira-Navarro et al. [2019] to include the effects of libration. The model developed here is readily applicable to the study of other moons and planets.

Tables

Table 1. Table 1: List of physical parameters

Name	Symbol	Core	Ocean	Crust
Outer radius (km)	R	192	[192 - 252]	252
Mass Density (kg/${\text{m}}^3$)	$\rho$	[2483 - 2357]*	1020	920
Shear modulus (GPa)	$\mu$	40	0	3.5

Equations

$${\mathcal{D}}_{\text{visc}}=2\text{Ek}{\int }_{\mathcal{V}}\stackrel{⌢}{\nabla \mathbf{v}}:\stackrel{⌢}{\nabla \mathbf{v}},$$

$${\nabla }^2\delta \varphi$$

$$-{\omega }^2\rho \mathbf{s}$$

$$\mathbf{\Delta }\mathbf{\sigma }\equiv -p\mathrm{g}+\mu \left(\nabla \mathbf{s}+\nabla {\mathbf{s}}^{\text{T}}\right),$$

$${\left[\delta \varphi \right]}_{-}^{+}$$

$${\left[\hat{\mathbf{n}}\cdot \left(\nabla \delta \varphi +4\pi G\rho \mathbf{s}\right)\right]}_{-}^{+}$$

$$[\hat{\mathbf{n}}\cdot \mathbf{\Delta }\mathbf{\sigma }{]}_{-}^{+}$$

$$[\hat{\mathbf{n}}\cdot \mathbf{s}{]}_{-}^{+}=0.$$

$$\delta \varphi =\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }\delta {\varphi }_{\ell ,m}{\text{Y}}_{\ell }^m(\theta ,\lambda ),$$

$$\delta {\varphi }_{2,-2}^{\text{tide}}$$

$$\delta {\varphi }_{2,0}^{\text{tide}}$$

$$\delta {\varphi }_{2,2}^{\text{tide}}$$

$$r(\theta ,\lambda )=R\{1+{\alpha }_0{\text{Y}}_2^0(\theta ,\lambda )+{\alpha }_2[{\text{Y}}_2^2(\theta ,\lambda )+{\text{Y}}_2^{-2}(\theta ,\lambda )]\},$$

$$\frac{\dot{r}}{R}=2i{\alpha }_2\dot{\delta }[{\text{Y}}_2^2(\theta ,\lambda )-{\text{Y}}_2^{-2}(\theta ,\lambda )].$$

$$\frac{\dot{r}}{R}=2{\alpha }_2ϵ(i\omega )[{\text{Y}}_2^2(\theta ,\lambda )-{\text{Y}}_2^{-2}(\theta ,\lambda )]\cos (\omega t).$$

$$i\omega \mathbf{v}+2\hat{\mathbf{z}}\times \mathbf{v}+\nabla p-\text{Ek}{\nabla }^2\mathbf{v}=0,$$

$$\text{Ek}\equiv \frac{\nu }{\Omega R_o^2},$$

$${\mathbf{v}}_{\text{tide/libration}}=i\omega {\mathbf{s}}_{\text{tide/libration}}.$$

$$\mathbf{v}={\mathbf{v}}_{\text{tide/libration}}.$$

$$\mathbf{v}=\nabla \times \nabla \times (P\mathbf{r})+\nabla \times (T\mathbf{r}),$$

$$\lambda \mathbf{u}+2\hat{\mathbf{\Omega }}\times \mathbf{u}+\nabla p-\text{Ek}{\nabla }^2\mathbf{u}=0,$$

$$\mathbf{v}\sim \sum _{\alpha }\frac{<\mathbf{f}|{\mathbf{u}}_{\alpha }>}{|{\lambda }_{\alpha }-i\omega |}{\mathbf{u}}_{\alpha },$$

$$\text{Le}=\frac{B_0}{\Omega R_o\sqrt{\mu \rho }}$$

$$m=-2\to$$

$$m=0\to$$

$$m=2\to$$

$$m=-2\to$$

$$m=0\to$$

$$m=2\to$$

$$m=-2\to$$

$$m=0\to$$

$$m=2\to$$

$$m=-2\to$$

$$m=0\to$$

$$m=2\to$$

$$\text{top}\to$$

$$\text{bottom}\to$$

Full text

1 Introduction

From the massive amount of information collected by Cassini, Saturn’s moon Enceladus currently appears as one of the most habitable moons in the Solar System. The presence of a subsurface reservoir of liquid water was inferred soon after the first few flybys (Porco et al., 2006). Measurements of the gravity field deduced from subsequent flybys helped establish the large spatial extent of this reservoir (Iess et al., 2014), but only libration observations could provide definitive evidence for a global-scale ocean rather than a regional sea (Thomas et al., 2016).

Tidal dissipation is the most likely power source for the observed geological activity (Nimmo et al., 2018). Current models, however, have so far been unable to account for the observed $\sim 10$ GW of endogenic heat flow emanating from the south polar region (Spencer et al., 2006; Howett et al., 2011) without invoking the presence of a highly porous (‘fluffy’) solid inner core (Roberts, 2015; Choblet et al., 2017) or a convecting icy crust (Běhounková et al., 2010, 2017). As suggested by Barr and McKinnon (2007), this latter hypothesis appears unlikely in view of the relatively small thickness of the crust ($\sim 20$ km) deduced from geodetic observations (Beuthe et al., 2016).

If dissipation is not concentrated in the core or in the shell, one possibility is that it takes place predominantly in the ocean. Most of the current models focusing specifically on the ocean layer rely on the solution of the Laplace Tidal Equations (LTE) whereby the fluid is modelled as a 2-dimensional thin layer (Tyler, 2014; Hay and Matsuyama, 2017; Beuthe et al., 2016; Matsuyama et al., 2018). However, gravity and topography data suggest that the average thickness of the ocean is not negligible, $\sim 38$ km if the crustal topography is isostatically supported, which is also consistent with the thin crust inferred from libration (Beuthe et al., 2016; Hemingway and Mittal, 2019). Such a thick ocean challenges the validity of the conclusions drawn from the 2-d models.

Additionally, the action of the Coriolis force on planetary oceans is known to support the existence of inertial waves within their volume (Poincaré, 1885). Morize et al. (2010) have proven that inertial waves can be excited through tidal deformation in their laboratory experiment. The role they play in planetary and astrophysical settings is still, however, far from clear. When viscosity is taken into account, the flow associated with these waves develops regions of intense shear within the fluid volume (Greenspan, 1968). These internal shear layers can significantly increase the total amount of dissipation in the ocean. The formula giving the total viscous dissipation for a flow of velocity field $\boldsymbol{v}$, is

[TABLE]

with $\mathord{\stackrel{{\scriptstyle\smallfrown}}{{\boldsymbol{\nabla}\boldsymbol{v}}}}\equiv\frac{1}{2}(\boldsymbol{\nabla}\boldsymbol{v}+\boldsymbol{\nabla}\boldsymbol{v}^{\text{T}})$ and where Ek denotes the Ekman number, which is typically very small in planetary applications ($\text{Ek}\sim 10^{-13}$ for Enceladus). In order to significantly contribute to the total dissipation, the gradient of velocity within the ocean volume in general and in the internal shear layers in particular, must be sufficiently strong to compensate the smallness of Ek. A study of the problem has recently been conducted by Rovira-Navarro et al. (2019) who concluded to the irrelevance of inertial waves in the energy budget. We reproduce their results in the present paper and extend them to include the contribution from libration which corresponds to a forcing with a much larger amplitude. We also take into account the presence of the icy crust. As the obliquity of Enceladus is very small (Baland et al., 2016), we focus our attention on the effects of eccentricity tides, which are likely to dominate over obliquity tides in the ocean layer (Chen and Nimmo, 2011; Tyler, 2014), even though, recent work by Hay and Matsuyama (2019) has shown that non-linear effects can alter this picture at least within the simplified model based on the LTE.

The present paper is structured as follows. In Sec. 2, we present our model including the details on the equilibrium state and the equations of motion and how these are used to compute the response to tidal and libration forcing. The results are presented in Sec. 3 and discussed in Sec. 4. We provide details about the interpretation of the forced response in terms of resonance with inertial waves in this last section and conclude with perspectives about the possible continuations of the present work.

The generality of the method described below makes it readily applicable to the study of other moons and planets.

2 Method

We can only reach the very low values of Ek in a sphere as it is not otherwise possible to model the very thin viscous layer at the ocean’s boundary. For this reason, we proceed in two steps. First we compute the tidal deformation at the ocean’s boundary of a non-rotating spherical model (Sec. 2.1). Second, we compute the resulting flow inside the ocean and with rotation taken into account (Sec. 2.2).

2.1 Deformation

2.1.1 Equilibrum state

We model the moon as a set of 3 homogeneous concentric spherical shells and we assume that each layer is incompressible. Table 1 gives the set of parameters used in the present paper. We set the size of the core to the central value 192 km obtained by Beuthe et al. (2016) with the reference datasets, and consider a range of interior models with variable ocean thickness.

In what follows, we use $\eta=\frac{R_{\text{core}}}{R_{\text{ocean}}}$ to represent the aspect ratio of the spherical shell.

2.1.2 Equations of deformation

The equation governing the oscillation around the equilibrium state are the Poisson equation and the conservation of momentum. Assuming incompressibility, these read :

[TABLE]

where $\phi_{0}$ is the equilibrium gravity potential, $\delta\phi$ denotes the (Eulerian) increments of gravity (i.e. self-gravitational + tidal) potential, $\boldsymbol{s}$ is the displacement vector field and $\boldsymbol{\Delta\sigma}$ denotes the (Lagrangian) increment of Cauchy stress (Dahlen and Tromp, 1998). The latter can be expressed in terms of the displacement and (Lagrangian) increment of pressure :

[TABLE]

where $\boldsymbol{\mathrm{g}}$ is the metric tensor. $\mu$ denotes the shear modulus, which is assumed to be constant in each layer, its values are listed in Table 1, as well as those for the density, $\rho$.

The gravity potential, displacement vector and stress tensor must satisfy the following conditions at each undeformed spherical boundary (Dahlen and Tromp, 1998) :

[TABLE]

where $\boldsymbol{\hat{n}}$ is the normal vector and the notation $\left[\cdot\right]^{+}_{-}$ denotes the difference between the values of the enclosed quantity on both sides of the boundary. The last constraint is that the normal component of the displacement be continuous across each internal boundary :

[TABLE]

The above does not hold at the free outermost boundary.

2.1.3 Tidal deformation

Enceladus gets deformed under the effect of tides caused by its slightly eccentric orbit around Saturn (eccentricity tides). The corresponding (Eulerian) increment of gravity potential can be decomposed in terms of spherical harmonics :

[TABLE]

where $\theta$ and $\lambda$ denote the colatitude and (east) longitude respectively.111We use the convention that $\text{Y}_{\ell}^{m}(\theta,\lambda)\equiv\sqrt{\frac{(\ell-m)!}{(\ell+m)!}}\text{P}_{\ell}^{m}(\cos\theta)e^{im\lambda}$

and $\text{P}_{\ell}^{m}(x)=\frac{(-1)^{m}}{2^{\ell}\ell!}(1-x^{2})^{\frac{m}{2}}\frac{d^{\ell+m}}{dx^{\ell+m}}(x^{2}-1)^{\ell}$, where $\text{P}_{\ell}$ are the Legendre polynomials. The only nonzero components for the contribution due to tides caused by the eccentricity of Enceladus’ orbit are

[TABLE]

where $\omega=\frac{2\pi}{1.37}~{}\text{days}^{-1}$ is the (diurnal) frequency of the forcing and $e=4.7\times 10^{-3}$ denotes the orbital eccentricty.

We compute the resulting displacement at the top and bottom of the ocean using Eq. (2) & (3) for different values of $\eta$. More details are given in Appendix A.1. Fig. 1 shows the spatial pattern for each tidal component for $\eta=0.838$ corresponding to a shell thickness of 23 km, as suggested by models of isostasy (Beuthe et al., 2016).

2.1.4 Libration forcing

Enceladus is not, in reality, spherical. And this fact influences its rotation. In particular, the equilibrium deformation caused by permanent tides induces a periodic motion of longitudinal libration under the perturbing effect of gravitational torques from external sources and the restoring effect of internal pressure and gravitational torques (inertial torques).

In the present paper, we treat the libration motion as a given and compute the resulting dissipation in the fluid layer of our simplified spherical setting. Our model of libration is based on the rigid-shell model of Baland and Van Hoolst (2010); we therefore neglect the limited decrease of libration amplitude introduced by elastic deformation (Van Hoolst et al., 2016). For the present, it is sufficient to observe that the libration of the ocean’s flattened boundaries can be represented as a superposition of several oscillations of fictitious spherical boundaries: a toroidal degree-1 oscillation where the librating spherical boundaries tend to viscously drag the ocean, and a spheroidal degree-2 oscillation where the librating flattened boundaries tend to push the ocean fluid, as if spherical boundaries underwent radial deformation.

The tangential displacement can be obtained directly from the libration amplitude as a function of $\eta$, the aspect ratio of the ocean (see Fig. 8b of Appendix A.2). It results in what is later referred to as the $m=0$ component.

The non-axisymmetric shape of Enceladus also causes a radial displacement at the top and bottom of the ocean, which can be estimated in the following manner. To first order, the radial coordinate of the triaxial boundary can be approximated as

[TABLE]

where $R$ is the mean radius and $\alpha_{0}$ and $\alpha_{2}$ are real numbers describing the polar and equatorial flattening, which we assume to be hydrostatic and compute using the first-order theory of figures. In the frame rotating at constant angular velocity $\boldsymbol{\Omega}$, libration can be accounted for by replacing : $\lambda\rightarrow\lambda+\delta(t)$ with $\delta\ll 1$. The first two terms in Eq. (13) do not depend on the azimuthal angle, $\lambda$, and so the velocity of the moving boundary comes out to be

[TABLE]

Now, if we set $\delta(t)=\epsilon\sin(\omega t)$, where $\epsilon$ and $\omega$ denote the amplitude and the (diurnal) frequency of libration respectively, we obtain

[TABLE]

The values of the (hydrostatic) flattening parameter $\alpha_{2}$ and the amplitude of libration $\epsilon$ both depend on the aspect ratio $\eta$. We provide simplified expressions of these in Appendix A.2. Hereafter, the radial displacement is referred to as the $m=\pm 2$ components.

2.2 Viscous dissipation

We use the (linearised) Navier-Stokes equation to model the motion of the ocean. In its dimensionless form and in the frame rotating with angular velocity $\boldsymbol{\Omega}~{}\equiv~{}\Omega~{}\boldsymbol{\hat{z}}$, this reads :

[TABLE]

where $\boldsymbol{v}$ denotes the velocity and $p$ denotes the reduced pressure. The Ekman number parametrises the balance between the viscous force and the Coriolis force. We define it as such

[TABLE]

where $\nu$ denotes the kinematic viscosity (here taken to be that of water: $\nu=10^{-6}~{}\text{m}^{2}/\text{s}$) and $R_{o}$ is the (mean) radius at the top of the ocean. In planetary settings, Ek is typically very small. For Enceladus, its value is of the order222There is actually a lot of uncertainty on the precise value of this parameter. Here we are only interested in its order of magnitude. $\text{Ek}\sim 3.5\times 10^{-13}$.

The motion inside the ocean is forced by the deformation of the boundary caused by tides and libration. We recover the velocity of the moving boundary from the displacement vector $\boldsymbol{s}$ via

[TABLE]

We use the no-slip boundary condition, which, in this case, amounts to impose that the velocity is continuous across the boundary (at the top and bottom of the ocean) :

[TABLE]

To first order, it is sufficient to enforce this condition at the boundary of the (spherical) equilibrium figure. Owing to the condition of incompressibility ($\boldsymbol{\nabla}\cdot\boldsymbol{v}=0$) we can write the velocity field as

[TABLE]

where $P$ and $T$ denote the poloidal and toroidal potentials respectively. These are then decomposed in terms of their spherical harmonics coefficients, $P_{\ell,m}$ and $T_{\ell,m}$. The two orthogonal projections of Eq. (16) that we use are given in Appendix C in terms of their spherical harmonics projections as Eqs. (39) and (40).

We solve these equations numerically using the techniques presented in Sec. 2 of Rekier et al. (2018) and we compute the dissipation using Eq. (1).

3 Results

3.1 Tides

Fig. 2 shows the dissipation for $m\in\{-2,0,2\}$ and the total dissipation for all these contributions (lower-right panel). There is a general trend towards larger dissipation as $\eta$ decreases. This is consistent with the results of Matsuyama et al. (2018) as this regime corresponds to a thicker ocean and a thinner icy crust in our model. The principal contribution to the total dissipation comes from the prograde tide ($m=-2$) as expected considering it is the dominant contribution to the tidal potential. The total contribution from the eccentricity tides lies somewhere in the range $10^{-6}<\mathcal{D}_{\text{visc}}<10^{-3}$ GW depending on the presence or absence of a peak.

The peaks in the dissipation profile are due to the augmented contribution caused by the internal shear layers inside the ocean. This can be observed from the plots on Fig. 3a & 3b, which show the density of kinetic energy for two different values of $\eta$ (both within the the range $0.82<\eta<0.85$ obtained by Beuthe et al. (2016)) one associated to larger dissipation and one to smaller dissipation. The corresponding values of $\mathcal{D}_{\text{visc}}$ can be read from Fig. 3c and differ by three orders of magnitude. The pattern of internal shear layers appears much sharper on Fig. 3b, corresponding to the higher dissipation while they have a lower intensity and seem to fade out more quickly on Fig. 3a.

3.2 Libration

Fig. 4 shows the dissipation for $m\in\{-2,0,2\}$ and the total dissipation for all these contributions (lower-right panel). We observe the same trend towards larger dissipation for decreasing $\eta$, similar to the situation with eccentricity tides.

The contributions $m=2$ and $m=-2$ lead to identical dissipation profiles, as expected given the symmetrical nature of the forcing. The contribution $m=0$ dominates for larger values of the viscosity but decreases rapidly as Ek goes down. This is the only contribution that shows no trace of peaks in its profile which indicates that most of the dissipation takes place inside the Ekman boundary layer. This picture is comforted by noting that the dissipation decreases with decreasing viscosity as $\mathcal{D}_{\text{visc}}\sim\text{Ek}^{1/2}$, as shown in Fig. 5 for $\eta=0.838$, a scaling law typical of Ekman boundary layers. Extrapolation of this curve to $\text{Ek}=3.5\times 10^{-13}$ gives $\mathcal{D}_{\text{visc}}\sim 0.001~{}\text{GW}$ for Enceladus. This is the dominant source of dissipation in the moon’s ocean if one disregards the other, more erratic, contributions. The total dissipation due to libration (Fig. 4, lower-right pannel) also shows how the three components may become commensurable at low Ekman number ($\text{Ek}\leq 10^{-10}$) at values of $\eta$ where there is a peak.

4 Discussion

4.1 Libration as the dominant source of dissipation

Toroidal libration ($m=0$) dominates over all other sources of dissipation for $\text{Ek}>10^{-10}$, as can be seen on Fig. 2 & 4, and becomes commensurable to the other components of libration for lower viscosities. Since the amount of dissipation induced by the toroidal drag scales nicely with the Ekman number, we can extrapolate the amount of dissipation to the Ekman number relevant to Enceladus from Fig. 5. The power scales as $\sim\text{Ek}^{1/2}$, which indicates the predominant role of the Ekman boundary layer. Extrapolating this power law to $\text{Ek}=3.5\times 10^{-13}$, the relevant value for Enceladus, we predict a value of $10^{-3}$ GW for the total dissipation. That is about 4 orders of magnitude below the $\sim 10$ GW observed by Cassini.

4.2 The role of inertial modes

Our results have revealed how the presence of internal shear layers can lead to a significant increase of the internal dissipation. Previous studies have attributed such peaks in the power profile to the existence of wave attractors (Rieutord et al., 2000; Ogilvie, 2013; Rovira-Navarro et al., 2019). We would like to complement this picture with our own interpretation based on the spectrum of free inertial modes.

In the absence of external body force, inertial modes are solutions to

[TABLE]

which is analogous to Eq. (16), except that $\lambda$ now denotes the complex eigenvalue. It is sensible to assume that the response to a body force (per unit volume), $\boldsymbol{f}$, of frequency $\omega$ can be decomposed, at least partially, onto the (infinite) set of inertial modes $\{\boldsymbol{u}_{\alpha}\}$ :

[TABLE]

where $<\cdot|\cdot>$ denotes the projection operator over pairs of vector fields. Ivers et al. (2014) and Backus and Rieutord (2017) have demonstrated the completeness of the above modal expansion for an inviscid fluid inside a spherical or ellipsoidal container (with no inner core). In analogy with those works, we define the projection operator as $<\boldsymbol{v}|\boldsymbol{w}>\equiv\text{Re}\int_{\mathcal{V}}\boldsymbol{v}^{*}\cdot\boldsymbol{w}$. A resonance takes place when the factor multiplying one or more $\boldsymbol{u}_{\alpha}$ becomes large, which happens when the distance between the forcing frequency and the eigenvalue approaches zero and/or $<\boldsymbol{f}|\boldsymbol{u}>$ is large. This typically happens when the spatial structure of the forcing shares some similarity with an eigenmode. This formalism can be applied to our case after we have found the relevant body force, $\boldsymbol{f}$. Details on how to do so are explained in Appendix B.

Fig. 6 illustrates the resonance with inertial waves. The upper and middle panels represent the evolution of the spectrum around $\lambda=i$ as a function of $\eta$ ($\text{Ek}=10^{-7}$). The colour scale is inherited from the middle plot with lighter colours corresponding to a stronger damping. The lower panel shows a superposition of the dissipation profile (in black) on top of the combination $\frac{<\boldsymbol{f}|\boldsymbol{u}_{\alpha}>}{|\lambda_{\alpha}-i|}$. There is a clear correlation between the two; the peaks in dissipation profile appear where the resonances with the eigenmodes are the strongest. The trend towards more dissipation in the deep ocean limit can be attributed to a large amount of weak resonances in this picture.

4.3 Conclusion and Future work

We have computed the viscous dissipation in Enceladus’s ocean and we have shown that the combined effect of eccentricity tides and libration is not sufficient to explain the heat flux observed by Cassini. We therefore confirm that the results of Matsuyama et al. (2018) remain valid even when the dissipation of the internal shear layers in the ocean’s bulk are taken into account. We show that librationally-induced dissipation dominates tidally-induced dissipation, but the overall power dissipated in the ocean is still negligible, as in Rovira-Navarro et al. (2019).

In light of our results, it is clear that the heat flux observed on Enceladus cannot be explained solely by viscous dissipation in the ocean alone, at least at the linear level. The source of heat must therefore be sought elsewhere. The ohmic dissipation in the shear layers due to the presence of Saturn’s magnetic field appears as a relevant candidate at first glance. However, the external magnetic field at Enceladus is quite small, of the order $\sim 0.6~{}\mu$T, which corresponds to a Lehnert number of $\text{Le}\sim 1.6\times 10^{-6}$. Lin and Ogilvie (2018) have shown that ohmic dissipation dominates over viscous dissipation only when $\text{Le}>\text{Ek}_{m}^{2/3}$, where $\text{Ek}_{m}$ is the magnetic Ekman number333The Lehnert and magnetic Ekman numbers are defined as

$\displaystyle\text{Le}=\frac{B_{0}}{\Omega R_{o}\sqrt{\mu\rho}}$ $\displaystyle\text{Ek}_{m}=\frac{\nu_{m}}{\Omega R_{o}^{2}}$

where $B_{0}$ is the background magnetic field, $\mu$ is the magnetic permeability and $\nu_{m}$ is the magnetic diffusivity of the fluid.. Recent estimates of the electrical conductivity inside the ocean (Vance et al., 2018) give $\text{Ek}_{m}\sim 1$. The ohmic dissipation is therefore unlikely to contribute significantly to the total.

Wilson and Kerswell (2018) have argued that libration could potentially generate turbulence inside the ocean and thus significantly increase the amount of dissipation.

Another possibility is that the present heat production does not balance the present heat flow, but was larger in the past during periods of larger eccentricities (Ojakangas and Stevenson, 1986).

The observed North-South dichotomy of the moon’s surface poses another important puzzle. The southern hemisphere is young and re-surfaced, while the northern hemisphere is old and cratered. This problem cannot be addressed in our simplified spherical model. The observed $J_{3}$ implies that the icy crust is thinner over the south polar terrain (Iess et al., 2014). We also know that the density of kinetic energy inside the shear layers is higher at the poles, close to the axis of rotation (Rieutord and Valdettaro, 1997), something that is visible on Fig. 3. It would be interesting to see how the dissipation at both poles changes when one takes the topography of the ocean into account. On the one hand, the internal shear layers tend to be more focused in the thin ocean limit. On the other hand, the overall dissipated power tends to be higher in the thick ocean limit. Both regimes are illustrated on the left and right sides of every dissipation profile in the present paper (e.g. Fig. 2 & 4).

Acknowledgement

The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Advanced Grant agreement No 670874).

Appendix A Tidal and libration forcings

A.1 Tidal displacements

The three spherical harmonics components of the radial and tangential displacements are shown on Fig. 7 for different values of the aspect ratio, $\eta$.

The analytical expressions for the displacement are quite long and impractical. Instead, we provide a set of simplified expressions obtained using a Pade approximant on the interval $\eta\in[\frac{16}{21},1]$. All expressions are in dimensionless units :

radial displacement at the top of the ocean

[TABLE]

radial displacement at the bottom of the ocean

[TABLE]

tangential displacement at the top of the ocean

[TABLE]

tangential displacement at the bottom of the ocean

[TABLE]

A.2 Libration forcing

The value of the parameters $\alpha_{2}$ and $\epsilon R$ at the top and bottom of the ocean are represented on Fig. 8. The analytical expressions for the amplitude $\epsilon R$ are quite long and impractical. And so we provide an approximation on the interval $\eta\in[\frac{16}{21},1]$ based on a Pade approximant (all expressions are in dimensionless units) :

Amplitude at the top and bottom of the ocean

[TABLE]

Appendix B Tidal dissipation with a body force

Some authors [Rieutord et al., 2000, Ogilvie, 2013, Rovira-Navarro et al., 2019] presented a slightly different formalism to compute the dissipation caused by tides. We now want to show the equivalence to ours. Their approach is based on the decomposition of the velocity field in two parts :

[TABLE]

The first part is attributed to the displacement caused by the equilibrium tides. Plugging the above ansatz into the momentum equation, leads to an equation for $\boldsymbol{v}^{\text{(d)}}$ :

[TABLE]

where $\boldsymbol{f}=-i\omega\boldsymbol{v}^{\text{(eq)}}-2\boldsymbol{\hat{\Omega}}\times\boldsymbol{v}^{\text{(eq)}}$ is treated as a body force per unit volume. Ogilvie [2013] argues that when $\boldsymbol{v}^{\text{(eq)}}$ is computed from the theory of equilibrium tides, one has $\boldsymbol{v}^{\text{(eq)}}=\boldsymbol{\nabla}X$, where $X$ is an harmonic scalar function ($\boldsymbol{\nabla}^{2}X=0$, owing to the condition of incompressibility). The harmonic component of degree $\ell$ therefore reads $X_{\ell,m}=(A_{\ell,m}r^{\ell}+\frac{B_{\ell,m}}{r^{\ell+1}})$, where $A_{\ell,m}$ and $B_{\ell,m}$ are scalar constants that depend on the deformation (see Ogilvie [2013] for details).

In our approach, we do not decompose the velocity field explicitly and rather solve for the full momentum equation, including the potential perturbation. This extra term, however, can be merged with the reduced pressure, which, in turn, disappears completely when one solves for the vorticity in the volume. The effect of the potential thus appears only in the expression of the boundary condition. This is the reason for which, the condition of continuity of the total flow, $\boldsymbol{v}$, across the physical boundaries is given by our Eq. (19) while Rovira-Navarro et al. [2019] simply have $\boldsymbol{v}^{(d)}=0$.

The two methods lead to results that are completely equivalent. The approach based on the boundary forcing is perhaps conceptually simpler and better suited to the usual language of fluid dynamicists. The description based on the body force $\boldsymbol{f}$, on the other hand, has the advantage to be directly usable in the discussion presented in Sec. 4.2.

Fig. 9 is a reproduction of Fig.6 (b) in Rovira-Navarro et al. [2019] using our settings. The discrepancy between the predicted amounts of dissipation in the deep ocean limit (left part of our plot, right part of theirs) is due to the presence of an icy crust in our model while it is neglected in their study. This leads to a decrease in the amount of energy dissipated in agreement with the results of Matsuyama et al. [2018]. Our curves also go up slightly in the thin ocean limit on the right part of the plot. This feature is absent in Rovira-Navarro et al. [2019] and is likely due to the fact that they use a purely radial forcing, thus neglecting tangential displacement at the boundary. When we take these into account, resonances appear in the thin ocean limit consistently with the results of Ogilvie [2009]. Appart from these differences, our curves bear a striking ressemblance with theirs, which illustrates the equivalence of our approaches.

Appendix C Navier-Stokes equation in spherical harmonics

Here below, we provide the expressions used to compute the flow inside the ocean. These are based on the curl of Eq. (16), i.e. the equation of vorticity, and the ansatz Eq. (20) :

[TABLE]