The Tides of Enceladus' Porous Core.

The inferred density of Enceladus' core, together with evidence of hydrothermal activity within the moon, suggests that the core is porous. Tidal dissipation in an unconsolidated core has been proposed as the main source of Enceladus' geological activity. However, the tidal response of its core has generally been modeled assuming it behaves viscoelastically rather than poroviscoelastically. In this work, we analyze the poroviscoelastic response to better constrain the distribution of tidal dissipation within Enceladus. A poroviscoelastic body has a different tidal response than a viscoelastic one; pressure within the pores alters the stress field and induces a Darcian porous flow. This flow represents an additional pathway for energy dissipation. Using Biot's theory of poroviscoelasticity, we develop a new framework to obtain the tidal response of a spherically symmetric, self-gravitating moon with porous layers and apply it to Enceladus. We show that the boundary conditions at the interface of the core and overlying ocean play a key role in the tidal response. The ocean hinders the development of a large-amplitude Darcian flow, making negligible the Darcian contribution to the dissipation budget. We therefore infer that Enceladus' core can be the source of its geological activity only if it has a low rigidity and a very low viscosity. A future mission to Enceladus could test this hypothesis by measuring the phase lags of tidally induced changes of gravitational potential and surface displacements.

Introduction

When Voyager 1 and Voyager 2 flew by Enceladus they revealed a surprisingly young and active world (Smith et al., 1982). More than two decades later, the Cassini spacecraft discovered water plumes erupting from Enceladus' South Polar Terrain (SPT) (Porco et al., 2006), showed that the SPT radiates ∼10 GW of energy to space (Howett et al., 2011), and demonstrated that Enceladus has a global subsurface ocean (Postberg et al., 2011; Thomas et al., 2016).

Explaining the high geological activity of Enceladus remains challenging (Nimmo et al., 2018). The moon's activity is linked to diurnal tides (e.g., Yoder, 1979). Enceladus is currently in a 2 : 1 mean‐motion orbital resonance with Dione that forces its orbital eccentricity, causing time‐varying tides that periodically deform the moon. As Enceladus is not perfectly elastic, part of the tidal energy is transformed into heat, a process known as tidal heating. Energy dissipation in Enceladus is ultimately dependent on tidal dissipation in Saturn. Dissipation within the planet causes a phase‐lag in Saturn's tidal bulge; consequently, rotational energy is transferred to Enceladus and Dione where part of it is dissipated (e.g., Nimmo et al., 2018).

Astrometric observations of the Saturnian system can be used to constrain the phase‐lag of Saturn's tidal bulge and estimate the amount of energy dissipated in the moons (Meyer & Wisdom, 2007). They suggest that Enceladus is in orbital and thermal equilibrium (Fuller et al., 2016; Lainey et al., 2012). If the moon is in thermal equilibrium, ice‐shell thickness estimates combined with measurements of the SPT thermal flux can be used to obtain the total energy produced within the moon, which adds up to ∼35 GW (Hemingway et al., 2018). However, explaining where and how this much energy is dissipated within Enceladus has been problematic (e.g., Nimmo et al., 2018), giving rise to Enceladus' energy puzzle.

Enceladus' ice shell is most likely brittle and conductive, limiting the amount of heat that can be dissipated within it to about 1 GW (Beuthe, 2019; Souček et al., 2019). Frictional heating along Enceladus' tiger stripes can contribute an additional 0.1–1 GW of energy dissipation (Pleiner Sládková et al., 2021), but, overall, tidal heating in the ice shell can only account for roughly 10% of the observed SPT thermal output. Ocean tides have been proposed as an additional heating mechanism (Tyler, 2011), but they only become important if Enceladus has an orbital obliquity two orders of magnitude higher than the expected value (Chen & Nimmo, 2011); or the ocean is unrealistically thin, radially stratified or turbulent (Chen et al., 2014; Hay & Matsuyama, 2019; Matsuyama, 2014; Rekier et al., 2019; Rovira‐Navarro et al., 2019, 2020; Tyler, 2020; Wilson & Kerswell, 2018). Modeling the core as a purely solid, viscoelastic body, tidal dissipation produced within it can only account for the observed thermal output if the core has a viscosity of η  < 1013 Pa s, much lower than that characteristic of rock. Because of this, substantial tidal dissipation in the core was first disregarded. However, Roberts (2015) and Choblet et al. (2017) recently suggested that a low enough viscosity can be attained if Enceladus' core is porous.

An Enceladan porous core is consistent with observations. The density of the core inferred from gravity data, 2.4 g cm−3 (e.g., Beuthe et al., 2016), is low compared to that of the minerals expected to form the bulk of the core (Choblet et al., 2017). Furthermore, the detection of salt‐rich particles (Postberg et al., 2009, 2011), silicon‐rich nanoparticles (Hsu et al., 2015) and molecular hydrogen (Waite et al., 2017) in material ejected by Enceladus' plumes suggests that the ocean interacts with the silicate core in hydrothermal systems. Taken together, these observations suggest that the moon's core is a porous, water‐saturated matrix of silicates or loosely packed rock pieces through which water can circulate.

Even though Roberts (2015) and Choblet et al. (2017) attributed Enceladus' activity to a porous core, they did not explicitly model how a porous core responds to tides. Instead, they modeled the core as a viscoelastic, rocky solid. The response of a porous, permeable, water‐saturated body to tidal forces differs from that of a pure solid. The deformation of the matrix induces a flow of water through the permeable interior, which in turns affects the response of the solid matrix and modifies the dissipation in the solid. Furthermore, the viscous flow of water through the pores adds an additional source of dissipation that may not be negligible.

While the tidal response of solid and liquid layers have been thoroughly examined (e.g., Beuthe, 2016; Chen et al., 2014; Jara‐Orué & Vermeersen, 2011; Kaula, 1964; Love, 1911; Matsuyama et al., 2018; Renaud & Henning, 2018; Rovira‐Navarro et al., 2019; Segatz et al., 1988; Tyler, 2008), the tidal response of bodies with porous layers has been subjected to much less scrutiny. Wang et al. (1999) estimated energy dissipation due to tidally induced flows in Earth's permeable seafloor and showed it to be negligible; Vance et al. (2007) applied the same approach to Enceladus' hydrothermal systems and reached a similar conclusion for the icy moon. Liao et al. (2020) developed a more complete approach based on Biot's theory of poroviscolasticity (Biot, 1941) and argued that the interaction between solid and liquid phases lead to a heat production that can easily exceed 10 GW and thus solve Enceladus' energy puzzle.

The model presented by Liao et al. (2020) included several simplifications that require further examination: (a) only the tidal response of the core was considered instead of that of the whole moon (core, ocean and ice shell); (b) the authors forced the problem via an imposed surface strain derived from viscoelastic models and only considered one component of the eccentricity tide instead of forcing the core with the complete tidal potential; (c) the authors neglected the effect of self‐gravity, the body force arising from the tidal deformation itself; and (d) they assumed that the displacement field was irrotational.

In this paper, we relax the assumptions of Liao et al. (2020) and develop a self‐consistent model to compute the tidal deformation and fluid flow of self‐gravitating bodies with porous layers (Section 2) that can be applied to Enceladus and other bodies with internal porous layers. The new approach is an extension of the standard theory of tides for self‐gravitating viscoelastic bodies (Love, 1906; Peltier, 1974; Sabadini et al., 2016; Saito, 1974; Takeuchi et al., 1962) to bodies with poroviscoelastic layers. We apply the new model to Enceladus (Section 3), and examine the circumstances under which sufficient tidal dissipation can be produced within the core to explain the thermal energy radiated by the moon. In Section 4 we conclude that Darcy dissipation is likely negligible in Enceladus' thermal budget, leaving a low‐viscosity, low‐rigidity core as the most plausible avenue for substantial tidal dissipation in the rocky core. We propose how future missions could test this hypothesis.

Methods

Our aim is to obtain the linear, periodic, tidal response of a body with internal porous layers to a tidal perturbation. We assume the moon is composed of N spherically symmetric layers of uniform properties. The boundary between a layer i and i + 1 is at radius r . We consider layers that are either purely liquid, purely solid, or a contiguous solid matrix with a permeable network of liquid‐filled pores. We use a viscoelastic model of the solid. The model can include as many layers as required to approximate the interior structure of the moon under consideration. Figure 1 shows the interior structure that is thought to be valid for Enceladus, consisting of an icy shell, a subsurface ocean and a porous core.

Figure 1

Interior structure of Enceladus consisting of three layers: a porous core, a subsurface ocean and an ice shell. The porous core boundary is assumed to be permeable.

Liao et al. (2020) applied Biot's theory of poroviscoelasticity to the tidal problem. However, the present work is the first derivation of a self‐consistent model of body tides in porous media; we therefore provide a detailed formulation of the problem, highlight key assumptions, and explain how it differs from the standard tidal theory for viscoelastic solids and the work of Liao et al. (2020). To this end, we start by presenting the governing equations for porous media in detail, show how they reduce to those of a pure solid, and discuss how internal liquid layers are modeled (Section 2.1); then we introduce the tidal potential and boundary conditions (Section 2.2); and, finally, sketch the solution method (Section 2.3), leaving full details of the mathematical formulation to the Appendices.

Governing Equations

We use the volume‐averaged mass and momentum conservation equations for a parcel of the moon. The parcel contains a solid and a liquid phase of densities ρ and ρ and volumes V and V , respectively. The parcel is at least one order of magnitude bigger than the typical grain size. This way, we can use continuum mechanics rather than explicitly modeling the microphysical interactions of grains within the parcel. Thus, the variables that follow should be understood as averages. The porosity of the parcel is simply defined as the ratio between the liquid (V ) and total (V) volumes:

The mass conservation equations for the liquid and solid phases can be written as (e.g., Ganesan & Poirier, 1990) Equations 2a and 2b can be added to obtain ρ is the bulk density, ρ = Φρ  + (1 − Φ)ρ ; and is the segregation flux given by the relative velocity of the liquid phase with respect to the solid phase,  = Φ(  −  ). We additionally introduce the variation of fluid content, defined as the amount of liquid entering the solid frame per unit of the solid frame (Cheng, 2016), with being a displacement vector.

The momentum equation in a frame rotating with the moon’s angular velocity (ω) is given by (McKenzie, 1984) is the stress tensor, ϕ is a potential that includes gravitational forces and the centrifugal force. is an interaction force between the solid and liquid phases given by p is the pore pressure and η is the liquid viscosity. κ is the matrix permeability, which depends on the geometry of the solid matrix. For a solid matrix made of uniform, spherical grains of size d , a commonly used expression is the Kozeny‐Carman law (Carman, 1997; Kaviany, 1995):

The inertial terms in the momentum equations (Equation (5a), (5b)) can be neglected. For the solid phase, the high viscosity and the long period of tidal forces as compared to seismic waves imply that the solid is in quasi‐equilibrium. In a porous medium, the interaction force is generally larger than inertial terms in Equation 5b, which in turn results in a small Reynolds number. Summing the two momentum equations, an equation for the bulk or total stress is obtained with the total stress tensor being We assume that deviatoric stresses in the liquid fluctuate on the pore scale and hence they volume‐average to zero, except for their contribution to the interaction force (Equation 6). The stress tensor of the liquid phase is thus isotropic and given by the pore pressure, with being the identity matrix. Using the interaction force expression, the liquid phase momentum equation reduces to Darcy's law, Note that if the porosity is zero, we recover the mass and momentum conservation equations for a solid.

A constitutive equation relating the stress tensor and pore pressure to kinematic variables is needed. We define the strain tensor as where is the volume average of the solid displacement ( ). The stress and strain tensors can be split into a mean and a deviatoric component: with ϵ  = tr( )/3 and P = −tr( )/3; tensile stress is taken to be positive.

In a compressible poroviscoelastic solid, deformation is associated with an effective stress ′ (Biot, 1941; Cheng, 2016), α is Biot's constant, the meaning of which will become evident later. We consider that the material is viscoelastic. Different rheological laws can be used to consider the behavior of a viscoelastic material (e.g., Renaud & Henning, 2018). We consider the Maxwell model, in which case effective stress and strain are related as or, in terms of total stress, Here, μ and η are the shear modulus and viscosity of the two‐phase aggregate, respectively. The Maxwell model does not capture the anelastic behavior of ices and silicates, which can become especially important when the forcing period is much smaller than the Maxwell time η/μ (Efroimsky, 2012). It is for this reason that recent studies have considered the Andrade rheological model (Andrade & Trouton, 1910) for ices (e.g., Castillo‐Rogez et al., 2011; Gevorgyan et al., 2020; Rambaux et al., 2010; Rhoden & Walker, 2022; Shoji et al., 2013) and silicates (e.g., Bierson & Nimmo, 2016; Efroimsky, 2012; Renaud & Henning, 2018; Rovira‐Navarro et al., 2021; Walterová & Běhounková, 2017). In Appendix F we discuss how more complex rheology models can be incorporated into our theory and demonstrate this using the Andrade model.

One further constitutive equation relates the pore pressure with the isotropic strain and the variation of fluid content as (Cheng, 2016), where K and K are the bulk modulus of the material in undrained (ζ = 0) and drained (p = 0) conditions, respectively. These are effective properties of the two‐phase medium; the drained modulus depends on the mechanical properties of the solid matrix; in contrast, the undrained modulus depends on both the properties of the liquid and the solid phases. If the material is microscopically homogeneous and isotropic, K, K and α can be obtained using the bulk modulus of the solid K and liquid K phases, and the bulk modulus of porosity K Φ, which measures the resistance to grain rearrangement, (e.g., Cheng, 2016), With this definition, it becomes apparent that the Biot parameter α relates the resistance to compression of the solid constituent and the porous matrix. A strong porous matrix (e.g., spherical holes) has a small α while an easily deformable matrix (i.e., slit or crack‐shaped pore spaces) has an α close to 1. If the solid is much less compressible than the frame, K Φ/K → 0 and therefore α → 1.

The previous set of constitutive equations define a material that, upon a stress perturbation, exhibits an elastic response and viscous creep. Moreover, a perturbation produces a pore pressure field that drives Darcian porous flow. Viscous creep and Darcian flow result in energy dissipation. The rate of volumetric tidal dissipation averaged over a tidal cycle due to these two processes is (e.g., Liao et al., 2020) respectively.

If the body has an internal liquid layer (i.e., ocean), alternate equations are required for that layer. We assume internal liquid layers are inviscid, incompressible and in hydrostatic equilibrium and the radial displacements follow equipotential surfaces, except at solid‐liquid interfaces, where this might be hindered (Jara‐Orué & Vermeersen, 2011). is the radial unit vector, and g the gravitational acceleration. Under these assumptions, surface (e.g., Hay & Matsuyama, 2017; Matsuyama, 2014; Rovira‐Navarro et al., 2020; Tyler, 2011) and internal waves (Rekier et al., 2019; Rovira‐Navarro et al., 2019) are excluded from the solution.

In all layers, the gravitational potential of the body can be computed using Poisson's equation,

Tidal Forcing and Boundary Conditions

We consider a synchronously rotating moon of radius R with an orbital frequency ω and eccentricity e. As the obliquity of Enceladus is expected to be very small (Chen & Nimmo, 2011), we focus on eccentricity tides and ignore obliquity tides. The tidal potential at a point with co‐latitude and longitude θ, φ located at radial distance r from the center of the moon is given by (e.g., Jara‐Orué & Vermeersen, 2011; Kaula, 1964) are normalized, real spherical harmonics of degree l and order m (Equation A2).

To solve the previous set of equations, boundary conditions at the moon's surface are required. The normal and shear stress at the surface are zero, the potential, ϕ, is continuous at the surface but its gradient is not. Using Poisson's equation and applying Gauss' theorem for an infinitesimal control volume surrounding the surface layer, we find We note that Liao et al. (2020) considered only the zonal component of the tidal potential (m = 0). Furthermore, the tidal forcing was imposed via a prescribed strain at the surface of the core instead of via the tidal potential ϕ as explained here, and the no‐stress boundary conditions (Equation 23) were not used. The distinct effect of each of these boundary conditions will be explored in Section 3.

Additional boundary conditions must be prescribed at internal boundaries. For the core–ocean and ocean–ice shell interfaces, we use the boundary conditions discussed in Jara‐Orué and Vermeersen (2011) and given in Appendix A. Nevertheless, an additional boundary condition should be provided at the porous layer interface (r ). Two different boundary conditions can be considered: no radial Darcy flux or force balance and continuity of fluid pressure. In the latter case, the ocean pressure at the core surface is balanced by the radial component of the stress tensor and the pore pressure equals the ocean pressure, which implies p(r ) + σ (r ) = 0. For Enceladus, we consider that the core–ocean boundary is permeable and thus use the second boundary condition.

Perturbation Theory and Solution Method

To make the analysis tractable, we linearize the equations of motion using perturbation theory. We split the density, stress tensor and potential force into a background and a perturbed component. The background component corresponds to the pre‐stressed, hydrostatic state arising from the self‐gravity of the body; the perturbed component is the result of the time‐dependent tidal potential, Here, ϕ 0 is the gravitational potential of the unperturbed body and ϕ Δ includes both the perturbing tidal potential and the potential arising from self‐gravitation of the perturbed body. In the unperturbed state, Equations 8 and 11 are given by: with ∇ϕ 0 = g .

Using the previous definitions and linearizing by assuming that the products of perturbation variables are negligible, the momentum equations can be written as:

The mass conservation equation can be obtained by linearizing Equation 3 and assuming that under small displacements the Lagrangian and Eulerian derivatives are approximately equal, Similarly, using the definition of the variation in liquid content (Equation 4), and the segregation flux, we obtain: Finally, the density change of the liquid phase can be obtained using the definition of the liquid bulk modulus,

The perturbed gravitational potential is obtained by solving the linearized Poisson's equation,

Equations (29a), (29b), (30), (31), (32), (33) reduce to those used in Liao et al. (2020) if the solid and liquid phases are assumed to be massless (ρ 0, ρ  = 0). If the layer is purely solid (Φ = 0, α = 0), we recover the classic equations used for a viscoelastic solid (Sabadini et al., 2016).

To obtain the tidal response of the body, the momentum Equation (29a), (29b), mass conservation Equations 30 and 31, and Poisson's Equation 33 together with constitutive Equations 16 and 17 should be solved under appropriate boundary conditions (Section 2.2). As the tidal forcing is periodic, we solve the equations of motion in the Fourier domain. We assume a solution proportional to exp(iωt) and transform the previous set of equations to the Fourier domain. Because of the symmetry of the problem, we solve the previous set of equations using spherical harmonics. We obtain stress and strain tensors, the pore pressure and Darcy flow, and we compute tidal dissipation in the solid and liquid phases using Equation (19a), (19b). Further details can be found in Appendices A, B, C, D.

Application to Enceladus

To understand how the predictions of this model differ from the previous treatment of Liao et al. (2020) we consider three different cases: (a) Enceladus' core forced via a prescribed surface strain as in Liao et al. (2020); (b) Enceladus' core with a free surface and forced with the tidal potential; and (c) a complete model of Enceladus consisting of a porous core, an ocean, and an ice shell, forced with the tidal potential. We begin by examining the simpler cases 1 and 2 to illustrate the effect of the boundary conditions used by Liao et al. (2020) (Section 3.1) and then move to the more complex, multilayered model to show how the ocean and ice shell affect the core's tidal response (Section 3.2).

We assume a core density of 2.4 g cm−3 consistent with gravity observations (Beuthe et al., 2016). Choblet et al. (2017) obtained a core porosity of 20–30% for realistic core compositions; we use a value of 20%. For a consolidated silicate core, the shear modulus is ∼1–10 GPa and the viscosity is ∼1020 Pa s or higher at low homologous temperature. However, if the core is unconsolidated, it can become weaker and the shear modulus and viscosity can be orders of magnitude lower than the typical values of silicates (Choblet et al., 2017; Goldreich & Sari, 2009; Nimmo et al., 2018). The parameters used are summarized in Table 1.

Table 1

Enceladus Physical and Mechanical Properties

Quantity	Symbol	Value	Units
Surface radius	R	252.1	km
Mass	M	1.08 ⋅ 1020	kg
Ocean thickness a	h ocean	38	km
Ice shell thickness a	h ice	23	km
Average core's density a	ρ core	2422	kg m−3
Ocean's density	ρ ocean	1000	kg m−3
Ice viscosity b	η ice	1 ⋅ 1018	Pa s
Ice shear modulus b	μ ice	3.3	GPa
Ice bulk modulus	K ice	33	GPa
Core shear modulus	μ s	0.01–10	GPa
Core solid phase bulk modulus c	K s	10	GPa
Core viscosity	η s	1010–1020	Pa s
Biot's constant	α	0–1	−
Water viscosity c	η l	1.9 ⋅ 10−3	Pa s
Water bulk modulus c	K l	2.2	GPa
Core permeability c	κ	10−8 − 10−4	m2
Core porosity	Φ	0.2	‐
Eccentricity	e	0.0047	‐
Orbital Period	T	33	h

Beuthe et al. (2016).

Hussmann and Spohn (2004).

Liao et al. (2020).

Core‐Only Model

Liao et al. (2020) studied the response of Enceladus' core to a prescribed radial strain imposed at the core's surface (in this section, R = r ). The strain was given by a degree‐2 order‐0 field of the form: . The amplitude of the strain was estimated from the theory of viscoelastic tides (Murray & Dermott, 2000). For a homogeneous body in which rigidity dominates over self‐gravity, the maximum radial strain attained at the poles is , where is the gravitational potential Love number and is the complex rigidity. Additionally, Liao et al. (2020) assumed the displacement field to be irrotational, ∇ ×  = 0, and considered the solid and liquid to be massless (ρ, ρ ) = 0. For Case 1 we make the same assumptions.

Liao et al. (2020) found that for Biot parameter α → 1, tidal dissipation is enhanced as compared to standard viscoelastic models. In what follows, we assess whether this still holds when the core is forced via the tidal potential rather than via a prescribed surface strain. We first take a fixed core permeability (κ = 10−8 m2), and compute tidal dissipation for different values of core viscosity and Biot parameter α. Changing α for a given porosity is equivalent to changing the ratio between the bulk modulus of porosity and the solid bulk modulus (Equation (18a), (18b)), as α approaches 1 the porous matrix becomes more compressible. Afterward, we study the role of the core's permeability. While we keep the porosity fixed, we note that variations in porosity between 0.2 and 0.3 result in changes of α and κ much smaller than the ranges explored below (Equations 7 and (18a), (18b), (18c)).

Figure 2 shows total tidal dissipation in the core for different values of Biot parameter α and core viscosity η for cases 1 and 2. For Case 1, we reproduce the results of Liao et al. (2020). For α = 0.95 tidal dissipation features two peaks, one at a core viscosity of ∼1011 Pa s, also characteristic of the viscoelastic response, and another one attained at higher core viscosity ∼1015 Pa s, only characteristic of the poroviscoelastic model. Around the two peaks, most of the energy dissipation occurs in the solid phase, as shown by the thin lines. As core viscosity increases further, dissipation in the solid decreases but the total dissipation remains high due to Darcy dissipation. This dissipation occurs in a shallow layer close to the core's surface (Figure 3c), where a strong pressure gradient develops that drives flows of up to 2 × 10−5 m s−1. The second dissipation peak is the result of the compressibility of the porous matrix. As α → 1, the drained bulk modulus decreases (Equation (7), (18a), (18b), (18c)), the second peak becomes more prominent, and Darcy dissipation also increases.

Figure 2

Tidal dissipation in Enceladus' core when forced with a prescribed surface strain (solid lines) or the tidal potential (dashed line) for different values of Biot coefficient α. The thin lines indicate tidal dissipation in the solid layer, the gray line shows the amount of Darcy dissipation for incompressible solid and liquid phases obtained using Equation (7), (18a), (18b), (18c).

Figure 3

Total tidal dissipation in an Enceladan core with a free surface for various values of core viscosity and permeability. In (a) the core is forced via a prescribed strain field of order degree 2 and order 0, in (b) via the tidal potential. The contour for GW is indicated in both plots. (c) and (d) show tidal dissipation in the liquid and solid phases for a meridional cut at longitude 0° for the two points indicated in (a) and (b). Both points have a viscosity of 1018 Pa s and result in the same amount of tidal dissipation (28 GW). For all cases we assume μ = 1 GPa, K  = 10 GPa, K  = 2.2 GPa and α = 0.95.

For Case 2, forcing by the tidal potential, the second dissipation peak is not present (heavy dashed line in Figure 2). Tidal dissipation reaches its maximum at the same core viscosity as in the non‐porous, viscoelastic case and then decreases as viscosity increase until Darcy dissipation becomes dominant. However, as opposed to (1), the amount of heat resulting from Darcy dissipation is less sensitive to the porous‐matrix compressibility and, more importantly, it is severely reduced. The prominent pressure gradients characteristic of Case 1 do not develop and the maximum flow velocities attained are reduced by two orders of magnitude. This suggests that explaining Enceladus' thermal budget in terms of poroviscoelastic dissipation may be more problematic. However, a highly permeable core may mitigate this to some extent.

The permeability dictates how easily water can flow through the core. Figures 3a and 3b show the total amount of internal heat production for cases 1 and 2 for different values of core permeability. As before, for both cases we observe the high dissipation band characteristic of a viscoelastic core with a low viscosity (∼1011 Pa s). As core viscosity increases Darcian dissipation becomes dominant and the total tidal dissipation becomes independent of core viscosity. In this regime, dissipation increases with permeability and the flow velocity is controlled by the dimensionless number Ω  = ωR 2 η /κμ, which can be understood as a ratio between the timescale of Darcy flow (R 2 η /κμ) and of the tidal perturbation (1/ω). When Ω ≪ 1, high flow velocities are attained (, Equation 11) which in turn results in high values of tidal dissipation (, Equation (7), (18a), (18b), (18c)). We note that when Ω becomes very high, some of the terms in the equations (e.g., the term A 87 in Equation (7), (18a), (18b), (18c)) can become very large, causing numerical problems. This limits the lowest value of permeability we can attain under our current formulation to ∼10−9 m2. Nevertheless, we derive an analytical expression for Darcian dissipation in the limit of an incompressible porous matrix and liquid (Appendix E) that presents good agreement with the numerical results (Figure 2) and shows that Darcian dissipation can be expected to further decrease for lower permeability values.

While dissipation increases with permeability in both cases, Darcy dissipation in the tidally forced case is lower than in the case where a surface strain is prescribed. In Case 1) a high surface stress follows from the imposed strain, producing a large gradient in pore pressure and driving high‐amplitude Darcian flow; in contrast, in Case 2 the no‐stress boundary conditions prevent this from occurring. To produce an amount of heat similar to that observed, a permeability of κ > 10−5 m2 is required (Figure 3b). These are high permeability values compared to the permeability of Earth's hydrothermal systems, which can reach values of about 10−8 m2 (Lauer et al., 2018). However, it is possible that Enceladus' core does not resemble such a system, but is instead akin to an unconsolidated rubble pile. In that case, Enceladus' core would be made up of loosely packed material through which water can easily circulate.

If Enceladus' core has a porosity of ∼0.2–0.3, a permeability of ∼10−5 m2 requires grain sizes of about ∼10–50 cm (Equation 7). This blocky structure could be the relic of a violent formation process such as Enceladus forming after a series of collisions of a previous generation of moons (Asphaug & Reufer, 2013; Ćuk et al., 2016).

Multi‐Layered Model

We now move to the more realistic multi‐layered model consisting of a porous core, a subsurface ocean and an ice shell (Case 3). Although Enceladus' ice shell is of variable thickness (e.g., Beuthe et al., 2016; Hemingway & Mittal, 2019; Čadek et al., 2016), we consider an ice shell of constant thickness equal to its average value. For the core radius, core density, and ice and ocean thicknesses we use values in agreement with gravity, shape and libration data (Beuthe et al., 2016). We keep these values and the rheological properties of the ice constant and vary the properties of Enceladus' core (Table 1).

As in the previous cases, we compute the tidal response of the moon for different values of core viscosity and permeability. The results are shown in Figure 4; they demonstrate that the presence of an ocean reduces the tidal response of the core. As in Cases 1 and 2, we find a peak in tidal dissipation for a core viscosity of ∼1011 Pa s. Around this value, tidal heating is compatible with Enceladus' thermal output. However, the high‐dissipation band is narrower than for Cases 1 and 2—the presence of an ocean and an ice crust reduces the tidal deformation of the core (Beuthe, 2015). Most importantly, the amount of Darcian dissipation is drastically reduced.

Figure 4

Total tidal dissipation in Enceladus' core for different values of core viscosity and permeability. We assume μ = 1 GPa, K  = 10 GPa, K  = 2.2 GPa and α = 0.95. The contour for GW is indicated.

As opposed to Case 2, Darcian dissipation is small even for a highly permeable core. The reduction in Darcian dissipation is due to the presence of the overlying ocean and ice shell. They impose a non‐zero pressure at the core surface that largely balances the forcing of the tidal potential. The Darcian flow q is driven by the modified pressure p* = p Δ + ρ ϕ Δ (Equation (29a), (29b)). Without an ocean, the pressure at the surface p Δ is 0 and the tidal force drives the Darcian flow, p* = ρ ϕ Δ. In contrast, if there is an overlying ocean and ice shell, the pressure at the core–ocean boundary is not zero. To understand this, consider a core that is just covered by an ocean. Under our assumptions, the ocean surface follows the equilibrium tide −ϕ Δ/g while the core–ocean boundary has a radial displacement of u. As a result, the pressure perturbation at the core–ocean boundary is p Δ = ρ g(−u − ϕ Δ/g). This ocean pressure partially compensates the driving tidal potential, leading to a modified pressure p* proportional to the the core radial displacement −ρ gu. Writing this in terms of h 2, the radial displacement Love number, p* = ρ h 2 ϕ Δ. Because of its rigidity, the radial displacement of the core is generally much smaller than the equilibrium tide (h 2 ≪ 1), explaining why Darcy dissipation is greatly reduced when considering an ocean‐covered core.

The presence of an ice shell above the ocean inhibits the ocean surface from following the equilibrium tide by imposing a pressure load. If the ice shell is completely rigid (no surface displacements), the pressure load compensates the equilibrium tide and is given by −ρ ϕ Δ. Hence the ice shell imposes a pressure equivalent to that imposed by the water column in an ice‐free ocean, leading to the same reduction in effective pressure at the core. Appendix E provides analytical solutions for an incompressible core with and without a free surface that further demonstrate the role of the pressure at the core–ocean boundary.

An Unconsolidated Core

In light of the reduced tidal dissipation in Enceladus' porous core in comparison to the findings of Liao et al. (2020), we reconsider the low viscosity band where dissipation is enhanced. In all cases, we found that tidal dissipation peaks at ∼1011 Pa s. Close to this viscosity value, the amount of heat dissipated within the core is compatible with Enceladus' observed thermal output. Such a low viscosity is incompatible with the viscosity characteristics of silicates (Roberts & Nimmo, 2008). However, Choblet et al. (2017) proposed that low effective viscosity values can be attained if Enceladus' core is unconsolidated. If this is the case, friction between grains can give rise to substantial dissipation. We next reconsider this hypothesis and discuss its plausibility and remaining unknowns.

If strains are small, a granular material essentially behaves as a monolith. The grains deform elastically and stresses are transmitted at the grain boundaries. Some viscous deformation can also occur mainly due to diffusion creep. Under higher strains, grain–grain sliding becomes important and, combined with inter‐granular friction, it can result in enhanced energy dissipation (Lambe & Whitman, 1969). Instead of using μ and η to characterize the rheology, a granular material is normally characterized in terms of an effective shear modulus (μ ) and a damping coefficient () (e.g., Choblet et al., 2017; Seed et al., 1986). These two variables are related to the shear modulus and the solid viscosity introduced before as

Figure 5 shows the amount of tidal dissipation in the core as function of μ and υ. For υ > 0.1 more than 10 GW can be generated in Enceladus' core, provided the shear modulus is low enough. The values of the damping coefficient and effective shear modulus depend on factors that include the amplitude of the deviatoric strain, the confining pressure, and the forcing frequency (Faul & Jackson, 2005; Lambe & Whitman, 1969; Seed et al., 1986).

Figure 5

Total tidal dissipation in Enceladus' core for different values of core effective shear modulus and damping coefficient. A viscoelastic core and ice shell are assumed. The red lines indicate the maximum deviatoric strain attained within the core.

Goldreich and Sari (2009) showed that the effective rigidity of an unconsolidated body is smaller than that of a monolith due the concentration of stresses in sharp contact points. They proposed that the shear modulus is controlled by the curvature radius at these contact points, which in turn depends on the yield strain of the material ϵ , and showed that the effective rigidity of an unconsolidated body can be estimated as , with R and g the surface radius and gravity and ρ the bulk density. This expression is consistent with an increase in shear modulus with confining pressure seen in the laboratory (Goddard, 1990). Using values representative of Enceladus's core (μ ∼ 1 GPa, ϵ  = 10−2), we find μ  ∼ 0.7 GPa—a value much higher than that required to attain high values of tidal dissipation for moderate values of damping coefficient υ (Figure 5).

Laboratory experiments can also be used to bound υ and μ . If the material experiences high deviatoric strains, υ increases. Laboratory data shows a pronounced increase of υ for ϵ > 0.01%; υ can reach values higher than 0.15 for ϵ > 0.1% (e.g., Seed et al., 1986; Rollins et al., 1998). For Enceladus' tidal amplitude, these high strains are only attained if the material has a low effective shear modulus (μ  ∼ 107 − 108 Pa) (Figure 5). The shear modulus increases with increasing overburden pressure and decreases with the amplitude of the deviatoric strain. For small strains (ϵ < 0.001%), the shear modulus of typical sand mixtures at Enceladus' core pressure (5–50 MPa) is on the order of 109 Pa (Seed et al., 1986), incompatible with enhanced dissipation. Strains on the order of 0.01%–0.1% can reduce the effective shear modulus by around 50% (Rollins et al., 1998; Seed et al., 1986). Furthermore, the effective shear modulus is expected to decrease and the damping coefficient to increase at lower forcing frequencies (Faul & Jackson, 2005). Unfortunately, laboratory data is only available for a frequency range (0.01–1 Hz) much higher than Enceladus' tidal frequency (∼10−5 Hz). It remains to be seen if such changes are sufficient to access the high‐dissipation region of Figure 5. Laboratory experiments at Enceladus‐like conditions (high confining pressure, low forcing frequency) are required to assess whether Enceladus' core is in the highly deformable state required for the generation of the observed thermal activity.

A more direct measurement of Enceladus' core viscosity could be provided by a future Enceladan mission. Due to the viscosity of the core and the ice shell, their tidal responses are characterized by phase lags with respect to the forcing. The ice‐shell and core phase lags affect the gravity field and the surface displacements of the moon in different ways. Because the ocean decouples the ice shell from the core, the phase lag of the moon's surface displacement is mostly dependent on the viscosity of the ice shell; in contrast, the phase lag of the moon's gravity field depends on both the viscosity of the core and the ice shell. A low‐viscosity core leads to a large gravity phase lag but has a much smaller effect on the surface‐displacement phase lag; in contrast, a low viscosity ice shell produces phase lags in the gravity field and surface displacements of similar magnitude. Therefore, by measuring the difference between gravity and surface displacement phase lags, we can distinguish between a low‐viscosity and a high‐viscosity core. Hussmann et al. (2016) proposed this strategy to constrain Europa's core viscosity and the viscosity of high‐pressure ice layers within Ganymede; a similar technique could be used for Enceladus (Marusiak et al., 2021).

The gravity and surface‐displacement phase lags are given by the phase lags of the gravitational and radial displacement Love numbers k 2 and h 2. Figure 6 shows the difference in gravity and surface‐displacement phase lags, . If the core has a low viscosity, phase‐lag differences up to 50° are attained. This also holds if the more complex Andrade rheology is considered (Appendix F). The gravity phase lag could be measured by precise tracking of a single or dual orbiter around Enceladus (e.g., Ermakov et al., 2021). The surface‐displacements phase lag could be measured using a laser or radar altimeter (Steinbrügge et al., 2015, 2018). Such measurement would help to constrain Enceladus' core viscosity and settle the long‐standing puzzle of where Enceladus' heat is coming from.

Figure 6

Difference in gravity and surface‐displacements phase lags for different values of core and ice shell viscosity. The thick and thin black contours show combinations of parameters for which tidal heating is 10 and 1 GW, respectively. A shear modulus of 1 GPa is assumed for the core; the properties of the ice shell are those given in Table 1, except for the viscosity, which we vary. Maxwell rheology is assumed for both ice shell and core.

Conclusions

Several observations suggest that Enceladus' core likely is a porous silicate matrix throughout which water can permeate. The tidal response of a porous core is expected to be markedly different from that of a non‐porous one: the presence of pores renders the body more compressible, pressure within the pores can alter the stress field, and water can flow within the core adding an additional dissipation mechanism. For a non‐porous core, tidal dissipation is only high if the core has a rigidity and a viscosity significantly lower than those expected for a monolithic silicate core. Recently, Liao et al. (2020) presented an analysis of the tidal response of Enceladus' porous core using Biot's theory of poroviscoelasticity and showed that poroviscoelastic effects can increase tidal dissipation for core properties compatible with those of silicates.

Liao et al. (2020) cautioned that their model relies in some assumptions that required further scrutiny: they considered the tidal response only of the core and ignored the overlying ocean and ice shell, neglected the effects of self‐gravity, and forced the problem via a prescribed surface strain rather than via the tidal potential. In this paper, we extended the model of Liao et al. (2020) by combining the theory of poroviscoelasticity (e.g., Biot, 1941; Cheng, 2016) with the theory commonly used to obtain the tidal deformation of viscoelastic, self‐gravitating bodies (e.g., Love, 1906; Sabadini et al., 2016).

With this theory, we assessed the relevance of the assumptions made by Liao et al. (2020), and studied the subspace of core properties for which tidal heating can explain Enceladus' thermal output. We started by considering a model of Enceladus' core without an overlying ocean and ice shell and showed that the boundary conditions at the core boundary play a central role. If the core is forced via a prescribed radial strain —as done in Liao et al. (2020)— tidal dissipation in both the solid and liquid phases can be orders of magnitude higher compared to that given by standard viscoelastic models. However, if a free surface is assumed (no‐stress) and the core is forced via a tidal body force, dissipation in both liquid and solid phases is reduced. In this scenario, Enceladus' thermal output can only be explained by tidal dissipation in the solid phase if the core viscosity is very low (η  ∼ 109–1012 Pa s) compared to that characteristic of silicates or via Darcian dissipation in the fluid if the core is highly permeable (κ > 10−5 m2). We then considered a more realistic multi‐layered model consisting of a porous core, a subsurface ocean and an ice‐shell. We showed that the presence of a hydrostatic ocean hinders the tidal response of the core. More importantly, tidal dissipation due to Darcy flow is severely reduced, making it difficult to reconcile tidal dissipation in Enceladus' rocky core with the observed thermal output.

Dissipation within the solid phase can still account for Enceladus' thermal output if the core is weak and has a low viscosity. This requirement appears to be incompatible with a monolithic, silicate core. Yet Choblet et al. (2017) ascribed the low shear modulus and viscosity to grain–grain friction in a fragmented core akin to a rubble pile. While it is true that a low viscosity can arise from this process, it only occurs if the tidal strain is sufficiently large, which requires a low shear modulus that is difficult to reconcile with laboratory experiments (Rollins et al., 1998; Seed et al., 1986). However, available laboratory data is not representative of Enceladus' core conditions, which points to the need for further laboratory work. Alternatively, we show that a future Enceladus mission could probe the core's viscosity by measuring the phase lag of tidally induced changes in the gravity field and surface deformation.

Other tidally active worlds might also have porous regions. Vigorous tidal heating can partially melt the mantle of a planet or a moon forming a porous sublayer filled with magma. Galileo's magnetometer data suggests that Io has a partially molten asthenosphere with a melt fraction of 20% (Khurana et al., 2011) (although alternative explanations have been proposed for the magnetometer data (Blöcker et al., 2018; Roth et al., 2017)). Io is the closest example of a magma‐rich world (Khurana et al., 2011; Peale et al., 1979; Spencer et al., 2020), but extrasolar worlds experiencing higher levels of tidal heating might also be common (e.g., Peters & Turner, 2013; Rovira‐Navarro et al., 2021). So far, attempts to compute the tidal response of bodies with a partially molten asthenosphere have either used the theory of viscoelasticity (e.g., Fischer & Spohn, 1990; Segatz et al., 1988), or the Laplace tidal equations commonly employed to model ocean tides (Hay et al., 2020; Tyler et al., 2015). Our model opens the door to study the tidal response of these worlds in a new light.