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

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

Mikael Beuthe
Royal Observatory of Belgium, Avenue Circulaire 3, 1180 Brussels, Belgium
E-mail: mikael.beuthe@observatoire.be

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. The micro approach is supplemented by a macro approach providing an independent determination of the core-shell partition of the total power. Unless the shell is incompressible, the assumption of a uniform Poisson’s ratio implies non-zero bulk dissipation. 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.

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 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 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.

Parameter Symbol Value Unit
Mean eccentricity -
Rotation rate
Bulk density
Surface gravity
Surface radius km
Reference thickness (if uniform shell) 23 km
Core radius 192 km
Density of ice and ocean 1000
Shear modulus of ice 3.5 GPa
Bulk modulus of ice 9.13 GPa
Poisson’s ratio of ice 0.33 -
Shear modulus of core (if non-porous) 40 GPa
Conductivity of ice (if uniform shell)
Activation energy
Melting ice temperature 273 K
Viscosity of ice at melting point Pa.s
Source given in Table 1 of Beuthe [2016].
Beuthe et al. [2016].
Helgerud et al. [2009].
Petrenko and Whitworth [1999] (Section 4 only); Eq. (63) is used in Section 5.
Table 1: Physical parameters used in this paper.

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:



  • is the radial displacement and is the auxiliary stress function;

  • and are the primary viscoelastic shell parameters (Table 2); is a secondary viscoelastic shell parameter, close to one (Table 2);

  • and 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 with standard matrix methods. Other variables, such as tangential displacements, strains, and stresses, can be computed from 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).

Symbol Name Definition
radial shell coordinate
depth parameter
viscosity of ice
complex shear modulus
Depth-integrated ()
th moment of
invariant second moment
invariant extensibility
invariant bending rigidity
bending rigidity
is the gas constant.
For Maxwell rheology, but other linear rheological models can be used.
Table 2: Viscoelastic shell parameters.

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 :


where is the degree-0 load (see Section 3.4 of Paper I) while the degree-1 load vanishes. If , the degree- load is given by


where is the weight of the tidal bulge in the unperturbed gravity field. The total perturbing potential depends on the primary tidal potential and on the secondary potential induced by the deformation of the whole body. In the thin shell approach, it is given by


where is the degree- density ratio,


The spherically symmetric structure below the shell (stratified ocean, viscoelastic core) enters into the problem via the nondimensional factor (),


where 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, and . In that case, the term in (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 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 ), I write the tidal potential as [Van Hoolst et al., 2013]


where for Enceladus [Thomas et al., 2016]. are the associated Legendre functions of degree and order depending on ( is the colatitude and is the longitude). The amplitudes of the optical and forced librations must be added because the forced libration and the tidal torque are 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: implies that the neglected corrections are of order . 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 . 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 equal to the mean motion. It is thus convenient to work with Fourier-transformed variables: . The dissipation power per unit volume averaged over one orbital period , in short the dissipation rate, reads (Appendix A of Beuthe [2013])


where the asterisk denotes complex conjugation. The tensors and 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 , and I drop the ‘tilde’ notation. According to the correspondence principle, linear viscoelasticity is introduced through complex moduli (shear modulus and bulk modulus ) in the frequency-domain stress-strain relation:


where denotes the trace of the 3D strain tensor and is the Kronecker delta. If Eq. (9) holds, the dissipation rate becomes


The term proportional to represents dissipation resulting from shear friction. The term proportional to 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’: or equivalently . 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 because it would require that Poisson’s ratio varies in tandem with the shear modulus (e.g. Appendix C of Beuthe [2015a]) according to the -- relation:


This would contradict the assumption of uniform 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: . This constraint together with Eq. (11) implies that


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


If there is no bulk dissipation, the factor should be replaced by . In the incompressible limit (, ), bulk dissipation always vanishes so that the conditions of ‘no bulk dissipation’ and ‘no Poisson dissipation’ are simultaneously satisfied. In that case, 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, is maximum at the forcing angular frequency if the viscosity of the material is equal to the critical viscosity , where is the elastic shear modulus. For Enceladus’s icy shell, the critical viscosity is equal to . But a conductive shell is not at all homogeneous: the viscosity, which controls , 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 , 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 , which has the advantage that the dependence on the tidal potential (eccentricity and libration) is factored out (Eq. (48)). If the shell thickness is in the range , 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 , but may reach 9% in the membrane limit if . 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’.

Figure 1: Impact of the dissipation condition on the total power dissipated in a thick shell. The total power is parameterized by the imaginary part of (see Eq. (48)). (A) if there is no bulk dissipation () for three values of the bottom viscosity ; the right-hand scale gives the shell power including the contribution of the forced libration. (B) Relative change of with respect to Panel A if there is no Poisson dissipation (black curves) or if the shell is incompressible (gray curves). Solid/dashed/dotted curves correspond to the bottom viscosities specified in Panel A. The vertical line corresponds to . See Section 3.1.

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)


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


where are 2D strain invariants:


If there is no bulk dissipation, the factor in Eq. (15) should be replaced by the factor ; both factors tend to 1 in the incompressible limit.

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


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


The shell dissipation rate depends on depth through and the depth parameter (Table 2). It does not depend on the choice of reference surface, since are invariant under a change of (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,


where the factor is associated with the integration measure in spherical coordinates. The shell dissipation flux can be expressed in terms of the variables and of the parameters (see Appendix C):


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


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: , , .

  • 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 (invariance under rigid displacements) and of (degree-1 gauge freedom), which belong to the null space of and .

  • It is inversely proportional to the area of the reference surface of the shell (as it should be), because is invariant under a change of while other quantities scale as , , and (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: .

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 :


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]):


where is the secondary potential due to the deformation of the body:


in which 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 :


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

If the core is elastic () and the ocean is inviscid, the total dissipated power is equal to the shell power. Setting into Eq. (26) and using Eqs. (3)-(4), I can write


which can be interpreted as follows: 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 are in phase if is real and (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 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


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  are represented in a solid layer by 6 radial response functions (). Since the tidal and pressure load solutions correspond to different boundary conditions, the full solution reads


where the tidal potential is evaluated at radius . The variables (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 (). 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 (radial displacement), (radial stress) and (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 and [Saito, 1974]. At the outer fluid surface, there is only one (inhomogeneous) boundary condition given by


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 ().

Defining the effective tidal potential for the core by


I can write the full solution in the core (Eq. (29)) as the tidal solution forced by :


The similar property () holding within the fluid leads to relations between pressure and tidal Love numbers: and , which were already noted as Eqs. (29)-(30) of Paper I (the superscript is omitted for tidal Love numbers).

For evaluation purposes, it is practical to write the effective tidal potential in terms of and the flexure solution :


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 (see Section 4.1), so that the second term in the brackets of Eq. (34) is small with respect to the first: (, 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 to (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


which results from Eqs. (3) and (28). Besides, is related to the effective tidal potential for the core by


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

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


The total power dissipated in the body thus reads




Similarly to Eq. (27), 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, vanishes because the tidal thin shell equations impose that are in phase.

If is the shell power in Eq. (40), 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, and can be expressed in terms of 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 ), 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




, , and correspond respectively to , , and in Beuthe [2013]. The differential operators and are defined in Appendix A. The functions are defined by Eqs. (33)-(34): is the radial displacement, is the potential for tangential displacement, and is the potential for radial-tangential stress (or strain). If the core and the ocean are both homogeneous, the functions appearing in are given by Eq. (D.5).

As the shell is laterally non-uniform, the functions 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 (or the related operator ), as can be seen for the operators and 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,


The total power dissipated in the core is obtained by integrating the flux over the core surface with surface element :


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 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 , 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 .

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 , three of which are needed here: (radial displacement), (potential for tangential displacement), and (potential for radial-tangential stress or strain).

Symbol Name Solution Example
Basic variables
radial displacement
stress function
total perturbing potential
tidal load
effective tidal potential
Spring Constants (SC)
membrane SC
bending SC
next-to-leading SC
thin shell SC
Tidal Love numbers (TLN)
fluid-crust radial TLN Eqs. (D.2)-(D.3)
fluid-crust gravitational TLN
radial TLN
gravitational TLN
are defined in Table 2; numerical values are given in Table 3 of Paper I.
is the degree- density ratio, see Eq. (5); , see Appendix A.
If , , ; and the core is elastic (other parameters given in Table 1).
Table 3: Degree- tidal deformation of a laterally uniform thin shell (see Section 4 of Paper I, except for which is given by Eq. (36)).

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


where are the radial weight functions while are angular functions representing the basic spatial patterns. The ratios (with ) 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 , and specifies the angular functions for degree- 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 :




The ratios (with ) measure the contributions of the patterns to the average flux (or to the total power).

Radial weight Angular function
Harmonic basis
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: (Eq. (7)). The symbols are the viscoelastic-gravitational solutions of Takeuchi and Saito [1972]. If the body is incompressible, and . is the weight for the angular average of the dissipation rate. denote the unnormalized associated Legendre functions of degree and order with argument .

Fig. 2 shows the basic dissipation patterns for degree-2 eccentricity tides combined with forced libration ( or ). 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). includes significant contributions from both degrees 2 and 4, has maxima along the equator and minima at middle latitudes. does not contribute in the thin shell approximation (but it contributes to core dissipation). is dominated by the degree-2 component of the squared tidal potential, is maximum at the poles, and zero along the tidal axis.

Figure 2: Basic patterns of tidal heating in a laterally uniform shell, due to degree-2 eccentricity tides and forced libration. The tidal axis goes through longitude. Each pattern repeats from to . The amplitude has been divided by . See Table 4 for analytical expressions.

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: and . In Appendix E, I substitute the analytical solutions for given in Table 3 into the general expression of the thin shell dissipation rate (Eqs. (18)-(19)). The resulting formulas for the radial weights and 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 , 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 , 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.

Figure 3: Dissipation rate in a laterally uniform conductive shell (, ): radial weight functions for (A) an incompressible shell, and (B) a compressible shell. Solid curves show the exact results for a thick shell without bulk dissipation. Dashed curves show the thin shell results without Poisson dissipation (Eq. (E.4)). Big dots indicate the values of the radial weights in the membrane limit. See Section 4.2.

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 (or less than 3% if ). For thinner shells, the contribution of Pattern C decreases, reaching a lower bound of to