Enceladus’s crust as a nonuniform thin shell:
II tidal dissipation
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 nonuniform icy shell thickness can cause the northsouth heating asymmetry by redistributing tidal heating either in the shell or in the core. Starting from the nonuniform 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 coreshell partition of the total power. Unless the shell is incompressible, the assumption of a uniform Poisson’s ratio implies nonzero 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 longwavelength 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. Nonsteady state models, though not investigated here, face similar difficulties in explaining the asymmetries of tidal heating and shell thickness.
Contents
 1 Introduction
 2 Tidal thin shell equations
 3 Dissipation in the thin shell approach
 4 Benchmarking against a laterally uniform thick shell
 5 Thermal equilibrium in a conductive shell
 6 Summary
 A Spherical differential operators
 B Dissipation rate in terms of and
 C Surface flux in terms of and
 D Deformation of the core
 E Dissipation in a laterally uniform thin shell
List of Figures
 1 Impact of dissipation condition
 2 Basic patterns of tidal heating in a laterally uniform shell
 3 Dissipation rate in a laterally uniform conductive shell
 4 Surface flux pattern in a laterally uniform conductive shell
 5 Core dissipation under a laterally uniform shell
 6 Total power dissipated in a laterally uniform shell
 7 Core dissipation under a laterally uniform shell
 8 Shell structure
 9 Conductive model
 10 Surface temperature
 11 Heat flux pattern in model ISOL
 12 Heat flux and rheology feedback if no core dissipation
 13 Partition and scaling of the shell dissipation flux
 14 Partition of the shell dissipation flux into membrane/mixed/bending contributions
 15 Surface heat flux if high dissipation within the core
 16 Surface heat flux if weak shell in SPT
List of Tables
 1 Physical parameters
 2 Viscoelastic shell parameters
 3 Tidal deformation of a laterally uniform thin shell
 4 Radial weights and angular functions factorizing the dissipation rate
 5 Pattern of surface shell flux
 6 Thermal equilibrium models
 7 Spherical differential operators: definitions
 8 Spherical differential operators: identities
 9 Nondimensional functions for core deformations
1 Introduction
Despite being small and cooling fast, Enceladus is warm. The first evidence for its presentday 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 nonuniform 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 SaturnEnceladusDione 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; RoviraNavarro 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 thermalorbital evolution cycle between phases of high and low dissipation. Most models of this sort postulate that Enceladus is now in a shortlived 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 presentday phase of low dissipation in a mediumthick 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 northsouth heating asymmetry if the shell is nonuniform, as suggested by Choblet et al. [2017]? Or else, is the nonsteady 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 nonuniform thickness. Here I will solve this problem with the nonuniform 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 nontrivial 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 gravitytopography 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 nonuniform shell induce a northsouth asymmetry in core dissipation?
In a previous study of the topic, Běhounková et al. [2017] used the finiteelement method (FEM) of Souček et al. [2016] to solve for the deformations of an elastic shell with nonuniform 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 temperaturedependent 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 nonzero 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 nonMaxwell rheological models (e.g. the pseudoAndrade rheology in Souček et al. [2019]); and that it neglects selfgravity (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 nonuniform 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 nonuniform 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 nonuniform 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 nonporous)  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. 
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 nonuniform thin spherical shell floating on a global ocean:
(1) 
where

is the radial displacement and is the auxiliary stress function;

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 
Depthdependent  
radial shell coordinate  
depth parameter  
viscosity of ice  
complex shear modulus  
Depthintegrated ()  
th moment of  
invariant second moment  
  
  
invariant extensibility  
invariant bending rigidity  
extensibility  
bending rigidity  
is the gas constant.  
For Maxwell rheology, but other linear rheological models can be used. 
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 :
(2) 
where is the degree0 load (see Section 3.4 of Paper I) while the degree1 load vanishes. If , the degree load is given by
(3) 
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
(4) 
where is the degree density ratio,
(5) 
The spherically symmetric structure below the shell (stratified ocean, viscoelastic core) enters into the problem via the nondimensional factor (),
(6) 
where is the tidal radial Love number of the fluidcrust 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 (selfattraction or selfgravity). 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 degree2 eccentricity tides plus the forced libration (denoted ), I write the tidal potential as [Van Hoolst et al., 2013]
(7) 
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 nonuniform 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 Fouriertransformed variables: . The dissipation power per unit volume averaged over one orbital period , in short the dissipation rate, reads (Appendix A of Beuthe [2013])
(8) 
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 frequencydomain 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 frequencydomain stressstrain relation:
(9) 
where denotes the trace of the 3D strain tensor and is the Kronecker delta. If Eq. (9) holds, the dissipation rate becomes
(10) 
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:
(11) 
This would contradict the assumption of uniform made in the theory of nonuniform 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
(12) 
If there is no Poisson dissipation, the dissipation rate thus reads
(13) 
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. CastilloRogez 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’.
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)
(14) 
Substituting this constraint into Eq. (13), I write the shell dissipation rate in absence of Poisson dissipation as
(15) 
where are 2D strain invariants:
(16)  
(17) 
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):
(18) 
If there is no Poisson dissipation, the terms in the RHS are given by
(19) 
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,
(20) 
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):
(21) 
If there is no Poisson dissipation, the terms in the RHS are given by
(22) 
The subscripts mem, mix, bend stand for membrane, mixed and bending contributions. This is analogous to the decomposition of the elastic energy in extensionalshearing, mixed, and bendingtwisting 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 depthintegrated 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 degree1 spherical harmonic components of (invariance under rigid displacements) and of (degree1 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 staticgeometric 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 :
(23)  
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 (ZschauPlatzman formula, see Eq. (7) of Platzman [1984]):
(24) 
where is the secondary potential due to the deformation of the body:
(25) 
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 :
(26) 
where the bracket notation denotes the angular average of (or degree0 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
(27) 
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 nontrivial check on the integrated dissipation rate given by Eq. (23).
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) viscoelasticgravitational equations from the center to the surface of the body before evaluating the dissipation rate at each point. Here, the laterally nonuniform shell breaks spherical symmetry. Nevertheless, radial viscoelasticgravitational 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
(28) 
which acts on the associated fluidcrust 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 fluidcrust body (Eq. (26) of Paper I).
The radial viscoelasticgravitational equations should now be solved in a 2layer body (core plus surface ocean) submitted to tidal and pressure loads. Beware that ‘2layer’ does not imply here uniform layers: density and rheology can be depthdependent. 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
(29) 
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 viscoelasticgravitational differential equations. At the coreocean 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 coreocean 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
(30)  
(31) 
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
(32) 
I can write the full solution in the core (Eq. (29)) as the tidal solution forced by :
(33) 
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 :
(34) 
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 Coreshell 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
(35) 
which results from Eqs. (3) and (28). Besides, is related to the effective tidal potential for the core by
(36) 
Using Eqs. (32), (35), and (36), I write the integrand of Eq. (24) as
(37)  
The total power dissipated in the body thus reads
(38) 
where
(39)  
(40) 
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 coreshell 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 selfconsistency 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 nonuniformity 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
(41) 
where
(42) 
, , 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 radialtangential 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 nonuniform, 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 coreocean boundary is equal to the dissipation rate integrated over the core radius,
(43) 
The total power dissipated in the core is obtained by integrating the flux over the core surface with surface element :
(44) 
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 radialangular 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 3layer 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 elasticgravitational 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 radialtangential 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  
nexttoleading SC  
thin shell SC  
Tidal Love numbers (TLN)  
fluidcrust radial TLN  Eqs. (D.2)(D.3)  
fluidcrust 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). 
If the internal structure is spherically symmetric, the dissipation rate can be factorized [Beuthe, 2013]:
(45) 
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 viscoelasticgravitational 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 :
(46) 
where
(47) 
The ratios (with ) measure the contributions of the patterns to the average flux (or to the total power).
Radial weight  Angular function  
Harmonic basis  
0  
2  
4 
Fig. 2 shows the basic dissipation patterns for degree2 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 degree2 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: 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.
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