Can a Combination of Convective and Magmatic Heat Transport in the Mantle Explain Io's Volcanic Pattern?

Tidal dissipation makes Jupiter's moon Io the most volcanically active body in the solar system. Most of the heat generated in the interior is lost through volcanic activity. In this study, we aim to answer the questions: Can convection and melt migration in the mantle explain the spatial characteristics of Io's observed volcanic pattern? And, if so, what constraints does this place on the viscosity and thickness of the convective layer? We examine three different spatial characteristics of Io's volcanic activity: (i) The presence of global volcanism, (ii) the presence of large‐scale variations in Io's volcanic activity, and (iii) the number of Io's volcanic systems. Our study relies on the assumptions that melt in the mantle controls Io's global volcanism, that the large‐scale variations of Io's volcanic activity are caused by nonuniform tidal heating, and that the spatial density of volcanoes correlates with the spatial density of convective anomalies in the mantle. The results show that the observed small and large‐scale characteristics of Io's volcanic pattern can be explained by sublithospheric anomalies influenced and caused by convective flow. Solutions that allow for active volcanism and Io's specific large‐scale variations in volcanic activity range from a thick mantle of a high viscosity ( 10 19 Pa s) to a thin asthenosphere of a low viscosity ( 10 12 Pa s). Provided that Io's volcanoes are induced by convective anomalies in the mantle, we find that more than 80% of Io's internal heat is transported by magmatic processes and that Io's upper mantle needs to be thicker than 50 km.

If sublithospheric processes are assumed to control the locations of volcanoes, the volcanic pattern at the surface can be used to constrain the sublithospheric convective system, such as its thickness, viscosity, and heating rate. A large number of numerical and laboratory experiments of convective systems show that convection signatures, such as the spatial density and the geometry of thermal anomalies, depend on the Rayleigh number and aspect ratio of the system (e.g., Galsa & Lenkey, 2007;Parmentier & Sotin, 2000;Vilella et al., 2018;Zhong, 2005;Zhou & Xia, 2010). If buoyant fluid is omnipresent below the crust, lithospheric anomalies, for example, variations in crustal strength, are the controlling factors. In this case, volcanic chains or regions with elevated heat output are likely to form as a result of fault propagation. The propagation rate and azimuth of the emerging faults are controlled by the stress field in the rocky or icy crust due to tidal deformation, self-gravitation, or other intrinsic sources. An example of this feature is Saturn's moon Enceladus (Crawford & Stevenson, 1985;Porco et al., 2006;Yin & Pappalardo, 2015).
For Io, it is unknown whether sublithospheric or lithospheric anomalies control its volcanic pattern. Observations show that Io's volcanic activity is, at first order, globally distributed with a higher concentration of volcanoes in equatorial regions Hamilton et al., 2013;Kirchoff et al., 2011;Veeder et al., 2015). This large-scale pattern can be seen in both Io's heat output and the spatial density of volcanoes. Many studies attempted to link Io's volcanic activity to interior processes (e.g., Bierson & Nimmo, 2016;Ross et al., 1990;Segatz et al., 1988;Shahnas et al., 2013;Tackley et al., 2001). These studies suggest that the variations in the volcanic activity at the surface are induced by heterogeneous tidal heating below Io's crust. Therefore, viscoelastic models or dynamic magma ocean models of Io's interior are generally evaluated based on the correlation between their resulting tidal heat production pattern and Io's observed volcanic activity pattern (e.g., Hamilton et al., 2013;Segatz et al., 1988;Tyler et al., 2015). However, Io's tidal heat production may be affected by thermal and chemical heterogeneities, which are neglected in spherically symmetric interior models. Thermal heterogeneities may even emerge self-consistently if the nonuniform nature of tidal dissipation is conserved in the temperature field of Io's dissipative layer (Steinke et al., 2020a). Combined with the fact that heat loss mechanisms and dissipative processes in Io's interior are not fully understood yet (Tyler et al., 2015) this could be a reason for the low correlation between derived dissipation patterns and surface observations. Statistical analysis of the spatial relation between adjacent surface features of volcanic origin reveals a clustered behavior (Hamilton et al., 2013). This indicates that volcanic features may share the same anomaly on scales of several tens of kilometers. However, the statistical analysis of only active hot spots, which represent the bulk of Io's heat output, shows uniform spreading (Hamilton et al., 2013). This uniform spreading can arise from: (i) obliteration of older randomly spread volcanic features due to the appearance of new volcanoes (Shoji & Hussmann, 2016), (ii) the competition of adjacent volcanic systems for ascending magma (Hamilton et al., 2013), or (iii) evenly spread thermal anomalies in a convective layer below Io's lithosphere . In any case, a uniform or random spreading implies that Io's largest volcano systems typically do not originate from the same anomaly. Global volcanic chains or bands of high heat output, such as present on Earth and Enceladus, do not seem to dominate Io's volcanic pattern (Kirchoff et al., 2011). The uniform spreading, in addition to Io's significant variations of volcanic activity, favors a sublithospheric control of Io's volcanic pattern. Previous studies have investigated Io's large-scale variations (Bierson & Nimmo, 2016;Hamilton et al., 2013;Tackley, 2001), Io's small-scale variations (Shahnas et al., 2013;Tackley et al., 2001), and Io's global heat output (e.g., Moore, 2001Moore, , 2003 separately from each other. However, as all these characteristics are controlled by the same properties of the sublithospheric layer, a joint analysis of conditions that simultaneously allow for all these characteristics could provide powerful constraints on Io's interior. The objective of this study is to test the hypothesis that convection and magmatic heat transport beneath Io's crust can explain the statistical spreading and the spatial density variations of Io's volcanic pattern. The novelty of our study is that three spatial characteristics of Io's volcanism are combined: (i) the presence of global volcanism, which requires the presence of magma, (ii) the large-scale variations in the volcanic density of spherical harmonic degree 4  , and (iii) the small-scale spreading, that is, the total number of Io's volcanoes. The aim of our analysis is to find properties of Io's upper mantle layer that can explain these three characteristics. Our analysis is based on three fundamental assumptions: (i) the STEINKE ET AL. 10.1029/2020JE006521 2 of 22 presence of volcanism and, therefore, the presence of melt constrain the mantle temperature, (ii) Io's large-scale volcanic activity variations are inherited from nonuniform tidal dissipation and the weakening of the heating pattern is caused by lateral convective flow, and (iii) the spatial frequency of thermal anomalies in the mantle controls the global number of hot-spot features at the surface. These three characteristics are strongly sensitive to our main model parameters, the viscosity, the thickness of the upper mantle layer, and the fraction between convective and magmatic heat transport. Depending on the values of these parameters, the mantle layer could be referred to as a homogeneous mantle ranging from the core-mantle boundary to the base of the lithosphere, a partially molten low-viscosity asthenosphere, or a thin mushy magma ocean.
Although we assume that the sublithospheric layer is convective, this does not exclude that a significant amount of heat could be transported by advecting magma to the crust. Even if melt advection is the dominant mechanism to carry heat, thermal anomalies caused by the underlying convective dynamics are assumed to control the location of melt generation and accumulation, and therefore the location of Io's volcanoes. In contrast to previous studies, which correlate dissipation patterns to Io's observed volcanic activity pattern (e.g., Hamilton et al., 2013;Ross et al., 1990;Tackley et al., 2001), our approach is observation-driven and therefore independent of tidal dissipation models. To identify combinations of model parameters that can explain the small-and large-scale characteristics of Io's volcanic pattern, we use the thermal model developed by Steinke et al. (2020a). This model was originally used to predict Io's mean mantle temperature and variations in Io's surface heat flux due to tidal dissipation and lateral heat flow of steady-state convection. Here, we use this model in an inverse manner by transferring it from the previously model-driven approach to an observation-driven approach. Furthermore, we use the scaling obtained by Vilella et al. (2018) to approximate the total amount of volcanic anomalies driven by convective dynamics in an internally heated system.
In Section 2, we discuss the volcanic observations used for our three-level analysis. The thermal model approximating the effect of the lateral heat flow as well as the scaling law approximating the average distance between volcanic features are introduced in Section 3. In Section 4, we present our results, that is, the global, large-scale and small-scale characteristics of Io's convective mantle, as a function of our model parameters. In Section 5, we discuss the influence of model assumptions on our resulting parameter space as well as possible lithospheric effects on Io's volcanic pattern.

Statistical Analysis of Hot-Spot Observations
In this section, the input for our observation-driven approach is presented. We require two physical parameters that we derive from surface observations: (i) The strength of large-scale variations of Io's magmatic surface heat flux, which is used to quantify the amount of lateral heat flow due to convection and (ii) the number of Io's volcanic anomalies, which will be compared to the number of thermal anomalies arising in the convective system.

Selection of Volcanic Data
Io's geologic surface features and infrared emission have been observed by Voyager, Galileo, and New Horizons (e.g., Davies et al., 2015;Schenk et al., 2001;Veeder et al., 2015;Williams et al., 2011b). Recent Jupiter InfraRed Auroral Mapper observations on-board the Juno mission (Mura et al., 2019) and continuous ground-based observations in the long wavelengths (de Kleer & de Pater, 2016;de Kleer et al., 2019) provide complementary information on the characteristics and location of Io's hot spots. Due to the diversity of these observations, a large number of different data sets of Ionian volcanic features and hot spots are available.
The localized behavior of Io's volcanoes makes the derivation of properties representing an average over space and time challenging. First, it is unknown how the activities of recently active and inactive hot spots evolve over geological time scales. In addition, clustered volcanoes can result from a single or multiple volcanic systems and the influence of intrusive volcanism on Io's heat flux pattern is not well constrained . Finally, variations in the observed volcanic density could originate from a nonuniform coverage of observations ( be able to quantify the influence of these uncertainties on our results, we set up and test two sets for Io's large-scale heat flux pattern, derived from different observations and filter options. For the number of Io's volcanic anomalies, we consider a minimum and a maximum value. Our standard observation-based input contains the locations of recently active hot spots discovered by Galileo, Juno, and Earth-based observation given in Davies et al. (2015), Mura et al. (2019), andde Kleer et al. (2019). We refer to this collection of volcano locations as the standard hot-spot set. The combination of these data sets leads to multiple entries. Our objective is therefore to remove duplicates resulting from uncertainties in the derived locations, while keeping volcanoes of volcanic clusters. As a first step, duplicates in the combined data sets of de Kleer et al. (2019) and Mura et al. (2019) are removed. Then, the indicated uncertainty ranges in the localization of volcanic centers of the de Kleer et al. (2019) and Mura et al. (2019) data set are used to remove any potential overlap with the data set of Davies et al. (2015). The filtering process results in a set of 383 volcanic centers, which are, at first order, globally distributed. The full set of volcanic centers is visualized in Figure 1a. The included Juno observation (Mura et al., 2019) reduces the observation bias of the previous observation sets toward lower latitudes (de Kleer & de Pater, 2016). However, it also introduces an additional observation bias toward Io's leading hemisphere. In order to quantify the effect of observation biases on our results, we do not include the Juno observations by Mura et al. (2019) in our second set.
A second set is adopted from Hamilton et al. (2013). It contains recently active volcanoes as well as regions of former activity, that is, patera floor units (Williams et al., 2011b). In total, 173 hot spots and 529 paterae are included (data sets 1 and 2 in the supplementary material of Hamilton et al., 2013). With 96 detected overlapping features between the hot spots and the paterae, this results in 606 locations. This data set is visualized in Appendix A and is used to show the robustness of our large-scale analysis results presented in Section 4.2.

Spherical Harmonic Analysis of the Volcanic Density
In this section, the volcanic data set is processed such that it can be used as an input for our geophysical model (Section 3). We use Io's spatial volcanic density of the previously introduced standard hot-spot set as a measure for Io's magmatic heat output. Note that we do not to take eruption intensities into account. On the one hand, different instruments and physical models are used for the approximation of the eruption temperatures, and subsequent eruption intensities. This reduces the comparability between different data sets. On the other hand, taking only the intensities of one observation set might also not represent Io's mean heat output. However, observations of Io's active volcanoes reveal that the eruption frequency and intensity of Io's volcanoes are not uniform and that only a small number of active hot spots, for example, Loki Patera, contribute to the bulk of Io's recent total heat flow (de Kleer et al., 2019;Veeder et al., 2012). Qualitative differences between Io's volcanic density pattern and Io's volcanic intensity pattern, and consequences for our results, are discussed in Section 5. STEINKE ET AL.   Spatial density of volcanic features, 10 -6 km -2 We conduct a spherical harmonic analysis in order to examine the large-scale density distribution of Io's volcanoes. Our analysis is based on the technique used by Kirchoff et al. (2011). The normalized spherical harmonic coefficients of degree l and order m are expressed as follows (Johnson & Richards, 2003): . obs N is the total number of volcanic features of our standard hot-spot set, lm P is the associated Legendre polynomial, and n  and n  are the co-latitude and the longitude of each volcanic feature. The coefficients can then be used to reconstruct the volcanic density field as a function of the co-latitude θ and the longitude ϕ (Kirchoff et al., 2011): We only consider spherical harmonics up to degree max 4 L  because the tidal dissipation pattern of a radially symmetric interior structure only contains contributions 4  (Beuthe, 2013). We assume that variations of higher degrees do not originate from the nonuniform tidal heating. Therefore, they are not included in our inversely applied lateral heat flow scaling (Section 3.2). The resulting field obs f is plotted in Figure 1b. The corresponding normalized spherical harmonic coefficients up to degree 4 are available in the data accompanying this paper (Steinke et al., 2020b).
The peak of the volcanic density is located 30-60 degrees east of the anti-sub-Jovian point as previously noted by Kirchoff et al. (2011) and Hamilton et al. (2013). Compared to the observation set adopted from Hamilton et al. (2013), the standard hot-spot set shows stronger contributions from nonsymmetric coefficients, such as 10 S . This is due to the included Juno data set, which mainly contains observations from Io's leading hemisphere (Mura et al., 2019). The resulting normalized coefficients are used as an input for our thermal model introduced in Section 3.2. A map of Io's volcanic density distribution up to degree 4 for the set from Hamilton et al. (2013) is presented in Appendix A.

Total Number of Io's Volcanic Anomalies
In addition to the large-scale variations in the volcanic activity pattern, the total number of anomalies causing Io's volcanism is an important parameter. To determine an upper and lower limit for the number of Io's volcanic anomalies, we need to account for the underlying physical processes affecting the link between the available surface observables and the source of volcanism in Io's interior.
The clustering behavior of volcanoes sharing one source of magma could lead to an overestimation of the total number of Io's volcanic anomalies. Therefore, we cannot simply take the number of volcanic features of our standard hot-spot set. Instead, following Davies et al. (2015), we use min 250 N  to approximate the minimum number of Io's volcanic anomalies. Small variations in this approximation do not affect our results significantly (Section 4.3). Investigating Io's magmatic heat transport through the crust, Spencer et al. (2020a) show that magmatic intrusions are a crucial component of Io's crustal heat balance. Consequently, not all of Io's magma is delivered to the surface. This could lead to an underestimation of the number of magmatic systems at the base of the lithosphere. To approximate the maximum number of anomalies in the upper part of the mantle, we use the number of features included in the Hamilton set (Appendix A) and assume that 80% of the magmatic systems do not reach the surface (Spencer et al., 2020a). This results in max 3030 N  . Altogether, we consider a conservative range of 250-3,030 for the number of Io's volcanic anomalies.

Modeling the Characteristics of Io's Convective Layer
Following our framework, we assume that the melt distribution below Io's lithosphere and Io's surface volcanism are closely linked. We investigate three different mechanisms of the mantle dynamics that can influence the melt distribution below Io's lithosphere: 1. Convection and magma transport regulate the mantle temperature, and therefore Io's melt fraction (Section 3.1). 2. Convection controls how the nonuniformly tidally generated heat is transported laterally, and therefore the extent of Io's large-scale variations of the melt distribution (Section 3.2). 3. Convection imposes regularly distributed cold downwellings and hot anomalies. The distance between these anomalies can be characterized by a specific wavelength whose value depends on the convective vigor. The resulting temperature variations lead to small-scale variations in the melt fraction, which are assumed to control the locations of Io's volcanic features (Section 3.3).
These characteristics of the convection pattern are dependent on parameters of Io's upper mantle, that is, the mantle viscosity, the mantle thickness, and the fraction between convective heat transport and magmatic heat transport. These three parameters are unknown for Io and the resulting parameter space shall be constrained in three levels. In each level, a physical model is applied, which relates the parameter space of the upper mantle layer to one of the three spatial characteristics described above. Furthermore, parameter combinations resulting in incompatible conditions of the thermal system are eliminated in the first-level analysis. This is due to the fact that the three parameters of Io's upper mantle given above, are not independent of each other, and in multiple ways connected to other properties of the thermal system. The aim of our three-level analysis is to progressively disregard all parts of the parameter space that are not able to explain Io's observed volcanic activity field. The fundamental relations, equations, and limitations of our three models are presented below.

Global Presence of Magma
For the first level, we assume that Io's convective layer has no lateral variations. Mantle temperatures above the solidus must be widespread, allowing for a widespread presence of melt, melt migration, and volcanism. This provides strong constraints on Io's interior. We assume that Io's interior is in thermal equilibrium (Lainey et al., 2009) meaning that all of the produced heat is transported to the surface. Thus, prod out , where out Q is the mean heat flux lost to space and prod Q is the heat produced in Io's interior and transported to Io's surface. Furthermore, we assume that heat generated by tidal dissipation within Io's interior is transported by melt advection, conduction, and solid-state convection (e.g., Elder, 2015;Steinke et al., 2020a). Consequently, Io's heat output can be divided into a heat flux cc Q characterized by relatively slow mechanisms, which are conduction and solid-state convection, and a heat flux mag Q characterized by relatively fast heat transport mechanisms via buoyant molten rocks: In the following, we refer to cc Q and mag Q as the convective and the magmatic heat flux, respectively. The ratio between the convective heat flux cc Q and the heat flux out Q is defined as in Steinke et al. (2020a): The determination of the heat flux fraction cc f requires the precise knowledge of the latent heat of the rising magma, the permeability of the mantle, and other parameters controlling the efficiency of magmatic heat transport (Bierson & Nimmo, 2016;Moore, 2001;Spencer et al., 2020a). As these parameters are challenging to obtain, we combine these parameters to one input parameter, namely cc f . Therefore, the above mentioned parameters are still implicitly included in our model. This approach has the crucial advantage to reduce the number of unconstrained parameters for our model. The value of cc f is unknown for Io and consequently one of the three parameters that we constrain in our parameter study.
We assume that due to compositional layering (Spencer et al., 2020b) or simply due to the presence of melt only in the upper part of Io's mantle, a large viscosity difference between the upper and the lower mantle emerges. We assume that the upper, low-viscosity layer of Io is convective (e.g., Monnereau & Dubuffet, 2002;Shahnas et al., 2013;Tackley, 2001;Tackley et al., 2001), and that Io's total heat output is produced in this layer. Depending on the presence of melt, this layer could also represent Io's whole mantle. In accordance with previous studies (e.g., Shahnas et al., 2013;Tackley et al., 2001), the influence of any deeper system is neglected. The thickness d and the viscosity η of the layer are unknown and form the remaining parameters. Any radial changes in the melt and temperature profile within the layer, which naturally arise in a system with magmatic and convective heat transport (e.g., Moore, 2001;Spencer et al., 2020a;Vilella & Kaminski, 2017), are neglected. Other properties within this sublithospheric layer, such as the density m  , the thermal expansivity m  , the conductivity m k and the diffusivity κ, are kept constant throughout the convective system and given in Table 1.
The central links between the convective and magmatic heat transport are the melt fraction and the mantle temperature of the system. These properties are controlled by both the efficiency of the magmatic heat transport (Moore, 2001(Moore, , 2003 and the efficiency of the convective heat transport. Note that even if magmatic segregation dominates the thermal state of the mantle (similar to Bierson & Nimmo, 2016;Moore, 2001;Spencer et al., 2020a), the given mantle conditions, such as the mantle viscosity and the temperature above the solidus temperature, still impose an additional amount of convective heat loss. In order to calculate the equilibrium mantle temperature of the system, we start with a global layer whose heat is partially removed by magmatic heat flux mag Q and check whether a certain parameter combination of  and d can produce the remaining convective heat flux cc cc out Q f Q  , while maintaining a mantle temperature between solidus and liquidus. That way we obtain conditions for which Equation 4 is valid. The Rayleigh-Roberts number of the global layer, which defines the ratio between convection driving and prohibiting forces, is given by: where 2 con con con con is the geometry factor accounting for spherical geometry (Deschamps et al., 2012), and con R is the radial distance from the center of Io to the bottom of the convective layer. We assume that melt does not interact with the convection pattern. The maximum temperature reached in the upper thermal boundary layer is a function of the three unknown parameters η, cc f , and d, and given by Steinke et al. (2020a): Here, a E is the activation energy, gas R the gas constant, and ,crit H Ra the critical Rayleigh-Roberts number.
Since we do not specify parameter values determining the efficiency of the magmatic heat transport, we obtain a solution volume in a three-dimensional parameter space instead of a single curve. The specific parameter combination η, cc f , and d at each point then defines these parameters, for example, the melt fraction, the melt segregation velocity, and the latent heat (Moore, 2003). However, constraints on these parameters imply limitations on the parameter space, which are discussed in Section 5.
Φ is the global melt volume fraction, that is, the porosity of the mantle, given by Hirschmann (2000): with the solidus temperature sol T and the liquidus temperature liq T . Minor nonlinear effects due to latent heat consumption of the different components of the mantle are neglected. The viscosity η of the system is related to the melt volume fraction of the convective system (Mei et al., 2002):  T , and liq T at 100-km depth used here are given in Table 1. We neglect any melting temperature reducing effect of sulfur or other volatiles as Io's high eruption temperatures indicate silicate volcanism (McEwen et al., 1998).

Large-Scale Anomalies
The aim of this second-level analysis is to find additional parameter sets of η, cc f , and d that can be disregarded as they would not be able to explain large-scale variations of Io's observed volcanic activity. Therefore, we assume that Io's volcanic density variations are caused by lateral melt fraction variations in the convective layer, which in turn are caused by large-scale variations in the tidal heat production pattern. The model for analyzing these large-scale variations in the volcanic density field is based on the approach by Steinke et al. (2020a). Our heat transport model is able to incorporate lateral differences in the heat production. In contrast to the model-driven approach by Steinke et al. (2020a), the present model is used in a reverse manner such that it is suited for an observation-driven approach and independent of the underlying heating processes. Note that we do not consider radial differences in the heat production within the convective layer, but only consider the projection of the produced heat Io's global heat production and output remains valid. Heat transport by magma facilitates a more radial heat flow direction because the time scale of rising magma is much shorter than that of solid-state convection. Our model presumes that melt, and the heat carried by melt, are immediately removed from the convective system. Melt migration in lateral direction, such as investigated for Earth Sparks & Parmentier, 1991;Turner et al., 2017), is assumed to not influence the large-scale pattern of the magmatic heat output. Even if melt is transported laterally over hundreds of kilometers to the nearest volcanic center, the investigated heat flux pattern, with wavelengths of 2,000 km and more, would be hardly affected.
A critical addition to the model introduced in Section 3.1 is the angle-wise use of Equations 6, 8, and 10. As discussed in Steinke et al. (2020a), we justify this regional approach with the good agreement between numerical results by Laneuville et al. (2013) and the analytical approach used by Vilella and Kaminski (2017; see further details in Section 4.1 of Vilella & Kaminski, 2017). In addition, we note that our results for the thermal anomalies (Section 4.3) indicate that the typical convection-length is significantly smaller than the wavelengths of the large-scale pattern considered here. Thus, the maximum temperature of the system (Equation 8), found in hot upwellings, can change regionally. Regional changes of the surface fraction between downwellings and hot upwellings are expected to remain small. Therefore, we argue that a regional application of the previous introduced scaling is reasonable.
As discussed in Section 2, we use the volcanic density fields to approximate Io's long-term heat output pattern   mag , Q   . We assume that the melt volume fraction Φ is proportional to the magmatic heat output mag Q . This requires that the percentage of intrusive-to-extrusive volcanism is uniform, implying no large differences in crustal and magma properties. From the pattern of Io's observed volcanic activity, we derive the following: where Φ is the mean melt volume fraction described by the integral over a sphere S of Io's radius Io R : In contrast to cc,out Q , the average resulting from the surface integral of   To obtain   cc,prod , Q   the filter is used reversely compared to Steinke et al. (2020a). Consequently, it has an amplifying effect. The tidally produced heat flux pattern, which is transported by convection, is calculated by: The total heat production pattern projected to Io's surface is obtained by: The inverse application of the blurring is only valid under the assumption that the present-day volcanic distribution mirrors the underlying tidal heating. This implies that the volcanic distribution is not a result of a random distribution or affected by additional interior or surface process other than lateral flow in convection cells with short wavelengths (see results in Section 4.3). The statistical significance of low-degree components of the volcanic density field (Kirchoff et al., 2011) and the good agreement between volcanic density and heat flow maps  for degree 2 indicate that Io's volcanic activity variations result from a physical process. We also include the less significant components of degrees 1, 3, and 4. However, compared to degree 2, we find that these components play only a minor role in constraining our parameter space (not shown). Provided that the assumptions above are valid, the inverse application of the blurring of all significant spectral contributions of the pattern is a unique process. However, since inherently lost components of high degrees cannot be included, the inverted solutions are not complete. An extensive discussion on how uncertainties in the pattern representing Io's heat flux affect our final results can be found in Section 5.
The aim of this procedure is to identify parts of the ηcc f -d parameter space for which this calculation   are disregarded. The procedure for the second-level analysis is visualized in Figure 2.

Number of Small-Scale Anomalies
For our third-level analysis, we estimate the total number of Io's volcanoes from the spatial density of thermal anomalies in a purely internally heated convective system. This estimation is based on the critical assumption that the volcanic features at the surface are related to locations of melt focusing or melt generation below Io's lithosphere. Both processes are induced by small compositional or thermal anomalies. Convection in the mantle with thermal up-and downwellings is a source of small-scale anomalies. The average distance between thermal anomalies of a system depends on the Rayleigh number, or the Rayleigh-Roberts number for an internally heated system, and can be analytically derived  or approximated by a scaling law (e.g., Parmentier & Sotin, 2000;Vilella et al., 2018;Zhong, 2005;Zhou & Xia, 2010).
It may be most intuitive to relate Io's volcanoes to evenly spread hot upwellings. However, volumetric heated systems, such as for Io, are cooled by active downwellings. Upwellings in these systems are usually broad and nonbuoyant (e.g., Parmentier & Sotin, 2000). However, small local temperature variations in the broad upwelling may lead to a significant difference in the melt fraction compared to cooler locations within the same upwelling. Thus, melt may be dominantly generated at the hottest points within these nonbuoyant upwellings (see yellow surrounded features in Figures A2a and A2b). The melt rises through vents to the surface (O'Reilly & Davies, 1981) and volcanoes form just above these hot anomalies.
In the other possible case, melt is focused at downwellings. This requires buoyant melt that has been produced in upwelling regions to migrate laterally over large distances. Magma transport of several hundreds of kilometers has been reported and studied for plume-ridge systems on Earth (Braun & Sohn, 2003). In this study, the authors find that a solid pressure gradient of an uprising plume can facilitate a lateral melt migration over distances of more than 100 km. With typical hot-spot distances of 100-400 km for Io, the lateral heat transport by magma does not affect the low spherical harmonic degrees of the heat flux pattern considered in the large-scale analysis (Section 3.2). Another potential link between surface volcanism and downwellings is delamination. Usually, downwellings lead to a thickening of the lithosphere. However, delamination of the lithosphere has been observed in numerical models of Io's asthenosphere (Schools & Montesi, 2019). Extreme basal heating reduces the lithosphere viscosity. As a result, the base of the lithosphere can more easily be dragged downwards, causing convective thinning of the lithosphere (Spohn & Schubert, 1983) and re-melting of the downwards dragged lithospheric material. The thinned lithosphere and the additional supply of lithospheric melt could facilitate Io's volcanism.   Based on high-resolution numerical simulations, Vilella et al. (2018) derived a scaling to characterize the convection pattern of an internally heated layer. An automatic detection procedure of downwellings at middepth of the convective layer results in the following relationship: where 2 d N is the dimensionless number of cold downwellings per unit area and N C is a scaling factor. The exponent N  is equal to 1 / 4 (Section 5.3 in Vilella et al., 2018). We convert 2 d N into the total number of downwellings: is the distance from the center of Io to the center of the convective layer. This scaling was initially developed for downwellings. Since for Io's heated mantle, the classical view of convective cells, in which every downwelling is associated with an upwelling, is not appropriate, it cannot be transferred for Io's broad upwelling regions. However, Figures 2a and 2b show that the hottest locations within upwelling regions, correlate with the locations of downwellings. A numerical investigation, presented in Appendix B, further shows that Equation 19 can be used to approximate the number of hot anomalies per unit area. As a consequence, Equation 19 is appropriate to approximate the number of volcanic anomalies, independently of whether these anomalies correlate with cold downwellings or hot anomalies within broad upwellings.
We use N C for a rigid boundary condition (Table 1) and assume an uncertainty range of 25  % in this coefficient, which we discuss in Appendix B. Since the exponential nature of the viscosity and the term 4 d in the Rayleigh-Roberts number dominate the resulting total number of volcanoes total N , the influence of the uncertainty in N C on the resulting parameter space is small (Section 4.3).

Results
Here, we show how our initial parameter space is reduced by disregarding combinations of parameters that cannot explain the observed spatial characteristics of Io's volcanic activity field. This is done in three levels. Each level represents a different spatial characteristic as presented in Section 3. While for the first level, that is, the calculation of Io's global mantle temperature, only Io-specific values given in Table 1 are required, the analysis of large-scale anomalies (second level) and small-scale anomalies (third level) requires Io's observed volcanic density distribution and the range of Io's volcanic centers presented in Section 2.

Global Presence of Magma
We use Equations 6 and 8 to calculate the mantle temperature and Rayleigh-Roberts number for our parameter space consisting of the viscosity η, the heat flux fraction cc f , and the thickness d of the convective layer (Table 1).
An overview of the three-dimensional parameter space with all remaining solutions is presented in Figure 3a. The domain forms a band from low viscosities and high heat flux fractions to high viscosities and low heat flux fractions (Figures 3b and 3c). Parameter combinations for which the average melt volume fraction is not compatible with the viscosity η are disregarded (Equation 11). For this melt-viscosity constraint, we use a conservative approach with a wide range for the initial melt-free viscosity 0  (Table 1). Additionally, models resulting in mantle temperatures above the liquidus temperature or below the solidus temperature are disregarded. This is because we only consider solutions with a heat flux value Constraints on the ηcc f -plane determining the thickness of the solution band (Figures 3b and 3c) result from these two restrictions. Note that at large d, small η, and large cc f the obtained Rayleigh-Roberts numbers are larger than Apart from our three main model parameters η, cc f , and d, uncertain properties of the convective system influence the mantle temperature and therefore the solution space. An important model uncertainty is the boundary condition between the convective layer and Io's lithosphere. In case a free-slip boundary instead of a rigid boundary is used, the critical Rayleigh-Roberts number decreases from H,crit 2772 Ra  to H,crit 868 Ra  (Vilella & Kaminski, 2017). In that case, the solution band in the ηcc f -plane is shifted by approximately half an order of magnitude toward higher mantle viscosities. The same effect occurs if 0  is reduced by 13% or if out Q is reduced by 30%. In addition, the liquidus and solidus temperatures of Io's mantle are uncertain. Changing these temperatures leads to a shift, thinning, or thickening of the valid parameter space in the ηcc f -plane. In particular, the parameter combinations with high cc f values are affected. For regimes in which heat is mainly extracted by melt ( cc f close to 0), the temperature is not as strongly affected by the viscosity as for a convection-dominated regime (Figures 3b and 3c). However, compared to the full parameter range investigated here, the shift would remain small. This indicates that our results are not strongly affected by uncertainties in these manlte properties.

Large-Scale Anomalies
For the large-scale analysis, we use the observation-driven approach presented in Section 3.2 and our standard hot-spot set presented in Section 2. Figure 4a shows the globally average mantle temperature ˆm T for the parameter space that corresponds to Io's large-scale volcanic density variations.
The cut at low viscosities ( This depends mostly on the Rayleigh-Roberts number, which varies between 12 10 and 13 10 at the cutoff, but also directly on the convective heat flux fraction cc f and the layer thickness d. We conclude that thick and highly turbulent convective systems cannot explain Io's volcanic distribution. STEINKE ET AL.  Additional thinning in the ηcc f -plane (Figures 4b and 4c) arises due to the laterally varying temperature profile, which locally exceeds or falls below the solidus-liquidus temperature range. In particular, parameter combinations with large peak-to-peak temperature variations and high average mantle temperatures are affected by this constraint. Furthermore, parameter combinations with a small layer thickness d, small cc f , and low viscosity η result in a Rayleigh-Robert number below the critical Rayleigh-Roberts number (Table 1). As convection would cease, this part of the parameter space is also disregarded.
The parameter space in agreement with the observations set adopted from Hamilton et al. (2013) is plotted in Appendix A. Both observation sets result in similar solutions. This indicates that the influence of data set choices, such as the filtering of volcano locations and the inclusion of Juno observations (Mura et al., 2019), is small.

Number of Small-Scale Anomalies
Finally, we focus on the number of convection-driven thermal anomalies, which we assume to correlate with Io's volcanic features (Section 3.3). In Figures 5a-5d, we present the numbers of thermal anomalies for the three-dimensional parameter space remaining from Section 4.2. We take the range of 250-3,030 volcanic centers derived for Io in Section 2 and add a 25% range to account for the uncertainty on the scaling discussed in Appendix B. We disregard all parameter combinations with smaller or larger values. Note that for models with approximately more than 5 10 thermal anomalies, a correlation between thermal anomalies at the base of the lithosphere and Io's volcanic surface features is not provided. Io's volcanoes could be located so close to each other that their magma chambers would form a continuous magma ocean. In this case, the volcano distribution is most likely driven by lithospheric anomalies. As this would violate our fundamental assumption that the volcanoes are correlated with thermal anomalies in the upper mantle, the third-level analysis cannot be applied to this part of the parameter space.
It can be seen that Io's range of volcanic centers can only be explained by a mantle thickness 50 d  km, high mantle viscosities, and a magma-dominated heat transport ( cc f 0.20). Solutions with a thickness below 50 km, either result in a Raleigh-Robert number below the critical Rayleigh-Roberts number, or the scaling results in more than 3,787 anomalies. A limit caused by too few volcanoes arises as well for models STEINKE ET AL. with extremely low heat flux fractions and thick convective layers. A large part of the remaining parameter space results in mantle temperatures just above the solidus temperature (compare Figures 4 and 5). The high uncertainties on Io's crustal thickness hardly affect our results. A crustal thickness of 60 km instead of 30 km changes the number of thermal anomalies by less than 5% (Equation 20).

Discussion
Our results place strong constraints on Io's upper mantle. However, several assumptions and simplifications have been made that require further justification.
The constraints provided by the second-level analysis (large-scale analysis) on the parameter space are controlled by the ratio between the largest amplitudes of Io's magmatic heat flux variations   mag , Q   and Io's mean heat output out Q . However, this ratio is affected by uncertainties. Observations and models indicate a significant amount of background heat flow, that is, heat originating from unaccounted sources (Spencer et al., 2020a;Veeder et al., 2015). The amount and distribution of this background heat flow could alter the ratio between the largest amplitudes of   mag , Q   and out Q and therefore the constraints on our three-dimensional parameter space. Accounting for an additional homogeneous background heat flux would allow for models with higher mantle temperatures and slightly higher Rayleigh-Roberts numbers. An additional uncertainty arises by how the available volcanic observations are used to derive   mag , Q   . We used Io's volcanic density pattern presented in Section 2.2 to characterize the variations in the magmatic heat output. However, Io's volcanoes show a large variation in their eruption intensities and frequencies, which we did not take into account. In order to estimate the effect of this simplification on the resulting parameter space, we perform an additional spherical harmonic analysis to obtain an alternative magmatic heat flux pattern, using the volcanic data set presented in Davies et al. (2015). For that, we follow the weighting approach of Kirchoff et al. (2011) and weight the hot-spot locations with their emitted power given in Davies et al. (2015). We obtain a strong degree 2 component compared to the average detected thermal emission (in agreement with Figure 3 in Davies et al., 2015). The ratio between degree 0 and the degree 2 contribution for the weighted pattern is more than two times larger than for the volcanic density patterns presented in Section 2 and Appendix A. Also other contributions (degrees 1, 3, and 4) show larger amplitudes. However, if we account for an additional homogeneous background heat flux, assuming that 60% of Io's magmatic heat output is not included in the detected emission (see discussion above), the results for the weighted magmatic heat flux pattern do not differ significantly from the results shown in Figures 4a and 1c. This shows that, although the magmatic heat flux pattern itself is affected, the incorporation of the hot-spot intensities has only a minor effect on the investigated parameter space. Only if we further minimize the amount of the background heat flow (degree 0) for the intensity-weighted pattern, we obtain stronger constrains: While the overall shape of the parameter space of Figure 4a remains, minor cutoffs along the concave boundary occur and additional parameter combinations resulting in high average mantle temperatures are disregarded. Therefore, the evenly spread volcanic density pattern can be considered as a conservative choice, while the effects of the above described uncertainties on our conclusions should be small.
Our results suggest that a significant part of Io's produced heat is transported by advecting magma. However, since Io's strongly heated mantle can maintain a sufficiently high Rayleigh-Roberts number to allow for convection, this result is not in contradiction with our hypothesis that the volcanoes are linked to the convection pattern. A magma-dominated heat transport for Io is in agreement with results by Bierson and Nimmo (2016), Moore (2001), andSpencer et al. (2020a). However, additional limitations on our parameter space need to be placed: At melt volume fractions above 30%-40%, fluid melt pockets enclose the solid particles (Moore, 2001). In this case, Io's upper mantle needs to be considered as a purely convective magma ocean with a very low viscosity. This is not in agreement with a partially magmatic and partially convective heat transport. In Section 4.2, all parameter combinations with melt volume fractions above 50% are disregarded due to other constraints. Eliminating solutions with melt volume fractions higher than 40% would further thin the solution space in the ηcc f -plane. However, since the critical threshold is uncertain and only has a small effect on the final parameter space, and therefore not on the conclusions, we do not apply this additional limitation to Figure 4. Furthermore, we note that the given melt fraction and cc f -fraction for each location of our solution space implies a specific melt segregation velocity (Moore et al., 2007). Because these velocities for Io are uncertain (Moore, 2001), we cannot eliminate any further parameter combina-tions. However, in case future work constrains the efficiency of melt to rise, the solution space could reduce further. Although our coupled model design does not violate fundamental thermal concepts within the given boundaries discussed above, the used scalings contain simplifications. Inaccuracies could arise as a result of strong interactions between melt and convective flow. Heat transport in a partially molten and convecting system has not received accurate quantitative treatment in the literature (see further discussion in Breuer & Moore, 2015). For example, magmatic emplacement and melt focusing were not initially included in our used scaling laws (Tackley, 2001;Vilella & Kaminski, 2017;Vilella et al., 2018). These interactions between solid and liquid phases in Io's tidally deformed interior could be driving mechanisms of Io's heat production and heat transport. Extensive lateral heat transport and interactions between tidal heating and rheology in mantle upwellings, shallow magma chambers, or decompacting boundary layers may explain Io's localized heat output anomalies such as Loki Patera. Assuming that these interactions only affect regional scales, we expect the effects on our solution space to be small. Nevertheless, it may be interesting to investigate the detailed effects of two-phase flow on Io's heat transport and heat production in relation to Io's focused heat output. This requires additional theoretical developments and is therefore beyond the scope of this work.
We find that the observed number of Io's volcanoes is orders of magnitudes smaller than the number of thermal anomalies resulting from a shallow and turbulent convective magma ocean with a low viscosity ( Figure 5). Based on our hypothesis that Io's volcanic distribution is controlled by convective motion, our results do not favor the presence of a turbulent convective magma ocean. Instead, our results favor a partially molten layer, with the majority of parameter combinations resulting in low to medium melt volume fractions. In these regimes, active downwellings and delamination (e.g., Schools & Montesi, 2019;Vilella & Kaminski, 2017) could pull lower parts of Io's lithosphere into Io's mantle. This could lead to a continuous recycling of the lithosphere material, required to explain Io's observed eruption temperatures (Keszthelyi et al., 2007;McEwen et al., 1998). Other reasons for Io's observed volcanic density, which are not based on the assumption that today's volcanic activity is controlled by Io's current convection pattern and do not preclude a magma ocean (Keszthelyi et al., 1999), need to be considered. For example, locations of the volcanoes could be a result of ancient impacts. Volcanoes could have formed in an earlier evolutionary stage and remained, while the interior kept evolving. However, the continuous and widespread resurfacing and observed changes in the volcanic activity in location and time (Geissler et al., 2004) makes the theory of frozen-in volcanoes unlikely. Another cause for Io's volcanic distribution could be crustal stresses that induce crustal weaknesses. Long faults accompanied by volcanic chains of global scale, which could indicate fault propagation due to tidal stresses, have not been observed. However, weaknesses in the crust due to thermal stresses, arising as a consequence of Io's continuous resurfacing, have been suggested (Kirchoff & McKinnon, 2009;Kirchoff et al., 2020). These stresses cause the elevation of Io's high mountains and may also cause volcanism. This mechanism would lead to continuous change of volcano locations, in agreement with Geissler et al. (2004). In the case magma forms evenly in the regions of diffuse upwellings, or the distance between adjacent thermal anomalies is sufficiently small for the formation of a continuous magma reservoir, competition for magma could also explain the uniform spreading of Io's hot spots (Hamilton et al., 2013). In order to make definitive conclusions on Io's interior, detailed modeling aiming to answer how magma can rise through Io's thick crust is essential. Stronger constraints on crustal properties and conditions at the base of the crust that invoke heat-piping could potentially rule out some of the possible theories for Io's volcanic activity pattern.

Conclusion
Io is covered by a large number of active hot spots and volcanic features (Figure 1). We investigated whether Io's volcanoes can be explained by convection processes below the lithosphere. Furthermore, we tested what constraints on the viscosity, the thickness, and the fraction between convective and total heat transport of Io's uppermost convective layer can be placed. This study is motivated by the observed uniform distribution of active hot spots (Hamilton et al., 2013), which fits the characteristics of uniformly distributed thermal anomalies in a convective system. Furthermore, global volcanic chains, which would indicate a strong relationship between volcanic activity and tectonics, have not been observed.
One key difference compared to previous studies is the observation-driven approach. In contrast to tidal dissipation modeling studies (e.g., Hamilton et al., 2013;Ross et al., 1990;Segatz et al., 1988;Tyler et al., 2015), our analysis is independent of the underlying tidal dissipation mechanism. The basis for our analysis is a set composed of the location of 383 recently active hot spots detected by satellite-based and Earth-based observations de Kleer et al., 2019;Mura et al., 2019). Another key difference of this study is that we investigated three different spatial characteristics simultaneously. These characteristics of Io's surface are assumed to be linked to corresponding spatial characteristics of convective flow in Io's interior. The three investigated characteristics act at different spatial wavelengths and are as follows: (i) the temperature of the upper mantle allowing for widespread volcanism, (ii) the large-scale (spherical harmonic degree 4  ) temperature variations originating from nonuniform tidal heating, and (iii) small-scale features arising from thermal anomalies in the upper mantle. Using steady-state scaling laws based on previous work by Steinke et al. (2020a) and Vilella et al. (2018), we quantified these three characteristics of convection processes as a function of the above given parameters of the convective layer. This allowed us to step-wise confine the parameter space in agreement with surface observations. Importantly, the three analysis levels are independent from each other, consequently, uncertainties and assumptions on one particular level do not affect the constraints on the parameter space given by other levels.
Our results show that, due to constraints from the mantle temperature and the relationship between temperature, melt, and viscosity, strong constraints can be placed on the three-dimensional parameter space. A thick and turbulent convective system with a low viscosity would erase any large-scale lateral temperature and melt variations in the interior. This is not in agreement with the low-degree variations in Io's volcanic activity pattern. Consequently, parameter combinations resulting in a high Rayleigh-Roberts number are disregarded (Figure 3a). If the total number of volcanic anomalies is related to the spatial density of downwellings or hot anomalies in Io's upper convective layer, additional constraints arise that rule out a fluid-phase dominated magma ocean. Based on observations and accounting for intrusive volcanism and clustered volcanoes, we approximated Io's total number of volcanoes to be between 250 and 3,030. This range of volcanic features can be obtained by a convective layer with a magma-dominated heat transport ( 80%), a mantle viscosity between 15 10 Pa s and 19 10 Pa s, and a upper mantle thickness 50 km.
We conclude that the global, large, and small-scale characteristic of Io's volcanic pattern can be explained by sublithospheric convection-driven anomalies in a system with convective and magmatic heat transport. We showed that our results are robust against parameter uncertainties and variations in the observation input (Section 4.1 and Appendix A). While this suggests that our working hypothesis is consistent with observations, it is important to keep in mind that other processes of lithospheric origin may also explain Io's distinct volcanic distribution (Section 5). More advanced numerical models studying the interactions between Io's asthenosphere, lithosphere, and crust, incorporating solid and liquid flow, are needed. Furthermore, observing Io's tidal deformation, rotation, topography, and gravity field with an Io-dedicated mission (Bierson & Nimmo, 2016;McEwen et al., 2019;Van Hoolst et al., 2020) could allow us to conclude with certainty whether the volcanic distribution is driven by lithospheric or sublithospheric anomalies.

Appendix A: Alternative Observation Set
To show that the obtained constraints on Io's interior are not purely a consequence of our chosen standard observation set, we here introduce a second observation set, which is adopted from Hamilton et al. (2013). It contains 173 active volcanoes as well as 529 volcanic features of former volcanic activity (data set 1 and 2 in the supplementary material of Hamilton et al., 2013). The 702 features are further reduced to 606, due to multiple entries. The locations of these features are presented in Figure A1a. All normalized spherical harmonic coefficients up to degree 4, resulting from the analysis presented in Section 2.2 are available under Steinke et al. (2020b). The adopted set also shows a peak in the volcanic density located east of the anti-sub-Jovian point. Compared to our standard hot-spot set, the set from Hamilton shows a strong- to the detection criteria, causing a large uncertainty on the parameterization of hot anomalies. To account for these uncertainties, we determine 2 d N for various values of the detection criteria and report these variations in Figure A2c. Variations of 2 d N are modest for low H Ra but increase with H Ra as the complexity of the convective flow increases (see Vilella et al., 2018, for more details). Interestingly, we find that the number of the hottest points within broad upwellings and the number of cold downwellings are somewhat equivalent ( Figure A2c). This could be seen in connection to the observation that the location of the hot anomalies are spatially correlated to the location of cold downwellings (Figures A2a and A2b). The only notable difference is the underestimation of 2 d N for large H Ra . However, in Section 4.3, the constraint on the number of volcanoes is used to exclude models with too many volcanoes. As such, this underestimation is not an issue as it implies more conservative results. We can therefore safely use Equation 19 to characterize the number of hottest points within upwelling regions per unit area.

Data Availability Statement
Used data can be found in the supplementary materials of Hamilton et al. (2013) and Davies et al. (2015), Table 5  (2:2). Blue and red colors correspond to material that is colder and hotter than the local horizontal average, respectively. The color scale is not the same in the two panels and was chosen to increase the visibility of the convective structures. Yellow contours illustrate the hot thermal anomalies that have been determined using the detection procedure described in Vilella et al. (2018). The change of the number of thermal anomalies per unit area ( 2 d N ) with the Rayleigh-Roberts number is reported in (c). The blue-shaded area corresponds to the number of hot thermal anomalies. The uncertainties due to the counting method are shown by the vertical extent of this area. The solid red curve is the scaling law for the number of cold anomalies proposed by Vilella et al. (2018), while the gray shaded area illustrates a 25% error on this count.