Magmatic intrusions control Io's crustal thickness

Io, the most volcanically active body in the solar system, loses heat through eruptions of hot lava. Heat is supplied by tidal heating and is thought to be transferred through the mantle by magmatic segregation, a mode of transport that sets it apart from convecting terrestrial planets. We present a model that couples magmatic transport of tidal heat to the volcanic system in the crust, in order to determine the controls on crustal thickness, magmatic intrusions, and eruption rates. We demonstrate that magmatic intrusions are a key component of Io's crustal heat balance; around 80% of the magma delivered to the base of the crust must be emplaced and frozen as plutons to match rough estimates of crustal thickness. As magma ascends from a partially molten mantle into the crust, a decompacting boundary layer forms, which can explain inferred observations of a high-melt-fraction region.


Introduction
Jupiter's moon Io is tidally heated by eccentricity forcing from its mean motion resonance with Europa and Ganymede (Lainey et al., 2009), resulting in extensive surface volcanism. This volcanism has lead to significant interest in understanding Io's internal structure and energy balance. The rate of tidal dissipation coincides closely with the rate of surface heat loss (Davies et al., 2015), implying that Io is close to a state of thermal equilibrium. Further, the spacing of volcanoes is approximately uniform across the surface (Kirchoff et al., 2011). These observations imply that Io's leading-order structure is spherically symmetric and roughly steady state. An understanding of this leading-order structure must serve as the foundation for investigations into spatial heterogeneity and temporal evolution.
Io's radial structure is determined by the heat and mass transport mechanisms operating in its interior.
Energy emission from the surface is concentrated at volcanic features, suggesting that volcanism, not conduction, is the primary heat transport mechanism in the crust. O'Reilly & Davies (1981) proposed that the export of heat across the crust by volcanic systems -a process known as heat piping -allows the growth of a thick crust, which limits the efficiency of conductive heat loss.
Heat transport in Io's mantle is more widely debated. The apparent thermal equilibrium requires that the rate of energy export matches the rate of tidal dissipation. Moore (2003) demonstrated that an equilibrium between convective heat transport and tidal dissipation would occur at melt fractions above disaggregation. However, the expected tidal heat production under these conditions is significantly less than the observed surface heat flux. This suggests that convection cannot be the primary mechanism for delivering heat to the crust. Alternatively, magmatic segregation is capable of transporting the observed tidal heat input at low melt fractions (Breuer & Moore, 2015). It is therefore plausible that Io's tidal heat is removed from the mantle by magmatic segregation and transported across the crust by a volcanic system. We adopt this basic hypothesis in the current study and address specifically what controls the total amount of magma produced, the amounts emplaced intrusively as plutons and extrusively as surface volcanism, and the controls these place on crustal thickness.
The thickness of Io's crust is determined by the depth to which it downwells before it is heated to its melting point. Previous work has not considered the dynamics of magma at the crust-mantle boundary or in the lower crust, but the emplacement of plutons introduces heat to the lower crust, which reduces the crust's thickness and modifies its thermal profile. We present a coupled model of crust and mantle dynamics that assesses these processes. The model is formulated to make predictions of elastic thickness, surface heat fluxes, and globally averaged eruption rates, predictions that can be readily tested by future missions to the Jupiter system. Our results indicate that the heat balance and melt transport mechanisms in the crust and at the crust-mantle boundary ultimately determine the thickness of Io's crust and the melt distribution below it.
Magnetic induction measurements -interpreted in terms of mantle electrical conductivity -have been used to infer the presence of a ∼ 50 km thick layer with ∼ 20% melt fraction (a "magma ocean") beneath Io's crust (Khurana et al., 2011). We note, however, that the interpretation of the induction measurements as a high-melt-fraction region is debated; Blcker et al. (2018) argue that interaction with Io's plasma environment is a better explanation of induction measurements than is a magma ocean.
Nonetheless, previous studies have proposed that this inferred high-melt-fraction layer could be a region of enhanced tidal dissipation (Hamilton et al., 2013;Bierson & Nimmo, 2016). Tidal dissipation theory predicts that for a homogeneous body, dissipation is highest at the center. A low viscosity layer underlying a rigid crust may allow the concentration of dissipation, but models that invoke it must explain how such a structure arises. We use our coupled dynamical model of the crust and mantle to investigate the feasibility of a high-melt-fraction layer occurring without the need for enhanced dissipation.
The manuscript is organised as follows. First we outline the physics of the model before presenting results showing the key controls on i ) crustal thickness, ii ) emplacement rates, and iii ) a high-meltfraction layer beneath the crust. We then discuss the implications of these results for interior structure and evolution.

Model Description
The model, shown schematically in figure 1, considers the dynamical effects of tidal dissipation on Io's crust and mantle. These are modelled as a continuum that is either solid (the crust) or partially molten (the mantle). We model melting, magmatic segregation, and compaction (the contraction of a solid matrix as melt is expelled) with a system of conservation equations for mass, momentum, and energy appropriate for a compacting two-phase medium (McKenzie, 1984). Magmatic flow can also occur in a volcanic plumbing system that stretches from the upper mantle to the surface, and that exchanges mass with the mantle and crust. In the upper mantle, melt can leave the pore space and enter the plumbing system and then, as melt rises in the plumbing system, it can form intrusions (freeze) in the crust, delivering mass and energy to the surroundings. The volcanic flux that reaches the surface (the eruptive flux) instantly cools and imparts a downward flux of cold surface material. The crust-mantle boundary is defined as the depth at which temperature equals the solidus of silicate rock; its location is determined as part of the model. Crustal thickness is thus defined as the distance over which the cold, downwelling surface material heats before it begins to re-melt.
The model invokes some simplifying assumptions. To leading order, Io's volcanoes are uniformly distributed across the surface (Kirchoff et al., 2011;Williams et al., 2011). Therefore we propose a spherically symmetric model, assuming that deviations from spherical symmetry are secondary effects imprinted on the leading-order radial structure. Tidal dissipation models that match Io's surface heat flux utilise either very low viscosities (Steinke et al., 2020) or empirically parameterised rheologies (Bierson & Nimmo, 2016;Renaud & Henning, 2018). Hence, to explore the leading-order dynamics without a dependence on poorly constrained parameters, we take tidal dissipation to be uniformly distributed. Below we assess the melt configurations this produces and discuss whether this may lead to significant radial partitioning of tidal heating. We assume one-component thermodynamics because prolonged and intense melting will have removed the more fusible minerals from the mantle to leave an olivine-dominated lithology (Keszthelyi & McEwen, 1997). We neglect the pressure-dependence of the melting temperature due to the small size of Io and hence the low pressures in the mantle. Finally we assume that melt is mobile in the partially molten mantle due to the large grain size, and hence large permeability, expected for an olivine-dominated, annealed mantle (Lichtenberg et al., 2019).
In the mantle (radii r b < r < r c ), temperature is at the melting point and heat transport occurs solely by magmatic transport of latent heat. Buoyancy causes the upward flow of magma, which is balanced by the downward flow of solid. In the crust (defined as having a spherically-averaged temperature below the melting point, r c < r < R), the temperature drop to the surface drives a conductive heat flux, while the downwelling solid crust transports the cold surface temperature inward (Schenk & Bulmer, 1998). The volcanic plumbing system continues to transport magma and latent heat upward through the crust.
We make minimal assumptions about what this plumbing system actually looks like, but require that it Figure 1: Schematic of the model for Io. Magma rises buoyantly through the mantle while the solid moves downward. At the top of the mantle, magma enters the volcanic plumbing system. Some of this magma is emplaced (intruded) into the ductile lower crust; the rest rises to the surface and fuels volcanic eruptions. The core is excluded from the model. interacts with the solid crust by emplacement of material (the formation of plutonic intrusions). This emplaced material is a crustal source of both mass and heat. The actual volume of the plumbing system is assumed to be negligible (consistent with the flow there being much faster than that of the solid and melt elsewhere).

Model equations
In the partially molten mantle, which has melt fraction φ(r, t), we represent the solid velocity by u and the relative liquid velocity (the Darcy segregation flux) by q. Darcy's law relates the segregation flux to pressure gradients and buoyancy where K 0 φ n is the permeability, n is the permeability exponent, ∆ρ is the density difference between solid and liquid, g = −gr is the gravity vector, η l is the liquid viscosity, and P = (1 − φ)(P liquid − P solid ) is the compaction pressure (Keller et al., 2013).
Transfer of material between the crust-mantle system and the volcanic plumbing system is considered in the conservation of mass equation. Making a Boussinesq approximation, conservation of solid and liquid mass require where Γ is the volume transfer rate of solid into liquid (the melting rate), E is the extraction rate to the plumbing system (a sink of liquid from the mantle), and M is the emplacement rate from the plumbing system (a source of solid to the crust). Adding equations (2) and (3) gives conservation of mass in the crust-mantle system where we have assumed that the total fraction of Io occupied by the plumbing system is negligible. The flux of material in the plumbing system q p increases when material is extracted from the mantle and decreases when material is emplaced back into the crust. Conservation of mass in the plumbing system is therefore given by We assume that the magma in the plumbing system is at the melting point T m and we parameterise the emplacement rate based on the temperature difference to the host material at temperature T , where c is the specific heat capacity, L is the latent heat, h is an emplacement rate constant (units s −1 ) and T e is the elastic-limit temperature. Regions of the crust colder than the elastic limit temperature are assumed to be brittle; magma propagates through these without permanent emplacement. This assumption is introduced to highlight and investigate the importance of the distribution of magmatic emplacement in controlling the crustal temperature profile. The detailed mechanisms of dike propagation and emplacement are subsumed in this parametrisation. Although h could plausibly be related to heat conduction in the vicinity of a dike, we treat it here as a free parameter and explore the model's behaviour for a wide range of values.
Extraction of liquid from the mantle into the plumbing system is expected to take place at the top of the mantle. We assume that this transfer is a function of liquid overpressure, where ν is an extraction rate constant (units s −1 Pa −1 ) and P c is a critical overpressure that liquid must attain in order to be extracted into the plumbing system. Liquid overpressure (compaction pressure) is related to the compaction rate ∇· u by the relationship (McKenzie, 1984) where ζ = η/φ is the compaction viscosity, related to the shear viscosity η. This form of the compaction viscosity is commonly assumed but other forms have been proposed with a weaker singularity as φ → 0 (Rudge, 2018).
We model heat transport in the crust and mantle of Io together, using an enthalpy method (Katz, 2008) so that no boundary conditions need be imposed an the crust-mantle boundary. Conservation of energy where bulk enthalpy is defined as H = T +Lφ/c, and κ is the thermal diffusivity. Changes in bulk enthalpy are caused by advection of sensible heat, advection of latent heat, diffusion, tidal heating, extraction of melt, and emplacement of melt, which are represented respectively by each term in equation (9). The integral of the tidal heating rate ψ over silicate Io gives the total tidal heating input Ψ.

Solution methods
The full model to be solved comprises the enthalpy equation (9) for H, the combination of equations (1) and (8) in a compaction equation for P and q, total mass conservation (4) for u, and conservation of mass in the plumbing system (5) for q p . The enthalpy solution gives φ and T through the definition of bulk enthalpy H = T + Lφ/c, by assuming that temperature is buffered to the melting point when enthalpy exceeds the melting temperature, while melt fraction is zero wherever temperature falls below the melting point. The system is scaled (see appendix A) and solved using the Portable, Extensible Toolkit for Scientific Computation (PETSc) (Balay et al., 2019a(Balay et al., ,b, 1997. Robust convergence is obtained by splitting the system of governing equations into three non-linear problems for enthalpy, pressure, and plumbing-system flux. These are solved iteratively at each timestep until the L 2 -norm of their residual vectors are all below a small tolerance (10 −7 ). We run the model to steady state and report the final, steady solutions. To facilitate exploration of the parameter space, an asymptotic approximation of these solutions is also employed. This is developed in appendix B.

Results
A representative solution of the model is plotted in figure 2. Panel (a) shows that melt upwells throughout the mantle and the solid correspondingly downwells. This process, driven by magmatic buoyancy, results in relatively low melt fractions (∼ 3%) except in a thin boundary layer beneath the crust that we discuss below. Magma that reaches the surface solidifies and cools to the surface temperature ( The lower crust (T > T e ) is heated by magmatic intrusions (emplacement), but the upper crust (T < T e ) is cold due to the downwelling cold surface material. c) Melt fractions are low throughout the mantle but increase in a thin boundary layer on the order of 50 km thick. d) Compaction occurs throughout most of mantle as the liquid is at low pressure, but beneath the crust P P c = 0.8 MPa, and downwelling solid is decompacted by liquid pressure. The elastic thickness is 80 km and the eruption rate is 1.1 cm/yr, with 99.5% of heat transport through the surface being volcanic. Parameter values can be found in table 1 in appendix A.
of the downwelling crustal material.
We now explore the behaviour of the model in relation to various parameters. We first investigate the effect of the parameters associated with the volcanic plumbing system, before turning to material and rheological parameters.

Dependence on volcanic plumbing system parameters
The main parameters that control the behaviour of the volcanic plumbing system are the emplacement rate constant h, and the elastic limit temperature T e . The behaviour of the system for three values of emplacement rate constant h and elastic-limit temperature T e is shown in figures 3-4. The solid lines are the full, numerical solutions to the model and the dashed lines are from the asymptotic approximations (see appendix B). Figure 3 demonstrates that the rate of magmatic emplacement h exerts a strong control on Io's crustal thickness. When the emplacement rate is zero ( fig. 3, h = 0), magma does not re-heat the sinking crust and therefore the cold surface material downwells far into the mantle before it is heated to its melting point (by basal conduction and tidal heating only), producing a > 600 km thick crust. Conversely, when emplacement is rapid, the crust is very quickly heated to its melting point and therefore is thin. Figure   3e shows how the volcanic plumbing flux changes through the crust. When there is no emplacement, the drop in volcanic plumbing flux is purely due to radial spreading. Figure 4 shows the effect of increasing the elastic-limit temperature T e . When T e is low, emplacement can take place throughout a large portion of the crust. This causes most of the crust to be hot (fig 4, T e = 825 K) and thus the elastic thickness is small. Increasing T e leads to the growth of a large, cold upper crust where no emplacement is taking place. This provides a large elastic thickness capable of supporting Io's mountains. Figure 4e shows that the total amount of material being emplaced at steady-state does not significantly change as T e is increased, and the thickness of the emplacement region (T > T e ) is roughly constant as T e increases. As discussed in section 4.1, the total amount of erupted and emplaced material is controlled by a global energy balance.
To comprehensively map the parameter space of T e and h we rely on our asymptotic approximation of the steady solution. Figures 3-4 show that there is good agreement between the full solutions and the asymptotic approximation. Figure 5 shows how (a) elastic thickness, (b) eruption rate, and (c) volcanic heat flux vary as a function of the emplacement rate constant h and elastic-limit temperature T e . Figure   5 confirms the trends seen in figures 3 and 4, and places them in the context of observable features. The blue region in figure 5 indicates the parameter space that gives reasonable elastic thicknesses (10−100 km) at reasonable brittle-ductile transition temperatures (homologous temperature 0.5 − 0.7 T m ). Figure 5a shows that the elastic thickness varies rapidly with relatively small changes in T e and h around the main solution (the central star, plotted in fig. 2). Elastic thickness is thus the most useful observation for constraining the characteristics of Io's volcanic plumbing system. Figure 5b shows that eruption rate reaches a maximum of ∼ 1.25 cm/yr over much of the parameter space, a value that is discussed further below. The conductive heat flux in this part of the parameter space is negligible ( fig. 5c) as virtually all of the input tidal heating is lost in eruptions. Figure 5b shows that if emplacement rate is very high and takes place through the majority of the crust, eruption rate goes to zero. All heat is lost by conduction through a thin lid in this case. When there is no magmatic emplacement (h = 0), the crust grows to be over 600 km thick due to the rapid downwelling of the cold surface temperature. As emplacement rate is increased, the heating that this provides to the crust can increasingly balance the cold downwelling surface temperature, resulting in smaller crustal thicknesses. T e = 960 K in these solutions. Liquid flux is the sum of the segregation flux q and plumbing-system flux q p . Parameters values can be found in table 1 in appendix A.

Dependence on material and rheological parameters
We now examine the effects of varying material and rheological parameters in the model. Figure 6 shows solutions of porosity and pressure for a range of critical compaction pressures P c , and figure 7 shows solutions of porosity for a range of matrix shear viscosities η and liquid visosities η l . We consider the critical overpressure P c to be a material parameter as it parameterises the strength of the downwelling crust.
Pressure differences between the solid and liquid in the partially molten rock causes the solid matrix to deform. When the solid pressure is higher, the matrix compacts, expelling liquid, and when the liquid pressure is higher, the solid decompacts and melt accumulates. Away from boundary layers, the melt fraction is almost entirely controlled by the buoyant segregation of magma. Figures 6 and 7a show that P c and η do not affect the melt fraction outside of boundary layers, while figure 7b shows that increasing the liquid viscosity η l increases melt fractions throughout the mantle. This is because the buoyancy-driven melt fraction is controlled by the permeability and liquid viscosity in Darcy's law (1).
As solid crust downwells, warms, and begins to melt, the high pressure of rising magma forces the solid matrix to decompact to accommodate infiltration of buoyant magma. This decompaction occurs over a region known as a decompacting boundary layer, in which the compaction pressure gradient term in Darcy's law (1) becomes important. Figure 6 shows that if a large compaction pressure (liquid overpressure) is required for magma to move out of the mantle pore-space into the plumbing system, large porosities build up beneath the crust. This is because the large liquid overpressure drives rapid decompaction of the downwelling solid material. High shear viscosities also lead to the development of a larger boundary layer with higher peak porosity ( fig. 7a). An increased shear viscosity (and so through ζ = η/φ, an increased bulk viscosity) causes greater resistance of the downwelling matrix to compaction, which means a larger length scale over which the compaction pressure gradient counteracts buoyancy, and in which the porosity must therefore increase to enable upward melt motion. The thickness of the boundary layer is on the order of the compaction length, which is an emergent length-scale governing the interaction of liquid and solid in a two-phase medium (see appendix B) (McKenzie, 1984).
The interaction between compaction pressure driving melt flow and compaction sets up an oscillation of the porosity and pressure (Spiegelman, 1993) that decays to match the buoyancy-dominated mantle region below. Figure 7 shows that high shear viscosities and low magma viscosities (both of which represent high bulk viscosities) damp porosity oscillations effectively so that porosity rapidly relaxes to the buoyancy-driven profile.
In figures 6 and 7 the location of the crust-mantle boundary and the flux of melt out of the mantle are not affected by the parameter changes (where porosity increases, melt velocity decreases and so the flux is unchanged). As such, the material parameters control the form of the decompacting boundary layer but do not affect the crust or volcanic plumbing system. In particular, and perhaps surprisingly, the value of P c , which encodes information about the strength of the crust in this model, does not significantly affect the crustal thickness.

Discussion
Our results demonstrate the importance of magmatic intrusions (emplacement) in controlling Io's crustal thickness and temperature profile. A significant proportion of magmas generated in the mantle must contribute to heating the cold downwelling crust -otherwise it would be extremely thick. Further, the intrusions must be concentrated in the lower crust if the upper crust is to retain significant elastic strength, which it requires to support Io's high mountains (McKinnon et al., 2001). The results also   suggest that a high-melt-fraction layer can develop beneath Io's crust due to decompaction, without requiring any radial partitioning of tidal heat. We now discuss these results further and consider their implications.

Magmatic intrusions
The foremost result of this work is captured in figure 5, which indicates that magmatic intrusions are the primary control on Io's crustal thickness. Models of heat piping that neglect heating from intrusions will result in steady-state crustal thicknesses of over 600 km. Models with a downwelling crust that fix the crustal thickness cannot be assumed to be in thermal steady state and may violate energy conservation.
The primary control on the temperature profile in Io's crust is the distribution of magmatic emplacement.
Unless intrusions are confined to the lower crust, the upper crust becomes hot and weak, so would be unable to support significant topography. With a sufficiently large elastic-limit temperature, colder regions (in the upper crust) have no emplacement and so remain strong. Kirchoff et al. (2020) note that the release of confining stress by crustal faults leads to extension in Io's upper crust, which manifest as rifts, pull-apart basins, and simple graben structures. This transition to an extensional regime may further explain the low level of emplacement that must be taking place in the upper crust. A more detailed description of heat piping might attempt to account for the mechanics and energetics of emplacement in the context of the crustal stress profile.
Unless the crust is unrealistically thin, our model shows that conductive heat loss at the surface is negligible ( fig. 5c), consistent with the conclusions of O'Reilly & Davies (1981). As shown in more detain in appendix C, the heat loss due to eruption is q s (ρL + ρc(T m − T s )) per unit surface area, where q s is the globally averaged eruption rate, ρL is the latent heat of erupted magma, and ρc(T m − T s ) is the sensible heat lost as erupted magma cools to the surface temperature T s . At steady state, with negligible conduction through the elastic crust, this surface heat flux must balance the total dissipation rate Ψ, implying a volcanic resurfacing rate For Io, equation (10) predicts a resurfacing rate of 1.25 cm/yr (Breuer & Moore, 2015); this is the maximum rate seen in figure 5b. Equation (10) also provides a means of estimating eruption rates for other tidally heated lava-worlds, utilising their tidal heating rate, size, and surface temperature, all of which are obtainable from transit observations (Bolmont et al., 2013). The predicted eruption rate can be compared to an estimate of the total melt production rate. Assuming that all tidal heating directly causes melting, the total melt production rate is Ψ/ρL. The predicted eruption rate is therefore less than the melt production rate by a factor of c(T m − T s )/(L + c(T m − T s )). For Io, this indicates that around 80% of magma must be emplaced into the crust. This analysis demonstrates that the eruption rate and the total amount of emplaced material are controlled by the tidal heating rate and the relative temperatures of the magma and the surface. This illuminates why these quantities remain almost unchanged in figures 3 and 4, despite significant differences in the thickness of the crust and the precise location at which emplacement occurs.

Figures 6 and 7 demonstrate that if a high-melt-fraction layer exists within Io, it does not necessarily
imply the presence of a magma ocean with melt fractions above disaggregation (Khurana et al., 2011;Tyler et al., 2015), nor was it necessarily formed by concentrated tidal heating (Moore, 2001;Bierson & Nimmo, 2016). Instead, our model suggests that a high-melt-fraction layer could form if a large liquid overpressure is required in order to inject dikes into the crust, or if the compaction length is large. McKinnon et al. (2001) and Kirchoff & McKinnon (2009) note that high compressive stresses must arise in Io's crust due to the downwelling of a spherical shell. If these compressive stresses are present at the crust-mantle boundary, they will have to be overcome by magma pressure to form dikes; this would indicate a high value for P c , promoting the formation of a high-melt-fraction region beneath the crust.
However, we have demonstrated that Io's lower crust must be hot due to extensive magmatic intrusion, and so it is likely that the stresses associated with downwelling crust would be accommodated by faulting in the upper crust and viscous creep in the lower crust. The relevant value for P c in the model is thus unclear; more investigation is needed into the strength and deformation mechanisms of Io's lower crust.
Nevertheless, we have demonstrated that a high-melt-fraction layer can in principle arise from a state of uniform tidal heating. Further work is needed to ascertain whether such a layer is likely to exist.
In particular, we must determine whether the required high shear viscosities are compatible with the observed tidal heating (Bierson & Nimmo, 2016;Renaud & Henning, 2018). If such a layer does not exist, our model suggests that melt fractions are relatively uniform within Io ( fig. 7), providing little drive for radial partitioning of tidal heating.
We emphasise that the presence and form of the decompacting boundary layer are independent of the model predictions for crustal thickness. In the context of our model, a decompacting boundary layer is a feature beneath the crust that does not affect the mass or heat fluxes out of the mantle, nor how energy is deposited in the crust; it is these other factors that control Io's crustal thickness.

Model limitations
A primary simplification made in formulating our model is the assumption of spherical symmetry. Other models suggest that tidal heating does not just vary with radius but also with latitude and longitude. Steinke et al. (2020) show that if dissipation is significantly concentrated in a layer beneath the crust, lateral variations in mantle temperature can exceed 100 K, perhaps impacting the usefulness of a onedimensional approach. More distributed heating leads to much lower lateral mantle temperature variations (∼ 1 K, Steinke et al. (2020)), but may still cause significant lateral variations in melting rate. The presence or absence of a high-melt-fraction layer is thus key to understanding Io's three-dimensional tidal heating distribution.
Convection of partially molten mantle is also a potentially important three-dimensional effect. If buoyant segregation is inefficient at transporting heat, convection of the two phase medium could occur and imprint deviations from the one-dimensional model presented here. However, we stress that convection cannot dominate Io's mantle heat transport due to its low efficiency at melt fractions relevant for Io (Breuer & Moore, 2015). Also neglected here is that compaction effects can cause the lateral migration and focusing of melt (Sparks & Parmentier, 1991;Turner et al., 2017), which could lead to channelisation of melt and may exert a control on Io's volcano distribution.
Another key simplification in this model is the assumption that Io is composed of a single chemical component. Continued melting in the interior is likely to have caused significant chemical stratification (Keszthelyi & McEwen, 1997). Melting of a polymineralic rock occurs over a range of temperatures, with more fusible minerals melting first. The upward migration and eruption of fusible melts plausibly depletes the deep mantle of fusible material and enriches the near-surface. Erupted fusible material would also melt at shallow depths upon burial, affecting the crustal thickness. The dynamics of this segregation is omitted here and may exert important controls on Io's volcanism and mantle melt distribution.

Conclusions
Io is a body that is complicated in its detail but observations suggest it has a simple structure at leading order. We have demonstrated that a coupled model of magmatic segregation, compaction, and heatpiping can explain this leading-order structure and the associated observations: globally averaged elastic thickness, eruption rate, and surface heat flux, as well as a possible high-melt-fraction layer beneath the crust. We have shown that magmatic intrusions into Io's crust are a fundamental control on its crustal thickness. Without the heating associated with the formation of magmatic intrusions, the crust would grow to be > 600 km thick. However, these intrusions must be confined to the lower crust if the upper crust is to retain sufficient strength to support Io's high mountains. We have also shown that an inferred high-melt-fraction region can be understood as a decompacting boundary layer if a process such as lateral compression makes it difficult for magma to migrate from the mantle into the crust. An extension of our model to include more elaborate chemical thermodynamics and three-dimensional flows will give insights into deviations from the spherically symmetric model developed here.

A Scaled model and parameter values
Here we non-dimensionalise the governing equations. Parameter definitions are given in table 1 and the scales and definitions of the non-dimensional parameters are given in table 2. We write for example u = u 0û , where u 0 is the velocity scale andû is the dimensionless velocity, insert similar expressions for all the variables into the equations, and finally drop the hats on the dimensionless quantities to arrive at a dimensionless model. For the temperature we write T = T s + T 0T with T 0 = T m − T s so that the non-dimensional temperature varies between 0 and 1. We also assume spherical symmetry, and write all quantities as a function of r, noting that ∇· u = r −2 ∂(r 2 u)/∂r where u is the radial component of the solid velocity.  q 0 q 0 = ψ 0 R/ρL 6.4 × 10 −9 m/s Solid velocity scale u 0 u 0 = q 0 6.4 × 10 −9 m/s Porosity scale Compaction parameter δ δ = ζ 0 K 0 φ n 0 /η l R 2 5.8 × 10 −3 The tidal heating scale ψ 0 is imposed, which gives the velocity scale q 0 which in turn gives the porosity scale φ 0 .
The non-dimensional equations for conservation of solid and liquid mass are where the emplacement rate M =ĥ(1 − T )I M . I M is an indicator function that equals 1 for T > T e and q p > 0, and equals zero otherwise. This ensures that emplacement only occurs above the elastic limit temperature, and provided there is melt present in the plumbing system to be emplaced. E =ν(P −P c )I E is the extraction rate, where the indicator function I E equals 1 for P > P c , and equals zero otherwise, ensuring that extraction only occurs in regions above the critical overpressure.
Total conservation of mass for the crust-mantle and plumbing system are where δ is a dimensionless parameter defined in table 2. This measures the typical size of the compaction pressure gradients relative to the buoyancy; it is expected to be relatively small. It can be related to the compaction length (McKenzie, 1984) l = ζ 0 K 0 φ n 0 /η l , by δ = l 2 /R 2 (so the square root of δ is the ratio of the compaction length to the radius of the planet).

B Asymptotic Approximation
To facilitate a rapid exploration of parameter space, we construct an approximation to the steady states of the model. This approximation makes use of the fact that the porosity scale φ 0 and the compaction parameter δ are both much less than unity. Neglecting them in the equations provides a good approximation over most of the crust and the mantle, apart from in the boundary layer just below the crust-mantle boundary (which is discussed below). Note that Pe −1 is also expected to be a small parameter, but we retain conduction in the equations because it is important in controlling the temperature distribution within the crust.

B.1 Mantle and Crust
In this approximation we consider the crust and mantle separately and solve for the position of the boundary r = r c between them. We assume that all extraction from the mantle occurs within the decompacting boundary layer, just below the base of the crust (as is verified by the full numerical solutions). The extraction term E is therefore non-zero only in a narrow region of thickness O(δ) and hence is neglected from the continuum equations. At the base of the crust r = r c , the entire flux q(r c ) (hereafter referred to as q c ) is transferred into the plumbing system. As q p = 0 in the mantle, emplacement M only appears in the crustal equations. Taking this into account, the combination of mass continuity equations (11) and (12) indicate that u = −q throughout the mantle. Since we also have T = 1 in the mantle, conservation of energy (eqn. 16) becomes simply so melting (or equivalently, the transport of latent heat by melt) balances tidal heating. Equation (17) can be integrated directly to give Darcy's law (eqn. 15a) becomes q = φ n so φ is deduced directly from q, and the compaction equation (15b) further indicates that P = −ψ/φ. In particular, these expressions give the values for the flux q c , porosity, and compaction pressure at the top of the mantle r = r c , which are used below to feed into the decompacting boundary layer, The flux q c transfers to the plumbing system at the top of the mantle for its continued transport through the crust, to which we now turn.
In the crust, q = 0 and the combination of mass equations (13) and (14) now requires u = −q p , the plumbing-system flux. Conservation of mass in the plumbing system (eqn. 14) becomes Conservation of energy in the crust (eqn. 16) becomes so advection of the solid balances conduction, tidal heating, and heating from intrusions (emplacement).
These two equations ( (20) and (21)) are solved together to determine the temperature profile T , the plumbing system flux q p (and hence the solid velocity u), and the position of the crust-mantle boundary r c . The required boundary conditions are q p = q c , T = 1, ∂T ∂r = 0, at r = r c , T = 0, at r = 1.
Although this solution to the crustal system involves a numerical integration, it is considerably more straightforward and faster than the solution to the full model.
From this approach we find a good approximation to the thickness of the crust, the temperature profile within the crust, the plumbing flux q p and emplacement rate (and hence the eruptive flux at the surface), as well as the porous melt flux and porosity in the majority of the mantle. A detail missing from the full solutions that is not yet captured by this asymptotic approximation is the high-porosity region -the decompacting boundary layer -just below the base of the crust. In the context of our model, apart from transferring melt from the porous mantle to the plumbing system, the details of this layer are unimportant in determining the large scale structure of the solutions (thickness and temperature distribution of the crust). However, in order to understand the dynamics further, we now analyse this region.

B.2 Decompacting Boundary Layer
The behaviour in the boundary layer is obtained by rescaling the equations locally to find a local approximation of the solution close to the crust-mantle boundary. A composite approximation valid over the whole domain can then be found by combining the two approximations (the 'outer' mantle and crust solution given above, and the 'inner' boundary layer solution).
As rising magma approaches the crust-mantle boundary at r = r c , it is no longer reasonable to neglect the compaction pressure gradient term in equation (15). Approaching the crust, rising magma in impeded by the low permeability of downwelling solid, which causes magma to accumulate. The accumulation of magma in this layer generates pressure that decompacts the low-porosity-solid that is downwelling from the crust (cf. Hewitt & Fowler (2008)). To understand what happens in this boundary layer, we must reintroduce the term proportional to δ in equation (15), and rescale lengths to consider the dynamics close to the crust-mantle boundary.
First we note that including extraction, the conservation of mass equations for the solid, liquid, and plumbing system in the mantle are We also note that the compaction pressure relation (eqn. 15b) can be combined with conservation of solid mass in the mantle (eqn. 23a), giving We then write r = r c −δZ, where Z is a boundary layer coordinate describing the rescaled distance beneath the crust-mantle boundary. Rewriting equations (23), (15a), and (24) in terms of this coordinate, we where O(φ 0 , δ) represents terms of order φ 0 or δ, which may be neglected (note that having rescaled into the boundary layer, the neglected terms are not the same as those neglected earlier, reflecting the different dominant physics in the boundary layer). We write φ 0 /δ = µ, and treat this as an O(1) parameter.
The first mass conservation equation (25a) indicates that u is approximately constant throughout this layer, and its value is given by matching the value in the mantle below, u = −q c , where q c is the liquid flux defined in equation (19). Moreover, adding together the second two mass conservation equations (25b-c) shows that q + q p is constant, and must be equal to q c to match with the mantle. The final two equations (25d-e) (which represent Darcy's law and the compaction relation) can therefore be written as Extraction occurs over the region 0 < Z < Z E (Z E is determined shortly). In this region we have E =νδ(P − P c ), and we take the limit of very large extraction rate constant ν (soνδ is large), such that P P c throughout this region. Equation (26a) can then be integrated to give and as P is approximately constant, equation (26b) gives q = φ n . The porosity and flux q therefore increase with Z through this extraction region and since q p = q c − q, the plumbing flux correspondingly decreases with Z until it reaches 0. This defines the position Z E at which extraction started. Substituting q c = q = φ n into equation (27) gives Turning to the region below extraction where E = 0, mass conservation equations (23b-c) indicate that q p = 0 and q = q c are now approximately constant. In this region equations (26a) and (26b) comprise a two-dimensional phase-plane problem for φ(Z) and P (Z). A solution is sought with P = P c and φ = q 1/n c at Z E (for continuity with the extraction region), and that matches the correct far-field behaviour as Z → ∞. The correct behaviour is that φ and P tend towards the values q 1/n c and −ψ/q 1/n c , given earlier in equation (19), to match with the rest of the mantle. Since this corresponds to a fixed point of the system that is a stable spiral or node, such a solution can be found. The solution involves decaying oscillations of both φ and P towards the far-field values, which are evident in figures 3 and 4, and even more so in figures 6 and 7.
It it worth pointing out that the 1/φ dependence of bulk viscosity is an important control on the porosity oscillation within this boundary layer. Other forms of bulk viscosity with a weaker singularity have also been suggested; for example with ζ ∼ − ln φ as φ → 0 (Rudge, 2018). The weaker dependence on porosity leads to greater oscillations, but the general form of the boundary layer is maintained.
To compare the asymptotic approximation with the full numerical solution in figures 3 and 4, we construct a combination of the 'outer' solution for the majority of the mantle (where q is given by equation (18), φ = q 1/n , and P = ψ/φ), together with the 'inner' solutions for the decompacting boundary layer.
Denoting the former solution as φ m (r) and the latter solution φ bl (Z), this composite solution is defined by with equivalent expressions for P and q, where the subtraction of the 'overlapping' value φ l is necessary to avoid double counting. The solutions so obtained are re-dimensionalised to produce the approximate solution that is shown as the dashed lines in figures 3 and 4, which is in good agreement with the full solutions.

C Analysis of heat flux and emplacement
Useful information can be obtained by integrating the energy equation (21) over the crust (from r = r c to r = 1). Substituting equation (20) into equation (21), recalling that u = −q p and integrating, gives where q c is the melt flux at the top of the mantle, which was defined in equation (19). Equation (30) where we recall that r b is the radius of the core. The left hand side here represents the heat loss due to eruption and conduction at the surface, and the right hand side is the total tidal input; this expression thus represents a global energy balance.
After re-dimensionalising the variables, this can be written as where k is the conductivity.
As shown in figure 5, for the parameter regime applicable to Io, the surface conductive heat flux is negligible. In this case we see that the resurfacing rate q s (that is, the value of q p at the surface) is given by q s = Ψ 4πR 2 (ρL + ρc(T m − T s )) .
close to R, that expression can be written dimensionally as q c ≈ Ψ 4πR 2 ρL .
As such, q s is approximately a fraction L/[L + c(T m − T s )] of q c , and the proportion of q c that must be emplaced is For Io, equation (33) gives a resurfacing rate of of 1.25 cm/yr. Compared to the flux into the base of the crust (eqn. 34) of 6.3 cm/yr, equation (35) predicts that 80% of magma produced inside Io is emplaced into the crust.