# Do tidally-generated inertial waves heat the subsurface oceans of Europa and Enceladus?

Marc Rovira-Navarro, Michel Rieutord, Theo Gerkema, Leo R. M. Maas, Wouter van Der Wal, Bert Vermeersen

► To cite this version:

Marc Rovira-Navarro, Michel Rieutord, Theo Gerkema, Leo R. M. Maas, Wouter van Der Wal, et al.. Do tidally-generated inertial waves heat the subsurface oceans of Europa and Enceladus?. Icarus, 2019, 321, pp.126-140. 10.1016/j.icarus.2018.11.010. insu-03674454

HAL Id: insu-03674454

https://insu.hal.science/insu-03674454v1

Submitted on 29 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Royal Netherlands Institute for Sea Research

This is a pre-copyedited, author-produced version of an article accepted for publication, following peer review.

Rovira-Navarro, M.; Rieutord, M.; Gerkema, T.; Maas, L.R.M.; van der Wal, W. & Vermeersen, B. (2019). Do tidally-generated inertial waves heat the subsurface oceans of Europa and Enceladus? *Icarus*, 321, 126-140

Published version: https://dx.doi.org/10.1016/j.icarus.2018.11.010

Link NIOZ Repository: http://www.vliz.be/nl/imis?module=ref&refid=308746

[Article begins on next page]

The NIOZ Repository gives free access to the digital collection of the work of the Royal Netherlands Institute for Sea Research. This archive is managed according to the principles of the [Open Access Movement](https://www.open Access movement.org), and the [Open Archive Initiative](https://www.openARCHIVE.org). Each publication should be cited to its original source - please use the reference as presented.
When using parts of, or whole publications in your own work, permission from the author(s) or copyright holder(s) is always needed.
# Do tidally-generated inertial waves heat the subsurface oceans of Europa and Enceladus?

Marc Rovira-Navarro<sup>a,b,*</sup>, Michel Rieutord<sup>c,d</sup>, Theo Gerkema<sup>e</sup>, Leo R.M. Maas<sup>f,a</sup>, Wouter van der Wal<sup>b</sup>, Bert Vermeersen<sup>b,e</sup>

<sup>a</sup>NIOZ Royal Netherlands Institute for Sea Research, Department of Ocean Systems OCS, and Utrecht University, P.O. Box 59, 1790 AB Den Burg, Texel, the Netherlands

<sup>b</sup>TU Delft, Faculty of Aerospace Engineering (AE), Building 62 Kluyverweg 1, P.O Box 5058, 2629 HS Delft / 2600 GB Delft, The Netherlands

<sup>c</sup>Université de Toulouse, UPS-OMP, IRAP, Toulouse, France

<sup>d</sup>CNRS, IRAP, 14, avenue Edouard Belin, 31400 Toulouse, France

<sup>e</sup>NIOZ Royal Netherlands Institute for Sea Research, Department of Estuarine and Delta Systems EDS, and Utrecht University, P.O. Box 140, 4400 AC Yerseke, the Netherlands

<sup>f</sup>Faculty of Science, Department of Physics, Institute for Marine and Atmospheric Research Utrecht, Utrecht University P.O Box 80.005, 3508 TA Utrecht

## Abstract

Some of the moons of the outer solar system harbour subsurface liquid oceans. Tidal dissipation plays an important role in preventing these oceans from freezing. In the past, most studies considered only tidal dissipation in the solid layers of these bodies (rock and ice). Recently, new studies considering tidal dissipation in the oceans of these moons have appeared. All of them make use of the shallow water approximation. However, the use of this approximation might not be adequate. Here we consider the linear non-hydrostatic three dimensional response of these oceans to tidal forcing with the full Coriolis force. To do so we consider an ocean of homogeneous density contained within a perfectly spherical shell and neglect the effect of the ice shell. We force the ocean with a time changing tidal potential and observe patterns of periodic inertial waves that take energy from the global tidal forcing and focus it along thin shear layers that propagate in the fluid. We focus on Europa and Enceladus, showing that inertial waves result in fluid flows of significant amplitude (a few cm/s). Nevertheless, we find that under the previously mentioned assumptions tidal dissipation due to inertial waves is several orders of magnitude smaller than Europa's radioactive heating and Enceladus' observed heat flux. Finally, we propose additional dissipation mechanisms that might play a relevant role in Europa and Enceladus and could be further investigated.

Keywords: Tides, Enceladus, Europa, Rotational dynamics

<footnote>*Corresponding author
Email address: marc.rovira@nioz.nl (Marc Rovira-Navarro)</footnote>
1. Introduction

Evidence for underground water reservoirs in some of the moons of the outer solar system has accumulated over the last few decades. The existence of a subsurface water body in Europa was first hypothesised by Cassen et al. (1979). They argued that radiogenic and tidal heating could melt Europa's interior and form an ocean. The hypothesis was confirmed when Khurana et al. (1998) and Kivelson et al. (2000) reported variations of the magnetic field in Europa associated with an induced magnetic field and showed that a layer of subsurface salty water is consistent with these observations. The existence of a subsurface ocean is also consistent with the complex geology of Europa's surface (e.g., Pappalardo et al., 1999).

The case of Enceladus is markedly different. This tiny moon of Saturn, its radius being just over 250 km, features vigorous geological activity. Cassini flybys revealed water plumes emanating from long parallel cracks (nicknamed tiger-stripes) on the South Pole (Porco et al., 2006; Spencer et al., 2006). The detection of salt-rich grains in the plumes (Postberg et al., 2009, 2011) and evidence of hydrothermal activity within the moon (Hsu et al., 2015) indicate that the plumes originate from an underground water reservoir. The observed forced libration of Enceladus' surface reveals that the ocean is not only restricted to the moon's southern pole but it completely decouples the ice shell from the moon's rocky mantle (Thomas et al., 2016).

Far away from the Sun, these oceans cannot be maintained by solar irradiation. Another heat source is therefore needed to prevent them from freezing; this heat source is most likely tidal dissipation. Due to orbital resonances with Io and Ganymede, Europa orbits Jupiter in an eccentric orbit. Similarly, in the Saturnian system, Enceladus' eccentricity is maintained by an orbital resonance with Dione. The orbital eccentricity results in a time-varying tidal potential that raises a prominent tide. The moons' obliquity results in additional latitudinal librations of the tidal bulge. Most studies focused on studying tidal dissipation in the solid layers of Europa and Enceladus (e.g., Ojakangas and Stevenson, 1989; Roberts and Nimmo, 2008). These studies show that the Europan ocean can be maintained by tidal and radiogenic heating; thermal models suggest that the ocean might be around 100 km thick (Hussmann et al., 2002). This is, however, not the case of Enceladus, where tidal and radiogenic heating in the solid parts of the moon are not sufficient to prevent a global ocean from freezing (Tobie et al., 2008; Bêhounková et al., 2017) unless Enceladus has an unconsolidated rocky core (Choblet et al., 2017).

By analogy with Earth, where most tidal energy dissipation occurs in the ocean, it has been suggested that tidal dissipation within the ocean plays a major role. Tyler (2008) was the first to propose that oceanic tidal currents could heat the moons of the outer solar system. Tyler (2008) considered the response of an ice-free shallow ocean of constant density. He solved the Laplace Tidal Equations (LTE) using the method of Longuet-Higgins (1968) for the different tidal constituents. He showed that Europa's obliquity excites planetary Rossby waves of considerable amplitude and suggested that this mechanism might close Enceladus' thermal energy budget (Tyler, 2009). However, it was later found that the obliquity of Enceladus is too small for this mechanism to generate enough tidal heating (Chen and Nimmo, 2011). Tyler (2014) and Hay and Matsuyama (2017) also showed that
the eccentricity tide can generate sufficient tidal heating but only for unrealistically shallow oceans. Subsequent studies considered the effect that the ice shell has on oceanic tidal dissipation and concluded that the ice crust further dampens the ocean response (Beuthe, 2016; Matsuyama et al., 2018).

A possible criticism to the previous studies is the assumption of an ocean of constant density. An unstratified ocean filters out internal gravity waves which have an important role in tidal dissipation and mixing in Earth's ocean (Garrett, 2003). There is little information about the stratification of Europa's and Enceladus' oceans. The ocean is heated from the bottom by tidal and radiogenic heat within the silicate mantle, likely creating a well-mixed convecting ocean (Goodman et al., 2004; Goodman and Lenferink, 2012; Soderlund et al., 2013). Nevertheless, under certain conditions the ocean might be stratified (Melosh et al., 2004; Vance and Brown, 2005). We do not challenge the assumption of an unstratified ocean and focus on the other assumptions underlying the LTE.

The main assumption of the LTE is that the ocean is in hydrostatic equilibrium, currents are predominantly horizontal. Because vertical currents are assumed to be small, the hydrostatic approximation is often used together with the traditional approximation, which consists in neglecting the terms of the Coriolis force linked to vertical motions in the ocean (Gerkema et al., 2008). For an unstratified ocean these assumptions hold as long as the ratio of the characteristic vertical and horizontal length scales is small (see Vallis (2006)). Using the ocean thickness and body's radius as a measure of vertical and horizontal length scales, respectively, higher ratios are obtained for Europa and Enceladus than for Earth (~0.06 and ~0.15 versus ~0.001). It is therefore expected that the neglected vertical velocity is of more relevance in these bodies, making the LTE incomplete to describe tidal currents in Europa's and Enceladus' subsurface oceans.

Without the hydrostatic and traditional approximations, new kinds of waves are possible in the unstratified oceans of the icy moons, the so-called internal inertial (gyroscopic) waves (Stern, 1963; Bretherton, 1964; Greenspan., 1969; Stewartson, 1971; Maas, 2001). These waves have properties that are markedly different from those of the more familiar surface waves (Maas, 2005) and have been suggested to be of importance for tidal dissipation in giant planets and binary stars (Rieutord and Valdettaro, 1997; Rieutord et al., 2001; Ogilvie and Lin, 2004; Rieutord and Valdettaro, 2010; Lainey et al., 2017). In this study, we want to take the young field of "planetary oceanography" one step further by exploring the relevance of inertial waves for tidal dissipation in the icy moons.

We consider an ocean contained within a deformable spherical shell and study tidally induced inertial waves for different ocean thicknesses. Our main goal is to quantify the amount of tidal heating that is generated by these waves to assess whether they are a significant component in the thermal energy budget of these moons. Additionally, we compare the flow amplitude of inertial waves to that of surface gravity waves obtained with the LTE and consider the footprint that they might have on the satellites' surface. To do so we solve the linearised incompressible Navier-Stokes equations for the different components of the eccentricity and obliquity tide using the spectral methods developed to study inertial waves in an astrophysical context (stars and giant planets) (e.g., Rieutord and Valdettaro, 1997; Ogilvie and Lin, 2004; Rieutord and Valdettaro, 2010).
The text is structured as follows: Section 2 introduces inertial waves and presents their main properties, in Section 3 we give the mathematical formulation of the problem, we introduce the tidal potential and then explain how the linearised Navier-Stokes equations are forced with this tidal forcing. Afterwards, the results are shown and discussed (Section 4) and conclusions are presented (Section 5)

## 2. Properties of Internal Inertial Waves

To understand the properties of internal inertial waves we start by considering a simplified situation. We consider an inviscid fluid of constant density ($\rho_0$) inside a container rotating with angular velocity $\Omega$. For this situation, the mass and momentum conservation equations, written in a co-rotating frame, are given by (e.g., Le Bars et al., 2015):

$$ \nabla \cdot \mathbf{u} = 0, \qquad (1a) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} + 2\Omega \times \mathbf{u} = -\nabla W. \qquad (1b) $$

$\mathbf{u}$ is the fluid's velocity, $W$ is a reduced pressure that contains the fluid pressure, the body's gravity potential and other possible conservative body forces. We neglect non-linear terms and seek plane wave solutions of the form:

$$ \mathbf{u} = Re \left[ \tilde{\mathbf{u}} e^{i(k \cdot x + \omega t)} \right], \qquad (2) $$

with $\mathbf{k}$ the wavevector, $\tilde{\mathbf{u}}$ the wave amplitude and $\omega$ the wave frequency. Introducing this trial solution to Equation (1) we obtain a dispersion relation of the form:

$$ \omega^2 = 4\Omega^2 \left( \frac{k_z^2}{k_x^2 + k_y^2 + k_z^2} \right), \qquad (3) $$

where $z$ is parallel to the rotational axis and $x, y$ normal to it. This dispersion relation is markedly different from that of the more familiar surface-gravity waves. While for surface waves the frequency of the wave only depends on the magnitude of the wavenumber, for internal inertial waves it only depends on the angle that the wavevector forms with the rotational axis. The group velocity, ($c_g = \frac{\partial \omega}{\partial \mathbf{k}}$) is perpendicular to the wavevector $\mathbf{k}$. Energy propagates along surfaces of constant slope which form an angle $\theta$ with the rotational axis:

$$ \theta = \arcsin \left( \frac{\omega}{2\Omega} \right), \qquad (4) $$

the so-called characteristics (see Figure 1a).

These properties lead to a fundamental difference in how inertial waves behave upon reflection as compared to surface waves (e.g., Maas, 2005). When a monochromatic surface wave packet encounters a surface it reflects specularly without changing its wavelength. In contrast, an internal inertial wave packet reflects keeping the wavevector inclination with respect to the rotational axis constant. Upon reflection, the wavelength can change.
<figure><img src="image_3.png" /><figcaption>Figure 1: Inertial wave propagation. The wave packet propagates with group velocity $c_g$ following characteristics that form an angle $\theta$ with the rotational axis. The wavevector $k$ is perpendicular to the group velocity(a). (b) Depending on the container's geometry a wave packet (blue) can converge towards a periodic trajectory (black).</figcaption></figure>

Depending on the container's geometry this can lead to focusing of the wave packet which becomes an efficient mechanism to transport energy from large to small scales.

The focusing properties of inertial waves can result in wave attractors. As its name indicates, wave attractors are trajectories in which energy accumulates. For ocean geometries that have focusing trajectories, two wave packets originating at different locations will converge towards the same trajectory (Figure 1b). Wave attractors have the peculiarity of focusing the energy of a large-scale forcing to smaller scales regardless of the nature and location of the excitation. This phenomenon has been observed in the laboratory, both for internal gravity waves (e.g., Maas et al., 1997; Brouzet et al., 2016) and internal inertial waves (e.g., Maas, 2001; Manders and Maas, 2003, 2004).

As it is of special relevance for astrophysical and geophysical applications, the propagation of inertial waves in spherical containers has been widely studied (e.g., Bryan, 1889; Stern, 1963; Bretherton, 1964; Greenspan, 1969; Stewartson, 1971, 1972; Rieutord and Valdettaro, 1997; Ogilvie, 2009). The response of an inviscid fluid inside an spherical container is given by well-behaved eigenmodes (Bryan, 1889). However when a nucleus is added, wave attractors can develop which lead to singularities; velocity increases without bound along the attractor (Stewartson, 1971). Something similar happens at the critical latitude ($\Theta_c$), the latitude at which a characteristic is tangent to the nucleus, where the velocity also develops a singularity (Rieutord et al., 2001). In reality, viscosity prevents the development of such singularities. The result is the development of prominent internal shear layers that propagate in the fluid domain following the characteristics where a significant amount of energy can be dissipated.

### 3. Problem Formulation

#### 3.1. The tidal potential

The obvious candidate for the generation of inertial waves in the oceans of the icy moons is the tidal potential caused by the obliquity and eccentricity of the moons. We can express the
tidal potential in terms of the orbital elements of the satellite (Jara-Orué and Vermeersen, 2011):

$$ \Psi = (\omega R)^2 \left(\frac{r}{R}\right)^2 \{\psi_0 + \psi_{ns} + \psi_{e0} + \psi_{e2} + \psi_{o1}\}, \quad (5) $$

where $\omega$ is the mean motion of the satellite given by $\omega^2 = GM_p/a^3$, $G$ is the universal gravitational constant, $M_p$ the mass of the planet and $a$ the semi-major axis of the satellite's orbit. As we are considering tidally locked satellites, $\omega$ is the same as the rotational frequency of the satellite ($\Omega$). $R$ is the satellite's radius and $r$ the distance from the center of the satellite to the considered point inside the satellite. The different components of the tidal potential are (Jara-Orué and Vermeersen, 2011):

$$ \psi_0 = -\frac{1}{2}P_{2,0}(\cos\theta) + \frac{1}{4}P_{2,2}(\cos\theta)\cos(2\phi), \quad (6a) $$

$$ \psi_{ns} = -\frac{1}{2}P_{2,2}(\cos\theta)\sin(2\phi + \Omega_{nst}\theta)\sin(\Omega_{nst}\theta), \quad (6b) $$

$$ \psi_{e0} = -\frac{3e}{2}P_{2,0}(\cos\theta)\cos(\omega t), \quad (6c) $$

$$ \psi_{e2} = \frac{e}{4}P_{2,2}(\cos\theta)[3\cos(2\phi)\cos(\omega t) + 4\sin(2\phi)\sin(\omega t)], \quad (6d) $$

$$ \psi_{o1} = P_{2,1}(\cos\theta)\sin(\epsilon)\cos(\phi)\sin(\Phi + \omega t). \quad (6e) $$

$\theta$ and $\phi$ are the co-latitude and longitude in the body-fixed reference frame, respectively. $\psi_0$ is the static component of the tidal potential, $\psi_{ns}$ arises due to non-synchronous rotation. The eccentricity tide is given by $\psi_{e0}$, which is caused by the variation in distance between the planet and the satellite; and $\psi_{e2}$, caused by the east/west libration of the position of the subplanet point on the moon's surface. On the other hand, the obliquity tide, $\psi_{o1}$, is the result of the latitudinal libration of the subplanet point due to the satellite's obliquity. $P_{l,m}$ are the associated Legendre polynomials of degree $l$ and order $m$ and $\Phi$, $e$ and $\epsilon$ are the argument of pericenter, the eccentricity and the obliquity of the moon, respectively. We only consider degree two terms of the tidal potential, as the amplitude of the other components rapidly decreases with increasing degree as $(r/a)^l$.

The non-axisymmetric term of the eccentricity tide ($\psi_{e2}$) and the obliquity tide ($\psi_{o1}$) can be further divided in a westward ($\psi_{e2w}$, $\psi_{o1w}$) and an eastward ($\psi_{e2e}$, $\psi_{o1e}$) propagating wave:

$$ \psi_{e2} = \psi_{e2e} + \psi_{e2w} = \frac{7e}{8}P_{2,2}(\cos\theta)\cos(2\phi - \omega t) - \frac{e}{8}P_{2,2}(\cos\theta)\cos(2\phi + \omega t), \quad (7) $$

$$ \psi_{o1} = \psi_{o1e} + \psi_{o1w} = \frac{\sin(\epsilon)}{2}P_{2,1}(\cos\theta)\sin(\phi - \omega t - \Phi) - \frac{\sin(\epsilon)}{2}P_{2,1}(\cos\theta)\sin(\phi + \omega t + \Phi). \quad (8) $$
<table><caption>Table 1: Physical parameters for Europa and Enceladus. Radius (R), mass (M), average density ($\rho_{av}$), surface gravity (g) and rotational rate (Ω) are obtained from: https://ssd.jpl.nasa.gov/?horizons and Chen et al. (2014). The obliquity (ε) of Europa and Enceladus are obtained assuming that the moons are in a Cassini state (Baland et al., 2012; Chen and Nimmo, 2011). The maximum amplitude of the different components of the equilibrium tide are computed considering a non-self gravitating ice-free ocean surround- ing an infinitely rigid core (Equation (17)) and using the maximum value of the satellites’ obliquity. The argument of pericenter (Φ) and ocean density ($\rho_o$) are assumed to be 0 and 1000 kg m⁻³, respectively.</caption><thead><tr><th></th><th>Europa</th><th>Enceladus</th></tr></thead><tbody><tr><td>R [km]</td><td>1565.0</td><td>252.3</td></tr><tr><td>M [kg]</td><td>4.797 &middot; 10<sup>22</sup></td><td>1.0805 &middot; 10<sup>20</sup></td></tr><tr><td>&rho;<sub>av</sub> [kg m<sup>-3</sup>]</td><td>2990</td><td>1606</td></tr><tr><td>g [m s<sup>-2</sup>]</td><td>1.31</td><td>0.11</td></tr><tr><td>&Omega; [rad s<sup>-1</sup>]</td><td>2.05 &middot; 10<sup>-5</sup></td><td>5.31 &middot; 10<sup>-5</sup></td></tr><tr><td>e [-]</td><td>0.0094</td><td>0.0047</td></tr><tr><td>&epsilon; [deg]</td><td>&lt; 0.1</td><td>&lt; 0.0015</td></tr><tr><td>z<sub>e0</sub><sup>(eq)</sup> [m]</td><td>11.08</td><td>11.50</td></tr><tr><td>z<sub>e2e</sub><sup>(eq)</sup> [m]</td><td>19.39</td><td>20.13</td></tr><tr><td>z<sub>e2w</sub><sup>(eq)</sup> [m]</td><td>2.77</td><td>2.88</td></tr><tr><td>z<sub>o2e</sub><sup>(eq)</sup> [m]</td><td>1.03</td><td>0.003</td></tr><tr><td>z<sub>o2w</sub><sup>(eq)</sup> [m]</td><td>1.03</td><td>0.003</td></tr><tr><td>&Phi; [deg]</td><td>0</td><td>0</td></tr><tr><td>&rho;<sub>o</sub> [kg m<sup>-3</sup>]</td><td>1000</td><td>1000</td></tr></tbody></table>

For this study we consider the different components of the eccentricity and obliquity tide and ignore the non-synchronous rotation of the satellite. We also ignore the static component of the tidal potential as it does not induce a dynamic ocean response. The relevant physical parameters of Europa and Enceladus are given in Table 1.

## 3.2. Governing equations and assumptions

In contrast to Section 2, we now consider the fluid to be viscous and expand the reduced pressure in its different components. The equations of motion can be written as:

$$ \nabla \cdot \mathbf{u} = 0, \qquad (9a) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} + 2\Omega \times \mathbf{u} = -\frac{1}{\rho_0} \nabla p' + \frac{1}{\rho_0} \nabla \cdot \boldsymbol{\tau} - \nabla \Psi - \nabla \Psi'. \qquad (9b) $$

The primed quantities denote deviations from hydrostatic equilibrium. $\Psi'$ represents the gravitational potential of the body, $p'$ is the fluid pressure and $\tau$ the stress tensor of the fluid, which we assume to follow Stokes’ constitutive law:

$$ \boldsymbol{\tau} = \mu(\nabla \mathbf{u} + \nabla \mathbf{u}^T), \qquad (10) $$

where $\mu$ is the dynamic molecular viscosity. We further use perturbation theory and consider the different quantities to be small. Under such consideration the non-linear advective term
(u · ∇u) is neglected. The validity of this assumption will be discussed later (Section 5). We introduce an effective viscosity μeff which is higher than the molecular viscosity of the fluid and accounts for non-modelled small-scale dissipation mechanisms (e.g., turbulent mixing, interactions of waves with turbulent convection, etc.).

Finally, we obtain the non-dimensional form of the equations of motion by using the inverse of the satellite's rotational frequency (2Ω)⁻¹ and the satellite radius (R) as time and length scales respectively:

$$ \hat{\nabla} \cdot \hat{\mathbf{u}} = 0, \qquad (11a) $$

$$ \frac{\partial \hat{\mathbf{u}}}{\partial \hat{t}} + \mathbf{e}_z \times \hat{\mathbf{u}} = -\hat{\nabla} \hat{W} + E \hat{\nabla}^2 \hat{\mathbf{u}}, \qquad (11b) $$

with e_z the direction of the satellite's rotational axis, W again the reduced pressure:

$$ W = \frac{p'}{\rho_o} + \Psi + \Psi', \qquad (12) $$

and the non-dimensional variables (indicated with a hat) and parameters defined as:

$$ \mathbf{u} = 2\Omega R \hat{\mathbf{u}}, \qquad (13a) $$

$$ W = 4\Omega^2 R^2 \hat{W}, \qquad (13b) $$

$$ E = \frac{\mu_{eff}}{2\rho\Omega R^2}, \qquad (14a) $$

$$ \hat{\omega} = \frac{\omega}{2\Omega}. \qquad (14b) $$

We have introduced the non-dimensional Ekman number (E), which gives the ratio of viscous to Coriolis forces; and the non-dimensional frequency (ω̂), which, since the problem is forced at the diurnal tidal frequency, equals 0.5. In the discussion that follows we use the non-dimensional equations; to avoid cumbersome notation we drop the hat from the variables.

As a starting point we assume that the moons have a free surface. If the surface wave speed (√gh, with g the satellite's surface gravity and h the ocean depth) is high enough so that the ocean adjusts quickly to forces varying at the tidal potential frequency, we can assume that the radial displacement of the ocean surface is given by the equilibrium tide (Tyler, 2008; Ogilvie, 2009). This assumption breaks down if the ocean is too shallow, in that case surface gravity waves dominate the ocean response and other kinds of resonances occur (Matsuyama et al., 2018). For the icy moons, the surface wave speed is high enough as long as the oceans are thicker than 0.78 km or 1.6 km for Europa and Enceladus, respectively. With ocean thicknesses of ~100 km (Hussmann et al., 2002) and ~38 km (Beuthe et al., 2016) for Europa and Enceladus respectively, the previous assumption seems reasonable for the problem at hand.
202

203

204

We follow the approach of Ogilvie (2005) and Rieutord and Valdettaro (2010) and split the fluid response into two components: the equilibrium tide $\mathbf{u}^{(\text{eq})}$ plus a dynamical tide $\mathbf{u}^{(\text{d})}$. By using this decomposition Equation (11) becomes:

$$ \frac{\partial \mathbf{u}^{(\text{d})}}{\partial t} + \mathbf{e}_z \times \mathbf{u}^{(\text{d})} = -\nabla W + E \nabla^2 \mathbf{u}^{(\text{d})} + \mathbf{f}^{(\text{eq})}. \quad (15) $$

205

$\mathbf{f}^{(\text{eq})}$ is an inertial force associated with the equilibrium tide given by:

$$ \mathbf{f}^{(\text{eq})} = -\frac{\partial \mathbf{u}^{(\text{eq})}}{\partial t} - \mathbf{e}_z \times \mathbf{u}^{(\text{eq})}. \quad (16) $$

206

207

208

209

Note that the dynamic tide is forced indirectly through the equilibrium tide as it does not satisfy the momentum equation (Ogilvie and Lin, 2004). If we consider a non self-gravitating ocean around an infinitely rigid core with a free-surface the equilibrium tide is simply given by:

$$ \mathbf{z}_{l,m}^{(\text{eq})} = \frac{\Psi_{l,m}}{g} \mathbf{e}_r, \quad (17) $$

210

211

212

213

where $\Psi_{l,m}$ is the degree $l$ order $m$ component of the tidal potential. The maximum amplitude of the equilibrium tide for Europa and Enceladus for the different tidal components is given in Table 1. If the effect of self-gravitation, the finite rigidity of the core or the presence of an ice-shell is considered the equilibrium tide can be obtained as:

$$ \mathbf{z}_{l,m}^{(\text{eq})} = Z_{l,m} \frac{\Psi_{l,m}}{g} \mathbf{e}_r, \quad (18) $$

214

215

216

217

218

219

220

with $Z_{l,m}$ being an admittance factor. Matsuyama et al. (2018) showed that an ice shell dampens the response of the ocean resulting in $Z_{l,m}$ to be smaller than or close to 1. In particular, for realistic ice shell thicknesses, the admittance is close to 1 for Europa while it can be more than one order of magnitude smaller for Enceladus. In this work, we assume the ocean to follow the equilibrium tide of an ice-free, non self-gravitating ocean in a satellite with an infinitely rigid core (Equation (17)). By doing so, we obtain an upper bound of tidal dissipation due to inertial waves in a spherical shell shaped ocean.

221

222

223

We assume the system to have a response with the same frequency as the forcing tidal potential and thus we consider the different fields to be proportional to $e^{i\omega t}$. By doing so we can eliminate time from the equations of motion and we are left with the system:

$$ \nabla \cdot \mathbf{u}^{(\text{d})} = 0, \quad (19a) $$

224

$$ i\omega\mathbf{u}^{(\text{d})} + \mathbf{e}_z \times \mathbf{u}^{(\text{d})} = -\nabla W + E\nabla^2\mathbf{u}^{(\text{d})} + \mathbf{f}^{(\text{eq})}, \quad (19b) $$

225

$$ \mathbf{f}^{(\text{eq})} = -i\omega\mathbf{u}^{(\text{eq})} - \mathbf{e}_z \times \mathbf{u}^{(\text{eq})}. \quad (19c) $$

226

227

228

We solve Equations (19) in a spherical shell, assuming no-slip boundary conditions at the fluid-solid interfaces ($\mathbf{u}_d = 0$).
## 3.3. Numerical Approach

We solve Equations (19) for the different constituents of the equilibrium tide (Equation (17)). As we did with the tidal potential, we expand the equilibrium tide into its different constituents: the order 0 eccentricity tide ($u_{e0}^{(eq)}$), and the eastward and westward components of the order 2 eccentricity tide and order 1 obliquity tide ($u_{e2e}^{(eq)}$, $u_{e2w}^{(eq)}$ and $u_{ole}^{(eq)}$, $u_{olw}^{(eq)}$) (see Appendix A). We then solve Equations (19) using $\omega = 0.5$ and $\omega = -0.5$ for the westward and eastward propagating components respectively.

We use the method of Rieutord and Valdettaro (1997, 2010), which is detailed in Appendix A. We make use of the spherical symmetry of the problem and expand the velocity and reduced pressure fields using L spherical harmonics in the horizontal direction. The resulting equations are discretised in the radial direction using Chebyshev polynomials on $N + 1$ Gauss-Lobatto collocation nodes. By doing so, Equation (19) and the boundary conditions result in an algebraic system of $(L - m + 1)(N + 1)$ linear equations. The associated matrix of this linear system is block-tridiagonal, and the system is solved using classical LU-factorization of a banded matrix. Except where otherwise indicated, the resolution, given by L and N, is chosen so that the truncation error is less than $\sim 10^{-4}$.

We are specially interested in computing tidal dissipation due to the tidally-induced flows. We compute the amplitude of tidal dissipation as (e.g., Ogilvie and Lin, 2004):

$$ \hat{D}_v = 2E \int_V c_{ij}^* c_{ij} dV, \qquad (20) $$
where $c_{ij}$ are the elements of the rate-of-strain tensor, $c_{ij} = \frac{1}{2}(\partial_i u_j + \partial_j u_i)$. We expand Equation (20) using spherical harmonics as indicated in Appendix B.

# 4. Results

## 4.1. Parameter Regime

The two parameters that control the fluid response are the ocean thickness and the Ekman number. We characterise the thickness of the ocean using the ratio of the inner to outer ocean radius, $\eta$. Europa's and Enceladus' ocean thicknesses are still not known. To understand the effect of ocean geometry on the propagation of inertial waves, we explore the range $\eta = 0.3-0.99$, which corresponds to ocean thicknesses ranging from 1095.5 to 15.7 km for Europa and 176.6 to 2.5 km for Enceladus. For this range we use a resolution of $\Delta\eta = 0.005$ which is equivalent to ocean thickness changes of 7825 m for Europa and 1262 m for Enceladus. However, current estimates suggest an $\eta$ of $\sim 0.93$ (Hussmann et al., 2002) for Europa and $\sim 0.85$ (Beuthe et al., 2016) for Enceladus. It is for this reason that we further explore the range $\eta = 0.8-0.99$ using a $\Delta\eta$ of 0.0005.

The Ekman number depends on the fluid viscosity (Equation (14)). If we compute the Ekman number using the molecular viscosity of water we obtain a value of $10^{-14}$ and $10^{-13}$ for Europa and Enceladus respectively. These low Ekman numbers require currently unattainable resolutions. In any case, when linearising the equations of motion we have introduced a new effective dynamic viscosity ($\mu_{eff}$) that accounts for small-scale dissipation
mechanisms. Given our limited knowledge about these oceans it is difficult to estimate the appropriate value of this parameter. It is however expected that the effective viscosity will be orders of magnitudes higher than the molecular viscosity (e.g., Ogilvie and Lin, 2004). In the following experiments we explore a wide range of Ekman numbers ($E = 10^{-4} - 10^{-10}$) to account for this uncertainty.

We first study the propagation of inertial waves for different ocean thicknesses and Ekman numbers and explain the main characteristics in terms of the properties defined in Section 2 such as wave attractors and the critical latitude singularity (Section 4.2). Afterwards, we focus on tidal dissipation and quantify the relevance of inertial waves for Europa's and Enceladus' thermal energy budget (Section 4.3).

## 4.2. Wave attractors and the critical latitude singularity

As explained in Section 2 the response of the fluid is dominated by the convergence of wave packets towards wave attractors and the critical latitude singularity. We start by studying the structure of wave attractors that are excited by tidal forcing. As the moons are phase-locked, the frequency of the forcing equals the rotational frequency of the moon, hence, the angle that the characteristics form with the vertical is 30° (Figure 1a). We use ray-tracing to study the propagation of energy in meridional planes of the fluid domain. The intersection of the characteristic surfaces with a meridional plane are straight lines (rays) along which energy propagates. We launch a wave packet from a point on the inner sphere without zonal velocity and follow its propagation along the characteristics until it converges towards an attractor in the same meridional plane. Note that when a wave packet is launched with a zonal component it may end up trapped in a meridional plane or escape meridional trapping and reflect endlessly around the domain (Rabitti and Maas, 2013, 2014). However, because of the symmetry of the tidal potential with respect to rotation around the z axis we focus on wave attractors in the meridional plane.

We launch a wave ray from a point equatorward of the critical latitude (5°) and one poleward (85°). The latitude of the last 100 inner boundary reflections are then noted and shown in a Poincaré plot in Figure 2a. Additionally, we compute the Lyapunov exponent, which measures the asymptotic rate at which two neighbouring rays converge:

$$ \Lambda = \lim_{N \to \infty} \sum_{n=1}^{N} \log \left| \frac{d\Theta_{n+1}}{d\Theta_n} \right| \quad (21) $$

$d\Theta$ is the angular distance between the reflection points of two neighbouring characteristics. High (in absolute value) Lyapunov exponents mean highly attracting wave attractors while a Lyapunov exponent equal to 0 indicates non-attracting trajectories. The Lyapunov exponent for different ocean geometries is displayed in Figure 2b.

We distinguish between two types of attractors: equatorial attractors, trapped equatorward of the critical latitude, and polar attractors, with reflections outside the previously mentioned interval. As can be seen in Figure 2a, the shape of the attractor depends on the ratio of inner to outer ocean radius ($\eta$). For some ocean geometries, attractors with few reflections (short attractors) and high (in absolute value) Lyapunov exponents exist. This is
the case of the band [0.5775, 0.6545] and [0.7235, 0.7420]. In the first case, both equatorial and polar attractors exist, while for the second only polar attractors appear. The peak for $\eta = 0.6545$ corresponds to an attractor with reflections infinitesimally close to the critical latitude, where there is nearly infinite focusing. For such attractors, the Lyapunov exponent goes to minus infinity as resolution in $\eta$ increases (Rieutord et al., 2001).

For $\eta < 0.57$, there are non-attracting periodic trajectories with a Lyapunov exponent close to 0. In fact, it can be shown that for $\eta \le 0.5$ all trajectories associated with the frequency $\omega/2\Omega = 0.5$ are strictly periodic and have a Lyapunov exponent equal to 0 (Rieutord et al., 2001). For these ocean geometries if a ray is launched from the inner sphere it will eventually return to the same point after some reflections, explaining the horizontal lines in Figure 2a.

The response of the fluid is characterised by the opposing effects of wave focusing and viscous diffusion. The focusing effect of the geometry tends to focus energy towards small scales ultimately leading to a singularity along the wave attractor. On the other hand, viscosity diffuses the velocity countering the geometrical focusing effect and the development of wave attractors. As we will now see, these two competing effects determine whether a wave attractor is excited for a particular ocean geometry or not.

We start by choosing two ocean geometries with a high (in absolute value) Lyapunov exponent, $\eta = 0.63$ and $\eta = 0.73$, and force the fluid with the axisymmetric eccentricity tidal component for two different Ekman numbers, $10^{-7}$ and $10^{-9}$. To achieve a truncation error less than $\sim 10^{-4}$, we use a resolution of $L = 500$ and $N = 200$ for $E = 10^{-7}$ and $N = 200$ and $L = 1200$ for $E = 10^{-9}$. A meridional cut of the sphere showing the amplitude of the kinetic energy and viscous dissipation is shown in Figure 3. For this tidal component the plots are independent of the meridional cut. Also superimposed, we show the wave attractors characteristic of these two ocean geometries.

As can be seen, in both cases wave attractors are generated. Internal shear layers are clearly observed along the attractors where both kinetic energy and viscous dissipation are enhanced. As in Rieutord and Valdettaro (2010) we observe the importance of the critical latitude singularity for producing the attractor. The attractor seems to be “fed” by the shear layer emanating from its critical latitude. As expected, as the Ekman decreases the internal shear layer becomes thinner.

We also explore the behaviour of the fluid for ocean geometries with a low (in absolute value) Lyapunov exponent. In such cases, the critical latitude singularity dominates the flow (Ogilvie and Lin, 2004; Ogilvie, 2009; Rieutord and Valdettaro, 2010). We distinguish two cases: ocean geometries where periodic non-focusing trajectories exist (e.g., $\eta = 0.35, 0.5$), $\Lambda = 0$, and geometries with long weakly-attracting wave attractors (e.g., $\eta = 0.75, 0.8815$), $\Lambda < 0$. For the first case, the shear layer emitted from the critical latitude propagates following the characteristics and widens slowly due to viscous dissipation and the lack of focusing (Figures 4a-4b). In contrast, for the second case we observe the focusing of the shear layer upon reflection (Figures 4c-4e).

For $\eta = 0.75$ and $0.8815$ we observe that as the ray approaches the polar regions both kinetic energy and dissipation are enhanced (Rieutord and Valdettaro, 1997). The effect of changing the Ekman number is similar to that reported for those cases where wave attractors
<figure><img src="image_4.png" /><figcaption>Figure 2: Poincaré plot (a) and Lyapunov exponent (b) for a forcing frequency of $\omega/(2\Omega) = 0.5$. In (a) the latitude of the last 100 reflection points on the inner sphere for rays launched from $\pm 5^\circ$ and $\pm 85^\circ$ latitude are indicated for different inner to outer ocean radius ratios ($\eta$). The shaded region indicates reflections equatorward of the critical latitude. Two specific ocean geometries for which short, highly attracting-wave attractors exist ($\eta = 0.63, 0.73$) are indicated in red and shown in (c) and (d), respectively. (b) gives the Lyapunov exponent for the polar attractors.</figcaption></figure>
347 appear. It is expected that if the Ekman number is sufficiently small the focusing effect will
348 prevail over viscous diffusion leading to the generation of long wave attractors for $\eta = 0.75$
349 and 0.8815 similar to those observed for $\eta = 0.63$ and 0.73. However, the low Ekman number
350 for which this transition occurs requires resolutions currently computationally unattainable.

351 The case of $\eta = 0.8815$ is especially interesting as it corresponds to an Enceladan ocean
352 thickness of around 30 km, which is close to current estimates for Enceladus (Thomas et al.,
353 2016; Čadek et al., 2016; Beuthe et al., 2016). For this ocean thickness we compute the
354 maximum velocity amplitude of inertial tides and study the possible manifestation of inertial
355 waves on Enceladus' surface. We find that the maximum velocity amplitude varies from 0.5
356 to 3 cm/s for Ekman numbers ranging from $10^{-7}$ to $10^{-10}$. Moreover, the maximum velocity
357 amplitude is always encountered in polar regions. These currents are stronger than those
358 induced by the eccentricity tide under the shallow water approximation, which for realistic
359 ocean thicknesses have a magnitude of around 1 mm/s (Tyler, 2009).

360 Vermeersen et al. (2013) suggested that wave attractors in a polar ocean basin could
361 be the origin of Enceladus’ tiger-stripes; here we investigate this hypothesis for the present
362 configuration (an unstratified global ocean of constant depth). Interaction between the fluid
363 and the ice shell are likely to occur at the inertial waves reflection points, where we observe
364 enhanced dissipation. For a 30 km thick ocean, we observe that the pressure at these points
365 is in the order of $\sim 1$ Pa. Moreover, we can estimate the melting rate of the ice shell at
366 these points by assuming that all the energy dissipated in a ray is focused there and does
367 not radiate outwards. We obtain a maximum melting rate in the order of 1 m every 10
368 kyr. Even though it is unclear how inertial waves would interact with the ice shell, these
369 numbers suggest that in an ocean of constant thickness it is unlikely that tidally-excited
370 inertial waves could be the origin of observed surface features, such as the tiger stripes.
371 More work is needed to study if other ocean geometries can result in enhanced energy
372 focusing and higher stresses.

### 4.3. Tidal dissipation

373 The main focus of this work is to study tidal dissipation within Europa and Enceladus. We
374 start by assuming an Ekman number of $10^{-7}$ and compute tidal dissipation for the different
375 tidal constituents and different ocean thickness ($\eta = 0.3 - 1$) as detailed in Appendix B
376 (Section 4.3.1). To strengthen our conclusions we analyse the effect of varying the Ekman
377 number on the degree-two order-zero tidal constituent (Section 4.3.2). For the different
378 shell geometries we first vary the Ekman number between $10^{-4} - 10^{-8}$. We then study the
379 asymptotic behaviour of dissipation with decreasing Ekman number for some relevant shell
380 geometries for which we use an Ekman number as low as $10^{-10}$.
381

#### 4.3.1. Tidal dissipation for the different tidal components

382 Dissipation due to tidally induced inertial waves in Europa and Enceladus for $E = 10^{-7}$,
383 is shown in Figure 5. As a reference, we also indicate the estimated value of radiogenic
384 heating in Europa (Hussmann et al., 2010) and the observed thermal output of Enceladus’
385 tiger stripes (Spencer et al., 2013).
386
<figure><img src="image_5.png" /><figcaption>Figure 3: Non-dimensional kinetic energy (left quadrant) and viscous dissipation (right quadrant) amplitude due to the degree two, order zero eccentricity tide for two ocean geometries where short wave attractors exist. The patterns are shown for different Ekman number (E) and inner to outer ocean radius ratio (η). Polar and equatorial wave attractors are superimposed in black and green, respectively, and the critical latitude is indicated in red. A logarithmic scale is used both for kinetic energy and viscous dissipation. The maximum values of the non-dimensional colour-scale correspond to a kinetic energy of 0.65 J m⁻³ and 4.69 J m⁻³; and a viscous dissipation of 0.26 µW m⁻³ and 4.97 µW m⁻³ for Europa and Enceladus, respectively.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Same as Figure 3 but for ocean geometries with a small Lyapunov exponent. Rays emanating from the critical latitude (and not wave attractors) dominate the flow response.</figcaption></figure>
It can be seen that the values of the tidal dissipation are well below the reference values for both Europa and Enceladus. Moreover, when we consider the region of the plot that corresponds to likely values of ocean thickness on Europa and Enceladus, $\eta > 0.9$ (Anderson et al., 1998; Hussmann et al., 2002) and $\eta > 0.8$ (Thomas et al., 2016; Čadek et al., 2016; Beuthe et al., 2016), respectively, we observe that the values of tidal dissipation are several orders of magnitude smaller than radiogenic heating in Europa and the observed heat flux in Enceladus. Nevertheless, it is interesting to study the contribution of the different tidal constituents. We see that the axisymmetric eccentricity and eastward tide dominate the fluid response. Moreover, the order one obliquity tide produces considerably lower tidal dissipation, especially in Enceladus. This fact follows from the small amplitude of this tidal component as compared with the others (see Table 1).

For most tidal constituents the tidal dissipation curve is markedly spiky (Figure 5). There are values of $\eta$ for which dissipation is enhanced by more than one order of magnitude. The complex system of internal-shear layers that arise due to the singularities of the problem results in different values of dissipation for different ocean geometries. This is not the case of the westward propagating obliquity tide whose dissipation is given by a much smoother curve. This is because for a degree two, order one forcing the equations of motion admit a purely toroidal solution (Rieutord and Valdettaro, 1997), the so called “spin-over” mode. This solution is an exact solution of the equations of motion in case stress-free boundary conditions are used.

However, when no-slip boundary conditions are used for the westward obliquity tide the “spin-over” mode disappears. The use of no-slip boundary conditions results in the development of an oscillatory Ekman layer at the solid-liquid boundary. The thickness of this Ekman layer scales with $\approx E^{1/2}$, however at the critical latitude it thickens and scales as $\sim E^{2/5}$ over a region of width $\sim E^{1/5}$ (Hollerbach and Kerswell, 1995; Kerswell, 1995). The thickening of the Ekman layer breaks the symmetry of the problem and launches inertial waves that propagates in the interior. However, the resulting internal shear layers are weaker in this case than those generated by the other tidal constituents, dissipation is dominated by the Ekman layer that forms in the solid-liquid boundary which is less dependent on the shell geometry.

### 4.3.2. Dependence of tidal dissipation on Ekman number

We study the dependence of tidal dissipation with Ekman number. For a given ocean thickness ($\eta$) we distinguish three different cases depending on the behaviour of viscous dissipation in the limit $E \to 0$. Dissipation can increase, decrease or become asymptotically constant. We call the first and second cases resonance and anti-resonance, respectively, after Rieutord and Valdettaro (2010). As we will see, the last case is associated with the appearance of a wave-attractor Note that this definition differs from the classic definition of a resonance, commonly associated with the excitation of eigenmodes.

We start by varying the Ekman number from $10^{-4}$ to $10^{-8}$ for the axisymmetric eccentricity tide (Figure 6). For low Ekman numbers ($E = 10^{-4}, 10^{-5}$) the flow is dominated by the Ekman boundary layer formed at the solid-fluid boundary and not internal shear layers and wave attractors (see Figure 4c). For such cases the dissipation curve is smooth, similarly
<figure><img src="image_7.png" /><figcaption>Figure 5: Tidal dissipation amplitude due to tidally excited inertial waves in Europa (a) and Enceladus (b). The ocean is assumed to be within a spherical shell of constant thickness (inner to outer ocean radius ratio, $\eta$). All values are computed using $E = 10^{-7}$. The axisymmetric and non-axisymmetric components of the eccentricity tide ($m = 0$ and $m = 2$) as well as the obliquity tide ($m = 1$) are considered. As reference, an estimate of Europa's radiogenic heating and the power radiated from Enceladus tiger stripes are indicated.</figcaption></figure>
429

430

431

to the westward order one obliquity tide. As the Ekman number decreases, internal shear layers and wave attractors become the dominant features of the fluid response and the curve becomes more spiky.

432

433

434

435

436

437

438

439

440

We see that for small values of $\eta$ dissipation decreases with Ekman number, an anti-resonance. The behaviour becomes more complicated for thinner oceans for which resonances occur for some values of $\eta$ (e.g., 0.75, 0.807, 0.8415, 0.8815), moreover, we see that these cases correspond to local maxima of tidal dissipation (in $\eta$). We select some representative cases of $\eta$ in which we observe wave attractors ($\eta$ = 0.63 and 0.73), anti-resonances ($\eta$ = 0.35) or resonances ($\eta$ = 0.76 and 0.8815) for further study. For these cases we further explore the change of dissipation as the Ekman number decreases (Figure 7). We use an Ekman number as low as $10^{-10}$, which requires enhanced resolution ($L$ = 2000, $N$ = 350) to achieve a truncation error of $\sim 10^{-3}$.

441

442

443

444

445

446

447

For a thick ocean ($\eta$ = 0.35), dissipation decreases with decreasing Ekman number. Dissipation approximately follows a power law of the form $\sim E^{0.35}$. As explained, for this ocean geometry, strictly non-attracting periodic trajectories exist. Rieutord and Valdettaro (2018) have shown that in such a situation the fluid responds to the excitation with a flow characterised by an ever increasing wavenumber that ultimately, when $E=0$, inhibits any response. This is actually the anti-resonance associated with the periodic trajectory $\omega = \sin(\pi/6)$.

448

449

450

451

452

453

454

455

456

457

458

We now consider the two cases with highly attracting short wave attractors, $\eta$ = 0.63 and $\eta$ = 0.73, depicted in Figure 3. Figure 7 shows that for these two cases dissipation tends towards an asymptotic limit. In a wave attractor there is a balance between focusing and diffusion. As the Ekman number is reduced, the thickness of the internal shear layer is reduced and the velocity gradient increased in such a way that dissipation remains constant. This situation is similar to that observed in Ogilvie (2005) who showed the asymptotic convergence of dissipation in wave attractors in a rectangular container. It is also interesting to note that the asymptotic limit is reached at higher Ekman numbers for $\eta$ = 0.63 than for $\eta$ = 0.73. Figure 3 shows that the $\eta$ = 0.73 attractor is longer and less attracting, viscous dissipation acts along a longer distance and prevents the attractor from developing for high Ekman number.

459

460

461

462

463

464

465

466

467

468

469

470

Finally, we focus on two cases where a resonance occurs ($\eta$ = 0.75 and 0.8815). These resonances are of special relevance if the Ekman number of Europa's and Enceladus' ocean is very low, for example close to that given by the molecular viscosity of water ($10^{-14} - 10^{-13}$). In such a case, if a resonant state is attained, dissipation could be considerably higher than the values computed for Ekman numbers several orders of magnitude higher. As seen in Figure 7, for these two cases dissipation increases with decreasing Ekman number until it reaches a maximum and then starts to decrease again. This behaviour can be understood by analysing the resonance peaks in Figure 6b. We see that as the Ekman number decreases the curve becomes more spiky, the resonance peak narrows. If the exact $\eta$ for which the resonance occurs is not chosen, dissipation will decrease with decreasing Ekman number as the resonance peak becomes narrower. As is evident from this plot, our current resolution does not allow us to resolve the exact value of $\eta$ at which dissipation reaches its maximum.

471

Higher resolution is needed in our explored parameter space to find the exact resonant
472

473

474

475

476

477

478

479

480

481

values of $\eta$. It is likely that at these exact resonant ocean geometries the actual dissipation for an Ekman number close to that given by the molecular viscosity of water ($E = 10^{-13}-10^{-14}$) is several orders of magnitude higher than what we find in Figure 6. However, we need to recall that as the Ekman number decreases the resonance peak becomes narrower, thus a resonance would only occur over a very narrow range of ocean thicknesses. For instance, we see that by using $\Delta\eta = 126$ m for Enceladus or 783 m for Europa the resonance peak for $\eta = 0.8815$ at an Ekman number of $10^{-9}$ is already missed. It seems impossible that the ocean of Europa or Enceladus is of the thickness precisely required for such resonances to occur, specially since this would require the sea-floor and ice shell to be devoid of topographical features of higher amplitude than the previously mentioned resolution.

482

## 5. Discussion and Conclusions

483

484

485

486

487

488

489

In this work we analyse for the first time tidal dissipation due to inertial waves in Europa and Enceladus. We consider tides caused by the eccentric orbit of these satellites as well as their obliquity. We consider an ocean contained within a spherical shell and use a spectral method to compute tidal dissipation for different ocean thicknesses and dynamic viscosity coefficients (given by the Ekman number). We neglect the effect of the ice shell and argue that our results represent an upper bound estimate for tidal dissipation for an ice-covered moon.

490

491

492

493

494

495

496

497

498

499

500

We find that dissipation depends strongly on ocean thickness and Ekman number. One of the challenges in computing tidal dissipation is the scarcity of information about Europa and Enceladus oceans. For instance, the effective viscosity, which depends on small-scale dissipative processes, is not known. It is for this reason that we analyse the behaviour of the fluid for a wide range of Ekman numbers. For a given ocean thickness we distinguish three different scenarios depending on the behaviour of dissipation with Ekman number. If a wave attractor is excited, an asymptotic limit is reached where dissipation is independent of viscosity. Otherwise, there can be an anti-resonance or resonance depending on the ocean geometry. Anti-resonant states, where we observe a decrease of dissipation with decreasing Ekman number, are attained for unrealistically thick oceans; while resonant states, where dissipation increases with decreasing Ekman number, are common for shallower oceans.

501

502

503

504

505

506

507

508

509

510

After analysing the effect of changing the ocean thickness and viscosity we conclude that under the aforementioned assumptions, tidal dissipation due to inertial waves is several orders of magnitude smaller than Europa's radiogenic heating and Enceladus' observed heat flux and thus does not play an important role in preventing these oceans from freezing. Still, we observe that the induced tidal currents can be one order of magnitude stronger than those obtained using the Laplace Tidal Equations (e.g., Tyler, 2008, 2009). We find that for a 30 km thick Enceladan ocean tidal currents of amplitude 3 cm/s are excited. We also consider the possible interaction of inertial waves and the ice shell and conclude that for global oceans of constant thickness it is unlikely that inertial waves could result in observable surface features.

511

512

The difficulty in explaining Enceladus' present state through tidal dissipation in its solid parts (ice and mantle) have led to a focus on tidal dissipation within the ocean in the past
<figure><img src="image_8.png" /><figcaption>Figure 6: (a) Tidal dissipation amplitude in an ocean contained within a spherical shell due to inertial waves excited by the degree-two order-zero eccentricity tide. The forcing frequency equals the rotational frequency of the body. Dissipation is given for varying ocean thickness, varying values of inner to outer ocean radius ratios ($\eta$) and Ekman number ranging from $10^{-4}$ to $10^{-8}$. (b) is a zoom-in of (a) for the most relevant range of ocean thicknesses for Europa and Enceladus ($\eta = 0.8 - 1$).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Non-dimensional tidal dissipation amplitude, as defined in Equation (B.4), as a function of Ekman number for different representative inner to outer ocean radius ratio ($\eta$) for the degree-two order-zero eccentricity tidal component. $\eta = 0.63, 0.73$, corresponds to an ocean with wave attractors with high (in absolute value) Lyapunov exponent; for $\eta = 0.35$ there is an antiresonance; and $\eta = 0.75, 0.8815$ correspond to two resonant states. The markers indicate the different numerical experiments performed.</figcaption></figure>

years. However, so far tidal dissipation estimations using the Laplace Tidal Equations have resulted in low values of tidal dissipation for Enceladus' estimated ocean thickness (Tyler, 2014; Hay and Matsuyama, 2017; Matsuyama et al., 2018). It was suggested that internal inertial waves might produce enough additional tidal dissipation. Here we show that the direct generation of inertial waves in a spherical shell does not provide significant additional heat to prevent the ocean from freezing either.

However, the discussion on tidal heating in the subsurface oceans of the icy moons is far from settled. In our study we have made several assumptions that need to be revisited. We have assumed the ocean to be (1) unstratified, (2) of constant-thickness (no topography), and (3) we have neglected non-linear effects. Each of these elements are discussed in the next paragraphs.

We are limited by the absence of direct observations of the extraterrestrial oceans. We do not have any information about their density structure, thus, we can only apply our knowledge of similar environments on Earth to make an educated guess. Following the suggestions of Goodman et al. (2004); Goodman and Lenferink (2012) we have considered the ocean to be convectively mixed. However, under some circumstances the ocean can be stratified (Melosh et al., 2004; Vance and Brown, 2005). In a stratified ocean, internal gravity waves, which play an important role in tidal dissipation on Earth (Garrett, 2003),
can also be generated, adding a new possible source of tidal dissipation.

In this study we consider an ocean contained within a spherical shell. However, it is
expected that Enceladus' and Europa's ocean basins deviate from this idealised shape. It
is known, for example, that Enceladus' ocean does not have a constant thickness. On the
contrary, the ice shell thickness varies from 7 km at the south pole to 24 km at the equator
(Beuthe et al., 2016; Čadek et al., 2016). Moreover, due to rotation and tidal deformations
both Europa and Enceladus have the shape of a triaxial ellipsoid (e.g., Nimmo et al., 2007;
Thomas, 2010). Wave attractors also exist in more general ocean basins. We expect that the
conversion of large-scale tides into small scale inertial waves will be facilitated by additional
topographic variations. Inertial waves for complex ocean geometries is an exciting topic for
future work.

Besides, the barotropic tide given by the LTE highly depends on the ocean's basin
geometry. The interaction of the barotropic horizontal currents with topography excites
internal waves. On Earth most tidal dissipation occurs on the shallow continental shelf where
barotropic currents have a higher amplitude, but an additional ~30% of tidal dissipation
occurs in the oceanic ridges due to the excitation of internal waves (Egbert and Ray, 2000,
2003). The effect of topography on both the barotropic and internal tide should be further
studied.

Finally, in this study we have used perturbation theory, we have ignored the non-linear
terms in the momentum equations. The use of non-linear terms would change the behaviour
of the fluid. For example, inertial waves could interact with the underling convective flow.
Another possible interaction is that of inertial waves with a mean-flow excited by libration
which can lead to the well-known elliptic instability (Kerswell, 2002). This flow instabil-
ity gives rise to fully three dimensional turbulence which might enhance dissipation and
change its geographical pattern. Both Europa and Enceladus might be subject to this in-
stability (Lemasquerier et al., 2017). Further experimental and numerical work is needed to
understand the relevance of this phenomenon in heating the icy moons and shaping their
surface.

559
**Acknowledgement**

This research has been financially supported by the Space Research User Support pro-
gram of the Netherlands Organization for Scientific Research (NWO) under contract number
ALW-GO/16-19. The authors would like to thank the editor F. Nimmo and two anony-
mous reviewers for their valuable comments and suggestions which helped to improve the
manuscript.

565
Appendix A. Numerical approach expanded

We solve Equations (19) using the spectral method of Rieutord and Valdettaro (1997).
We expand the different fields using vector spherical harmonics:
$$ \mathbf{u}^{(\mathbf{d})} = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \left[ u_{m}^{l}(r) \mathbf{R}_{l}^{\mathbf{m}} + v_{m}^{l}(r) \mathbf{S}_{l}^{\mathbf{m}} + w_{m}^{l}(r) \mathbf{T}_{l}^{\mathbf{m}} \right] \qquad (A.1a) $$

568

$$ W = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} W_{m}^{l}(r) Y_{l}^{m}(\theta, \phi) \qquad (A.1b) $$

569

570 with:

$$ \mathbf{R}_{l}^{\mathbf{m}} = Y_{l}^{\mathbf{m}} e_{\mathbf{r}} \qquad \mathbf{S}_{l}^{\mathbf{m}} = \nabla Y_{l}^{\mathbf{m}} \qquad \mathbf{T}_{l}^{\mathbf{m}} = \nabla \times \mathbf{R}_{l}^{\mathbf{m}} \qquad (A.2) $$

571

572 where $Y_l^m$ are normalized spherical harmonics which satisfy the following orthogonal relations:

$$ \int_{4\pi} \mathbf{R}_{l}^{\mathbf{m}} \cdot \mathbf{S}_{l'}^{\mathbf{m}'} d\Omega = \int_{4\pi} \mathbf{S}_{l}^{\mathbf{m}} \cdot \mathbf{T}_{l'}^{\mathbf{m}'} d\Omega = \int_{4\pi} \mathbf{T}_{l}^{\mathbf{m}} \cdot \mathbf{R}_{l'}^{\mathbf{m}'} d\Omega = 0 \qquad (A.3a) $$

573

$$ \int_{4\pi} \mathbf{R}_{l}^{\mathbf{m}} \cdot \mathbf{R}_{l'}^{\mathbf{m}'} d\Omega = \delta_{ll'} \delta_{mm'} \qquad (A.3b) $$

574

$$ \int_{4\pi} \mathbf{S}_{l}^{\mathbf{m}} \cdot \mathbf{S}_{l'}^{\mathbf{m}'} d\Omega = \int_{4\pi} \mathbf{T}_{l}^{\mathbf{m}} \cdot \mathbf{T}_{l'}^{\mathbf{m}'} d\Omega = l(l+1)\delta_{ll'}\delta_{mm'} \qquad (A.3c) $$

575

Using the previous expansion the continuity equation becomes:

$$ v_m^l = \frac{1}{l(l+1)r} \frac{\partial r^2 u_m^l}{\partial r}. \qquad (A.4) $$

576

577

Taking the curl of Equation (15) to eliminate the potential and introducing Equation (A.4) we get the following Equations for the $\mathbf{R}_l^m$, and $\mathbf{T}_l^m$ velocity components:

$$ E \Delta_l w_m^l + \left( \frac{im}{l(l+1)} - i\omega \right) w_m^l = -A(l,m) r^{l-1} \frac{d}{dr} \left( \frac{u_m^{l-1}}{r^{l-2}} \right) - A(l+1,m) r^{-l-2} \frac{d}{dr} \left( r^{l+3} u_m^{l+1} \right) + f_R, \qquad (A.5a) $$

$$ E \Delta_l \Delta_l (r u_m^l) + \left( \frac{im}{l(l+1)} - i\omega \right) \Delta_l (r u_m^l) = B(l,m) r^{l-1} \frac{d}{dr} \left( \frac{w_m^{l-1}}{r^{l-1}} \right) + B(l+1,m) r^{-l-2} \frac{d}{dr} \left( r^{l+2} w_m^{l+1} \right) + f_T; \qquad (A.5b) $$

578

579 with:

$$ A(l, m) = \frac{1}{l^2} \sqrt{\frac{l^2 - m^2}{4l^2 - 1}}, \quad B(l, m) = l^2(l^2 - 1)A(l, m); \quad \Delta_l = \frac{1}{r} \frac{d^2}{dr^2} r - \frac{l(l+1)}{r^2}. \qquad (A.6) $$
580

The forcing terms $f_R$ and $f_T$ are given by:

$$f_R = -\frac{r}{l(l+1)} \nabla \times \mathbf{f}^{(\text{eq})} \cdot \mathbf{R}_l^m \qquad (\text{A.7a})$$

$$f_T = l(l+1)\nabla \times \mathbf{f}^{(\text{eq})} \cdot \mathbf{T}_l^m \qquad (\text{A.7b})$$

581

582

583

584

585

586

The values of $f_T$ and $f_R$ are obtained for the different tidal constituents. We need to obtain the non-dimensional components of the different constituents of the equilibrium tide $\mathbf{u}_{eq}$ in terms of spherical harmonics. We compute their amplitude by taking the time derivative of Equation (6) and using the definition of the equilibrium tide given in Equation (16). We write the resulting expressions using the previously defined normalised spherical harmonics $Y_{l,m}$:

$$\mathbf{u}_{e0}^{(\text{eq})} = -i\frac{3}{2}\sqrt{\frac{\pi}{5}}\frac{\omega^2 e R}{g}r^2 \mathbf{R}_2^0 = k_{e0}i r^2 \mathbf{R}_2^0 \qquad (\text{A.8a})$$

$$\mathbf{u}_{o1e}^{(\text{eq})} = -\sqrt{\frac{3\pi}{5}} \frac{\omega^2 R \sin \epsilon}{2g} r^2 \mathbf{R}_2^1 = -k_{o1e} r^2 \mathbf{R}_2^1 \qquad (\text{A.8b})$$

$$\mathbf{u}_{o1w}^{(\text{eq})} = \sqrt{\frac{3\pi}{5}} \frac{\omega^2 R \sin \epsilon}{2g} r^2 \mathbf{R}_2^1 = k_{o1w} r^2 \mathbf{R}_2^1 \qquad (\text{A.8c})$$

$$\mathbf{u}_{e2e}^{(\text{eq})} = i\frac{7}{4}\sqrt{\frac{6\pi}{5}}\frac{\omega^2 e R}{g}r^2 \mathbf{R}_2^2 = k_{e2e}i r^2 \mathbf{R}_2^2 \qquad (\text{A.8d})$$

$$\mathbf{u}_{e2w}^{(\text{eq})} = -i\frac{1}{4}\sqrt{\frac{6\pi}{5}}\frac{\omega^2 e R}{g}r^2 \mathbf{R}_2^2 = k_{e2w}i r^2 \mathbf{R}_2^2. \qquad (\text{A.8e})$$

587

588

Introducing Equation (A.8) into Equation (A.7) we can get the values of $f_R$ and $f_T$ for $m=0, m=1$ and $m=2$ components, ($f_R^0, f_T^0, f_R^1, f_T^1$ and $f_R^2, f_T^2$).

$$f_R^0 = k_{e0}i \frac{r^2}{\sqrt{35}} \delta_{l,3} - k_{e0}i \frac{r^2}{\sqrt{15}} \delta_{l,1}, \quad f_T^0 = 6k_{e0}\omega r \delta_{l,2} \qquad (\text{A.9a})$$

589

$$f_R^1 = k_{o1} \frac{2r^2}{3} \sqrt{\frac{2}{35}} \delta_{l,3} - k_{o1} \frac{r^2}{2} \sqrt{\frac{3}{15}} \delta_{l,1}, \quad f_T^1 = -6k_{o1}i r(\omega + \frac{1}{2})\delta_{l,2} \qquad (\text{A.9b})$$

590

$$f_R^2 = k_{e2}i \frac{r^2}{\sqrt{63}} \delta_{l,3}, \quad f_T^2 = 6k_{e2}\omega r \delta_{2,l} + 6k_{e2}r \delta_{l,2}. \qquad (\text{A.9c})$$

591

592

593

594

595

$k_{e2e}$ ($k_{o1e}$) or $k_{e2w}$ ($k_{o1w}$) should be used in the place of $k_{e2}$ ($k_{o1}$) depending if we solve for the westward or eastward wave. Note that the constant $k$ contains information about the physical characteristics of each moon. Due to the fact that the PDEs are linear we solve the system for $k=1$ and multiply the solution by the corresponding value of $k$ to obtain the solution for each configuration.
596

597

598

599

By inspecting Equations A.5 it is clear that the different orders are decoupled, but there exists a coupling between degrees. However, not all the degrees are coupled, symmetric ($\mathbf{X}_m^+$) and antisymmetric ($\mathbf{X}_m^-$) modes with respect to the equator are decoupled and can be solved independently.

$$ \mathbf{X}_m^+ = \begin{bmatrix} u_m^m \\ w_m^{m+1} \\ u_m^{m+2} \\ w_m^{m+3} \\ \vdots \end{bmatrix}, \qquad \mathbf{X}_m^- = \begin{bmatrix} w_m^m \\ u_m^{m+1} \\ w_m^{m+2} \\ u_m^{m+3} \\ \vdots \end{bmatrix}. \qquad (A.10) $$

600

601

Note that for the $m = 0$ and $m = 2$ tidal components the $\mathbf{X}_0^+$ and the $\mathbf{X}_2^+$ solutions are excited, respectively, while $m = 1$ excites the $\mathbf{X}_1^-$ solution.

602

603

604

The partial differential equations A.5 are discretised in the radial direction using the Chebyshev polynomials on the Gauss-Lobatto collocation nodes (see Section 3.3). More details about this method can be found in Rieutord and Valdettaro (1997).

605

## Appendix B. Tidal dissipation computation

606

We expand Equation (20) as:

$$ \hat{D}_v = E \int_{\eta}^{1} d(r) r^2 dr, \qquad (B.1) $$

607

608

where $d(r)$ is the dissipation in a spherical layer with radius $r$. We can write the previous equation using the spherical harmonic decomposition of the velocity (Equation (A.1)):

$$ d(r) = \sum_{l=0}^{L} \sum_{m=-L}^{L} 3 \left| \frac{\partial u_m^l}{\partial r} \right|^2 + l(l+1)(|s_m^l|^2 + |t_m^l|^2) + (l-1)l(l+1)(l+2) \frac{|v_m^l|^2 + |w_m^l|^2}{r^2}, \quad (B.2) $$

609

with:

$$ s_m^l = \frac{\partial v_m^l}{\partial r} + \frac{u_m^l - v_m^l}{r} \qquad (B.3a) $$

610

$$ t_m^l = r \frac{\partial}{\partial r} \frac{w_m^l}{r}. \qquad (B.3b) $$

611

612

e.g. Rieutord (1987). The total value of $\hat{D}_v$ is obtained by numerically integrating Equation (B.1). Finally, we obtain the dimensional value as:

$$ D_v = 8k^2 \rho_o \omega^3 R^5 \hat{D}_v \qquad (B.4) $$

613

where $k$ depends on the tidal component considered and is given by Equation (A.8).
614 References

615 Anderson, J. D., Schubert, G., Jacobson, R. A., Lau, E. L., Moore, W. B., Sjogren, W. L., 1998. Europa's differentiated internal structure: Inferences from four galileo encounters. *Science* 281 (5385), 2019–2022.

616 URL https://doi.org/10.1126/science.281.5385.2019

618 Baland, R.-M., Yseboodt, M., Van Hoolst, T., 2012. Obliquity of the Galilean satellites: The influence of a global internal liquid layer. *Icarus* 220 (2), 435–448.

620 URL https://doi.org/10.1016/j.icarus.2012.05.020

621 Bêhounková, M., Souček, O., Jaroslav, H., Čadek, O., 2017. Plume activity and tidal deformation on Enceladus influenced by faults and variable ice shell thickness. *Astrobiology* 17 (9), 941–954.

623 URL https://doi.org/10.1089/ast.2016.1629

624 Beuthe, M., 2016. Crustal control of dissipative ocean tides in Enceladus and other icy moons. *Icarus* 280, 278–299.

626 URL https://doi.org/10.1016/j.icarus.2016.08.009

627 Beuthe, M., Rivoldini, A., Trinh, A., 2016. Enceladus's and Dione's floating ice shells supported by minimum stress isostasy. *Geophys. Res. Lett.* 43 (19), 10,088–10,096.

629 URL http://dx.doi.org/10.1002/2016GL070650

630 Bretherton, F. P., 1964. Low frequency oscillations trapped near the equator. *Tellus* 16 (2), 181–185.
631 URL http://dx.doi.org/10.1111/j.2153-3490.1964.tb00159.x

632 Brouzet, C., Ermanyuk, E. V., Joubaud, S., Sibgatullin, I., Dauxois, T., 2016. Energy cascade in internal-wave attractors. *Europhys. Lett.* 113 (4), 44001.

634 URL https://doi.org/10.1209/0295-5075/113/44001

635 Bryan, G., 1889. The waves on a rotating liquid spheroid of finite ellipticity. *Phil. Trans. R. Soc. A* 180, 187–219.

637 URL http://doi.org/10.1098/rsta.1889.0006

638 Čadek, O., Tobie, G., Van Hoolst, T., Mass, M., Choblet, G., Lefèvre, A., Mitri, G., Baland, R.-M., Bêhounková, M., Bourgeois, O., Trinh, A., 2016. Enceladus's internal ocean and ice shell constrained from Cassini gravity, shape, and libration data. *Geophys. Res. Lett.* 43 (11), 5653–5660.

641 URL http://dx.doi.org/10.1002/2016GL068634

642 Cassen, P., Reynolds, R. T., Peale, S. J., 1979. Is there liquid water on Europa? *Geophys. Res. Lett.* 6 (9), 731–734.

644 URL http://dx.doi.org/10.1029/GL006i009p00731

645 Chen, E., Nimmo, F., Glatzmaier, G., 2014. Tidal heating in icy satellite oceans. *Icarus* 229, 11–30.

646 URL https://doi.org/10.1016/j.icarus.2013.10.024

647 Chen, E. M. A., Nimmo, F., 2011. Obliquity tides do not significantly heat Enceladus. *Icarus* 214 (2), 779–781.

649 URL https://doi.org/10.1016/j.icarus.2011.06.0078

650 Choblet, G., Tobie, G., Sotin, C., Bêhounková, M., Čadek, O., Postberg, F., Souček, O., 2017. Powering prolonged hydrothermal activity inside Enceladus. *Nat. Astron.* 1 (12), 841–847.

652 URL https://doi.org/10.1038/s41550-017-0289-8

653 Egbert, G. D., Ray, R. D., 2000. Significant dissipation of tidal energy in the deep ocean inferred from satellite altimeter data. *Nature* 405, 775.

655 URL http://dx.doi.org/10.1038/35015531 http://10.0.4.14/35015531

656 Egbert, G. D., Ray, R. D., 2003. Semidiurnal and diurnal tidal dissipation from TOPEX/Poseidon altimetry. *Geophys. Res. Lett.* 30 (17).

658 URL https://doi.org/10.1029/2003GL017676

659 Garrett, C., 2003. Internal tides and ocean mixing. *Science* 301 (5641), 1858–1859.

660 URL http://doi.org/10.1126/science.1090002

661 Gerkema, T., Zimmerman, J. T. F., Maas, L. R. M., van Haren, H., 2008. Geophysical and astrophysical fluid dynamics beyond the traditional approximation. *Rev. Geophys.* 46 (2), 1–33.

663 URL http://dx.doi.org/10.1029/2006RG000220

664 Goodman, J. C., Collins, G. C., Marshall, J., Pierrehumbert, R. T., 2004. Hydrothermal plume dynamics
665 on Europa: Implications for chaos formation. J. Geophys. Res. Planets 109 (E3).

666 URL http://dx.doi.org/10.1029/2003JE002073

667 Goodman, J. C., Lenferink, E., 2012. Numerical simulations of marine hydrothermal plumes for Europa and other icy worlds. Icarus 221 (2), 970–983.

668 URL https://doi.org/10.1016/j.icarus.2012.08.027

669 Greenspan., H. P., 1969. The Theory of Rotating Fluids. Cambridge University press, Cambridge Monographs on Mechanics and Applied Mathematics, New York.

670 Hay, H. C., Matsuyama, I., 2017. Numerically modelling tidal dissipation with bottom drag in the oceans of Titan and Enceladus. Icarus 281, 342 – 356.

671 URL https://doi.org/10.1016/j.icarus.2016.09.022

672 Hollerbach, R., Kerswell, R. R., 1995. Oscillatory internal shear layers in rotating and precessing flows. J. Fluid Mech. 298, 327–339.

673 URL https://doi.org/10.1017/S0022112095003338

674 Hsu, H.-W., Postberg, F., Sekine, Y., Shibuya, T., Kempf, S., Horányi, M., Juhász, A., Altobelli, N., Suzuki, K., Masaki, Y., Kuwatani, T., Tachibana, S., Sirono, S., Moragas-Klostermeyer, G., Srama, R., 2015. Ongoing hydrothermal activities within Enceladus. Nature 519, 207–210.

675 URL http://dx.doi.org/10.1038/nature14262

676 Hussmann, H., Choblet, G., Lainey, V., Matson, D. L., Sotin, C., Tobie, G., Van Hoolst, T., 2010. Implications of rotation, orbital states, energy sources, and heat transport for internal processes in icy satellites. Space Sci. Rev. 153 (1), 317–348.

677 URL https://doi.org/10.1007/s11214-010-9636-0

678 Hussmann, H., Spohn, T., Wieczerkowski, K., 2002. Thermal equilibrium states of europa's ice shell: Implications for internal ocean thickness and surface heat flow. Icarus 156 (1), 143 – 151.

679 URL https://doi.org/10.1006/icar.2001.6776

680 Jara-Orué, H. M., Vermeersen, B. L., 2011. Effects of low-viscous layers and a non-zero obliquity on surface stresses induced by diurnal tides and non-synchronous rotation: The case of Europa. Icarus 215 (1), 417 – 438.

681 URL https://doi.org/10.1016/j.icarus.2011.05.034

682 Kerswell, R. R., 1995. On the internal shear layers spawned by the critical regions in oscillatory ekman boundary layers. J. Fluid Mech. 298, 311–325.

683 URL http://doi.org/10.1017/S0022112095003326

684 Kerswell, R. R., 2002. Elliptical instability. Annu. Rev. Fluid Mech. 34 (1), 83–113.

685 URL https://doi.org/10.1146/annurev.fluid.34.081701.171829

686 Khurana, K. K., Kivelson, M. G., Stevenson, D. J., Schubert, G., Russell, C. T., Walker, R. J., Polanskey, C., 1998. Induced magnetic fields as evidence for subsurface oceans in Europa and Callisto. Nature 395 (6704), 777–780.

687 URL http://dx.doi.org/10.1038/27394

688 Kivelson, M. G., Khurana, K. K., Russell, C. T., Volwerk, M., Walker, R. J., Zimmer, C., 2000. Galileo magnetometer measurements: A stronger case for a subsurface ocean at Europa. Science 289 (5483), 1340–1343.

689 URL http://doi.org/10.1126/science.289.5483.1340

700 Lainey, V., Jacobson, R. A., Tajeddine, R., Cooper, N. J., Murray, C., Robert, V., Tobie, G., Guillot, T., Mathis, S., Remus, F., Desmars, J., Arlot, J.-E., Cuyper, J.-P. D., Dehant, V., Pascu, D., Thuillot, W., Poncin-Lafitte, C. L., Zahn, J.-P., 2017. New constraints on Saturn's interior from Cassini astrometric data. Icarus 281, 286 – 296.

701 URL https://doi.org/10.1016/j.icarus.2016.07.014

702 Le Bars, M., Cbron, D., Le Gal, P., 2015. Flows driven by libration, precession, and tides. Annual Review of Fluid Mechanics 47 (1), 163–193.

703 URL https://doi.org/10.1146/annurev-fluid-010814-014556

704 Lemasquerier, D., Grannan, A. M., Vidal, J., D., C., Favier, B., Le Bars, M., Aurnou, J. M., 2017. Librationdriven flows in ellipsoidal shells. J. Geophys. Res. Planets 122 (9), 1926–1950.
716 URL https://doi.org/10.1002/2017JE005340

717 Longuet-Higgins, M. S., 1968. The eigenfunctions of Laplace's tidal equations over a sphere. Phil. Trans. R. Soc. A 262 (1132), 511-607.

718 URL http://doi.org/10.1098/rsta.1968.0003

719 Maas, L. R. M., 2001. Wave focusing and ensuing mean flow due to symmetry breaking in rotating fluids. J. Fluid Mech. 437, 13-28.

720 URL http://doi.org/10.1017/S0022112001004074

721 Maas, L. R. M., 2005. Wave attractors: Linear yet nonlinear. Int. J. Bifurcat. Chaos 15 (09), 2757-2782.

722 URL http://doi.org/10.1142/S0218127405013733

723 Maas, L. R. M., Benielli, D., Sommeria, J., Lam, F.-P. A., Aug. 1997. Observation of an internal wave attractor in a confined, stably stratified fluid. Nature 388 (6642), 557-561.

724 URL https://doi.org/10.1038/41509

725 Manders, A. M. M., Maas, L. R. M., 2003. Observations of inertial waves in a rectangular basin with one sloping boundary. J. Fluid Mech. 493, 59-88.

726 URL http://doi.org/10.1017/S0022112003005998

727 Manders, A. M. M., Maas, L. R. M., 2004. On the three-dimensional structure of the inertial wave field in a rectangular basin with one sloping boundary. Fluid Dyn. Res. 5 (1), 1-21.

728 URL https://doi.org/10.1016/j.fluiddyn.2004.03.004

729 Matsuyama, I., Beuthe, M., Hay, H. C., Nimmo, F., Kamata, S., 2018. Ocean tidal heating in icy satellites with solid shells. Icarus 312, 208-230.

730 URL https://doi.org/10.1016/j.icarus.2018.04.013

731 Melosh, H., Ekholm, A., Showman, A., Lorenz, R., 2004. The temperature of Europa's subsurface water ocean. Icarus 168 (2), 498-502.

732 URL https://doi.org/10.1016/j.icarus.2003.11.021

733 Nimmo, F., Thomas, P., Pappalardo, R., Moore, W., 2007. The global shape of europa: Constraints on lateral shell thickness variations. Icarus 191 (1), 183-192.

734 URL https://doi.org/10.1016/j.icarus.2007.04.021

735 Ogilvie, G. I., 2005. Wave attractors and the asymptotic dissipation rate of tidal disturbances. J. Fluid Mech. 543, 19-44.

736 URL http://doi.org/10.1017/S0022112005006580

737 Ogilvie, G. I., 2009. Tidal dissipation in rotating fluid bodies: a simplified model. Mon. Not. R. Astron. Soc. 396 (2), 794-806.

738 URL http://dx.doi.org/10.1111/j.1365-2966.2009.14814.x

739 Ogilvie, G. I., Lin, D. N. C., 2004. Tidal dissipation in rotating giant planets. Astrophys. J. 610 (1), 477-509. URL http://doi.org/10.1086/421454

740 Ojakangas, G. W., Stevenson, D. J., 1989. Thermal state of an ice shell on Europa. Icarus 81 (2), 220-241. URL https://doi.org/10.1016/0019-1035(89)90052-33

741 Pappalardo, R. T., Belton, M. J. S., Breneman, H. H., Carr, M. H., Chapman, C. R., Collins, G. C., Denk, T., Fagents, S., Geissler, P. E., Giese, B., Greeley, R., Greenberg, R., Head, J. W., Helfenstein, P., Hoppa, G., Kadel, S. D., Klaasen, K. P., Klemaszewski, J. E., Magee, K., McEwen, A. S., Moore, J. M., Moore, W. B., Neukum, G., Phillips, C. B., Prockter, L. M., Schubert, G., Senske, D. A., Sullivan, R. J., Tufts, B. R., Turtle, E. P., Wagner, R., Williams, K. K., 1999. Does Europa have a subsurface ocean? Evaluation of the geological evidence. J. Geophys. Res. Planets 104 (E10), 24015-24055.

742 URL http://dx.doi.org/10.1029/1998JE000628

743 Porco, C. C., Helfenstein, P., Thomas, P. C., Ingersoll, A. P., Wisdom, J., West, R., Neukum, G., Denk, T., Wagner, R., Roatsch, T., Kieffer, S., Turtle, E., McEwen, A., Johnson, T. V., Rathbun, J., Veverka, J., Wilson, D., Perry, J., Spitale, J., Brahic, A., Burns, J. A., DelGenio, A. D., Dones, L., Murray, C. D., Squyres, S., 2006. Cassini observes the active South Pole of Enceladus. Science 311 (5766), 1393-1401.

744 URL http://doi.org/10.1126/science.1123013

745 Postberg, F., Kempf, S., Schmidt, J., Brilliantov, N., Beinsen, A., Abel, B., Buck, U., Srama, R., 2009. Sodium salts in E-ring ice grains from an ocean below the surface of Enceladus. Nature 459, 1098-1101.

746
767 URL http://dx.doi.org/10.1038/nature08046

768 Postberg, F., Schmidt, J., Hillier, J., Kempf, S., Srama, R., 2011. A salt-water reservoir as the source of a compositionally stratified plume on Enceladus. Nature 474, 620–622.

769 URL http://dx.doi.org/10.1038/nature10175

771 Rabitti, A., Maas, L. R., 2014. Inertial wave rays in rotating spherical fluid domains. J. Fluid Mech. 758, 621–654.

773 URL https://doi.org/10.1017/jfm.2014.551

774 Rabitti, A., Maas, L. R. M., 2013. Meridional trapping and zonal propagation of inertial waves in a rotating fluid shell. J. Fluid Mech. 729, 445–470.

776 URL https://doi.org/10.1017/jfm.2013.310

777 Rieutord, M., 1987. Linear theory of rotating fluids using spherical harmonics part I: Steady flows. Geophys. Astro. Fluid 39 (3), 163–182.

779 URL http://dx.doi.org/10.1080/03091928708208811

780 Rieutord, M., Georgeot, B., Valdettaro, L., 2001. Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum. J. Fluid Mech. 435, 103–144.

782 URL http://doi.org/10.1017/S0022112001003718

783 Rieutord, M., Valdettaro, L., 1997. Inertial waves in a rotating spherical shell. J. Fluid Mech. 341, 77–99.
URL http://doi.org/10.1017/S0022112097005491

785 Rieutord, M., Valdettaro, L., 2010. Viscous dissipation by tidally forced inertial modes in a rotating spherical shell. J. Fluid Mech. 643, 363–394.

787 URL http://doi.org/10.1017/S002211200999214X

788 Rieutord, M., Valdettaro, L., Jun. 2018. Axisymmetric inertial modes in a spherical shell at low Ekman numbers. J. Fluid Mech. 844, 597–634.

790 URL http://doi.org/10.1017/jfm.2018.201

791 Roberts, J. H., Nimmo, F., 2008. Tidal heating and the long-term stability of a subsurface ocean on Enceladus. Icarus 194 (2), 675 – 689.

793 URL https://doi.org/10.1016/j.icarus.2007.11.010

794 Soderlund, K. M., Schmidt, B. E., Wicht, J., Blankenship, D. D., 2013. Ocean-driven heating of Europa's icy shell at low latitudes. Nat. Geosci. 7, 16–19.

796 URL http://dx.doi.org/10.1038/ngeo2021

797 Spencer, J. R., Howett, C. J. A., Verbiser, A., Hurford, T. A., Segura, M., Spencer, D. C., Sep. 2013. Enceladus heat flow from high spatial resolution thermal emission observations. European Planetary Science Congress Abstracts 8.

800 Spencer, J. R., Pearl, J. C., Segura, M., Flasar, F. M., Mamoutkine, A., Romani, P., Buratti, B. J., Hendrix, A. R., Spilker, L. J., Lopes, R. M. C., 2006. Cassini encounters Enceladus: Background and the discovery of a south polar hot spot. Science 311 (5766), 1401–1405.

803 URL http://doi.org/10.1126/science.1121661

804 Stern, M. E., 1963. Trapping of low frequency oscillations in an equatorial boundary layer. Tellus 15 (3), 246–250.

806 URL http://dx.doi.org/10.1111/j.2153-3490.1963.tb01384.x

807 Stewartson, K., 1971. On trapped oscillations of a rotating fluid in a thin spherical shell. Tellus 23 (6), 506–510.

809 URL http://dx.doi.org/10.1111/j.2153-3490.1971.tb00598.x

810 Stewartson, K., 1972. On trapped oscillations of a rotating fluid in a thin spherical shell ii. Tellus 24 (4), 283–287.

812 URL http://dx.doi.org/10.1111/j.2153-3490.1972.tb01555.x

813 Thomas, P., 2010. Sizes, shapes, and derived properties of the saturnian satellites after the cassini nominal mission. Icarus 208 (1), 395 – 401.

815 URL https://doi.org/10.1016/j.icarus.2010.01.025

816 Thomas, P., Tajeddine, R., Tiscareno, M., Burns, J., Joseph, J., Loredo, T., Helfenstein, P., Porco, C., 2016. Enceladus's measured physical libration requires a global subsurface ocean. Icarus 264, 37 – 47.

817
818 URL https://doi.org/10.1016/j.icarus.2015.08.037

819 Tobie, G., Čadek, O., Sotin, C., 2008. Solid tidal friction above a liquid water reservoir as the origin of the south pole hotspot on Enceladus. Icarus 196 (2), 642–652.

820 URL https://doi.org/10.1016/j.icarus.2008.03.008

821 Tyler, R. H., 2008. Strong ocean tidal flow and heating on moons of the outer planets. Nature 456 (7223), 770–772.

822 URL http://dx.doi.org/10.1038/nature07571

823 Tyler, R. H., 2009. Ocean tides heat Enceladus. Geophys. Res. Lett. 36 (15), L15205.

824 URL http://dx.doi.org/10.1029/2009GL038300

825 Tyler, R. H., 2014. Comparative estimates of the heat generated by ocean tides on icy satellites in the outer solar system. Icarus 243 (Supplement C), 358 – 385.

826 URL https://doi.org/10.1016/j.icarus.2014.08.037

827 Vallis, G. K., 2006. Atmospheric and Oceanic Fluid Dynamics. Cambridge University Press, Cambridge, U.K.

828 Vance, S., Brown, J., 2005. Layering and double-diffusion style convection in Europa's ocean. Icarus 177 (2), 506 – 514.

829 URL https://doi.org/10.1016/j.icarus.2005.06.005

830 Vermeersen, B. L., Maas, L. R., van Oers, S., Rabitti, A., Jara-Orue, H., 2013. Tidal-Induced Ocean Dynamics as Cause of Enceladus' Tiger Stripe Pattern. AGU Fall Meeting Abstracts, P53B–1848.