Tidal Heating in a Subsurface Magma Ocean on Io Revisited

We investigate the tidal dissipation in Io's hypothetical fluid magma ocean using a new approach based on the solution of the 3D Navier‐Stokes equations. Our results indicate that the presence of a shallow magma ocean on top of a solid, partially molten layer leads to an order of magnitude increase in dissipation at low latitudes. Tidal heating in Io's magma ocean does not correlate with the distribution of hot spots, and is maximum for an ocean thickness of about 1 km and a viscosity of less than 104 Pa s. Due to the Coriolis effect, the k2 Love number can depend on the harmonic order. We show that the analysis of k2 may not reveal the presence of a fluid magma ocean if the ocean thickness is less than 2 km. If the fluid layer is thicker than 2 km, k20 ≈ k22/2 ≈ 0.7.


Introduction
Jupiter's moon Io is the most active volcanic body in the Solar System, with more than 400 known volcanoes, 150 of which are erupting at any given time (e.g., Radebaugh et al., 2001;Schenk et al., 2001;Veeder et al., 2012).Io's volcanic activity is driven by the heat generated by tidal friction caused by its orbital resonance with Europa and Ganymede (Peale et al., 1979).The average endogenous heat production is estimated to be of the order of 100 TW (e.g., Lainey et al., 2009;J. R. Spencer et al., 2000;Veeder et al., 1994), which is significantly more than the heat output of the Earth.
The dissipative response of the body to tidal forcing is determined by its internal structure, size and the frequency of forcing.Dissipation of tidal energy can occur in different ways, in both the solid and liquid regions of the body.In a solid material, the dissipative properties depend on the elastic moduli, the viscosity and the anelastic response, while the dissipation in a liquid is controlled by a single parameter, viscosity, varying with temperature and composition.
Physical models of solid body tides on Io usually assume that most of the heat is generated in a partially molten layer beneath Io's lithosphere or in a deeper, potentially dissipative mantle (e.g., Bierson & Nimmo, 2016;Hamilton et al., 2013;Kervazo et al., 2021;Renaud & Henning, 2018;Ross & Schubert, 1985;Segatz et al., 1988;Steinke et al., 2020).The presence of a partially molten layer in Io's upper mantle is predicted by the models of magmatic heat transfer (e.g., Moore, 2001Moore, , 2003;;D. C. Spencer et al., 2020;Steinke et al., 2020) and is consistent with the estimates of eruption temperatures indicating that a substantial portion of Io's mantle is partially molten, with porosity between 3% and 25% (Keszthelyi et al., 2007).
An alternative model to explain Io's heat output has been proposed by Tyler et al. (2015).The model assumes that the tidal heating is concentrated in a hypothetical magma ocean that is approximated by a fluid layer.The main difference between the solid and fluid tides is that the tidal deformation of a fluid layer is affected by the Coriolis force, an effect that is negligible in solid tide models.Tyler et al. (2015) shows that fluid-tide models predict different patterns of tidal heating than solid-tide models and can explain Io's heat production over a wide range of plausible parameters.
The existence of a magma ocean on Io was predicted by Peale et al. (1979), shortly before the Voyager 1 mission discovered Io's active volcanism (Smith et al., 1979).At the same time, the mission revealed mountains with elevations of ∼10 km, suggesting that Io must have a cold lithosphere (O'Reilly & Davies, 1981) that is significantly thicker than that proposed by Peale et al. (1979).Although the concept of a mushy magma ocean was supported by geological analysis (Keszthelyi et al., 1999), some scientists remained skeptical and argued that heat transport by melt segregation would lead to rapid cooling, reducing the melt fraction and preventing the formation of a magma ocean (e.g., Moore, 2001).The question of whether Io has a magma ocean was reopened in 2011 when Khurana and co-workers analyzed the magnetometer data collected by the Galileo spacecraft near Io and showed that the data were consistent with the presence of a global conductive layer.Taking into account the electrical properties of partially molten rocks, Khurana et al. (2011) interpreted this layer as a magma ocean with a thickness exceeding 50 km and a rock melt fraction of a few tens of percent.However, this model was challenged by Roth et al. (2017) and Blöcker et al. (2018) who argued that the interaction of the Jovian magnetosphere with Io's plasma environment is a more likely explanation than a magma ocean.Recently, Miyazaki and Stevenson (2022) have explored the steady state of a solid layer with a high melt fraction ("magmatic sponge" with porosity >0.2).They showed that for a wide range of parameters such a layer would be unstable and it would swiftly separate into two phases, leading to the formation of a subsurface magma ocean.Such an ocean would likely contain some amount of crystals but it would behave rheologically as a liquid.The existence of a magma ocean does not contradict the results of Roth et al. (2017) and Blöcker et al. (2018) because the magma layer may be relatively thin (∼1-10 km) and the magnetic induction signal from Io's interior may be weak compared to the magnetic field perturbations caused by the plasma interaction with Io's asymmetric atmosphere.The presence of a magma ocean on Io has recently been supported by new observations by the Juno spacecraft, which show that Io's hot spots are more or less evenly distributed over the surface, but polar volcanoes emit less energy than volcanoes at lower latitudes (Davies et al., 2024).Such a heat flux distribution cannot be explained by tidal dissipation in the deep interior, but it is consistent with the presence of a global magma ocean and/or shallow asthenospheric heating (e.g., Matsuyama et al., 2022;Tyler et al., 2015).
In this study, we investigate the tidal dissipation in Io's hypothetical magma ocean using a new approach based on the solution of the three-dimensional Navier-Stokes equations.Unlike the study of Tyler et al. (2015), where the mechanical coupling between the ocean and the solid parts of the moon was neglected, the flow in the ocean is calculated simultaneously with the deformation of the lithosphere and the sub-oceanic mantle.The resulting maps of tidal dissipation are compared with the geological evidence and the possible role of a magma ocean in Io's thermal evolution is discussed.

Method
The models presented in this paper were obtained by solving the following set of equations: where σ is the incremental stress tensor, ρ is the density, ω is the angular frequency, V t is the tidal potential, V g is the gravitational potential due to the deformation, t is the time, v is the velocity, σ d is the deviatoric part of σ, η is the dynamic viscosity, μ is the shear modulus and • T denotes the transpose of a tensor.Equation 1 is the momentum equation including the Coriolis force ( 2ρω × v) and the time-varying tidal potential (Kaula, 1964), Here, r, θ, and ϕ are the spherical coordinates, P 20 and P 22 are the associated Legendre functions, and Geophysical Research Letters where e is the eccentricity and ω is the angular speed.The effect of the obliquity tides has been neglected (Baland et al., 2012).Equation 2is the continuity equation for an incompressible flow.Finally, Equation 3 is the Maxwell constitutive law for a viscoelastic body expressed in terms of velocity v.In the case of μ → ∞, the equation reduces to a constitutive relation for a Newtonian fluid.We assume that the surface of the moon is free to move, the velocity and traction vectors are continuous at the internal interfaces and the core is in hydrostatic equilibrium.
Equations 1-3 are solved in the time domain using the semi-spectral method developed by Aygün and Čadek (2023b).The spherical harmonic expansions are truncated at degree 20-200 depending on the viscosity of the magma ocean, while the spherical harmonic coefficients are discretized in 800 unevenly spaced radial points.
The radial resolution ranges from 1 to 50 m in the magma ocean and from 100 m to 4 km in the rest of the mantle.Tidal dissipation (heat power per unit volume) is calculated using the formula where P is the rotation period, t 0 is an arbitrary time and the symbol : denotes the Frobenius inner product (σ d ij σ d ij in the Cartesian components).Equation 6can be used to calculate tidal dissipation in both liquid and solid phases, provided that the solid is represented by a Maxwell viscoelastic model (Souček et al., 2016, see Section S5 in Supporting Information S1).The total dissipation or the total heat production, H, is calculated as the integral of h over the volume of the magma ocean.The spatial distribution of dissipation in Io's magma ocean is presented in the form of maps showing the heating h integrated over the thickness of the ocean ("tidal heat flux"), where R m and d are the outer radius and the thickness of the magma ocean, respectively.
The modeling approach used here is different from that used by Tyler et al. (2015) in two main respects.First, the tidal response of the magma ocean is calculated by solving the three-dimensional (3D) Navier-Stokes equations, while Tyler et al. (2015) used the Laplace tidal equations (LTE) where the tidal flow is described as a barotropic two-dimensional sheet flow.Second, the response of Io to tidal loading is calculated not only in the magma ocean but also in the lithosphere and the sub-oceanic mantle, which allows us to precisely quantify the mechanical and gravitational coupling between the three layers.
Our method is also more general than the methods used in studies investigating the tidal response of water ocean worlds (e.g., Beuthe, 2016;Matsuyama et al., 2018;Rekier et al., 2019;Rovira-Navarro et al., 2019).To couple the flow in the ocean with the deformation of the crust, Beuthe (2016) and Matsuyama et al. (2018) proposed to solve the LTE together with the equations governing the viscoelastic deformation of the overlying shell.Although their method is similar at first glance to our approach, it differs from it in that the boundaries of the ocean are treated as free-slip surfaces and the flow velocity does not change with radius.As recently shown by Aygün and Čadek (2023b), the method by Beuthe (2016) 2019) determine the dissipation rate in the ocean by using the 3D Navier-Stokes equations, but assume that the deformation of the crust is not affected by the flow in the ocean and can therefore be imposed as a boundary condition at the surface of the ocean.However, this assumption is valid only if the thickness of the ocean layer is greater than about 0.01R m ≈ 15 km, that is, outside the thickness range considered in the present study (Aygün & Čadek, 2023b).
The density structure of Io is chosen to match the total mass (8.9319 • 10 22 kg) and MoI factor (0.37685, Anderson et al., 2001) and to satisfy constraints on its material composition (Anderson et al., 1996).The upper boundary of the magma ocean is set to a depth of 30 km and the thickness of the ocean is varied from 100 m to 10 km.We assume that the magma behaves as a Newtonian liquid and its viscosity is constant throughout the ocean.The

Geophysical Research Letters
10.1029/2023GL107869 viscosity of the magma ranges from 100 Pa s (hot mafic magma) to 10 7 Pa s (low-temperature magma with solid crystals suspended in the liquid phase, see, e.g., Philpotts and Ague (2009)).For simplicity, we assume that the density of the magma is the same as the density of the mantle (for details, see Table S1 in Supporting Information S1).
The lithosphere and the mantle below the ocean are assumed to behave as a Maxwell viscoelastic solid.Since the main focus of this study is to examine tidal heating in a hypothetical magma ocean, the material parameters of the lithosphere and the sub-ocean mantle are chosen so that the tidal heat production outside the magma ocean is much smaller than Io's current heat output.

Results
As illustrated in Figure 1, the tidal flow in the magma ocean produces a wide variety of heating patterns and even small changes in ocean thickness can lead to order of magnitude changes in the total heat production.The heat flux distributions are symmetric about the equator and most of them, but not all, also about the tidal axis.The heat flux patterns are dominated by dissipation at low latitudes, typical of fluid models.The highest heat production (>10 4 TW) is found in the case where tidal heating is concentrated in an equatorial zone at latitudes below 30°.This equatorial zone is clearly separated from the low-dissipation regions at higher latitudes and its position shows a remarkable correspondence with Io's yellow bright plains (Williams et al., 2011, , 2015), except for a study including the moon-moon tidal interaction (Hay et al., 2020, based on the communication with Hamish Hay).Although the zonal distribution of dissipation is obtained for models with a thin ocean (d ≈ 1 km), the velocity field has a strong radial component (v r /v θ ≈ v r /v ϕ ≈ 0.2) that cannot be found by solving the LTE.The dependence of tidal heating on the thickness of the ocean and the viscosity of magma is shown in Figure 2a.Inspection of the figure shows that the heat power of ∼100 TW (Io's observed heat output) is exceeded over a wide range of magma viscosities and ocean thicknesses, with maximum values achieved for d ≈ 1 km and η between 10 2 and 10 4 Pa s.When d is significantly larger than 1 km, the 100 TW limit can only be achieved for η > 10 4 Pa s.The results in Figure 2a are qualitatively similar to those obtained from the solution of the LTE (Figure 4a in Tyler et al. (2015)) but they differ in three minor respects: First, the ocean thickness for which the maximum dissipation is reached is about twice the value predicted by Tyler et al. (2015), second, no dissipation maxima are found for d < 1 km and third, the dissipation in the ocean vanishes when d → 0. As  S1 in Davies et al. (2024).

Geophysical Research Letters
10.1029/2023GL107869 illustrated in Figures 2b and 2c, the models of magma ocean tides cannot explain the observed distribution of Io's hot spots.Most of the models satisfactorily explain the hot spots at latitudes below 50°but do not account for hot spots in the polar regions.
The tidal heating models discussed above were obtained under the assumption that Io's solid layers are only weakly dissipative.In Supporting Information S1, we also show the results obtained for the case where the magma ocean is underlain by a 100-km thick, low viscosity "magmatic sponge" layer.Comparison of Figure 1 with Figure S1 in Supporting Information S1 indicates that the total dissipation in the ocean is only weakly sensitive to the viscosity of the sub-oceanic mantle.
The origin and the distribution of tidal heating within Io's interior remain a subject of debate.The question could, in principle, be answered by a dedicated mission making close flybys of Io and providing new information about Io's gravity signature.Bierson and Nimmo (2016) have demonstrated that the tidal Love number k 2 of Io is highly sensitive to the presence of a fluid layer beneath the surface.The value of k 2 should be about 0.5 for Io with a fluid magma ocean and 0.1 if Io's mantle behaves as a solid (see also Kervazo et al., 2022).The first value should be regarded as a rough estimate because it was obtained under the assumption that the dynamic effect of the tidal flow in the ocean can be neglected.
In the absence of a fluid magma ocean, the gravitational response of Io to tidal forcing can be expressed as , where a is Io's radius, V t is the tidal potential (Equation 4) and Δt 2 is the time lag.The response is described by only two parameters, k 2 and Δt 2 .If Io has a fluid magma ocean, the tidal deformation generates a degree-2 flow in the ocean, which is further modulated by the Coriolis effect.Depending on the parameters of the model, the resulting flow can deform the surface and internal density interfaces, generating a gravitational signal that is much more complex than in the case of a solid body.The gravitational response at degree 2 varies with the harmonic order and is described by 6 parameters, k 20 , k c 22 , k s 22 , Δt 20 , Δt c 22 , and Δt s 22 , which can be determined from the following equations (cf.Equation 4).
V g,20 (a,t) = k 20 V t,20 (a,t Δt 20 ), (8) where V g,20 , V c g,22 , and V s g,22 are the coefficients of the gravitational potential induced by tidal potential coefficients V t,20 , V c t,22 and V s t,22 , respectively.
In the case of solid-body tides, k 20 = k c 22 = k s 22 = k 2 and Δt 20 = Δt c 22 = Δt s 22 = Δt 2 .On the other hand, if Io has a magma ocean, the degree-2 Love numbers and time delays can significantly vary with the order (m = 0, 2).While k 20 increases with the increasing ocean thickness, reaching a maximum of 0.85 for d = 10 km and η = 3 • 10 6 Pa s (Figure 3a), k c 22 and k s 22 (Figures 3b and 3c) are strongly affected by the Coriolis effect and correlate with the total heat production (cf. Figure 2a).The maximum values of k c 22 and k s 22 are about 10 times greater than the maximum value of k 20 .Large differences are also found between Δt 20 on one side and Δt c 22 and Δt s 22 on the other.The values of the tidal Love numbers and time lags corresponding to Io's current dissipative power (≈100 TW) are shown by the blue lines.If the ocean thickness is small and/or the magma viscosity is high, k 20 ≈ k c 22 ≈ k s 22 < 0.1 and Δt 20 ≈ Δt c 22 ≈ Δt s 22 < 2 hr, suggesting that, in this case, the presence of a magma ocean has little effect on the large-scale deformation of the lithosphere.

Discussion and Conclusions
Our results confirm the conclusion of Tyler et al. (2015) that the tidal heating in a fluid magma ocean can explain Io's observed heat production over a broad range of magma viscosities and ocean thicknesses.Compared to the studies using the LTE method (Matsuyama et al., 2022;Tyler et al., 2015), the solution of the Navier-Stokes equations is characterized by a greater variety of heat flux patterns, which vary depending on the parameters of the ocean.The highest heat production is found for models with a thin (∼1 km thick) fluid layer and low (<10 4 Pa s) viscosity.The total heat production decreases with the ocean thickness (Figure 2a) suggesting that the thickness of the ocean (if it exists) is likely to be less than ∼10 km.An extrapolation of the results indicates that a thick (∼50 km) magma ocean is only weakly dissipative unless the viscosity of magma is higher than 10 8 Pa s.
It is tempting to speculate that the formation of a magma ocean may have been triggered by a temporary increase in Io's eccentricity (Hussmann & Spohn, 2004), resulting in the enhancement of tidal heating and an increase of porosity in a magmatic sponge.The subsequent formation of a shallow (∼0.1-1 km thick) magma ocean (Miyazaki & Stevenson, 2022) led to a sharp increase of the tidal dissipation rate in Io's interior, especially at low (<30°) latitudes (Figure 2c).Due to a positive feedback between temperature and tidal dissipation, the magma production rate was faster than the rate of magma extraction, leading to an increase in the magma ocean thickness and a gradual change in the dissipation pattern (Figure 1).Since dissipation in the ocean decreases with the ocean thickness (Figure 2a), the total heat production gradually stabilized at the present level.A subsequent decrease in eccentricity did not necessarily lead to the disappearance of the ocean because once the ocean is formed, the heat production in the magma layer is high enough to maintain it in a liquid state.
The question of whether present-day Io has a magma ocean or not is difficult to answer.It is usually assumed that the Io's volcanic activity should be correlated with the distribution of tidal heating.Magma ocean models (Figure 1) predict enhanced tidal heating at low latitudes and low tidal heating in the polar regions.However, new data from the Juno mission indicate that the density of hot spots does not decrease toward the poles (Davies et al., 2024;Zambon et al., 2022), in contrast to previous data sets where most of the hot spots were located at low latitudes (e.g., Davies et al., 2015).The low correlation between the hot spots and the predicted tidal heating at high latitudes does not necessarily mean that there is no magma ocean on Io.If the ocean is global (i.e., if the magma is stored in a global continuous reservoir), the melt is present everywhere beneath the lithosphere, and the amount of magma that gets to the surface then depends on the local conditions, rather than on the distribution of tidal heating.In other words, volcanic eruptions can occur at any location where conditions in the lithosphere are favorable for the ascent of magma (Crisp, 1984;de Kleer et al., 2019;Jaeger et al., 2003).New observations of Io's thermal emission by the Juno spacecraft Jovian Infrared Auroral Mapper suggest that polar hot spots emit less energy than hot spots at lower latitudes (Davies et al., 2024).This finding is consistent with the presence of a magma ocean, although enhanced tidal heating at low latitudes can also be predicted by models with a dissipative asthenosphere (e.g., Kervazo et al., 2021;Matsuyama et al., 2022;Tyler et al., 2015).Another possible explanation for hot spots in polar regions is the tidal heating in the solid sub-oceanic mantle (Tyler et al., 2015).Unlike dissipation in a liquid ocean, which is concentrated at low latitudes, dissipation in the deep mantle mainly occurs near the poles (e.g., Beuthe, 2013) and, therefore, it can compensate for the decrease in ocean tidal heating at high latitudes.If this is the case, the current distribution of hot spots on Io results from the combined effect of solid and fluid tidal dissipation (see also Figure S2 in Supporting Information S1).However, the question is whether the heat flux obtained in this way is physically meaningful because the temperature field in Io's deep mantle is likely to be affected by convection (Tackley, 2001;Tackley et al., 2001) and the heat transfer between the solid mantle and the liquid magma ocean can be modulated by the solid-liquid phase transition (e.g., Labrosse et al., 2018).As shown by Steinke et al. (2020), mantle convection tends to significantly reduce the lateral temperature signal, depending on the thickness of the convective layer, the viscosity and the mode of heat transfer.Finally, it is possible that the magma ocean has a variable thickness and may even be absent in some areas.Dissipative behavior of such an ocean is difficult to predict but it is likely that tidal heating would strongly vary laterally and would be affected by regional resonance effects.Bierson and Nimmo (2016) suggested that the presence of a fluid magma ocean on Io could be detected by measuring the tidal Love numbers.Our results indicate that if tidal heating preferentially occurs in a magma ocean and the total heat production is about 100 TW, then the degree-2 Love numbers are either less than 0.1 if d < 2 km (i.e., about the same as in the case of solid tides) or greater than 0.7 if d > 2 km.While in the former case, k c 22 ≈ k s 22 ≈ k 20 , in the latter case, k c 22 and k s 22 are twice as big as k 20 .In both cases the time lag is less than 2 hr.The fact that the tidal Love numbers are not sensitive to the presence of a liquid magma ocean if the ocean thickness is small needs to be taken into account in future analyses of Io's gravity signature.

•
Due to the Coriolis effect, the degree-2Love numbers for models with a magma ocean can depend on the harmonic order• The tidal Love numbers are not sensitive to the presence of a fluid magma ocean if the thickness of the fluid layer is less than 2 km Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.Distribution of tidal dissipation in the magma ocean, Equation7, evaluated for different thicknesses and viscosities of the ocean (Mollweide equal-area projection centered on 180°W longitude).The numbers above each map represent the total heat production of the magma ocean in TW.

Figure 2 .
Figure 2. (a) Total heat production in the magma ocean as a function of ocean thickness and magma viscosity.The contours plotted in blue correspond to Io's present-day heat production (100 TW).(b) Heat flux map calculated for η = 10 5 Pa s and d = 10 km.While the total heat production (96 TW) is compatible with the observation, the dissipation pattern shows no relationship with the current distribution of hot spots (blue circles), as only a half of the hot spots are located in a region of higher-thenaverage heat flux.(c) Heat flux map corresponding to the model with the highest dissipation (2.1 • 10 4 TW).In this case, only about 35% of hot spots are located in high dissipation regions.(d) Geological map of Io(Williams et al., 2011).For the sake of comparison with the geological map, the tidal heat flux in panels (b) and (c) is shown in the equidistant cylindrical projection.The color scale is the same as in Figure1.The locations of hot spots are taken from TableS1inDavies et al. (2024).

Figure 3 .
Figure 3. Tidal Love numbers (a-c) and time lags (d-f) as functions of the ocean thickness and magma viscosity.The blue lines correspond to models with a heat output of 100 TW.The contour interval in panels (a-c) is 0.1 if the Love number is less than 1.5 and 0.5 if the Love number is greater than or equal to 1.5.The contour interval in panels (d-f) is 2 hr.