Mathilde Kervazo

Gabriel Tobie

Gaël Choblet

Caroline Dumoulin

Marie Běhounková

Inferring Io's interior from tidal monitoring

Keywords: Io, Tides, solid body, Satellites, dynamics

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The volcanically active moon of Jupiter, Io, emits a spectacular amount of volcanic plumes and lavas, corresponding to an average heat flux esti-mated to be of the order of 2.5 W/m 2 (e.g. Veeder et al., 1994;[START_REF] Lainey | Strong tidal dissipation in io and jupiter from astrometric observations[END_REF]. This means surface heat flux is almost thirty times larger than the one estimated for the Earth (Turcotte and Schubert, 2002). However, in contrast to terrestrial planets with recent or active silicate volcanism, Io's present heat output cannot be explained by radiogenic sources and remnant cooling, but requires an additional source, most likely extreme dissipation of tidal energy owing to its relatively close distance to Jupiter (e.g. [START_REF] Peale | Melting of io by tidal dissipation[END_REF][START_REF] Segatz | Tidal dissipation, surface heat flow, and figure of viscoelastic models of io[END_REF][START_REF] Ross | Internal structure of io and the global distribution of its topography[END_REF][START_REF] Beuthe | Spatial patterns of tidal heating[END_REF][START_REF] Bierson | A test for io's magma ocean: Modeling tidal dissipation with a partially molten mantle[END_REF]Steinke et al., 2020a). The heat production rate, estimated to more than 100 times the radiogenic heat rate, is sufficient to melt portions of Io's interior, feeding hundreds of continually erupting volcanoes (e.g. [START_REF] Lopes-Gautier | Active volcanism on io: Global distribution and variations in activity[END_REF][START_REF] Lopes | Lava lakes on io: Observations of io's volcanic activity from galileo nims during the 2001 fly-bys[END_REF][START_REF] Davies | Map of io's volcanic heat flow[END_REF][START_REF] Cantrall | Variability and geologic associations of volcanic activity on io in 2001-2016[END_REF][START_REF] Mura | Infrared observations of io from juno[END_REF].

While the idea of substantial melting in Io's mantle is broadly consistent with the interpretation of several observations, the localization of melt production in the interior and the way it is redistributed to the surface remains debated and poorly constrained. High eruption temperatures of Io's silicate volcanism indicate that the interior, at least the upper mantle beneath the crust, is largely molten [START_REF] Keszthelyi | Revisiting the hypothesis of a mushy global magma ocean in io[END_REF][START_REF] Keszthelyi | A post-galileo view of io's interior[END_REF], with melt fraction estimated to about 20%-30% [START_REF] Keszthelyi | New estimates for io eruption temperatures: Implications for the interior[END_REF]. High concentration of melts below Io's near surface is also consistent with models describing the release of interior heat via melt extraction [START_REF] Moore | Tidal heating and convection in io[END_REF][START_REF] Bierson | A test for io's magma ocean: Modeling tidal dissipation with a partially molten mantle[END_REF]Steinke et al., 2020a;[START_REF] Spencer | Magmatic intrusions control io's crustal thickness[END_REF], but a deep melt source cannot be ruled out (e.g [START_REF] Monnereau | Is io's mantle really molten?[END_REF]. Furthermore, the presence of a high melt fraction (∼ 20 -30%) layer has been argued to be consistent with Galileo magnetic induction measurements [START_REF] Khurana | Evidence of a global magma ocean in io's interior[END_REF], although the measured perturbations could also result from plasma interactions with the atmosphere [START_REF] Blöcker | Mhd modeling of the plasma interaction with io's asymmetric atmosphere[END_REF]. Higher melt fractions (> 30%) and the possible existence of a fully liquid magma ocean cannot be ruled out with the available Galileo magnetic induction measurements but seems to be at odds with the phase of the auroral spot oscillations [START_REF] Roth | Constraints on io's interior from auroral spot oscillations[END_REF]. In summary, even if there is some kind of consensus regarding the existence of a partially molten zone in Io's interior, the depth and extent of this zone still remains unconstrained.

The total amount of heat produced by tidal friction and its distribution in the interior is intimately linked to the structure and thermal state of Io's interior, especially the distribution of temperature and melt fraction [START_REF] Bierson | A test for io's magma ocean: Modeling tidal dissipation with a partially molten mantle[END_REF]Steinke et al., 2020a;[START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF]. The presence of melts is expected to strongly impact the viscoelastic properties of rocks (e.g. [START_REF] Budiansky | Elastic moduli of a cracked solid[END_REF][START_REF] Mavko | Velocity and attenuation in partially molten rocks[END_REF][START_REF] Takei | Constitutive mechanical relations of solid-liquid composites in terms of grain-boundary contiguity[END_REF]Hirth and Kohlstedt, 1995a,b;[START_REF] Kohlstedt | Rheology of partially molten rocks[END_REF][START_REF] Scott | The effect of large melt fraction on the deformation behavior of peridotite[END_REF], and, in return, the amplification of tidal heating in molten areas can create a positive feedback R7 on melt generation. Describing the mechanical response of partially molten rocks for a wide range of melt fraction is thus essential to correctly describe the tidal friction in Io's hot interior. Most of previous studies varied the elastic and viscous properties of the mantle and of the partially molten layer in an arbitrary manner in order to match the observed heat output (e.g. [START_REF] Segatz | Tidal dissipation, surface heat flow, and figure of viscoelastic models of io[END_REF][START_REF] Renaud | Increased tidal dissipation using advanced rheological models: Implications for io and tidally active exoplanets[END_REF]Steinke et al., 2020a). To our knowledge, the only models accounting for the combined evolution of elastic and viscous properties as a function of melt fraction were proposed by [START_REF] Bierson | A test for io's magma ocean: Modeling tidal dissipation with a partially molten mantle[END_REF] and [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF]. In the latter study, we showed that bulk dissipation, neglected in all previous studies on Io, can contribute significantly to the tidal heat budget and impact the tidal heating pattern, especially when a thin and highly molten layer is considered. Beyond Io, better understanding the link between melt production and tidal dissipation has implications for the eraly history of rocky planets in our solar system and in extrasolar systems (e.g [START_REF] Moore | Heat-pipe planets[END_REF][START_REF] Renaud | Increased tidal dissipation using advanced rheological models: Implications for io and tidally active exoplanets[END_REF]. R8 Clues to the nature of tidal dissipation inside Io have long been sought from surface heat fluxes (Veeder et al., 2012) and volcano distributions [START_REF] Ross | Internal structure of io and the global distribution of its topography[END_REF][START_REF] Kirchoff | Global distribution of volcanic centers and mountains on io: Control by asthenospheric heating and implications for mountain formation[END_REF][START_REF] Hamilton | Spatial distribution of volcanoes on io: Implications for tidal heating and magma ascent[END_REF][START_REF] Rathbun | The global distribution of active ionian volcanoes and implications for tidal heating models[END_REF], an exploration complicated by the transient nature of volcanic events as well as the apparent dearth of polar observations [START_REF] Rathbun | The global distribution of active ionian volcanoes and implications for tidal heating models[END_REF][START_REF] Mura | Infrared observations of io from juno[END_REF]. Lithospheric thickness and topography have also been proposed as useful proxy for long-timescale heat flux [START_REF] Ross | Internal structure of io and the global distribution of its topography[END_REF]Steinke et al., 2020a;[START_REF] Spencer | Tidal controls on the lithospheric thickness and topography of io from magmatic segregation and volcanism modelling[END_REF]. Such observations provide some insights on the lateral variations of heat production and mechanical properties, but cannot be used to really constrain the internal structure. A complementary way to access key information about the average internal structure of planetary bodies is the measurement of libration (Van Hoolst et al., 2020) and tidal deformation from spacecraft tracking, altimetry and high-resolution imaging (e.g. [START_REF] Mazarico | Simulated recovery of europa's global shape and tidal love numbers from altimetry and radio tracking during a dedicated flyby tour[END_REF][START_REF] Dumoulin | Tidal constraints on the interior of venus[END_REF][START_REF] Steinbrügge | Viscoelastic tides of mercury and the determination of its inner core size[END_REF][START_REF] Park | Advanced pointing imaging camera (apic) for planetary science and mission opportunities[END_REF]. Such measurements allow the determination of the Love numbers k 2 , h 2 and l 2 , which quantify the gravitational potential and the surface radial and horizontal displacement associated to tidal motions, respectively. Concept of R1 future missions dedicated to Io, as the 'Io Volcano Observer' (IVO) [START_REF] Mcewen | Follow the heat: Io volcano observer[END_REF], one of the Discovery finalists, although not selected by NASA R1 , may provide the first estimate of these quantities, essential to better understand how heat and magma is generated and transported to the surface.

In this context, the goal of the present study is to compute the tidal response of Io's interior for various distributions of internal melt consistent with the estimated heat output and to predict the expected gravimetric, altimetric and geodetic signatures and heat flow pattern that might be compared to future measurements. A coherent melt profile between the subsurface partially molten layer and the underlying deep mantle is considered following petrological and two-phase flow arguments [START_REF] Moore | The thermal state of io[END_REF][START_REF] Keszthelyi | New estimates for io eruption temperatures: Implications for the interior[END_REF][START_REF] Bierson | A test for io's magma ocean: Modeling tidal dissipation with a partially molten mantle[END_REF][START_REF] Spencer | Magmatic intrusions control io's crustal thickness[END_REF] R9 . We vary the melt fraction R10 from a few percent R10 melt fraction to values just above the so-called rheological critical melt fraction (> 20 -30%). This corresponds to a sharp transition from the solid-dominated behavior to the liquid-dominated behavior. The rheological parameterization developed by [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF] is used to consistently take into account the role of melt fraction on the elastic and viscous parameters of Io's partially molten interior. The computation of the tidal response described by both shear and bulk viscoelastic rheology is performed for a wide range of internal parameters and is analyzed in terms of predicted tidal Love number k 2 , h 2 and l 2 as well as in terms of polar heat flux and heat flow patterns.

Method

2.1. Rheological model and properties of Io's interior

Existing geophysical constraints

The main observational constraints about Io's interior come from the gravity measurements [START_REF] Anderson | Io's gravity field and interior structure[END_REF], magnetic induction signals [START_REF] Khurana | Evidence of a global magma ocean in io's interior[END_REF] as well as heat budget assessment from infrared remote sensing (e.g. Veeder et al., 1994;[START_REF] Spencer | Discovery of gaseous s2 in io's pele plume[END_REF] and astrometric measurements [START_REF] Lainey | Strong tidal dissipation in io and jupiter from astrometric observations[END_REF]. Galileo gravity data put relatively good constraints on the mean density (3527.8 kg m -3 ) and the moment of inertia (0.37685 R2 , [START_REF] Anderson | Io's gravity field and interior structure[END_REF]. Magnetic induction signals [START_REF] Khurana | Evidence of a global magma ocean in io's interior[END_REF] are consistent with a partially molten layer beneath the crust, but do not really constrain its melt content or depth. Based on various estimates, the average heat flux is evaluated to 2.24 ± 0.45 Wm -2 (e.g. Veeder et al., 1994;[START_REF] Spencer | Discovery of gaseous s2 in io's pele plume[END_REF][START_REF] Lainey | Strong tidal dissipation in io and jupiter from astrometric observations[END_REF], corresponding to a total power ranging between 65 and 125 TW. In the following we choose 100 TW as a reference value and consider a range between 65 and 125 TW when computing the total heat production in Io's interior (Table 1).

Interior model parameters

Based on geophysical constraints, the interior model considered here consists of (from surface to center) a silicate crust, a melt-rich upper mantle layer, called asthenosphere hereafter, a melt-poor silicate deep mantle and a liquid metallic core (Figure 1a). For simplicity, the density of each layer is assumed to be uniform. In the reference model (Table 1), we set the crustal thickness to 30 km (also corresponding to the asthenospheric depth in our model) with a density of 3000 kg/m 3 . The core density is set to 5165 kg m -3 , corresponding to an eutectic Fe-FeS core composition (e.g. [START_REF] Segatz | Tidal dissipation, surface heat flow, and figure of viscoelastic models of io[END_REF][START_REF] Anderson | Galileo gravity results and the internal structure of io[END_REF]. This assumption results in a core radius of 955 km and a mantle density of 3263 kg/m 3 in order to satisfy the average density and moment of inertia. To assess the sensitivity of tidal deformation to internal structure parameters, we explore values as large as 90 km for the crustal thickness and 8000 kg/m 3 for the core density, corresponding to a pure iron core. For the dense core end-member, the size of the core reduces to 665 km.

All internal layers (except the liquid iron core) are assumed to behave as viscoelastic solid media at tidal frequencies. To take into account the effect of melt on the viscoelastic response, a specific rheological model including the description of melt, inspired from the law proposed by [START_REF] Costa | A model for the rheology of particle-bearing suspensions and partially molten rocks[END_REF], is considered (see [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF], and further described in section 2.1.3).

To account for the unknown composition and thermal state of Io's silicate mantle, we consider viscosity values η sol for the solid rock matrix ranging from 10 16 to 10 20 Pa s. A viscosity of the mantle as low as 10 16 Pa s is chosen as a lower bound, allowing to produce 100 TW in Io's interior without invoking the presence of partial melt in our model. We note that recent work shows a drastic drop in viscosity when a very small amount of melt is present in silicate rocks [START_REF] Holtzman | Questions on the existence, persistence, and mechanical effects of a very small melt fraction in the asthenosphere[END_REF]. It is therefore wise to study the effect of a very low viscosity in a mantle as hot as that of Io. For the reference value, we choose the typical value expected for the solid-state viscosity for dry olivine-dominated rocks near their melting point, 10 19 Pa s (e.g. [START_REF] Karato | Rheology of the upper mantle: A synthesis[END_REF], and consider a value of 10 20 Pa s, as an upper bound. For the reference model, we choose a value of the reference shear modulus µ mant = 60 GPa, justified for olivine [START_REF] Tan | High-temperature viscoelasticity of fine-grained polycrystalline olivine[END_REF][START_REF] Jackson | Grainsize-sensitive viscoelastic relaxation in olivine: Towards a robust laboratory-based model for seismological application[END_REF].

Based on these experimental studies, the shear modulus of olivine is assumed to range typically between 50 and 70 GPa (Table 1).

We divide the crust into a top rigid crust (with a shear modulus µ crust = 65 GPa and a high viscosity η crust = 10 23 Pa s) above a less rigid and less viscous crust which rheological properties are those of the mantle (that is µ mant and η sol ). The thickness of the elastic part of the crust is set to 2/3 of the total crustal thickness, in an arbitrary manner. Sensitivity tests showed that the relative thickness of the elastic/viscoelastic parts of the crust has a negligible influence on the results reported hereafter. The bulk modulus K for the crust, mantle and core is set to 200 GPa. Note that the effect of temperature and pressure change with depth on the silicate bulk modulus is only moderate: values typically range between 150 and 250 GPa for Io' s mantle pressure conditions (e.g. [START_REF] Jackson | Elasticity, composition and temperature of the earth's lower mantle: a reappraisal[END_REF]. Varying the bulk modulus of the solid phase has a neglibile effect on the results [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF] and therefore a single value of K=200 GPa is considered for both the crust and the mantle. The iron core is assumed to be fully liquid and inviscid, and therefore its shear modulus and viscosity are set to zero.

Melt-based rheological model

We consider that the whole mantle (either the deep mantle or the asthenosphere) can undergo partial melting, in agreement with the simulations of [START_REF] Spencer | Magmatic intrusions control io's crustal thickness[END_REF]. Their study of Io's internal dynamics in a two-phase flow framework showed that a higher melt fraction region, corresponding to a decompaction boundary layer owing to maximal liquid pressure beneath the crust, is expected to surmount a deeper mantle with less melt present.

Compared to their simulations, our assumed melt profiles are simplified, with melt fraction considered uniform in the two distinct layers corresponding to the asthenosphere and the rest of the mantle underneath (Figure 1b). Melt fractions are referred to φ ast and φ mant respectively, with φ mant related to φ ast by a factor typically ranging between 1/10 (in magenta, Figure 1b) and 1/3 (in black, Figure 1b) following the results of [START_REF] Spencer | Magmatic intrusions control io's crustal thickness[END_REF].

While magnetic induction measurements as well as petrological indices favor significant amounts of near-surface melting (about 20-30%, e.g. [START_REF] Keszthelyi | New estimates for io eruption temperatures: Implications for the interior[END_REF][START_REF] Khurana | Evidence of a global magma ocean in io's interior[END_REF], the amount of melt in the deep mantle is not constrained. We therefore investigate a relatively large range for φ ast (from 0 to 35%, Table 1) and vary the mantle melt fraction according to the asthenospheric melt fraction. We set the thickness of the melt-rich asthenosphere to 100 km in our reference model, and vary it from 50 to 200 km (Table 1) which is in the range of what is commonly assumed in the literature for this layer (e.g. [START_REF] Segatz | Tidal dissipation, surface heat flow, and figure of viscoelastic models of io[END_REF][START_REF] Khurana | Evidence of a global magma ocean in io's interior[END_REF]Steinke et al., 2020a).

[Table 1

about here.]

We developed a synthetic law, inspired from the melt-based viscosity law of [START_REF] Costa | A model for the rheology of particle-bearing suspensions and partially molten rocks[END_REF], to account for the melt dependence of elastic shear and bulk moduli, and shear and bulk viscosity, consistent with existing experimental and theoretical constraints (Figure 1c, see [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF] for more details). R5

The rheological law is the following:

•(φ) = • l 1 + Θ δ [1 -F (Θ, ξ, γ)] B(1-φ * ) , (1)
where rheological parameter • is either the shear viscosity η, the shear modulus µ or the bulk modulus K. Two auxiliary functions are introduced:

Θ = (1 -φ)/(1 -φ * ), (2)
F = (1 -ξ) erf √ π 2(1 -ξ) Θ(1 + Θ γ ) . (3)
Besides parameter B, the Einstein coefficient (set to 2.5), all other parameters are tuned to reproduce the available constraints on the specific rheological parameter, from the solid-state endmember • s to the fully liquid state endmember • l . The values of δ, ξ, γ and φ * also depend on the specific rheological parameters considered. These are listed in Table 2. For the bulk viscosity, the simplified relationship ζ = η φ is considered (see [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF], for more details).

This formulation allows us to predict both viscous and elastic properties on the entire range of melt fraction, below or above the critical melt fraction, φ c , using a single expression. The latter reproduces the rheological transition associated with a sharp decrease from the solid behaviour to the liquid behaviour (e.g. [START_REF] Renner | On the rheologically critical melt fraction[END_REF]. For φ < φ c , the parameters used in the rheological law (Eq.1) have been derived to reproduce experimental and theoretical constraints published in the literature on the viscous and elastic properties of terrestrial partially molten rocks (see Kervazo et al., 2021, Appendix A). The formulation of [START_REF] Costa | A model for the rheology of particle-bearing suspensions and partially molten rocks[END_REF] was derived to reproduce the shear viscosity drop occurring at elevated melt fraction, which is relatively well documented in the literature. For shear and bulk elastic moduli as well as bulk viscosity, for which no experimental constraints exist at elevated melt fraction, we assume the same kind of behaviour when φ > φ c R5 . Note that in [START_REF] Costa | A model for the rheology of particle-bearing suspensions and partially molten rocks[END_REF] the viscosity is given as a function of solid volume fraction (φ S ), while here we use the melt fraction (φ S = 1 -φ). In Table 3 of [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF], the last parameter was incorrectly labeled φ * , whereas

it should read (1 -φ * ). R *
In this approach, the rheological parameters, µ, K, η, ζ are considered as the effective parameters of the melt-rock mixture, and the viscoelastic deformation of this mixture is treated as a single media. This does not consider dissipation associated with differential motions between the liquid melt and the solid matrix, usually called melt squirt. Efficient dissipation due to this process is expected to occur when the fluid phase oscillates in the solid matrix at very high frequencies, several orders of magnitude larger than tidal frequencies (e.g. [START_REF] Hammond | Upper mantle seismic wave attenuation: Effects of realistic partial melt distribution[END_REF][START_REF] Faul | Shear wave attenuation and dispersion in melt-bearing olivine polycrystals: 2. microstructural interpretation and seismological implications[END_REF][START_REF] Carcione | Differential form and numerical implementation of biot's poroelasticity equations with squirt dissipation[END_REF]. R6

[Table 2 about here.]

The rheological critical melt fraction φ c is set to 30 % (red lines in Figures 1 b andc). This value marks the transition between the solid-dominated behavior and the liquid-dominated behavior. The assumed value of φ c has no significant effect on the results, what matters is the considered value of the melt fraction relative to this critical value [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF], i.e. whether the melt fraction is above or below φ c . In the following, the results will thus be presented relative to φ ast /φ c . In the context of the present study, focused on solid-body tides, we consider melt fractions up to φ c and slightly above in order to mimic the significant drop in strength, covering the whole range of mantle and asthenospheric partial melting needed to match Io's heat budget.

We note that what we call asthenosphere here is a layer substantially less viscous than the deeper regions (and that the asthenosphere of the Earth) but still considerably more viscous than a fully liquid magma ocean.

[Figure 1 about here.]

Computation of viscoelastic tidal deformation

The viscoelastic deformation of Io under the action of periodic tidal forces is computed following the method of Tobie et al. (2005), adapted to the particular case of Io where the viscoelatic response accounts for shear and bulk dissipative processes in partially molten layers in [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF].

The only difference relative to the approach of Tobie et al. (2005) is the consideration of a non-zero imaginary part for the complex bulk modulus, related to a bulk viscosity using a Maxwell viscoelastic model (see Kervazo et al., 2021, for details).

Using the density profile and rheological properties described in Section 2.1, the Poisson equation and the equation of motions are solved for small perturbations in the frequency domain assuming a compressible viscoelastic media using the problem formulation based on radial functions as defined by [START_REF] Takeuchi | Seismic surface waves[END_REF] for the equivalent elastic problem.

The complex Love numbers, k 2 , h 2 and l 2 , characterizing the potential perturbation and surface radial and tangential displacements, respectively, are computed by integrating the radial functions associated with the radial and tangential displacements (y 1 and y 3 , respectively), the radial and tangential stresses (y 2 and y 4 ), and the gravitational potential (y 5 ), as defined by [START_REF] Takeuchi | Seismic surface waves[END_REF] and adapted to the viscoelastic case (Tobie et al., 2005;[START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF]. The full set of equations, the boundary conditions (center, liquid-solid interface, surface) as well as the numerical scheme to solve them in detail are provided in the appendix of [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF].

Following the approach proposed in [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF], two different viscoelastic models are used for the solid layers depending on the degree of partial melting. For partially molten rocks with high melt fraction (>20%), the estimated Maxwell time approaches Io's tidal period. A Maxwell model thus provides reasonable estimate of the dissipation rate when partial melting occurs, which we consider. In our previous study, we showed that bulk dissipation becomes significant for melt fraction exceeding the critical melt fraction and must be consistently considered when computing the tidal response. In this case, both bulk and shear components of the viscous and elastic parameters (respectively shear and bulk viscosity η and ζ and shear and bulk modulus µ and K), affected by the presence of melt (see Section 2.1.3), are taken into account when computing tidal deformation. In our model, this applies to the asthenosphere.

For rocks with low to moderate melt fraction (φ < 20 %), bulk dissipation becomes negligible, but shear dissipation using the Maxwell model is underestimated. We thus use the more appropriate Andrade model [START_REF] Castillo-Rogez | The tidal history of iapetus: Spin dynamics in the light of a refined dissipation model[END_REF][START_REF] Efroimsky | Tidal dissipation compared to seismic dissipation: In small bodies, earths, and super-earths[END_REF][START_REF] Bierson | A test for io's magma ocean: Modeling tidal dissipation with a partially molten mantle[END_REF] for which the magnitude and time decay of anelastic deformation is described by the shear modulus, µ, the shear viscosity η and two parameters α and β, characterizing the transient viscoelastic response. The typical value of the parameter α required to explain the Q factor of the Earth's mantle [START_REF] Ray | Constraints on energy dissipation in the earth's body tide from satellite tracking and altimetry[END_REF] is estimated between 0.23 and 0.28 [START_REF] Tobie | Tidal response of rocky and ice-rich exoplanets[END_REF]. Varying the Andrade parameter α from 0.2 to 0.3 has only a small effect for viscosity values ranging between 10 16 and 10 20 Pa s, we therefore consider a single value for α set to 0.3 for all calculations shown hereafter. Following Castillo-Rogez et al. ( 2011), we assume that β µ α-1 η -α , which is justified for olivine minerals [START_REF] Tan | High-temperature viscoelasticity of fine-grained polycrystalline olivine[END_REF][START_REF] Jackson | Grain-size-sensitive seismic wave attenuation in polycrystalline olivine[END_REF]. The Andrade parameter β is automatically varied as the melt fraction varies because of its dependency on µ and η. Typical values range between 5.7 × 10 -14 and 7.2 × 10 -12 , in agreement with [START_REF] Jackson | Sound velocities and elasticity of aluminous mgsio3 perovskite: implications for aluminum heterogeneity in earth's lower mantle[END_REF]. In our model, the Andrade rheology is considered for the deep mantle underneath the asthenosphere as well as for the crust.

Matching Io's heat budget for various melt distribution

To test various distributions of melt between Io's deep mantle and asthenosphere, we vary systematically the viscosity of the solid matrix η sol and the mantle and asthenospheric melt fractions, φ mant and φ ast respectively, over a large range of possible values (10 16 -10 20 Pa s for the solid-state viscosity and 0-35% for the whole range of asthenospheric and mantle melt fractions). For a given viscosity of the solid matrix and configuration of the asthenosphere (depth, thickness, φ mant /φ ast ratio), we determine iteratively the asthenospheric melt fraction (see typical melt profile in Figure 1b) that match Io's heat budget, typically ranging between 65 and 125 TW (e.g. [START_REF] Lainey | Strong tidal dissipation in io and jupiter from astrometric observations[END_REF].

Figure 2 illustrates the solutions in terms of φ ast /φ c and η sol (Figure 2a) able to produce the appropriate range for Io's heat budget. Varying η sol and φ mant /φ ast also means exploring different ratios of the tidal power produced in the asthenosphere P ast relative to the total tidal power within Io P glob (Figure 2b). For the example displayed in Figure 2, two ratios φ mant /φ ast are compared: 1/10 (in magenta) and 1/3 (in black), both producing a total power P glob of 100 TW. Whatever the assumed ratio φ mant /φ ast , two distinct regimes can be distinguished. As long as φ ast < φ c , tidal dissipation in the mantle is predominant. Whatever the assumed melt fraction in the manlte (φ mant = φ ast /10 or φ ast /3), the dissipation in the asthenosphere becomes predominant for almost the same value of φ ast , about 1 point above φ c . This transition between mantle-dominated and asthenosphere-dominated regime correspond to different solid-state viscosity values: 1.6 × 10 17 Pa s for φ mant /φ ast = 1/10 and 1.5 × 10 18 Pa s for φ mant /φ ast = 1/3. The same behavior is observed whatever the total tidal power P glob (65, 100, 125 TW)

with only a shift in the location of the transition. In the following, we consider as a reference model φ mant /φ ast =1/3 and P glob =100 TW.

[Figure 2 about here.]

Tidal heating patterns

In order to better understand the difference in terms of heat production between the mantle-dominated and asthenosphere-dominated regimes, the radial and lateral distributions of the tidal heat produced in Io's interior are shown as a function of the solid-state viscosity in Figure 3. Low viscosities (< 1.5×10 18 Pa s) allow for a power produced in the mantle greater than 50 TW (Figure 3a, orange curve). The contribution of the crust (blue line in Figure 3a) is negligible (< 2 TW over the whole range of explored parameters).

For the reference case displayed on Figure 3, the relative contribution of the asthenospheric tidal power becomes larger than the mantle counterpart at viscosities larger than 1.5×10 18 Pa s. Among these calculations, we select five representative models to illustrate the transition from one regime to another in terms of spatial distribution of tidal heating rate H tide and resulting surface pattern of radially integrated tidal heat flux q tide (Figure 3c, bottom and top respectively). For that purpose, the local tidal heating rate per unit of volume H tide (r, θ, φ) averaged over one cycle is evaluated at any point inside the body from the complex stress and strain tensors, determined from the radial functions (see Tobie et al. (2005) and appendices in Kervazo et al.

(2021) for details). By integrating radially the volumetric heating rate H tide over the viscoelastic layers, we then compute the tidal surface heat flux q tide . Such a simplified representation does not take into account the complexity of internal heat transfer (convection and melt extraction) which may affect the surface heat flux pattern, so that comparison with real surface data should be done with care. However, such an approach is similar to what has been classically done in the literature (e.g. [START_REF] Segatz | Tidal dissipation, surface heat flow, and figure of viscoelastic models of io[END_REF][START_REF] Beuthe | Spatial patterns of tidal heating[END_REF][START_REF] Hamilton | Spatial distribution of volcanoes on io: Implications for tidal heating and magma ascent[END_REF]Steinke et al., 2020a), so it makes sense for comparison with previous studies. All the results displayed on Figure 3 [Figure 3 about here.]

The degree-two shape of the tidal potential results in a modulation of tidal heating with latitude and longitude. As already shown by previous studies [START_REF] Segatz | Tidal dissipation, surface heat flow, and figure of viscoelastic models of io[END_REF][START_REF] Beuthe | Spatial patterns of tidal heating[END_REF][START_REF] Hamilton | Spatial distribution of volcanoes on io: Implications for tidal heating and magma ascent[END_REF], these modulations are expressed differently in terms of the spatial pattern of tidal heat flux depending on the contribution of the asthenosphere versus the mantle (Figure 3c top). The obtained heat flux pattern can be interpreted as the linear combination of two end-member patterns:

• a mantle-dominated pattern (pattern I on Figure 3c) corresponding to maxima at the poles (∼ 4.5 W/m 2 , Figure 3b) and minima along the equator at the sub-and anti-Jovian points (∼ 1 W/m 2 , Figure 3b); an averaged value of ∼ 2.2 W/m 2 is observed between -45 and +45 o latitude on Figure 3b;

• an asthenosphere-dominated pattern (pattern V on Figure 3c) corresponding to maxima along the equator at proximity of the sub-and anti-Jovian points (∼ 3 W/m 2 , Figure 3b) and minima at the poles (∼ 1.8 W/m 2 , Figure 3b).

The pattern obtained for P ast P mant = 50%P glob (Pattern III) corresponds to a combination of these two end-member patterns resulting in moderate heat flux variations, corresponding to less than 0.6 W/m 2 of variation between the maxima (3 W/m 2 ) and the minima (2.4 W/m 2 ). Moderate changes of asthenosphere/mantle heat repartition (± 5 TW, patterns II and IV) can, however, result in significant changes in heat flow pattern, with a progressive decay of the polar contribution (see the significant change of polar heat flux between model II and IV Figure 3b).

We note that the asthenosphere-dominated pattern obtained here is dif-ferent that the one classically published in previous studies [START_REF] Segatz | Tidal dissipation, surface heat flow, and figure of viscoelastic models of io[END_REF][START_REF] Beuthe | Spatial patterns of tidal heating[END_REF][START_REF] Hamilton | Spatial distribution of volcanoes on io: Implications for tidal heating and magma ascent[END_REF]Steinke et al., 2020a), which exhibits maxima at mid latitude near the sub-and anti-jovian points. This results from the incorporation of bulk dissipation which strongly affects the response of the partially molten asthenosphere when φ ast > φ c [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF]. Our obtained asthenospheric pattern is similar to Pattern A predicted by [START_REF] Beuthe | Spatial patterns of tidal heating[END_REF], which is expected to be the pattern associated to bulk dissipation (see Eq. ( 23) in [START_REF] Beuthe | Spatial patterns of tidal heating[END_REF]. When bulk dissipation is considered, the heat flow pattern in the asthenosphere is characterized by a preponderance of the σ rr rr components, resulting in maxima at the equator where the radial displacement and energy associated with change of volume are maximal [START_REF] Beuthe | Spatial patterns of tidal heating[END_REF][START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF]. Note that the radial-angular decomposition approach proposed by Beuthe ( 2013) is fully equivalent to our approach where we directly compute the patterns from the full stress and strain tensors [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF]. As shown by [START_REF] Beuthe | Spatial patterns of tidal heating[END_REF], in absence of bulk dissipation, asthenospheric dissipation should be charaterized by a strong reduction of heat flow in polar regions and heat flow maxima around 30 • latitude (pattern B in [START_REF] Beuthe | Spatial patterns of tidal heating[END_REF]. The incorporation of bulk dissipation limits the heat flow reduction in polar region and enhances heat flow in the equatorial band. It is difficult to determine which pattern is more consistent with existing data [START_REF] Davies | Map of io's volcanic heat flow[END_REF], in particular due to the lack of observation in the polar regions, but future measurements may be more conclusive and may confirm if bulk dissipation play a role in the heat budget of Io.

In summary, with our modeling approach, we can determine consistently the repartition of dissipation between the asthenosphere and the mantle for a wide range of mantle viscosity and melt fraction. R3 The asthenospheredominated regime occurs only if the melt fraction in the asthenosphere, φ ast , exceeds the critical melt fraction φ c (=30% in the case considered here). In this regime, in order to match the total power of 100 TW, the viscosity in the mantle η sol must be larger than 7 × 10 16 Pa s -4 × 10 17 Pa s, depending on the melt profile. In absence of melt, a mantle viscosity as low as 10 16

Pa s is required to produce a total power of 100 TW. The mantle-dominated regime is characterized by deep heat production in the polar regions with volumetric heating rate as high as 10 -5 W/m 3 and always larger than 10 -6

W/m 3 (Figure 3c, bottom), which may strongly impact the largescale convective dynamics of the mantle (e.g. [START_REF] Tackley | Three-dimensional simulations of mantle convection in io[END_REF]. The asthenospheredominated regime is characterized by even stronger heating rate, maximum at the mantle/asthenosphere interface and further enhanced in the equatorial regions. Such a strong heating rate may also strongly impact the convective heat transfer through the asthenosphere (e.g [START_REF] Tackley | Three-dimensional simulations of mantle convection in io[END_REF], which in return may influence the resulting surface heat. Future modeling will be needed to assess the resulting surface heat flow pattern, but we can already anticipate that the difference of average heat flow between the polar and equatorial regions may be a good diagnostic of the dissipative regime.

Predicted Love numbers

In this section, we explore the influence of melt disbtribution in the deep mantle and the asthenosphere on the Love numbers, k 2 , h 2 and l 2 , which characterize the global response of the interior to tidal forcing in terms of in-duced potential, radial and horizontal displacements, respectively, and which could be assessed by future missions. We show how the Love numbers vary as a function of another independent observable, the polar heat flux, which, as we will show, is a good proxy of the heat repartition between the asthenosphere and the mantle R4,R11 . Before presenting the predicted range of Love numbers as a function of heat repartition in the interior, we first evaluate the role of bulk rheology on the global response. As discussed in Section 4 and shown in [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF], bulk dissipation can have a major effect when melt fractions exceed φ c , which is possibly expected in the asthenosphere.

Contribution of bulk viscoelasticity in the asthenosphere

The general trends of the Love numbers and the ratio P ast /P glob R4,R11

for a global heat budget of 100 TW as a function of the polar heat flux F pole for various asthenospheric melt fractions R4,R11 are displayed in Figure 4. In order to quantify the role of bulk dissipation in this regime, we specifically report results obtained for k 2 , h 2 , l 2 and P ast /P glob with (colored symbols)

and without (dashed black lines) bulk viscosity. Four asthenosphere configu- and two values of the crust-asthenosphere interface depth (30 km (circles)

and 60 km (triangles)). As expected, the contribution of bulk dissipation is negligible in the mantle-dominated regime (i.e. for φ ast /φ c <1, blue to orange color in Figure 4). The two solutions (with or without bulk viscosity) are very close for F pole > 3.2 -3.8 R4,R11 , and diverge below R4,R11 this value (corresponding to red colors in Figure 4), highlighting the rheological transition (φ ast /φ c =1, depicted as black circles).

In the asthenosphere-dominated regime, results obtained without considering bulk viscoelasticity for k 2 , h 2 and l 2 generally present a progressive increase with increasing F pole , slightly more pronounced for b ast = 50 km than for b ast = 100 km in the case of k 2 (Figure 4). Once bulk dissipation is included, the trend with increasing F pole is inverted for the radial and horizontal displacement Love numbers, h 2 and l 2 , with a strong amplification for the case with b ast = 50 km. Bulk viscoelasticity also impacts the k 2

value, but to a lesser extent. We interpret this particular influence on h 2 as a result of the impact of bulk viscoelasticity on the radial component of the strain and stress tensor [START_REF] Kervazo | Solid tides in io's partially molten interior: Contribution of bulk dissipation[END_REF], which influence is even more pronounced for a thin asthenosphere. The depth of the crust-asthenosphere interface has however no significant effect on the results, leading to increased values of Love number less than 10% when a depth of 60 km is considered instead of 30 km.

As discussed in Section 4, the polar heat flux is high for the mantledominated regime (> 3.5 W/m 2 ), compared to the asthenospheric one (comprised between 0.5 and 3.5 W/m 2 ), whatever the rheological model assumption (i.e. with or without bulk dissipation). The polar heat flux appears to be inversely proportionnal to P ast /P glob , and therefore may be used as a good indicator of the heat repartition between the mantle and the asthenosphere R4,R11 . Considering bulk dissipation in addition to shear dissipation leads to a polar heat flux up to 50% higher than for cases taking into account shear dissipation only, but the linear trend remains R4,R11 . Compared to the Love numbers, decreasing the asthenospheric thickness is not associated to a pronounced enhancement of F pole .

[Figure 4 For h 2 and l 2 (Figure 5, middle lines), more contrasting variations can be observed between the two regimes, especially for a thin asthenosphere as already discussed in section 5.1. For asthenosphere thicknesses of 100-200 km, h 2 and l 2 in the asthenosphere-dominated regime are always larger by about 30% and 50%, respectively, than the ones in the mantle-dominated regime.

For a 50-km asthenosphere, a strong amplification is observed between the two regimes, h 2 and l 2 being up to four times and seven times, respectively, larger in the asthenosphere-dominated regime. Moreover, in this particular case of a thin layer, the effect of the global tidal power is modest for the mantle-dominated dissipation regime (30% for h 2 and l 2 ) but considerable in the asthenospheric-dominated regime (70% for h 2 and 65% for l 2 ).

[Figure 5 about here.]

Synthesis

Besides the two main factors discussed before, the global power (between 65 and 125 TW) and the asthenospheric thickness (between 50 and 200 km), several interior parameters, whose values are still uncertain, may affect the Love numbers k 2 , h 2 and l 2 and the polar heat flux, all of them being observable values R4,R11 . Figure 6, which summarizes the expected range of variations for different interior parameters, shows that k 2 and to a lesser extent, h 2 , are most sensitive to the core size and the reference (solid-state) value of the shear modulus. By contrast, l 2 and to a lesser extent, F pole , are mostly sensitive to the dissipation regime.

Even if k 2 is measured with high accuracy, it appears difficult to determine the mantle/asthenosphere structure from the sole measurement of k 2 as there is a strong overlap between the two dissipation regimes (mantledominated (blue lines) and asthenosphere-dominated (red lines), Figure 6a) and as the interpretation of k 2 will depend on the core size and reference shear modulus. For h 2 (Figure 6b), even if there is a clear separation between mantle-dominated and asthenosphere-dominated solutions, controls by the core size and reference shear modulus would also make the interpretation ambiguous.

As shown on Figure 6c, l 2 would be much more discriminating regarding the dissipation regimes. Love number l 2 is much less sensitive to core size and shear modulus, and mostly depends on the melt fraction in the asthenosphere.

The same tendency is obtained for the polar heat flux (Figure 6d), even if there is a less clear distinction between the two regimes, especially as there is an uncertainty on the global heat power. Measuring l 2 with a precision of 0.01 would distinguish whether Io's heat budget is dominated by mantle or asthenosphere dissipation and would then put constraint on the melt fraction.

As noted earlier, the value of the polar heat flux F pole also mostly reflects whether Io lies in the mantle-or asthenosphere-dominated heating regimes and is roughly independent of other parameters. A precise evaluation may provide an independent complementary constraint and enable to discriminate between these. Additionally, F pole is globally more sensitive to Io's actual heat budget and may help refine its location in the range 65-125 TW.

[Figure 6 about here.]

Discussion and conclusion

In the present study, we focused on Io's interior models including varying percents R3 of partial melting in the deep mantle and asthenosphere (a layer substantially less viscous than the deeper regions but still considerably more viscous than a fully liquid magma ocean). We quantified how the tidal response varies as a function of R3 the viscosity of the solid matrix and the asthenospheric melt fraction relative to the rheological critical melt fraction, φ c (which value may vary between 20 and 40% in the literature). For a melt fraction in the asthenosphere below φ c and viscosity below 10 17 -10 18 Pa s, we show that R3 the dissipation is dominated by mantle dissipation, while dissipation in the asthenosphere becomes dominant above this critical value.

Consistent with previous studies (e.g. [START_REF] Segatz | Tidal dissipation, surface heat flow, and figure of viscoelastic models of io[END_REF][START_REF] Hamilton | Spatial distribution of volcanoes on io: Implications for tidal heating and magma ascent[END_REF][START_REF] Beuthe | Spatial patterns of tidal heating[END_REF]Steinke et al., 2020a), our calculations show that R3 the transition between mantle-dominated and asthenosphere-dominated regimes results in significant changes in terms of heat flux patterns, with strong reduction of the polar heat flux with increasing contribution of the asthenospheric dissipation. The main novelty in our study is to quantify the role of bulk dissipation and to determine the consequences for tidal heat patterns for a wide range of melt fractions R3 . For asthenospheric melt fraction above φ c , the bulk viscoelastic response is predicted to affect significantly the dissipation process, resulting in major changes of tidal heating distribution compared to classical approaches neglecting bulk dissipation (e.g. [START_REF] Segatz | Tidal dissipation, surface heat flow, and figure of viscoelastic models of io[END_REF][START_REF] Beuthe | Spatial patterns of tidal heating[END_REF][START_REF] Hamilton | Spatial distribution of volcanoes on io: Implications for tidal heating and magma ascent[END_REF]Steinke et al., 2020a).

Interestingly, for the asthenosphere-dominated regime (Patterns IV and V on Figure 3), the maxima of heat production along the equator are shifted by ±25 -30 • relative to the sub-jovian and anti-jovian points, which appears consistent with the maxima of volcano density reported by [START_REF] Hamilton | Spatial distribution of volcanoes on io: Implications for tidal heating and magma ascent[END_REF]. We note, however, that we obtained four maxima at about 30 the comparison between modeled heat flow pattern and observed volcanic activity pattern should be made with caution, as tidal heat production may be affected by unconstrained thermal and chemical heterogeneities. In spite of these limitations, the apparent correlation might constitute the first clue that bulk dissipation plays a significant role in Io's tidal dissipation processes.

This also provides a solution alternative to the longitudinal offset attributed to dissipation in a fully liquid magma ocean (Tyler et al., 2015).

Our computation showed that the polar heat flux provides a good proxy of heat repartition in the interior (mantle-dominated vs. asthenospheredominated regime). A future determination of this value with an accuracy < 1W/m 2 would be sufficient to distinguish the two regimes, for a given assumed heat budget R4,R11 . Despite the drastic changes in terms of tidal heating distribution and polar heat flux R 4,11 between mantle-dominated and asthenosphere-dominated regimes, the change in terms of gravitational Love number, k 2 , remains very subtle. Consistent with the results of [START_REF] Bierson | A test for io's magma ocean: Modeling tidal dissipation with a partially molten mantle[END_REF], we found that, for a partially molten interior with φ < 35%, the Love number k 2 should range between 0.07 and 0.11, in the case of a lowdensity metallic core. However, even if a low k 2 value (0.07-0.08) would be more consistent with the mantle-dominated regime, there is no clear separation between the two regimes as several other interior parameters (e.g. shear modulus, core density) may affect the amplitude of tides. In the case of a partially molten interior, at least for melt fractions lower than ∼ 35%, the k 2 Love number is most sensitive to the density of the metallic core and to a lesser extent to the shear modulus of the solid matrix. A k 2 Love number below 0.06-0.07 will suggest a high density core while a value significantly above 0.1 will suggest a much higher melt fraction. For a fully liquid magma ocean, we predict that the Love number depends mostly on the crust thickness and should range between 0.6 and 0.9 for crust thicknesses ranging between 90 and 30 km, consistent with the results of [START_REF] Bierson | A test for io's magma ocean: Modeling tidal dissipation with a partially molten mantle[END_REF] and Van Hoolst et al. (2020).

A multiple-flyby mission to Io, similar to R1 the Io Volcano Observer (IVO, McEwen et al., 2019), would have R1 the capability to provide quantifications of k 2 from radio tracking with an expected accuracy of ∼0.05 (Ryan Park, pers. comm. in Van Hoolst et al., 2020). Such a precision will be sufficient to discriminate univocally between a partially molten interior and a fully liquid magma ocean, but it will not be sufficient to provide constraints on the interior structure and melt distribution if k 2 < 0.1. A dedicated orbiter mission may provide much better accuracy, comparable to what is expected on Ganymede (0.0001) with the ESA JUpiter ICy moon Explorer (JUICE, [START_REF] Cappuccio | Ganymede's gravity, tides and rotational state from juice's 3gm experiment simulation[END_REF]. However, based on our results, k 2 alone will be sufficient to constrain the melt distribution. Before any hypothetical selection of a future mission R8 , close flybys of Io during the extended phase of the JUNO mission [START_REF] Bolton | The juno mission[END_REF], may provide first assessment of k 2 or at least put some upper limit on the tidal amplitude.

The displacement Love number, radial h 2 and horizontal l 2 , may be more conclusive, as they show a more clear distinction between mantle-dominated and asthenosphere-dominated regimes. Values of h 2 smaller than 0.15 would indicate a mantle-dominated regime, but similarly to k 2 some ambiguity would remain as h 2 varies significantly with the core density and shear modulus of the solid mantle. Only a value smaller than 0.1 will be a confirmation that dissipation is mostly located in the deep mantle. Love number l 2 is more conclusive as there is a clear separation between mantle-dominated and asthenosphere-dominated regimes and as l 2 is less sensitive to the interior parameters.

The main challenge is to reach a sufficient accuracy to detect the subtle difference in h 2 and l 2 (±0.01) between the two regimes. [START_REF] Park | Advanced pointing imaging camera (apic) for planetary science and mission opportunities[END_REF] indicated that high resolution imaging may allow the determination of the Love numbers h 2 and l 2 with accuracy of 0.09 and 0.07 after ten flybys.

Such an accuracy will be sufficient to determine whether Io has a global magma ocean or only a partially molten mantle, but it will remain insufficient to distinguish between the two dissipation regimes. Nevertheless, the first estimation of [START_REF] Park | Advanced pointing imaging camera (apic) for planetary science and mission opportunities[END_REF] indicates that the required accuracy may be reached with a mission with higher number of flybys or directly in orbit around Io. Accurate measurements of horizontal displacement will also be essential to constrain the libration (Van Hoolst et al., 2020). Combined determination of Love numbers [START_REF] Park | Advanced pointing imaging camera (apic) for planetary science and mission opportunities[END_REF], libration (Van Hoolst et al., 2020) and magnetic induction [START_REF] Khurana | Evidence of a global magma ocean in io's interior[END_REF][START_REF] De Kleer | Tidal heating: Lessons from io and the jovian system-final report[END_REF] will be required to fully reconstruct the melt profile in Io's interior.

Another important result of our study concerns the possible amplification of the Love numbers, h 2 and l 2 , due to the bulk viscoelastic response in a thin asthenosphere when the melt fraction is above the critical value. In this case, h 2 and l 2 can be amplified by a factor of two, while k 2 remains below 0.1. The combined detection of high h 2 and l 2 values and a low k 2 value by future missions will confirm that bulk dissipation plays a crucial role in the heat budget of Io and will provide constraints on the thickness and melt fraction of the partially molten asthenosphere.

As already mentioned in section 2.1.1, significant uncertainties remain on the determination of the global heat budget and the lateral variations of heat flux. The interpretation of the Love numbers in term of interior structure depends on the total heat budget estimate. Our calculations show that, for similar interior structures and melt distributions, the current uncertainty on the global heat budget (65-125 TW) results in a difference of about 0.01 in k 2 and up to 0.05 in h 2 . Without improvement in the total heat budget, the interpretation of the Love numbers will remain ambiguous. Future astrometric measurements (e.g. [START_REF] Lainey | Strong tidal dissipation in io and jupiter from astrometric observations[END_REF] and volcanic thermal emission mapping (e.g. Veeder et al., 1994Veeder et al., , 2012) ) will be essential to reduce the uncertainty on Io's global heat budget. Precise determination of the orbital dynamics of Io can provide accurate estimates of its mean motion rate change and hence of the total energy dissipated in its interior (e.g. [START_REF] Lainey | Strong tidal dissipation in io and jupiter from astrometric observations[END_REF]. The signature of tidal dissipation in Io is observable not only in Io's orbital dynamics but also in Europa's and Ganymede's ones as the three moons are coupled through the Laplace resonance (e.g. [START_REF] Dirkx | Dynamical modelling of the galilean moons for the juice mission[END_REF]. Spacecraft radio tracking during close flybys of the Galilean moons (JUNO at Io, Europa Clipper and JUICE at Europa, Ganymede and Callisto) combined with optical and UV astrometric data will allow the quantification of the orbital dynamics of the Galilean moons with unprecedented accuracy and a precise determination of total dissipated power independently from thermal emission mapping.

Thermal emission and volcanic landforms mapping are also essential to determine the lateral variations of heat production and transport and possible correlation with the different dissipation regimes. Even if previous mapping suggests both latitudinal and longitudinal variations in volcanic activity [START_REF] Hamilton | Spatial distribution of volcanoes on io: Implications for tidal heating and magma ascent[END_REF][START_REF] Davies | Map of io's volcanic heat flow[END_REF], a lack of complete surface coverage, in particular in the polar regions, makes the comparison between observation and model prediction not conclusive. The general consensus is that more polar observations are needed to fully address this question [START_REF] Rathbun | The global distribution of active ionian volcanoes and implications for tidal heating models[END_REF][START_REF] De Kleer | Tidal heating: Lessons from io and the jovian system-final report[END_REF]. Our calculations confirm that determining the polar heat flux is essential to constrain the dissipative regime of Io R4,R11 . to the contribution of bulk dissipation. Our calculations also showed that, once above the rheological critical melt fraction, small changes in melt fraction can have a large impact on the dissipation rate. Strong retroaction between melt-induced tidal heating and melt production may lead to significant lateral variations in both melt content in the asthenosphere and crustal thickness, which may be amplified by bulk viscoelstic response. Convective transport may also strongly affect melt and heat transport [START_REF] Tackley | Three-dimensional simulations of mantle convection in io[END_REF]Steinke et al., 2020b), making the surface expression of melt and heat production rather complicated [START_REF] De Pater | A 2020 observational perspective of io[END_REF]. Future modeling efforts are required to take into account lateral variations on tidal dissipation rate, including both shear and bulk dissipation and their consequences on melt and heat transport to the surface. Long-wavelength topography and gravity data combined with magnetic induction measurements by future missions will be essential to constrain the lateral variations in depth and melt content of the asthenosphere. Comparison with maps of heat flow and volcanic edifices may reveal the complex interplay between tidal heating, melt generation and extraction, which are likely the main drivers of Io's evolution.

  reproduce the same average heat flux over Io's surface (total power equal to 100 TW), but they result in different dissipation patterns and local values of the surface heat flux owing to change in dissipation regimes.

  rations are considered corresponding to two different values of asthenospheric thickness (b ast = 50 km (left column) and b ast = 100 km (right column))

  Figure 5. Whatever the assumed asthenosphere thickness and heat budget (65-125 TW), a comparable change in k 2 amplitude is obtained between mantle-dominated and asthenosphere-dominated regime (Figure 5, top line).

  The volcanic activity is by nature sporadic (e.g. deKleer et al., 2019b), which creates another complication in the estimation of averaged surface heat flux, representative of the total energy dissipated inside Io. First observations from distant flybys obtained by the instrument JIRAM onboard JUNO have already offered a unique opportunity to map the polar regions[START_REF] Mura | Infrared observations of io from juno[END_REF]. Future close flybys by JUNO may provide a full coverage of the polar region, which is crucial to determine the dominant dissipative regime in Io's interior. Moreover, prior to a dedicated mission to Io, the ESA's JUpiter ICy moon Explorer (JUICE) mission[START_REF] Grasset | Jupiter icy moons explorer (juice): An esa mission to orbit ganymede and to characterise the jupiter system[END_REF]) and NASA's Europa Clipper[START_REF] Howell | Nasa's europa clipper-a mission to a potentially habitable ocean world[END_REF], may also allow complement mapping of hotspot activity from distant flybys, bringing additional constraint on the temporal change of volcanic activity and the heat flow mapping.All calculations presented here neglect lateral variations in interior properties: the crust thickness and the melt fraction are assumed uniform over the entire globe. As indicated bySteinke et al. (2020a) and[START_REF] Spencer | Tidal controls on the lithospheric thickness and topography of io from magmatic segregation and volcanism modelling[END_REF], lateral variations in heat production and transport may result in significant variations in melt content and crust thickness. Based on our computation performed for two different crust thicknesses of 30 and 60 km (Figure4), we can anticipate that lateral variations of crust thickness of the order of a few kilometers as predicted bySteinke et al. (2020a) should have a minor effect on the dissipation rate in the asthenosphere, and almost zero effect on the dissipation rate in the mantle. For larger lateral variations in crust thickness (∼ 20 -30 km) as predicted by[START_REF] Spencer | Tidal controls on the lithospheric thickness and topography of io from magmatic segregation and volcanism modelling[END_REF], the dissipation rate may significantly change in the asthenosphere compared to predictions with constant thickness, especially for a thin asthenosphere (∼ 50 km) due

Figure 2 :

 2 Figure2: Configurations of φ ast and η sol matching Io's heat budget (a) and corresponding fraction of heating in the asthenosphere (P ast /P glob ) (b) for φ mant /φ ast =1/3 (in black) and φ mant /φ ast =1/10 (in magenta). (a) Required asthenospheric melt fraction φ ast /φ c as a function of the solid-state viscosity η sol in order to reproduce Io's heat budget (100 TW is indicated by the solid lines, and the range 65 to 125 TW is displayed by dashed lines and dashed dotted lines respectively). (b) Asthenospheric melt fraction φ ast /φ c as a function of the relative power produced in the asthenosphere (P ast /P glob ). The transition at the critical value φ ast /φ c =1 is depicted in red on both panels.

Figure 5 :Figure 6 :

 56 Figure 5: Influence of Io's tidal heat budget P glob on the tidal Love numbers k 2 , h 2 and l 2 and P ast /P glob R4,11 . These are displayed as a function of the polar heat flux F pole R4,11 . The value of φ ast /φ c is indicated with the color scheme for the reference case P glob =100 TW. The two black lines correspond to endmembers for Io's heat budget: 65 TW (dashed) and 125 TW (dashed dotted). Three asthenospheric thicknesses are investigated: b ast = 50 km, b ast = 100 km and b ast = 200 km. On each curve, the black circle denotes the critical melt fraction.

Table 2 :

 2 Parameter values employed for the various rheological parameters (equations 1-3). ×10 -7 -2.3 × 10 -10 6.39 ×10 -7 -8 × 10 -7 0.102

			η	µ	K
	• s		10 16 -10 20 P as	50-70 GPa	200 GPa
	• l		1 Pa s	10 Pa	1 GPa
	δ		25.7	2.10	2.62
	ξ 1.5 γ 5	5	5
	1 -φ *	R *	0.569	0.597	0.712

Acknowledgments

This research received funding from the French "Agence Nationale de Recherche" ANR (OASIS project, ANR-16-CE31-0023-01) and from CNES for the preparation of the ESA JUICE and NASA Europa Clipper missions. M.B. acknowledges Czech Science Foundation grant no. 19-10809S.

 Tobie, G., Mocquet, A., Sotin, C., 2005. Tidal dissipation within large icy satellites: Applications to Europa and Titan. Icarus 177, 534-549. Turcotte, D.L., Schubert, G., 2002. Geodynamics. Cambridge university press.

Tyler, R.H., Henning, W.G., Hamilton, C.W., 2015. Tidal heating in a magma ocean within jupiter's moon io. The Astrophysical Journal Supplement Series 218, 22.

Van Hoolst, T., Baland, R.M., Trinh, A., Yseboodt, M., Nimmo, F., 2020. (2) the asthenosphere, described by a Maxwell rheology including both shear and bulk viscoelastic rheology; (3) the inviscid liquid metallic core. In the case of the asthenosphere and deep mantle, the effect of melt on the viscoelastic parameters is accounted for. (b) Typical melt profile in Io's mantle. The ratio between the melt fraction in the deep mantle (φ mant ) and in the asthenosphere (φ ast ) is 1/10 (in magenta) and 1/3 (in black) and the reference viscosity η sol is set to 10 19 Pa s. (c) Effect of melt fraction φ on the viscoelastic parameters: shear viscosity η and modulus µ (solid lines), bulk viscosity ζ and modulus K (dashed lines), over the range of interest (from 0 to 35% for the melt fraction). The two endmembers corresponding to η sol = 10 16 and 10 20 Pa s are shown for viscosities. The red vertical line denotes the rheological critical melt fraction (φ c ) where the transition between solid-state and liquid-state behavior occurs.