Enceladus's crust as a non-uniform thin shell: II tidal dissipation

TL;DR

This study investigates how non-uniform icy shell thickness on Enceladus influences tidal heating distribution, finding that shell variations alone cannot explain observed south pole heating, suggesting additional core dissipation or non-steady processes.

Contribution

It introduces a non-uniform thin shell model to analyze tidal dissipation, highlighting the limited role of shell heterogeneity in explaining Enceladus's heating asymmetry.

Findings

Shell dissipation flux increases where the shell is thinner.

Nominal shell dissipation is insufficient to match observed heat flux.

Core dissipation alone cannot produce the observed asymmetry.

Abstract

Tidal heating is the prime suspect behind Enceladus's south polar heating anomaly and global subsurface ocean. No model of internal tidal dissipation, however, can explain at the same time the total heat budget and the focusing of the energy at the south pole. I study here whether the non-uniform icy shell thickness can cause the north-south heating asymmetry by redistributing tidal heating either in the shell or in the core. Starting from the non-uniform tidal thin shell equations, I compute the volumetric rate, surface flux, and total power generated by tidal dissipation in shell and core. If the shell is laterally uniform, the thin shell approach predicts shell dissipation with a few percent error while the error on core dissipation is negligible. Variations in shell thickness strongly increase the shell dissipation flux where the shell is thinner. For a hard shell with long-wavelength variations, the shell dissipation flux can be predicted by scaling with the inverse local thickness the flux for a laterally uniform shell. If Enceladus's shell is in conductive thermal equilibrium with isostatic thickness variations, the nominal shell dissipation flux at the south pole is about three times its value for a shell of uniform thickness, which remains negligible compared to the observed flux. The shell dissipation rate should be ten times higher than nominal in order to account for the spatial variations of the observed flux. Dissipation in an unconsolidated core can provide the missing power, but does not generate any significant heating asymmetry as long as the core is homogeneous. Non-steady state models, though not investigated here, face similar difficulties in explaining the asymmetries of tidal heating and shell thickness.

Tables

Table 1. Table 1: Physical parameters used in this paper.

Parameter	Symbol	Value	Unit
Mean eccentricitya	$e$	$0.0047$	-
Rotation ratea	$\omega$	$5.307\times {10}^{-5}$	$s^{-1}$
Bulk densitya	${\rho }_b$	$1610$	$\mathrm{kg}{m}^{-3}$
Surface gravitya	$g$	$0.1135$	$m{s}^{-2}$
Surface radiusa	$R$	$252.1$	km
Reference thicknessb (if uniform shell)	$d$	23	km
Core radiusb	$R_c$	192	km
Density of ice and ocean	$\rho$	1000	$\mathrm{kg}/m^3$
Shear modulus of icec	${\mu }_e$	3.5	GPa
Bulk modulus of icec	$K_e$	9.13	GPa
Poisson’s ratio of icec	${\nu }_e$	0.33	-
Shear modulus of core (if non-porous)	${\mu }_{\mathrm{ce}}$	40	GPa
Conductivity of iced (if uniform shell)	$k_{ice}$	$651/T$	${\mathrm{Wm}}^{-1}{K}^{-1}$
Activation energy	$E_a$	$59.4$	$\mathrm{kJ}{\mathrm{mol}}^{-1}$
Melting ice temperature	$T_m$	273	K
Viscosity of ice at melting point	${\eta }_m$	${10}^{13}-{10}^{15}$	Pa.s
a Source given in Table 1 of Beuthe [2016].
b Beuthe et al. [2016].
c Helgerud et al. [2009].
d Petrenko and Whitworth [1999] (Section 4 only); Eq. (63) is used in Section 5.

Table 2. Table 2: Viscoelastic shell parameters.

Symbol	Name	Definition
Depth-dependent
$\zeta$	radial shell coordinate	$r-R$
$z$	depth parameter	$\zeta /\left(R+\zeta \right)-\left(1-\chi \right)$
$\eta$	viscosity of icea	${\eta }_m\exp \left(E_a\left(T^{-1}-T_m^{-1}\right)/R_g\right)$
$\mu$	complex shear modulusb	${\mu }_e/\left(1-i{\mu }_e/\omega \eta \right)$
Depth-integrated ($\epsilon =d/R$)
${\mu }_p$	$p$th moment of $\mu$	${\left(1/d\right)}^{p+1}{\int }_d\mu {\zeta }^p𝑑\zeta$
${\mu }_{\mathrm{inv}}$	invariant second moment	${\mu }_2-{\left({\mu }_1\right)}^2/{\mu }_0$
$\chi$	-	$\left({\mu }_0+\epsilon {\mu }_1\right)/\left({\mu }_0+2\epsilon {\mu }_1+{\epsilon }^2{\mu }_2\right)$
$\psi$	-	${\mu }_0/\left({\mu }_0+\epsilon {\mu }_1\right)$
${\alpha }_{\mathrm{inv}}$	invariant extensibility	${\left(2\left(1+\nu \right){\mu }_{0}d\right)}^{-1}$
$D_{\mathrm{inv}}$	invariant bending rigidity	$2{\mu }_{\mathrm{inv}}{d}^3/\left(1-\nu \right)$
$\alpha$	extensibility	$\chi \psi {\alpha }_{\mathrm{inv}}$
$D$	bending rigidity	$\chi \psi D_{\mathrm{inv}}$
a $R_g=8.314J{K}^{-1}{\mathrm{mol}}^{-1}$ is the gas constant.
b For Maxwell rheology, but other linear rheological models can be used.

Table 3. Table 3: Degree- n 𝑛 n tidal deformation of a laterally uniform thin shell (see Section 4 of Paper I, except for U n ∘ superscript subscript 𝑈 𝑛 U_{n}^{\circ} which is given by Eq. ( 36 )).

Symbol	Name	Solutiona,b	Examplec
Basic variables
$w_n$	radial displacement	$h_n{U}_n^T/g$
$F_n$	stress function	$\left(\chi /\alpha \right)/\left({\delta }_n^{\prime }-1-\nu \right)\left(w_n/R\right)$
${\Gamma }_n$	total perturbing potential	$\left(k_n+1\right)U_n^T$
$q_n$	tidal load	$\rho g{\Lambda }_n{w}_n$
$U_n^{\circ }$	effective tidal potential	$\left(\left(k_n+1\right)/\left(k_n^{\circ }+1\right)\right)U_n^T$
Spring Constants (SC)
${\Lambda }_n^M$	membrane SC	${\delta }_n^{\prime }/\left({\delta }_n^{\prime }-1-\nu \right){\left(\rho g{R}^2{\alpha }_{\mathrm{inv}}\right)}^{-1}$	$21.267+0.730i$
${\Lambda }_n^B$	bending SC	${\delta }_n^{\prime }\left({\delta }_n^{\prime }-1+\nu \right)D_{\mathrm{inv}}/\left(\rho g{R}^4\right)$	$0.378+0.038i$
${\Lambda }_n^{corr}$	next-to-leading SC	$\left(\chi /\psi -1\right){\Lambda }_n^M+\left(\chi \psi -1\right){\Lambda }_n^B$	$0.020+0.003i$
${\Lambda }_n$	thin shell SC	${\Lambda }_n^M+{\Lambda }_n^B+{\Lambda }_n^{corr}$	$21.665+0.771i$
Tidal Love numbers (TLN)
$h_n^{\circ }$	fluid-crust radial TLN	Eqs. (D.2)-(D.3)	$1.594$
$k_n^{\circ }$	fluid-crust gravitational TLN	$h_n^{\circ }-1$	$0.594$
$h_n$	radial TLN	$h_n^{\circ }/\left(1+\left(1+{\xi }_n{h}_n^{\circ }\right){\Lambda }_n\right)$	$0.0448-0.0015i$
$k_n$	gravitational TLN	$\left(1+{\Lambda }_n\right)h_n-1$	$0.0167-0.0006i$
a $(\alpha ,{\alpha }_{\mathrm{inv}},D_{\mathrm{inv}},\chi ,\psi )$ are defined in Table 2; numerical values are given in Table 3 of Paper I.
b ${\xi }_n$ is the degree-$n$ density ratio, see Eq. (5); ${\delta }_n^{\prime }=-\left(n-1\right)\left(n+2\right)$, see Appendix A.
c If $n=2$, $d=23\mathrm{km}$, ${\eta }_m={10}^{13}\mathrm{Pa}.s$; and the core is elastic (other parameters given in Table 1).

Table 4. Table 4: Radial weights and angular functions factorizing the dissipation rate in a laterally uniform body forced by degree-2 tides (Eq. ( 45 )). The harmonic basis is given for eccentricity tides combined with the forced libration: γ 1 = 1 + γ 0 / 2 ​ e subscript 𝛾 1 1 subscript 𝛾 0 2 𝑒 \gamma_{1}=1+\gamma_{0}/2e (Eq. ( 7 )). The symbols y i subscript 𝑦 𝑖 y_{i} are the viscoelastic-gravitational solutions of Takeuchi and Saito [ 1972 ] . If the body is incompressible, r ​ ∂ r y 1 = − 2 ​ y 1 + 6 ​ y 3 𝑟 subscript 𝑟 subscript 𝑦 1 2 subscript 𝑦 1 6 subscript 𝑦 3 r\partial_{r}y_{1}=-2y_{1}+6y_{3} and f K = 0 subscript 𝑓 𝐾 0 f_{K}=0 . f T subscript 𝑓 𝑇 f_{T} is the weight for the angular average of the dissipation rate. P n ​ m subscript 𝑃 𝑛 𝑚 P_{nm} denote the unnormalized associated Legendre functions of degree n 𝑛 n and order m 𝑚 m with argument cos ⁡ θ 𝜃 \cos\theta .

$J$	Radial weight $f_J$	Angular function ${\Psi }_J$
$A$	$\left(4/3\right){\left|r{\partial }_r{y}_1-y_1+3y_3\right|}^2$	${\Psi }_0+{\Psi }_2+{\Psi }_4$
$B$	$6{\left|r{y}_4/\mu \right|}^2$	${\Psi }_0+\left(1/2\right){\Psi }_2-\left(2/3\right){\Psi }_4$
$C$	$24{\left|y_3\right|}^2$	${\Psi }_0-{\Psi }_2+\left(1/6\right){\Psi }_4$
$K$	$\frac{\mathrm{Im}\left(K\right)}{\mathrm{Im}\left(\mu \right)}{\left|r{\partial }_r{y}_1+2y_1-6y_3\right|}^2$	${\Psi }_0+{\Psi }_2+{\Psi }_4$
$T$	$f_A+f_B+f_C+f_K$	${\Psi }_0$
$n$	Harmonic basis ${\Psi }_n$
0	$\frac{3}{5}\left(3+4{\gamma }_1^2\right)e^2$
2	$-\frac{3}{7}\left(3+8{\gamma }_1^2\right)e^2{P}_{20}+\frac{9}{14}{e}^2{P}_{22}\cos 2\phi$
4	$\frac{9}{140}\left(27+16{\gamma }_1^2\right)e^2{P}_{40}-\frac{27}{140}{e}^2{P}_{42}\cos 2\phi +\frac{3}{1120}\left(9-16{\gamma }_1^2\right)e^2{P}_{44}\cos 4\phi$

Table 5. Table 5: Surface shell flux of a laterally uniform shell: percentage contribution of Pattern C to the average surface flux. Pattern A contributes the remainder. The shell is 23 ​ km 23 km 23\rm\,km thick and the bottom viscosity is either 10 13 superscript 10 13 10^{13} or 10 14 ​ Pa . s formulae-sequence superscript 10 14 Pa s 10^{14}\rm\,Pa.s (results separated by | | | ). See Section  4.2 .

Dissipation constraint	Thick shell	Thin shell	Membranea
No Poisson dissipation	92.2 $|$ 93.1	91.7 $|$ 91.8	75.0 $|$ 75.0
No bulk dissipation	92.0 $|$ 94.0	-	76.0 $|$ 80.9
Incompressible	96.3 $|$ 96.5	96.5 $|$ 96.6	81.8 $|$ 81.8
a Eq. (E.7) with $\nu =0.33$, $\overline{\kappa }=0.42|\mathrm{\hspace{0.17em}0.06}$, and $\nu =0.5$, respectively.

Table 6. Table 6: Thermal equilibrium models: characteristics, total power (conductive/shell/core), partition of the shell power and feedback effect. UNI/ISO/THIN specify the shell thickness model; L/H denote the low/high level of shell dissipation; C indicates that core dissipation occurs (shell and core structures are detailed in Section  5.2 ).

Models	UNI-L	ISO-L	UNI-H	ISO-H	ISO-LC	ISO-HC	THIN-LC
Shell thickness	uniform	isostatic	uniform	isostatic	isostatic	isostatic	thin SPT
Shell dissipation	low	low	high	high	low	high	low
Core dissipation	–	–	–	–	high	high	high
Power (GW)
Conductive power	31.0	34.5	31.8	35.4	34.5	35.5	40.3
Shell power	0.96	1.16	12.9	15.3	1.22	15.9	2.99
Core power	0	0	0	0	33.3	19.5	37.3
Shell partition (%)
Membrane	86.7	89.8	85.1	89.2	89.8	89.2	57.4
Mixing	7.8	2.0	9.7	3.2	2.0	3.3	-3.0
Bending	5.6	8.2	5.2	7.6	8.2	7.5	45.6
Rheology feedback (%)	2	2	37	35	2	37	5

Table 7. Table 7: Spherical differential operators: definitions. See Appendices A and B of Paper I and Appendix B of Beuthe [ 2013 ] . The null space is characterized by the harmonic degree of the functions belonging to it.

Notation	Definition	Null space
Tensorial operators
${\overline{𝒪}}_1={𝒪}_1-1$	${\partial }_{\theta }^2$	$n=0$
${\overline{𝒪}}_2={𝒪}_2-1$	${\left(\sin \theta \right)}^{-2}{\partial }_{\phi }^2+\cot \theta {\partial }_{\theta }$	$n=0$
${\overline{𝒪}}_3={𝒪}_3$	${\left(\sin \theta \right)}^{-1}\left({\partial }_{\theta }{\partial }_{\phi }-\cot \theta {\partial }_{\phi }\right)$	$n=0,1$
Scalar operators
$\Delta$	${\overline{𝒪}}_1+{\overline{𝒪}}_2$	$n=0$
${𝒟}_2(a;b)$	$\left({\partial }_{\theta }a\right)\left({\partial }_{\theta }b\right)+{\left(\sin \theta \right)}^{-2}\left({\partial }_{\phi }a\right)\left({\partial }_{\phi }b\right)$	$n=0$
${𝒟}_4(a;b)$	$\left({\overline{𝒪}}_{1}a\right)\left({\overline{𝒪}}_{1}b\right)+\left({\overline{𝒪}}_{2}a\right)\left({\overline{𝒪}}_{2}b\right)+2\left({\overline{𝒪}}_{3}a\right)\left({\overline{𝒪}}_{3}b\right)$	$n=0$
${\Delta }^{\prime }$	${𝒪}_1+{𝒪}_2$	$n=1$
${𝒟}_4^{\prime }(a;b)$	$\left({𝒪}_{1}a\right)\left({𝒪}_{1}b\right)+\left({𝒪}_{2}a\right)\left({𝒪}_{2}b\right)+2\left({𝒪}_{3}a\right)\left({𝒪}_{3}b\right)$	$n=1$
$𝒞(a;b)$	${\Delta }^{\prime }\left(a{\Delta }^{\prime }b\right)$	$n=1$
$𝒜(a;b)$	$\left({𝒪}_{1}a\right)\left({𝒪}_{2}b\right)+\left({𝒪}_{2}a\right)\left({𝒪}_{1}b\right)-2\left({𝒪}_{3}a\right)\left({𝒪}_{3}b\right)$	$n=1$

Table 8. Table 8: Spherical differential operators: identities. See Appendices A and B of Paper I and Appendix B of Beuthe [ 2013 ] .

Identity				Remark
$a$ and $b$ are arbitrary functions of $(\theta ,\phi )$, $a_{\mathrm{uni}}$ is uniform:
(a)	${\Delta }^{\prime }$	=	$\Delta +2$	–
(b)	${𝒟}_2(a;b)$	=	$\frac{1}{2}\left[\Delta \left(ab\right)-\left(\Delta a\right)b-a\left(\Delta b\right)\right]$	–
(c)	${𝒟}_4(a;b)$	=	$\left(\Delta a\right)\left(\Delta b\right)+\left(\Delta a\right)b+a\left(\Delta b\right)+2ab-𝒜(a;b)$	–
(d)	${𝒟}_4^{\prime }(a;b)$	=	$\left({\Delta }^{\prime }a\right)\left({\Delta }^{\prime }b\right)-𝒜(a;b)$	–
(e)	$𝒜(a;b)$	=	$\frac{1}{4}[-{\Delta }^{\prime }{\Delta }^{\prime }\left(ab\right)-\left({\Delta }^{\prime }{\Delta }^{\prime }a\right)b-a\left({\Delta }^{\prime }{\Delta }^{\prime }b\right)$	–
			$+\mathrm{\hspace{0.17em}2}\left({\Delta }^{\prime }a\right)\left({\Delta }^{\prime }b\right)+\mathrm{\hspace{0.17em}2}{\Delta }^{\prime }\left(\left({\Delta }^{\prime }a\right)b+a\left({\Delta }^{\prime }b\right)\right)$
			$-\mathrm{\hspace{0.17em}2}({\Delta }^{\prime }\left(ab\right)+\left({\Delta }^{\prime }a\right)b+a\left({\Delta }^{\prime }b\right))+8ab]$
(f)	$𝒜(a_{\mathrm{uni}};b)$	=	$a_{\mathrm{uni}}{\Delta }^{\prime }b$	–
$a_n$ is a spherical harmonic of degree $n$:
(g)	$\Delta a_n$	=	${\delta }_{n}a$	${\delta }_n=-n\left(n+1\right)$
(h)	${\Delta }^{\prime }{a}_n$	=	${\delta }_n^{\prime }a$	${\delta }_n^{\prime }=-\left(n-1\right)\left(n+2\right)$
(i)	$𝒜(a_n;a_n^{\ast })$	=	$\frac{1}{4}\left[-\Delta \Delta +\left(4{\delta }_n^{\prime }-6\right)\Delta +4{\delta }_n^{\prime }\right]{\left|a_n\right|}^2$
(j)	$⟨𝒜(a_n;a_n^{\ast })⟩$	=	${\delta }_n^{\prime }⟨{\left|a_n\right|}^2⟩$	$⟨x⟩={\left(4\pi \right)}^{-1}{\int }_{S}x𝑑\Omega$

Table 9. Table 9: Non-dimensional functions for core deformations. Variables are the harmonic degree n 𝑛 n , the density ratio ξ = ρ / ρ b 𝜉 𝜌 subscript 𝜌 𝑏 \xi=\rho/\rho_{b} , the relative core radius y = R c / R 𝑦 subscript 𝑅 𝑐 𝑅 y=R_{c}/R , and the reduced radius r ^ = r / R c ^ 𝑟 𝑟 subscript 𝑅 𝑐 \hat{r}=r/R_{c} .

$f_n$	=	$n{\left[2\left(n-1\right)\left(3+4n+2n^2\right)\right]}^{-1}$
$p_A$	=	$\left(2\left(n-1\right)+3y^{2n+1}\right)\left(1-\xi \right)+\left(2n+1\right)y^3\xi$
$p_B$	=	$\left(2n+1-3\xi \right)\left[2\left(n-1\right)\left(1-\xi \right)+\left(2n+1\right)y^3\xi \right]-9\left(1-\xi \right)y^{2n+1}\xi$
$A_n$	=	$f_n\left(2n+1\right)\left(1-\xi \right)p_A$
$B_n$	=	$f_n\left(1-\xi \right)p_B$
$f_{n1}\left(\hat{r}\right)$	=	$\frac{1}{2n+1}\left(n\left(n+2\right)-\left(n^2-1\right){\hat{r}}^2\right){\hat{r}}^{n-1}$
$f_{n3}\left(\hat{r}\right)$	=	$\frac{1}{n\left(2n+1\right)}\left(n\left(n+2\right)-\left(n-1\right)\left(n+3\right){\hat{r}}^2\right){\hat{r}}^{n-1}$
$f_{n4}\left(\hat{r}\right)$	=	$\frac{2}{2n+1}\left(n-1\right)\left(n+2\right)\left(1-{\hat{r}}^2\right){\hat{r}}^{n-1}$

Equations

$$\mathcal{C}(D;w)-(1-\nu )\mathcal{A}(D;w)+R^3\mathcal{A}(\chi ;F)$$

$$\mathcal{C}(\alpha ;F)-(1+\nu )\mathcal{A}(\alpha ;F)-R^{-1}\mathcal{A}(\chi ;w)$$

$$q=q_0+\sum _{n\ge 2}{q}_n,$$

$$q_n=-\rho g\left(w_n-{\Gamma }_n/g\right),$$

$${\Gamma }_n={\upsilon }_n\left(U_n^T+g{\xi }_n{w}_n\right),$$

$${\xi }_n=\frac{3}{2n+1}\frac{\rho }{{\rho }_b}.$$

$${\upsilon }_n=\frac{h_n^{\circ }}{1+{\xi }_n{h}_n^{\circ }},$$

$$U(t,\theta ,\phi )=(\omega R{)}^{2}e\left(-\frac{3}{2}{P}_{20}\cos \omega t+P_{22}\left(\frac{3}{4}\cos 2\phi \cos \omega t+\left(1+\frac{{\gamma }_0}{2e}\right)\sin 2\phi \sin \omega t\right)\right),$$

$$P(r,\theta ,\phi )=\frac{1}{T}{\int }_0^T{\sigma }_{ij}(t){\dot{ϵ}}_{ij}(t)dt=\frac{\omega }{2}\mathrm{Im}\left({\tilde{\sigma }}_{ij}(\omega ){\widetilde{ϵ}}_{ij}^{\ast }(\omega )\right),$$

$${\sigma }_{ij}=2\mu {ϵ}_{ij}+\left(K-2\mu /3\right)ϵ{\delta }_{ij},$$

$$P(r,\theta ,\phi )=\omega \mathrm{Im}(\mu )\left({ϵ}_{ij}{ϵ}_{ij}^{\ast }-\frac{1}{3}{\left|ϵ\right|}^2\right)+\frac{\omega }{2}\mathrm{Im}(K){\left|ϵ\right|}^2.$$

$$K=\frac{2}{3}\frac{1+\nu }{1-2\nu }\mu .$$

$${\rm Im}(K)=\frac{2}{3}\,\frac{1+\nu}{1-2\nu}\,{\rm Im}(\mu)\hskip 14.22636pt\mbox{($\nu$ is real)}\,.$$

$$P(r,\theta ,\phi )=\omega \mathrm{Im}(\mu )\left({ϵ}_{ij}{ϵ}_{ij}^{\ast }+\frac{\nu }{1-2\nu }{\left|\mathrm{Tr}{ϵ}_{ij}\right|}^2\right).$$

$${ϵ}_{\zeta \zeta }=-\frac{\nu }{1-\nu }\left({ϵ}_{\theta \theta }+{ϵ}_{\phi \phi }\right).$$

$$P_{shell}(\zeta ,\theta ,\phi )=\omega \mathrm{Im}(\mu )\left({\mathcal{E}}_2+\frac{\nu }{1-\nu }{\mathcal{E}}_{\mathrm{tr}}\right),$$

$${\mathcal{E}}_2$$

$${\mathcal{E}}_{\mathrm{tr}}$$

$$P_{shell}(\zeta ,\theta ,\phi )=\omega \mathrm{Im}(\mu )\left({\mathcal{E}}^{FF}+{\mathcal{E}}^{Fw}+{\mathcal{E}}^{ww}\right).$$

$${\mathcal{E}}^{FF}$$

$${\mathcal{E}}^{Fw}$$

$${\mathcal{E}}^{ww}$$

$${\mathcal{F}}_{shell}(\theta ,\phi )={\int }_d{P}_{shell}(\zeta ,\theta ,\phi ){\left(1+\zeta /R\right)}^{2}d\zeta ,$$

$${\mathcal{F}}_{shell}(\theta ,\phi )={\mathcal{F}}_{mem}+{\mathcal{F}}_{mix}+{\mathcal{F}}_{bend}.$$

$${\mathcal{F}}_{mem}$$

$${\mathcal{F}}_{mix}$$

$${\mathcal{F}}_{bend}$$

$${\dot{E}}_{shell}$$

$${\dot{E}}_{tot}=\frac{\omega }{2}\frac{1}{4\pi GR}\sum _n\left(2n+1\right){\int }_S\mathrm{Im}\left(U_n^T{U_n^{\prime }}^{\ast }\right)dS,$$

$$U_n^{\prime }={\Gamma }_n-U_n^T,$$

$${\dot{E}}_{tot}=-\frac{5}{2}\frac{\omega R}{G}\mathrm{Im}({\upsilon }_2)⟨|U_2^T{|}^2⟩+2\pi R^2\omega \rho \mathrm{Im}⟨U_2^T{\upsilon }_2^{\ast }{w}_2^{\ast }⟩,$$

$${\dot{E}}_{tot}={\dot{E}}_{shell}=\frac{\omega }{2}\mathrm{Im}{\int }_S{q}_2{w}_2^{\ast }dS,$$

$$U_n^P={\xi }_n{q}_n/\rho ,$$

Full text

1. 1 Introduction

2. 2 Tidal thin shell equations

3. 2.1 Flexure equations

4. 2.2 Tidal forcing

5. 3 Dissipation in the thin shell approach

6. 3.1 Shear, bulk and Poisson dissipation

7. 3.2 Dissipation inside the shell

8. 3.2.1 Shell dissipation rate

9. 3.2.2 Shell dissipation flux

10. 3.2.3 Shell power

11. 3.3 Dissipation inside the core

12. 3.3.1 Effective tidal potential for the core

13. 3.3.2 Core-shell partition of total power

14. 3.3.3 Core dissipation rate

15. 4 Benchmarking against a laterally uniform thick shell

16. 4.1 Factorization of dissipation rate

17. 4.2 Patterns in a laterally uniform thin shell

18. 4.3 Total power in core and shell

19. 4.4 Dissipation in the core

20. 5 Thermal equilibrium in a conductive shell

21. 5.1 Coupling dissipation to heat transfer

22. 5.2 Shell structure and core rheology

23. 5.3 Conductive model

24. 5.4 Surface temperature

25. 5.5 Results and discussion

26. 5.5.1 Shell dissipation

27. 5.5.2 Partition and scaling

28. 5.5.3 Core dissipation

29. 6 Summary

30. A Spherical differential operators

31. B Dissipation rate in terms of $F$ and $w$

32. C Surface flux in terms of $F$ and $w$

33. D Deformation of the core

34. E Dissipation in a laterally uniform thin shell

1. 1 Impact of dissipation condition
2. 2 Basic patterns of tidal heating in a laterally uniform shell
3. 3 Dissipation rate in a laterally uniform conductive shell
4. 4 Surface flux pattern in a laterally uniform conductive shell
5. 5 Core dissipation under a laterally uniform shell
6. 6 Total power dissipated in a laterally uniform shell
7. 7 Core dissipation under a laterally uniform shell
8. 8 Shell structure
9. 9 Conductive model
10. 10 Surface temperature
11. 11 Heat flux pattern in model ISO-L
12. 12 Heat flux and rheology feedback if no core dissipation
13. 13 Partition and scaling of the shell dissipation flux
14. 14 Partition of the shell dissipation flux into membrane/mixed/bending contributions
15. 15 Surface heat flux if high dissipation within the core
16. 16 Surface heat flux if weak shell in SPT

1. 1 Physical parameters
2. 2 Viscoelastic shell parameters
3. 3 Tidal deformation of a laterally uniform thin shell
4. 4 Radial weights and angular functions factorizing the dissipation rate
5. 5 Pattern of surface shell flux
6. 6 Thermal equilibrium models
7. 7 Spherical differential operators: definitions
8. 8 Spherical differential operators: identities
9. 9 Non-dimensional functions for core deformations

1 Introduction

Despite being small and cooling fast, Enceladus is warm. The first evidence for its present-day warm state was the detection by Cassini of a south polar hot spot [Spencer et al., 2006; Howett et al., 2011], with temperature reaching $200\rm\,K$ along the south polar faults [Goguen et al., 2013; Spencer et al., 2018]. A second line of evidence arises from the presence of underground liquid water, first inferred at the south pole from the composition of Saturn’s E ring [Postberg et al., 2009] and Enceladus’s plume [Waite et al., 2009; Postberg et al., 2011]. Geodesy measurements (gravity, topography, and libration) then proved beyond doubt the existence of a global ocean beneath a thin shell of non-uniform thickness [Iess et al., 2014; Thomas et al., 2016; Beuthe et al., 2016].

Tidal heating is the only source of energy sufficient to keep Enceladus warm [Squyres et al., 1983; Ross and Schubert, 1989], and might be sufficient to keep it in a steady state [Fuller et al., 2016; Nimmo et al., 2018]. Energy conservation in the Saturn-Enceladus-Dione system, however, does not constrain where and how energy is dissipated inside Enceladus. At the present time, viscoelastic dissipation within the shell is too small to produce the observed heat flux and to maintain thermal equilibrium [Squyres et al., 1983; Roberts and Nimmo, 2008], while dissipation within the ocean is completely negligible [Beuthe, 2016; Rovira-Navarro et al., 2019]. Theoretically, dissipation within an unconsolidated core can be large enough [Ross and Schubert, 1989], but it implies unverified assumptions about the viscoelastic behaviour of porous matter [Roberts, 2015; Choblet et al., 2017]. Instead of being in a steady state, Enceladus could oscillate in a thermal-orbital evolution cycle between phases of high and low dissipation. Most models of this sort postulate that Enceladus is now in a short-lived stage of high dissipation, leading to the same problem as above regarding heat production, besides raising the question of Enceladus’s special status at the present time (see review in Nimmo et al. [2018]). Alternatively, Enceladus could oscillate between two phases of nearly equal duration: the present-day phase of low dissipation in a medium-thick shell and another phase of high dissipation in a very thin shell [Luan and Goldreich, 2017].

The south polar localization of the thermal anomaly and the spatial variations of shell thickness are often studied separately from the problem of the global heat budget. This means throwing away a strong constraint on dissipation models. For example, does a very dissipative core produce a north-south heating asymmetry if the shell is non-uniform, as suggested by Choblet et al. [2017]? Or else, is the non-steady state model compatible with the generation of shell thickness variations? Before addressing these questions, we must be able to compute tidal dissipation in a body made of a viscoelastic core, a global ocean, and a viscoelastic shell of non-uniform thickness. Here I will solve this problem with the non-uniform thin shell approach of Beuthe [2018] (hereafter named Paper I). The approach will be benchmarked against the solution for dissipation in a laterally uniform thin shell, which has not appeared before in the literature (dissipation in a laterally uniform membrane was studied in Beuthe [2015a]). Besides being a non-trivial check of the consistency of the thin shell approach, this step provides an estimate of the error due to the thin shell approximation.

Without prejudging whether Enceladus is truly in a steady state, it seems natural to apply first the method to the simplest case, i.e. thermal equilibrium between tidal heating and conductive cooling. As a concrete example, I will assume that Enceladus’s shell varies in thickness according to the isostatic interpretation of gravity-topography data [Beuthe et al., 2016]. Given that previous studies made it clear that dissipation within the shell does not provide enough power at the present time, I will push the envelope of allowable parameters [Wolfe, 1979] in order to address the following questions:

• •

Does the variation of shell thickness increase the total dissipated power?

• •

How much does shell thinning increase the dissipation flux at the south pole?

• •

Is there a general relation between local shell thickness and local tidal dissipation?

• •

Can we tune core and shell parameters so as to match both the total conductive flux and its spatial distribution, assuming thermal equilibrium?

• •

Does the non-uniform shell induce a north-south asymmetry in core dissipation?

In a previous study of the topic, Běhounková et al. [2017] used the finite-element method (FEM) of Souček et al. [2016] to solve for the deformations of an elastic shell with non-uniform thickness, and predicted several tens of GW of tidal heating by assuming a shell of uniform viscosity. This last assumption is not realistic for a conductive shell because viscosity increases by orders of magnitude from the bottom of the shell to the surface; the total power dissipated in the shell is actually much lower. Older 3D studies (reviewed in Běhounková et al. [2017]) investigated the effect of laterally varying rheology on dissipation within a very thick convecting shell, now disfavoured by Cassini data. Souček et al. [2019] recently improved the FEM approach by taking into account temperature-dependent viscosity and now predict less than $2\rm\,GW$ of dissipation in Enceladus’s shell. The FEM approach has the advantages of including faults (modelled as frictionless open slots) and of avoiding the errors intrinsic to the thin shell approach (mainly the thin shell approximation and the non-zero bulk dissipation). Its disadvantages are that it does not include core dissipation; that it becomes unstable at high viscosity contrasts (need of viscosity cutoff); that it is not readily adaptable to non-Maxwell rheological models (e.g. the pseudo-Andrade rheology in Souček et al. [2019]); and that it neglects self-gravity (this approximation is justified for Enceladus but is not valid for large satellites). Preliminary benchmarking against the thin shell approach shows that the two methods agree within a few percent error [Běhounková et al., 2018]. Thus, the technical differences mentioned above are not problematic. The choice of the method rather depends on whether one wants to include faulting or core dissipation in the model.

The rest of the paper is made of four parts: (1) summary of the tidal thin shell equations from Paper I, (2) thin shell approach to dissipation in a non-uniform shell and in a spherically symmetric core, (3) benchmark against a laterally uniform thick shell, and (4) tidal heating in a body with a dissipative core and a non-uniform conductive shell in thermal equilibrium.

2 Tidal thin shell equations

2.1 Flexure equations

In Paper I, I used thin shell theory to compute the tidal deformations of an icy satellite with a viscoelastic core, a global ocean and a laterally non-uniform viscoelastic shell. This approach rests on the following assumptions about the shell: thickness less than 10 to 20% of the surface radius, uniform Poisson’s ratio, homogeneous density, no density contrast with the top layer of the ocean, linear viscoelasticity. Enceladus’s internal structure approximately satisfies these criteria (Table 1). In addition, the core structure must be spherically symmetric.

The thin shell approach to tidal deformations is neatly summarized by the tidal thin shell equations. These two differential equations govern the tidal deformations of a non-uniform thin spherical shell floating on a global ocean:

[TABLE]

where

• •

$w$ is the radial displacement and $F$ is the auxiliary stress function;

• •

$\alpha$ and $D$ are the primary viscoelastic shell parameters (Table 2); $\chi$ is a secondary viscoelastic shell parameter, close to one (Table 2);

• •

${\cal C}$ and ${\cal A}$ are spherical differential operators of order 4 (Appendix A).

In a spherical harmonic basis, the tidal thin shell equations become a system of coupled linear equations which can be solved for $(F,w)$ with standard matrix methods. Other variables, such as tangential displacements, strains, and stresses, can be computed from $(F,w)$ with the pseudospectral transform method: angular derivatives are applied in the spectral domain whereas products of fields are computed in the spatial domain (see Paper I).

2.2 Tidal forcing

The RHS of Eq. (1) represents the tidal loading which can be expressed as a sum of spherical harmonics of degree $n$:

[TABLE]

where $q_{0}$ is the degree-0 load (see Section 3.4 of Paper I) while the degree-1 load vanishes. If $n\geq 2$, the degree-$n$ load is given by

[TABLE]

where $\rho{}gw_{n}$ is the weight of the tidal bulge in the unperturbed gravity field. The total perturbing potential $\Gamma_{n}$ depends on the primary tidal potential $U_{n}^{T}$ and on the secondary potential induced by the deformation of the whole body. In the thin shell approach, it is given by

[TABLE]

where $\xi_{n}$ is the degree-$n$ density ratio,

[TABLE]

The spherically symmetric structure below the shell (stratified ocean, viscoelastic core) enters into the problem via the nondimensional factor $\upsilon_{n}$ ($\upsilon_{n}\geq 1$),

[TABLE]

where $h_{n}^{\circ}$ is the tidal radial Love number of the fluid-crust body (i.e. the same body except that the icy shell behaves as a fluid). If the core is not deformable, $h_{n}^{\circ}=(1-\xi_{n})^{-1}$ and $\upsilon_{n}=1$. In that case, the term $\xi_{n}w_{n}$ in $\Gamma_{n}$ (Eq. (4)) can be interpreted as the geoid perturbation due to the tidal deformation (self-attraction or self-gravity). In the simplest model with a deformable core, the core is viscoelastic, homogeneous, and incompressible, and the ocean is homogeneous, in which case $h_{n}^{\circ}$ can be computed with Eq. (D.2).

Since Enceladus is in a synchronous orbit with negligible obliquity, tidal deformations are mainly due to eccentricity tides of degree 2. These tides can be expressed as the sum of a radial tide (due to the varying distance to Saturn) and a librational tide (due to the optical libration, i.e. the varying direction of Saturn in the frame fixed to the satellite) [Murray and Dermott, 1999]. The latter is enhanced by the 1:1 forced (or physical) libration. Including degree-2 eccentricity tides plus the forced libration (denoted $\gamma=-\gamma_{0}\sin\omega t$), I write the tidal potential as [Van Hoolst et al., 2013]

[TABLE]

where $\gamma_{0}=0.12^{\circ}$ for Enceladus [Thomas et al., 2016]. $P_{2m}$ are the associated Legendre functions of degree $2$ and order $m$ depending on $\cos\theta$ ($\theta$ is the colatitude and $\varphi$ is the longitude). The amplitudes of the optical and forced librations must be added because the forced libration and the tidal torque are $180^{\circ}$ out of phase (e.g. Fig. 5 of Hemingway et al. [2018]). Thus, including the 1:1 forced libration increases tidal dissipation (by about 28%, see Eq. (49)), as already shown for a homogeneous body [Wisdom, 2004], instead of decreasing it as concluded by Běhounková et al. [2017].

Introducing the forced libration into the tidal potential, as above, makes sense for a completely solid body, in which the differential rotation between layers is negligible. This procedure is however problematic if the shell and the solid core are decoupled by a global ocean. In that case, the core has a smaller libration than the shell (by a factor of 10 in amplitude, A. Trinh, private comm.), so that the tidal potential (7) is only valid for the shell. Moreover, the differential rotation of the core and shell induces gravitational and pressure couplings between the shell and core which should be taken into account in an additional forcing term. Here I assume that the shell and core have no differential rotation so that they are forced by the same tidal potential (7). For the shell, this approximation is reasonable as long as the forced libration is significantly smaller than the optical libration: $\gamma_{0}/2e=0.223$ implies that the neglected corrections are of order $(\gamma_{0}/2e)^{2}\ll 1$. For the core, overestimating dissipation does not pose a problem because core rheology is adjusted so that the total heat budget is satisfied.

3 Dissipation in the thin shell approach

In this section, I set forth the full methodology required to compute dissipation in the non-uniform shell and the internally spherically symmetric core. First, I explain how the assumption of a uniform Poisson’s ratio affects dissipation. Next, I derive formulas for dissipation in the shell (rate, flux, and power) in terms of the basic thin shell variables $(F,w)$. Finally, I show how to compute dissipation in the core by the way of the effective tidal potential. This method is applied to partition the total power into core and shell contributions and the compute the spatially dependent core dissipation rate.

3.1 Shear, bulk and Poisson dissipation

In the micro approach, dissipation is computed locally by multiplying at each point the microscopic stress by the strain rate. As Enceladus is rotating synchronously with its mean motion, tidal deformations due to the eccentric orbit are periodic with an angular frequency $\omega$ equal to the mean motion. It is thus convenient to work with Fourier-transformed variables: $V(t)={\rm Re}(\tilde{V}(\omega)\,e^{i\omega t})$. The dissipation power per unit volume averaged over one orbital period $T$, in short the dissipation rate, reads (Appendix A of Beuthe [2013])

[TABLE]

where the asterisk denotes complex conjugation. The tensors $\sigma_{ij}$ and $\epsilon_{ij}$ denote stress and strain in the time domain (without tilde) or frequency domain (with tilde). Henceforth I work with frequency-domain variables which depend implicitly on $\omega$, and I drop the ‘tilde’ notation. According to the correspondence principle, linear viscoelasticity is introduced through complex moduli (shear modulus $\mu$ and bulk modulus $K$) in the frequency-domain stress-strain relation:

[TABLE]

where $\epsilon=\epsilon_{rr}+\epsilon_{\theta\theta}+\epsilon_{\varphi\varphi}$ denotes the trace of the 3D strain tensor and $\delta_{ij}$ is the Kronecker delta. If Eq. (9) holds, the dissipation rate becomes

[TABLE]

The term proportional to ${\rm Im}(\mu)$ represents dissipation resulting from shear friction. The term proportional to ${\rm Im}(K)$ represents bulk dissipation, which is poorly constrained although seismic data suggest that it is much smaller than shear dissipation, at least at seismic frequencies [Durek and Ekström, 1995; Resovsky et al., 2005]. Thus, one usually assumes the condition of ‘no bulk dissipation’: ${\rm Im}(K)=0$ or equivalently $K=K_{e}$. This assumption has been criticized on theoretical grounds [Morozov, 2015] and might be invalid in the presence of melt [Takei and Holtzman, 2009]. Recently, Ricard et al. [2014] argued that seismic attenuation could be due to the laminated structure of the mantle, in which case intrinsic dissipation would not be constrained at all. The occurrence of bulk dissipation is thus an open question. Here, I cannot impose that $K=K_{e}$ because it would require that Poisson’s ratio $\nu$ varies in tandem with the shear modulus $\mu$ (e.g. Appendix C of Beuthe [2015a]) according to the $K$-$\mu$-$\nu$ relation:

[TABLE]

This would contradict the assumption of uniform $\nu$ made in the theory of non-uniform thin shells (see Section 2). I impose instead the condition of ‘no Poisson dissipation’, i.e. that Poisson’s ratio remains equal to its elastic value which is real and uniform: $\nu=\nu_{\rm e}$. This constraint together with Eq. (11) implies that

[TABLE]

If there is no Poisson dissipation, the dissipation rate thus reads

[TABLE]

If there is no bulk dissipation, the factor $\nu/(1-2\nu)$ should be replaced by $-1/3$. In the incompressible limit ($K\rightarrow\infty$, $\nu\rightarrow 1/2$), bulk dissipation always vanishes so that the conditions of ‘no bulk dissipation’ and ‘no Poisson dissipation’ are simultaneously satisfied. In that case, ${\rm Tr}\,\epsilon_{ij}=0$ and the dissipation rate reduces to the first term of Eq. (13).

Before going further, I need to specify the rheology. The simplest linear viscoelastic model is the one of Maxwell (see Table 2). In this model, ${\rm Im}(\mu)$ is maximum at the forcing angular frequency $\omega$ if the viscosity of the material is equal to the critical viscosity $\eta_{\rm crit}=\mu_{\rm e}/\omega$, where $\mu_{e}$ is the elastic shear modulus. For Enceladus’s icy shell, the critical viscosity is equal to $6.6\times 10^{13}\rm\,Pa.s$. But a conductive shell is not at all homogeneous: the viscosity, which controls ${\rm Im}(\mu)$, varies by orders of magnitude between the cold surface and the bottom of the shell at the melting point (the strain varies much less because it is controlled by the upper and colder part of the shell). Thus, the dissipation rate is maximum in a thin layer where viscosity is closest to the critical viscosity. By the same logic, the total power dissipated in the shell is maximum if the viscosity at the bottom of the shell is lower than the critical viscosity. Since I am interested to maximize tidal dissipation within the shell, I generally assume that the bottom viscosity is equal to $10^{13}\rm\,Pa.s$, which the lower bound for the viscosity of ice in the low stress regime [Tobie et al., 2003; Barr and Showman, 2009]. Nonetheless, I will also consider higher viscosities for benchmarking. Although Andrade rheology (e.g. Castillo-Rogez et al. [2011]) is more realistic than Maxwell and more dissipative at high viscosity, it makes little difference if the bottom viscosity is lower than the critical viscosity (see example in Section 4.3). There is no problem, however, to switch in the thin shell approach from Maxwell to Andrade rheology.

Fig. 1 shows how the dissipation condition affects the total power dissipated in a thick shell (‘thick’ meaning no thin shell approximation) which is laterally uniform and conductive (computational details are given in Section 4.1). Instead of the power itself, the figure shows the imaginary part of the gravitational Love number $k_{2}$, which has the advantage that the dependence on the tidal potential (eccentricity and libration) is factored out (Eq. (48)). If the shell thickness $d$ is in the range $20-50\rm\,km$, assuming ‘no Poisson dissipation’ changes the total power by less than 2% whereas assuming an incompressible shell changes it by -5 to -7% (the precise number depends on the bottom viscosity). For thinner shells, the effect of the ‘no Poisson dissipation’ condition depends a lot on the bottom viscosity: it is below 1% if $\eta_{\rm m}=10^{13}\rm\,Pa.s$, but may reach 9% in the membrane limit if $\eta_{\rm m}=10^{15}\rm\,Pa.s$. In comparison, the incompressible assumption changes the total power by 9 to 19% in the membrane limit. In conclusion, assuming an incompressible shell, as is done in the propagator matrix approach, generally leads to a larger error than the condition of ‘no Poisson dissipation’.

3.2 Dissipation inside the shell

3.2.1 Shell dissipation rate

Apart from the assumption of no Poisson dissipation, the dissipation rate given by Eq. (13) is completely general and valid anywhere in the body. I will now restrict it to the shell and work in the thin shell limit. The plane stress approximation underlying thin shell theory implies that (Eq. (C.1) of Paper I)

[TABLE]

Substituting this constraint into Eq. (13), I write the shell dissipation rate in absence of Poisson dissipation as

[TABLE]

where $({\cal E}_{2},{\cal E}_{\rm tr})$ are 2D strain invariants:

[TABLE]

If there is no bulk dissipation, the factor $\nu/(1-\nu)$ in Eq. (15) should be replaced by the factor $(4{\rm Re}(\nu)-1-|\nu|^{2})/(3|1-\nu|^{2})$; both factors tend to 1 in the incompressible limit.

If the tidal thin shell equations are solved for $(F,w)$, the shell dissipation rate can be written as a bilinear form in these variables and their complex conjugates (see Appendix B):

[TABLE]

If there is no Poisson dissipation, the terms in the RHS are given by

[TABLE]

The shell dissipation rate depends on depth through ${\rm Im}(\mu)$ and the depth parameter $z$ (Table 2). It does not depend on the choice of reference surface, since $(w,\alpha{}F,z/R)$ are invariant under a change of $R$ (see Appendix G of Paper I and Eq. (J.3) of Paper I).

3.2.2 Shell dissipation flux

The shell dissipation flux is the energy flux, due to dissipation within the shell, through the reference surface of the shell (chosen here to be the outer surface of the body). If the heat transfer is radial, the shell dissipation flux is equal to the dissipation rate integrated over the shell thickness,

[TABLE]

where the factor $(1+\zeta/R)^{2}$ is associated with the integration measure in spherical coordinates. The shell dissipation flux can be expressed in terms of the variables $(F,w)$ and of the parameters $(\nu,\alpha,D,\chi)$ (see Appendix C):

[TABLE]

If there is no Poisson dissipation, the terms in the RHS are given by

[TABLE]

The subscripts mem, mix, bend stand for membrane, mixed and bending contributions. This is analogous to the decomposition of the elastic energy in extensional-shearing, mixed, and bending-twisting terms [Novozhilov, 1964; Axelrad, 1987]. Membrane and bending contributions are always positive whereas the mixed contribution can be negative.

The shell dissipation flux has several nice properties:

• •

It depends on rheology through depth-integrated shell parameters: ${\rm Im}(\alpha)$, ${\rm Im}(\chi)$, ${\rm Im}(D)$.

• •

It depends on scalar quantities and scalar differential operators, making it easy to compute it with the pseudospectral transform method.

• •

It does not depend on the degree-1 spherical harmonic components of $w$ (invariance under rigid displacements) and of $F$ (degree-1 gauge freedom), which belong to the null space of $\Delta^{\prime}$ and ${\cal A}$.

• •

It is inversely proportional to the area of the reference surface of the shell (as it should be), because $w$ is invariant under a change of $R$ while other quantities scale as $F\sim{}R^{-2}$, $\alpha\sim{}D\sim{}R^{2}$, and $\chi\sim{}R$ (see Appendix G of Paper I).

• •

It satisfies the static-geometric duality (Eq. (17) of Paper I), a transformation exchanging the LHS of the governing equations: $(w,D,\nu,\chi)\leftrightarrow(R^{2}F,-\alpha,-\nu,\chi)$.

3.2.3 Shell power

The total power dissipated in the shell, or shell power, is obtained by integrating the shell dissipation flux over the reference surface with surface element $dS=R^{2}\sin\theta\,d\theta\,d\varphi$:

[TABLE]

In general, this integral must be evaluated numerically, but it can be done analytically if the shell is laterally uniform (see Appendix E). Eq. (23) embodies the micro approach to tidal dissipation, in which the total power dissipated in the shell is evaluated by integrating the microscopic dissipation rate over the volume of the shell.

An alternative approach to tidal dissipation (macro approach) consists in computing the total power dissipated in the body from the work done by the tidal potential (Zschau-Platzman formula, see Eq. (7) of Platzman [1984]):

[TABLE]

where $U_{n}^{\prime}$ is the secondary potential due to the deformation of the body:

[TABLE]

in which $\Gamma_{n}$ is the total perturbing potential (Eq. (3)). If the primary tidal potential is of degree two, substituting Eq. (4) into Eq. (24) yields the total power in terms of $(U_{2}^{T},w_{2})$:

[TABLE]

where the bracket notation $\langle x\rangle$ denotes the angular average of $x$ (or degree-0 spherical harmonic coefficient of $x$). Once the tidal thin shell equations have been solved for $(F,w)$, Eq. (26) immediately yields the total dissipated power.

If the core is elastic (${\rm Im}(\upsilon_{2})=0$) and the ocean is inviscid, the total dissipated power is equal to the shell power. Setting ${\rm Im}(\upsilon_{2})=0$ into Eq. (26) and using Eqs. (3)-(4), I can write

[TABLE]

which can be interpreted as follows: $\dot{E}_{shell}$ is equal to the dissipative part of the power developed by the bottom load acting on the shell (see Appendix E of Beuthe [2015a]). Dissipation only occurs at degree 2 because $(q_{n},w_{n})$ are in phase if $\upsilon_{n}$ is real and $U_{n}^{T}=0$ (see Eqs. (3)-(4)). Eq. (27) provides a non-trivial check on the integrated dissipation rate given by Eq. (23).

If the core is viscoelastic, it is tempting – but not correct – to interpret the two terms of the RHS of Eq. (26) as core and shell contributions. In Section 3.3.2, I will explain how to split $\dot{E}_{tot}$ between core and shell using the effective tidal potential.

3.3 Dissipation inside the core

3.3.1 Effective tidal potential for the core

Another source of dissipation arises from the tidal deformations of the viscoelastic core. If the whole body has a spherically symmetric structure, one first solves the (radial) viscoelastic-gravitational equations from the center to the surface of the body before evaluating the dissipation rate at each point. Here, the laterally non-uniform shell breaks spherical symmetry. Nevertheless, radial viscoelastic-gravitational equations are still valid in the core (and ocean) as long as the interior structure is spherically symmetric under the shell. The effect of the viscoelastic shell can be represented as a pressure load

[TABLE]

which acts on the associated fluid-crust body (defined as the body equivalent to the original one except that the shell has no rigidity; see Section 3.3 of Paper I). Besides this pressure load, tides deform the fluid-crust body (Eq. (26) of Paper I).

The radial viscoelastic-gravitational equations should now be solved in a 2-layer body (core plus surface ocean) submitted to tidal and pressure loads. Beware that ‘2-layer’ does not imply here uniform layers: density and rheology can be depth-dependent. In the standard approach [Takeuchi and Saito, 1972], the displacements, stresses and gravitational perturbations due to a forcing of degree $n$ are represented in a solid layer by 6 radial response functions $y_{in}(r)$ ($i=1...6$). Since the tidal and pressure load solutions correspond to different boundary conditions, the full solution reads

[TABLE]

where the tidal potential is evaluated at radius $R$. The variables $y_{in}^{\circ J}$ (J=T,P) are propagated from the center (where they depend on 3 unknown constants) to the surface of the solid core by solving 6 viscoelastic-gravitational differential equations. At the core-ocean boundary, a first boundary condition is given by the condition of zero shear stress ($y_{4n}^{\circ J}=0$). In the limit of static deformations [Saito, 1974], a second homogeneous boundary condition is provided by the fluid constraint (e.g. Beuthe [2015b]), which relates $y_{1n}^{\circ J}$ (radial displacement), $y_{2n}^{\circ J}$ (radial stress) and $y_{5n}^{\circ J}$ (gravitational potential perturbation). As there is no third boundary condition at the core-ocean boundary, one must go on and solve the differential equations in the fluid layer. Under the assumption of static deformations, the gravitational potential decouples from displacements: it becomes sufficient to propagate the variables $y_{5n}^{\circ J}$ and $y_{7n}^{\circ J}=y_{6n}^{\circ J}+(4\pi{}G/g)y_{2n}^{\circ J}$ [Saito, 1974]. At the outer fluid surface, there is only one (inhomogeneous) boundary condition given by

[TABLE]

which result from the usual boundary conditions for tidal forcing and pressure loading (Eqs. (C.5) and (E.3) of Beuthe [2016]). Thus, the pressure and tidal load solutions are related within the core by $y_{in}^{\circ P}=-y_{in}^{\circ T}$ ($i=1...6$).

Defining the effective tidal potential for the core by

[TABLE]

I can write the full solution in the core (Eq. (29)) as the tidal solution forced by $U_{n}^{\circ}$:

[TABLE]

The similar property $y_{in}^{\circ P}=-y_{in}^{\circ T}$ ($i=5,7$) holding within the fluid leads to relations between pressure and tidal Love numbers: $k_{n}^{\circ P}=-h_{n}^{\circ}$ and $h_{n}^{\circ P}=-1/\xi_{n}-h_{n}^{\circ}$, which were already noted as Eqs. (29)-(30) of Paper I (the superscript $T$ is omitted for tidal Love numbers).

For evaluation purposes, it is practical to write the effective tidal potential in terms of $U_{n}^{T}$ and the flexure solution $w_{n}$:

[TABLE]

This formula gives us a preliminary estimate of how much lateral variations of the shell structure influence the core dissipation pattern. If the shell is laterally uniform, the radial displacement is related to the tidal potential by $w_{n}=h_{n}U_{n}^{T}/g$ (see Section 4.1), so that the second term in the brackets of Eq. (34) is small with respect to the first: $\xi_{n}\,h_{n}\,U_{n}^{T}\ll{}U_{n}^{T}$ ($\xi_{2}h_{2}\sim 10^{-2}$, see Section 4.1). Lateral variations in the shell result in small deviations from this prediction. The maximum deviation occurs at the south pole, where the radial displacement at the south pole changes from $-0.5\rm\,m$ to $-0.7\rm\,m$ (see Fig. 8 of Paper I). Thus, lateral variations in the shell change the effective tidal potential by less than 1% and cause negligible deviations in the core dissipation pattern with respect to the case of a laterally uniform shell.

3.3.2 Core-shell partition of total power

In Section 3.2.3, I used the macro approach in order to compute the total power dissipated in the body, and showed that it is equal to the shell power if the core is elastic. If the core is viscoelastic, the total power can be split into core and shell contributions with the help of the effective tidal potential. Note first that the total perturbing potential can be written as

[TABLE]

which results from Eqs. (3) and (28). Besides, $\Gamma_{n}$ is related to the effective tidal potential for the core by

[TABLE]

which results from Eqs. (4) and (34).

Using Eqs. (32), (35), and (36), I write the integrand of Eq. (24) as

[TABLE]

The total power dissipated in the body thus reads

[TABLE]

where

[TABLE]

Similarly to Eq. (27), $\dot{E}_{shell}$ can be identified as the dissipative part of the power developed by the bottom load acting on the massless shell. It is thus equal to the shell power, the difference with Eq. (27) being that dissipation now occurs at all harmonic degrees. If the shell is elastic, $\dot{E}_{shell}$ vanishes because the tidal thin shell equations impose that $(q_{n},w_{n})$ are in phase.

If $\dot{E}_{shell}$ is the shell power in Eq. (40), $\dot{E}_{core}$ must be equal by subtraction to the power dissipated in the core. This claim is confirmed in all generality by integrating the dissipation rate over the volume below the shell, using energy conservation (Eq. (179) of Takeuchi and Saito [1972]) and the surface boundary conditions for tidal and pressure loading solutions. As a caveat, beware that the above core-shell partition relies on the assumption of no differential rotation between shell and core (see Section 2.2).

For a laterally varying shell, $\dot{E}_{core}$ and $\dot{E}_{shell}$ can be expressed in terms of $(w_{n},U_{n}^{T})$ by substituting Eq. (3)-(4) and Eq. (34) into Eqs. (39)-(40). Once the tidal thin shell equations have been solved, it is straightforward to evaluate the core and shell contributions to the total dissipated power. Comparing the numerical values yielded by the micro and macro formulas (Eqs. (23) and (40)) for the shell power provides a self-consistency check for numerical codes.

3.3.3 Core dissipation rate

Spatial patterns of core dissipation can only be computed in the micro approach. If the forcing is of degree 2 (or more generally of a given degree $n$), the dissipation rate in a spherically symmetric layer can be factorized into radial and angular parts which depend on the internal structure and on the square of the forcing potential, respectively (see Beuthe [2013]). The non-uniformity of the shell, however, introduces harmonic degrees other than degree 2 in the forcing potential (Eq. (32)), which interfere when computing the squared forcing potential. Thus, it is not sufficient to consider the patterns corresponding a tidal forcing of a given harmonic degree. Nonetheless, it is advantageous to express the dissipation rate in terms of spherical differential operators acting on scalar radial functions, instead of stresses and strains, because derivatives of scalar fields can be computed efficiently with the pseudospectral transform method.

Following Section 2.2 of Beuthe [2013], I write the dissipation rate within the incompressible core as

[TABLE]

where

[TABLE]

$E_{A}$, $E_{B}$, and $E_{C}$ correspond respectively to $E_{radial}+E_{tan2}-E_{dilat}/3$, $E_{shear}$, and $E_{tan1}$ in Beuthe [2013]. The differential operators ${\cal D}_{2}$ and ${\cal D}_{4}$ are defined in Appendix A. The functions $\Upsilon_{i}$ are defined by Eqs. (33)-(34): $\Upsilon_{1}$ is the radial displacement, $\Upsilon_{3}$ is the potential for tangential displacement, and $\Upsilon_{4}$ is the potential for radial-tangential stress (or strain). If the core and the ocean are both homogeneous, the functions $y_{in}^{T}$ appearing in $\Upsilon_{i}$ are given by Eq. (D.5).

As the shell is laterally non-uniform, the functions $\Upsilon_{i}$ are superpositions of spherical harmonics of different degrees. In that case, the various terms appearing in the RHS of Eq. (42) can be evaluated with the transform method mentioned in Section 2.1. Derivatives only appear through the spherical Laplacian $\Delta$ (or the related operator $\Delta^{\prime}$), as can be seen for the operators ${\cal D}_{2}$ and ${\cal D}_{4}$ with identities (b), (c), and (e) of Appendix A.

Under the assumption of radial heat transfer, the heat flux at the core-ocean boundary is equal to the dissipation rate integrated over the core radius,

[TABLE]

The total power dissipated in the core is obtained by integrating the flux over the core surface with surface element $dS_{c}=R_{c}^{2}\sin\theta\,d\theta\,d\varphi$:

[TABLE]

This expression should be equal to the power obtained in the macro approach (Eq. (39)). If the core and the ocean are both homogeneous, this integral can be done analytically with the solution of Appendix D, yielding Eq. (39) in which $k_{n}^{\circ}$ is given by Eq. (D.2). Note that interfering harmonic degrees do not contribute to the angular integral. The equivalence between Eqs. (39) and (44) is another example of the consistency between micro and macro approaches to tidal dissipation.

4 Benchmarking against a laterally uniform thick shell

In this section, I benchmark the thin shell solution for tidal dissipation against the thick shell solution for a laterally uniform shell. The primary aim is to quantify the impact of the thin shell approximation on dissipation patterns and on the total power dissipated in the shell and core. Dissipation patterns will be analyzed with the radial-angular factorization method. Regarding the total power, it can be expressed in terms of the imaginary part of the Love number $k_{2}$, making it easy to study the error due to the thin shell approximation.

4.1 Factorization of dissipation rate

For this benchmark, Enceladus is modelled as a 3-layer body made of a homogeneous core, a homogeneous ocean, and a conductive icy shell. The internal structure is spherically symmetric. Model parameters are given in Table 1. The rheology of the shell is modelled as in Section 4.2.3 of Paper I: the rheology is described with the Maxwell model; the viscosity depends on temperature through an Arrhenius relation (Table 2); the temperature profile is the solution of the Cartesian 1D heat equation without internal source (Eq. (56) of Paper I); the surface temperature is uniform and set to $59\rm\,K$.

Before computing dissipation, one should solve for tidal deformations. In the thin shell approach, the tidal thin shell equations can be solved analytically in the spherical harmonic basis (see Table 3). The exact solution for a thick shell is obtained by integrating the elastic-gravitational equations for the spherically symmetric problem in the static limit [Takeuchi and Saito, 1972; Saito, 1974], as done for the core in Section 3.3.1. The solution is a set of six radial functions $y_{i}$, three of which are needed here: $y_{1}$ (radial displacement), $y_{3}$ (potential for tangential displacement), and $y_{4}$ (potential for radial-tangential stress or strain).

If the internal structure is spherically symmetric, the dissipation rate can be factorized [Beuthe, 2013]:

[TABLE]

where $f_{J}$ are the radial weight functions while $\Psi_{J}$ are angular functions representing the basic spatial patterns. The ratios $(f_{A}+f_{K},f_{B},f_{C})/f_{T}$ (with $f_{T}=\sum_{J}f_{J}$) measure the contributions of the patterns to the angular average of the dissipation rate. Table 4 gives the weight functions in terms of the viscoelastic-gravitational solutions $y_{i}$, and specifies the angular functions for degree-$2$ eccentricity tides combined with the forced libration (see Eq. (7) and Appendix D of Beuthe [2013]).

In a similar fashion, the shell surface flux can be expressed as a weighted sum of the angular functions $\Psi_{J}$:

[TABLE]

where

[TABLE]

The ratios $({\cal F}_{A}+{\cal F}_{K},{\cal F}_{B},{\cal F}_{C})/{\cal F}_{T}$ (with ${\cal F}_{T}=\sum_{J}{\cal F}_{J}$) measure the contributions of the patterns to the average flux (or to the total power).

Fig. 2 shows the basic dissipation patterns $(\Psi_{A},\Psi_{B},\Psi_{C})$ for degree-2 eccentricity tides combined with forced libration ($\gamma_{0}=0.12^{\circ}$ or $\gamma_{1}=1.223$). These patterns do not differ much from the patterns without forced libration (see Fig. 1 of Beuthe [2013]; beware that this figure shows patterns with zero mean and unit standard deviation). $\Psi_{A}$ includes significant contributions from both degrees 2 and 4, has maxima along the equator and minima at middle latitudes. $\Psi_{B}$ does not contribute in the thin shell approximation (but it contributes to core dissipation). $\Psi_{C}$ is dominated by the degree-2 component of the squared tidal potential, is maximum at the poles, and zero along the tidal axis.

4.2 Patterns in a laterally uniform thin shell

If the shell is laterally uniform, the thin shell dissipation rate must be factorizable in radial weights and angular functions. Thin shell assumptions (i.e. purely tangential stress) imply that only two patterns contribute: $\Psi_{A}$ and $\Psi_{C}$. In Appendix E, I substitute the analytical solutions for $(F,w)$ given in Table 3 into the general expression of the thin shell dissipation rate (Eqs. (18)-(19)). The resulting formulas for the radial weights $f_{J}$ and ${\cal F}_{J}$ are given by Eqs. (E.4) and (E.6).

Regarding the dissipation rate, Fig. 3 shows that the thin shell approximation is very accurate for the dominant weight function $f_{C}$, with or without compressibility. It is less accurate for the coefficient of Pattern A: the difference is in part due to the error on the weight function $f_{A}$, and in part to bulk dissipation (if the shell is compressible). However, the error is very small near the bottom of the shell where dissipation is highest.

Regarding the surface flux pattern, Table 5 gives the contribution of the dominant Pattern C to the average surface flux for different types of solutions (thick shell/thin shell/membrane), and for different dissipation conditions. The thin shell approach predicts much better the dissipation pattern than the membrane approximation: the error is less than 1% if $\eta_{\rm m}=10^{13}\rm\,Pa.s$ (or less than 3% if $\eta_{\rm m}=10^{14}\rm\,Pa.s$). For thinner shells, the contribution of Pattern C decreases, reaching a lower bound of $75\%$ to $82\%$ in the membrane limit (more details are given in Appendix E). Fig. 4 shows the surface flux pattern predicted in the membrane approximation and with the thin shell approach (the latter is indistinguishable from the thick shell solution if $\eta_{\rm m}=10^{13}\rm\,Pa.s$). All in all, the surface flux pattern remains similar as shell thickness and bottom viscosity vary, but the dissipation contrast between poles and tidal axis is much stronger if the shell is rather thick. The min/max dissipation contrast in a thin shell is much higher than in a solid body (or in the core, see Section 4.4), for which the surface flux is about twice as large at the poles than along the tidal axis.

4.3 Total power in core and shell

In the macro approach, the total dissipated power is given by a surface integral, the integrand of which is the product of the primary and secondary tidal potentials (Eq. (24)). Substituting $\Gamma_{n}=(k_{n}+1)U_{n}^{T}$ (Table 3) in the Zschau-Platzman formula, I get back the well-known formula for the total power dissipated in a body with a spherically symmetric structure:

[TABLE]

The effect of the internal structure is hidden in the imaginary part of the gravitational Love number, while the external forcing appears as the averaged squared tidal potential [Zschau, 1978; Platzman, 1984]. For degree-2 eccentricity tides combined with the forced libration, the latter is given by (see Table 4)

[TABLE]

The amplification due to the forced libration is the same as for a homogeneous body [Wisdom, 2004], because of the (arbitrary) assumption that the core and shell librate with the same angle.

If the core is elastic, the total power $\dot{E}_{tot}$ is equal to the shell power $\dot{E}_{shell}$. Eq. (49) shows that the forced libration ($\gamma_{0}=0.12^{\circ}$) increases the total power dissipated in the shell by 28%. Fig. 5 shows the total power dissipated in a laterally uniform shell as a function of the shell thickness and the viscosity at the bottom of the conductive shell. Rheology is either Maxwell or Andrade, the latter with parameters equal to $\alpha_{A}=0.25$ and $\beta_{A}\approx\mu_{\rm e}^{\alpha_{A}-1}\eta^{-\alpha_{A}}$ [Castillo-Rogez et al., 2011]. If the shell is thicker than a few km, dissipation increases as the bottom viscosity decreases from $10^{16}$ to $10^{13}\rm\,Pas.s$. Lowering the bottom viscosity below that value does not increase dissipation much further. Andrade rheology leads to more dissipation than Maxwell if the bottom viscosity is larger than $10^{13}\rm\,Pa.s$, but makes little difference below that threshold.

In Fig. 6, I compare the thin shell approximation of ${\rm Im}(k_{2})$ with the exact thick shell results for different values of the shell thickness and bottom viscosity. Fig. 6A confirms that the thin shell result is a well-behaved approximation of the exact result: the error with respect to the thick shell power with no Poisson dissipation is below 4% if $d<50\rm\,km$ (or below 2% if $d<20\rm\,km$) and tends to zero with decreasing shell thickness. Fig. 6B shows that the error with respect to the thick shell power with no bulk dissipation is below 4% if $20<d<50\rm\,km$. For thinner shells, the error increases again if the bottom viscosity is higher than $10^{13}\rm\,Pa.s$, but the error is approximately bounded by the curve for $\eta_{\rm m}=10^{15}\rm\,Pa.s$. In the membrane limit ($d\rightarrow 0\rm\,km$) the error ranges from 0 to 9%, depending on the bottom viscosity (the membrane error is exactly $(1-2\bar{\kappa})/(11+2\bar{\kappa})$ where $\bar{\kappa}$ is the effective bulk dissipation varying between 0 and 0.5; see Eq. (E.7)).

Suppose now that the core is viscoelastic. The partition of the total power (Eqs. (39)-(40)) becomes, after substituting the expressions for $(w_{n},q_{n},U_{n}^{\circ})$ from Table 3,

[TABLE]

In the membrane limit ($\Lambda_{n}\rightarrow\Lambda_{n}^{M}$), these formulas are identical to Eqs. (98)-(102) of Beuthe [2015a]. All quantities in the RHS can be numerically evaluated. In particular, the tidal Love numbers $(k_{n}^{\circ},k_{n},h_{n})$ and the thin shell spring constant $\Lambda_{n}$ can be computed with the analytical formulas of Table 3, which remain valid if the core is viscoelastic (see Section 4.4).

The macro formula for $\dot{E}_{shell}$ (Eq. (51)) is identical to the volume-integrated dissipation rate (Eq. (E.11)): this is again an example of the micro-macro equivalence in tidal dissipation [Beuthe, 2013, 2015a]. The substitutions $\Lambda_{n}\rightarrow\Lambda_{n}^{M}$ and $\Lambda_{n}\rightarrow\Lambda_{n}^{B}$ in Eq. (51) yield to first order the sum of the membrane and mixed contributions, and the bending contribution, respectively (Eq. (E.11) again). The bending contribution is about 5% of the total power for a laterally uniform shell with $d=23\rm\,km$ and $\eta_{\rm m}=10^{13}\rm\,Pa.s$ (it increases if the shell is laterally non-uniform, see Section 5.5). Finally, energy conservation ($\dot{E}_{core}+\dot{E}_{shell}=\dot{E}_{tot}$) is guaranteed by the identity (E.13).

4.4 Dissipation in the core

In this section, I will first show that the thin shell approximation causes a very small error on core dissipation, before studying the conditions required for high core dissipation and examining the core dissipation pattern.

First, what is the error on core dissipation due to the thin shell approach? The formula for core dissipation (Eq. (50)) is actually valid beyond the thin shell approximation if the shell density is homogeneous and there is no density contrast at the shell-ocean boundary (the ocean density can increase with depth). In that case, the full solution in the core can be obtained from the fluid-crust solution in the core by gravity scaling (see Appendix F of Beuthe [2015a]):

[TABLE]

As dissipation depends on the product of stress and strain, this procedure accounts for the factor $|(k_{n}+1)/(k_{n}^{\circ}+1)|^{2}$ in Eq. (50). For thin shells, gravity scaling is equivalent to the effective tidal potential trick, as is seen by substituting $U_{n}^{\circ}$ from Table 3 into the core solution given by Eq. (33):

[TABLE]

Using the property that Eq. (50) is applicable to both thin and thick shells (with the restrictions on the shell density mentioned above), I can estimate the thin shell error on $\dot{E}_{core}$ from the thin shell error on $|k_{2}+1|^{2}$ which is less than 0.1% ($k_{2}\sim 10^{-2}$ and the thin shell error on $k_{2}$ is a few percent). Core dissipation thus mainly depends on internal structure through the factor ${\rm Im}(k_{n}^{\circ})/|k_{n}^{\circ}+1|^{2}$.

Tidal heating reaches several tens of GW if the unconsolidated core is modelled as a very soft viscoelastic material [Choblet et al., 2017] (Roberts [2015] studied before the enhancement of tidal heating in a fluffy core, but without global ocean and with the shear modulus of the core larger than the shear modulus of ice). The complex shear modulus of the homogeneous core is parameterized in terms of the elastic shear modulus $\mu_{c\rm e}$ and a nondimensional parameter $\delta$ (zero if the core is elastic, otherwise positive):

[TABLE]

For Maxwell rheology, $\delta$ is related to the core viscosity $\eta_{c}$ by $\delta=\mu_{c\rm e}/(\omega\eta_{c})$, but the above expression is generally valid for any linear rheological model. Choblet et al. [2017] parameterize core rheology with the effective shear modulus $\mu_{\rm eff}=|\mu_{c}|$ and the dissipation function $Q_{\mu}^{-1}=\delta/\sqrt{1+\delta^{2}}$. The latter ranges from 0.2 to 0.8, corresponding to $\delta$ values between 0.5 and 2.

The total power dissipated in the core is maximum if ${\rm Im}(k_{n}^{\circ})/|k_{n}^{\circ}+1|^{2}$ is maximum (see above). At constant $\delta$, this occurs for a homogeneous core if

[TABLE]

where $A_{2}\cong 0.24$ (see Appendix D). Fig. 7A shows the total power dissipated in the core as a function of $|\mu_{c}|$ and $\delta$. The core is homogeneous and incompressible, the ocean is homogeneous and the shell is conductive. The assumption of no differential rotation (Section 2.2) results in overestimating core dissipation by about 25%, but it does not matter since the core rheology is unknown. In Section 5, I set $\delta=1$ and adjust $|\mu_{c}|$ so that the core power is equal to the difference between the conductive power and the shell power.

Although the spatial pattern of dissipation could be obtained from the core dissipation rate under a non-uniform shell (Eq. (41)), it is simpler to use the radial-angular factorization method. After substituting Eq. (52) and Eq. (D.5) into the formulas of Table 4, I can write the radial weights within the core as

[TABLE]

where $h_{2}^{\circ c}$ is the fluid-crust radial Love number at the core-ocean boundary (Eq. (D.4)). The reduced radial weights $\bar{f}_{J}$ are functions of the reduced radius $\hat{r}=r/R_{c}$:

[TABLE]

Apart from a global scaling factor, the radial weights are identical to those for a homogeneous body of radius $R_{C}$ (Eq. (55) of Beuthe [2013]). Therefore, the dissipation pattern in a homogeneous core is identical to the one within a homogeneous body. If heat is transported radially, the flux patterns $(A,B,C)$ at the core surface are weighted by

[TABLE]

Patterns A, B, and C thus contribute respectively 13, 31, and 56 % of the average flux.

The core dissipation flux (Fig. 7B) has nearly no degree-4 harmonic component. The inclusion of the forced libration does not change much the pattern: the flux enhancement due to libration varies between 23% (along the leading-trailing axis) and 30% (at the poles and along the tidal axis). Similarly to shell dissipation, core dissipation is higher at the poles than along the tidal axis. The maximum dissipation contrast, however, is only about a factor of two, whereas it can be larger by an order of magnitude for a conductive shell. Thus, core dissipation hardly explains the spatial variations of the observed surface flux.

5 Thermal equilibrium in a conductive shell

In this section, I study the conditions under which Enceladus’s non-uniform shell is in thermal equilibrium between tidal heating and conductive cooling. I assume here that thermal equilibrium implies a shell in a steady state, but this is not necessarily true because shell thickness variations are progressively destroyed by viscous relaxation at the shell-ocean boundary. This mechanism must be dynamically compensated by ocean freezing or ice melting at the shell-ocean boundary [Čadek et al., 2019].

5.1 Coupling dissipation to heat transfer

Ojakangas and Stevenson [1989] were the first to compute lateral variations of shell thickness by balancing heat production (due to tidal dissipation within the shell and heat flow from the core) against conductive heat transfer. Their model, however, does not take into account the lateral variations of shell thickness and rheology when computing tidal dissipation. Conversely, one should in principle include the effect of tidal dissipation when computing the local temperature profile which determines the rheology of ice. Therefore, tidal dissipation and heat transfer should be solved as a coupled system.

The conductive equilibrium solution is found by iteration. The viscoelastic shell parameters $(\alpha,D,\chi)$ are initially evaluated for the no-dissipation temperature profile (Eq. (64)). One iteration consists of the following three steps:

solving the tidal thin shell equations for the stress function $F$ and the radial displacement $w$ (as in Paper I). 2. 2.

computing the shell dissipation rate $P_{shell}$ (Eqs. (18)-(19)) and the shell dissipation flux ${\cal F}_{shell}$ (Eqs. (21)-(22)). 3. 3.

solving numerically the heat equation with as a source term. The viscoelastic shell parameters corresponding to this new temperature profile are simultaneously evaluated.

The procedure is reiterated until the value of the dissipated power stabilizes, which normally happens after a few iterations. The relative increase in shell power between the no-dissipation solution and the iterated solution will be called rheology feedback. More iterations are needed if the rheology feedback is large.

5.2 Shell structure and core rheology

On the basis of gravity, topography, and libration data, Enceladus is thought to be made of a large silicate core, surrounded by a deep ocean and a thin icy shell. Gravity-topography data combined with the hypothesis of isostasy result in the following model of shell thickness variations (see discussion in Section 5.1 of Paper I):

[TABLE]

where $(d_{00},d_{20},d_{22},d_{30})=(22.8,-12.1,1.3,3.7)\rm\,km$ with $1\sigma$ errors of $(4,2.4,0.3,0.7)\rm\,km$ (uncertainties are ignored below). $P_{nm}$ are the unnormalized associated Legendre functions. Contrary to Paper I, our models include non-zonal variations of shell thickness. The resulting shell thickness is $14.4$ and $7\rm\,km$ at the north and south poles, respectively, and varies between $24.95$ and $32.75\rm\,km$ along the equator (the thickest shell is along the tidal axis). Such models are denoted ‘ISO’ (solid curves in Fig. 8A) whereas models with uniform thickness (equal to $22.8\rm\,km$) are denoted ‘UNI’ (dashed curve in Fig. 8A). When studying asymmetric core dissipation, I also consider the model ‘THIN’ in which the shell is very thin at the south pole. It is parameterized by Eq. (59) in which $(d_{00},d_{20},d_{22},d_{30})=(20.8,-12.1,0,5.7)\rm\,km$ (dotted curve in Fig. 8A). The resulting shell thickness is $14.4$, $26.85$, and $3\rm\,km$ at the north pole, equator, and south pole, respectively.

The icy shell responds to deformations as a linear viscoelastic material with Maxwell rheology. Elastic and viscoelastic parameters are given in Table 1. The viscosity depends on temperature through an Arrhenius relation (Table 2). In this paper, I assume that the bottom viscosity is $10^{13}\rm\,Pa.s$ in order to maximize shell dissipation (see Section 4.3). This value is at the lower end of the range usually considered for melting ice [Tobie et al., 2003; Barr and Showman, 2009] and could lead to fast viscous ice flow destroying the topography at the bottom of the shell. Balancing viscous flow against ice-water phase change, Čadek et al. [2019] argue that the bottom viscosity should be larger than $3\times 10^{14}\rm\,Pa.s$.

Newly published laboratory studies of the anelastic response of ice at tidal frequencies suggest that dissipation could be an order of magnitude higher than predicted by Andrade rheology. [McCarthy and Cooper, 2016]. This phenomenon is simulated here by increasing the nominal dissipation rate by a factor of 10. These models are denoted by the letter ‘H’ (for ‘High’) whereas the models with nominal dissipation are denoted by the letter ‘L’ (for ‘Low’).

In the model ‘THIN’, the shell is further weakened at the south pole (because of faulting) by multiplying the elastic shear modulus by

[TABLE]

where $a=0.1$, $b=5$, and $\theta_{0}=130^{\circ}$ (see Fig. 8B).

If the core is non-porous and elastic, the shear modulus of the core is set to $\mu_{\rm ce}=40\rm\,GPa$. If the core is porous and viscoelastic, core rheology is parameterized by Eq. (54) with $\delta=1$, and the elastic shear modulus is adjusted so that the core power and shell power sum to the conductive power. The shear modulus must be about 1000 times smaller than its elastic value for a non-porous silicate core: $\mu_{\rm ce}=34$, $68.5$, and $32.4\rm\,MPa$ in models ISO-LC, ISO-HC, and THIN-LC, respectively (or $|\mu_{\rm c}|=24.0$, $48.4$, and $22.9\rm\,MPa$, respectively).

5.3 Conductive model

Given that the shell thickness is smaller than $40\rm\,km$, the shell is most likely in a conductive state [Barr and McKinnon, 2007; Mitri and Showman, 2008]. If the shell is radially in thermal equilibrium and the lateral heat transfer is negligible, the temperature $T(r,\theta,\varphi)$ satisfies the radial heat equation:

[TABLE]

The conductive flux denotes the conductive heat flux at the surface of the shell:

[TABLE]

The conductivity of ice decreases as $1/T$ in the range $(40,175)\rm\,K$ but falls below this line at high temperature [Slack, 1980; Petrenko and Whitworth, 1999]. In planetology, the conductivity of ice is often approximated by $k_{ice}=a/T$ with either $a=567{\rm\,W/m}$ [Klinger, 1980] or $a=651{\rm\,W/m}$ [Petrenko and Whitworth, 1999]. The former relation underestimates the conductivity over the whole range but is only 3% too low at the melting temperature, whereas the latter fits well the data up to $200\rm\,K$ but is 11% too high at the melting temperature (see Fig. 9A). When solving numerically the heat equation, I adopt the more accurate fit of Andersson and Inaba [2005],

[TABLE]

with an estimated error of 5% in the range $(40,273\rm\,K)$. Note that the data reviewed by Slack [1980] pertain to single crystals of pure ice: the conductivity is certainly modified by polycrystalline anisotropy and salt contamination. Another factor neglected here is the insulating effect of a $100\rm\,m$-thick snow cover which, if present, lowers the near-surface conductivity by 1 or 2 orders of magnitudes (Fig. 4 of Travis and Schubert [2015]).

Boundary conditions are given by fixing the temperature at the top and bottom of the shell. The surface temperature $T_{\rm s}$ is mainly determined by radiative equilibrium with the annual solar insolation, with a small contribution due to internal heating (see Section 5.4). The temperature at the bottom of the shell (of radius $R_{o}$) is equal to the melting temperature $T_{\rm m}=273\rm\,K$, because the temperature of the ocean is expected to be only slightly less than the freezing point of pure water [Glein et al., 2018].

The heat equation in presence of a source (Eq. (61)) is solved with the BVP solver for boundary value ordinary differential equations [Shampine et al., 2006]. The starting guess is the analytical solution obtained by neglecting dissipation and approximating the conductivity with $k_{ice}=a/T$:

[TABLE]

This profile predicts slightly lower temperatures than the Cartesian profile (identical to Eq. (64) except for the factors $R_{o}/R$ and $R/r$ in the exponents; see Eq. (56) of Paper I). The conductive flux associated with Eq. (64) reads

[TABLE]

where the factor $R_{o}/R$ is the correction due to spherical geometry.

Fig. 9B shows the temperature profile in the shell if there is no dissipation. The top curve (dotted straight line) is associated with constant conductivity in Cartesian geometry. The two intermediate curves result from a conductivity inversely proportional to temperature (Eq. (64), either in Cartesian (solid gray) or in spherical geometry (dashed black). The lowest curve (solid black) is the solution of the heat equation in spherical geometry for the best-fitting conductivity (Eq. (63)). The three curves obtained with a variable conductivity do not differ much between themselves, but decrease much more steeply than the linear profile at the bottom of the shell. Thus, models with constant conductivity overestimate the thickness of the most dissipative layer and the resulting dissipation.

5.4 Surface temperature

The surface temperature depends on solar insolation, albedo, and internal heating sources, all of which vary with latitude and possibly with longitude too. Enceladus’s albedo $A$ is high and varies between $0.74$ and $0.81$ north of $60^{\circ}\rm\,S$ [Howett et al., 2010], while $A=0.80$ fits well the data close to the South Pole [Howett et al., 2011]. Without much error, we can assume a uniform albedo of $A=0.81$ as in Spencer et al. [2006].

Solar insolation is globally proportional to the solar irradiance at Saturn (${\cal F}_{\rm sat}=14.8\rm\,Wm^{-2}$) and varies locally with latitude. Roberts and Nimmo [2008] took the latter factor into account with the approximate formula of Ojakangas and Stevenson [1989]. Neglecting internal heating and assuming unit emissivity, they predicted surface temperatures between 61 and $80\rm\,K$. This range, however, results from setting the average equatorial temperature to $80\rm\,K$, above the subsolar temperature of $76\rm\,K$ [Spencer et al., 2006]. Such values are higher than the near-surface temperature used in models of tidal dissipation, i.e. the approximately constant temperature below the penetration depth of diurnal and seasonal temperature oscillations (respectively about $1\rm\,cm$ and $1\rm\,m$, see Howett et al. [2010, 2011]). Using the formula of Ojakangas and Stevenson [1989] with an albedo $A=0.81$ yields instead an equatorial temperature of $63\rm\,K$ (Fig. 10). A second problem is that Ojakangas and Stevenson [1989] rightly assume that Europa’s obliquity with respect to Jupiter’s orbital plane is small ($i\sim 3^{\circ}$), whereas Enceladus’s obliquity with respect to Saturn’s orbital plane is large ($i\sim 27^{\circ}$). In the latter case, the insolation formula of Nadeau and McGehee [2017] gives a better fit of the mean annual insolation with the additional advantage of being continuous:

[TABLE]

where ${\cal F}_{\rm sat}/4$ is the global annual average insolation. The distribution function $s(\cos\theta,\cos i)$ is approximated by a 6th-order expansion in Legendre polynomials (with unit average on the sphere). As an aside, note that the approximation of Nadeau and McGehee [2017] overestimates the polar insolation if $i<12^{\circ}$; the formula of Ojakangas and Stevenson [1989] actually gives a better fit near the poles if $i<6^{\circ}$ (A. Nadeau, private comm.).

At the poles, the conductive flux becomes comparable to the radiative flux. For example, Eq. (65) yields ${\cal F}_{\rm cond}=150\rm\,mW/m^{2}$ if $d=7\rm\,km$ at the south pole, i.e. 40% of the radiative flux in equilibrium with solar insolation ($380\rm\,mW/m^{2}$ at the same location). For a black body, the equilibrium between solar insolation, reemitted radiation, and internal heating reads

[TABLE]

where $\sigma=5.67\times 10^{-8}\rm\,Wm^{-2}K^{-4}$ is the Stefan-Boltzmann constant. If internal heating is ignored, the mean annual near-surface temperature varies between $51\rm\,K$ (at the poles) and $62\rm\,K$ (at the equator), as shown in Fig. 10. If the crust is nowhere much thinner than $7\rm\,km$, the conductive flux can be approximated by Eq. (65) in which $T_{\rm s}$ is the surface temperature before the correction and $a=632\rm\,W/m$. With this correction, the surface temperature in the Model ISO increases by 2 and 4 degrees at the north and south poles, respectively, while the increase is only about half a degree along the equator, with a very small longitudinal variation. In a model with a very thin crust at the south pole (e.g. $d=2\rm\,km$), the correction is larger and Eq. (67) should be solved self-consistently.

5.5 Results and discussion

5.5.1 Shell dissipation

In good approximation, the conductive flux is inversely proportional to the shell thickness, with marginal influences of surface temperature and in-shell dissipation (Eq. (65)). If the shell thickness is isostatic, the conductive flux is highest at the south pole where it reaches $150\rm\,mW/m^{2}$ (Fig. 11A), with an average over the SPT (below $55^{\circ}\rm\,S$ latitude) of $85\rm\,mW/m^{2}$. The conductive power emitted by the SPT is then $6\rm\,GW$, which is comparable to the SPT power deduced from Cassini infrared data (between 4 and $20\rm\,GW$, see review by Spencer et al. [2018]). It is thus reasonable to treat the conductive power of the isostatic shell as an observational constraint, to be matched by the power generated within Enceladus.

For nominal (i.e. low) dissipation, the shell dissipation flux is only a few percent of the conductive flux (Fig. 11B), as found by Souček et al. [2019]. While the non-uniformity of the shell only slightly increases the shell power (by 20%, see Table 6), it has a major effect on the pattern of the shell dissipation flux. Dissipation within a floating thin shell is typically highest at the poles and lowest along the tidal axis. This is already true for a laterally uniform shell (Fig. 4) but these contrasts are further enhanced by isostatic variations of shell thickness (Fig. 11B and Fig. 12). In that case, the flux is amplified by nearly a factor of 3 at the south pole, where the shell is thinner by a factor of 3 (with respect to the average thickness).

If shell dissipation is ten times higher, the shell power increases more than ten times because of the significant rheology feedback (35% instead of 2% on average, see Table 6). For the same reason, a non-zero forced libration has a larger effect on the shell power if shell dissipation is high (42% increase instead of 29% in the nominal case). In the end, the shell power of the model ISO-H makes up 43% of the conductive power. Furthermore, rheology feedback is much stronger at the poles than at the equator; in particular, it significantly enhances the flux at the south pole (Fig. 12D). Thus, high dissipation within the shell can generate the major part of the conductive flux pattern (Fig. 12C). Nevertheless, it cannot explain the whole conductive flux, whatever the dissipation enhancement, because shell dissipation within the shell remains close to zero along the tidal axis. Thermal equilibrium requires a non-zero heat flux from below the shell, either from core or ocean dissipation.

5.5.2 Partition and scaling

Before looking at models with core dissipation, it is instructive to examine the partition of shell dissipation into membrane/mixed/bending contributions. The shell power is clearly dominated by the membrane contribution (between 85 and 90%, see Table 6). This power partition is slightly misleading because the mixed term switches sign between the poles and the equator (Fig. 13). At the equator, the mixed contribution largely cancels the (always positive) membrane and bending contributions, whereas it makes up a significant portion of the south polar flux (23% and 12% in models UNI-L and ISO-L, respectively). It is thus important to include mixed and bending contributions, and also to keep both of them: neglecting the latter may result in a negative surface flux at the equator.

The membrane term is mainly responsible for amplifying the flux in areas with thinner shell (Fig. 13). The decomposition of the shell dissipation flux (Eq. (22)) explain this behaviour: membrane, mixed and bending terms term are proportional to ${\rm Im}(\alpha)\sim 1/d$, ${\rm Im}(\chi)\sim{}d$, and ${\rm Im}(D)\sim{}d^{3}$, respectively. The factors of Eq. (22) that depend on $F$ and $w$ are less sensitive, in the membrane approximation, to the local shell thickness $d$. On the one hand, the stress function – in the membrane limit of a hard shell – does not depend on the shell properties (see Eq. (80) of Paper I; a shell is hard if $(\mu_{e}/\rho{}gR)(d/R)\gg 1$, see Section 4.3.2. of Paper I). On the other, the shell structure weakly affects the radial displacement of the shell (Fig. 8 of Paper I). Thus, the shell dissipation flux in a laterally non-uniform shell can be approximated by the following (full) scaling rule:

[TABLE]

where superscripts ‘uni’ and ‘iso’ refer to the laterally uniform (UNI) and non-uniform (ISO) models. The terms in the RHS can be computed with analytical formulas (Eq. (E.5)). The amplification factor is defined by

[TABLE]

Thus, if the rheology is approximately laterally uniform, the dissipation flux is inversely proportional to the shell thickness. Beware that the scaling rule (69) only holds if lateral variations of shell properties are of long wavelength and if the shell is everywhere hard, as in models ISO-L and ISO-H. It breaks down, for example, in the model THIN-LC in which bending effects are important and the shell is locally soft at the south pole (Fig. 14).

In Paper I, I showed that the surface stress (and strain) of Enceladus scales with the factor $|\alpha|\sim 1/d$ (see Eq. (88) of Paper I). Why does the shell dissipation flux scale in the same way? By definition, it is proportional to the time average of stress times strain rate, integrated over the shell thickness. In a conductive shell, dissipation mostly occurs in a thin layer at the bottom of the shell. If dissipation is too low to perturb the temperature profile, the thickness of the thin dissipative layer is proportional to the local shell thickness. The amplitudes of the stress and strain rate at the bottom of the shell both scale with $\alpha$. Therefore, the shell dissipation flux scales as $\alpha^{2}d\sim 1/d$ for a hard shell in the membrane limit, which gives back Eq. (69).

5.5.3 Core dissipation

Since shell dissipation cannot account for the full conductive flux, we now turn to models with high core dissipation. Two pieces of (indirect) evidence for a highly porous and dissipative core are, first, its low density (between 2300 and $2600\rm\,kg/m^{3}$, see McKinnon [2015]; Beuthe et al. [2016]) and, second, the hydrothermal activity inferred from silica nanoparticles and molecular hydrogen observed in the plume [Hsu et al., 2015; Sekine et al., 2015; Waite et al., 2017; Glein et al., 2018]. As already seen with the laterally uniform model (Fig. 7), it is possible to tune the viscoelastic core parameters so that core dissipation provides the missing part of the conductive power; it can be similarly done in the laterally non-uniform model by successive adjustments. Unfortunately, laboratory experiments on porous materials have only been done at conditions very different (higher frequency, lower pressure) of those to which the core is submitted during tidal loading. Thus, we don’t know whether the effective shear modulus of the core $|\mu_{c}|$ can really be 1000 times smaller than its elastic value for non-porous silicates.

Although the core can potentially provide whatever is needed for the total heat budget, it faces a more difficult job in explaining the north-south asymmetry of the observed heat flux. The problem is probably made worse by the ocean circulation which partially averages the core heat flux, but I will ignore this complication here. In models with isostatic thickness variation, the core dissipation flux varies only by a factor of two between the poles and the equator (Fig. 15 and Fig. 16B). In model ISO-LC, shell dissipation provides only 5% of the south polar flux while the core contribution reaches 40%, leaving more than half of the conductive flux unaccounted for. If shell dissipation is high (model ISO-HC), shell and core dissipation provide respectively 76% and 22% of the south polar flux, the total matching well the conductive flux. Alhough there is a mismatch of about 25% at the north pole and at the equator, one can probably find a model of lateral variation in thickness and rheology for which the fit is perfect.

Choblet et al. [2017] floated the idea that the shell structure induces a north-south asymmetry in heat production within the core. I will test their proposal with the model THIN-LC, which has a shell structure similar to the one proposed in the Supplementary Information of their paper: the shell is very thin ($3\rm\,km$ thick) and very soft ($\mu_{\rm e}=0.35\rm\,GPa$) in the SPT (Fig. 8). An immediate problem with this model is that the conductive flux more than doubles at the south pole and cannot be matched by internal dissipation. Leaving this problem aside, we observe that core dissipation differs by only a few percent between the north and south poles (Fig. 16), confirming thus the preliminary analysis made after Eq. (34). Therefore the very asymmetric shell structure does not translate into a large asymmetry in core dissipation.

Although hydrothermal flow within the core focuses the tidally dissipated heat towards the poles [Travis and Schubert, 2015; Choblet et al., 2017], it cannot by itself produce a strong north-south asymmetry unless the core is very inhomogeneous. An asymmetrical core structure could be due to its conditions of formation. This impact of this assumption on tidal dissipation and hydrothermal flow remains to be investigated.

6 Summary

Enceladus’s high heat output is attributed to dissipative eccentricity tides, but more by default than by a true understanding of its modus operandi. So far, geodynamical models have been unable to account for the anomalous heating at the south pole both in magnitude and in localization. The north-south asymmetry in geophysical activity is particularly telling, because it negates the usual assumption of a spherically symmetric internal structure which entails equal dissipation at the north and south poles. The only other information we have about an asymmetric interior resides in gravity and topography data: in an isostatic framework, they imply that the shell is thinnest at the south pole, medium thick at the north pole, and thickest at the equator. Thus, the north-south asymmetries observed at the surface and inferred in the interior make it pressing to compute tidal dissipation in a body with a floating shell of variable thickness. It is also likely that the shell is softer in the SPT, so that lateral variations of shell rheology should also be included.

In Paper I, I developed the theory of non-uniform viscoelastic thin shells, coupled it to tides, and solved the resulting tidal thin shell equations in terms of the stress function $F$ and the radial displacement $w$. The ultimate purpose was the prediction of tidal stresses within the shell. Here, I show how to use the same $(F,w)$ solution to predict tidal dissipation in the core and shell. The former must be internally spherically symmetric, but the latter can be fully non-uniform in thickness and rheology. The shell must however satisfy the following requirements (approximately satisfied in Enceladus): thickness less than 10 to 20% of the surface radius, homogeneous shell density, negligible density contrast between the shell and the top layer of the ocean, uniform Poisson’s ratio, and linear viscoelasticity. The uniformity of Poisson’s ratio deserves special attention as it implies a new dissipation constraint, called ‘no Poisson dissipation’ in contrast with the more usual condition of ‘no bulk dissipation’ related to the uniformity of the bulk modulus. Dissipation predictions under the two constraints differ by a few percent, which is generally smaller than the error due to the common incompressible assumption, except for the thinnest shells (Fig. 1). While the uniformity of Poisson’s ratio can be seen as a drawback of the thin shell approach, it is not certain either that bulk dissipation vanishes in tidal heating processes.

The thin shell approach is basically a 2D theory because its primary variables ($F$ and $w$) do not depend on depth. Nevertheless, thin shell theory does not assume that all variables are constant with depth. In particular, the strain varies quasi linearly across the shell thickness, while the stress varies in proportion to the local strain and local shear modulus. Thus, the thin shell approach predicts 3D quantities such as the volumetric dissipation rate (Eqs. (18)-(19)). Integrating the dissipation rate then yields the shell dissipation flux at the surface and the total shell power (Eqs. (21)-(23)). In that ‘micro’ framework, the shell dissipation rate, flux, and power are expressed as sums of membrane ($FF$), mixed ($Fw$), and bending ($ww$) contributions, the first one being by far the largest.

Regarding core dissipation, I compute core deformations by treating the non-uniform shell as a pressure load (effective tidal potential method), before computing the core dissipation rate, flux, and power, which are nearly insensitive to lateral variations of shell properties (Eqs. (41)-(44)). If one is only interested by the total power produced in the core, it is easier to partition the total power produced in the body into core and shell contributions (Eqs. (39)-(40)): this is the ‘macro’ approach to tidal dissipation, to be contrasted with the previous ‘micro’ approach. Comparing the core power and shell power obtained in the micro and macro approaches provides a self-consistency check on the correctness of numerical codes.

Benchmarking against a laterally uniform shell provides a good way to estimate errors due to the thin shell approximation. The impact of a non-zero shell thickness on the shell power is less than 2 and 4% (for a laterally uniform shell) if the shell thickness is less than 20 and $50\rm\,km$, respectively. These errors are computed by comparing apples with apples, i.e. thin and thick shell models having the same uniform Poisson’s ratio (Fig. 6A). If the thin shell is compared to the thick shell with no bulk dissipation, the difference is below 4% if the shell thickness is between 20 and $50\rm\,km$, but can climb to 9% for very thin shells (Fig. 6B). Spatial patterns of the shell dissipation rate and flux are well predicted by the thin shell approach. The benchmarking of a laterally non-uniform thin shell against the finite-element method of Souček et al. [2019] is outside the scope of this paper, but preliminary comparisons with FEM results show good agreement [Běhounková et al., 2018].

As an illustration, I compute tidal dissipation within Enceladus’s shell assuming that it is conductive, in thermal equilibrium, and with thickness variations predicted by isostasy. The conductive model is as realistic as possible regarding ice conductivity (dependent on temperature) and surface temperature (latitudinal variation). Dissipation and heat transfer are solved as a coupled system, because viscosity depends strongly on temperature which in turn depends on internal heating, although rather weakly. If dissipation is low, the temperature profile is mostly determined by the surface and melting temperatures, but high dissipation results in a large rheology feedback (about 35%): the rheology of the bottom of the shell becomes softer because of locally higher temperatures. The addition of the forced libration to eccentricity tides increases dissipation by about 30%. Variations in shell thickness only slightly increase the total shell power (here by 20%), but they have a major effect on the shell dissipation pattern: dissipation is highest where Enceladus’s shell is thinnest. In particular, the shell dissipation flux at Enceladus’s south pole is three times larger in the isostatic model ($d_{SP}=7\rm\,km$) than in the uniform thickness model ($d=23\rm\,km$). If the shell is hard and the lateral variations of the shell properties are of long wavelength, one can actually predict the shell dissipation flux by scaling the flux for a laterally uniform shell (Eqs. (68)-(69)). The scaling factor depends on the inverse shell thickness (at least if rheology is nearly laterally uniform) and can be understood in the thin shell approach by the dominance of the membrane contribution.

The shell power for nominal values of viscoelastic parameters is only a few percent of the conductive power, in agreement with Souček et al. [2019]. If the shell dissipation rate is ten times higher, as suggested by recent laboratory experiments, it contributes to nearly half the conductive power and accounts for most of the spatial variations of the conductive flux. The same effect is obtained with an eccentricity three times larger than the present one. It does not make sense to explore higher rates of dissipation within the shell, because the shell dissipation flux remains too small along the tidal axis to balance conductive cooling. In a steady state model, dissipation within an unconsolidated core must contribute the rest. Given that core or ocean dissipation cannot be avoided, it is logical to ask whether one could entirely dispense with shell dissipation. Core dissipation, however, remains nearly north-south symmetric even in a model with an extremely thin and soft shell at the south pole. Therefore, dissipation in a homogeneous core cannot be responsible for the flux asymmetry observed at the surface.

In conclusion, explaining Enceladus’s heat anomaly in magnitude and localization requires pushing the envelope very far: dissipation must be simultaneously high in the shell and in the core. To avoid this non-parsimonious solution, one should investigate whether dissipation within a inhomogeneous core is sufficient to maintain the non-uniform shell in thermal equilibrium. Alternatively, non-steady state models – which maybe solve the problem of the total power without recourse to core dissipation – must demonstrate that they can predict the surface heat pattern and the variations of shell thickness. Models in which Enceladus’s shell is currently thinning face the same difficulties as thermal equilibrium models in explaining the origin of the emitted energy and the non-uniformity of tidal heating and shell thickness. If Enceladus’s shell is currently thickening [Luan and Goldreich, 2017], the higher dissipation flux in the past could have provided enough heat to keep the ocean liquid, but present-day shell thickness variations are more difficult to explain: thickness variations generated in the past are not only smaller (because they are bounded by the smaller average thickness) but they also tend to be averaged out by faster thickening where the shell is thinner.

Acknowledgments

I thank Attilio Rivoldini for his help with the BVP solver, Antony Trinh for information about the libration of the core, and Alice Nadeau for explanations about equilibrium surface temperatures. This work is financially supported by the Belgian Federal Science Policy Office through the Brain Pioneer contract BR/314/PI/LOTIDE and the PRODEX grant No. 4000120791.

Appendix A Spherical differential operators

Table 7 gives the definitions of the spherical differential operators appearing in this paper. Table 8 gives their expressions in terms of scalar operators.

Appendix B Dissipation rate in terms of $F$ and $w$

The shell dissipation rate (Eq. (15)) is a linear combination of the strain invariants ${\cal E}_{2}$ and ${\cal E}_{\rm tr}$ (Eqs. (16)-(17)). The tangential strains are given by Eq. (J.1) of Paper I, here corrected for a typo:

[TABLE]

where $z$ and $\chi$ are defined in Table 2. The potentials for the tangential load $\Omega$ and its moment $\Omega_{M}$ (Eq. (E.5) of Paper I) vanish here ($\Omega=\Omega_{M}=0$) because tangential loads are neglected in the tidal coupling of the thin shell. In Paper I, the factor $(1-\nu)$ was missing in front of the terms $(\Omega+\Omega_{M})$, but other equations of Paper I are not affected by this typo. Another typo without consequences is that the definitions of the elastic extensibility and bending rigidity are exchanged in Eq. (4) of Paper I.

Computing the trace invariant ${\cal E}_{\rm tr}$ is straightforward. After substituting Eq. (B.1) into Eq. (17), I apply the identity ${\cal O}_{1}+{\cal O}_{2}=\Delta^{\prime}$ (Table 7). I write the result as

[TABLE]

where

[TABLE]

and ${\cal E}_{\rm tr}^{F^{*}w}=({\cal E}_{\rm tr}^{Fw^{*}})^{*}$.

Computing the strain invariant ${\cal E}_{2}$ is a bit more involved. After substituting Eq. (B.1) into Eq. (16), I combine the operators $({\cal O}_{1},{\cal O}_{2},{\cal O}_{3})$ into $\Delta^{\prime}$, ${\cal A}$, and ${\cal D}_{4}^{\prime}$ using the definitions of Table 7. Next, I express ${\cal D}_{4}^{\prime}$ in terms of $\Delta^{\prime}$ and ${\cal A}$ with the identity (d) of Table 8. I write the result as

[TABLE]

where

[TABLE]

and ${\cal E}_{2}^{F^{*}w}=({\cal E}_{2}^{Fw^{*}})^{*}$. Substituting Eqs. (B.2) to (B.5) into Eq. (15), I can write the dissipation rate within the shell in the form of Eqs. (18)-(19).

Appendix C Surface flux in terms of $F$ and $w$

The shell surface flux (Eq. (20)) is given by

[TABLE]

where $({\cal E}^{FF},{\cal E}^{Fw},{\cal E}^{ww})$ are given by Eq. (19) and $\bar{\zeta}=\zeta/R$.

In more compact notation, I write

[TABLE]

where

[TABLE]

First, the integral over ${\cal E}^{FF}$ is proportional to

[TABLE]

Now,

[TABLE]

so that

[TABLE]

Second, the integral over ${\cal E}^{Fw}$ includes a term proportional to

[TABLE]

Now,

[TABLE]

so that

[TABLE]

The contribution of the other term (proportional to $\alpha{}z^{*}$) gives the same result.

Third, the integral over ${\cal E}^{ww}$ is proportional to

[TABLE]

Now,

[TABLE]

so that

[TABLE]

Substituting Eqs. (C.6), (C.9), and (C.12) into Eq. (C.1), I can write the shell surface flux as Eqs. (21)-(22).

Appendix D Deformation of the core

In the thin shell approach, the tidal deformation of the core is obtained by solving the viscoelastic-gravitational problem for the associated fluid-crust model forced by an effective tidal potential (Section 3.3.1). If the core is homogeneous and incompressible and the ocean is homogeneous, the propagator matrix method [Sabadini and Vermeersen, 2004] yields analytical formulas for the Love numbers and the functions $y_{in}^{T}$.

Non-dimensional parameters are defined by

[TABLE]

The fluid-crust tidal Love numbers are (Appendix C.2 of Beuthe [2015b])

[TABLE]

where $(A_{n},B_{n})$ are polynomials in $(y,\xi)$ defined in Table 9. If the core is very rigid ($\mu_{c}>0.1\rm\,GPa$), these Love numbers become independent of the core parameters:

[TABLE]

The radial Love number at the core-ocean boundary reads

[TABLE]

The three radial functions required for the dissipation rate read

[TABLE]

where $\hat{r}=r/R_{c}$ and ${\rm f}_{ni}(\hat{r})$ are non-dimensional functions defined by Table 9. This solution tends to the well-known solution for a homogeneous body if one takes the limits $y\rightarrow 1$ and $\xi\rightarrow 0$ (e.g. Eq. (D.3) of Beuthe [2016]). Since Eq. (D.5) depends on $(y,\xi)$ through the common factor $h_{n}^{\circ c}$, the strains in the core of a 3-layer body are scaled down by a common factor from the strains in a homogeneous body.

Appendix E Dissipation in a laterally uniform thin shell

Dissipation rate

If the shell is laterally uniform and the tidal potential is of degree $2$, it is possible to factorize the thin shell dissipation rate $P_{shell}$ (Eqs. (18)-(19)) in radial and angular functions as in Eq. (45). One just needs to know that:

The variables $(F,w)$ are proportional to the tidal potential $U_{2}$ (see Table 3). 2. 2.

The operator $|\Delta^{\prime}U_{2}|^{2}$ depends on the angular function $\Psi_{A}$:

[TABLE] 3. 3.

The operator ${\cal A}(U_{2};U_{2}^{*})$ is a linear combination of the angular functions $(\Psi_{A},\Psi_{C})$:

[TABLE]

where the first equality results from the identity (i) of Table 8 and the second one from Eq. (22) of Beuthe [2013].

The three components of $P_{shell}$ are thus equal to

[TABLE]

where $f_{2}=(R|h_{2}|/g)^{2}\omega^{4}$. Equating $P_{shell}$ to the factorized dissipation rate (Eq. (45)), I obtain the dissipation weight functions for the thin shell:

[TABLE]

where $r=R+\zeta$. The three terms within the brackets correspond to membrane, mixed, and bending contributions, respectively.

Surface flux

Similarly to the dissipation rate, the surface flux of a laterally uniform thin shell can be written as a weighted sum of the angular functions $\Psi_{J}$. Equating ${\cal F}_{shell}$ (Eqs. (21)-(22)) to the factorized surface flux (Eq. (46)), one can show that

[TABLE]

in which $f_{2}=(R|h_{2}|/g)^{2}\omega^{4}$ as above. The surface flux weights (Eq. (47)) thus read

[TABLE]

In the membrane limit, the contribution of Pattern C to the average surface flux is

[TABLE]

The first line can be deduced from Eq. (E.6) while the second line results from Eq. (94) of Beuthe [2015a]. The parameter $\bar{\kappa}$ quantifies the effective bulk dissipation of the thin shell and varies between 0 (incompressible limit) and 0.5 for a conductive shell with $\nu_{e}=0.33$. If the membrane is incompressible, the ${\cal F}_{C}/{\cal F}_{T}$ ratio is equal to $9/11$ (whatever the dissipation condition) and reproduces (if there is no forced libration) the pattern obtained by Ojakangas and Stevenson [1989] (see Eq. (61) of Beuthe [2013]). If the membrane is compressible, the ${\cal F}_{C}/{\cal F}_{T}$ ratio is always equal to $3/4$ ($\nu_{e}=1/3$) if there is no Poisson dissipation but varies between $3/4$ and $9/11$ in there is no bulk dissipation, the precise value depending on the bottom viscosity (see Table 5).

Thin shell power

In the micro approach, the total power dissipated in the shell is obtained by integrating the surface flux (Eq. (23)). If the shell is laterally uniform, the integral can be done analytically by substituting the degree-$n$ solution of Table 3 into the surface flux equations (Eqs. (21)-(22)), replacing the operator $\Delta^{\prime}$ by its eigenvalue $\delta_{n}^{\prime}$ and integrating the operator ${\cal A}$ with the identity (j) of Appendix A. The result reads

[TABLE]

where the membrane, mixed, and bending terms are given by

[TABLE]

in which $c_{n}=(2\pi\omega R^{2}\rho/g)|h_{n}|^{2}\langle|U_{n}^{T}|^{2}\rangle$. For the reference conductive shell ($d=23\rm\,km$, $\eta_{\rm m}=10^{13}\rm\,Pa.s$), the membrane/mixed/bending contributions to the degree-2 total power are

[TABLE]

The membrane and bending terms can be summed with the identity ${\rm Im}(a^{*}b)+2\,{\rm Im}(a){\rm Re}(b)={\rm Im}(ab)$ (valid for any complex numbers $a$ and $b$), where $a=\chi$ and $b=\Lambda_{n}^{M}/\psi$. The total power can thus be written as

[TABLE]

where $\Lambda_{n}^{corr}$ includes the next-to-leading contributions (see Table 3). For the same reference conductive shell as above, the membrane/bending/next-to-leading contributions to the degree-2 total power are

[TABLE]

The mixing contribution is now included in the membrane term, whereas it was separate in Eq. (E.10). Both equations predict the same bending contribution of 5% (differences are of second order in the thin shell approximation; $\Lambda_{2}^{corr}$ is also of second order, see Appendix I of Paper I).

Core-shell partition

The partition of the total power (Eq. (48)) into core and shell contributions is equivalent to the decomposition of ${\rm Im}(k_{n})$ into a term proportional to ${\rm Im}(k_{n}^{\circ})$, associated with the interior below the shell, and a term proportional to ${\rm Im}(\Lambda_{n})$, associated with the shell itself:

[TABLE]

This identity was already proved for a membrane (see Appendix H of Beuthe [2015a]) and is also valid for a thin shell.