Probing the icy shell structure of ocean worlds with gravity-topography admittance

The structure of the icy shells of ocean worlds is important for understanding the stability of their underlying oceans as it controls the rate at which heat can be transported outward and radiated to space. Future spacecraft exploration of the ocean worlds (e.g., by NASA's Europa Clipper mission) will allow for higher-resolution measurements of gravity and shape than currently available. In this paper, we study the sensitivity of gravity-topography admittance to the structure of icy shells in preparation for future data analysis. An analytical viscous relaxation model is used to predict admittance spectra given different shell structures determined by the temperature-dependent viscosity of a tidally heated, conductive shell. We apply these methods to the ocean worlds of Europa and Enceladus. We find that admittance is sensitive to the mechanisms of topography support at different wavelengths and estimate the required gravity performance to resolve transitions between these mechanisms. We find that the Airy isostatic model is unable to accurately describe admittance universally across all wavelengths when the shell thickness is a significant fraction of body's radius. Our models suggest that measurements of admittance at low spherical harmonic degrees are more sensitive to thick shells with high tidal dissipation, and may complement ice-penetrating radar measurements in constraining shell thickness. Finally, we find that admittance may be used to constrain the tidal dissipation within the icy shell, which would be complementary to a more demanding measurement of the tidal phase lag.


INTRODUCTION
The existence of subsurface oceans within icy satellites is of great interest due to their potential habitability. In this paper, we focus on two such ocean worlds: Jupiter's moon Europa and Saturn's moon Enceladus. The presence of a global subsurface ocean on Europa has been inferred from observations by Galileo of an induced magnetic field. A global salty ocean is believed to be the conductive fluid causing the induced magnetic field (Khurana et al. 1998). The presence of liquid water in Enceladus was inferred from observations of water vapor plumes (Porco et al. 2006). Thomas et al. (2016) inferred the global nature of the ocean from observations of the large amplitude of physical libration.
Considering life as we know it, the presence of large amounts of liquid water, essential elements from chondritic material at the base of the ocean, energy from tides and radiogenic sources as well as chemical gradients make ocean worlds promising astrobiology targets (Hand et al. 2009). A major consideration for assessing habitability is the persistence and stability over geologic timescales of a liquid water ocean to allow for developing and sustaining life. The structure of the overlying icy shell, its thickness, heat transport mechanism, and the spatial distribution of tidal heating are particularly important for understanding the overall heat budget and the long-term survivability of the ocean.
For Europa, the structure of the icy shell has been studied using classical Airy isostasy model and observations of topography (e.g., Kadel et al. 2000), flexural analysis (e.g., Billings & Kattenhorn 2005;Nimmo et al. 2011), and tidal-convective equilibrium model (e.g., Moore 2006), but large uncertainty remains for the mean shell thickness. For Enceladus, shell thickness was constrained by gravity, shape and libration data (e.g., Hemingway & Mittal 2019). However, for both Europa and Enceladus the properties of the icy shell that govern heat transport, such as viscosity and temperature profile, remain poorly constrained.
In this paper, we explore the use of gravity-topography admittance to study the icy shells of ocean worlds Europa and Enceladus. Gravity-topography admittance Z n is defined as a wavelength-dependent ratio of gravity to topography spectral amplitudes: Z n = gravity amplitude shape amplitude at spherical harmonic degree n [mGal/km], where gravity and shape data are obtained from spacecraft observations. Gravity-topography admittance can be computed from the spherical harmonic expansion coefficients of the gravity and shape up to the degree of the lowerresolution data set. Typically, the accuracy and resolution of admittance are limited by the quality of the gravity data set. Current observations of gravity and shape of Europa and Enceladus from the Galileo and Cassini spacecraft, respectively, are limited to long-wavelength measurements. In terms of spherical harmonic expansions, gravity fields have been measured up to spherical harmonic degree 2 for Europa (spatial scale of 4900 km) (Anderson et al. 1998;Casajus et al. 2020) and up to degree 3 for Enceladus (Iess et al. 2014) (spatial scale of 528 km). Nimmo et al. (2007) made ellipsoidal fits of Europa's shape, and Enceladus' shape has been mapped up to spherical harmonic degree 16 (Tajeddine et al. 2017) (spatial scale of 100 km). Future missions such as the upcoming NASA's Europa Clipper will deliver higher-resolution measurements of the gravity and shape allowing determination of gravity-topography admittance over a larger range of spatial scales. Our goal is to prepare for these future higher-resolution data by exploring the sensitivity of gravity-topography admittance to various aspects of the icy shell structure.
Gravity-topography admittance bears clues to the icy shell structure, specifically to topography support mechanisms. The dominant topography support mechanism could vary depending on the spatial scale. The Airy isostasy model, used extensively for the Earth (e.g., Watts 2001), describes surface topography supported by buoyancy force arising due to a topography of a density-discontinuity interface (i.e., a crustal root) located at a certain depth called "depth of compensation". Topographic loads can also be supported by elastic and viscous stresses within the body. Since the interface between the icy shell and ocean is a phase boundary, melting and freezing can produce non-hydrostatic basal topography inducing a flow throughout the shell (Čadek et al. 2019). These different topography support mechanisms influence the amplitudes of topography at different wavelengths at the surface and the base of the shell, which affects the moon's gravity field and, therefore, is reflected in the admittance spectrum.
In summary, the goals of the paper are: 1. To provide an algorithm for computing gravity-topography admittance suitable for icy shells with large gradients of viscosity.
2. To explore the sensitivity of gravity-topography admittance to viscosity profile, shell tidal heating, and shell thickness.

METHODS
In our exploration of the sensitivity of admittance to icy shell structure, we follow the process laid out in Fig. 1. We model Europa and Enceladus as three spherically symmetric layers, placing an icy shell on top of a liquid water ocean that overlies a solid mantle. For each choice of shell thickness, the density of the mantle and thickness of the ocean are computed to satisfy the total mass, radius, and moment of inertia factor of the moon as given in Table 1. In addition, viscoelastic parameters shown in Table 1 are assigned to each layer and are used in the tidal heating model. Starting with this icy shell structure, a temperature profile is found by solving the conductive heat equation with tidal heating as a source term. Since tidal heating depends on a temperature-dependent ice rheology, we compute the temperature profile iteratively by alternating calculations of tidal heating and heat conduction to arrive at a steady-state solution. We use this converged temperature profile to determine the shell's viscosity profile. The viscosity profile is then used in a viscous relaxation model to find the shape amplitudes at the surface and the base of the icy shell, treating the ice-ocean boundary as either a material or a phase boundary. Finally, the relation between the shell's surface and base amplitude at each spherical harmonic degree is used to obtain the expected gravity-topography admittance spectrum. We will now proceed with a detailed description of these steps.  The rheologic model of ice governs its response to applied stresses. We assume that the icy shell flows viscously on geologic timescales. The properties of ices under conditions relevant to icy moons are poorly understood, largely due to the difficulty of reproducing low stress, frequency, and temperature conditions with laboratory experiments (≤ 0.1 MPa, Tobie et al. (2003)). The flow mechanisms of diffusion creep and grain boundary sliding (GBS) are most relevant for the low stress and strain rates, low temperatures, and small grain sizes expected in icy moons (Goldsby & Kohlstedt 2001). The diffusion creep flow mechanism, extrapolated to these low stresses from experimental data and described by Goldsby & Kohlstedt (2001), results in a Newtonian flow, in which strain rate is proportional to stress. GBS has lower activation energy compared to diffusion creep and, unlike diffusion creep, is non-Newtonian, which introduces a stress dependence of viscosity. Diffusion creep likely dominates at the conditions and grain sizes in the icy shell (Moore 2006). Thus, we assume viscosity is independent of stress and the flow is Newtonian described by the diffusion creep deformation mechanism (e.g., Showman & Han 2004;Tobie et al. 2003;Mitri & Showman 2005). Under these assumptions, the temperature-dependent viscosity of pure water ice is described by:

Assumptions of ice rheology
where E a is the activation energy of diffusion creep taken to be 59.4 kJ mol −1 (Goldsby & Kohlstedt 2001), T m is the melting temperature of ice taken to be 273 K, and R = 8.314 J mol −1 K −1 is the universal gas constant. The grain size dependency is included within the viscosity at the melting point η melt . The grain size is estimated to be in the range from 0.1 to 1 mm (McKinnon 1999;Kirk & Stevenson 1987), corresponding to η melt between 10 13 Pa s and 10 15 Pa s. There are several simplifications in our relaxation modeling. We do not model solid-state convection within the shell. Convection, if it occurs, is most likely confined to the bottom part of the shell and is less likely to occur in thinner shells (e.g., McKinnon 1999). The inferred large amplitude of basal shell topography at Enceladus argues against convection as convection would likely rapidly relax this topography (Hemingway & Mittal 2019). In addition, the rheology of ice at the base of the shell near the melting point can be influenced by premelting and partial melting that would reduce viscosity enhancing the flow (Tobie et al. 2003). Additionally, impurities within the ice, such as salts and silicates, can affect grain sizes and increase viscosity (Barr & Showman 2009).
At the surface, applying Eq. 2 at the equilibrium surface temperatures of Europa (92 K) or Enceladus (59 K) implies that ice near the surface has a viscosity on the order of 10 30 and 10 50 Pa s, respectively. Such high viscosities will prevent the relaxation of surface topography on geologic timescales. Passey (1983) and Bland et al. (2012) find that craters on Enceladus are highly relaxed despite these low surface temperatures. Passey (1983) determine a range of surface viscosities on Enceladus between 10 24 and 10 25 Pa s to explain the relaxation of the observed craters, and suggest that an insulating layer at the surface could increase the effective surface temperature and thus decrease viscosity. On Europa, Showman & Han (2004) argue that fractures in ice and brittle deformation could affect ice rheology at the surface and parametrize this more complex rheology by also imposing an upper viscosity bound. We follow this simplification of a bounded viscosity used by Showman & Han (2004) andČadek et al. (2017 for our viscous relaxation model and consider values η bound = 10 22 − 10 25 Pa s. We compute tidal heating (Section 2.2) and the temperature structure of the shell (Section 2.3) using unbounded viscosity, and the resulting viscosity profile will subsequently be set to η bound where η(r) > η bound to compute the relaxation of the shell (Section 2.4).

The effect of tidal heating
As an icy moon orbits its parent planet, tidal forces cause a periodic deformation of the moon resulting in heat production. In this paper, we introduce tidal heating within the icy shell through its influence on the temperature profile and, thus, viscosity profile through Eq. 2. We use the Maxwell rheology to model viscoelastic deformation of a spherically symmetric, layered body under the influence of tidal potential. For a chosen shell thickness and viscosity profile, our layered model moon is described by the viscoelastic and orbital parameters listed in Table 1. While the mantle and ocean are treated as single layers with constant viscoelastic parameters, the icy shell is broken up into 80 layers prescribed by the viscosity profile. The degree 2 tidal potential expanded to first order in eccentricity for a tidally-locked satellite is used to set the surface boundary condition. This potential drives the deformation of the viscoelastic shell causing tidal heating. We follow Tobie et al. (2005) to compute a radial profile of surface averaged volumetric tidal dissipation rate, h tide (r). We use h tide (r) as a source term in the heat conduction equation in the following Section 2.3. A full description of the tidal dissipation model and the computation of h tide (r) is given in Appendix A.
Tidal dissipation is maximized in the portion of the shell where Maxwell time, defined as the ratio of viscosity to elastic shear modulus, is close to the forcing period (e.g., Tobie et al. 2003). The viscosity for which this occurs can be expressed as η Maxwell = µ/ω. For tidal forcing at Europa's orbital frequency ω = 2.0 × 10 −5 rad s −1 and taking elastic shear modulus µ = 3.3 GPa, maximum dissipation occurs at viscosity η Maxwell = 1.6 × 10 14 Pa s. For Enceladus, a higher orbital frequency ω = 5.3 × 10 −5 rad s −1 leads to maximum dissipation at a lower viscosity η Maxwell = 6.2 × 10 13 Pa s. As discussed in Section 2.1, the viscosity at the base of the shell η melt is taken as a range from 10 13 to 10 15 Pa s covering a range of grain sizes. Since η Maxwell is close to η melt , we expect that our tidal heating model will produce maximum tidal dissipation near the base of the shell. In summary, the amount of tidal heating and the shell temperature profile depends strongly on the assumed melt viscosity.

Temperature structure
Conductive heat transport through a tidally heated shell defines its temperature profile, which controls its viscosity profile. Since tidal heating depends on the viscosity profile, we need to solve for tidal dissipation and temperature profile simultaneously. In a steady state, assuming that heat transport is dominated in the radial direction, the temperature and tidal dissipation satisfy the 1D conductive heat equation in spherical coordinates: where the source term h tide (r) is the surface averaged volumetric tidal dissipation rate found in Section 2.2 and thermal conductivity for pure water ice is given by κ(T ) = 0.4685 + 488.12/T in W m −1 K −1 (Hobbs 2010). Once the tidal dissipation term h tide (r) is known, Eq. 3 can be solved numerically as a boundary value problem with temperatures set at the surface and base of the shell. The temperature at the base is assumed to be T b = 273 K and the surface temperature T s is calculated as an equilibrium temperature using a fast rotator approximation, which gives T s = 92 K for Europa and T s = 59 K for Enceladus (see Appendix B). For the range of viscosity values η melt considered, Tobie et al. (2003) and Roberts & Nimmo (2008) show for Europa and Enceladus, respectively, that the tidal heating rate at the base is larger than the heat flux from radiogenic heating in the silicate mantle. However, the tidal heat production in the rocky mantles can be significant if the rocky mantle is loosely consolidated (Roberts 2015). In our modeling, however, we neglect mantle dissipation and the heat flux from the ocean for simplicity, considering only the tidal heat produced within the icy shell. Thus, the mantle is assumed to have an infinite viscosity and, thus, deforms only elastically.
To solve for a steady-state solution of the heat conduction equation, we converge the tidal dissipation and temperature profiles with the following iterative approach: 1. The solution to the boundary value problem of Eq. 3 excluding the tidal heating term (h tide (r) = 0) is set as an initial guess for the temperature profile.
2. The temperature profile obtained is used to find the viscosity profile with Eq. 2, and tidal dissipation h tide (r) is computed.
3. Then, the full Eq. 3 is used with h tide (r) and a new temperature profile is found as the solution to the boundary value problem.
Steps 2 and 3 are repeated until the temperature profile has converged or the temperature at any point within the shell has exceeded 273 K. Temperatures above the melting point in portions of the shell would cause different flow mechanisms to dominate due to partial melting, which would invalidate our ice rheology assumption (Eq. 2). Thus, we exclude such temperature profiles from the relaxation modeling.

Viscous relaxation model
Given a shell structure, we model the viscous relaxation of the icy shell as an incompressible fluid in a spherical shell with self-gravitation following the methodology of Hager & Clayton (1989). The icy shell is described by a sequence of layers of constant density and viscosity as prescribed by the viscosity profile. We have tested how the relaxation study results depend on the number of layers and found that 80 layers are needed to arrive at a converged solution. Incompressible Stokes flow along with Poisson's equation for gravitational potential are formulated as a system of six linear differential equations for a spherically symmetric shell (Hager & Clayton 1989). The state vector is given as radial functions of radial and poloidal velocities, radial normal stress, and poloidal shear stress as coefficients of vector spherical harmonics. Gravitational potential perturbation and its derivative are given as coefficients of scalar spherical harmonics. Given boundary conditions at the surface and base of the shell, the system of equations can be solved with a propagator matrix method (Gantmacher & Brenner 2005). Thus, our layered model can be expressed as a single linear system that relates the state vector at the ocean-ice interface to the state vector at the ice-outer space interface. A full description of the viscous relaxation model and the implementation of boundary conditions (Section 2.4.1) are given in Appendix C.

Boundary conditions
The two boundaries of the icy shell are the ice-outer space interface at the surface of the moon and the ice-ocean interface at the base of the shell. Under the formulation of Hager & Clayton (1989), the shapes of the interfaces are represented as loads applied at the artificial spherical boundaries r = R and r = R − D, for the mean radius R and shell thickness D. The amplitude of this topography is assumed to be small compared to its wavelength 2πR/ n(n + 1) ≈ 2πR/n for a spherical harmonic degree n. The loads generated by topography create perturbations in radial stress and gravity potential. The surface ice-outer space boundary is a material boundary. Thus, the radial velocity at the boundary directly affects the surface topography. Finally, free-slip boundary conditions are imposed at both boundaries, setting shear stress to zero. At the bottom ice-ocean boundary, two types of boundary conditions are considered as illustrated in Fig. 2. A bottom material boundary condition, which we will refer to as "bottom material boundary", assumes that no material crosses the water-ice interface through phase transitions. A bottom material boundary is similar to the boundary condition described by Hager & O'Connell (1979) for an isostatic rebound problem, and "constant load" described by Beuthe (2020). With this material condition, the radial component of the flow at the base of the icy shell will create or remove topography. The two interfaces will approach hydrostatic equilibrium over time. On the other hand, a bottom phase boundary condition, or "bottom phase boundary", allows for a flow across the ice-ocean interface. Under this boundary condition, described byČadek et al. (2019) to model the shell of Enceladus, the surface non-hydrostatic topography is maintained in a dynamic equilibrium by the flow at the base of the shell that is balanced by phase transitions (melting or freezing):n wheren is a unit vector normal to the interface, v is the velocity of ice flow, L is the latent heat of fusion for ice, ρ 1 is the density of the shell, and q 1 , q 2 are heat fluxes from the icy shell and ocean respectively. Our bottom phase boundary is the same as "constant shape" described by Beuthe (2020).
We solve for the Stokes flow using each of these boundary conditions with the goal of obtaining the asymptotic shape of the icy shell at spherical harmonic degree n for the calculation of gravity-topography admittance (Section 2.5). The shape of the icy shell at the surface and base are expanded in spherical harmonics at radii R and R − D, and their spherical harmonic coefficients h t and h b are referred to as shape coefficients or amplitudes, where the dependence on n is implied. We note that the term "shape", when used in describing our computations, is distinct from "topography", which refers to the elevation difference between shape and the equipotential surface. However, we use the term "topography" liberally to describe "shape" throughout the paper when describing features at the boundaries or topography support mechanisms. The shape of the ice-ocean boundary is described by the shape ratio w n = −h b /h t . This is different from the topographic ratio used inČadek et al. (2019), which describes the ratio of topographies. The negative sign in the definition of w n is a matter of convention and we use it to make w n positive for the conventional Airy isostasy case.
Solving for the Stokes flow under the described boundary conditions requires two different approaches. For the bottom material boundary, the shapes of the interfaces are functions of time and eventually relax to hydrostatic equilibrium. Applying the bottom material boundary, the propagator matrix solution of the Stokes flow equations (Hager & Clayton 1989) can be rearranged as a 2 × 2 system of ordinary differential equations for radial velocities dht dt and dh b dt (see Appendix C for details). The solution for the shape amplitudes T is given as a sum of two decaying exponentials describing the time evolution of the interfaces: where T are the eigenvectors of the system, and τ 1 = 1/γ 1 , τ 2 = 1/γ 2 are the characteristic decay times. We set τ 1 < τ 2 , so that the second term in Eq. 5 corresponds to the longer decay time. Two modes arise from Eq. 5, which are described by Hager & O'Connell (1979) as symmetric and antisymmetric modes of relaxation.
In the symmetric mode, the two interfaces move in the same direction. For example, if a mountain were placed on the surface and allowed to relax, a crustal root would initially form at the base of the icy shell. Conversely, in the antisymmetric mode, the two interfaces move in opposite directions. In this mode of relaxation, the height of the mountain and the crustal root will decrease with time, approaching hydrostatic equilibrium. We obtain the asymptotic state of the shell from the eigenvector [B s , B b ] T of Eq. 5 corresponding to the longer decay time τ 2 . The degree-dependent shape ratio is given by For the bottom phase boundary, there is no time-dependent solution of the interfaces. Instead, we solve for h t and h b as constants for the dynamic equilibrium between ice flow and phase transitions at the base. In this state, the radial velocity of the surface interface is set to zero to maintain the shape of the surface. However, the radial velocity at the base can be nonzero as long as the flow across the bottom interface is balanced by phase transitions satisfying Eq. 4. By setting surface shape h t to unity, the propagator matrix solution to the system of Hager & Clayton (1989) can be rearranged into a system of four equations with four unknowns: bottom shape h b , bottom radial velocity, bottom poloidal velocity, and surface poloidal velocity. Thus, the shape ratio is equal to w n = −h b . We note that with an observed value of h t , assuming that the shell is in this dynamic equilibrium, we can find bottom radial velocity, and thus, heat flux from Eq. 4. Our method of applying the bottom phase boundary is similar to the spectral method used byČadek et al. (2019).

Gravity-topography admittance
Finally, we compute the gravity-topography admittance using the shape amplitudes of viscously relaxed icy shells. We follow an approach of Ermakov et al. (2017) used to model the admittance of Ceres. Assuming shape amplitude is small relative to its wavelength, we use a mass-sheet approximation to get a linear relationship between gravity and shape. Summing the contributions to gravity from the shape of the surface and bottom boundaries of the shell, and dividing by surface shape, we find an expression for admittance Z n at degree n for shape ratio w n : whereρ is the mean density of the moon, ρ 1 is the density of the icy shell, ρ 2 is the ocean density, and ∆ρ = ρ 2 − ρ 1 is the ice-ocean density contrast. See Appendix D for derivation. Thus, with a degree-dependent shape ratio w n obtained from the viscous relaxation model, we find the admittance spectrum for a given shell structure.
We have two points of comparison for admittance. First, in the classical case of Airy isostatic compensation and assuming Cartesian geometry, the amplitude of the bottom shape is related to the amplitude at the surface by the ratio of the shell density to the density contrast, leading to w n = ρ 1 /∆ρ. Airy compensated values of admittance can inform us of topography supported by buoyancy forces acting on the bottom topography, or isostatic roots. The second limiting case is that of uncompensated surface topography, for which there is no corresponding bottom topography h b = 0. Therefore, w n = 0 for uncompensated surface topography. The admittance for the uncompensated topography provides an upper bound for admittance for a given density of a shell. Uncompensated values can indicate that topography is supported by stresses within the shell. We note, however, that Airy admittance does not necessarily provide a lower bound on admittance.

The effect of the bottom boundary condition
We begin by studying the effects of bottom material boundary and bottom phase boundary on admittance. A model of Europa with a 30-km shell is created for a range of viscosity profiles. Viscosity is assumed to decrease exponentially with depth. We vary the viscosity contrast between the top and the bottom of the shell (η s /η b ). Fig. 3 compares admittance spectra for both boundary conditions. We observe that for the uniform viscosity case (η s /η b = 1) shown in purple, the two boundary conditions produce admittance spectra that are nearly identical to the Airy isostasy case shown as the dotted curve. However, as the viscosity gradient steepens, the two boundary conditions start to diverge at higher spherical harmonic degrees.
The bottom phase boundary produces lower admittance values than for bottom material boundary. Lower admittance values correspond to larger values of shape ration w n . Thus, in this case, the shell has larger amplitude bottom topography relative to its surface topography compared to the bottom material boundary case. This is expected as the dynamic equilibrium between phase transitions and ice flow sustains more bottom topography compared to a material boundary. However, large amplitude topography at the surface could indicate that an unrealistic heat flux is needed to sustain melting and freezing. Bottom phase boundary can be justified if a sufficient heat source is present at the appropriate wavelength and magnitude to support the observed surface topography. This would be difficult to reconcile with tidal heating. The tidal potential is dominated by the degree 2 component. The higher-degree terms decrease quickly in amplitude. For example, for Europa, the degree 3 term is 1/430 the degree 2 term (Sabadini et al. 2016). Therefore, tidal heating could drive phase transitions only at low degrees, unless there are significant heterogeneity within the shell or in the heat flow from the ocean.

The effect of the viscosity gradient within the shell
Here, we explore the effect of the shell viscosity structure on admittance using exponential viscosity profiles of variable steepness. To illustrate the effect of the viscosity profile, Fig. 4 shows the velocity field comparing the ice flow within a uniform and gradient-viscosity shell for the bottom material and phase boundaries. The flow within the icy shell of Enceladus is shown for illustrative purposes, as the Enceladus' shell is thicker compared to its radius relative to Europa's shell. The flow pattern for the four shown cases is qualitatively similar for Europa.
For uniform-viscosity shells, the flow is quasi-uniform in amplitude throughout the shell beneath the topographic load. (Fig. 4A and B). On the other hand, a viscosity gradient in the shell causes a strong lateral flow concentrated at the low-viscosity base for both types of bottom boundary conditions ( Fig. 4D and C). For bottom material boundary, this lateral flow rapidly relaxes basal topography compared to the relaxation of the high-viscosity near-surface layers. For bottom phase boundary, there is a net freezing at the base of the shell associated with the topographic load at the pole. This freezing is balanced by upward radial flow such that the bottom interface remains stationary.
We observe that as the viscosity gradient steepens, admittance values approach high, uncompensated values at progressively lower degrees for bottom material boundary (Fig. 3A). Uncompensated values indicate that there is little basal topography relative to surface topography. Thus, surface topography is not supported by buoyancy as in the Airy isostasy case. Instead, viscous stresses within the shell provide the dominant topography support mechanism. For very large viscosity contrasts (the red bounded and unbounded curves for the bottom material boundary, Fig. 3A), the high relative viscosity at the surface impedes the relaxation of surface topography even at large scales while the bottom topography relaxes, leading to uncompensated admittance values at low degrees. For the bottom phase boundary, we find that admittance spectra are not very sensitive to the steepness of the viscosity gradient (Fig. 3B), especially at the . Admittance spectra of Europa with a 30 km icy shell. We show admittance spectra with both (A) bottom material and (B) bottom phase boundary conditions for a range of viscosity gradients. Viscosity is assumed to decrease exponentially with depth, with the viscosity contrast from the surface to base ηs/η b where η b = 10 14 Pa s is fixed. For the two high viscosity contrasts, admittance spectra for bounded viscosity η bound = 10 22 Pa s is shown. The uniform viscosity case (ηs/η b = 1) is plotted in purple. Admittance spectra assuming uncompensated topography and Airy isostasy are shown by the black dashed and dotted curves, respectively. Admittance values approach uncompensated values at lower spherical harmonic degrees with steepening viscosity gradient.
lowest degrees. This observation also holds for admittance spectra of Enceladus, and is consistent with the observation of similar shape amplitudes of the bottom interfaces seen in Fig. 4B and D.
We note that admittance spectra are only affected by the viscosity contrast η s /η b and not the actual values of η s or η b , which, instead, affect the timescales of relaxation, with higher viscosities increasing the time it takes to reach the asymptotic state of topography. The choice of the top viscosity bound affects the admittance spectra, with higher values for η bound increasing the viscosity contrast, causing admittance to approach uncompensated values at lower degrees.

Tidal heating influence on viscosity
We proceed from simple viscosity gradients to more realistic viscosity profiles by including a temperature-dependent viscosity and tidal heating. We incorporate the effect of tidal heating on a conductive temperature profile of the icy shell using a steady-state solution obtained from the iterative method described in Section 2.3. Fig. 5 shows the effect of including tidal heating on the temperature (Fig. 5A), viscosity (Fig. 5B), and volumetric tidal dissipation rate (Fig. 5C) for two choices of η melt . As mentioned in Section 2.2, tidal dissipation is maximized if the forcing period is near the Maxwell time of the material being deformed. The viscosity for which this occurs for Europa (η Maxwell ≈ 1.6 × 10 14 Pa s) is shown as the dotted line in Fig. 5B. The viscosity profiles of the 30-km shells modeled with two values of η melt cross η Maxwell near the base of the shell. For η melt = 10 13 Pa s, this intersection occurs at ≈78% of the shell's depth, which corresponds to a maximum in tidal dissipation as expected. For the case with η melt = 10 14 Pa s, maximum dissipation occurs closer to the base of the shell. This causes tidal dissipation to be lower throughout the shell, which is seen in the temperature profile in Fig. 5A, where η melt = 10 14 Pa s causes the temperature to be closer to the no tidal heating case. Although not shown in Fig. 5, η melt = 10 15 Pa s further decreases the influence of tidal heating, and the temperature profile becomes nearly identical to the case without tidal heating. In summary, the choice of η melt , and thus the assumption of the grain size of ice, has a large effect on tidal heating and, therefore, the temperature structure of the shell.

Shell thickness
Finally, we use our complete model with temperature-dependent viscosity and tidal dissipation to study the sensitivity of admittance to the shell thickness. Predictions of Europa's shell thickness vary from local estimates of several km to more than 30 km as a global average (Billings & Kattenhorn 2005). We consider mean thicknesses between D = 8 and 35 km. Depending on the viscosity at the base η melt , the amount of tidal heating for thick shells causes temperature profiles to not converge with our iterative method described in Section 2.3. For η melt = 10 14 Pa s, a temperature profile did not converge above a shell thickness of 32 km due to temperatures exceeding the melting point in our iterative method.
The shell thickness has a strong effect on the admittance spectrum, which is illustrated in Fig. 6. The bottom material boundary condition is used in Fig. 6A and the bottom phase boundary condition is used in Fig. 6B. We observed that increasing shell thickness lowers viscosity at relative depths (Fig. 6C) as thicker shells generate more heat from tidal dissipation. To understand the behavior of admittance and topography support mechanisms at different wavelengths, we apply Jeffreys' theorem. It states that the minimum stress difference required to support a surface load is ≈ 1/3 of the load and is concentrated in a region with dimensions comparable to the width of the load (Melosh 2011). Thus, a surface topographic load samples the shell beneath it to a depth comparable to its wavelength. A short-wavelength load feels only the near-surface viscosity value. A load with wavelength comparable to the shell thickness samples deeper into the shell and is sensitive to its viscosity profile. For a load with a wavelength much longer than the shell thickness, the shell acts essentially as a membrane. Thus, such a load is not sensitive to the viscosity profile.
As shell thickness increases, admittance for bottom material boundary (Fig. 6A) approaches uncompensated values at progressively lower degrees or longer wavelengths. For thicker shells, it is difficult to propagate buoyancy-produced stresses from the base of the shell to support surface topography, causing topography at progressively longer wavelengths to be supported by stresses within the shell. Additionally, wavelengths much longer than the thickness of the shell do not sense the viscosity variation throughout the shell, thus yielding an Airy-like admittance. This is seen for the 8-km shell at low degrees (n < 15) and at the lowest degrees for a 16-km shell (Fig. 6A). For shorter wavelengths, admittance becomes sensitive to the steep viscosity profile resulting from a tidally heated conductive shell and approaches uncompensated values.
For bottom phase boundary (Fig. 6B), admittances are generally closer to Airy isostasy admittance values. As the shell thickness increases, the degree at which admittance deviates from the Airy isostasy becomes progressively lower, however at low degrees (n < 15), values of admittance are similar to that of Airy isostasy. For thicker shells, the divergence between boundary conditions occurs at lower degrees. In a scenario where the shell is thick, measurements showing a dip in admittance at a high degree-characteristic of bottom phase boundary-could indicate the presence of a heat source pattern sufficient to maintain the dynamic equilibrium of bottom phase boundary.

Topography support mechanism
Admittance can provide insight about topography support mechanism. High, uncompensated values of admittance can indicate that surface topography that is supported by viscous stresses within the shell. Low admittance could indicate that topography is supported by buoyancy, as predicted from the Airy isostasy model. Alternatively, low admittance could indicate that surface topography is supported by flow arising from phase transitions at the base of the shell. In this case, the heat flux pattern controls the rate of basal melting and freezing. Further study on the distribution of basal heat flux is needed to understand the wavelengths at which bottom phase boundary could be applicable. If tidal heating in the shell is small, for example due to high values of η melt , the influence of radiogenic and tidal heating from the mantle and ocean (Roberts 2015;Rekier et al. 2019; Rovira-Navarro et al. 2019) may become more important. Ocean circulation patterns resulting from mantle heat sources and salinity gradients could affect the distribution of heat at the base of the shell (Kang et al. 2020(Kang et al. , 2021.

Ice rheology
Temperature-dependent viscosity for a tidally heated conductive shell yields a steep viscosity gradient within the shell. The depth at which Maxwell time is crossed influences the distribution of tidal dissipation, which, in turn, determines the extent of the high-temperature, low-viscosity region at the base of the shell. On Europa, tidal dissipation is concentrated at the base of the shell, where viscosity is low and is close to η Maxwell . We find that the low viscosity at the base causes rapid relaxation of bottom topography by lateral flow, leading to uncompensated surface topography. Thus, the Airy isostasy model is likely not suitable for interpreting admittance except at the longest wavelengths, which are not sensitive to the steep viscosity gradient, or where the bottom phase boundary is applicable.
Ice deformation mechanisms are difficult to study in the laboratory, and different non-Newtonian flow mechanisms such as grain boundary sliding may be more appropriate when considering solid-state convection (Barr & Showman 2009). Although our analytical viscous relaxation model does not allow for convection, we expect that convection at the base of the shell would effectively relax shorter wavelength bottom topography resulting in uncompensated values of admittance.
In addition to the flow mechanism, the grain size of ice has a large effect on the viscosity profile, tidal dissipation, and thus the temperature structure of the shell. The viscosity of ice at the base of the shell at the melting temperature depends on this poorly constrained grain size of ice. For Europa and Enceladus, the viscosity at which Maxwell time is crossed, η Maxwell , lies within the range of melt viscosities η melt = 10 13 − 10 15 Pa s typically considered. We find that tidal heating has little influence on temperature structure when η melt > η Maxwell . If η melt < η Maxwell , the influence of η melt on temperature profile is amplified depending on the extent of the low viscosity region and the relative depth at which η Maxwell is crossed. The melting temperature T melt changes with pressure and affects the value of η melt and the temperature gradient. However, changes in grain size likely have a larger effect on η melt , thus we can approximate different melting temperatures with the range of η melt considered.
The diffusion creep temperature-dependent viscosity predicts high viscosities at the surfaces of Europa and Enceladus, implying very slow relaxation of surface topography leading to uncompensated values of admittance. Properties of shallow sub-surface ice such as fracturing, porosity, and impurities affect deformation mechanisms, viscosity, and thermal conductivity. An insulating layer, considered by Bland et al. (2012) in their study of crater relaxation on Enceladus, could raise the shallow sub-surface temperature and, thus, lower viscosities near the surface. In addition, isolated periods of warming may have occurred in the past (Bland et al. 2012). Such warming episodes could have facilitated viscous relaxation affecting admittance. Thus, measuring admittance and comparing the inferred viscosity profile to predictions from the crater relaxation study could validate the notion of past heating episodes on Enceladus. We now apply our methods to a model of Enceladus and compare it to Europa. We assume a shell thickness of 21 km for Enceladus taking the central value from Hemingway & Mittal (2019), and compare it to Europa with a 30 km shell using melt viscosity η melt = 10 14 Pa s for both ocean worlds. Enceladus has a thicker shell in proportion to its radius compared to Europa. We show a comparison of admittance spectra in Fig. 7. Since the wavelength is inversely proportional to degree, the same degree on Enceladus corresponds to a shorter wavelength than on Europa. Thus, we expect admittance spectra of bottom material boundary for Enceladus to approach uncompensated values at a lower spherical harmonic degree than for Europa due to shorter wavelengths at low degrees for Enceladus sampling the viscosity gradient of the shell. Comparing admittance of bottom phase boundary, we see that admittance for Enceladus deviates from the Airy isostasy case at low degrees more than for Europa, perhaps making it easier to distinguish topography support mechanism from future observations. For tidal heating, the key differences between the two moons are orbital frequency and surface temperature. The higher orbital frequency of Enceladus causes η Maxwell = 6.2 × 10 13 Pa s to be lower than that of Europa. Thus, in Fig. 7 with η melt = 10 14 Pa s, η melt > η Maxwell for Enceladus while η melt < η Maxwell for Europa. A lower temperature of Enceladus leads to a steeper viscosity gradient within Enceladus' shell. This reduces the thickness of the layer where viscosity close to η Maxwell , thus causing an overall reduction of the shell tidal heating. As a result, tidal heating does not affect the shell temperature profile on Enceladus as strongly as it does for Europa.

Comparing Europa to Enceladus
For the bottom phase boundary, Fig. 7 shows lower values of admittance compared to the bottom material boundary at low degrees for our Enceladus model. In general, Airy-compensated values of admittance only hold for low spherical harmonic degrees, or long wavelength is much greater than the shell thickness. Thus, we reach the same conclusion aš Cadek et al. (2019) that the Airy isostasy model can be applied for long wavelengths but may be inaccurate at short wavelengths for the bottom phase boundary. 4.4. Sensitivity of future gravity and shape data to the icy shell structure In this subsection, we characterize what can be achieved with future data at Europa and Enceladus by comparing various shell structure endmembers to expected mission performance. To access the accuracy of admittance recovery two items are needed. First, one needs to estimate the accuracy of gravity coefficients determination, which can be done either with a simplified approach of Bills & Ermakov (2019) or with a full mission simulation using covariance analysis. Second, a global shape model or, at least, an estimate of the shape power spectrum is required. A global shape model is available for Enceladus (Tajeddine et al. 2017) but is not currently available for Europa to our knowledge. Thus, we can only approximately judge about admittance recovery at Europa based on the expected resolution of the gravity field. Park et al. (2011) conducted a covariance analysis study simulating Europa Clipper flyby mission and found that Europa's gravity field may be resolved up to degree 10. A transition between low (Airy compensation-like) admittance values to higher, uncompensated values marks the shift of the dominant topography support mechanism from buoyancy to viscous stresses within the shell. Thus, we can pose a question: under what conditions does this transition occur at n < 10? We find the transition is captured in admittance spectra at n < 10 for shell thicknesses over ≈24 km for η melt = 10 14 Pa s (Fig. 6A). If a lower melt viscosity is assumed (η melt = 10 13 Pa s, not shown in Fig. 6A), the tidal heating is stronger within the shell and the admittance transition is captured for shell thicknesses over ≈15 km. Thus, we conclude, that admittance is more sensitive to warm and thick shells. Extending the capability of capturing the transition of topography support mechanism to higher degrees will improve the sensitivity of admittance to thinner shells. Ermakov et al. (2021) presented a covariance analysis for Enceladus orbiter missions. Ermakov et al. (2021) studied the recovery of Enceladus' gravity field and tides simulating one month of continuous radio-tracking for a single orbiter with radio-tracking to the Earth and a GRAIL-like dual spacecraft with inter-satellite tracking. The ranging accuracy was assumed 10 −7 km s −1 for the single spacecraft case and 10 −9 km s −1 for the dual spacecraft case. These accuracies are typical for X-band and Ka-band ranging, respectively. We used the covariance analysis by Ermakov et al. (2021) and Enceladus shape model by Tajeddine et al. (2017) to estimate the performance of two orbiter mission configurations in recovering admittance. The gravity error RMS spectra for the two orbiter configurations are found from the diagonal elements (i.e., variances) of the gravity covariance matrix σ 2 Cnm and σ 2 Snm as: The gravity error RMS spectra for the two orbiter configurations are shown in Fig. 8B and correspond to blue and yellow curves in Fig. 6 in Ermakov et al. (2021). In order to match Ermakov et al. (2021), we present the error in radial gravitational acceleration coefficient M gg n g(n + 1), where g is surface gravity. For comparison, we show gravity error RMS from two currently available gravity field models by Iess et al. (2014) for degrees 2 and 3. The variance of admittance at degree n is given by: where D n is a vector of partial derivatives of degree-n admittance with respect to gravity coefficients: Here, V tt n is the shape variance spectrum (see Eq. D.3 in Appendix D for definition),h n = [Ā n0 ,Ā n1 ,B n1 , ...,Ā nn ,B nn ] is a vector of normalized shape coefficients from the shape model of Tajeddine et al. (2017), and C n is the degree n submatrix of the covariance matrix from Ermakov et al. (2021). Note that Tajeddine et al. (2017) provided unnormalized shape coefficients. Shape is assumed to be exact, thus the error comes solely from gravity. The admittance error σ Zn is shown for the single and dual spacecraft configurations in Fig. 8A. Excluding covariances by looking at the diagonal elements of C n , we find the correlations between recovered coefficients contribute little to the admittance error for both orbiter configurations. Similarly, admittance error is shown for the gravity error for the two gravity field models from Iess et al. (2014).
The ability to distinguish between different endmembers of the shell structure using admittance can be used to place a requirement on the gravity error. A gravity error RMS spectum is estimated from an admittance error requirement σ Zn by: where we assume gravity error uniform across all orders for each degree and there are no correlations between coefficients of different orders. A range of admittance errors from 1 to 30 mGal km −1 are shown as gravity error RMS spectra in Fig. 8B. An admittance error requirement is satisfied if the gravity error RMS of a given spacecraft mission configuration is below that for a desired admittance error represented by colored curves Fig. 8B). For the model of Enceladus shown in Fig. 7, the difference in admittance due to the boundary conditions are from ≈20 to 30 mGal km −1 at degrees up Figure 8. Comparison of different admittance error thresholds to simulated observational error from two Enceladus orbiter missions from Ermakov et al. (2021) and current gravity error from Iess et al. (2014). Admittance error (A) and gravity error RMS spectra (B) are shown for a single spacecraft with X-band tracking with 10 −7 km s −1 accuracy and a GRAIL-like dual spacecraft with inter-satellite Ka-band tracking with 10 −9 km s −1 accuracy (corresponding to Fig. 6 in Ermakov et al. (2021)). Current observational gravity error and admittance error from two gravity field models (SOL1 and SOL2) are shown for degrees 2 and 3 (Iess et al. 2014). A set of gravity performance curves is shown for a given admittance error threshold with factor g(n + 1) with gravitational acceleration of Enceladus g for units of mGal (B).
to 10. As can be seen in Fig. 8B, both mission configurations studied in (Ermakov et al. 2021) would yield admittance accuracy smaller than the expected difference between the two boundary conditions. Thus, such mission configurations would provide sufficient accuracy to distinguish between two types of behaviour of the shell-ocean interface. The use of admittance may augment ice-penetrating radar in constraining shell thickness. One of the goals of Europa Clipper's REASON (Radar for Europa Assessment and Sounding: Ocean to Near-surface) instrument is to detect the ice-ocean interface. However, the attenuation of radio waves within ice limits the detection of the ice-ocean interface for the case of a thick shell. Kalousová et al. (2017) find that direct detection of Europa's ocean with REASON may be possible up to 15 km for a conductive shell or in an area of cold downwelling within a convective shell. In the case that the icy shell of Europa is thick, admittance, which we find is more sensitive to thicker shells, may augment radar in the constraining of shell thickness.

Sensitivity of admittance to tidal dissipation
A measurement of the tidal phase lag, or equivalently, the imaginary part of the tidal Love number k 2 can provide a constraint on the total dissipation within the satellite (Peale et al. 1979). However, this measurement is challenging and so far has been achieved only for the Earth (Ray et al. 2001), the Moon (Williams et al. 2014) and Mars (Bills et al. 2005), which have abundant data sets. Early future data for ocean worlds will likely be more limited. Park et al. (2011) estimated that Europa Clipper will enable determination of k 2 of Europa with an uncertainty of 0.009. A more recent study by Verma & Margot (2018) is more pessimistic and predicts that Europa's k 2 can be determined with an uncertainty of ≈ 0.05 − 0.06 depending on the Europa Clipper trajectory and the selection of the groundbased radio-tracking assets. We computed the real and imaginary parts of Europa's Love numbers (See Appendix A for details) and explored their dependence on the shell thickness and viscosity profile. Fig. 9 shows how Re(k 2 ) and Im(k 2 ) depend on the shell thickness and melt viscosity η melt for a model of Europa with viscosity bound η bound = 10 24 Pa s. We observe that Re(k 2 ) depends primarily on the shell thickness. On the other hand Im(k 2 ), and therefore total dissipation, varies strongly depending on both shell thickness and η melt . In addition, Fig. 9 shows gravity-topography admittance at degree 3 and 10 as filled contours. It can be seen that the admittance contours are, in general, not parallel to the contours of Re(k 2 ) or Im(k 2 ). Thus, measuring admittance would provide non-degenerate information about the dissipation within the icy shell.
We note that the values of Im(k 2 ) are smaller than the estimated uncertainty of k 2 recovery (Park et al. 2011;Verma & Margot 2018) assuming Im(k 2 ) is measured to the same precision as Re(k 2 ). Thus, directly constraining total tidal dissipation within Europa will be challenging from the Europa Clipper data. Park et al. (2011) showed that the gravity field of Europa can be estimated to degree 10. More recently, Verma & Margot (2018) showed that only degree 3 and 4 of the gravity field can be estimated. The global shape can be estimated from fitting limb profiles or from building a geodetic network of reference points, which has been done previously for Enceladus using the Cassini flyby data Combinations of thickness and η melt where a temperature profile did not converge according to our criteria (Section 2.3) are shown in grey. (Nimmo et al. 2011;Tajeddine et al. 2017). Thus, since measuring the gravity field, the shape and, therefore, gravitytopography admittance, at low degrees is less demanding than measuring Im(k 2 ), admittance measurements would likely allow constraining tidal dissipation within Europa's shell prior to the Im(k 2 ) measurement becoming available.
Finally, if tidal dissipation is estimated from both Im(k 2 ) and admittance, the potential disagreement between these estimates could indicate either a violation of our modeling assumptions such as: the choice of η bound and the rheology of ice near the cold surface; the choice of the bottom boundary condition; the assumptions of uniformity of tidal heating distribution and conductive heat transport; or that significant heating occurs in the rocky mantle or the ocean.

CONCLUSIONS
We simulated the viscous relaxation of the icy shells of ocean worlds to test the sensitivity of gravity-topography admittance to the structure of the shell. We find that admittance is sensitive to the topography support mechanism. Our models show that the behavior of the bottom interface-whether it is primarily a material or phase interfacestrongly influences the admittance spectrum. The base of the icy shell may behave as a phase boundary at low degrees where tidal heating can drive melting and freezing. At higher degrees, a lack of a heat source would render a material boundary condition more applicable. If the base is treated as a material boundary, low viscosities cause the relaxation of basal topography by lateral flow faster than the relaxation of higher-viscosity surface topography, leading to uncompensated topography. This makes the Airy isostasy model not suitable for interpreting admittance at degrees for which bottom material boundary is applicable except at the longest scales, where topographic loads with wavelengths much larger than shell thickness do not sense the viscosity variation throughout the shell.
A transition between topography support mechanisms may be resolved with future measurements, allowing us to probe the shell structure of ocean worlds. The viscosity structure of the shell is controlled by the rheology of ice. In particular, the poorly constrained grain size of ice influences the tidal heat production and the extent of the low viscosity region at the base of the shell affecting the temperature profile. At low degrees, admittance spectra are more sensitive to temperature and viscosity structures for thick shells with high tidal dissipation. Measurements of higher-degree admittance would allow for studying a wider range of potential shell structures. The higher sensitivity of admittance to the thickness of thicker shells may augment constraining the shell thickness by ice-penetrating radar, which is more sensitive to thinner shells. Thus, a combination of radar and gravity measurements in future ocean world missions may improve the robustness of the measurement strategy.
We also find that admittance measurement can be used to constrain the tidal dissipation within the icy shell. Such a measurement would be complementary to a demanding measurement of the imaginary part of the tidal Love number, which is proportional to the total tidal dissipation. Measuring admittance in addition to Im(k 2 ) can be used to separate the shell heating from the mantle and ocean heating. Knowledge of heating distribution within ocean worlds would allow for a better understanding of the heat budget and stability of their oceans, which is critical for their long-term habitability.

ACKNOWLEDGMENTS
This work sprung form conversations with Roger Fu. The authors thank Ryan Park for providing covariance matrices for gravity measurements from simulations. The authors also thank Peter James, Doug Hemingway, Mikael Beuthe, and Shunichi Kamata for their helpful insights. RA was supported by the National Science Foundation-Department of Energy (DOE) partnership for plasma science and engineering (grant DE-SC0016248). AE received support through a Cassini Data Analysis grant (NNX16AI43G). BM acknowledges support from the U.S. Department of Energy and National Nuclear Security Administration (grant DE-NA0003842).

A. Tidal heating
We compute the volumetric tidal heat production rate using a spherically symmetric internal structure model and a Maxwell model of viscoelasticity. Our modeling follows the steps of Takeuchi & Saito (1972), who define six radial functions y n i (r) to describe tidal flow field. These functions are radial multipliers of the corresponding fields expanded in spherical harmonics. Index n refers to a spherical harmonic degree. We restrict our analysis to degree 2 tides and, for now on, we will omit the dependence on n. The functions y i describe radial and tangential displacements (y 1 and y 3 , respectively) and radial and tangential stresses (y 2 and y 4 ) expanded in vector spherical harmonics. Gravitational potential (y 5 ) is expanded in scalar spherical harmonics. y 6 contains the radial derivative of gravitational potential and is formulated in the following way by Takeuchi & Saito (1972) to simplify the surface boundary condition: where ρ is the density of the layer. Note that definition for y 6 (r) is different from another commonly used notation of Alterman et al. (1959). These six radial functions are found by solving a system of linear differential equations within layers of constant density and viscoelastic parameters: where A tidal is a 6×6 matrix given by Eq. 82 in Takeuchi & Saito (1972). The surface boundary condition at r = R for surface radius R are given by vanishing stresses (y 2 (R) = y 4 (R) = 0) and y 6 (R) = (2n + 1)/R describing potential and its derivative continuity. Across the liquid-solid interfaces, which occurs at the ice-ocean and ocean-mantle boundaries for our three-layer model, additional boundary conditions are required. At the center of the planet r = 0, displacements are zero (y 1 (0) = y 3 (0) = 0) and gravitational potential is also zero (y 5 (0) = 0). The system of differential equations becomes singular at r = 0. An analytical solution at the surface of a small homogeneous sphere is used as a starting solution of radial functions y i (r) (see Takeuchi & Saito (1972) and Martens (2016) for more detail). This starting solution at r 0 R, consisting of three linearly independent solutions, is then propagated through successive layers by solving Eq. A.2 with an eighth order Runge-Kutta numerical integrator scheme using variable normalization from Martens (2016). Finally, the three independent solutions at the surface are linearly weighted to satisfy the surface boundary condition and combined to form the six radial functions y i (r). Complex-valued tidal Love number k 2 is found as y 5 (R) − 1.
To find the radial distribution of tidal heating, we follow Tobie et al. (2005). Tidal heating is driven by the periodic tidal potential. The degree 2 tidal potential for a tidally-locked satellite on an eccentric orbit, computed to the first order in eccentricity, and evaluated on the surface of the satellite is given by Moore & Schubert (2000): P 2,0 (cos θ) cos ωt + 1 4 P 2,2 (cos θ) (3 cos ωt cos 2φ + 4 sin ωt sin 2φ) , where ω is the orbital (and rotational) frequency, e is the eccentricity, θ is colatitude, φ is longitude, t is time, P 2,0 and P 2,2 are unnormalized associated Legendre functions. Our internal structure model consists of three main layers. The solid icy shell and mantle are described by the elastic shear modulus µ, Poisson's ratio ν, and viscosity η, while the ocean is described by bulk modulus κ following the values given in Table 1. The viscosity profile in the shell is broken up into 80 constant viscosity layers. The viscoelastic Maxwell rheology moduli used in solving for y i (r) are given by complex-valued Lamê's first parameterλ and complex shear modulusμ in terms of their real-valued counterparts along with the bulk modulus κ and viscosity η: We then use the radial functions y i (r) to find the sensitivity parameter H µ introduced by Tobie et al. (2005) to describe the radial sensitivity to the shear modulus (Eq. 33 of Tobie et al. (2005)), as it enables computing the radial distribution of the volumetric dissipation rate. A simplified version of H µ is given by Eq. 25 in Beuthe (2013)). Finally, the radial profile of surface averaged volumetric tidal dissipation rate averaged over moon's orbit is computed with Eq. 37 of Tobie et al. (2005) accounting for a sign correction made by Beuthe (2013): h tide (r) = 21 10 ω 5 R 4 e 2 r 2 H µ Im(μ). (A.6) Our tidal dissipation code was benchmarked by comparing our y i (r) profiles with values from Fig. C1 and Fig. C2 of Kamata et al. (2015).

B. Surface temperature
The surface temperature of our spherically symmetric moon is given by the uniform equilibrium temperature using a fast rotator approximation: where F is solar irradiance, A is the Bond albedo, and σ is the Stefan-Boltzmann constant. For Europa, we use solar irradiance at Jupiter of 50.26 W m −2 and Bond albedo 0.68 (Grundy et al. 2007) and find surface temperature T s ≈ 92 K. For Enceladus, we use solar irradiance at Saturn 14.82 W m −2 and a Bond albedo of 0.81 (Spencer et al. 2006) and find surface temperature T s ≈ 59 K.

C. Stokes flow with self-gravitation
Viscous relaxation of the icy shell is modeled by solving for the Stokes flow of an incompressible fluid in a spherical shell with self-gravitation. Similarly, to the tidal heating computation, a spherically symmetric icy shell is broken up into layers of constant density ρ and viscosity η, where the viscosity is set by the temperature profile found in Section 2.3. The governing equations are: where v is the velocity, V is the potential, p is the pressure, G is the universal gravitational constant and F is the external force, equal to gravity in our case. To solve for the Stokes flow, we follow Hager & Clayton (1989), who define six radial functions y n i (r) where n is spherical harmonic degree (Eq. 4.11 -4.18). These functions are not to be confused with the y functions discussed above for the tidal heating problem. Radial and poloidal velocities (y n 1 and y n 2 ), radial normal stress (y n 3 ), and poloidal shear stress (y n 4 ) are the coefficients of the corresponding vector fields expanded in vector spherical harmonics. Gravitational potential perturbation (y n 5 ) and its radial derivative (y n 6 ) are coefficients of scalar spherical harmonics. The decoupling of the system of equations for the velocities and stresses from that of the potential is achieved by introducing a new set of variables u n i and v n i . Again, we will drop the dependence on degree n for convenience. The components of the state vectors u and v are given in terms of the y i variables as: where ρ is local density, ρ 0 is a reference density, and η 0 is a reference viscosity chosen close to η bound . The radial derivatives for u i and v i can be written in matrix form: where ν = log(r/R), R is the mean radius of the moon, A Stokes is a 4 × 4 matrix and B Stokes is a 2 × 2 matrix given by Eq. 4.33 and 4.34 in Hager & Clayton (1989). These matrices depend only on viscosity η within a layer and spherical harmonic degree n. We note that the same system of equations C.1 is solved by Hager & O'Connell (1979) as a nearly identical coupled 6 × 6 system. With a choice of boundary conditions at the surface and base of the icy shell, a solution can be propagated from the base to the surface analytically using the propagator matrix method (Eq. 4.43 for the u-system and Eq. 4.44 for the v-system in Hager & Clayton 1989).
A free-slip boundary condition is imposed at both boundaries, setting poloidal shear stresses and, correspondingly, u 4t and u 4b to zero. To account for self-gravity in the boundary conditions for u 3 , we need to solve the second equation of Eqs. C.3 for v 1 at the surface and the base of the shell. We solve the v-system given the shape amplitudes of the top and bottom interface: h t and h b , respectively. We get the following expressions for v 1t and v 1b : where r b = R − D for shell thickness D; ρ 1 and ρ 2 are the densities of the icy shell and the ocean, respectively, with a density contrast of ∆ρ = ρ 2 − ρ 1 . Using the solution for v 1t and v 1b , we find the boundary conditions for u 3 at the surface (u 3t ) and base (u 3b ) using Eq. 4.53 and 4.54 in Hager & Clayton (1989). The boundary conditions for u 3t and u 3b are given as functions of shape amplitudes h t and h b as: where g 1 and g 2 are the gravitational accelerations at the top and bottom boundary, respectively. We can write these two equations in a matrix form: The propagator matrix solution for the u-system in Eqs. C.3 is: where u t and u b are the state vectors at the surface and the base, respectively. P is the propagator matrix, which is found by a product of layer-wise propagator matrices for each layer of uniform viscosity between bottom (k-th) boundary r b = r k = R − D, and the surface r 1 = R.
p 1,1 p 1,2 p 1,3 p 1,4 p 2,1 p 2,2 p 2,3 p 2,4 p 3,1 p 3,2 p 3,3 p 3,4 p 4,1 p 4,2 p 4,3 p 4,4      = P 1 (r 1 , r 2 ) · P 2 (r 2 , r 3 ) · ... · P k−1 (r k−1 , r k ) (C.8) The propagator matrix for the i-th layer is given by the following matrix exponential: For the bottom material boundary, we solve for a time-evolving system illustrated in Fig. 2. Thus, the radial flow velocities u 1t and u 1b are unknowns. We rearrange Eq. C.7 and use the boundary condition for u 3 from Eq. C.5 to express a system of linear differential equations for the flow velocities u 1 and u 2 at the top and the base of the shell: are the elements of the u 3 matrix from Eq. C.6. From Eq. C.10, we get the time evolution of h as a sum of two decaying exponentials, resulting in Eq. 5. The degree-dependent shape ratio w n = −h b /h t needed to compute admittance converges to the ratio of the elements of the eigen-vector corresponding to the longer relaxation time.
The bottom phase boundary condition is different in that we solve for a dynamic equilibrium illustrated in Fig. 2. Thus, the topography amplitude vector h is constant in time. This gives us the condition that radial velocity at the surface (u 1t ) is zero. However, although topography at the base h b is constant, radial velocity (u 1b ) can be nonzero because the flow across the boundary is assumed to be balanced by phase transitions (melting or freezing). Similar to bottom material boundary, u 3t and u 3b are taken as functions of h t and h b using Eq. C.6 and a free-slip boundary condition gives u 4t = u 4b = 0. We again rearrange the terms in the propagator matrix solution (Eq. C.7), expressing a vector of unknowns in terms of the known quantities:      p 1,3 B b p 1,1 0 p 1,2 p 2,3 B b p 2,1 −1 p 2,2 p 3,3 B b − B t p 3,1 0 p 3,2 p 4,3 B b p 4,1 0 p 4,2 We set h t = 1 to get four equations with four unknowns and solve for h b , u 1b , u 2t , and u 2b . Finally, we find the shape ratio w n = −h b /h t = −h b , which we use in computing the gravity-topography admittance.

D. Gravity-topography admittance
In this section, we present how gravity-topography admittance is 1) computed from the shape and gravity data, and 2) modeled using the degree-dependent shape ratio w n found in the previous subsection. Our approach is similar to Ermakov et al. (2017) in their study of Ceres' admittance. The measured shape of the icy moon can be expressed as a spherical harmonic expansion (Wieczorek 2007): r(θ, φ) = R 1 + ∞ n=1 n m=0 (Ā nm cos(mφ) +B nm sin(mφ))P nm (cos(θ)) , (D.1) where θ is colatitude, φ is longitude, R is the mean radius,Ā nm andB nm are normalized coefficients of the spherical harmonic expansion, andP nm are normalized associated Legendre functions for spherical harmonic degree n and order m. Similarly, the gravitational potential is also expanded in spherical harmonics: U (r, θ, φ) = GM r 1 + ∞ n=2 n m=0 R 0 r n (C nm cos(mφ) +S nm sin(mφ))P nm (cos(θ)) , where R 0 is the reference radius, r is the observation radius, andC nm andS nm are normalized spherical harmonic coefficients of the gravitational potential. Note that the reference radius is not necessarily equal to the mean radius R. The variance spectrum of the shape V tt n , gravity V gg n and the cross-variance spectrum of gravity and shape V gt n are found from the spherical harmonic coefficients in Eq. D.1 and D.2 as: We are interested in finding the gravity-topography admittance between radial gravitational acceleration and the shape. Thus, we need to differentiate Eq. D.2 with respect to r, which results in simply multiplying each term in Eq. D.2 by (n + 1)/r. Gravity-topography admittance Z n in terms of gravity-shape cross variance and shape variance is given as: This expression is known to produce a more accurate admittance estimate in the case when gravity data is more noisy than the shape data (McKenzie 1994), which is typically the case for ocean worlds. It is fully equivalent to: Z n = V gg n V tt n R gt n · GM R 3 (n + 1), (D.5) where R gt n = V gt n / V gg n V tt n is the gravity-shape correlation coefficient. We note that in the viscous relaxation models presented in this paper, the gravity and shape are always "in phase". Thus, the correlation is either 1 or -1. Since our model concerns small-amplitude topography, it is spectrally pure. That is the topography harmonic of degree n and order m produces the gravity perturbation only of the same degree and order. Thus, Eq. D.5 can be simply written in terms of shape and gravity coefficients: where σ = C nm ,S nm and h = Ā nm ,B nm are degree-n coefficients of gravity and shape, respectively. In order to compute the model admittance, we need to find the relationship between the gravity and shape spherical harmonic coefficients. Assuming a homogeneous spherical body with topography small in amplitude relative to its wavelength (a "mass-sheet" approximation), we can use the first-order term of Eq. 10 in Wieczorek & Phillips (1998) to get the linear relationship between gravity and shape coefficients: where R vol is the radius of the volume-equivalent sphere. For a body consisting of multiple layers, topography at the interfaces associated with density changes generates perturbation in gravity. We only consider topography at the surface and the base of the icy shell assuming the ocean floor has no topography. Thus, the total gravity includes contributions only from the upper two interfaces. The gravity contribution of each layer needs to be weighted by fractional mass of that layer (ρ i /ρ)(r i /R) 3 , whereρ is the mean density of the body. Since we seek the gravity coefficients referenced to the mean radius of the body (R 0 = R in Eq. D.2), the contributions to gravity from i-th layer must be upward propagated by (r i /R) n to find the contribution to gravity at the surface. Finally, we reference the topography coefficients to the mean radius R of the body in Eq. D.1. Thus, the topography at the base of the shell h b needs to be upscaled by R/(R − D). In summary, we find the total gravity coefficient at degree n as: σ = 3 2n + 1 · h t · ρ 1 ρ fractional mass In this derived expression, we ignored the factor (R/R vol ) 3 , since the difference between R vol and R is small for a small-amplitude topography, which is one of our assumptions anyway. Plugging in h b = −w n h t , dividing by the surface shape coefficients h t and multiplying by GM R −3 (n + 1) to find the ratio of radial gravitation acceleration to shape coefficients (in mGal km −1 ), the gravity-topography admittance simplifies to: The shape ratio w n is found given a viscosity profile and boundary conditions described in Appendix C. Plugging in w n = ρ 1 /∆ρ into Eq. D.9, we recover admittance for the standard Cartesian Airy compensation model: (D.10) If w n = 0, we recover admittance for an uncompensated topography: Z n = GM R 3 · 3(n + 1) (2n + 1) · ρ 1 ρ . (D.11)