# Spatial patterns of tidal heating

###### Abstract

In a body periodically strained by tides, heating produced by viscous friction is far from homogeneous. The spatial distribution of tidal heating depends in a complicated way on the tidal potential and on the internal structure of the body. I show here that the distribution of the dissipated power within a spherically stratified body is a linear combination of three angular functions. These angular functions depend only on the tidal potential whereas the radial weights are specified by the internal structure of the body. The 3D problem of predicting spatial patterns of dissipation at all radii is thus reduced to the 1D problem of computing weight functions. I compute spatial patterns in various toy models without assuming a specific rheology: a viscoelastic thin shell stratified in conductive and convective layers, an incompressible homogeneous body and a two-layer model of uniform density with a liquid or rigid core. For a body in synchronous rotation undergoing eccentricity tides, dissipation in a mantle surrounding a liquid core is highest at the poles. Within a soft layer (or asthenosphere) in contact with a more rigid layer, the same tides generate maximum heating in the equatorial region with a significant degree-four structure if the soft layer is thin. The asthenosphere can be a layer of partial melting in the upper mantle or, very differently, an icy layer in contact with a silicate mantle or solid core. Tidal heating patterns are thus of three main types: mantle dissipation (with the icy shell above an ocean as a particular case), dissipation in a thin soft layer and dissipation in a thick soft layer. Finally, I show that the toy models predict well patterns of dissipation in Europa, Titan and Io. The formalism described in this paper applies to dissipation within solid layers of planets and satellites for which internal spherical symmetry and viscoelastic linear rheology are good approximations.

Keywords:
Tides, Solid body - Planetary dynamics - Europa - Io -Titan

The final version of this preprint is published in Icarus (www.elsevier.com/locate/icarus)

doi:10.1016/j.icarus.2012.11.020

## 1 Introduction

Fifty years ago, William Kaula found that the heat dissipated by tidal friction within the Moon is extremely nonuniform both radially and laterally [Kaula, 1963, 1964]. At that time, nonuniform tidal heating was a rather academic subject as these variations were not observable: tidal heating in the Moon is indeed much smaller than radiogenic heating. Spatial variations of tidal heating were actually a byproduct of computing the total power dissipated in the body, a key factor in modeling orbital evolution. Kaula’s calculations were based on the microscopic (or micro) approach to tidal dissipation, which starts with the computation of viscoelastic tidal strains. A bit earlier, Munk and MacDonald [1960] had given an approximate formula for the total power dissipated by tides with a macroscopic (or macro) approach. This alternative approach does not require the computation of tidal strains: the tidal bulge lags the tidal forcing by an angle parameterizing the viscous response. Their formula, however, was limited to an incompressible homogeneous body, in contrast with the micro approach which is applicable to any model of internal structure.

New theoretical developments had to wait until 1978, when two papers significantly improved tidal heating computations. First, Peale and Cassen [1978] reconsidered both macro and micro approaches to tidal heating, correcting several errors in the various formulas. They also mapped tidal heating variations for eccentricity tides and obliquity tides, showing that tidal dissipation in a homogeneous Moon is maximum at the poles in the former case and at the equator in the latter. Furthermore they showed that the presence of a large liquid core enhances dissipation in the mantle. Second, Zschau [1978] derived a simple formula for the total dissipated power in a spherically stratified compressible body, in which the influence of the body’s internal structure appears through , the imaginary part of the gravity tidal Love number . Zschau’s formula was the first step toward reconciling the macro and micro approaches, since can be computed if the internal structure of the body is known. Tobie et al. [2005b] later applied the variational method to relate to the imaginary part of the volume-integrated strain power. They also computed the radial distribution of the dissipated power in terms of deformation functions that can be evaluated with standard methods for any spherically stratified model. I will show that the formulas of Zschau [1978] and Tobie et al. [2005b] are recovered by spatially averaging the local power obtained in the micro approach, thus bridging the last gap between the macro and micro approaches. My paper, however, is primarily about the angular distribution of the dissipated power.

Spatial variations in tidal heating became relevant when spacecraft sent back incredible surface data showing that the Galilean satellites Io and Europa undergo strong tidal deformations and heating. Tidal heating was the only explanation for Io’s volcanism [Peale and Cassen, 1978] but what was going on beneath the surface was a mystery. Maybe the distribution of volcanoes could be used in order to constrain the internal structure of the satellite? Segatz et al. [1988] suggested that dissipation occurred either in the whole mantle or mostly in a thin asthenosphere close to the surface. These two models predict completely different patterns of surface heat flux with maximum dissipation at the poles in the former case and at the equator in the latter (a mix of the two is of course possible). Galileo data have often been interpreted as favoring the asthenospheric dissipation model though the issue remains controversial [Lopes-Gautier et al., 1999; Tackley et al., 2001; Kirchoff et al., 2011; Veeder et al., 2012; Hamilton et al., 2012]. Another idea consists in using long wavelength topography as an indirect measure of the heat flux (valid if the topography is isostatically compensated), but the promising analysis of Voyager data [Ross et al., 1990] was not confirmed by Galileo measurements [Thomas et al., 1998]. Lateral variations of surface heat flux however also depend on the heat transport mechanism which could be determinant.

On icy satellites like Europa, Titan and Enceladus, spatial variations of tidal heating are interesting for other reasons. In these satellites, tidal heating can melt the ice at depth and create a global subsurface ocean. The covering icy shell varies in thickness because of spatial variations in tidal heating and solar insolation. Ojakangas and Stevenson [1989a, b] computed icy shell thickness variations on Europa with the aim of predicting nonsynchronous rotation and polar wander. Nimmo et al. [2007] and Nimmo and Bills [2010] used the same method to predict long wavelength topography on Europa and on Titan, respectively. Variations of the thickness of the icy shell (or more generally the lithosphere) also influence surface tectonics [Beuthe, 2010].

Finally, spatial variations of tidal heating are important because of their strong coupling to convection. Several convection models have used as input tidal heating predicted by spherically stratified models [Tackley et al., 2001; Tobie et al., 2003; Roberts and Nimmo, 2008]. An important limitation of this approach is the neglect of lateral viscosity variations on tidal heat production, which can be taken into account by giving up the assumption of spherical symmetry and solving simultaneously for convection and tidal dissipation [Běhounková et al., 2010; Han and Showman, 2010]. Chaos terrain on Europa could be a visible result of spatially varying tidal heating enhanced by a local drop in viscosity.

Until now, predicting spatial patterns of tidal dissipation meant computing the dissipated power at every point within the body. This laborious procedure obscures the link between the internal structure and the resulting pattern, and makes it difficult to look for all possible patterns generated by realistic internal structures. In this paper, I show that the dissipated power at a given radius within a spherically stratified body is the linear combination of three basic patterns that depend only on the tidal potential. The coefficients weighting the patterns depend are radial functions which can be computed with standard methods developed for tidal deformation problems once the internal structure of the body has been specified.

I study the influence of the internal structure on dissipation patterns by computing the dissipation weight functions in various toy models without assuming a specific rheology. The toy models are the thin icy shell above an ocean, the homogeneous body (relevant to a completely solid body with little stratification or to a solid core) and the incompressible two-layer body, either with a liquid core or with rigid core, the latter case being relevant to dissipation in a soft layer (such as an asthenosphere) above a more rigid layer. It is well-known that dissipation patterns are completely different if dissipation occurs in the deep mantle and in a thin asthenosphere. Besides these two classes of patterns, I show that dissipation in a thick asthenosphere leads to a third type of dissipation pattern, with maximum heating at the equator as in asthenospheric dissipation but with a lower content in harmonic degree four. Toy models predict well dissipation patterns within real bodies though not their magnitude. I will give three examples of this by computing dissipation weight functions for realistic internal structures of Europa, Titan and Io.

## 2 Dissipated power

### 2.1 Power, strains and tidal potential

In the micro approach to dissipation, the dissipated power within the planet or satellite is expressed in terms of tidal strains (this formula is derived in Appendix A and compared with other expressions found in the literature). If tides operate at only one angular frequency , the dissipated power per unit volume averaged over one orbital period is given by

(1) |

where (resp. ) is the complex shear (resp. bulk) modulus, is the Fourier transform of the strain tensor and is the trace of . All quantities implicitly depend on the frequency and on the point within the planet where the power is evaluated. In this paper, the tilde on viscoelastic parameters indicates that they are complex and frequency-dependent (the tilde is dropped if the parameters are purely elastic).

If there are several tidal frequencies, the total power averaged over time is a sum over these frequencies (interferences vanish, see Eq. (68)) so that it becomes essential to know the frequency dependence of the viscoelastic parameters (i.e. the rheology). Unfortunately the rheology of planetary bodies is poorly constrained [Jackson, 2007; Karato, 2008]. Earth’s mantle has been mainly studied at frequencies that are much higher (laboratory experiments), moderately higher (seismic attenuation and seismic anisotropy) or much lower (Chandler wobble and postglacial rebound) than tidal frequencies [Karato, 2010]. It is indeed difficult to determine Earth’s viscous response at tidal frequencies [e.g. Benjamin et al., 2006; Nakada and Karato, 2012] and even more so for other bodies. Maxwell rheology is often used in studies of tidal dissipation because it is the simplest model in which the response changes from elastic to viscous as the frequency decreases. It is however not clear how the Maxwell viscosity is related to the true viscosity of the material [Ross and Schubert, 1986; Bills et al., 2005; Sotin et al., 2009]. Rheological models depending on more parameters such as the Andrade model [Castillo-Rogez et al., 2011] or the extended Burgers model [Nimmo et al., 2012] could be more realistic. Moreover there has been a long-standing debate on whether viscous deformations in Earth’s mantle are mainly due to diffusion creep or to dislocation creep, corresponding to a linear or nonlinear rheology, respectively [Karato and Wu, 1993]. In this paper, I assume that the rheology is linear without being more specific about it except in applications to real bodies for which I use the Maxwell model. Besides I consider only the dominant tidal frequency; this restriction is appropriate for Solar System bodies, but it should be lifted for exoplanets with short orbital periods.

Dissipation within the Earth is generally much larger in shear than in uniform compression (), at least in the seismic frequency band [Anderson, 1989]. For this reason, most models of planetary dissipation assume that the bulk modulus is real and independent of frequency [Ranalli, 1995, p. 142], though it is not necessarily true in presence of a high fraction of partial melt [Schmeling, 1985]. In Earth’s asthenosphere, the ratio of bulk to shear dissipation could possibly reach 30% [Durek and Ekström, 1995]. Bulk dissipation should thus be kept in mind when studying bodies (such as Io) where asthenospheric dissipation could be the dominant mechanism.

The strains appearing in Eq. (1) are induced by tidal deformations which are in turn related to the tidal potential. The tidal potential at the surface of the deformed body can be expanded in spherical harmonics of degree and order [Kaula, 1964], each component being the superposition of several terms of angular frequencies ,

(2) |

where means real part of , is the colatitude and is the longitude in a frame attached to the body. The coefficients are defined by this equation if you know the tidal potential from the formulas in Kaula [1964]; they are not the same as the complex coefficients of the Fourier series of since the frequencies are not all different. If the orbital eccentricity and the obliquity of the body are not zero, an infinite number of terms (indexed by ) contribute to the potential at a given degree and order . I neglect all tidal perturbations except the dominant contribution of degree and I sum over the order . The component of the tidal potential at frequency can be written as

(3) | |||||

where means complex conjugate.

Strains are related to the tidal potential through displacements. If the internal structure of the body is spherically symmetric, displacements due to tides can be written as [Alterman et al., 1959; Takeuchi and Saito, 1972; Saito, 1974]

(4) |

where ) are the Fourier components of the displacement defined as in Eqs. (62).

The two radial functions and depend on the internal structure of the body and have the dimension of inverse acceleration. These functions are part of a set of six radial functions that solve six coupled linear differential equations of first order (I drop from now on the explicit dependence on ). Strains depend on , and their derivatives. If derivatives are eliminated with the help of the differential equations, strains also depend on the functions and respectively associated with the stresses and (or ). The remaining functions and are related to the gravity potential. The derivatives are related to by Eqs. (82) of Takeuchi and Saito [1972]:

(5) |

The quantity actually measures the shear strain (see Eqs. (8) below). In an incompressible medium, the first equation reduces to

(6) |

Since stresses with a radial component vanish at the surface , the boundary conditions for and are

(7) |

Beware that the definitions of vary between authors (I follow here Takeuchi and Saito [1972]). In particular, the functions of Sabadini and Vermeersen [2004] are equivalent to the functions of Takeuchi and Saito [1972].

Kaula [1963, 1964] already used the functions in his pioneering papers on tidal dissipation; his example was recently followed by Tobie et al. [2005b] and Roberts and Nimmo [2008]. This formalism has the advantage that it is widely used in geophysics to compute the deformation of Earth’s surface for complicated internal structures [Spada et al., 2011]. Peale and Cassen [1978] unfortunately returned to the method of Love [1911] which is difficult to apply correctly to structures more complicated than two layers of the same density.

In spherical coordinates, I divide strain components in three classes: (1) the purely radial component is the radial strain, (2) the purely angular components are the tangential strains, and (3) the mixed components are the radial-tangential shear strains. I now insert Eqs. (4) into the strain-displacement equations in spherical coordinates [Takeuchi and Saito, 1972, Eqs. (20)]. Note that the non-diagonal strains in Takeuchi and Saito [1972] must be multiplied by if they are to be the components of the strain tensor. The radial and shear strains read

(8) |

Doing the same for the tangential strains, I get

(9) |

where are differential operators of degree two on the sphere defined by Eqs. (80)-(81).

Though Eqs. (8)-(9) are well-known, their tensorial structure with respect to rotations is always ignored in geophysics: is a scalar, are the components of a vector and are the components of a second-order symmetric tensor (Greek indices denote or ). The tensorial character of strains under rotations is the key to simplifying the formula for the dissipated power.

### 2.2 Strain invariants

From a technical point of view, this is the key section of the paper. In order to evaluate the dissipated power in terms of the tidal potential, I need to compute in Eq. (1) the quantities and which are invariant under arbitrary coordinate transformations [Malvern, 1969]. This general invariance is not obvious once the strains have been expressed in spherical coordinates. Nevertheless invariance under rotations should remain obvious in that basis. It is indeed easy to construct rotationally invariant terms: take the product of two scalars (such as ), take the scalar product of two vectors (such as ), contract two tensors (such as ), or take the trace of a tensor.

Let us define

(10) |

with the subscript dilat denoting dilatation, since is the infinitesimal change of volume.

In spherical coordinates, reads

(11) |

The terms (scalar) and (trace of a tensor) are each invariant under rotations but it is easier to keep them together in Eq. (11). is the sum of three terms that are each invariant under rotations:

(12) |

where

(13) |

, and are invariant because they are respectively the product of two scalars, the scalar product of two vectors and the contraction of two tensors.

Inserting Eqs. (8)-(9) and the first Eqs. (82) into Eq. (11), I write the first invariant as

(14) |

where is the spherical Laplacian defined by Eq. (76). Inserting Eqs. (8) and (79) into the first and second Eqs. (13), I write the radial and shear invariants as

(15) |

Inserting Eqs. (9) and (82) into the third Eqs. (13), I write the tangential invariant as the sum of two terms invariant under rotations:

(16) |

where

(17) |

### 2.3 Factorized power

The last step consists in rewriting the dissipated power given by Eq. (1) in terms of the various strain invariants of Section 2.2:

(19) |

Examining the factorized strain invariants (Eqs. (18)), we see that the power at a given radius depends on the angles through a linear combination of three factors,

(20) |

showing that there are at most three independent spatial patterns.

The three invariants without derivatives of can be combined as

(21) |

The three functions (20) are not the best choice of basic angular functions because of compensations between the three terms. Instead I will maintain as much as possible the distinction between the radial, radial-tangential shear and tangential dissipative terms. These contributions have indeed very different magnitudes in very different interior models (see Section 4.3). Besides , I thus propose to use as basic angular functions the combinations appearing in and (see Eqs. (18)):

(22) |

These functions are dimensionless since the factor absorbs the dimension of the squared tidal potential ( is the mean motion and is the surface radius, see Eq. (89)). The numbers 12, 22 and 48 result from setting in Eqs. (18). The patterns and are normalized so that the terms without derivatives are the same in the three functions. In the following, the index collectively refers to .

Inserting Eqs. (18), (21) and (22) into Eq. (19), I write the dissipated power as

(23) |

where the weight functions and are positive and depend on the internal radial functions :

(24) |

The weights functions have the dimension of squared inverse acceleration. When averaging over angles, I will need the sum of the weight functions denoted as

(25) |

In practice, is evaluated from , and (see Eqs. (5)). The terms , and are associated with radial-tangential shear, tangential and bulk dissipation, respectively, whereas the term is a combination of radial and tangential dissipation. The terms and depend on the same angular function and can be grouped into one term. Formulas (22), (23) and (24) are the central results of this paper.

At the surface, weight functions can be related to the displacement Love numbers which characterize how the body responds to tidal forcing. The displacements and are indeed proportional to the tidal Love numbers and , while and vanish because of Eqs. (7):

(26) |

where is the surface gravity. Moreover, the first of Eqs. (5) evaluated at the surface gives

(27) |

where I replaced and by Young’s modulus and Poisson’s ratio :

(28) |

I compute the weight functions at the surface with Eqs. (24) to (27):

(29) |

The amplitude of tidal deformation (quantified by ) varies by orders of magnitude depending on the global structure of the body. By contrast, the ratio typically ranges from 0.2 to 0.3 [Beuthe, 2010, Fig. 10]. Eqs. (29) show that ratios of weight functions at the surface are insensitive to the magnitude of the viscoelastic response . In good approximation, this property also holds inside the body. It is thus useful to define dimensionless rescaled weight functions by

(30) |

For example, rescaled weight functions do not depend on the rheology if the body is incompressible and homogeneous (Section 4.2) or if it is consists of an incompressible viscoelastic mantle above a liquid core of the same density (Section 4.3).

Weight functions have a few properties that do not depend on a specific internal structure. First, always vanishes at the surface (see Eqs. (29)) and at internal fluid/solid interfaces (where a condition similar to Eq. (7) holds). Pattern B can thus be disregarded close to such boundaries, for example in a thin icy shell above an ocean. Second, and vanish in the incompressible limit at the boundary of an infinitely rigid layer because there. Dissipation thus approximately follows Pattern B close to rigid/viscoelastic boundaries, for example at the bottom of an icy shell in contact with a silicate mantle. Third, I can relate to at the surface with Eqs. (29):

(31) |

This ratio is less than one-half for silicates () and less than one-fifth for ice (). If the material is incompressible, it is possible to go further by applying Eq. (6):

(32) |

so that is expected to be much smaller than at all radii if incompressibility is a good approximation. Since dissipation also receives contributions from terms B and C, the contribution of bulk dissipation to total dissipation remains in general minor even if (see discussion in Section 2).

## 3 Spatial patterns, global power and surface flux

### 3.1 Harmonic content of angular functions

The angular functions defined by Eqs. (22) can be computed from the expansion of the squared norm of the tidal potential in spherical harmonics. Since the tidal potential is of degree two, the expansion of its squared norm only includes functions of degrees zero, two and four:

(33) |

where

(34) |

The constant is the spatial average of the squared norm of the nondimensional tidal potential,

(35) |

while and represent spatial variations with zero average. Table 1 gives the numerical values of the coefficients for a satellite in synchronous rotation with non-zero eccentricity and obliquity.

eccentricity tides () | ||||||
---|---|---|---|---|---|---|

obliquity tides () | 0 | |||||

interference () |

In terms of the harmonic functions , the angular functions read

(36) |

The angular dependence cancels in the following combination:

(37) |

Eqs. (36) give us a qualitative understanding of the spatial patterns. and mainly differ by the sign of the degree-four term. By contrast, the term of degree two in is of opposite sign with respect to (or ) and the degree-four contribution to is much smaller, even though contains derivatives of the fourth order. The harmonic content of the angular functions can be measured with the variance (see Eq. (94)) if one knows the tidal potential. Table 2 shows that the content of the variance in degree four is much larger for obliquity tides than for eccentricity tides though is still mainly of degree two. Fig. 1 shows the basic patterns of tidal heating (angular functions ) for a body in synchronous rotation with either eccentricity tides or obliquity tides, drawn from Eqs. (36) in which I substitute Eq. (34) and the values of Table 1.

unit | ||||
---|---|---|---|---|

eccentricity tides | 30 | 44 | 1 | % |

obliquity tides | 64 | 76 | 5 | % |

The relative importance of the three patterns can be determined from the comparison of the weight functions but one should keep in mind that the do not have the same variance. I thus define equal-variance weight functions through

(38) |

The standard deviations are tabulated in Table 3 for synchronous rotation (eccentricity tides only or obliquity tides only). For eccentricity tides, and . For obliquity tides, .

eccentricity tides () | |||
---|---|---|---|

obliquity tides () |

Instead of , I could use as basic angular functions and express the spatially varying part of the power as a linear combination of only two functions ( and ). The radial weights associated with and , however, are not always positive. For example, the degree-two component of the power has a minimum either at the poles or at the equator, depending on the sign of the radial function. Physically, these two cases are interpreted as different spatial patterns.

### 3.2 Global power

I show here that the dissipated power integrated over the volume of the body is equal to the global power computed with the macro approach to tidal heating discussed in the Introduction. The angular average of the local power at radius can be read from Eqs. (23), (25) and (36):

(39) |

and have the same dependence on the functions as the radial sensitivity functions and in Tobie et al. [2005b] (my formulas are more compact but strictly equivalent). For a satellite in synchronous rotation undergoing eccentricity tides, (see Table 1) and , in which case is identical to the function of Tobie et al. [2005b] (their Eq. (37)) except for a difference in sign explained below. After having shown that the global power is equal to , these authors correctly interpret as the radial distribution of the dissipation rate per unit volume averaged over angles.

Using variational principles, Tobie et al. [2005b] prove that conservation of energy relates the sum of the strain and kinetic energies to the potential energy, the former two being integrated throughout the body whereas the latter is computable at the surface. The dissipated part of this energy balance results in

(40) |

where is the radius of the tidally perturbed body, is the tidal gravity Love number of degree two and is the gravitational constant. I obtain the correct sign in the right-hand side of Eq. (40) by correcting Eqs. (34)-(35) of Tobie et al. [2005b] as follows: must be replaced by its complex conjugate in their Eq. (34) so that the left-hand side of their Eq. (35) has an additional minus sign.

The local power integrated over the volume of the body is thus given by

(41) | |||||

This equation agrees with the formula that Zschau [1978] obtained for the energy dissipated during one tidal cycle by considering the dephasing between the tidal potential and the potential induced by tidal deformations (the formulation in terms of was introduced by Platzman [1984]). Therefore Eqs. (39), (40) and (41) prove the equivalence between the micro and macro approaches to tidal heating.

Substituting the values of Table 1 into Eq. (41) and setting , I get the well-known formula for the total dissipated power within a satellite in synchronous rotation [Cassen et al., 1980; Segatz et al., 1988; Chyba et al., 1989],

(42) |

which is valid up to second order in eccentricity and obliquity (terms proportional to vanish when averaged over angles). In the macro approach to tidal heating, the total power is often expressed in terms of the quantity which Segatz et al. [1988] define through

(43) |

I inserted a minus sign in this equation because is usually assumed to be positive whereas is negative. There is however no simple relation between this and the specific dissipation function of the same name [Munk and MacDonald, 1960], as discussed by Zschau [1978] and Segatz et al. [1988]. Contrary to what is sometimes said in the literature, Eq. (42) is valid when the body is inhomogeneous (but spherically stratified) and compressible. The macroscopic derivation of this equation does not even require a linear rheology. Note that existing extensions of Eq. (42) to arbitrary eccentricity and obliquity [Wisdom, 2008; Levrard, 2008] assume that is proportional to frequency. While this assumption makes computations easy, it does not correspond to any plausible rheology [Efroimsky and Lainey, 2007]. Eq. (42) thus remains the best available formula for the total dissipated power within a satellite in synchronous rotation if the eccentricity and obliquity are not too large.

In the Solar System, the obliquity of large satellites other than the Moon has not yet been observed except for Titan [Stiles et al., 2008, 2010] but it is theoretically predicted to be very small [Gladman et al., 1996; Peale, 1999; Bills, 2005; Baland et al., 2012]. Thus obliquity tides are at present negligible for large satellites except the Moon where they contribute about 40% of the tidal heating. Heating patterns are however not observable on the Moon because radiogenic heating is dominant. Though Moon’s obliquity may have been larger in the past, tidal heating due to obliquity tides was always smaller than radiogenic heating [Peale and Cassen, 1978; Wisdom, 2006]. Regarding exoplanets, Winn and Holman [2005] suggested that obliquity tides account for bloated ‘hot Jupiters’ but this idea has been criticized on the grounds that the corresponding Cassini state is unstable [Fabrycky et al., 2007; Levrard et al., 2007; Peale, 2008]. Tides could also be caused by forced librations in longitude [Wisdom, 2004], the effect of which can be computed by adding new terms proportional to in the tidal potential (Eq. (88)).

### 3.3 Surface flux

Observations from space give indications on the distribution of surface heat flux but not on the dissipated power at depth. Computing angular variations of the surface heat flux thus involves assumptions about heat transport. Let me first define as the fraction of the global power due to the dissipation term ,

(44) |

with a similar formula for ( and ). By definition, .

If the heat is radially transported to the surface, the surface heat flux due to tidal dissipation is given by

(45) |

where is the spatially-averaged heat flux at the surface,

(46) |

Under the assumption of radial heat transport, the coefficients are the weights for the three basic angular patterns of the surface flux. Non-radial transport mechanisms within the planet tend to average lateral variations in heating. On the other hand, an analysis including a more realistic modeling of heat transport within the planet (for example by convection) should also take into account lateral variations in viscosity due to tidal heating itself. I further discuss this topic in the Conclusions. Finally the heat passing through the crust can be channelled into heat pipes (for Io) or fractures (Enceladus’ tiger stripes), thus further modifying the distribution of surface heat flux.

## 4 Interior structure and spatial patterns: toy models

### 4.1 Thin shell

#### 4.1.1 Rheology constant with depth

If the body contains an ocean beneath a thin icy shell (or an asthenosphere beneath a lithosphere), displacements do not vary much with depth within the shell and can be approximated by their value at the surface which depend on the displacement Love numbers and . In the membrane limit of thin shell theory, these Love numbers are related by

(47) |

where is the viscoelastic counterpart of Poisson’s ratio. This relation was previously derived for an elastic thin shell [Beuthe, 2010, Eq. (D.12)]. Using the correspondence principle, I apply it here to a viscoelastic thin shell. The membrane limit of thin shell theory is appropriate for loads of large wavelength with respect to the shell thickness : the harmonic degree of the load should satisfy [Beuthe, 2008, Eq. (83)], that is for tidal loads of degree two. This condition is satisfied in thin shell theory, in which one usually assumes that the shell thickness is less than 10% of the radius.

When does the thin shell assumption break down for Eq. (47)? Consider an incompressible () two-layer body of uniform density in the liquid core limit for which the ratio is given by Eq. (114) with . Deviations from the thin shell limit () are smaller than 10% if . This constraint coincides with the usual requirement of thin shell theory mentioned above and is satisfied in most models of icy satellites with a subsurface ocean.

Eq. (47) is not very sensitive to rheology. For materials constituting planetary crusts, ranges between 0.25 (silicates) and 0.325 (elastic ice, Gammon et al. [1983]). The value of depends on the rheological model. If the rheology is Maxwell, the real part of varies between the elastic value (high viscosity) and the incompressible value 0.5 (low viscosity), while its imaginary part is always much smaller than its real part. Neglecting the imaginary part, we see that changes from to (less than 10%) as changes from to .

I compute weight functions in the thin shell approximation with Eqs. (25), (29) and (47):

(48) |

In absence of bulk dissipation, the relative weights of the patterns do not depend on the rheology: Pattern C dominates in the icy shell ( for eccentricity tides) whatever the value of . Of course, the amount of tidal heating depends on the rheology of ice through the common multiplying factor depending on and . Therefore, the dissipated power in a thin shell above a fluid layer (e.g. an icy crust above an ocean) is well approximated by

(49) |

where is related to the ratio of bulk dissipation to shear dissipation:

(50) |

Fig. 2 shows the resulting spatial pattern for eccentricity tides and obliquity tides if there is no bulk dissipation (see Section 5.4 for a comparison with the literature).

#### 4.1.2 Love numbers

In the formula for the dissipated power (Eq. (49)), the Love number does not affect the spatial pattern but regulates the total amount of tidal heating. It can be computed with an interior model, for example the thin shell limit of the incompressible body of uniform density (see Eqs. (119)). I give here slightly more general formulas (obtained with the same method) valid for an incompressible body of nonuniform density having a fully rigid mantle surrounded by an ocean and a thin icy shell. The density and elastic structure below the mantle do not matter: a solid or liquid core of higher density may be present. The ocean and the icy shell have the same density. In this thin shell/rigid mantle approximation, the Love numbers are given by

(51) |

The parameters , and are the shell thickness, the density of the ice (or ocean) and the average density of the body, respectively. The factor is the value of in the limit of vanishing shell thickness:

(52) |

Eqs. (51) expanded at first order in agree with Eqs. (3)-(6) of Wahr et al. [2006]. It is however important not to expand Eqs. (51) if one wants to compute good estimates of the imaginary part of Love numbers. Though Wahr et al. [2006]’s formulas take into account the density difference between ice and ocean, they are not well-suited to a viscoelastic interpretation because is not necessarily of the same order as . Using Eqs. (49) and (51), I can check that the integration over the volume of the shell of the dissipated power yields the global power (Eq. (42)) as computed with the macro approach (equivalently, one can verify that Eq. (40) is satisfied).

#### 4.1.3 Depth-dependent rheology

The above formulas implicitly assumed that the icy shell has a homogeneous rheology characterized by . Actually, the assumption of depth-independent rheology is not realistic. The viscosity of ice indeed depends on the local temperature of the ice and thus varies a lot between the top ( for Europa) and bottom () of the icy shell, where it is at its melting point. This variation has a huge effect on the shear modulus , as can be seen for example in a Maxwell model [Wahr et al., 2009]. However, Poisson’s ratio varies much less as I already argued in Section 4.1.1. In practice, one divides the shell into an outer conductive layer and an inner convective layer, the former being mainly elastic and the latter being viscoelastic [McKinnon, 1999]. Dissipation mainly occurs in the convective layer whereas the conductive layer influences deformations.

What is the effect of depth-dependent rheology on the thin shell formula for the dissipated power? There is no problem is letting be depth-dependent in Eq. (49), but is a constant so that cannot vary with depth in the formulas for Love numbers, Eqs. (47) and (51). In this last equation, we see that is constant () and that is multiplied by the thickness of the shell . In thin shell theory, it can be shown that the product in Eqs. (51) arises from the integration of the shear modulus over the thickness of the shell [Beuthe, 2008]:

(53) |

In good approximation, (resp. ) is determined by the conductive (resp. convective) layer where (resp. ) is highest. Thus the real (resp. imaginary) part of the Love numbers, corresponding to deformation (resp. dissipation), is mainly determined by the conductive (resp. convective) layer. The substitution preserves energy conservation: one can either check that Eq. (40) is satisfied, or more explicitly that the integration over the volume of the shell of the dissipated power with depth-dependent yields the global dissipated power (Eq. (42)). Incompressibility () is the simplest assumption for compatible both with depth-dependent rheology and with the absence of bulk dissipation. In contrast with Eq. (53), it would not make sense to simply replace by its average in Eq. (47). Note that the averaging procedure described by Eq. (