Tidal Heating Did Not Dry out Io and Europa

The Galilean satellites exhibit a clear trend in composition from the rocky Io, close to Jupiter, to the icy and distant Callisto. Proposed causes of this trend can be roughly divided by when the trend developed—either as a result of the material accreted, as a byproduct of the accretion process, or due to the subsequent diverging evolution of the moons. While the first two options have been heavily favored in the existing literature and were therefore studied previously, in this work, we directly address the last of these possibilities. To do so, we determine the range of plausible tidal heating experienced by these moons and how efficiently that energy could be converted into mass loss. We find that while the total tidal energy does exceed the energy required to lose an ice shell, the loss process would have to be highly energy efficient. Examining a range of loss processes, we find that only in extreme cases could enough mass be lost from Europa and no cases where enough mass could be lost from Io. We conclude that this compositional gradient must have been in place by the end of accretion.


Introduction
One of the most striking features of the Galilean satellites is the monotonic trend in composition with distance from Jupiter, from the rocky Io to the icy Callisto. Despite this trend being known for over 40 years, no consensus explanation has been found. Explanations for this compositional gradient can be roughly separated into three categories based on when the gradient was established: first, a strong "ice line" effect, wherein icy material sublimated before being accreted onto Io (Lunine & Stevenson 1982;Ronnet et al. 2017). Second, the water ice could have been lost from the surface during the accretion process by impacts or atmospheric processes (Mosqueira et al. 2010;Dwyer et al. 2013;Bierson & Nimmo 2020). Third, the ice could have been lost over solar system history due to the extreme tidal heating experienced by Io and Europa. While this final possibility has been mentioned repeatedly in the literature (Canup & Ward 2009;Dwyer et al. 2013), very little has been done to quantify its plausibility. In this work, we aim to directly assess if tidal heating could have significantly altered the bulk composition of Io and Europa after their formation.
In this work, we consider that all the Galilean satellites formed with the same bulk density of 1800 kg m −3 (ice mass fraction of 34%). This would imply that the initial radii of Io and Europa were significantly larger than they are today ( Table 1). For convenience, we will refer to the hypothetical early state of these moons as Icy Io and Icy Europa.
Tidal heating is the driver for the geology observed on both Io and Europa today. On Io, a heat flux of roughly 2 W m −2 is accommodated by advecting large amounts of molten rock to the surface in "heat-pipe" volcanism (Moore et al. 2007;Lainey et al. 2009;Veeder et al. 2012). Europa's heat flux is less well known but likely between 0.02-0.1 W/m 2 (Hussmann et al. 2002;Nimmo & Manga 2002;Ruiz & Tejero 2003). For both Io and Europa, these heat fluxes are orders of magnitude larger than radiogenic heating ( Table 2). It is also important to note that these heat fluxes are fairly insignificant compared to the solar heating of 54 W m −2 . Consequently, the average surface temperatures of Io and Europa are very similar except where the heat flow originating from the interior is locally concentrated at the surface, like at volcanic centers.
In Table 2, we compute the minimum energy needed to transition the bulk density of Io and Europa from 1800 kg m −3 to their present-day values. The two large energy sinks are the energy required to move material from the satellite surface to Jupiter's Hill radius (E G , Appendix B) and the latent heat of vaporization for that ice (E L ). For these worlds, the gravitational energy required to remove this mass is about 10 times larger than the latent heat of vaporization. We can take modernday estimates for the tidal heating at Io and Europa and assume that heating was constant over solar system history for a first estimate of the total energy deposited. Ratioing this energy input with the energy needed (E G /E tide ) provides the efficiency with which that energy would have to be used to remove the initial mass of ice. For these simple assumptions, we find these satellites could have lost an initial ice shell if the average energy efficiency of mass loss was 3% for Io and 30% for Europa over solar system history ( Table 2).
In this work, we leverage this general energy efficiency approach. In Section 2, we consider if the average tidal heating rate differed significantly from their modern values and the impacts of tidal heating variations. This is done to more thoroughly estimate the energy available to drive mass loss. In Section 3, we estimate the energy efficiency of different escape processes including geysers and atmospheric escape. Section 4 discusses the implications of these results and places them into the larger context of icy satellites in the solar system.

Tidal Heating over Solar System History
To establish a constraint for the necessary energy efficiency, we need to estimate how much total energy has been dissipated in both Io and Europa. For a tidally locked satellite, the amount of energy dissipated by tidal heating scales as where ω is the orbital frequency, R is the satellite radius, e is the eccentricity, and k Im 2 ( ) is the imaginary part of the tidal potential Love number (e.g., Peale & Cassen 1978;Peale et al. 1979).
k Im 2 ( ) describes the efficiency with which the body dissipates energy and depends on the internal structure.
Tidal dissipation is a strong function of orbital eccentricity. In a two-body system, tidal heating takes energy from the orbit, thus decreasing the eccentricity. Io, Europa, and Ganymede however are in a Laplace resonance. In this configuration, satellite-satellite interactions increase their eccentricity. This competition between tidal heating and satellite interactions allows the nonzero eccentricities to be maintained over solar system history. In this case, the total heating rate in the satellites does not depend on their k Im 2 ( ), but on the k Im 2 ( ) of Jupiter (Meyer & Wisdom 2007). How the heating is distributed between the satellites in resonance does still depend on their internal structures; however, the inner satellites will always experience the most heating due to the strong dependence on ω. The k Im 2 ( ) of Jupiter is also poorly understood and may itself be a function of time and orbital frequency (Fuller et al. 2016). Even with uncertainty in Jupiter's k Im 2 ( ), over gigayear timescales, the total dissipated energy in Io and Europa is not dependent on the satellite structure. Therefore, we expect the average heating rate to be similar, with a factor of a few, to their modern values.
In addition to the physical argument above, we can use numerical models to estimate the mean and variations in tidal heating these moons may have experienced. Equation (1) establishes two potential causes for different amounts of tidal heating in the past, changes in the internal structure, and changes in the orbit. We will now address each of these in turn.
Physical arguments have been used to suggest that the Laplace resonance between Io, Europa, and Ganymede was established during or shortly after their formation (Peale & Lee 2002). We will adopt this assumption because any delay in the formation of the Laplace resonance would lead to less total energy dissipation. Because the tidal heating rate both depends on, and actively changes, the internal structure and orbital eccentricity, a cyclic behavior can arise (Ojakangas & Stevenson 1986). This would have led both Io and Europa to have experienced periods of high and low eccentricity over their histories. Integrations of the coupled thermal-orbital histories suggest that Io's heating rate may have varied by roughly an order of magnitude while Europa's only varied by a factor of 3 (Hussmann & Spohn 2004). These models include the variations in both ω and e as the orbits evolve. Hussmann & Spohn (2004) find that short spikes in tidal heating are followed by much more prolonged periods of low tidal heating. Hence, we infer that the long-term mean dissipation is close to or less than the modern value.
The other key controlling parameters in Equation (1) are R and k Im 2 ( ). We forward-calculate k Im 2 ( ) for a variety of iceshell thicknesses assuming that little dissipation occurs in the ocean layer (Chen et al. 2014;Hay & Matsuyama 2019). We calculate k Im 2 ( ) assuming a Maxwell rheology using the matrix propagation method of Segatz et al. (1988). We use a simple two-layer model for the present-day Io with a core surrounded by a low-viscosity mantle, allowing us to reproduce the observed heat flux of 2 W m −2 ( Table 3). The assumed viscosity of 10 15 Pa s may be unrealistically low but compensates for the fact that the Maxwell rheology underestimates dissipation at high viscosity (Bierson & Nimmo 2016;Renaud & Henning 2018). We then expand this model to a radius of 2600 km by adding an additional 780 km H 2 O ice layer (Icy Io from Table 1). By choosing the viscosity of the ice Note. Rock mass fractions calculated assuming an ice density of 920 and rock density of 3500 kg m −3 . Icy Io and Icy Europa are terms used throughout this work for the hypothetical properties of these moons if they formed with a bulk density of 1800 kg m −3 . Note. For the estimates of ice-shell thickness, the thermal conductivity is assumed to be 3 W m −1 K −1 . Calculations for the gravitational energy required include the change in mass over time and the impact of Jupiter's gravity (details in Appendix B). L v is the latent heat of vaporization of water (2.5 × 10 6 J kg −1 ).

Table 3
Interior Structure used for Calculating the Im(k 2 ) Results Shown in Figure 1 (Equation (1)) Core 800 7500 50 × 10 9 10 30 Mantle 1820 3300 2 × 10 15 50 × 10 9 Ocean Varied 1000 0 0 Ice shell 2600 1000 1.5 × 10 14 3.3 × 10 9 Note. The core and mantle values replicate Io's observed bulk density and tidal dissipation. The two bottom rows are only used in the Icy Io case. The outer ocean radius is chosen dynamically as a function of the ice thickness studied.
such that we match the Maxwell time to the forcing frequency, 1.5 × 10 14 Pa s, we can maximize the tidal dissipation.
The results of this modeling are shown in Figure 1. When the H 2 O layer is entirely frozen (no ocean), the resulting dissipation would be above 32 W m −2 . However, in steady state, this dissipated energy has to be transported to the surface. For modern icy worlds, this is achieved by conductive or convective transport through the ice shell (Moore 2006;Robuchon & Nimmo 2011). The very high dissipation resulting from this setup would quickly melt the ice and leave Io with a very thin equilibrium ice shell of around 400 m. In this new state, 97% of the dissipation still occurs in the silicate part of the body, due to the small volume of the frozen ice shell. In this configuration, we obtain a decrease in k Im 2 ( ) from 1.4 × 10 −2 to 2.5 × 10 −3 compared to Io without the H 2 O layer. As a result, the total energy dissipation is nearly the same in our Io and Icy Io cases (10 14 W) but due to the increased radius, the surface heat flux decreases to 1.2 W m −2 .
While the modeling shown in Figure 1 is focused on an Io analog, the same principles apply to Europa. On modern Europa, most of the heat dissipation is expected to occur in the ice shell (Moore 2006;Sotin et al. 2009). The ice-shell thickness is itself dependent on tidal heating, not the rock to ice ratio. If Europa formed with a high ice fraction the increased radius (Table 1) would cause the volume of the ice shell (and therefore total energy dissipated) to increase by less than a factor of 2.
In summary, even considering dramatic changes in internal structure, we find the mean tidal heating rate for Io and Europa over the solar system history is unlikely to be significantly different from their modern values.

Energy Efficiency of Mass-loss Processes
Generally, there are two pathways that mass could have been removed from Io and Europa over solar system history. First, the mass could have been removed directly through some eruptive process. We will use the term geysers to describe this process although we are agnostic to the exact physical mechanism causing the eruptions. Second, a water-rich atmosphere could form and the mass would be lost via atmospheric escape processes. These are each address in detail below.
Regardless of the mechanism, the minimum escape rate that could remove a primordial ice shell is simply given by the mass of the ice shell over the age of the solar system. This comes out to 10 13 kg yr −1 for Io and 4 × 10 12 kg yr −1 for Europa. Because these escape rate estimates assume that the process operates over all of solar system history, if escape only occurred during the heating spikes modeled by Hussmann & Spohn (2004), the required escape rate may need to be at least an order of magnitude or greater.

Geysers
The first potential pathway for losing material from Io and Europa is eruptive plumes. If Io formed with a thick ice shell, the volcanic plumes we observe today would occur at the base of the ocean. Therefore, we will only consider water geysers erupting through an ice shell. The case example for geysers on any icy world is Enceladus. The plume activity at Enceladus is driven by tidal heating ) and releases enough material to form Saturn's E-ring. Some observations have suggested that larger geysers may be active on Europa (Roth et al. 2014).
Because the geysers at Enceladus are the best-studied ones, we will structure our discussion around those. Enceladus is estimated to have a globally averaged heat flux of ∼200 mW m −2 (Spencer & Nimmo 2013), much higher than that of Europa. Ingersoll & Ewald (2011) used images of Enceladus's geysers to estimate the velocity distribution of erupted material. They found that both Lorentzian and exponential distributions well fit the observed data. Ingersoll & Ewald (2011) then numerically integrated the trajectories of erupted material to determine the mass fraction lost. These dynamical integrations provide a similar result to simply determining the fraction of the velocity distribution that exceeds the escape velocity. For this discussion, we will estimate the mass loss by the fraction of the velocity distribution that exceeds the escape velocity. For full details, see Appendix C.
At Enceladus, 10%-20% of the erupted material is lost forming the E-ring. If eruptive geysers had a similar characteristic eruption velocity (90 m s −1 ) on our Icy Io or Europa the fraction lost would be =5%. With escape fractions so small, the energy efficiency depends on what occurs to the material that does not escape. The bound material will likely fall back on the surface and radiate its thermal energy to space. So unless the gas portion of the erupted material could escape via other processes (discussed in Section 3.2), the energy efficiency of this process is extremely poor. In order for eruptive geysers to have notable mass loss and energy efficiency, the characteristic eruption velocity needs to be 1 km s −1 .
Given this, the only plausible way for geysers to support mass loss is for them to sustain an atmosphere where other processes drive the mass loss (Section 3.2). In this case, the erupted gas mass must at least match the atmospheric escape rate. If the gas to solid ratio in the eruptions is similar to Enceladus, ∼50% (Ingersoll & Ewald 2011), then the eruption mass would need to be twice the atmospheric loss rate. While we do not directly calculate the energy loss associated with  Table 3. The red curve shows the Im(k 2 ) on the left axis and the resulting heat flow given Io's modern orbit along the right axis. The blue curve is the conductive heat flow (right axis) of an ice shell with the thermal conductivity of ice of 3 W m −1 K −1 . Tidal dissipation and conductive heat flux balance with an ice-shell thickness of 400 m. this, it is another potential process that will reduce the energy efficiency calculated below.

Atmospheric Loss/Jeans Escape
At present, Io loses about 10 3 kg/s of material due to iondriven atmospheric sputtering (Wilson et al. 2002). To go from our assumed Icy Io state to Io's current bulk composition requires a mass-loss rate of about 3 × 10 6 kg/s (M L /4.5 Gyr; Table 2) over all of solar system history (half that for Europa). Io's ion sputtering rate could have been higher in the past if Io had a thick atmosphere with a high exobase (Smyth & Combi 1988;Johnson 2004); however, the exobase would need to be several Io radii in height before the flux would be large enough to notably change the bulk composition. This is consistent with Johnson (2004), who estimated that Io could have lost as much as 10 19 kg over solar system history (three orders of magnitude less than the 4×10 22 kg needed). Ion escape is even less plausible for Europa which is estimated to have lost only 6 × 10 17 kg (Johnson 2004), four orders of magnitude less than required.
When considering these ion-driven escape rates in the context of losing a massive ice shell, there are two other important factors. The first is that Io's current atmosphere is primarily SO and SO 2 produced by volcanic outgassing. If Io had a thick ice shell, these volcanic volatiles would likely not reach the surface. In this case, Io's atmospheric pressure could be significantly lower than the modern value, ∼10 −4 Pa, maybe becoming closer to that of Europa, ∼10 −7 Pa (Bagenal & Dols 2020). The other factor is the feedback between the atmosphere itself and the plasma ions that cause sputtering (Bagenal & Dols 2020). While modeling these interactions and feedbacks is beyond the scope of this work, we are not aware of any work suggesting the sputtering rate would have been three to four orders of magnitude faster for most of solar system history. Bierson & Nimmo (2020) suggested that atmospheric escape could drive mass loss during the formation of the Galilean satellites. They suggested that the dominant atmospheric massloss mechanism was hydrodynamic escape because the gravity of the satellites was lower (as they had not finished accretion) and atmospheric temperatures were kept high by a large energy flux. This energy flux was a combination of radiation from the opaque proto-Jovian disk and the accretion of new material. When considering atmospheric loss over solar system history, none of these factors are in place.
Given the above, the only atmospheric loss mechanism that could have any chance of removing this much mass is Jean's escape. Jean's escape is atmospheric mass loss by the fraction of molecules that are moving faster than the escape velocity above the exobase. The escape rate (#/s) for Jean's escape is given by where n is the number density, v esc is the escape velocity, and v T is the most probable speed in a Maxwell-Boltzmann distribution. All values (n, v T , λ) are at the atmospheric exobase. The λ term is the Jean's parameter and is given by where r is the radius of the exobase. Because these moons are deep in Jupiter's gravity well, it is important to calculate the escape velocity to the Hill radius of the satellite with respect to Jupiter (not infinity). We do so as where r H is the Hill radius.
Applying the minimum escape rates (Section 3), at any given temperature, we can calculate the exobase number density (or equivalently atmospheric pressure) needed to sustain that escape rate (Figure 2). For Io and Europa's modern surface temperature of 100 K, the required pressures are preposterously high, requiring more than all of the initial ice mass to be in the atmosphere. If, however, the atmosphere was supported by geysers (temperatures between 250 and 300 K) Io would require pressures of ∼500 Pa and Europa of ∼0.5 Pa. In both cases, those pressures are subsaturated and so the atmosphere would not be actively condensing onto the surface until it cooled. The central challenge with a geyser-supported atmosphere is that the blackbody radiation for those temperatures is in excess of 200 W m −2 , which is unsustainable by tidal or solar heating. The only alternative way that such an atmosphere could persist is for it to radiate inefficiently by having a low infrared emissivity.
The energy fraction going to escape can be calculated by comparing the energy-loss rate of atmospheric escape with that of radiation: x f Here, f is the mass-loss flux in kg s −1 m −2 . If the atmospheric emissivity, ò, were constant, ξ would strongly decrease at higher temperatures as more energy would be radiated. ò is itself a function of temperature and pressure. We use the HITRAN database (Appendix D) to calculate emissivities for all temperature and pressure conditions shown in Figure 2. Within our parameter space, the dominant control on emissivity is the change in atmospheric pressure needed to maintain escape. Figure 3 shows the resulting energy efficiency for Io and Europa compared to the thresholds derived in Section 2. For both satellites, only at high temperatures (corresponding to pressures 100 Pa) would escape be energy efficient. To achieve an efficiency greater than 3% for Io, this atmosphere would need a temperature in excess of 300 K, significantly above the temperature of a presumed subsurface ocean. For this reason, we conclude that there is no plausible mechanism rage of parameters where escape would be energy efficient enough for Io.

Impact of Assumed Temperature and Pressure Structure
In the above calculations, we assumed an isothermal temperature structure for the atmosphere. A more realistic temperature structure would decrease with altitude (either adiabatically or following the saturated vapor pressure via a moist adiabat). This will have two effects. First, the radiating temperature would drop. However, for any plausible scenario atmospheric radiation needs to already be negligible (achieved by being optically thin). Therefore, this will not impact the plausibility of any escape pathway. The second effect of that lower temperature is that it will reduce the escape rate, thus requiring higher temperatures to remove enough mass.
In addition, it was assumed that the minimum exobase pressure needed to maintain escape was the surface pressure. This is an extremely generous lower limit for two reasons. First, the escape rate scales linearly with atmospheric density (Equation (2)), and thus, if the mass was lost over a shorter period of time, the density would need to be correspondingly larger. Second, for the low-temperature high-pressure cases, the exobase would have to be significantly removed from the surface. This would imply surface pressures orders of magnitude larger than those in Figure 2, enhancing the concerns of radiative cooling.
Importantly, both of these assumptions greatly favor higher energy efficiencies. Therefore, we consider the energy estimates in Figure 1 to be generous upper limits.

Discussion and Conclusion
In this work, we made generous assumptions to determine if there was any plausible mechanism by which Io and Europa could have evolved from an initial icy moon into the rocky worlds observed today. For both moons, the total tidal dissipation is larger than the energy required to remove the mass of an ice shell; however, the mass-loss process would need to be highly energy efficient. The gravity of these moons is too large for geysers to cause any significant mass loss. For Io, we also find no plausible atmosphere that could drive mass loss. For Europa there is a narrow sliver of parameter space, where such evolution is possible, exists. Such an evolutionary pathway would require the existence of a water vapor atmosphere with a pressure of a few Pascals and a temperature near 270 K over most of Europa's history. Such an atmosphere would need to have been supported by constant eruptions of geysers. If the geysers contained the same mass estimated for the Roth et al. (2014) observation (10 6 kg), millions of eruptions per year for all of solar system history would be required. Then, the conditions driving these constant eruptions must have ceased long enough ago that the ring of ejected material around Jupiter had been dynamically lost.
Given that an early water-loss mechanism is necessarily required for Io and the incredible specificity required for longterm loss of water for Europa, we propose that it is far more likely that both satellites have changed their bulk composition very little since their formation.
We thank two anonymous reviewers for their helpful feedback. Table 4 presents all variables used throughout this text. ( ) ) is constant as 40 mW m −2 for Io and 20 mW m −2 for Europa. The efficiency is higher at higher temperatures because the atmospheric pressure needed to maintain the minimum escape rate is lower ( Figure 2). This lower pressure leads to lower atmospheric emissivity ( Figure 5) and therefore less energy lost to radiation.

Appendix A Variables Used
Here, ρ i is the density of ice (920 kg m −3 ), ρ 0 is the modern satellite density, and r 0 is the modern satellite radius. Equation (B4) was numerically integrated from the ice-rich radii to the modern satellite radii (Table 1).

Appendix C Geyser Calculation Details
Ingersoll & Ewald (2011) found that the velocity distribution of particles in Enceladus's plumes are best fit by either a Lorentzian or Exponential distribution. We use these same velocity distributions to estimate the mass loss and energy efficiency per eruption Table 5. We find that integrating the fraction of the velocity distribution above the escape velocity provides comparable estimates for the mass fraction lost as the dynamical integration of Ingersoll & Ewald (2011). Because the escape velocity is so much higher on Io and Europa compared with Enceladus, a significantly lower fraction of material is lost. For a comparable mass fraction to Enceladus to be achieved, the characteristic eruption velocity, v 0 , has to approach 1 km s −1 .
To fully determine the energy efficiency would require knowledge of the solid to vapor fraction and temperature of those components. We can however assess the energy efficiency related to the velocity distribution itself. This is done by integrating over the kinetic energy in the entire plume and comparing that with the minimum kinetic energy required to lose the escaping fraction of the plume, i.e., where ξ is the energy efficiency, f lost is the mass fraction of escaped material, and v esc is the escape velocity. This integration unfortunately cannot be done for the Lorentzian distribution because the term v dp dv 2 does not go to zero as the velocity approaches infinity (see Figure 4).
Using the exponential velocity distribution we find that for v 0 = 250 m s −1 , more than three times the velocity of the Enceladus plumes, the energy efficiency is 1%, well below the required 10% efficiency. For v 0 = 1000 m s −1 , the energy efficiency does get to 27%. We also stress that this is the upper bound on the energy efficiency as we are not including the thermal energy or latent heat. From this, we find the the plumes would need to have characteristic eruption velocities approaching 1000 m s −1 to significantly account for Europa's mass loss. As the escape velocity of Io is similar (2300 m s −1 ), these same results hold for Io as well.