Thin-shell Tidal Dynamics of Ocean Worlds

We consider an ocean of total thickness H_ocean that is partitioned into N layers of densities ρ_i , where i stands for the index of the layer. Additionally, we consider that the ocean might be covered by a solid shell of thickness H_shell. A sketch of the system for a two-layer ocean is shown in Figure 1. The thickness of each layer is given by h_i , which we further split into its unperturbed average thickness H_i and a perturbation,

$$\begin{eqnarray}&&{h}_{i}={H}_{i}+{\eta }_{i}-{\eta }_{i+1}.\end{eqnarray} \tag{ 1 }$$

For each layer, we can write the mass conservation equation as

$$\begin{eqnarray}&&{\partial }_{t}({\eta }_{i}-{\eta }_{i+1})+{\rm{\nabla }}\cdot {{\boldsymbol{s}}}_{i}=0,\end{eqnarray} \tag{ 2 }$$

where s _i is the momentum flux, defined as

$$\begin{eqnarray}&&{{\boldsymbol{s}}}_{i}={\int }_{{h}_{i}}{{\boldsymbol{u}}}_{i}{dr},\end{eqnarray} \tag{ 3 }$$

and u is the tangential flow velocity.

Following the thin-shell approximation, we consider that the radial momentum equation is simply given by hydrostatic equilibrium and obtain the momentum equation of each layer by vertically averaging the other components of the momentum equations

$$\begin{eqnarray}&&{\partial }_{t}{{\boldsymbol{s}}}_{i}-f{{\boldsymbol{s}}}_{i}\times {{\bf{1}}}_{r}=-\displaystyle \frac{{h}_{i}}{{\rho }_{i}}{\rm{\nabla }}{p}_{i}-{h}_{i}{\rm{\nabla }}\phi -D({h}_{i},{{\boldsymbol{s}}}_{i}),\end{eqnarray} \tag{ 4 }$$

where _pi is the pressure within the layer. D(h_i , s _i ) is a function that parameterizes dissipation within the layer, ϕ is the gravitational potential, 1_r is the radial unit vector, and f is the Coriolis parameter

$$\begin{eqnarray}&&f=2{\rm{\Omega }}\cos \theta ,\end{eqnarray} \tag{ 5 }$$

with θ being the colatitude and Ω the body's rotational frequency. We have assumed that perturbations are small and linearized the equations of motion. This makes the problem analytically tractable and allows us to obtain wave solutions, but it implies that nonlinear effects that can lead to instabilities and turbulence are ignored (see discussion in Section 3.3).

The pressure within each layer can be obtained using hydrostatic balance

$$\begin{eqnarray}\begin{array}{rcl}{p}_{i} & = & q+g\left[\Space{0ex}{3.55ex}{0ex}{\rho }_{\mathrm{shell}}\left({H}_{\mathrm{shell}}+{\eta }_{0}-{\eta }_{1}\right)\right.\\ & & +\,\displaystyle \sum _{k=1}^{i-1}{\rho }_{k}({H}_{k}+{\eta }_{k}-{\eta }_{k+1})+{\rho }_{i}\\ & & \times \left.\left(-z+{\eta }_{i}-\displaystyle \sum _{k=1}^{i-1}{H}_{k}-{H}_{\mathrm{shell}}\right)\right],\end{array}\end{eqnarray} \tag{ 6 }$$

where q is the dynamic pressure exerted by the solid shell on the ocean, η₀ is the displacement of the solid surface, and z is measured from the surface and is positive upward. From the previous expression it is evident that even though there might be no perturbations (η = 0), a pressure gradient driving ocean flows can arise if layers are not constant in thickness. We consider that this is compensated by a transport mechanism that maintains the stratification structure. As we use the linearized LTEs, we assume that this flow does not interact with tidally driven currents and only consider flows induced by gradients in η.

The dissipative function D can take different forms depending on the leading dissipation mechanism (e.g., Hay & Matsuyama 2017)

$$\begin{eqnarray}&&D({h}_{i},{{\boldsymbol{u}}}_{i})=\alpha {h}_{i}{{\boldsymbol{u}}}_{i}+{c}_{D}|{{\boldsymbol{u}}}_{i}|{{\boldsymbol{u}}}_{i}+{h}_{i}\nu {\rm{\Delta }}{{\boldsymbol{u}}}_{i}.\end{eqnarray} \tag{ 7 }$$

The first term is a linear dissipation term; ocean dissipation is proportional to the velocity via the Rayleigh parameter α. It is akin to the viscous term of a damped oscillator and considers that dissipative processes have a characteristic timescale equal to the inverse of the Rayleigh parameter for all spatial and temporal scales. The second term is proportional to the square of the velocity via the bottom drag coefficient c_D and is often used to account for dissipation in fluid-solid boundaries due to friction. The third term is a diffusive term; ν can be either the dynamic or an eddy viscosity. On Earth tidal dissipation occurs mostly as a result of bottom drag and is well captured by the second term. There exist empirical constraints on the value of c_D from Earth studies; because of this, several studies of tides in icy moons have considered bottom drag as the main dissipative source (Chen & Nimmo 2016; Hay & Matsuyama 2017, 2019). However, this has the disadvantage of introducing a nonlinear term in the equations of motion and might not capture other dissipative processes at play in subsurface oceans. For these reasons many studies have used a linear drag term (Tyler 2011, 2014; Matsuyama 2014; Beuthe 2016; Matsuyama et al. 2018; Hay et al. 2020; Rovira-Navarro et al. 2020). We will assume a linear dissipation term that can account for different dissipative processes and vary α over various orders of magnitude. It is useful to note that scaling analysis can be used to obtain an effective Rayleigh coefficient α for a given value of bottom drag coefficient c_D . Matsuyama et al. (2018) showed α to be between 10⁻⁵ and 10⁻¹¹ s⁻¹ for Europa and Enceladus.

Plugging Equation (6) into Equation (4) and considering Rayleigh dissipation, we find the linearized thin-shell equations

$$\begin{eqnarray}&&{\partial }_{t}{\eta }_{i}-{\partial }_{t}{\eta }_{i+1}+{\rm{\nabla }}\cdot ({{\boldsymbol{s}}}_{i})=0,\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&\begin{array}{c}{{\rm{\partial }}}_{t}{{\boldsymbol{s}}}_{i}+\alpha {{\boldsymbol{s}}}_{i}-f{{\boldsymbol{s}}}_{{\boldsymbol{i}}}\times {{\bf{1}}}_{r}\\ \text{}=\,-{c}_{i}^{2}{\rm{\nabla }}\left({\eta }_{i}-{\eta }_{i+1}+\displaystyle \sum _{k\gt i}({\eta }_{k}-{\eta }_{k+1})\right.\\ \text{}\text{}+\,\left.\displaystyle \sum _{k\lt i}\displaystyle \frac{{\rho }_{k}}{{\rho }_{i}}({\eta }_{k}-{\eta }_{k+1})+\displaystyle \frac{q}{g{\rho }_{i}}+\displaystyle \frac{\phi }{g}\right).\end{array}\end{eqnarray} \tag{ 8b }$$

Here c_i ² is the square of the surface wave speed $\sqrt{{{gh}}_{i}}$; we will assume that displacements η are smaller than the layer thickness and hence use ${c}_{i}^{2}={{gH}}_{i}$.

Viscous processes within the ocean cause energy dissipation. In our approach, dissipative processes are parameterized via the Rayleigh coefficient α. The average amount of tidal dissipation in layer i for a period T is given by

$$\begin{eqnarray}&&{\dot{E}}_{i}=\displaystyle \frac{1}{T}{\int }_{T}{\int }_{V}{\rho }_{i}\alpha {{\boldsymbol{u}}}_{i}\cdot {{\boldsymbol{u}}}_{i}{dV}.\end{eqnarray} \tag{ 9 }$$

Alternatively, the total tidal dissipation can be computed by considering the work done by the tidal force in the ocean (Chen et al. 2014),

$$\begin{eqnarray}&&\dot{E}=-\displaystyle \frac{1}{T}{\int }_{T}{\int }_{V}\rho {\rm{\nabla }}{\phi }^{T}\cdot {\boldsymbol{u}}{dV}.\end{eqnarray} \tag{ 10 }$$

The pressure term q couples ocean and shell dynamics. q is related to the displacement of the ocean surface, η₁. As mentioned above, we use thin-shell theory for shells of variable thickness (Beuthe 2008, 2018). For shells of constant thickness, the pressure load results in a normal load at the bottom of the shell. In contrast, for shells of variable thickness a lateral load proportional to the slope of the shell−ocean interface also arises. As we consider long-wavelength thickness variations, which are characterized by small slopes, we neglect this component and assume that the ocean pressure exerts only a normal load. We use the massless approach (Beuthe 2015a), which assumes that the density contrast between the ocean and solid shell is zero and adds the shell mass to the ocean. This implies that we neglect the gravitational perturbation that results from the ocean−shell density contrast, which is small owing to the small density contrast and hence results in a small error (e.g., 2% in the case of Enceladus; Beuthe 2018). The solid shell dynamics is then described by

$$\begin{eqnarray}&&\begin{array}{c}{ \mathcal M }(\psi \chi {D}^{(b)},{\eta }_{0})-(1-\nu ){ \mathcal N }(\psi \chi {D}^{(b)},{\eta }_{0})\\ \quad +\,{R}^{3}{ \mathcal N }(\chi ,F)={R}^{4}q\end{array}\end{eqnarray} \tag{ 11a }$$

$$\begin{eqnarray}&&\begin{array}{c}{ \mathcal M }(\psi \chi {\alpha }^{(e)},F)-(1+\nu ){ \mathcal N }(\psi \chi {\alpha }^{(e)},F)\\ \,-\,{R}^{-1}{ \mathcal N }(\chi ,{\eta }_{0})=0.\end{array}\end{eqnarray} \tag{ 11b }$$

F is an auxiliary variable called the stress function, and η₀ is the radial surface displacement. From the thin-shell approximation it follows that the transverse displacement is constant within the shell, and hence η₁ = η₀. ${ \mathcal M }$ and ${ \mathcal N }$ are two differential operators defined in Appendix F. α^(e) and D^(b) are the inverse of the extensional rigidity and the bending rigidity:

$$\begin{eqnarray}&&{\alpha }^{(e)}=\frac{1}{2(1+\nu ){\mu }_{0}{H}_{\mathrm{shell}}},\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&{D}^{(b)}=\frac{2{\mu }_{\mathrm{inv}}{H}_{\mathrm{shell}}^{3}}{1-\nu }.\end{eqnarray} \tag{ 12b }$$

Here ν is the Poisson ratio and μ₀ is the vertically averaged shear modulus, or more generally the p = 0 moment of the shear modulus,

$$\begin{eqnarray}&&{\mu }_{p}=\displaystyle \frac{1}{{H}_{\mathrm{shell}}^{p+1}}{\int }_{{H}_{\mathrm{shell}}}\mu {h}^{p}{dh},\end{eqnarray} \tag{ 13 }$$

with h being the distance from a reference surface within the shell. Parameter μ_inv is given by

$$\begin{eqnarray}&&{\mu }_{\mathrm{inv}}={\mu }_{2}-\displaystyle \frac{{\mu }_{1}}{{\mu }_{0}}.\end{eqnarray} \tag{ 14 }$$

χ and ψ are two purely geometrical parameters given by

$$\begin{eqnarray}&&\left(\chi ,\psi \right)=\left(\displaystyle \frac{{\mu }_{0}+\delta {\mu }_{1}}{{\mu }_{0}+2\delta {\mu }_{1}+{\delta }^{2}{\mu }_{2}},\displaystyle \frac{{\mu }_{0}}{{\mu }_{0}+\delta {\mu }_{1}}\right),\end{eqnarray} \tag{ 15 }$$

with δ = H_shell/R. Apart from ocean currents and shell deformations, the tidal potential also deforms the mantle. Additionally, the mantle deforms in response to the ocean load and shell pressure q. As we will see, the effects of a nonrigid mantle can be accounted for using the tidal, load, and pressure Love numbers of the mantle (Matsuyama 2014; Beuthe 2016; Matsuyama et al. 2018).

The tidal response is obtained by solving Equations (8) and (11b). We use a semianalytical method based on spherical harmonics. The method is similar to that previously used to solve the LTEs for oceans of constant thickness and density (Longuet-Higgins 1968; Tyler 2008; Chen et al. 2014; Matsuyama 2014; Beuthe 2016) but is extended to accommodate both ocean and solid shell thickness variations, as well as ocean stratification.

We first split the gravitational potential into two components, one corresponding to the forcing tidal potential ϕ^T, and the other, ϕ^SG, arising from the deformation of the body and its self-gravity. We also rewrite the forcing term using

$$\begin{eqnarray}&&{\eta }^{{\rm{T}}}=-\displaystyle \frac{{\phi }^{{\rm{T}}}}{g},\end{eqnarray} \tag{ 16 }$$

where η^T is the displacement of the ocean surface if ocean dynamics, self-gravity, and the presence of a shell are ignored. We decompose s _i using a Helmholtz decomposition:

$$\begin{eqnarray}&&{{\boldsymbol{s}}}_{i}={\rm{\nabla }}{{\rm{\Phi }}}_{i}+{\rm{\nabla }}\times ({{\rm{\Psi }}}_{i}{{\bf{1}}}_{r}).\end{eqnarray} \tag{ 17 }$$

Using these definitions and taking ∇ · and 1_r · ∇ × of Equation (8b), Equation (8) can be rewritten in terms of Φ_i and Ψ_i instead of s _i (see Appendix A).

The resulting equations are solved in the spectral domain using spherical harmonics. The forcing term η^T is decomposed as

$$\begin{eqnarray}&&{\eta }^{{\rm{T}}}=\displaystyle \sum _{l,m}{\eta }_{l,m}^{{\rm{T}}}{Y}_{l}^{m}\exp (i\omega t),\end{eqnarray} \tag{ 18 }$$

where ω is the frequency of the forcing and Y_l ^m is a spherical harmonic of degree l and order m,

$$\begin{eqnarray}&&{Y}_{l}^{m}={\left(-1\right)}^{m}\sqrt{(2l+1)\displaystyle \frac{(l-m)!}{(l+m)!}}{P}_{l,m}(\cos \theta )\exp ({im}\varphi ).\end{eqnarray} \tag{ 19 }$$

P_l,m are associated Legendre functions (e.g., Ivers & Phillips 2008; Arfken et al. 2013).

The solution is similarly expanded as

$$\begin{eqnarray}\begin{array}{rcl}({\eta }_{i},{{\rm{\Phi }}}_{i},{{\rm{\Psi }}}_{i},q,F) & = & \displaystyle \sum _{l,m}({\eta }_{i,l,m},{{\rm{\Phi }}}_{i,l,m},{{\rm{\Psi }}}_{i,l,m},{q}_{l,m},{F}_{l,m})\\ & & \times \,{Y}_{l}^{m}\exp (i\omega t).\end{array}\end{eqnarray} \tag{ 20 }$$

We split the parameters c_i ², H_shell, χ, (ψ χ α^(e)), and (ψ χ D^(b)) into a uniform and a variable component:

$$\begin{eqnarray}&&\begin{array}{c}({c}_{i}^{2},{H}_{\mathrm{shell}},\chi ,\chi \psi {\alpha }^{(e)},\chi \psi {D}^{(b)})\\ \,=\,({c}_{i,0}^{2},{H}_{\mathrm{shell},0},{\chi }_{0},{(\chi \psi {\alpha }^{(e)})}_{0},{(\chi \psi {D}^{(b)})}_{0})\\ \,+\,\displaystyle \sum _{l\ne 0,m}({c}_{i,l,m}^{2},{H}_{\mathrm{shell},{\rm{l}},{\rm{m}}},{\chi }_{l,m},{(\chi \psi {\alpha }^{(e)})}_{l,m},{(\chi \psi {D}^{(b)})}_{l,m}){Y}_{l}^{m}.\end{array}\end{eqnarray} \tag{ 21 }$$

As lateral variations must be real, we have that if there are sectoral or tesseral variations (m ≠ 0), the expansion coefficients A_l,m must satisfy ${A}_{l,-m}={\left(-\right)}^{m}{A}_{l,m}^{* }$, where * indicates the complex conjugate. Similarly, the Coriolis term—Equation (5)—can be written using spherical harmonics as

$$\begin{eqnarray}&&f=\displaystyle \frac{2{\rm{\Omega }}}{\sqrt{3}}{Y}_{1}^{0}.\end{eqnarray} \tag{ 22 }$$

The gravitational potential arising from the deformation of the body (ϕ^SG) consists of a component resulting from the deformation of the outer layers (ϕ^L) and a component arising as a result of the deformation of the mantle. If the mantle is spherically symmetric and we use the thin-shell approximation, ϕ^SG can be written in terms of Love numbers as (Beuthe 2016)

$$\begin{eqnarray}&&{\phi }_{l,m}^{{\rm{SG}}}={k}_{l}^{{\rm{T}}}{\phi }_{l,m}^{{\rm{T}}}+(1+{k}_{l}^{{\rm{L}}}){\phi }_{l,m}^{{\rm{L}}}+{k}_{l}^{{\rm{P}}}{\phi }_{l,m}^{{\rm{P}}}.\end{eqnarray} \tag{ 23 }$$

Here k_l ^T, k_l ^L, k_l ^P are the gravitational tidal, load, and pressure Love numbers, respectively. ${\phi }_{l,m}^{{\rm{L}}}$ and ${\phi }_{l,m}^{{\rm{P}}}$ are the load and pressure potentials, which, assuming that density changes within the ocean are small Δρ/ρ₁ ≪ 1, are given by

$$\begin{eqnarray}&&{\phi }_{l,m}^{{\rm{L}}}=-\frac{3}{2l+1}\frac{{\rho }_{1}}{{\rho }_{b}}g{\eta }_{l,m}=-{\xi }_{l}g{\eta }_{l,m},\,\,{\phi }_{l,m}^{{\rm{P}}}=-{\xi }_{l}\frac{{q}_{l,m}}{{\rho }_{1}}.\end{eqnarray} \tag{ 24 }$$

Here ρ_b is the bulk density of the body and η_l,m is the difference between the displacement of the top and the bottom of the ocean η_1,l,m − η_N+1,l,m . The deformation of the bottom of the ocean is similarly given by

$$\begin{eqnarray}&&{\eta }_{N+1,l,m}=-\frac{{h}_{l}^{{\rm{T}}}}{g}{\phi }_{l,m}^{{\rm{T}}}-\frac{{h}_{l}^{{\rm{L}}}}{g}{\phi }_{l,m}^{{\rm{L}}}-\frac{{h}_{l}^{{\rm{P}}}}{g}{\phi }_{l,m}^{{\rm{P}}},\end{eqnarray} \tag{ 25 }$$

where h_l ^T, h_l ^L, h_l ^P are the radial Love numbers. If the mantle is rigid, the Love numbers are 0; Appendix D provides analytical expressions of the Love numbers for a homogeneous mantle. Finally, it is useful to nondimensionalize the equations of motion. We use the definitions

$$\begin{eqnarray}&&\begin{array}{c}(\hat{\eta },\hat{{\rm{\Phi }}},\hat{{\rm{\Psi }}},\hat{\phi },\hat{q},\hat{F},\hat{\omega },\hat{\alpha },\hat{c},{\hat{\alpha }}^{\left(e\right)},{\hat{D}}^{\left(b\right)})\\ \,=\,(\frac{\eta }{R},\frac{{\rm{\Phi }}}{{R}^{3}{\rm{\Omega }}},\frac{{\rm{\Psi }}}{{R}^{3}{\rm{\Omega }}},\frac{\phi }{{R}^{2}{{\rm{\Omega }}}^{2}},\frac{q}{{\rho }_{1}{gR}},{\alpha }_{0}^{\left(e\right)}F,\\ \,\,\frac{\omega }{{\rm{\Omega }}},\frac{\alpha }{{\rm{\Omega }}},\frac{c}{{c}_{0}},\frac{{\alpha }^{\left(e\right)}}{{\alpha }_{0}^{\left(e\right)}},\frac{{\alpha }_{0}^{\left(e\right)}}{{R}^{2}}{D}^{\left(b\right)}),\end{array}\end{eqnarray} \tag{ 26 }$$

with ${c}_{0}=\sqrt{{{gH}}_{\mathrm{ocean}}}$ and ${\alpha }_{0}^{(e)}=1/2(1+\nu ){\mu }_{0}{H}_{\mathrm{shell},0}$.

Using Equations (20) and (21) and projecting Equations (11b) and (A1) onto spherical harmonics employing the properties listed in Appendix G, we can get an algebraic set of differential equations. First, from mass conservation in the ocean (Equation (A1a)), we relate ${\hat{\eta }}_{i}$ and ${\hat{{\rm{\Phi }}}}_{i}$ as

$$\begin{eqnarray}&&{\hat{\eta }}_{i,l,m}={\hat{\eta }}_{i+1,l,m}-i\displaystyle \frac{l(l+1)}{\hat{\omega }}{\hat{{\rm{\Phi }}}}_{i,l,m}.\end{eqnarray} \tag{ 27 }$$

The previous expression can be used recursively to obtain ${\hat{\eta }}_{l,m}$,

$$\begin{eqnarray}&&{\hat{\eta }}_{l,m}=-\displaystyle \sum _{i}i\displaystyle \frac{l(l+1)}{\hat{\omega }}{\hat{{\rm{\Phi }}}}_{i,l,m}.\end{eqnarray} \tag{ 28 }$$

Using Equations (23), (25), and (27), we write Equations (11b), (A1b), and (A1c) as

$$\begin{eqnarray}&&\begin{array}{c}\left[i\hat{\omega }+\hat{\alpha }+{{ \mathcal C }}_{l,m}^{\left(1\right)}-\frac{i{{ \mathcal D }}_{l}{\hat{c}}_{i,0}^{2}}{\hat{\omega }\epsilon }\left(1-{\xi }_{l}\left(1+{\upsilon }_{l}^{L}\right)\right)\right]{\hat{{\rm{\Phi }}}}_{i,l,m}+{{ \mathcal C }}_{l,m}^{\left(2\right)}{\hat{{\rm{\Psi }}}}_{i,l-1,m}+{{ \mathcal C }}_{l,m}^{\left(3\right)}{\hat{{\rm{\Psi }}}}_{i,l+1,m}+\mathop{\underbrace{\frac{{\hat{c}}_{i,0}^{2}}{\epsilon }\left(\frac{{\rho }_{1}}{{\rho }_{i}}-{\xi }_{l}{\upsilon }_{l}^{{\rm{P}}}\right){\hat{q}}_{l,m}}}\limits_{{\rm{s}}{\rm{h}}{\rm{e}}{\rm{l}}{\rm{l}}}\\ +\,\displaystyle \mathop{\displaystyle \underbrace{{\sum }_{\mathop{{l}_{1},{m}_{1}}\limits_{{l}_{2},{m}_{2}}}{{{ \mathcal T }}^{\left(1\right)}}_{{l}_{1},{m}_{1},{l}_{2},{m}_{2}}^{l,m}\frac{{\hat{c}}_{i,{l}_{1},{m}_{1}}^{2}}{\epsilon }\left[-\frac{i{{ \mathcal D }}_{{l}_{2}}}{\hat{\omega }}\left(\left(1-{\xi }_{{l}_{2}}\left(1+{\upsilon }_{{l}_{2}}^{{\rm{L}}}\right)\right){\hat{{\rm{\Phi }}}}_{i,{l}_{2},{m}_{2}}\displaystyle \mathop{\displaystyle \underbrace{-{\sum }_{k\ne i}{\xi }_{{l}_{2}}(1+{\upsilon }_{{l}_{2}}^{{\rm{L}}}){\hat{{\rm{\Phi }}}}_{k,{l}_{2},{m}_{2}}+{\sum }_{k\gt i}{\hat{{\rm{\Phi }}}}_{k,{l}_{2},{m}_{2}}+{\sum }_{k\lt i}\frac{{\rho }_{k}}{{\rho }_{i}}{\hat{{\rm{\Phi }}}}_{k,{l}_{2},{m}_{2}}}}\limits_{{\rm{s}}{\rm{t}}{\rm{r}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{f}}{\rm{i}}{\rm{c}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}}\right)+\mathop{\underbrace{\left(\frac{{\rho }_{1}}{{\rho }_{i}}-{\xi }_{{l}_{2}}{\upsilon }_{{l}_{2}}^{{\rm{P}}}\right){\hat{q}}_{{l}_{2},{m}_{2}}}}\limits_{{\rm{s}}{\rm{h}}{\rm{e}}{\rm{l}}{\rm{l}}}\right]}}\limits_{{\rm{o}}{\rm{c}}{\rm{e}}{\rm{a}}{\rm{n}}\text{}{\rm{t}}{\rm{h}}{\rm{i}}{\rm{c}}{\rm{k}}{\rm{n}}{\rm{e}}{\rm{s}}{\rm{s}}\text{}{\rm{v}}{\rm{a}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}{\rm{s}}}\\ -\,\displaystyle \mathop{\displaystyle \underbrace{i{{ \mathcal D }}_{l}\frac{{{\hat{c}}^{2}}_{i,0}}{\epsilon \hat{\omega }}\left[-{\sum }_{k\ne i}{\xi }_{l}(1+{\upsilon }_{l}^{{\rm{L}}}){\hat{{\rm{\Phi }}}}_{k,l,m}+{\sum }_{k\gt i}{\hat{{\rm{\Phi }}}}_{k,l,m}+{\sum }_{k\lt i}\frac{{\rho }_{k}}{{\rho }_{i}}{\hat{{\rm{\Phi }}}}_{k,l,m}\right]}}\limits_{{\rm{s}}{\rm{t}}{\rm{r}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{f}}{\rm{i}}{\rm{c}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}}=\frac{{\hat{c}}_{i,0}^{2}}{\epsilon }(1+{\upsilon }_{l}^{{\rm{T}}}){\hat{\eta }}_{l,m}^{{\rm{T}}}+\displaystyle \mathop{\displaystyle \underbrace{{\sum }_{\mathop{{l}_{1},{m}_{1}}\limits_{{l}_{2},{m}_{2}}}{{{ \mathcal T }}^{\left(1\right)}}_{{l}_{1},{m}_{1},{l}_{2},{m}_{2}}^{l,m}\frac{{\hat{c}}_{i,{l}_{1},{m}_{1}}^{2}}{\epsilon }(1+{\upsilon }_{{l}_{2}}^{T}){\hat{\eta }}_{{l}_{2},{m}_{2}}^{{\rm{T}}}}}\limits_{{\rm{o}}{\rm{c}}{\rm{e}}{\rm{a}}{\rm{n}}\text{}{\rm{t}}{\rm{h}}{\rm{i}}{\rm{c}}{\rm{k}}{\rm{n}}{\rm{e}}{\rm{s}}{\rm{s}}\text{}{\rm{v}}{\rm{a}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}{\rm{s}}},\end{array}\end{eqnarray} \tag{ 29a }$$

$$\begin{eqnarray}&&\begin{array}{c}\left[i\hat{\omega }+\hat{\alpha }+{{ \mathcal C }}_{l,m}^{\left(1\right)}\right]{\hat{{\rm{\Psi }}}}_{i,l,m}-{{ \mathcal C }}_{l,m}^{\left(2\right)}{\hat{{\rm{\Phi }}}}_{i,l-1,m}-{{ \mathcal C }}_{l,m}^{\left(3\right)}{\hat{{\rm{\Phi }}}}_{i,l+1,m}\\ -\,\displaystyle \mathop{\displaystyle \underbrace{{\sum }_{\mathop{{l}_{1},{m}_{1}}\limits_{{l}_{2},{m}_{2}}}{{{ \mathcal T }}^{\left(2\right)}}_{{l}_{1},{m}_{1},{l}_{2},{m}_{2}}^{l,m}\frac{{\hat{c}}_{i,{l}_{1},{m}_{1}}^{2}}{\epsilon }\left[-\frac{i{{ \mathcal D }}_{{l}_{2}}}{\hat{\omega }}\left(\left(1-{\xi }_{{l}_{2}}\left(1+{\upsilon }_{{l}_{2}}^{{\rm{L}}}\right)\right){\hat{{\rm{\Phi }}}}_{i,{l}_{2},{m}_{2}}\displaystyle \mathop{\displaystyle \underbrace{-{\sum }_{k\ne i}{\xi }_{{l}_{2}}(1+{\upsilon }_{{l}_{2}}^{{\rm{L}}}){\hat{{\rm{\Phi }}}}_{k,{l}_{2},{m}_{2}}+{\sum }_{k\gt i}{\hat{{\rm{\Phi }}}}_{k,{l}_{2},{m}_{2}}+{\sum }_{k\lt i}\frac{{\rho }_{k}}{{\rho }_{i}}{\hat{{\rm{\Phi }}}}_{k,{l}_{2},{m}_{2}}}}\limits_{{\rm{s}}{\rm{t}}{\rm{r}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{f}}{\rm{i}}{\rm{c}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}}\right)+\mathop{\underbrace{\left(\frac{{\rho }_{1}}{{\rho }_{i}}-{\xi }_{{l}_{2}}{\upsilon }_{{l}_{2}}^{{\rm{P}}}\right){\hat{q}}_{{l}_{2},{m}_{2}}}}\limits_{{\rm{s}}{\rm{h}}{\rm{e}}{\rm{l}}{\rm{l}}}\right]}}\limits_{{\rm{o}}{\rm{c}}{\rm{e}}{\rm{a}}{\rm{n}}\text{}{\rm{t}}{\rm{h}}{\rm{i}}{\rm{c}}{\rm{k}}{\rm{n}}{\rm{e}}{\rm{s}}{\rm{s}}\text{}{\rm{v}}{\rm{a}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}{\rm{s}}}\\ =\,-\displaystyle \mathop{\displaystyle \underbrace{{\sum }_{\mathop{{l}_{1},{m}_{1}}\limits_{{l}_{2},{m}_{2}}}{{{ \mathcal T }}^{\left(2\right)}}_{{l}_{1},{m}_{1},{l}_{2},{m}_{2}}^{l,m}\frac{{\hat{c}}_{i,{l}_{1},{m}_{1}}^{2}}{\epsilon }(1+{\upsilon }_{{l}_{2}}^{T}){\hat{\eta }}_{{l}_{2},{m}_{2}}^{{\rm{T}}}}}\limits_{{\rm{o}}{\rm{c}}{\rm{e}}{\rm{a}}{\rm{n}}\text{}{\rm{t}}{\rm{h}}{\rm{i}}{\rm{c}}{\rm{k}}{\rm{n}}{\rm{e}}{\rm{s}}{\rm{s}}\text{}{\rm{v}}{\rm{a}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}{\rm{s}}};\end{array}\end{eqnarray} \tag{ 29b }$$

$$\begin{eqnarray}&&\begin{array}{c}-{il}(l+1)\left[{\left(l+2\right)}^{2}{\left(l-1\right)}^{2}+(1-\nu )(l+2)(l-1)\right]\frac{{\left(\chi \psi {\hat{D}}^{(b)}\right)}_{0}}{\hat{\omega }}{\sum }_{k}^{N}{\hat{{\rm{\Phi }}}}_{k,l,m}-{\chi }_{0}(l+2)(l-1){\hat{F}}_{l,m}\\ +\,\displaystyle \mathop{\displaystyle \underbrace{{\sum }_{\mathop{{l}_{1},{m}_{1}}\limits_{{l}_{2},{m}_{2}}}\left[\frac{-i}{\hat{\omega }}({l}_{2}+1){l}_{2}\left({{ \mathcal M }}_{{l}_{1},{m}_{1},{l}_{2},{m}_{2}}^{l,m}-(1-\nu ){{ \mathcal N }}_{{l}_{1},{m}_{1},{l}_{2},{m}_{2}}^{l,m}\right){\left(\chi \psi {\hat{D}}^{(b)}\right)}_{{l}_{1},{m}_{1}}{\sum }_{k}^{N}{\hat{{\rm{\Phi }}}}_{k,{l}_{2},{m}_{2}}+{{ \mathcal N }}_{{l}_{1},{m}_{1},{l}_{2},{m}_{2}}^{l,m}{\chi }_{{l}_{1},{m}_{1}}{F}_{{l}_{2},{m}_{2}}\right]}}\limits_{\mathrm{shell}\ \mathrm{thickness}\ \mathrm{variations}}={\gamma }^{-1}{\hat{q}}_{{lm}}\end{array}\end{eqnarray} \tag{ 29c }$$

$$\begin{eqnarray}&&\begin{array}{c}\left[{\left(l+2\right)}^{2}{\left(l-1\right)}^{2}+(1+\nu )(l+2)(l-1)\right]{\left(\chi \psi {\hat{\alpha }}^{(e)}\right)}_{0}{\hat{F}}_{l,m}-i\frac{1}{\hat{\omega }}(l+2)(l+1)l(l-1){\chi }_{0}{\sum }_{k}^{N}{\hat{{\rm{\Phi }}}}_{k,l,m}\\ +\,\displaystyle \mathop{\displaystyle \underbrace{{\sum }_{\mathop{{l}_{1},{m}_{1}}\limits_{{l}_{2},{m}_{2}}}\left[\left({{ \mathcal M }}_{{l}_{1},{m}_{1},{l}_{2},{m}_{2}}^{l,m}-(1+\nu ){{ \mathcal N }}_{{l}_{1},{m}_{1},{l}_{2},{m}_{2}}^{l,m}\right){\left(\chi \psi {\hat{\alpha }}^{(e)}\right)}_{{l}_{1},{m}_{1}}{\hat{F}}_{{l}_{2},{m}_{2}}+\frac{i}{\hat{\omega }}{l}_{2}({l}_{2}+1){{ \mathcal N }}_{{l}_{1},{m}_{1},{l}_{2},{m}_{2}}^{l,m}{\chi }_{{l}_{1},{m}_{1}}{\sum }_{k}^{N}{{\rm{\Phi }}}_{k,{l}_{2},{m}_{2}}\right]}}\limits_{\mathrm{shell}\ \mathrm{thickness}\ \mathrm{variations}}=0.\end{array}\end{eqnarray} \tag{ 29d }$$

The operators ${ \mathcal D }$, ${ \mathcal C }$, ${ \mathcal T }$, ${ \mathcal N }$, and ${ \mathcal M }$ are given in Appendix F, and υ^T, υ^L and υ^P are a set of coefficients related to the deformation of the mantle:

$$\begin{eqnarray}&&{\upsilon }_{l}^{{\rm{T}}}={k}_{l}^{{\rm{T}}}-{h}_{l}^{{\rm{T}}},\,\,{\upsilon }_{l}^{{\rm{L}}}={k}_{l}^{{\rm{L}}}-{h}_{l}^{{\rm{L}}},\,\,{\upsilon }_{l}^{{\rm{P}}}={k}_{l}^{{\rm{P}}}-{h}_{l}^{{\rm{P}}}.\end{eqnarray} \tag{ 30 }$$

The main nondimensional parameters controlling the response of the moon are given in Table 1. $\hat{\omega }$ is the nondimensional forcing frequency. $\hat{\alpha }$ is the inverse of the characteristic ocean damping timescale; a high value indicates an overdamped system. The Lamb parameter  is given by the ratio between the velocity at which a perturbation of frequency Ω and wavelength R propagates and the surface wave speed in a shell-free ocean. When $| \epsilon \hat{\omega }| \ll 1$, the terms on the left-hand side of Equation (8b) are much smaller than those on the right-hand side. Ocean dynamics is then negligible, and the ocean follows the equilibrium tide, which for a nongravitating shell-free body is simply given by ${\hat{\eta }}^{{\rm{T}}}$. For synchronous satellites experiencing tides $\hat{\omega }=1$, and a look at Table 2 shows that the Lamb parameter is small for all bodies, indicating that surface waves do not play a significant role in their response to planet tides. This is not the case for tides raised by other moons in the system, for which the forcing frequency $| \hat{\omega }|$ might be much higher than 1 and hence $| \epsilon \hat{\omega }| \ll 1$ might not hold (Hay et al. 2020, 2022). The solid shell adds two additional nondimensional parameters: γ and ${\hat{D}}_{0}^{(b)}$. γ is the effective rigidity and indicates the importance of the shell on the tidal response. In bodies with a high effective rigidity (hard-shell bodies) the shell hinders ocean surface displacements—e.g., the amplitude of the equilibrium tide is reduced. In contrast, in bodies with low effective rigidity (soft-shell bodies) the shell does not hinder ocean surface displacements—the equilibrium tide is barely affected by the presence of the shell (Beuthe 2018; Matsuyama et al. 2018). From Table 2, we see that small bodies generally behave as hard-shell bodies and bigger ones as soft-shell bodies. ${\hat{D}}_{0}^{(b)}$ is related to the importance of bending effects. The role of bending effects depends on the wavelength (see Equation (B2)); bending effects become more relevant as the wavelength decreases and start to become dominant when $l\gtrsim 1/{({\hat{D}}_{0}^{(b)})}^{1/4}$. ${\hat{D}}_{0}^{(b)}$ is small in all cases, indicating that bending effects are small for the long wavelengths characteristic of tides. ξ encapsulates the importance of ocean loading and ocean self-gravitation, if self-gravitation is ignored, ξ = 0. The amplitude of mantle deformations depends on the mantle's effecftive rigidity ${\hat{\mu }}_{m}$. As shown in Appendix D, the Love numbers of the mantle depend on this parameter as $\sim 1/{\hat{\mu }}_{m}$. The effective rigidity ${\hat{\mu }}_{m}$ is generally large for the bodies considered here, which implies that the Love numbers and hence υ are small. Because of this, in what follows we will neglect mantle deformations.

Table 1. List of the Main Nondimensional Parameters Governing Tidal Dynamics

Uniform Density and Thickness
Name	Symbol	Definition	Meaning
Nondimensional frequency	$\hat{\omega }$	$\tfrac{\omega }{{\rm{\Omega }}}$	Ratio of forcing to rotating frequencies
Nondimensional damping frequency	$\hat{\alpha }$	$\tfrac{\alpha }{{\rm{\Omega }}}$	Relevance of ocean dissipation processes
Ocean loading coefficient	ξ	$\tfrac{{\rho }_{1}}{{\rho }_{b}}$	Relevance of ocean loading and self-gravity
Lamb parameter		$\tfrac{{{\rm{\Omega }}}^{2}{R}^{2}}{{{gH}}_{\mathrm{ocean}}}$	Relevance of surface gravity waves
Effective rigidity	γ	$\tfrac{2{\mu }_{0}(1+\nu ){H}_{\mathrm{shell}}}{g{\rho }_{1}{R}^{2}}$	Relevance of the solid shell
Effective bending rigidity	${\hat{D}}_{0}^{(b)}$	$\tfrac{1}{12(1-{\nu }^{2})}{\left(\tfrac{{H}_{\mathrm{shell}}}{R}\right)}^{2}$	Relevance of bending effects
Mantle effective rigidity	${\hat{\mu }}_{m}$	$\tfrac{{\mu }_{m}}{{\rho }_{m}{g}_{m}{R}_{m}}$	Role of mantle deformation
Variable Density
Relative layer thickness	${\hat{c}}_{i}^{2}$	$\tfrac{{H}_{i}}{{H}_{\mathrm{ocean}}}$	Layer-to-total ocean thickness ratio
Density ratio	ri	$\tfrac{{\rho }_{1}}{{\rho }_{i}}$	Density contrast within the ocean
Variable Thickness
Relative ocean thickness variations	${\hat{c}}_{l,m}^{2}$	$\tfrac{{H}_{l,m}}{{H}_{\mathrm{ocean},0}}$	Amplitude of ocean thickness variations
Relative ice thickness variations	${\hat{H}}_{\mathrm{shell},l,m}$	$\tfrac{{H}_{\mathrm{shell},l,m}}{{H}_{\mathrm{shell},0}}$	Amplitude of shell thickness variations

Download table as:  ASCIITypeset image

Table 2. Physical Parameters of a Selection of Tidally Active Bodies of the Solar System

	Moon	Io	Europa	Ganymede	Callisto	Mimas	Enceladus	Dione	Titan	Triton
Radius (km)	1737.4	1821.5	1560.8	2631.2	2410.3	198.2	252.1	561.7	2574.7	1353.4
Mass (1022 kg)	7.3	8.9	4.8	14.8	10.8	0.0037	0.011	0.11	13.5	2.1
Average density (103 kg m−3)	3.34	3.53	3.01	1.94	1.84	1.15	1.61	1.48	1.88	2.06
Surface gravity (m s−2)	1.62	1.80	1.32	1.43	1.25	0.064	0.11	0.23	1.35	0.78
Shell thickness (km)	50	30	30	125	115	28	23	99	100	150
Ocean thickness (km)	50	50	100	150	80	44	38	65	255	150
Ocean thickness variations (km)	⋯	⋯	7	⋯	⋯	⋯	22	⋯	12	⋯
Tidal period (days)	27.32	1.77	3.55	7.16	16.69	0.94	1.37	2.74	15.95	5.88
Eccentricity (10−3)	55.4	4.1	9.4	1.3	7.4	19.6	4.7	2.2	28.8	⋯
Inclination (deg)	5.16	0.036	0.466	0.177	0.192	1.57	0.003	0.028	0.306	156.87
Obliquity (deg)	6.68	0.0021	0.053	0.033	0.24	0.041	0.00045	0.0020	0.32	0.35
Lamb parameter ()	7.2 × 10−3	6.2 × 10−2	7.8 × 10−3	3.3 × 10−3	1.1 × 10−3	8.3 × 10−2	4.2 × 10−2	1.5 × 10−2	4.0 × 10−4	2.4 × 10−3
Effective rigidity (γ)	6.2 × 10−1	3.1 × 10−1	8.2 × 10−2	1.1 × 10−1	1.4 × 10−1	9.67 × 10	2.8 × 10	1.2 × 10	9.8 × 10−2	9.2 × 10−1
Bending (${\hat{D}}^{(b)}$)	7.7 × 10−5	2.5 × 10−5	3.45 × 10−5	2.1 × 10−4	2.1 × 10−4	1.8×10−3	7.8 × 10−4	2.9 × 10−3	1.4 × 10−4	1.1 × 10−3
Amp. Equation eccentricity tide (m)	2.3	40.9	23.3	2.2	2.1	229	23.7	6.8	9.4	0
Amp. Equation obliquity tide (m)	2.21	0.17	1.05	0.04	0.54	3.83	0.018	0.05	0.83	3.2
Surface heat flux (W m−2)	⋯	2.24	0.02	⋯	⋯	⋯	0.02–0.04	⋯	0.005	⋯

Note. Mass, radius, surface gravity, orbital period, eccentricity, and inclination are from NASA/JPL's Solar System Dynamics database (https://ssd.jpl.nasa.gov). Only the obliquity of the Moon and Titan correspond to measured values (Stiles et al. 2008; Lang 2011). For the remaining satellites the obliquity is computed by assuming that the satellite is in a Cassini state; see Chen et al. (2014) and Baland et al. (2016) for the icy moons and Baland et al. (2012) for Io. Representative shell and ocean thickness values: Io (Jaeger et al. 2003; Khurana et al. 2011); Europa (Anderson et al. 1998; Hussmann et al. 2002); Ganymede (Kivelson et al. 2002; Vance et al. 2018); Callisto and Titan (Vance et al. 2018); Enceladus and Dione (Beuthe et al. 2016); Mimas (Rhoden & Walker 2022); Triton (Nimmo & Spencer 2015). For the Moon, ocean thickness, shell thickness, and Lamb parameter are representative of values attained during the freezing process (Chen & Nimmo 2016). Ocean thickness variations for Europa, Enceladus, and Titan are from Nimmo et al. (2007), Beuthe et al. (2016), and Lefevre et al. (2014). Parameters , γ, and ${\hat{D}}^{(b)}$ are obtained assuming a Young modulus of 8.8 and 170 GPa for ice and rock, respectively, and a Poisson ratioof ν = 0.33. For Io the heat flux is estimated from astrometric measurements (Lainey et al. 2009); for Enceladus, it is obtained from Cassini's infrared observations of the south pole (Howett et al. 2011) and heat flux estimates obtained assuming that the moon is in thermal equilibrium (Hemingway et al. 2018); for Europa, it is estimated assuming that Europa is in thermal equilibrium (Hussmann et al. 2002); and for Titan, it is inferred from topography (Nimmo & Bills 2010).

Download table as:  ASCIITypeset image

Lateral ocean or shell thickness variations add new nondimensional parameters: ${\hat{c}}_{l,m}^{2}$, which depends on ocean thickness variations, and ${(\chi \psi {\hat{\alpha }}^{(e)})}_{l,m}$ and $(\chi \psi \hat{D}{)}_{l,m}^{(b)}$, which follow from shell thickness variations. Moreover, each ocean density interface adds two new nondimensional parameters: the density ratio between the two layers r, and the ratio between the thickness of the new two layers resulting from the new density interface, ${\hat{c}}_{i}^{2}$. As we will see, these two parameters control the excitation of internal waves.

We can get some insight into the problem by examining Equation (29d). The extra terms that arise from thickness variations and stratification are indicated. The terms with ${ \mathcal D }$ are related to gravity waves in a shell-free ocean. Terms proportional to q arise from the presence of a shell, which, as we will later see, effectively increases the surface gravity wave speed. Terms with ${ \mathcal C }$ are the result of the Coriolis force. Without these terms, equations of different degree l would be decoupled. The terms containing ${ \mathcal T }$ are the result of ocean thickness variations. They couple equations of different degree and order even in the absence of rotation. Similarly, the terms ${ \mathcal M }$ and ${ \mathcal N }$ arise as a result of lateral variations in shell properties. If there are only meridional thickness variations m₁ = 0, equations with different orders are uncoupled as ${{ \mathcal T }}_{{l}_{1},0,{l}_{2},{m}_{2}}^{l,m}={{ \mathcal M }}_{{l}_{1},0,{l}_{2},{m}_{2}}^{l,m}={{ \mathcal N }}_{{l}_{1},0,{l}_{2},{m}_{2}}^{l,m}=0$ if m ≠ m₂. If the ocean is unstratified (N = 1), it is of uniform thickness (${\hat{c}}_{l,m}^{2}=0$ for l ≠ 0), and we ignore the presence of a solid shell (γ = 0), we recover the classic LTEs used by Longuet-Higgins (1968) to study the response of a terrestrial global ocean and first employed by Tyler (2008, 2009) to study subsurface oceans. If we consider a spherically symmetric solid shell, we obtain the equation used in Beuthe (2016) (see Appendix B).

Equation (29d) constitutes a linear system of 2(1 + N)(L + 1)² equations of the form

$$\begin{eqnarray}&&{\bf{A}}{\boldsymbol{x}}={\boldsymbol{f}}.\end{eqnarray} \tag{ 31 }$$

N is the number of layers, and L is the degree at which we cut the system of equations. A follows from Equations (29d), x is the vector containing the unknowns

$$\begin{eqnarray}\begin{array}{ccc}{\boldsymbol{x}} & = & \left[{\hat{{\rm{\Phi }}}}_{i,l,m}\quad {\hat{{\rm{\Psi }}}}_{i,l,m}\quad {\hat{q}}_{l,m}\quad {\hat{F}}_{l,m}...\quad {{\rm{\Phi }}}_{N,L,L}\quad {\hat{{\rm{\Psi }}}}_{N,L,L}\right.\\ & & \,{\left.{\hat{q}}_{L,L}\quad {\hat{F}}_{L,L}\right]}^{\dagger },\end{array}\end{eqnarray} \tag{ 32 }$$

and f is the forcing term

$$\begin{eqnarray}&&{\boldsymbol{f}}={\left[{f}_{i,l,m,{\rm{\Phi }}}\,{f}_{i,l,m,{\rm{\Psi }}}\,0\,0\,...\,{f}_{N,L,L,{\rm{\Phi }}}\,{f}_{N,L,L,{\rm{\Psi }}}\,0\,0\right]}^{\dagger },\end{eqnarray} \tag{ 33 }$$

with f_i,l,m,Φ and f_i,l,m,Ψ being the forcing terms on the right-hand side of Equations (29a) and (29b). The superscript † indicates the transpose.

If we only consider a forcing of a given degree and order and restrict the problem to meridional ocean and shell thickness variations, the system reduces to 2(1 + N)(1 + L) equations. Moreover, if the shell is spherically symmetric, q and F can be written in terms of Φ and the system reduces to 2N(1 + L) (see Appendix B).

The approach presented above allows us to obtain the tidal response due to a general tidal perturbation. Here we will focus on the tides experienced by a synchronous satellite of nonzero eccentricity e and obliquity θ_i . In this case $| \hat{\omega }| =1$. The tidal potential has two components, the eccentricity (ϕ_e ) and the obliquity tide (ϕ_o ),

$$\begin{eqnarray}&&{\hat{\phi }}^{{\rm{T}}}={\hat{\phi }}_{e}+{\hat{\phi }}_{o}+O({e}^{2})+O({\sin }^{2}{\theta }_{i}).\end{eqnarray} \tag{ 34 }$$

The tidal response of the ocean depends on the direction in which the tidal potential propagates. Because of this, the tidal potential is further split into westward- and eastward-propagating components (e.g., Chen et al. 2014; Beuthe et al. 2016),

$$\begin{eqnarray}\begin{array}{rcl}{\hat{\phi }}_{e} & = & {\hat{\phi }}_{2,0}^{W}{Y}_{2}^{0}{e}^{i\hat{t}}+{\hat{\phi }}_{2,0}^{E}{Y}_{2}^{0}{e}^{-i\hat{t}}+{\hat{\phi }}_{2,2}^{W}{Y}_{2}^{2}{e}^{i\hat{t}}\\ & & +{\hat{\phi }}_{2,2}^{E}{Y}_{2}^{2}{e}^{-i\hat{t}}+{\rm{c.c}}\end{array}\end{eqnarray} \tag{ 35a }$$

$$\begin{eqnarray}&&{\hat{\phi }}_{o}={\hat{\phi }}_{2,1}^{W}{Y}_{2}^{1}{e}^{i\hat{t}}+{\hat{\phi }}_{2,1}^{E}{Y}_{2}^{1}{e}^{-i\hat{t}}+{\rm{c.c}}.\end{eqnarray} \tag{ 35b }$$

The W and E superscripts indicate the westward and eastward components, respectively. The zonal component is a standing wave—it does not propagate in the longitudinal direction—however, splitting it into $+\hat{\omega }$ and $-\hat{\omega }$ components is convenient. The amplitude of each component of the tidal potential is given in Table 3. The solution is similarly expanded into westward- and eastward-propagating components

$$\begin{eqnarray}&&{\hat{{\rm{\Phi }}}}_{i}^{W}=\displaystyle \sum _{l,m}\left[{\hat{{\rm{\Phi }}}}_{i,l,m}^{W}{Y}_{l}^{m}{e}^{i\hat{t}}+{\hat{{\rm{\Phi }}}}_{i,l,-m}^{W}{Y}_{l}^{-m}{e}^{-i\hat{t}}\right],\end{eqnarray} \tag{ 36a }$$

$$\begin{eqnarray}&&{\hat{{\rm{\Phi }}}}_{i}^{E}=\displaystyle \sum _{l,m}\left[{\hat{{\rm{\Phi }}}}_{i,l,m}^{E}{Y}_{l}^{m}{e}^{-i\hat{t}}+{\hat{{\rm{\Phi }}}}_{i,l,-m}^{E}{Y}_{l}^{-m}{e}^{+i\hat{t}}\right],\end{eqnarray} \tag{ 36b }$$

and the same for the remaining variables—${\hat{{\rm{\Psi }}}}_{i}$, $\hat{F}$, and $\hat{q}$. Since the different variables should be real, it follows that ${\hat{{\rm{\Phi }}}}_{i,l,m}={\left(-\right)}^{m}{\hat{{\rm{\Phi }}}}_{i,l,-m}^{* }$, and similarly for the other fields.

Table 3. Amplitude of the Components of the Eccentricity and Obliquity Tides

Tidal Component	Amplitude
${\hat{\phi }}_{2,0}^{W}$	$\tfrac{3}{8}\sqrt{\tfrac{1}{5}}e$
${\hat{\phi }}_{2,0}^{E}$	$\tfrac{3}{8}\sqrt{\tfrac{1}{5}}e$
${\hat{\phi }}_{2,2}^{W}$	$\tfrac{1}{8}\sqrt{\tfrac{6}{5}}e$
${\hat{\phi }}_{2,2}^{E}$	$-\tfrac{7}{8}\sqrt{\tfrac{6}{5}}e$
${\hat{\phi }}_{2,1}^{W}$	$\tfrac{1}{4}\sqrt{\tfrac{6}{5}}\sin {\theta }_{i}$
${\hat{\phi }}_{2,1}^{E}$	$\tfrac{1}{4}\sqrt{\tfrac{6}{5}}\sin {\theta }_{i}$

Download table as:  ASCIITypeset image

Because the equations of motion have been linearized, we can compute the tidal response of the ocean for each component of the tide considering a unit forcing of ${\hat{\eta }}_{2,m}^{T}=1$ by solving Equation (29d) and then combining the solutions using the weights in Table 3. However, tidal dissipation is not linear (Equations (9) and (10)), and hence it should be generally computed once all the components of the tidal response have been obtained. Using the definition of s (Equation (3)), plugging Equations (35) and (36) into Equation (10), and using the orthogonality properties of vector spherical harmonics (Equation (G16b)), we find the total nondimensional energy dissipation

$$\begin{eqnarray}\begin{array}{rcl}\hat{\dot{E}} & = & -48\pi \displaystyle \sum _{i}^{N}\displaystyle \frac{{\rho }_{i}}{{\rho }_{1}}\left[{\hat{\phi }}_{2,0}^{W}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{i,2,0}^{W})+{\hat{\phi }}_{2,0}^{E}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{i,2,0}^{E})\right.\\ & & +\,{\hat{\phi }}_{2,2}^{W}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{i,2,2}^{W})+{\hat{\phi }}_{2,2}^{E}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{i,2,2}^{E})\\ & & +\left.\,{\hat{\phi }}_{2,1}^{W}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{i,2,1}^{W})+{\hat{\phi }}_{2,1}^{E}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{i,2,1}^{E})\right].\end{array}\end{eqnarray} \tag{ 37 }$$

The dimensional tidal dissipation is given by

$$\begin{eqnarray}&&\dot{E}=\rho {{\rm{\Omega }}}^{3}{R}^{5}\left(\displaystyle \frac{R{{\rm{\Omega }}}^{2}}{g}\right)\hat{\dot{E}},\end{eqnarray} \tag{ 38 }$$

where the nondimensional term in parentheses accounts for the fact that while we force the problem with ${\hat{\eta }}^{T}=1$, the magnitude of ${\hat{\eta }}^{T}$ experienced by each satellite depends on its orbital frequency, radius, and surface gravity.

As noted above, examining Equation (37), we see that we cannot always separate the contribution of the different components of the tide to the energy budget. For example, if the eccentricity tide produces an m = 1 ocean tidal response (${\hat{{\rm{\Phi }}}}_{2,1}^{X})$, a term arising from the combination of this term with the obliquity tide (${\hat{\phi }}_{2,1}^{X}$) of amplitude $\propto e\sin {\theta }_{i}$ appears. Nevertheless, in some special cases the contribution of tidal dissipation due to each component of the tide can be separated. If there are not lateral variations with odd order m, the contribution of eccentricity and obliquity tides can be separated. Similarly, if there are not lateral variations with even order m, the contribution of the m = 0 and m = 2 components can be separated. In the special case of only zonal variations, m = 0, the contribution of all tidal components to tidal dissipation can be separated. Below we will focus on this case and distinguish between tidal dissipation due to the eccentricity and the obliquity tide. Furthermore, to make the nondimensional tidal response independent of the eccentricity and obliquity of a specific body, we will obtain the tidal response for e = 1 and $\sin {\theta }_{i}=1$, which can be easily scaled to a particular body by multiplying by its eccentricity and obliquity.

The tidal response of unstratified oceans of constant thickness and density has been widely studied in the past decade (e.g., Tyler 2008, 2009, 2011; Chen et al. 2014; Matsuyama 2014; Beuthe 2016; Hay & Matsuyama 2017, 2019; Matsuyama et al. 2018). If the ocean is unstratified and both ocean and solid shell are of constant thickness, Equation (29d) reduces to

$$\begin{eqnarray}&&\begin{array}{c}\left(i\hat{\omega }+\hat{\alpha }+{{ \mathcal C }}_{l,m}^{\left(1\right)}-\displaystyle \frac{i{{ \mathcal D }}_{l}}{\epsilon \hat{\omega }}\left(1-{\xi }_{l}+{{\rm{\Lambda }}}_{l}\right)\right)\\ {\hat{{\rm{\Phi }}}}_{l,m}\,+{{ \mathcal C }}_{l,m}^{\left(2\right)}{\hat{{\rm{\Psi }}}}_{l-1,m}+{{ \mathcal C }}_{l,m}^{\left(3\right)}{\hat{{\rm{\Psi }}}}_{l+1,m}=\displaystyle \frac{{\hat{\eta }}_{l,m}^{{\rm{T}}}}{\epsilon }\end{array}\end{eqnarray} \tag{ 39a }$$

$$\begin{eqnarray}&&\left(i\hat{\omega }+\hat{\alpha }+{{ \mathcal C }}_{l,m}^{(1)}\right){\hat{{\rm{\Psi }}}}_{l,m}-{{ \mathcal C }}_{l,m}^{(2)}{\hat{{\rm{\Phi }}}}_{l-1,m}-{{ \mathcal C }}_{l,m}^{(3)}{\hat{{\rm{\Phi }}}}_{l+1,m}=0,\end{eqnarray} \tag{ 39b }$$

where Λ_l is a nondimensional shell spring constant proportional to the effective rigidity γ (Beuthe 2015a; see Appendix B).

Equation (39) can be cast into a matrix form and the solution obtained by inverting A. When A is close to singular, a resonance can occur. The frequencies $\hat{\omega }$ for which this happens are the eigenfrequencies of the system. The eigenfrequencies cannot be analytically obtained except if the moon is not rotating or the Coriolis terms are neglected (${ \mathcal C }=0$), in which case A becomes diagonal. Note that in the nonrotating limit the rotational frequency of the body Ω is 0; instead of using the rotational frequency to nondimensionalize the equations of motion, we use the orbital frequency, ∣ω∣ (Equation (26)). Provided that the damping is small $\hat{\alpha }\ll 1$, the eigenfrequencies are

$$\begin{eqnarray}&&{\hat{\omega }}_{l}^{2}=\displaystyle \frac{l(l+1)}{\epsilon }\left(1+{{\rm{\Lambda }}}_{l}-{\xi }_{l}\right).\end{eqnarray} \tag{ 40 }$$

Using the previous expression, we can obtain the values of , for which the ocean resonates at the tidal frequency (${\hat{\omega }}_{l}=\pm 1$),

$$\begin{eqnarray}&&{\epsilon }_{l}^{\mathrm{res}}=l(l+1)(1+{{\rm{\Lambda }}}_{l}-{\xi }_{l}),\end{eqnarray} \tag{ 41 }$$

which is equivalent to the expressions given in Beuthe (2015b) and Beuthe et al. (2016). These resonances correspond to surface gravity waves. Examining Equation (41), it becomes evident that the resonant increases with increasing degree or, equivalently, the resonant ocean thickness decreases with increasing degree. Moreover, we see that the presence of an elastic shell shifts the resonant ocean thickness toward lower values. From the values of and γ characteristic of the moons of the outer solar system (Table 2), we see that for all of them $\epsilon \ll {\epsilon }_{2}^{\mathrm{res}}$, which implies that planet tides do not excite resonant surface gravity waves. Using $\epsilon \ll {\epsilon }_{2}^{\mathrm{res}}$ and Equation (27), we find that in such cases the ocean surface displacement is well approximated by

$$\begin{eqnarray}&&{\hat{\eta }}_{l,m}\approx \displaystyle \frac{{\hat{\eta }}_{l,m}^{{\rm{T}}}}{1+{{\rm{\Lambda }}}_{l}-{\xi }_{l}};\end{eqnarray} \tag{ 42 }$$

the ocean surface follows the equilibrium tide.

As discussed above, the inclusion of rotation couples equations of different degrees l, changing the resonant frequencies of the system. Moreover, even though the tidal forcing is restricted to degree l = 2, modes containing shorter wavelengths can be excited, and, as opposed to the nonrotating case, the eigenfrequencies of the system depend on the order m and on whether the tidal wave propagates westward or eastward. Nevertheless, equatorial symmetric and antisymmetric modes remain decoupled. An equatorial symmetric forcing (e.g., ${\hat{\eta }}_{2,2}^{{\rm{T}}}$) results in a symmetric flow field (${{\boldsymbol{x}}}^{\mathrm{sym}}={[{\hat{{\rm{\Psi }}}}_{\mathrm{1,2}},{\hat{{\rm{\Phi }}}}_{\mathrm{2,2}},{\hat{{\rm{\Psi }}}}_{\mathrm{3,2}},{\hat{{\rm{\Phi }}}}_{\mathrm{4,2}}...]}^{\dagger }$), while an antisymmetric forcing (e.g., ${\hat{\eta }}_{2,1}^{{\rm{T}}}$) results in an antisymmetric flow field (${{\boldsymbol{x}}}^{\mathrm{asym}}={[{\hat{{\rm{\Psi }}}}_{\mathrm{1,1}},{\hat{{\rm{\Phi }}}}_{\mathrm{2,1}},{\hat{{\rm{\Psi }}}}_{\mathrm{3,1}},{\hat{{\rm{\Phi }}}}_{\mathrm{4,1}}...]}^{\dagger }$).

The inclusion of rotation also introduces Rossby waves into the system. Tyler (2008) showed that the westward component of the obliquity tide ($l=2,m=1,\hat{\omega }=1$) can excite strong Rossby waves in thick subsurface oceans. The Rossby wave solution can be well approximated by cutting off Equation (39) at degree 2. Doing so, an analytical solution for Rossby waves can be obtained (Beuthe 2016),

$$\begin{eqnarray}&&{\hat{{\rm{\Phi }}}}_{2,1}^{W}\approx \displaystyle \frac{\sqrt{5}}{3}\hat{\alpha }{\hat{{\rm{\Psi }}}}_{1,1}^{W},\end{eqnarray} \tag{ 43a }$$

$$\begin{eqnarray}&&{\hat{{\rm{\Psi }}}}_{1,1}^{W}\approx \displaystyle \frac{\sqrt{5}}{1+\tfrac{5}{3}\hat{\alpha }\left(\hat{\alpha }+i\left(\tfrac{2}{3}-\tfrac{6}{\epsilon }\left(1+{{\rm{\Lambda }}}_{2}-{\xi }_{2}\right)\right)\right)}\displaystyle \frac{1}{\epsilon }{\hat{\eta }}_{2,1}^{{\rm{T}}}.\end{eqnarray} \tag{ 43b }$$

If there is no energy dissipation ($\hat{\alpha }=0$), the solution only contains the term ${\hat{{\rm{\Psi }}}}_{2,1}^{W}$ and is equivalent to the solution derived in Tyler (2008). Moreover, as noted in Beuthe (2016), if $\hat{\alpha }\ll 1$ and $10\hat{\alpha }(1+{{\rm{\Lambda }}}_{2}-{\xi }_{2})/\epsilon \ll 1$, the flow velocity, ∇ × (Ψ_1,1 1_r /H), is independent of ocean thickness and tidal dissipation increases linearly with ocean thickness.

The main characteristics of the tidal response of an unstratified ocean of constant thickness can be observed in Figure 2. For both eccentricity and obliquity tides, surface gravity waves occur at various values of Lamb parameter . Nevertheless, this occurs far from the region characteristic of the moons considered. Additionally, the excitation of Rossby waves is evident by an increase of tidal dissipation with decreasing for the obliquity tide. As expected from Equations (43), the amount of tidal dissipation due to Rossby waves depends on $\hat{\alpha }$ (see Figure 2(d)). The lower $\hat{\alpha }$ is, the thicker the ocean for which tidal dissipation due to obliquity attains its maximum. This dependence disappears if tidal dissipation is modeled using nonlinear bottom drag rather than a Rayleigh dissipation coefficient, in which case dissipation becomes independent of ocean thickness (Hay & Matsuyama 2017, 2019). As opposed to surface gravity waves, Rossby waves can result in a significant amount of tidal dissipation for satellites with a high obliquity, most notably Triton (see Figure 2(e)), provided that the assumption of constant ocean thickness holds.

Zoom In Zoom Out Reset image size

In this section we consider an ocean of uniform density but variable thickness covered by a shell of either constant or variable thickness. Rovira-Navarro et al. (2020) already investigated the tidal response of a global ocean of variable thickness using finite element modeling (FEM). However, the semianalytical method presented here allows a better exploration of the parameter space. In Appendix H, we show that, as expected, both methods lead to the same results.

We will first discuss how ocean thickness variations alter the eigenfrequencies of the ocean (Section 3.2.1). Afterward, we will discuss the important role of thickness variations for the excitation of Rossby waves and the damping of a moon's inclination (Section 3.2.2). Finally, we will consider the specific case of Enceladus in more detail (Section 3.2.3).

3.2.1. The Eigenmodes and Eigenfrequencies of an Ocean of Variable Thickness

As opposed to the case of an ocean and shell of constant thickness, if the ocean and shell are not spherically symmetric, the terms containing ${ \mathcal T }$, ${ \mathcal M }$, and ${ \mathcal N }$ in Equation (29d) are not zero. The terms arising as a result of ocean thickness variations resemble the Coriolis term; they couple equations of different degrees. Moreover, while the Coriolis term does not couple equations of different orders, sectoral and tesseral thickness variations (m_t ≠ 0) couple equations of different orders. An ocean with thickness variations of order m_t experiencing a forcing of order m_f responds with modes of orders m_f + nm_t , with n being an integer, some of which might resonate with the tidal force.

Even though sectoral and tesseral shell or ocean thickness variations are possible in subsurface oceans, we focus on meridional variations (m_t = 0, l_t ≠ 0). This reduces the system of equations from 2(1 + L)² to 2(1 + L). Contrary to the case of an ocean of uniform thickness, equatorially symmetric ( x ^sym) and antisymmetric ocean modes ( x ^asym) are not necessarily decoupled. Nevertheless, if ocean thickness variations are symmetric with respect to the equator, ${{{ \mathcal T }}^{(1)}}_{{l}_{t},0,{l}_{2},{m}_{2}}^{l,m}=0$ when l and l₂ have different parity and ${{{ \mathcal T }}^{(2)}}_{{l}_{t},0,{l}_{2},{m}_{2}}^{l,m}=0$ when both l and l₂ have the same parity. This implies that if the ocean has symmetric ocean thickness variations, symmetric and antisymmetric modes are decoupled, as was the case for a constant-thickness ocean. In contrast, antisymmetric ocean thickness variations couple symmetric and antisymmetric modes, making it possible for a symmetric forcing to excite an antisymmetric ocean response and vice versa.

To illustrate the effects of ocean thickness, we consider long-wavelength ocean thickness variations of degree 2 and 3 (${\hat{c}}_{2,0}^{2}\ne 0$ or ${\hat{c}}_{3,0}^{2}\ne 0$). An ocean with ${\hat{c}}_{2,0}^{2}\gt 0$ is symmetric with respect to the equator and thicker at the poles than at the equator, and one with ${\hat{c}}_{3,0}^{2}\gt 0$ is antisymmetric with respect to the equator and thicker at the north pole than at the south pole.

The effects of ocean thickness variations in the tidal response of a subsurface ocean are illustrated in Figure 3. As in Rovira-Navarro et al. (2020), we observe that degree 2 ocean thickness variations decrease ${\epsilon }^{\mathrm{res}}$. Degree 3 ocean thickness variations do not alter the resonant ocean thickness for modes already excited in a constant-thickness ocean but introduce new modes, which appear as new bands of intense dissipation in Figures 3(c) and (d). These are antisymmetric modes excited by the symmetric eccentricity tide and symmetric modes excited by the antisymmetric obliquity tide.

Zoom In Zoom Out Reset image size

The most dramatic effect is the change of tidal dissipation due to obliquity tides for low values of Lamb parameter (see Figures 3(b) and (d)). Ocean thickness variations inhibit the excitation of Rossby waves. This is evident as a decrease of energy dissipation for thick oceans ( ≪ 1) as ocean thickness variations become more prominent in Figures 3(b) and (d). For positive degree 2 topography (${\hat{c}}_{2,0}^{2}\gt 0$) Rossby waves are not resonant with the tidal frequency; on the other hand, degree 3 topography shifts Rossby wave resonances to thinner oceans (higher ).

When ocean thickness variations are included, the Rossby wave solution given by Equation (43) does not approximate the tidal response of the ocean anymore. We can still find an approximated solution cutting Equation (B3b) at degree 2, but it is less transparent than that for oceans of constant thickness (see Appendix C). As ${\hat{c}}_{2,0}^{2}$ increases, ${\mathfrak{R}}({{\rm{\Phi }}}_{\mathrm{2,1}})$ decreases, which leads to a decrease in energy dissipation (Equation (37)) for the westward obliquity tide. As noted in Rovira-Navarro et al. (2020), this can be understood as ocean thickness variations detuning the ocean from the excitation of Rossby waves. This is illustrated in Figure 4, where tidal dissipation due to a westward forcing is shown as a function of ocean thickness variation and forcing frequency. When ocean topography is added, the maximum amount of energy dissipation due to Rossby waves no longer occurs at the tidal frequency.

Zoom In Zoom Out Reset image size

3.2.2. Consequences for the Thermal−Orbital Evolution of Satellites

We have seen that zonal degree 2 and 3 ocean thickness variations hinder the excitation of Rossby waves by the obliquity tide. It has been suggested that Rossby waves play a prominent role in the thermal budget and orbital evolution of some moons of the outer solar system (Tyler 2008, 2009; Chen et al. 2014; Nimmo & Spencer 2015; Downey et al. 2020), as well as during the early evolution of the Moon (Chen & Nimmo 2016). However, the amount of tidal dissipation resulting from obliquity tides shown in Figure 2(c) can be reduced by several orders of magnitude if there are ocean thickness variations (Figures 3(b) and (d)). Are there any grounds to expect long-wavelength ocean thickness variations in subsurface oceans?

Ocean thickness variations can arise as a result of internal dynamic processes. A clear example is tidal dissipation. Tidal dissipation in the solid layers of a body (i.e., mantle and shell) is not spherically symmetric; the tidal dissipation pattern characteristic of both mantle and shell is mainly degree 2 (Beuthe 2013). We can suspect that such patterns translate into shell thickness variations of the same degree. Viscous flow within the shell tends to homogenize the shell thickness, but shell−ocean interactions can preserve shell thickness variations (Ashkenazy et al. 2018; Čadek et al. 2019; Lobo et al. 2021).

Except for Titan, the obliquities of the outer solar system moons have not been measured. Assuming that the moons are in a Cassini state, Chen et al. (2014) computed the obliquity of candidate ocean worlds and showed that, except for Triton, obliquity ocean tides are unlikely to be a major contributor to the moons' thermal budget. Nimmo & Spencer (2015) showed that ocean tidal dissipation due to obliquity tides can maintain a subsurface ocean in Triton. Our results demonstrate that this is only the case if the ocean is of uniform thickness.

A decrease in obliquity tidal dissipation also has consequences for the orbital evolution. Obliquity ocean tides dampen the moons' inclination (Chyba et al. 1989):

$$\begin{eqnarray}&&\tan i\displaystyle \frac{{di}}{{dt}}\approx -\displaystyle \frac{a}{{GMm}}{\dot{E}}_{\mathrm{obl}},\end{eqnarray} \tag{ 44 }$$

where m and M are the mass of the moon and the planet, respectively, and a is the moon's semimajor axis. If the moon's obliquity is small and the inclination is close to 0 or π, the inclination damping timescale is given by

$$\begin{eqnarray}&&{\tau }_{i}\approx \displaystyle \frac{{GMmg}}{a{\rho }_{1}{R}^{6}{{\rm{\Omega }}}^{5}}{\hat{\dot{E}}}_{\mathrm{obl}}^{-1}\left\{\begin{array}{l}{\left(\displaystyle \frac{i}{{\theta }_{i}}\right)}^{2}\ \mathrm{if}\ i\lt \pi /2,\\ {\left(\displaystyle \frac{i-\pi }{{\theta }_{i}}\right)}^{2}\ \mathrm{if}\ i\gt \pi /2,\end{array}\right.\end{eqnarray} \tag{ 45 }$$

where we have used Equation (38). As illustrated by the previous formula, the inclination timescale depends on the physical and orbital parameters of the moon and the inverse of the nondimensional obliquity tidal dissipation.

Figure 5(a) shows the inclination damping timescale for several moons assuming a constant-thickness ocean for different values of ocean dissipation parameter $\hat{\alpha }$. Except for the Moon, where we use a semimajor axis of 20 Earth radii, the damping timescales are computed for the present semimajor axis of the satellites. This provides an upper bound of the damping timescale throughout the evolution of the moon, as for a synchronous satellite τ_i ∝ a^13/2. Except for Ganymede, Triton, and Dione, damping timescales below the age of the solar system can be attained for a wide range of ocean dissipation properties. Moons with an inclination damping timescale smaller than the age of the solar system are expected to have a small inclination, or else their inclinations have been recently excited or they only had short-lived subsurface oceans. However, ocean thickness variations provide a way to reconcile a high orbital inclination of primordial origin and with long-lived subsurface oceans.

Zoom In Zoom Out Reset image size

The case of the Moon is especially interesting. Chen & Nimmo (2016) showed that obliquity tides in a Lunar magma ocean can quickly dampen the Moon's obliquity, making it difficult to reconcile the present Moon's obliquity with the existence of a magma ocean for more than 10 Myr. This contrasts with geochronological constraints that predict life spans up to 200 Myr (see Chen & Nimmo 2016, and references therein). A possible answer to this puzzle can come from ocean thickness variations. Departures from a constant-thickness ocean can increase the inclination damping timescale by several orders of magnitude; for a 50 km thick magma ocean, a thickness difference between pole and equator of just 3 km—of the same magnitude as the Moon's polar flattening—increases the inclination damping timescale from ∼1 to ∼100 Myr (Figure 5(b)). But can we expect the Moon's magma ocean to present such thickness variations?

The Lunar magma ocean solidified in two distinct phases. During the first phase, the magma reached the lunar surface and crystallization occurred from the bottom up. In the second phase, a floating crust formed, slowing down the cooling process. During the first stage, 80% of the magma ocean solidified in about 1000 yr; the remaining 20% solidified during the second stage (Meyer et al. 2010; Elkins-Tanton et al. 2011). Magma ocean variations can arise during both stages. During the first one, rotation can result in latitude-dependent crystallization timescale (Maas & Hansen 2019) and produce bottom thickness variations. Once the crust forms, spatial differences in tidal dissipation within the crust can be expected to produce shell thickness variations.

Similarly, ocean thickness variations can reconcile the present inclination of the outer solar system moons with long-lived oceans. To illustrate this, we consider the cases of Callisto and Titan (Figure 5(b)). The case of Titan is especially interesting, as there is evidence of kilometric shell thickness variations (Lefevre et al. 2014). Equator-to-pole thickness variations of just 1 km, ${\hat{c}}_{20}^{2}\sim 1\times {10}^{-3}$, are already sufficient to increase Titan's inclination damping timescale from ∼300 Myr to a value higher than the age of the solar system. Similarly, ocean thickness variations of the same order of magnitude can increase the damping timescale of Callisto, providing an alternative explanation to the recent crossing of the 2:1 mean motion orbital resonance with Ganymede proposed in Downey et al. (2020).

3.2.3. The Tidal Response of an Enceladan Ocean

The tidal response of an Enceladan ice shell of variable thickness has been studied using the thin-shell equations (Beuthe 2018, 2019) and using FEM (Souček et al. 2019). In both cases, the effect of ocean dynamics was ignored. Here we will consider the full tidal response (ice shell and ocean) of Enceladus.

We consider the ice-shell and ocean thicknesses to be (Beuthe et al. 2016)

$$\begin{eqnarray}&&{H}_{\mathrm{ice}}=23(1-0.2326{Y}_{2}^{0}+0.0643{Y}_{3}^{0})\quad [\mathrm{km}]\end{eqnarray} \tag{ 46a }$$

$$\begin{eqnarray}&&{H}_{\mathrm{ocean}}=38(1+0.1408{Y}_{2}^{0}-0.0389{Y}_{3}^{0})\quad [\mathrm{km}].\end{eqnarray} \tag{ 46b }$$

To obtain the tidal response of the ocean and ice shell, we first expand ${\hat{c}}^{2}$, χ, $(\psi \chi {\hat{\alpha }}^{(e)})$, and $(\psi \chi {\hat{D}}^{(b)})$ in spherical harmonics and afterward solve Equation (29d).

Figure 6 shows the resulting surface displacements. As found in Beuthe (2018) and Souček et al. (2019), ice-shell thickness variations lead to higher-amplitude surface displacements at both the equator and the pole. Moreover, the asymmetry introduced by the degree 3 ice-shell thickness variations translates into a hemispherical asymmetry. The south pole experiences higher-amplitude tides.

Zoom In Zoom Out Reset image size

Our results are very similar to the conductive ice-shell case of Beuthe (2018). This is because the surface displacements obtained considering ocean dynamics and ignoring them ( → 0) are almost identical (Figures 6(b) and (c)), which justifies ignoring ocean dynamics to study the tidal deformation of Enceladus's ice shell if the ocean is unstratified.

We now turn our attention to stratified oceans. As opposed to unstratified oceans, stratified oceans can support internal gravity waves. Internal gravity waves account for around 25% of tidal dissipation on Earth (e.g., Munk & Wunsch 1997). Could they also be of relevance in subsurface oceans?

Tyler (2020, 2011) considered the role of internal waves using the "one-and-a-half layer model," which assumes that a stratified ocean can be modeled by using the equations characteristic of an unstratified ocean provided that an effective ocean depth related to the stratification profile is used (Käse & Tyler 2001). Below we will consider in which cases this approximation can be used.

Our method allows us to compute the tidal response of an ocean with an arbitrary number of layers and an arbitrary density profile. Before proceeding, we need to consider which stratification profile is appropriate for subsurface oceans. Each layer adds two new parameters to the problem (see Table 1). As there are few constraints on the stratification of subsurface oceans, the simpler the model the better. Based on this argument, a two-layer model seems appropriate, but can such a simple model capture the behavior of a subsurface ocean?

It is illustrative to first consider Earth's ocean. The stratification profile of Earth's ocean is normally modeled using a three-layer model that consists of a well-mixed layer at the surface; a thin, strongly stratified thermocline; and the deep ocean, which is weakly stratified (Gerkema 2019). The system can be further simplified if the deep ocean is assumed to be well mixed and the thermocline taken to be very thin. Then, Earth's ocean can be approximated as a two-layer system with a density ratio ρ₁/ρ₂ ≈ 0.999.

On Earth, solar heating and freezing and melting at the poles drive the global overturning circulation. Together with mixing due to tides, this results in the stratification structure mentioned above. The dynamics of subsurface oceans is very different from Earth's. Yet several mechanisms can lead to ocean stratification. Depending on pressure and salinity, the maximum density of water is not reached at the freezing point but slightly above. This can result in the formation of a cold "stratosphere" beneath the ice shell. Below the cold stratosphere a well-mixed warmer layer can be present. Melosh et al. (2004) estimate that the thickness of the cold stratosphere for Europa's ocean could be a few hundred meters and ρ₁/ρ₂ ≈ 0.9999. Zeng & Jansen (2021) proposed that the same process can lead to the formation of a stratified layer on top of a well-mixed layer on Enceladus provided that the ocean salinity is low.

On the other end of the spectrum, oceans close to saturation composition can also become stratified. The uptake of salt by buoyant parcels at the seafloor causes upwelling parcels to lose buoyancy before reaching the ice and can produce an ocean divided into a well-mixed warm salty layer below a cold fresh layer of water. Both salinity and temperature gradients decrease upward, and the system can enter the double diffusion convection regime (Vance & Brown 2005; Vance & Goodman 2009).

Finally, a process similar to Earth's overturning circulation might also take place in subsurface oceans. Enceladus ice-shell geometry results in the viscous flow of ice from the equator to high latitudes (Čadek et al. 2019). To maintain the current ice-shell geometry, melting and freezing should occur at the poles and equator, respectively. This can lead to global circulation patterns and a complex lateral and radial stratification profile (Lobo et al. 2021). Lobo et al. (2021) show that these circulation patterns can produce a low-density layer below the ice shell whose thickness changes from 1 to 0 km from pole to equator. The density ratio between the upper and lower layer is ≈0.998. Ashkenazy & Tziperman (2021) showed that the interactions of Europa's ocean with its ice shell can also result in a complex radial and lateral stratification profile characterized by weakly stratified polar regions (r ∼ 0.998) and unstable convective regions in the tropics.

Based on the previous discussion, we will adopt the simple two-layer model. We will first consider the general characteristics of the tidal response of a stratified ocean using the simple two-layer system (Section 3.3.1), and then we will examine the potential contribution of internal tides to Enceladus's heat budget and how it compares with its observed endogenic activity (Section 3.3.2).

3.3.1. Resonances in the Two-layer System

As we did for the one-layer system, we start by considering the nonrotating case. This allows us to find analytical solutions to the equations of motion and study the general characteristics of the ocean response (Appendix E). The eigenfrequencies of the nonrotating system are given by

$$\begin{eqnarray}\begin{array}{rcl}{\hat{\omega }}_{l}^{2} & = & \displaystyle \frac{l(l+1)}{\epsilon }\left(1+{{\rm{\Lambda }}}_{l}-{\xi }_{l}\right)\\ & & \times \,\left[\displaystyle \frac{1}{2}\pm \sqrt{\displaystyle \frac{1}{4}+\displaystyle \frac{r-1}{1+{{\rm{\Lambda }}}_{l}-{\xi }_{l}}{\hat{c}}_{1}^{2}{\hat{c}}_{2}^{2}}\right],\end{array}\end{eqnarray} \tag{ 47 }$$

where r is the density ratio, r ≡ ρ₁/ρ₂, and to simplify notation we have written ${\hat{c}}_{i,0}^{2}$ as ${\hat{c}}_{i}^{2}$. Note that ${\hat{c}}_{i}^{2}$ is simply the ratio between the thickness of layer i and the total ocean thickness (Table 1).

If the ocean is unstratified, r = 1, we recover the eigenfrequencies of the single-layer system—Equation (40). If more layers were added to the system, more internal modes would be added. We can get further insight into the eigenfrequencies of the system by using that, for the moons under consideration, the Lamb parameter is small, → 0. Moreover, we use that the density ratio between the two layers is close to 1, r → 1, but (1 − r)/ is finite. Under these assumptions, Equation (47) reduces to

$$\begin{eqnarray}&&{\hat{\omega }}_{1,l}^{2}=\displaystyle \frac{l(l+1)}{\epsilon }\left(1+{{\rm{\Lambda }}}_{l}-{\xi }_{l}\right),\end{eqnarray} \tag{ 48a }$$

$$\begin{eqnarray}&&{\hat{\omega }}_{2,l}^{2}=\displaystyle \frac{l(l+1)(1-r)}{\epsilon }{\hat{c}}_{1}^{2}{\hat{c}}_{2}^{2}.\end{eqnarray} \tag{ 48b }$$

The first set of eigenmodes (${\hat{\omega }}_{1,l}$) corresponds to the surface gravity waves discussed for unstratified oceans (Section 3.1); the second set (ω_2,l ) corresponds to internal gravity waves. As opposed to surface gravity modes, the eigenfrequencies of internal modes are independent of the effective rigidity of the moon and hence the properties of the solid shell. Internal modes are characterized by the difference between the displacements of the density interface and that of the surface, η_d = η₂ − η₁. This difference is given by (Appendix E)

$$\begin{eqnarray}\begin{array}{rcl}{\hat{\eta }}_{d,l,m} & \approx & {\left[1-{\left(\displaystyle \frac{\hat{\omega }}{{\hat{\omega }}_{2,l}}\right)}^{2}+\displaystyle \frac{i\hat{\alpha }\hat{\omega }}{{\hat{\omega }}_{2,l}^{2}}\right]}^{-1}\\ & & \times \left[\displaystyle \frac{{{\rm{\Lambda }}}_{l}}{1+{{\rm{\Lambda }}}_{l}-{\xi }_{l}}+\displaystyle \frac{{\hat{\omega }}^{2}-i\hat{\alpha }\hat{\omega }}{\left.{\hat{\omega }}_{2,l}^{2}(1+{{\rm{\Lambda }}}_{l}-{\xi }_{l}\right)}{\hat{c}}_{1}^{2}\right]{\hat{\eta }}_{l,m}^{{\rm{T}}}.\end{array}\end{eqnarray} \tag{ 49 }$$

If the ocean is unstratified, ${\hat{\omega }}_{2,l}=0$, we recover the equilibrium tidal response (Equation (42)). The difference between the surface displacement and that of the interface is proportional to ${\hat{c}}_{1}^{2}$. When interface and surface coincide, ${\hat{c}}_{1}^{2}=0$, the difference is zero; in contrast, if the interface is set at the ocean bottom, ${\hat{c}}_{1}^{2}=1$, we recover a zero displacement for it. When the ocean is stratified and close to a resonance, ${\hat{\omega }}_{2,l}^{2}\approx {\hat{\omega }}^{2}$ and ${\hat{\eta }}_{d,l,m}$ becomes very large. The strength of the resonance depends on the moon's rigidity. We can distinguish between two scenarios. If the moon has a low effective rigidity, Λ_l ≪ 1, the term proportional to ${\hat{c}}_{1}^{2}$ is higher than the one proportional to Λ_l . The strength of the resonance is then proportional to ${\hat{c}}_{1}^{2};$ the thinner the upper low-density layer is, the weaker the resonance. In contrast, in the rigid-lid limit, Λ_l ≫ 1, the strength of the resonance is independent of ${\hat{c}}_{1}^{2}$,

$$\begin{eqnarray}&&{\hat{\eta }}_{d,l,m}\approx {\left[1-{\left(\displaystyle \frac{\hat{\omega }}{{\hat{\omega }}_{2,l}}\right)}^{2}+\displaystyle \frac{i\hat{\alpha }\hat{\omega }}{{\hat{\omega }}_{2,l}^{2}}\right]}^{-1}{\hat{\eta }}_{l,m}^{{\rm{T}}}.\end{eqnarray} \tag{ 50 }$$

In this case, the upper layer remains still and the velocities of the upper and lower layers are opposite in direction,

$$\begin{eqnarray}&&{\hat{{\rm{\Phi }}}}_{2,l,m}=-{\hat{{\rm{\Phi }}}}_{1,l,m}=\displaystyle \frac{i\hat{\omega }}{l(l+1)}{\hat{\eta }}_{d,l,m},\end{eqnarray} \tag{ 51 }$$

producing a strong shear. Moreover, the ratio between the magnitude of the velocity of the upper and lower layer is given by the ratio between their thicknesses, ${\hat{c}}_{2}^{2}/{\hat{c}}_{1}^{2}$. Plugging Equation (51) into Equation (38), we can obtain tidal dissipation in the two-layer system in the rigid-lid limit

$$\begin{eqnarray}\begin{array}{rcl}\hat{\dot{E}} & = & -48\pi \left(\displaystyle \frac{{\rho }_{2}-{\rho }_{1}}{{\rho }_{1}}\right)\left[{\hat{\phi }}_{2,0}^{W}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{2,2,0}^{W})\right.\\ & & +\,{\hat{\phi }}_{2,0}^{E}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{2,2,0}^{E})+{\hat{\phi }}_{2,2}^{W}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{2,2,2}^{W})+{\hat{\phi }}_{2,2}^{E}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{2,2,2}^{E})\\ & & +\left.\,{\hat{\phi }}_{2,1}^{W}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{2,2,1}^{W})+{\hat{\phi }}_{2,1}^{E}{\mathfrak{R}}({\hat{{\rm{\Phi }}}}_{2,2,1}^{E})\right],\end{array}\end{eqnarray} \tag{ 52 }$$

showing that the amount of tidal dissipation in the rigid-lid limit is proportional to the density contrast.

We can recover the one-and-a-half layer system used by Tyler (2011, 2020) if we consider that either layer is much thicker than the other, ${\hat{c}}_{2}^{2}\approx 1$ or ${\hat{c}}_{1}^{2}\approx 1$. Equation (48) reduces to

$$\begin{eqnarray}&&{\hat{\omega }}_{2,l}^{2}=\displaystyle \frac{l(l+1)}{{\epsilon }_{r}},\end{eqnarray} \tag{ 53 }$$

with _r = Ω² R²/(1 − r)gH₁ and _r = Ω² R²/(1 − r)gH₂ for ${\hat{c}}_{2}^{2}\approx 1$ and ${\hat{c}}_{1}^{2}\approx 1$, respectively. Comparing Equations (53) and (40), it becomes apparent that the eigenfrequencies of the one-and-a-half layer system are the same as those of the one-layer ocean without a shell and self-gravity with Lamb parameter _r . Furthermore, if there is a rigid lid, Λ_l ≫ 1, the expression for ${\hat{\eta }}_{d,l,m}$ (Equation (50)) is the same as what one would obtain in the one-layer system except that ${\hat{\omega }}_{2}$ appears in the denominator instead of ${\hat{\omega }}_{1}$. Therefore, the two-layer system behaves as the surface-free one-layer system but with = _r if the lower or upper layer is much thicker than the upper layer (${\hat{c}}_{2}^{2}\approx 1$ or ${\hat{c}}_{1}^{2}\approx 1$) and the moon has a high effective rigidity (e.g., Enceladus and Dione). In this case the one-and-a-half layer model of Tyler (2011, 2020) is appropriate to describe the tidal response of a stratified ocean. Solutions of the one-layer system can be used to estimate dissipation in the two-layer system provided that the reduced gravity $g^{\prime} =g({\rho }_{2}-{\rho }_{1})/{\rho }_{2}$ is used instead of the surface gravity to obtain the Lamb parameter and dissipation is scaled with the factor (ρ₂ − ρ₁)/ρ₂.

Resonances with the tidal force occur when $| \hat{\omega }| =1$. As we saw for an unstratified ocean of constant thickness, the eigenfrequencies of surface gravity waves are much higher than the diurnal forcing frequency, ${\hat{\omega }}_{1,l}^{2}\gg 1$, and hence surface gravity wave resonances do not occur. However, this is not necessarily the case for internal gravity waves. Using Equation (48b), we find the condition for internal gravity wave resonances to be

$$\begin{eqnarray}&&{\hat{c}}_{1}^{2}=\displaystyle \frac{1}{2}\pm \sqrt{\displaystyle \frac{1}{4}-\displaystyle \frac{\epsilon }{l(l+1)(1-r)}},\end{eqnarray} \tag{ 54 }$$

from which it follows that there exist two different configurations for which a mode of degree l can resonate. From Equation (54) it also follows that resonances of degree l can be excited as far as

$$\begin{eqnarray}&&1-r\gt 4\epsilon /l(l+1).\end{eqnarray} \tag{ 55 }$$

The small value of Lamb parameter characteristic of subsurface oceans (Table 2) makes it possible for internal wave resonance to be excited for small density contrasts, r close to 1. In the limit of 1 − r ≫ 4/l(l + 1), Equation (54) reduces to ${\hat{c}}_{1}^{2}=\epsilon /l(l+1)(1-r)$ and ${\hat{c}}_{1}^{2}=1\,-\epsilon /l(l+1)(1-r)$. These two resonant thicknesses correspond to the limiting cases of a ${\hat{c}}_{2}^{2}\approx 1$ and ${\hat{c}}_{1}^{2}\approx 1$ discussed above in the context of the one-and-a-half layer system. In this limit there exist two resonances, occurring when either the lower layer or the upper layer has _r = l(l + 1).

The previous discussion holds for nonrotating moons. As was the case for unstratified oceans, in nonrotating moons modes of different degrees are not coupled. Therefore, the tidal potential can only excite the degree 2 modes. In contrast, when the Coriolis force is included, a forcing of degree 2 can excite modes of higher degree.

The main characteristics of the tidal response of a two-layer ocean are shown in Figure 7, which shows tidal dissipation in a stratified ocean due to the eccentricity tide for different stratification profiles (r, ${\hat{c}}_{1}^{2}$), Lamb parameter (), and rigidity (γ). In all cases, bands of intense tidal dissipation are observed. Each of the bands corresponds to a different resonant mode. As expected from the simple nonrotating case (Equation (54)), given a density contrast, there are two values of ${\hat{c}}_{1}^{2}$ for which the mode can resonate. We can also observe that when the lower layer is much thicker than the upper layer (${\hat{c}}_{1}^{2}\ll 1$), the resonant bands have a slope of ≈−1 in the ${\rm{log}}(1-r)$–${\rm{log}}({\hat{c}}_{1}^{2})$ space, as expected from the one-and-a-half layer approximation in the nonrotating case (Equation (53)).

Zoom In Zoom Out Reset image size

The shape and intensity of these bands depend on the Lamb parameter and the effective rigidity γ. Figures 7(a) and (c) illustrate the effect that the Lamb parameter has on the excitation of internal modes. The Lamb parameter alters the shape of the resonant bands; as decreases, the resonant bands shift to smaller density contrasts. For a given density contrast 1 − r, more resonances can be excited the smaller is. This is expected from the nonrotating model, Equation (55), from which it follows that the lower is, the smaller the density contrast needs to be for a resonance of degree l to be excited.

Figures 7(a) and (b) show how the effective rigidity of the moon also affects the excitation of internal waves. First, we notice that the background tidal dissipation is higher the lower the effective rigidity γ, which follows from the fact that a lower rigidity results in more intense surface gravity waves (see also Figure 2). Second, we observe that the locations of the resonant bands are the same but their intensities decrease, which follows from Equations (48b) and (49). For a hard-shell moon, the strength of the resonances is independent of both effective rigidity and the ${\hat{c}}_{1}^{2}$ for which it occurs. In contrast, for soft-shell moons, the intensity of the resonance depends on ${\hat{c}}_{1}^{2}$. As ${\hat{c}}_{1}^{2}$ decreases, the amplitude of the resonance decreases.

We also observe that tidal dissipation depends on the density contrast 1 − r (Equation (52)). This is more evident when the moon has a high effective rigidity and low Lamb parameter (Figure 7(c)). We observe that outside resonances dissipation decreases as 1 − r decreases until it becomes independent of 1 − r when internal waves have a smaller contribution to the energy budget than surface waves. Similarly, the intensity of a given resonance band also decreases with 1 − r.

3.3.2. Tidal Dissipation in a Stratified Enceladan Ocean

In the previous section we considered the tidal response of stratified oceans and showed that internal wave resonances can be excited. We saw that internal wave resonances are stronger for hard-shell moons such as Enceladus and Dione. With this in mind, we will consider the specific case of Enceladus because of its relatively thick ice shell to assess whether internal wave resonances might be behind the moon's thermal activity.

Figure 8 shows tidal dissipation for different stratification profiles and Rayleigh dissipation parameters α and r. For some combinations of ocean stratification and Rayleigh parameter, resonances can result in more tidal heating than Enceladus's thermal output (Howett et al. 2011). For a strongly stratified ocean (1 − r = 10⁻²), this occurs for α = 10⁻⁷ s⁻¹. As the ocean becomes less stratified, dissipation decreases (Equation (52)), and a lower α is required to reconcile tidal heating with Enceladus's thermal output. For a density contrast similar to that used in Lobo et al. (2021) and similar to Earth's, 1 − r = 10⁻³, an α as small as 10⁻¹⁰ s⁻¹ is required. As α decreases, the dissipation peak becomes sharper, and hence the region of the parameter space for which high values of tidal dissipation are attained narrows. While energy dissipation from internal waves can be close to Enceladus's thermal output, the geographical patterns are characterized by dissipation at low latitudes (Figure 9), which is at odds with the structure of Enceladus's ice shell and the location of the highest endogenic activity at the south pole. Though our model produces a tidal heating distribution that is inconsistent with these observations, it does not preclude the possibility that internal waves are a significant contributor to Enceladus's total thermal budget. The disagreement between predicted and observed heating patterns might indicate that the constant-thickness two-layer model employed is too simplistic. Furthermore, it is uncertain how ocean circulation might redistribute the heat generated as a result of tides.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In resonances, strong ocean currents develop. For example, for 1 − r = 10⁻³, ${\hat{c}}_{1}^{2}=6.573\times {10}^{-1}$, and α = 10⁻⁹ s⁻¹, velocities of ∼1 m s⁻¹ are attained. When this occurs, nonlinear effects can become important (Hay et al. 2020). This is especially aggravated for internal wave resonances, as a shear develops at the interface between the two layers (Equation (51)), which can lead to Kelvin–Helmholtz instabilities, turbulence, and mixing (Stewart & Dellar 2013). On Earth's oceans the Kelvin–Helmholtz instability plays an important role in ocean mixing and stratification (Smyth & Moum 2012). A similar process might occur in subsurface oceans.