Tidal Heating at Europa Using the Multifrequency Analysis of Tidal Heating Toolkit

We evaluate the thermal response of Europa’s ice shell to the gravity perturbation of Jupiter using a comprehensive toolkit (dubbed Multifrequency Analysis of Tidal Heating, MATH) that evaluates tidal heating and heat flow in planetary bodies. The tidal heating source is radially resolved and can capture the temperature-dependent pattern of heat production with depth. We use MATH to examine the steady-state thermal profiles through a conductive Europan ice shell, considering the tidal effects of long-term eccentricity variations, obliquity, libration, and nonsynchronous rotation (NSR). In each case, we vary the shell thickness, calculate the combination of tidal heating generated within the shell and heat flux into the shell base that would maintain that thickness, and track the resulting surface heat flux. We find that Europa’s ice shell should be, on average, no more than 12–17 km thick, due to long-term eccentricity variations, and could be even thinner if the basal heat flux is nonnegligible. These results are more consistent with inferences from Europa’s surface geology than previous tidal heating studies that used simplified interior models, a globally averaged tidal heating rate, and/or constant eccentricity. We also find that, for a given stable shell thickness, the surface heat flow remains fairly consistent even as other parameters are varied, perhaps providing an additional method for constraining ice shell thicknesses on ocean-bearing moons. Although Europa’s tidal heat budget and shell thickness seem relatively insensitive to constant obliquity, libration, and NSR, moons in closer-in orbits may be more sensitive to these other effects and should be further explored.


Introduction
Although it has long been accepted that Europa, one of Jupiter's moons, maintains a thin ice shell (relative to its radius) over a thick liquid water ocean, the thickness of the shell has been heavily debated. In particular, there appears to be a substantial discrepancy between ice shell thicknesses inferred from models of observed geologic features and those derived from models of Europa's heat budget (Billings & Kattenhorn 2005;Nimmo & Manga 2009). We begin with an overview of studies that have estimated Europa's ice shell thickness using a variety of approaches. This overview is intended to highlight the challenges faced in constraining Europa's ice shell thickness; it is not intended to be an exhaustive list of such studies. The thickness of Europa's ice shell affects the geologic processes within the shell, which, in turn, control the propensity for material transport between the ocean and surface. Europa's habitability, the plausibility of sampling oceanic material at the surface, and the feasibility of sampling the ocean directly all depend on the overall ice shell thickness, making robust estimates particularly important for our continued exploration of Europa. Billings & Kattenhorn (2005) compiled estimates of ice shell thickness within the literature (see also Nimmo & Manga 2009). The thinnest estimates of hundreds of meters, but as much as several kilometers in some studies, come from methods in which the "flexure" of the ice shell is computed; these values only apply to the elastic portion of the ice shell. Other approaches to matching observed geologic features (e.g., the heights of rafts within chaos features) probe the full depth of the ice shell and tend to result in ice shell thickness estimates of 1-6 km, although two studies stated a few to tens of kilometers were plausible.
Europa's craters provide another method of estimating the ice shell thickness, as some impactors clearly penetrated to the ocean while others did not. Across the 12 crater-based studies compiled by Billings & Kattenhorn (2005), seven find a total ice shell thickness between 2.5 and 10 km, several suggest 10-15 km, and only one concluded that the ice shell is >15 km thick. Schenk & Turtle (2009) report a range of 4-20 km, with a preference for the upper end of the range, although they point out that a combination of temporal variations in ice shell thickness, uncertainties in the thermal and mechanical properties of the ice, and the small number of impacts on Europa contribute to a large uncertainty. More recent hydrocode modeling of impacts into a conductive ice shell over an ocean suggests that Europa's craters are most consistent with a 7 km ice shell (Bray et al. 2014). Another crater modeling study found that the morphological transitions between crater types on Europa could be matched by either an 8 km thick ice shell over an ocean or a 5-7 km conductive lid over a thick, convecting ice shell (Silber & Johnson 2017).
Thermal models have been used to deduce the thickness at which Europa's heat generation is matched by its heat output (i.e., thermal equilibrium). Europa's ice shell thickness is strongly affected by tidal heating produced in its ice shell and, to a lesser extent, by tidal heating in the ocean and heat produced within the rocky core, which enter the ice shell from below (e.g., Nimmo & Manga 2009 and references therein). Tidal heating mainly results from Europa's slightly elliptical (i.e., eccentric) orbit around Jupiter. Although dissipation typically leads to circularization of a moon's orbit, Europa's mean motion resonances with Io and Ganymede allow it to both dissipate energy and maintain its eccentricity (Peale 1999).
However, the eccentricity is not constant: it oscillates around a long-term average (Hussmann & Spohn 2004;Nimmo & Manga 2009). Hence, Europa's interior will transition between states of higher-than-average tidal heating, during periods of high eccentricity, and lower-than-average heating during low eccentricity epochs. Heating variations may also lead to temporal changes in ice shell thickness (e.g., Hussmann & Spohn 2004).
Europa's tidal heating budget can also be affected by obliquity (a tilt of the spin pole), physical libration (a wobble about the sub-Jupiter point), and nonsynchronous rotation (NSR; rotations of the ice shell relative to the direction of Jupiter). Theoretical analysis of Europa's rotation state (Bills et al. 2009), as well as Europa's fracture patterns (e.g., Rhoden et al. 2010Rhoden et al. , 2021, strongly suggest the presence of a nonzero obliquity. Long-period NSR cannot be ruled out, but there is no upper limit on the NSR period, and the shell may instead be synchronously locked (Hoppa et al. 2001;Bills et al. 2009). Detailed studies of tidally driven fractures, called cycloids, suggest that libration may also be nonzero (Rhoden et al. 2010), but there are currently no measurements that constrain its amplitude or phase. Librations can strongly affect the tidal heating budget of a moon (e.g., Rhoden & Walker 2022).
Past thermal models have typically included only heating from Europa's eccentricity (see Section 2). Billings & Kattenhorn (2005) report a range of thicknesses from these models of 13-38 km, which is substantially thicker than estimates based on model-based interpretations of geologic features. More recent thermal modeling studies have reported similar results of 30 km or more (e.g., Nimmo & Manga 2009;Hussmann et al. 2015 and references therein). However, all of these results are susceptible to the effects of simplifications in the modeling approach, as well as uncertainty in many model parameters.
Thermal modeling results also depend on the amount of heat available at the base of the ice shell (i.e., basal heating), which can be estimated (e.g., Tobie et al. 2005;Hussmann et al. 2015), albeit with large uncertainties. Recent work suggests that the ocean itself may generate tidal heating that would enter the ice shell from below (Tyler 2014;Beuthe 2016;Matsuyama et al. 2018;Hay & Matsuyama 2019;Rovira-Navarro et al. 2019). Heat transport at the rock-ocean interface, due to tidal heating in the silicate interior, as well as radiogenic heating, could also be substantial (e.g., Běhounková et al. 2021), but how efficiently that heat can be transported through the ocean is not well known (e.g., Soderlund et al. 2020 and references therein).
A common assumption in thermal models is that Europa's ice shell undergoes convection (e.g., Moore 2006;Hussmann et al. 2015). Convection is a very efficient heat transport mechanism, so much more heat can be generated and transported without causing the ice shell to melt. As a consequence, thermal models that assume convection provide the largest estimates of ice shell thickness reported by Billings & Kattenhorn (2005). Whether the shell is convecting is unknown, particularly due to lack of information as to the mechanical and physical properties of Europa's ice that control convection (Barr & Showman 2009). As a result, the thickness range over which convection might occur spans the range of thicknesses implied by other methods.
To summarize, there is a gap between ice shell thicknesses derived from model fits to observed geologic features, which are typically 7-15 km, and those from thermal models that lack observational benchmarks, which are typically >20 km. Here, we attempt to reconcile these inferences by adopting more sophisticated techniques for computing tidal heating and its transport through the ice shell, and then by considering variations in Europa's eccentricity and additional components of the tide, which affect the total amount of tidal heating over time. We determine the ranges of stable ice shell thicknesses that correspond to tidal heating due to eccentricity variations, obliquity, physical libration, and NSR, most of which have not been included in past tidal heating studies of Europa. Like the eccentricity, obliquity is not constant-both the magnitude of the obliquity and the direction of the spin pole change with time (Bills et al. 2009). However, we address only simple cases of constant obliquity in this work. We find that long-term eccentricity variations have the largest effect on Europa's equilibrium ice shell thickness and substantially reduce the gap between thickness estimates from geologic features and tidal heating models.

Comparison with Previous Tidal Heating Studies
Many authors have computed tidal heating within Europa and estimated the corresponding equilibrium ice shell thickness (e.g., Ojakangas & Stevenson 1986;Hussmann et al. 2002;Nimmo et al. 2003;Tobie et al. 2005;Moore 2006;Beuthe 2016;Hussmann et al. 2015 and references therein;Howell 2021). We have identified several common assumptions or simplifications used in tidal heating studies: (1) assuming constant material parameters for an ice shell that is represented by only one or two layers; (2) applying a globally averaged tidal heating rate; and/or (3) assuming convection is occurring within the ice shell as an initial condition. We describe the effects of these assumptions, which were often required to reach any solution due to mathematical complexity and the limitations of computational resources when the studies were conducted. Furthermore, early studies were often addressing whether any ocean was possible, so conservative estimates and assumptions that might result in thicker shells provided valuable upper limits; the exact thickness of the ice shell was of lesser importance.
Dissipation depends, in part, on the material parameters of the ice shell. Temperature-dependent parameters, such as viscosity, can have a particularly large effect on the total tidal heating and its distribution within the ice shell. For a convecting ice shell, it is reasonable to treat the convecting portion as isothermal, and, therefore, constant viscosity. However, if the ice shell is transporting heat via conduction, the viscosity changes with depth from >10 21 Pa s at the surface to the melting-point viscosity around 10 13 Pa s at the base (Collins et al. 2009). If the ice shell is only represented by one or two sets of material parameters, the variation in tidal heating with depth is lost and tidal heating may be underestimated, particularly if the specified viscosity is not resonant with the mean motion frequency. In addition, calculation of the strain caused by tidal deformation (which is a key component of the dissipation calculation) depends on the rigidity of the ice shell, which is also temperature dependent by way of the viscosity. Several tidal heating papers (e.g., Ojakangas & Stevenson 1986;Nimmo et al. 2003) invoke the strain equations of Peale & Cassen (1978), in which they assume a thin shell with constant rigidity, thereby neglecting the variations in dissipation due to the rigidity variations with depth.
Once the tidal heating is calculated, the equilibrium shell thickness is usually determined by taking a global average of the heat generated within the ice shell, as though all of the heat is evenly distributed throughout. However, the viscosity dependence of tidal dissipation means that the ice shell will not be uniformly heated with depth. In fact, the lower part of the ice shell (where viscosity is low and resonant with Europa's mean motion frequency) will produce much more heat than other parts of the shell. Because the ice near the base of the shell is already close to the melting point, adding heat there can lead to melting and thinning of the ice shell. Global averaging reduces the amount of heat near the base of the shell and increases it in colder parts of the shell such that predictions of melting may be reduced.
Our method deviates from previous work in several ways. We retain additional terms in the tidal potential that are relevant for forcing frequencies other than the mean motion frequency. We also divide the ice shell into a sufficient number of layers that a conductive thermal profile can be reproduced. We calculate heat generation throughout the ice shell and track how that heat is transported across the shell. Rather than prescribing a basal heat flux, we vary the ice shell thickness and calculate the basal heat flux that can maintain each value. The maximum possible ice shell thickness would correspond to a minimum basal heat flux. The minimum ice shell thickness would be zero (i.e., fully melted), although the basal heat flux required to achieve any near-zero thickness would be implausibly high. Given the large uncertainty in the ability of Europa's ice shell to convect (Barr & Showman 2009), particularly for thin shells, we restrict the present study to purely conductive thermal profiles. However, we note that convecting ice shells would be associated with higher surface heat flux signatures than we predict here, and future effort should go to modeling convection, as well.

Methods
Here, we describe the aspects of tidal forcing we explore and provide the model parameters we use for Europa. Details of our mathematical formulation for tidal heating and heat transfer are provided in the Appendix. Our approach is based on the methods of Sabadini & Vermeersen (2004), and, in a broad sense, follows the approaches used in previous tidal heating studies (see Section 2), with improvements to the tide and rheology implementations. Although we will, in the next section, show that our present-day eccentricity model provides a similar ice shell thickness to previous work, we note that the vast differences in the methodologies, reproducibility, and applicability of tidal heating models in the literature make detailed benchmarking very challenging. Any future effort to reconcile these different approaches would be highly valuable to the community.
We test models that include a constant eccentricity at the present-day value, quasi-static eccentricity at values other than present-day, a periodically varying eccentricity, a nonzero constant obliquity, physical libration with the present eccentricity, and NSR at several periods. While our numerical model can also address time-varying obliquity and spin pole precession, the extensive parameter space for combinations of time-varying eccentricity and time-varying obliquity warrants an independent future study.
Our model Europa has a 167 km hydrosphere over a 1394 km rocky mantle/core. We use a density of 1000 kg m −3 for both phases of water (e.g., Moore & Schubert 2000;Tobie et al. 2005), so that phase changes will not affect total mass, and set the interior density at 3824 kg m −3 . This parameterization predicts a total Europa mass of 4.80 × 10 22 kg and 0.346 for the mean moment of inertia. We vary the portion of the hydrosphere that is frozen (i.e., the ice shell thickness). The associated ocean thickness adjusts to preserve the hydrosphere thickness. Physical, thermal, and orbital parameters are provided in Table 1.
We evaluate tidal heating using both Maxwell and Andrade rheologies (see Appendix A.1), with Andrade parameters of α = 1/3 and ζ = 1 (e.g., Castillo-Rogez et al. 2011;Shoji et al. 2013). To set the initial thermal/viscosity structure within the ice shell, we assume a melting-point viscosity of 1 × 10 13 Pa s at the base of the ice shell and a surface temperature of 100 K. We then divide the ice shell into 50 layers, assigning each layer a viscosity between the surface and base values based on a conductive thermal profile (see Appendix A.1). Previous analysis of Mimas showed that tidal heating and ice shell thickness are sensitive to both of these parameters, with higher melting-point viscosity leading to thicker shells and a lower surface temperature leading to thinner shells (Rhoden & Walker 2022). Given the wide range of tidal forcing we are exploring, we have chosen not to vary either parameter in these models.
For Europa, and other worlds with ice shells that are thin relative to their radii, ice shells approximated by two layers are generally sufficient for calculating tidal heating due to eccentricity. However, for bodies with thicker shells relative to their sizes, such as Enceladus or Mimas, additional layers are required to adequately capture the temperature-dependent dissipation with depth. In addition, longer-period forcing will produce maximum dissipation at different temperatures, and thus depths, in the ice shell, which may require additional layers to fully resolve. Here, we use the largest number of layers to illustrate the most generic form of our model. We start by using an eccentricity of 0.0094, which is the current average eccentricity determined from JPL Horizons ephemerides. However, the results of previous work indicate there are scenarios in which the eccentricity may vary periodically over timescales of ∼140 Myr (Hussmann & Spohn 2004). In cases where the eccentricity changes very slowly compared to the thermal equilibrium timescale, we can ignore the inertial forces and treat the eccentricity as a quasi-static constant, varying the exact value we use for the eccentricity.
Because we are assuming conduction, we estimate the diffusion timescale as τ = L 2 /α, where L gives the length over which the heat conducts and α = k/ρC is the diffusivity, which depends on thermal conductivity (k = ∼3 W m −1 K −1 ), density (ρ = 1000 kg m −3 ), and specific heat (C = 2000 J kg −1 K).
Using these values for Europa, we determine the timescale as a function of shell thickness, L (Figure 1), for different ice shell thicknesses. For the sample case we describe in Section 4.1 (35 km ice shell, present-day eccentricity, Andrade rheology), the timescale to achieve thermal equilibrium is roughly 25.5 Myr, which is considerably lower than the ∼140 Myr cycle for eccentricity variations reported by Hussmann & Spohn (2004). We also find that the simulation time to achieve steady state in our thickest shells is about 10 Myr, in rough agreement with the timescale analysis. Hence, we can treat variable eccentricity in a quasi-static way by testing different values of (constant) eccentricity.
We also explore how periodic eccentricity variations affect the conditions for thermal equilibrium by modulating the tidal source term based on the variable eccentricity model of Hussmann & Spohn (2004). In Figure 2, we show their eccentricity evolution along with our approximations of the periodicity. We fit a sinusoid of the form e(t) = e 0 + e A Sin(e n t) p to their eccentricity model, setting the amplitude, frequency, and base eccentricity as e A = 0.02, e n = 7 × 10 −16 s −1 , and e 0 = 0.007, respectively. This frequency corresponds to a period of roughly 140 Myr in the even-powered sinusoid. To adjust the fit, we used two powers: one to examine more broad peaks in eccentricity (p = 2) and another for sharper, more rapidly changing eccentricity peaks (p = 4).
In both cases, this simple model reasonably captures the eccentricity variation found short of about 5 billion years and represents the current epoch. The fits do become too wide (in a) or are too low (in a and b) beyond 5 billion years, but they capture the general variation we aim to examine fairly well. Applying this eccentricity model gives us a long-term average of the thermal environment, about which we expect the instantaneous solutions to oscillate. It is important to note that eccentricity evolution studies are parameter dependent, and different patterns of eccentricity evolution can emerge as these parameters are varied. Here, we explore the behavior of the ice shell in response to periodic eccentricity variations in a generic way, using Hussman & Spohn 2004 and a guide. More robust studies of the Europan eccentricity states would be beneficial for more detailed estimates of Europa's ice shell thickness.
Our analysis includes an investigation into the effects of a constant obliquity, in which the spin pole is not oriented normal to the orbital plane. We define the obliquity as the angle between the orbit normal and the spin pole. We explore obliquity values ranging from 0 to 2°degrees, which should cover the range of possible Europan obliquity values (Bills et al. 2009). We further consider that the spin rate of the body may be influenced by a low-amplitude forced libration in longitude. To examine this effect, we apply a sinusoid to the position of Jupiter in Europa's sky (f s ; see Appendix A.5 for details) with f s → f s − L A Sin(L n t + L p ). Here, L A is the amplitude of the libration, L n gives the rate (which we assume is equal to the orbital mean motion), and L p is the phase of the libration relative to the sub-Jovian location, although tidal heating is independent of this parameter. We evaluate libration amplitudes that range from 0 to 1°, which should be observable if occurring at Europa. We also examine the heating in Europa's ice shell assuming a small departure from synchronous rotation (Ω NSR ). We assume the rotation rate is very slightly greater than the mean motion, given by Ω = n + Ω NSR . We consider NSR periods of 10, 100 kyr, 1, 10, and 100 Myr. Periods at the shorter end of this range are not consistent with patterns of tectonic features and would create enough tidal stress to dominate over the eccentricity-driven stress (Kattenhorn & Hurford 2009;Rhoden et al. 2021), but we include them here for illustrative purposes.

Overview of Tidal Heating Response
To illustrate the effects of heat generation and transport in the ice shell, in Figure 3 we show the temperature, tidal heating, and heat flux profiles for a sample case of a 35 km ice shell, which can achieve steady state with negligible basal heating. The initial temperature and flux profiles (blue lines) represent the equilibrium state of the shell without the tidal source term. In this initial phase, the temperature profile follows our prescribed conductive profile (see Equation (A4) in Appendix A.2). To maintain the specified thickness (here, 35 km) without tidal heating requires a large basal heat flux, which is then transported through and out of the ice shell as a surface heat flux. When the tidal source, at current eccentricity, is applied, the temperature increases near the shell base, as shown by the orange and green curves in the figures, which correspond to Maxwell and Andrade rheology models, respectively. As a result, the thermal gradient becomes less steep at the base, and the basal heat flux is reduced to a value that can maintain the prescribed shell thickness. The increase in temperature also causes the near-surface temperature gradient and associated flux to increase by a small amount. In the center panel of Figure 3, we see the tidal heat with depth for a 35 km ice shell in steady state. The Maxwell rheology produces a higher and narrower distribution of tidal heat, deposited deeper in the ice shell, whereas the Andrade rheology generates a lower but wider distribution of heat, with the peak higher in the ice shell.
In the rest of this section, we report the set of steady-state ice shell thicknesses associated with a suite of tidal effects at Europa. First, we examine the range of thicknesses that are stable for different values of constant/quasi-static eccentricity, and then we explore a sinusoidally varying eccentricity. We also address the effects of constant obliquity, a forced longitudinal libration of the sub-Jupiter point, and, finally, apply a NSR (very slightly faster than synchronous) to examine the effects of a complete long-period rotation of the shell over the ocean. We have found that our results are largely insensitive to the choice of Andrade versus Maxwell rheology for the ice shell thicknesses of interest. Hence, we focus on our Andrade results and provide a short summary of the differences between the two models. Table 2 provides a comparison of the surface heat fluxes which result from a 10 to 25 km ice shell.

Present-day Eccentricity
We find a range of ice shell thicknesses that can achieve steady state at Europa's present-day eccentricity, depending on the basal heat flux. Thin shells can move heat quickly, so do not increase temperature as readily as thick shells do. Thus, thin shells require a very high basal heat, where thick shells can preserve steady state at lower basal flux. We show the range of allowable ice shell thicknesses, and their dependence on the basal heat flux, in Figure 4, which we have dubbed a "zipper plot." The strength of the tide determines the separation between the basal and surface flux curves. Stronger tides, either locational or due to increased eccentricity, will separate (unzip) the basal and surface flux curves, as shown here. Weaker tides, perhaps at a location of low tidal signal or due to near-zero eccentricity, will cause the basal and surface fluxes to collapse (zip) onto each other to coincide with the sourceless solution (dashed lines). The plot also shows the surface flux that each solution predicts.
For Europa's present-day eccentricity, we find that the maximum possible steady-state shell thickness is 36 km using an Andrade rheology. Tidal heating produced within the shell will sustain this thickness with the minimum heat entering the shell from below. By increasing the basal heat flux, we find steady-state solutions for thicknesses between 1 and 36 km. We do not test thicknesses less than 1 km to avoid situations in which the actively melting/freezing zone at the base of the ice shell is a significant fraction of the total ice shell.
It is important to note that, while the tide can draw down the basal flux required to maintain a particular thickness, it has a very minimal effect on the surface flux. Thus, making predictions of surface flux becomes an easy analytic task in most cases (see Appendix A.2). The detail-oriented reader will note that, in the sourceless condition, the basal flux is greater than the surface flux, a result of the increasing surface area with radius.
In general, the two rheology models predict very similar thermal structures in the ice, with Andrade resulting in slightly reduced basal flux, slightly increased surface flux, and a smoother transition between the two regions for a particular ice shell thickness. The maximum shell thickness is 35.5 km with the Maxwell model versus 36 km with Andrade. Overall, we find that Andrade does generate slightly more heat than Maxwell; our steady-state 35 km Maxwell shell globally integrated heat rate is 0.60 TW, while Andrade was 0.65 TW. This may come as a surprise, Andrade producing more heat globally yet providing thinner shells. Locally, however, the tidal heat rate is at least one-half the Maxwell value. Andrade distributes the heat out over the colder temperatures and so heat balance is locally attainable near the shell base without an increase in temperature over the melting temperature. With the very localized Maxwell heat spike at over twice the strength of Andrade, the balance is more difficult to achieve and the temperatures can rise over the melt temperature more easily, and the shell melts back. It is very dependent on where the heat is deposited rather than it is on how much heat is deposited.

Quasi-static Eccentricity
Treating the eccentricity as a constant, we generated zipper plots ( Figure 5) at eccentricities ranging from near-zero (e min = 0.001) to over twice the current value (e max = 0.025). In most cases, the Andrade and Maxwell models predict very similar behavior. However, when the eccentricity is low and the shell is very thick, the Andrade model requires nearly an order of magnitude less basal heating than is found with Maxwell. Shells thicker than 100 km will likely undergo convection, which we do not model here, so these thick shell cases are not realistic. We find that higher-eccentricity cases reduce the range of ice shell thicknesses that can be maintained to overall thinner shells because tidal heating increases with eccentricity. For a particular ice shell thickness, the basal heat flux required to maintain that thickness is reduced as eccentricity increases.
In Figure 6, we have compiled the maximum stable ice shell thickness with eccentricity over the range explored. With eccentricities lower than about 0.0075, the maximum thickness shells become thick enough that convection may be a plausible mechanism for heat transport, in conflict with our assumption of conductive heat transport. Hence, future work examining thermal equilibrium in ices undergoing convection would be warranted.

Oscillatory Eccentricity
With the oscillatory eccentricity models (Figure 2), for which zipper plots are shown in Figure 7, we find that the ice shells cannot be as thick as with the present-day eccentricity model (Section 4.1). The maximum shell thicknesses are 17.5 and 20 km for the "broad peaks" (p = 2) and "sharp peaks" (p = 4) eccentricity models, respectively, using an Andrade rheology. We note that the p = 4 case can tolerate slightly thicker shells, due to the spiky nature of the eccentricity function, which spends less time at high eccentricity than the p = 2 case does. Thus, we infer that on long timescales (e.g., 100s of Myr), the average shell thickness should be no more than 17-20 km, assuming the minimum basal heat flux. If the basal heat flux is consistently high, the average shell thickness will be thinner.
If we compute the average eccentricity of each of our periodic eccentricity curves, we can compare the resulting shell thicknesses to those determined using a constant eccentricity at that average value (e.g., the quasi-static cases in Section 4.2). The average eccentricity values are 0.017 for p = 2 and 0.0145 for p = 4. Using these eccentricity values, we find that the ice shells are close to the thicknesses found using the periodic approach. Because our functions for eccentricity variation are slightly skewed, such that the eccentricity spends more than half the time below the average eccentricity, the shell thicknesses using the constant average eccentricities are 1-2 km thinner than the values we find using the periodic functions.  Note. N/A indicates the shell cannot achieve stability at that thickness.
Assuming that Europa is undergoing periodic eccentricity variations, we find that its average ice shell thickness is controlled by the long-term average of the eccentricity. Hence, Europa's ice shell can be substantially thinner, on average, than the ice shell would be with a constant eccentricity at the present-day value. On shorter timescales, the thickness of the shell will oscillate around the average, seeking the quasi-static solution (Section 4.2), with a value depending on the strength of the basal heat flux and instantaneous eccentricity. Again, we find very little difference between the rheology model results, and we observe the same reduced basal flux and enhanced surface flux in Andrade relative to Maxwell.

Obliquity
We next examine the effects of obliquity on stable ice shell thicknesses. We assume the standard eccentricity tide, at the current value, and apply a hypothetical obliquity ranging from 0°to 2°. As expected, obliquity allows thinner ice shells to achieve steady state, although very large obliquity values would be required to have a substantial effect on Europa's overall heat budget. As shown in Figure 8, we find that obliquity values less than 1°have a small effect, resulting in a reduction in maximum shell thickness of <5 km relative to the eccentricity-only models. The maximum thickness quickly recedes from 35 km with zero obliquity to about 20 km when a 2°obliquity is applied. Larger obliquity values would thin the shell even more, but such large obliquity values are less likely as they should have been detectable if actually present at the system.

Physical Libration
Here, we introduce a hypothetical forced libration in longitude and examine how it affects the thermal structure of Europa's ice. For simplicity, we use Europa's current eccentricity and assume the libration happens at the same rate as the orbital mean motion. For tidal heating, the phase of the libration does not contribute to the result, but the phase would be important for tidal stresses. The important parameter here is the libration amplitude, for which we examine up to 1°in increments of 0.25°. We choose these values for demonstration purposes only; they do not represent observed libration amplitudes for Europa. We include them to examine the effect of libration on Europa, noting the applicability at other ocean worlds like Mimas (Rhoden & Walker 2022). As shown in Figure 9, adding libration reduces the basal heat flux required to maintain relatively thick ice shells (>20 km) but, as with obliquity, large values are required for substantial effects. To achieve a 10 km ice shell requires ∼50 mW m −2 of basal heating with all of the libration amplitudes we tested, assuming constant eccentricity at the present-day value and an Andrade rheology.

Nonsynchronous Rotation
We find that, in all cases, NSR does not affect the steadystate thermal environment within Europa's ice shell, although NSR can have a significant effect on tidal stresses. For the same basal heat fluxes, the shell thicknesses are comparable to runs that did not include NSR (see Section 4.1). The most interesting thing to note is the presence of the "NSR heat bump," which is a low-magnitude signal observed surfaceward of the typical peak in tidal heating. In Maxwell models (Figure 10, left), this bump is obvious in log space, with a position and magnitude that depends on the NSR period. Longperiod NSR has lower-magnitude bumps that appear nearer to the surface, whereas shorter-period NSRs have higher-   magnitude bumps closer to the shell base. With an Andrade rheology (Figure 10, right), the effects of NSR period differences are much less pronounced, and all NSR rates generate nearly identical tidal heating. In the Andrade model, heat is more evenly distributed throughout the shell, so the heat rates are already elevated in colder ices as compared with the Maxwell model, making the effects of additional NSR heat slightly stronger in all cases, however less dramatic in terms of frequency.

Surface Heat Flux
In Table 2, we list the surface heat fluxes for a 10 and 25 km ice shell for the various models discussed in this work. As can be seen, there is very little variation between the different scenarios we model, although there is a small increase in surface flux tied to increasing tidal signal. In general, however, the values are more strongly related to shell thickness. Even when we remove tidal heating, the surface heat flux associated with a particular thickness remains close to the value without tidal heating, but, in that case, the heat comes in the form of basal heating.

Discussion
Our motivation for applying MATH to Europa was to determine whether a more sophisticated rheologic model and/ or more complex tidal forcing could reconcile estimates of ice shell thicknesses derived from Europa's geology with those computed from tidal heating models. We find that the largest effect on ice shell thickness estimates is the treatment of eccentricity, whereas the adoption of an Andrade rheology had a minimal effect. The choice of rheological model could have a more substantial effect on moons with closer-in orbits because the mean motion frequency will control the position of the tidal heating peak within the ice shell.
We find that the periodically varying eccentricity model substantially reduces the maximum ice shell thickness, as compared with constant eccentricity at the present-day value, from 35 to 17 km. The ice shell could be thinner with a larger basal heat flux (see below). Our 35 km result is consistent with past estimates of ice shell thickness from tidal heating models that assumed a constant eccentricity (e.g., Tobie et al. 2005;Hussmann et al. 2006). Neither physical libration at the mean motion frequency, obliquity, nor NSR have substantial effects on the range of possible steady-state ice shell thicknesses or the basal heat flux requirements to achieve a particular ice shell thickness.
The range of plausible basal heat fluxes for Europa is large and highly uncertain. Estimates within the literature include values of 6-7 mW m −2 (Hussmann et al. 2006;Nimmo & Manga 2009) to >25 mW m −2 (Běhounková et al. 2021;Howell 2021), and all of these studies point out that the true value is unknown. The lower end of these basal heat flux estimates would be sufficient to maintain a 15 km thick shell   using the variable eccentricity model, whereas the upper end would correspond to a ∼12 km thick ice shell. The basal heat flux would have to be at least double the highest estimates in order to maintain a 10 km thick shell at Europa. Reproducing as much heat as we observe at tidally heated moons, such as Enceladus, has proved challenging (e.g., Roberts 2015 and references therein), and much is still unknown about heat generated within, and transported through, the oceans of icy moons (e.g., Soderlund et al. 2020). Hence, the strongest conclusion we can draw is that Europa's ice shell ought to be less than 17 km (as a long-term average), as opposed to the thicker estimates found in previous tidal heating studies, and a thickness of 10-15 km is plausible. Thus, we have been able to largely close the gap between geologically inferred thickness estimates and tidal heating analyses, assuming those features were emplaced during high-eccentricity excursions.
We have reported the average ice shell thickness that results from eccentricity variations. However, there will be times when Europa's eccentricity is higher than the present-day value and tidal heating in both the ice shell and the deeper interior are enhanced. Europa's surface age is thought to be 10-200 million years (Bierhaus et al. 2009), whereas the eccentricity variations found in previous works can have a ∼140 Myr cycle (Hussmann & Spohn 2004). Therefore, the shell could record one to several complete oscillations in eccentricity. The surface may, thus, display a combination of features that formed at different shell thicknesses and in both cooling/freezing and warming/melting epochs. Our analysis of quasi-steady-state eccentricity suggests that an eccentricity of ∼0.02 or larger (approximately double the present-day eccentricity) could maintain an ice shell thinner than 10 km with low to moderate basal heating. Therefore, we suggest that the crater population from which a 7 km thick ice shell has been inferred represents the most recent epoch of high eccentricity and tidal heating.
It is also possible that an effect we did not include enhances the tidal heating budget, leading to the thinner average ice shell that is reflected in geologic features. In this study, we have only examined the effects of constant obliquity due to the large (and largely unconstrained) parameter space. However, obliquity variations will affect Europa's tidal heating budget, perhaps enabling thinner ice shells without high basal heating. As more information as to the likely value of Europa's obliquity, the period of obliquity variations relative to eccentricity variations, and the precession rate become available, variable obliquity should be incorporated into tidal heating studies of ice shell thickness. Similarly, we only tested physical libration at the mean motion frequency. Longer-period librations would deposit tidal heating at different depths within the ice shell than eccentricity, which may lead to a different equilibrium solution. Again, more constraints are needed to narrow the plausible parameter space ahead of a detailed libration study.
Our model has high radial resolution but does not preserve global spatial variations in tidal heating (as in Ojakangas & Stevenson 1989) or variations in surface temperature due to, e.g., solar isolation. Tidal heating should lead to regions that are thinner than the averages we compute, and we have found that equilibrium shell thickness is sensitive to the surface temperature, which may result in additional thickness variations (e.g., Rhoden & Walker 2022). However, the warm ice near the base of the shell may be able to flow and smooth out those variations (Nimmo & Manga 2009). Hence, a more sophisticated 3D model is required to track spatial variations in tidal heating, lateral heat flow, and the flow of warm ice in order to determine the extent to which lateral thickness variations can persist. It would then be possible to address whether a thinner ice shell is achievable locally, if not globally.
The amount of heat that can be transported out of the shell will govern whether melting or freezing occurs, and thus set the equilibrium shell thickness. Hence, across all of our models, we find that the surface heat flow associated with a particular ice shell thickness remains fairly consistent, although the proportion of heat coming from tidal heating versus basal heating can vary. An important implication of these results is that surface heat flux measurements could provide precise constraints on Europa's ice shell thickness. Based on our current suite of simulations, we find that the surface heat flux varies by only ∼4-8 mW m −2 for a given ice shell thickness, across many unknown tidal factors (e.g., libration amplitude or NSR). In contrast, the difference in surface heat flux between a 10 and 25 km ice shell is >30 mW m −2 , suggesting that measurements with a precision of approximately 10 mW m −2 or better could be used to constrain ice shell thickness. Although outside the scope of this work, it would be useful to evaluate how sensitive surface heat flow is to other parameters, such as surface temperature, and further refine the measurement requirements needed to determine the ice shell thickness to within a few kilometers. Figure 10. Tidal heat rate in 25 km Europa shells using Maxwell rheology (left) and Andrade (right) for different nonsynchronous rotation rates. While the "heat bumps" are observed in the tidal heating profile, they do not affect the final thicknesses of the shells or substantially alter their temperature profile.

Conclusions
We have presented a sophisticated model of tidal heating that more accurately captures the rheological behavior of the ice shell and treats complex tidal forcing. We find that our new toolkit finds ice shell thicknesses that are more consistent with estimates from Europa's geologic features of ∼7-15 km. An additional outcome of our model development is that MATH can treat a more diverse set of moons than past models, including Enceladus, which has both eccentricity and measured physical libration (Rambaux et al. 2010), and Triton, a zeroeccentricity moon with high obliquity.
When applied to Europa, we find that, at the present-day eccentricity, very high basal heating of ∼80 mW m −2 would be required to maintain an ice shell that is 7 km thick, the thickness inferred from crater modeling (Schenk & Turtle 2009;Bray et al. 2014;Silber & Johnson 2017). For more modest basal heating rates (<25 mW m −2 ), the present-day eccentricity would sustain ice shells ∼18-35 km thick, consistent with previous studies. Accounting for Europa's periodically varying eccentricity substantially reduces the range of stable ice shell thicknesses for modest basal heat flows to 12-17 km; a 7 km ice shell requires 70 mW m −2 of basal heating in that case. These thicknesses correspond to the long-term average thickness that results from the eccentricity variations. However, during higher-eccentricity departures, sub-10 km ice shells become much more feasible. Given the timescales for periodic eccentricity evolution versus Europa's surface age, it is plausible that some surface features formed in thinner shells, during a recent epoch of higher eccentricity. As more constraints on Europa's rotation state become available, we may find that additional tidal heating is produced by obliquity variations and/or additional librations, whereas NSR, constant obliquity, and librations at the mean motion frequency do not substantially thin the shell.
Finally, we have found that surface heat flux is quite sensitive to ice shell thickness, and relatively insensitive to other model parameters that we vary. While more work is needed to fully assess the diagnostic power of heat flow measurements, these results suggest an additional technique by which ice shell thickness, and even basal heat flow, could be constrained through remote-sensing techniques. This work was supported by the NASA Solar System Workings and SESAME programs (grant Nos. 80NSSC19K1026 and 80NSSC19K0613, respectively). The results reported in this paper can be reproduced using the parameter values given in the text, tables, and supporting online materials, along with the equations and descriptions provided in the cited literature. The authors have no competing interests to report.

Appendix
Here, we begin with a general overview of viscoelasticity and how rheological models predict heating at various frequencies. We then discuss our heat-transfer methods and tidal source term. Finally, we describe the tidal model in enough detail to reproduce the results or to apply MATH to other bodies of interest.

A.1. Rheological Considerations
Elastic deformation, which stores energy, becomes viscoelastic deformation, which dissipates energy, with the application of a particular rheology model. These models are generally reminiscent of circuit diagrams, where various elastic elements (springs and dashpots) are combined in series and/or parallel to represent a mass element. This modifies the rigidity, μ, into a complex-valued material parameter through the correspondence principle (e.g., Peltier 1974). In these cases, the real component of the rigidity describes the restorable elastic energy that is related to stress while the imaginary part describes an irrecoverable energy converted to heat through viscous dissipation.
The specific function of viscosity that the rigidity will take depends on the rheology model chosen, but in all cases this introduces the parameters of viscosity, η, and forcing frequency, ω. There are a variety of models that could be applied to Europa's ice. The Maxwell model is the most frequently employed, but recent studies argue that the Andrade rheology is more appropriate because it can capture the anelastic behavior expected in Europa's ice (e g., Castillo-Rogez et al. 2011 where i is the unit imaginary number. The Andrade formulation uses the gamma function, Γ, and the Andrade exponential, α. Here we assume the Andrade parameter ζ = 1. The rigidity, thus, becomes complex as μ = μ r + iμ i , where the real part describes the storage of elastic energy and the imaginary part corresponds to energy dissipated as heat. In both rheology models (with α = 1/3 and ζ = 1), the peak in the imaginary component that controls heat is found when (μ/ωη) = 1. This means that for every order of magnitude increase in forcing frequency, ω, the peak is found at a corresponding lower viscosity, η. That is, different tidal frequencies will produce a peak in tidal heating in different locations, tied to the value of the local viscosity. Viscosity is a dramatically sensitive function of temperature and, therefore, so is the complex rigidity. In Europa's ice, and in planetary bodies in general, there is a very large range of temperatures to consider. For example, Europa's average surface temperature is about 100 K, but the base of the ice shell is at or around the melting temperature of water, 273 K. Typically, an Arrhenius relationship is used to convert between temperature and viscosity (e.g., Shoji et al. 2013). For this work, we use a power-law representation, as described by Ojakangas & Stevenson (1989), but with a power L = 27 (as was used in Hussmann & Spohn 2004), which is appropriate for diffusion creep. Specifically, we use where T M is the ice melting temperature of 273 K, and η M is the ice viscosity at the melt temperature, which we assign a value of 10 13 Pa s by default. The results are sensitive to the choice of melting-point viscosity.
We can now write the imaginary rigidity as a function of temperature, μ i (T), in a simple analytic form. Figure A1 shows a comparison of this form of the rigidity between the two rheological models. In both cases, the peak of rigidity, which is associated with heat generation at Europa's mean motion frequency, is at 246 K-near to, but not quite at, the melting temperature. Thus, we expect that most tidal heating will be generated near the base of the ice shell for forcing at Europa's orbital frequency. Furthermore, the peak of the Maxwell model is twice the value of the Andrade model, but Andrade's peak is broader, generating more heat in colder ices. The parameters T M and η M cause horizontal shifts to this profile. Manipulating α causes similar shifts, but also a change in magnitude. For example, setting α = 1 (with our choice of ζ = 1) shifts the peak to about 240 K with a magnitude equal to the Maxwell peak. In this work, we do not explore the parameter space for T M , η M , or α and use only the values described above.

A.2. Heat Transfer
The method of heat transfer is assumed to be conduction. We use the standard heat equation, which says that a mass element with density, ρ, and specific heat capacity, C, will experience a change in temperature, T, if there is any divergence in the heat flux (k ∇T). In the case of volumetric tidal heat sources, Q, a temperature change, can also result if there is any imbalance in flux divergence and the source. We use an inversely temperature-dependent conductivity such that k = A/T, where A is a constant. Thus, our heat equation is given by where ∂ t indicates the time derivative. We assume that both ice and liquid water have the same density of ρ = 1000 kg m −3 to simplify the model such that phase changes will not affect the total mass of Europa. We use the specific heat capacity of ice C = 2000 J kg −1 K −1 . Our choice of the constant A = 567 W m −1 comes from Nimmo & Manga (2009), who cites Klinger (1980). Our heat-transfer model is in spherical geometry and onedimensional in the radial direction, so the temperature we solve for is a function of both radius and time, as T(r, t). Generally, the volumetric tidal heat source, Q, in units of W m −3 , is a function of all three spatial dimensions and time; however, we average it in colatitude and longitude to simplify the numerical integration, as described in Appendices A.3 and A.4. Ultimately, Q is also a function of radius and time due to its dependence on the local temperature, as in Q(T(r, t)).
We take boundary conditions in temperature at the shell base T(r b , t) = T b = 273 K and surface T(R, t) = T R = 100 K. Here, R = 1651 km describes the body radius of Europa. The radius of the shell base, r b = R − d, depends on the thickness of the ice shell, d, which is the free parameter in this model. It can take a value anywhere from d = 0 km (ocean surface) to d = 167 km (grounded ice shell). For an initial condition in temperature, we use the temperature solution for the heat equation without a source. That is, we solve, ρC∂ t T = ∇ · A(∇T/T) in a spherical geometry to find where, again, T R = 100 K and T b = 273 K. We alternatively refer to this initial condition as the sourceless condition. As always, the temperature gradient is proportional to heat flux.

A.3. Tidal Heat Source
Here we derive the source term, Q, from elastic theory. For an incompressible body, a mass element experiences a volumetric elastic energy density Ξ = με ij ε ij * , where μ is the material rigidity and ε ij is the ith-jth component of the local strain tensor. Summation over all nine components is assumed and the * indicates the complex conjugate operation. When strain can be represented as a Fourier series, like e e = w e ij ij t i (and it can, as we will show below), then the time derivatives are simple, and we can write the volumetric elastic energy We again encounter the imaginary number, i, complex rigidity, μ, and forcing frequency, ω. Notice that, due to the absolute value of a complex exponential, the time dependence of the tidal strain falls out of the solution. Thus, this expression represents the time average of the strain over one cycle in that frequency. When multiple frequencies are present in the tidal strain, this expression is computed individually at each frequency and then summed together. This represents the total elastic energy density rate due to some strains. For elastic materials, the energy is completely recoverable (μ is real). Because we seek dissipation rates of viscoelastic materials, in which the rigidity is complex (μ = μ r + iμ i ; see Appendix A.1), we isolate the imaginary component of the rigidity to write the energy density form for tidal heating (also in W m −3 ) as wm e e = *˜( ) Q 2 .
A 6 ij ij i

A.4. Tidal Strain
The complicated part of the tidal heating model is the strain. Many simplifying assumptions can be made to arrive at analytic results or to speed-up numerical computation. Yet, however grueling, we desire to keep this formulation as general as possible.
Before finding the strain, we define displacement as which assumes that displacements are separable into an angular factor that depends on the tidal potential, Ψ(θ, f), and radial factors, y 1 (r) and y 2 (r). The tidal potential is the applied potential (see Appendix A.5) at a single frequency. To determine radial y-functions, we use the normal mode and propagator matrix methods of Sabadini & Vermeersen (2004).
In this method we can relate the total potential, Φ, stress tensor, σ, and displacement vector, u, using the differential equations Here, (i) is an expression for force balance, (ii) asserts material incompressibility, and (iii) requires that the total perturbation to the forcing potential satisfies the Laplace condition. As above, ρ gives the material density and g is the gravitational acceleration. The stress term introduces the rigidity, μ, and is related to the strain with σ = λ(∇ · u) I + 2 με. Here λ represents the bulk Lamé parameter and I is the three-dimensional identity matrix. Sabadini & Vermeersen (2004) describe the method to solve these differential equations in a medium with homogenous rigidity and density, and then they describe the boundary conditions at internal interfaces that allow for multilayer solutions. We choose a total of 50 layers in the ice shell to accurately capture the thermal gradients expected, which modify the rigidity and viscosity (see Appendix A.1). To include subsurface oceans, as is necessary for Europa and other ocean worlds, we use the boundary conditions noted in Jara- Orué & Vermeersen (2011).
Once the displacement has been found, the simple relationship from displacement to strain is written e  = + ( u 1 2  ( ) ) u T , where the T indicates the transpose operation.

A.5. Tidal Potential
The gravitational potential of a source mass M, at a body-tobody center distance, d, away from another mass of radius, R, is written in spherical harmonics,  Y ℓ, , as where the * indicates the complex conjugate, (θ, f) is the spherical location of interest in body coordinates, and (θ s , f s ) gives the coordinates to M in the sky. This is the standard form for the tidal potential, and all that is needed is to relate the source location, (d, θ s , f s ), to orbital parameters. The distance between bodies is given by the equation of an ellipse, with a being the semimajor, e the eccentricity, and f the true anomaly (discussed later), such that We use a series of rotations to find that the source mass position in the sky is given by where o is the spin obliquity and ζ describes the "longitude" of the spin pole (with zero defined such that the spin pole points to lead the source at periapse, and 90°points the pole toward M).
Here, Ω is the spin rate and the true anomaly, f, is a function of time, t. The true anomaly is a function of time and orbital mean motion (n = (GM/a 3 ) 1/2 ). However, it is not analytic and numerical methods are required. We gloss over the details of this textbook problem but describe the result. By using Newton's method, we can represent the true anomaly as a series in eccentricity up to any desired order and precision. Typically, the eccentricities are small enough that expansion to first order is sufficient for most studies, and we have used first order for our analysis of Europa here; however, we present the result to third order for completeness: The final step is to carry though the expansion in eccentricity for the true anomaly through to the full tidal result. Either an additional series expansion in eccentricity or a Fourier series will produce the same result, which enables the tide to be written as where e is the exponential function, i the imaginary number, and Θ is the coefficient for a particular frequency, ω, in the Fourier series. The sum is over all the frequencies in the tide. While this series expansion completes the eccentricity formulation, we seek to incorporate the effects of additional time-varying processes, a few of which we describe briefly. Libration, NSR, and spin pole precession all affect the tidal pattern relative to the latitude-longitude grid applied to Europa, as well as introduce additional energy to the system.
A forced libration in longitude can be included in two equivalent ways. The first is to include an additional term on the spin rate angle, such that Ωt becomes Ωt-L A Sin(L n t + L p ), where L A is the amplitude of the libration (assumed to be dependent on eccentricity), L n is the libration rate, and L p is the libration phase. The other method is to simply modify the position of M in the sky with f s becoming f s − L A Sin(L n t + L p ). This must be applied before the expansion in eccentricity and L A must also be a (typically linear) function of eccentricity. W - NSR is a process by which a very-long-period (relative to the spin period) drift of the surface around the spin pole makes the spin rate just slightly different from synchronous (Ω = n when synchronous). This long-period process corresponds to a lowfrequency perturbation to the spin rate such that Ω becomes n + Ω NSR .
Finally, we can consider a precession of the spin pole by letting the "longitude" of the spin pole, ζ, be a function of time.
Here we use ζ = ζ 0 + ζ n t, where ζ 0 tells us the initial orientation of the axis and ζ n describes its rate of change in time. Note that this considers precession of the axis at a fixed obliquity value.
Eccentricity, obliquity, and variations in both quantities affect the magnitude and lateral excursions of the tidal bulge.
The precession above can be combined with a sinusoid in obliquity so that both directions of orientation of a spin pole can be varied. By representing obliquity as , the spin pole can be oriented in multiple directions at once to track the movement of the pole. The subscripts 0, A, and n indicate the initial obliquity, amplitude, and time rate, respectively. The p power on the sine function can be used to control the shape of the function.
Orbital eccentricity can also be variable, represented as + ( ) e e e t Sin A n p 0 . In this case, the transformation must be applied after the expansion in eccentricity. The subscripts 0, A, and n are again initial eccentricity, amplitude, and time rate, while the power p controls the function shape.

A.5.1. Specific Interesting Examples
The number of frequencies and terms that appear in the tide can be quite large when all the considerations above are included. Most problems will not require all of these, and it will be fine to consider one or two processes at a time. In this section, we build up a few of the most common configurations of the tide, starting with the simplest and building up from there. In all these examples, we choose the degree two tide and expand to first order in eccentricity.
Included for most cases, we present a "Tidal Table." This table consists of rows that depict the coefficient, Θ, frequency, ω, as well as the average power, J, contained within each term. The power is important to note as the rate of tidal heating is dependent on this quantity. We compute the tidal power by averaging each term's absolute value over latitude and longitude using . In all cases, the coefficients for terms at opposite frequency will always be complex conjugates of one another. Thus, the tidal power for terms with opposite frequency are always equal. We also report, in the first row of each result, the static tide at zero frequency. This represents the time-independent component of the tide, and therefore is used to identify the average tide and body figure but does not contribute to stress or heat.

A.5.2. A Classic Synchronous Rotator
The first, and most simple case, to consider is a classic synchronous rotator. Here, the spin rate has synchronized to be equal to the mean motion (Ω = n). There are no obliquity effects or forced libration, and the eccentricity is assumed fixed to some constant value. For this configuration, the tide includes three terms, as shown in Table A1.

A.5.3. An Oscillating Eccentricity
When we consider a tide where the eccentricity varies with time, the form for the classic synchronous rotator undergoes spectral splitting. As described in Sections 3 and 4.3, we employ a sinusoidal form for the eccentricity to emulate the results of Hussmann & Spohn (2004). Specifically, and symbolically, the eccentricity is modeled with e(t) = e 0 + e A Sin(e n t) p . The specific values the parameters take are described in Section 3 of the text, whereas here we keep these parameters symbolic. When p = 2 is selected, for the broader eccentricity peaks, the tide splits as described in Table A2. When the p = 4 form is used, the result is as shown in Table A3.  A synchronous rotator with some nonzero obliquity will contain at least 11 frequencies in its tidal spectrum, {5n, 4n, 3n, 2n, n, 0, −n, −2n, −3n, −4n, −5n}. When precession is included, the tide contains 27 frequencies, and, when an asynchronous rotator has obliquity and precession, the tide splits into 45 frequencies. Given the large number of terms, and a strong dependence on the value of obliquity, we examine these tides graphically rather than in table form.
We present a synchronous rotator with obliquity in Figure A2. The plots contain the power of each "positive" tidal frequency over the range of possible obliquity values. The "negative" frequencies all have the same power as their "positive" counterpart, so those are not displayed. What is most evident in comparing Figure A2, top and bottom, is the very muted effect on the tidal power because of (a rather large) eccentricity. Furthermore, it can be seen that, except for the static tide, the power is very small at low obliquity.
This trend is recovered in the results (Section 4.4) where we see that the small values of obliquity explored in this work do not affect the thermal stability of Europa's ice shell significantly. Note that the gold curve corresponding to frequency n, with the apparent near-zero valued tidal power at o = 0, is in fact the standard tide of Appendix A.5.2 that was used for Section 4.1. Thus, it does still have a value, but one that is very close to that with no obliquity. We can see, however, that large values of obliquity will cause tremendous tidal signals. At o > π/2, for example, these additional frequencies in the tide are very large. In future work, we desire to conduct a similar study with obliquity tides at Triton as a result of these findings.

A.5.5. A Forced Librator
When the effects of libration are considered, the tide picks up two additional frequencies (Table A4), at plus and minus the libration rate, L n . In this case, it is assumed that the libration amplitude, L A , is a linear function of the orbital eccentricity. Hence, the inclusion of these new terms is the only change from the classic synchronous body. Note that when the libration amplitude is zero, L A = 0, those rows vanish, and the synchronous result is recovered.
rates (also identified as eccentricity-dependent terms) condense into the mean motion frequency, and the ±(2n − 2Ω) frequencies (with no eccentricity dependence) collapse back to the zero frequency. When terms with like frequency are summed, the result of Table A1 is achieved.

A.5.7. A Nonsynchronous Rotator
If the spin rate is close to synchronized, but just slightly greater than the mean motion, a slow drift reorients the surface relative to the subprimary location. We set the spin rate to n + Ω NSR ; the resulting tidal terms are shown in Table A6. The tidal terms have the same coefficients and powers as the asynchronous rotator (Table A5), but the frequencies are different. Here, because Ω NSR = n, most of the frequencies are at about the mean motion rate. The most important difference is in the ±2Ω NSR frequencies, which are very-high-power terms in the tidal spectrum, due to their independence from eccentricity. This means that the NSR terms deliver a lot of power from the tide and should be of general interest to tidal studies.