Compositional layering in Io driven by magmatic segregation and volcanism

Magmatic segregation and volcanic eruptions transport tidal heat from Io's interior to its surface. Several observed eruptions appear to be extremely high temperature ($\geq$ 1600 K), suggesting either very high degrees of melting, refractory source regions, or large amounts of viscous heating on ascent. To address this ambiguity, we develop a model that couples crust and mantle dynamics to a simple compositional system. We analyse the model to investigate chemical structure and evolution. We demonstrate that magmatic segregation and volcanic eruptions lead to differentiation of the mantle, the extent of which depends on how easily high temperature melts from the more refractory lower mantle can migrate upwards. We propose that Io's highest temperature eruptions originate from this lower mantle region, and that such eruptions act to limit the degree of compositional differentiation.


Introduction
Jupiter's moon Io is the most volcanically active body in the solar system. Its volcanism is a result of tidal heating from its mean motion resonance with Europa and Ganymede, which causes widespread melting in its interior (Peale et al., 1979;O'Reilly & Davies, 1981). Despite its long history of study, it is not well known to what extent melting and volcanism control Io's interior structure and evolution and, in particular, if these processes create compositional layering within the mantle. The best constraints on interior structure would be provided by measurements of the composition and temperature of erupted lavas. To keep pace with recent improvements in observational techniques (e.g. Davies et al. (2016Davies et al. ( , 2017; de Kleer, de Pater, et al. (2019); de Kleer, Nimmo, and Kite (2019)), interior evolution models that are predictive of eruption temperatures and compositions are increasingly required. Keszthelyi and McEwen (1997) presented an initial attempt to estimate the geochemical and petrological structure of Io's interior that would arise from the extensive volcanism. They predicted that the crust would be dominated by felsic lavas rich in incompatible elements and that the mantle would be dominantly a forsterite-rich dunite. When the initial Galileo observations suggested widespread eruption of ultramafic lavas and constrained the temperature of the Pillan eruption to 1870±25 K (McEwen et al., 1998, this model was abandoned. It was replaced by a model that called upon a region with ∼ 50% partial melting at the base of the crust. This configuration hypothetically allowed efficient recycling of the erupted lavas back into the mantle (Keszthelyi et al., , 2004. This magma-ocean model was supported by Galileo magnetometer results (Khurana et al., 2011) and is consistent with the suggestion of magnesian orthopyroxenes in Ionian lavas (Geissler et al., 1999). The magma-ocean model predicts a well-mixed and geochemically homogeneous mantle (Keszthelyi et al., 2004); erupted lavas would be largely uniform in temperature and composition, most likely similar to terrestrial komatiites (Williams et al., 2000).
However, there were significant challenges to the magma-ocean model as proposed in Keszthelyi et al. (2004). For example, once partial melting exceeds ∼ 20%, the shear modulus drops to the point that tidal dissipation cannot match the surface heat flow (Moore, 2003;Bierson & Nimmo, 2016;Renaud & Henning, 2018), limiting the possible thickness of such a high-melt-fraction layer. Furthermore, applying a different thermal model to the Pillan eruption, its temperature was revised down to ∼ 1600 K (Keszthelyi et al., 2007). Indeed, even the initial McEwen et al. (1998) results showed most eruptions being consistent with ∼ 1300 K (i.e., basaltic) temperatures. Spectroscopic constraints on the mineralogy of Io's lavas were always known to be weak because the Galileo camera did not observe far enough into the infrared to reliably detect other key minerals such as olivine (Geissler et al., 1999). These issues led to a revised magma-ocean model with the maximum degree of mantle partial melting only reaching ∼ 25% and decreasing rapidly with depth (Keszthelyi et al., 2007). Auroral hotspot oscillations have been used as evidence against a magma-ocean (Roth et al., 2017), and reanalysis of the magnetometer results suggests that plasma interactions with the atmosphere provide an alternative explanation to a magma ocean (Blöcker et al., 2018;de Kleer, McEwen, & Park, 2019). More recently,  showed that high melt fractions can arise within a decompacting boundary layer at the top of a low-melt-fraction mantle. Indeed, the distinction between a magma-ocean model and a low-melt-fraction model has significantly reduced since Keszthelyi and McEwen (1997) and McEwen et al. (1998); at this point, the hypothesis that Io is a largely solid body that has undergone significant magmatic differentiation needs to be investigated.
In this work we present a fluid dynamical model of crust and mantle dynamics that builds on the recent work of  by including compositional evolution. The compositional model is in the form of a two-component phase diagram between hypothetical refractory and fusible components. We use this simplified theory to investigate the effect of magmatic segregation and volcanic eruptions on leading-order chemical structure. Our results show that magmatic segregation causes a rapid differentiation of the mantle, with fusible material in the upper mantle and crust, and refractory material at depth. Magma forms in both the upper and lower mantle and, importantly, magma must be able to leave the lower mantle in order to facilitate heat loss. The model exhibits two distinct modes of behaviour, depending on the fate of magma produced in the lower mantle. If lower mantle melts stall within the upper mantle, high temperature eruptions should not occur. However, if these refractory melts migrate to the surface, they can provide an explanation for the highest temperature eruptions observed on Io.
The manuscript is organised as follows. First we outline the physics of the model before presenting results showing the two distinct modes of behaviour. We demonstrate the time evolution of both modes, and investigate the effect of bulk composition on the system. We then discuss these results in the context of present and potential future observations.

Model description
The model, shown schematically in figure 1, considers the evolution and dynamics of a tidally heated body composed of a mixture of two chemical components. It is an extension of that described in  using the same equations. Here it is extended to consider conservation of chemical species and the effect of composition on melting behaviour, using a phase diagram described below. We consider the crust and mantle to be a continuum that can either be entirely solid or partially molten, depending on the local energy content, and solve a system of conservation equations for mass, momentum, energy, and chemical species.
Alongside the continuum, we model a magmatic plumbing system that provides a means of upward magma transport distinct from magmatic segregation. Keszthelyi and McEwen (1997) proposed that deep, refractory magmas may sometimes ascend to the surface from great depth, but a mechanism to allow this has not been explored. We assume that anywhere magma reaches high overpressure, it enters into a magmatic plumbing system and migrates upward; this system can be present in both the mantle and the crust. Possible physical interpretations of this plumbing system will be considered in the discussion section. When magma enters the plumbing system, it transports the local melt composition and temperature upward into the upper mantle and crust. The flux of melt that reaches the surface is the erupted flux and its composition sets the composition of the newly resurfaced crust.
We revisit the thermochemical melting models that have been used to predict the segregation of Io's mantle into an upper fusible layer and a deep layer of almost-pure olivine (Keszthelyi & McEwen, 1997). Our approach is to simplify the compositional model to two representative end-members, aiding their incorporation into a dynamical framework. We consider Io to be composed of a mixture of these two components, the melting behaviour of which is described by the two-component phase diagram shown in figure 2. The presence of fusible material (component A) significantly reduces the melting point of the refractory component (component B), and so upon heating, fusible melts are produced until component A is almost entirely removed from the system. These types of compositional model have proven fruitful in studies of mantle melting at mid-ocean ridges (Katz, 2010;Katz & Weatherley, 2012).
As in , we assume spherical symmetry motivated by the global distribution of Io's volcanoes (Kirchoff et al., 2011;Williams et al., 2011). We focus our analysis on the chemical evolution of the system, and therefore take tidal dissipation to be uniform, avoiding dependence on poorly constrained rheological parameters (Bierson & Nimmo, 2016;Renaud & Henning, 2018). We neglect the pressure-dependence of the melting temperature due to the small size of Io and hence the low pressures in the mantle. However, for more detailed petrologic modeling, it may be important to include this effect (e.g., Keszthelyi et al. (2007)).
Our model considers the time-dependent evolution of the interior structure and composition, and explores the evolution to a steady state. We also develop a reduced model to elucidate key features of the dynamics predicted by the full model. The reduced model is formulated at steady state and its structure is motivated by solutions obtained to the full model; it is detailed in Appendix C.

Model equations
We consider a generic refractory component B and a fusible component A, and the phase diagram shown in figure 2. The concentration of the fusible component A in phase i (solid s or liquid l) is denoted c i , and that of the refractory component is 1−c i . The solidus temperature T s is given by and the liquidus temperature T l is given by where T B is the melting point of the refractory component, T A is the melting point of the fusible component, and γ > 0 is a parameter that controls the amount of fusible material that is incorporated in a solid solution with component B. We allow this small degree of solid solution simply because it provides a smoothed solidus curve, which facilitates our numerical method (the effect of smoothing the solidus is small, and is dis--4-manuscript submitted to JGR: Planets Figure 1. Schematic of the model. Magma rises buoyantly in the mantle while the solid moves downwards. If a critical overpressure is exceeded, magma is extracted to a magmatic plumbing system. It freezes (is emplaced) from the plumbing system back into the continuum at a rate defined in equation (10). Some magma reaches the surface, fueling volcanic eruptions and burying the crust. The composition of erupted magma determines the composition of the crust.
The core is excluded from this model. cussed in Appendix C). As γ → 0, the smoothed solidus approaches that of pure refractory component B. The chosen form for the solidus should not be interpreted as representative of any underlying thermodynamics.
The model of  is described by conservation equations for mass, momentum, and energy in a compacting two-phase medium and conservation of mass and energy equations in the magmatic plumbing system. These are 1 ρC where u is the solid velocity, q = φ(v liquid −u) is the Darcy segregation flux, E is the extraction rate to the plumbing system, and M is the emplacement rate from the plumbing system. Porosity is denoted by φ, and K 0 φ n is the permeability, in which n is the permeability exponent. In addition, ∆ρ is the density difference between solid and liquid, g = −gr is the gravity vector, η l is the liquid viscosity, P = (1 − φ)(P liquid − P solid ) is the compaction pressure, and ζ = η/φ is the compaction viscosity related to shear viscosity η. Bulk enthalpy is defined as H = ρCT +ρLφ, T is temperature, L is the latent heat, C is the specific heat capacity, ρ is the density, ψ is the volumetric tidal heating rate, κ is the thermal diffusivity, T p is temperature in the plumbing system, and q p is the plumbing system flux.
Conservation of mass (3) tells us that material leaves the crust-mantle system by extraction to the plumbing system and enters the crust-mantle system by emplacement from the plumbing system back into the continuum. We note that "emplacement" may have different interpretations in other works, but here it simply means the arrest and freezing of rising plumbing-system melts within the interior. Conservation of momentum is formulated by the combination of Darcy's law (4a), which tells us that fluid flow is driven by buoyancy and compaction pressure gradients, with the compaction relation (4b), which relates the liquid overpressure to the compaction rate ∇· u (McKenzie, 1984). Equation (4b) includes magmatic emplacement because we assume that emplacement does not cause fluid pressurisation. Conservation of energy (5) tells us that changes in bulk enthalpy occur by the advection of sensible and latent heat, diffusion of sensible heat, tidal heating, the energy removed by extraction, and the energy delivered by emplacement. We note that in  bulk enthalpy was normalised by the volumetric heat capacity ρC. Conservation of mass (6) in the plumbing system tells us that the plumbing system flux increases when material is extracted from the mantle and decreases when material is emplaced back into the continuum. Equation (7) represents conservation of energy in the plumbing system. There are no time derivatives in equations (6)-(7) because the plumbing system is assumed to occupy negligible volume.
To the equations above, we add an equation that tracks the composition of the system ∂c ∂t where c = φc l +(1−φ)c s is the phase averaged composition and c p is the composition of material in the plumbing system. This equation tells us that changes in phase averaged composition occur through advection of the liquid composition, advection of the solid composition, extraction of the liquid to the plumbing system, and emplacement of the plumbing system material. We neglect compositional diffusion due to the large advective velocities compared to chemical diffusivity. The composition of plumbing system material is given by a conservation of chemical mass equation where the plumbing system composition can only change by the addition of melts from the crust-mantle system of a different composition.
As in  we assume that the emplacement rate of magma from the plumbing system to the continuum is proportional to the temperature difference between the plumbing system material and the local continuum where T e is an elastic limit temperature below which no emplacement occurs . The emplacement rate constant h is discussed at length in , but here we propose that it may have different values in the mantle h M and the crust h C (the crust is where T < T A ). The mechanisms by which magma propagates through a partially-molten medium are likely to be very different to those in a solid, and so would be expected to have a different efficiency of magma transport. In this work, h C is directly analogous to h in Spencer,  and the behaviour with different values of h M will be explored.
Extraction of liquid from the mantle into the plumbing system is treated in the same way as in ; the transfer is taken to be a function of liquid overpressure, where ν is an extraction rate constant (units s −1 Pa −1 ), and P c is a critical overpressure that the liquid must exceed in order to be extracted into the plumbing system. We recall that P is the overpressure relative to the lithostatic pressure P solid , not the absolute liquid pressure P liquid .
The full model to be solved comprises equations (3)-(11), which govern the time evolution of temperature, porosity, and composition, as well as the magma and solid velocities. The phase averaged composition c and the bulk enthalpy H uniquely define the temperature, porosity, and liquid and solid compositions through the solidus and liquidus equations (1)-(2), the definition of bulk enthalpy, and the definition of phase averaged composition. The boundary conditions state that there is zero solid and liquid velocity and zero heat flux at the base of the mantle (r m in figure 1), and that there is a prescribed surface temperature T s . The composition at the surface is set by the erupted composition, which together with the zero basal fluxes, conserves the bulk composition. The bulk composition is therefore effectively set by the initial conditions.
Parameter values and definitions are given in table 1. The system is scaled (see Appendix A) and spherical symmetry is assumed so that all variables are a function of only radial position r and time. The system is solved using the Portable, Extensible Toolkit for Scientific computation (PETSc) (Balay et al., 1997(Balay et al., , 2019(Balay et al., , 2020Katz et al., 2007). Details of the implementation are given in Appendix B. The system is benchmarked against the single-chemical-component model in Spencer, Katz, and Hewitt (2020).

Results
The steady-state behaviour of the model across parameter space can be broadly divided into two distinct modes. This division is on the basis of the transport of refrac-tory melts that form in the lower mantle, which is controlled by the value of the mantle emplacement constant h M . The results in this section are framed to exhibit the contrasting behaviour of these two modes; the implications of each mode will be discussed further below. In mode 1, rising refractory magma in the magmatic plumbing system interacts and exchanges substantial energy with the lower-temperature partially-molten upper mantle. This drives all plumbing-system magmas to freeze within the upper mantle and, as a result, refractory melts to not reach the crust. In mode 2, refractory plumbingsystem magmas rise through the upper mantle with little to no interaction. These melts reach the base of the crust, combine with more fusible melts, and are erupted to the surface. Figures 3 and 4 show steady-state solutions for the full model for each of the two modes. Figure 5 shows the evolution of the model from an initial uniform state, again for each of the two modes. Finally, in figure 6 we summarise the behaviour of the model as a function of the bulk composition of the body, demonstrating the transition between the two modes. These figures are discussed further below.
In this paper we do not explore the parameter space of the crustal emplacement constant h C , the elastic limit temperature T e , nor the critical extraction pressure P c . The effect of variation in these parameters was considered by  and their effects here are the same. The crustal emplacement constant h C and the elastic limit temperature T e control the thickness and temperature distribution in the crust, and the critical extraction pressure P c affects the melt fraction in decompacting boundary layers that occur where magma is extracted to the plumbing system. In the results presented here, we choose values of h C and T e that give reasonable crustal thicknesses and temperature distributions. We take P c = 0 and explore whether compositional effects also exert a control on melt fractions. Figure 3 shows temperature, porosity, fluxes, and compositions at steady state for two representative values of h M . Refractory magmas that form in the lower mantle are transferred to the magmatic plumbing system at the top of the lower mantle, enabling their continued rise. As they rise through the upper mantle, they are emplaced at a rate proportional to h M , and it is the size of this parameter that distinguishes the two modes. Mode 1 arises when h M is sufficiently large that all the melt from the lower mantle is emplaced into the mid-and upper mantle. Mode 2 arises when some of the melt extracted from the lower mantle reaches the crust, which occurs if h M is sufficiently small. Solid lines in figure 3 are steady-state solutions to the full model; dashed lines are solutions to the reduced model (see Appendix C).

Two Modes of Magmatism
The two modes share various features that can be identified from figure 3. We discuss these similarities before considering their differences. Some features are similar to those in the one-component case of , which we cover only briefly here. The radial porosity profiles in figure 3b,f show that the uniform tidal heating causes melt to form throughout the mantle. Figure 3c,g shows that these melts rise buoyantly while the solid correspondingly sinks. Where melt reaches high pressure it is extracted into the plumbing system, through which it continues to rise. The crustal plumbing system carries melt to the surface where it erupts. The globally-averaged eruption rate is the surface plumbing-system flux in figure 3c,g. Over long timescales and given the negligible surface conduction, this global eruption rate must extract heat at the same rate that it is input to the interior by tidal heating. The upward flux of melt through the crustal magmatic plumbing system is balanced by downwelling of the solid crust. This recycles erupted material back into the mantle.
At steady state in both modes, the mantle has segregated into three layers: a refractory lower mantle with T = T B , a low-melt-fraction mid-mantle with T A < T < T B , and a fusible upper mantle with T ≈ T A . As crustal solid downwells through the -9-manuscript submitted to JGR: Planets upper mantle, tidal heating causes the formation of fusible melts, which buffers the temperature close to T A . With continued melting and the buoyant segregation of fusible melts, material downwelling out of the upper mantle is almost exhausted in fusible material and so its solidus temperature has increased according to the phase diagram. In this midmantle region, tidal heating primarily acts to raise the temperature of the solid. As a result, melting rate and porosity are low in the mid-mantle, as seen in both modes in figure 3b,f. Further, the Darcy flux in the mid-mantle is approximately zero (figure 3c,g), so heat transport across this region occurs only by conduction, advection in the plumbing system, and downward solid advection, a result that we discuss below. Continued heating as the solid downwells through the mid-mantle melts out the remaining small amount of fusible material, and the solid is raised to the refractory melting point T B . Melting rate and thus porosity increase in the lower mantle because, as in the upper mantle, all imparted tidal heating directly causes melting.
Magma rising through a two-phase medium cannot pass into impermeable regions. Such regions act as barriers to flow, causing an increase in magma pressure, which forces the solid to decompact and produces higher melt fractions (figure 3b,f). The crust represents such an impermeable barrier to melts rising through the upper mantle, and similarly, the mid-mantle region acts as an essentially impermeable barrier to melts rising from the lower mantle. The high liquid pressure below these layers causes melts to be extracted into the magmatic plumbing system. Magma extracted from the lower mantle is composed entirely of the refractory component and is at temperature T B . Flow through the plumbing system enables these refractory melts to migrate from the lower mantle into the colder overlying mantle and crust. The differences between the two modes are then a consequence of what happens to this melt. The mid-and upper mantle are below the melting point of the refractory component, and it may be expected that these lower temperatures causes refractory plumbing system material to be emplaced during ascent.
In mode 1 (figure 3a-d), this emplacement is significant -it acts to exhaust the plumbing system of refractory material before it reaches the crust. As refractory melts are emplaced they release their latent heat to the upper mantle, providing additional heat to melt surrounding fusible material. This is reflected in the rapid increase of Darcy flux in the lower part of the upper mantle in figure 3c. The emplacement of refractory melts into the upper mantle eventually exhausts the material in the plumbing system, as shown by the plumbing system flux in 3c. Where the plumbing system material runs out, the melting rate in the upper mantle decreases to just that produced by tidal heating, which causes the change in gradient of the Darcy flux in the upper mantle in figure 3c. The change in melting rate caused by the cessation of emplacement means that downwelling solid must suddenly decompact, creating a high-porosity decompacting layer in the upper mantle, which can be seen in figure 3b.
Mode 2 (figure 3e-h) is the case where at least some of the melt that is extracted from the lower mantle makes it all the way to the surface. The end-member shown in figure 3 is when h M = 0, in which case there is no emplacement in the upper mantle at all. The plumbing-system flux still decreases in figure 3g, but only due to radial spreading in a spherical coordinate system, and so the total volume of melt extracted from the lower mantle reaches the top of the upper mantle. Fusible magmas extracted at the top of the upper mantle combine with refractory plumbing system melts rising from below, producing crustal plumbing-system material with a volumetrically averaged temperature and composition. This crustal plumbing-system material describes either an average of non-interacting melts of different temperatures and compositions, or a mixture with an intermediate composition; we assume that the effect is the same on the long timescales considered here. The crustal plumbing system melts are emplaced into the crust at a rate determined by h C and the temperature of the melt, and with a distribution determined by T e . Material that erupts onto the surface in mode 2 is at a higher temperature than in mode 1, and so serves as a more efficient heat-loss mechanism. This increased heat-loss efficiency results in a lower eruption rate and a thinner crust (see below). Figure 4 shows a schematic of temperature, mass transport, and the phase diagram. Colours in figure 4 denote composition according to the phase diagram in panel c. Mode 1 is characterised by a strong segregation of fusible and refractory material; refractory material does not erupt, instead it is cycled between the lower mantle and the deep parts of the upper mantle, whilst fusible material is cycled between the upper mantle and the crust. In mode 2, refractory material is cycled from the lower mantle to the surface, and fusible material is cycled from the upper mantle to the surface. In both modes, the lower mantle is composed purely of refractory material, and the mid-mantle spans compositions corresponding to the steep section of the solidus in figure 4c. In mode 1 there is a transition from almost pure refractory to pure fusible material above the region of the upper mantle where emplacement takes place (figure 3d). In mode 2, the segregation of the mantle is much less complete, as shown by figure 3h. The lack of mantle emplacement means that refractory melts rise all the way to the surface. The intermediate-composition erupted material is buried down through the crust and upper mantle, and its composition gradually changes due to the melting of the fusible material by tidal heating. Figure 5 shows how both modes of the model evolve to steady state, presenting results for eruption rate, temperature, porosity, and composition. We assume an initially homogeneous body with a bulk composition of 50% fusible material that is initially on its solidus throughout. Other initial conditions, for example starting uniformly cold, or with a cold lithosphere, result in the same broad behaviour, but starting on the solidus removes the spin-up time required to heat the mantle. Thus, despite not knowing the precise 'initial condition', various distinctive behaviours can be found that may have im-portant implications for the evolution of Io and other heat-pipe bodies. The left column of figure 5 shows the evolution of mode 1, and the right column shows the evolution of mode 2. Note that steady state is reached much more rapidly in mode 1 and so the time axis of mode 2 is significantly expanded. The final steady states are those shown previously in figure 3.

Time-Evolution to Steady State
The early (t ≤ 5 Myr) evolution of the model is the same for both modes. Fusible (pure-A) melts are produced throughout the mantle and rise upward. They are erupted onto the surface and so a cold fusible crust begins to grow. The upper mantle is being continually resupplied with fusible material as it is buried though the crust and remelted at its base. There is no such resupply of fusible material to the deep mantle, which becomes increasingly refractory. After ∼ 5 Myr, about 20% of Io's volume has been erupted and reburied; the lower mantle is almost completely depleted in fusible material. As a result, melting rate there drops and the solid starts to climb the solidus toward T = T B (figure 2). Panels a and d in figure 5 show that the decreased melting rate in the lower mantle reduces the eruption rate to almost zero. This reduction in eruption rate causes the crust to thin, increasing conductive heat loss from the surface. Once the lower mantle has been heated to T B , the 3-layer mantle structure described above in the steadystate solution emerges. From this point in the evolution onward, the mid-mantle is acting as an impermeable barrier to refractory melts formed in the lower mantle. The presence of this barrier causes melt to accumulate at the top of the lower mantle, as shown by the bright region at ∼ 1300 km in figure 5c. The accumulation of melt at the top of the lower mantle increases liquid overpressure, which initiates the extraction of refractory melt to the magmatic plumbing system. It is at this point, after around 15 Myr, that the evolution of the two modes diverge.
In mode 1, the emplacement of the refractory melts into the upper mantle creates a band of intermediate composition there, but the top of the upper mantle and the crust remain purely composed of the fusible material. Steady state is reached after ∼ 30 Myr, coinciding with the attainment of thermal equilibrium, where heat loss from eruptions equals that input by tidal heating. In mode 2, the deep refractory melts make it to the surface, and the crust -initially composed of purely fusible material -becomes of intermediate composition. As there is little to no emplacement in the upper mantle the downwelling crust maintains its composition, which results in cyclic behaviour where the composition of new crust depends on the downwelling composition of the crust a few Myr previously. For example, the initial, purely fusible crust creates a pulse of fusible melt at ∼ 40 Myr, which produces a new pulse of erupta, more fusible than that in the intervening period. This cycle continues with a decreasing amplitude of differences between erupta compositions until eventually a steady state is reached after around ∼ 200 Myr. Thermal equilibrium is reached after ∼ 100 Myr, which can be seen by the constant eruption rate after ∼ 100 Myr in figure 5e. bulk composition is fusible enough) then it is all emplaced before reaching the surface and the erupted composition is purely fusible (mode 1). For a given value of h M there is a critical bulk composition that divides mode 2 from mode 1 (figure 6d). Equivalently, -14-manuscript submitted to JGR: Planets for a given bulk composition there is a critical h M which divides mode 2 (low h M ) from mode 1 (high h M ).

Bulk Composition and Mantle Emplacement Rate
A prominent feature of figure 6 is that the crustal thickness and eruption rate both decrease at more refractory bulk compositions. When hot, refractory melt reaches the surface, the eruption rate and crustal thickness drop. The drop in crustal thickness is due to increased emplacement and the higher temperature of the material that is emplaced . The implications of this are discussed below.

Discussion
Our results demonstrate that magmatic segregation and volcanic eruptions lead to a rapid differentiation of the mantle. Fusible material is cycled in the upper mantle and crust, and its depletion at depth generates a refractory lower mantle that rises to its melting point. The fate of high temperature refractory magmas formed in the lower mantle controls the degree of chemical differentiation and the composition and temperature of erupted products. If high-temperature refractory melts freeze in the upper mantle (mode 1), no refractory lavas will be observed at the surface and the mantle will be fully differentiated. Alternatively, if refractory melts can migrate to the surface (mode 2), refractory eruptions will be observed and the mantle will not be fully differentiated.
We first discuss the differentiation caused by magmatic segregation and the mantle structure it produces. Next we discuss the key results from each mode, analysing their successes and shortcomings in explaining present observations, and their predictions for future observations. We then consider how lower mantle extraction and the migration of deep refractory melts could be interpreted physically, before finally discussing the limitations and future directions of this work.

Differentiation by Magmatic Segregation
The formation of a pure-refractory lower mantle at steady state is a necessary consequence of magmatic segregation in our model. Magmas that form in the lower mantle rise toward the upper mantle, leaving behind an increasingly refractory residuum, a feature shown in the time evolution plots in figure 5. The composition of the lower mantle only reaches steady state when all fusible material has been removed. Compositional stratification in our model can be best understood by noting that solids are continually moving downward (see solid flux in figure 3c,g), and are continually heated as they downwell. Continued heating of intermediate compositions produces fusible melts that segregate buoyantly upward, leading to increasingly refractory compositions with depth.
The structure of the mid-and upper mantle depends on both the phase diagram and the fate of refractory magmas produced in the lower mantle. For our simple twocomponent phase diagram, the upper mantle is at the fusible melting temperature T A , and the mid-mantle must span the temperature range between T A and the temperature T B of the pure-refractory region below. The reduced model, formulated in Appendix C, shows that the thickness of the mid-mantle (T A < T < T B ) is determined by the rate at which downwelling solid is heated from T A to T B , which is slowest (and thus the midmantle is thickest) when no emplacement takes place there. If emplacement of the lower refractory magma there is very efficient (see the largest value of h M in figure 6b) the midmantle is thin and there is almost complete segregation between a pure refractory lower mantle and a pure fusible upper mantle. On the other hand, if refractory melts migrate far into the upper mantle, differentiation is less complete. The upward migration reduces the thickness of the pure refractory lower mantle, and increases the thickness of intermediatecomposition upper mantle.
With a more detailed phase diagram, we would expect a general structure similar to that proposed here but with greater complexity. In particular, the chemistry of the crust and uppermost mantle would likely be much more complex, with layering controlled by melting temperature, and potentially influenced by near-surface sulphur cycling. Sul-phur may be acting as a volatile that reduces melting temperatures (Battaglia et al., 2014). The formation of a lowermost olivine layer is expected to be a feature of any relevent silicate phase diagram, and so our prediction of the formation of high temperature refractory melts is expected to hold. Any temperature range in the mantle over which there is not significant melting would be present as a low-melt-fraction layer that acts as a barrier to melts rising from below, potentially leading to magma overpressure and, in the context of our model, transfer to a plumbing system.

Implications of the Two Magmatic Modes
In this section we discuss the specific results and implications of each mode, analysing the degree to which each mode can explain current observations, and the predictions they make of future observations.
In mode 1, high-temperature refractory magmas formed in the lower mantle migrate into the upper mantle and freeze, delivering their latent heat to the fusible surroundings. The additional melting this emplacement causes can manifest as a high-melt-fraction decompacting layer, as seen in figure 3b. Magnetic induction models have been used to infer the presence of a ≥ 50 km region of ≥ 20% melt fraction beneath Io's crust (Khurana et al., 2011). This has been previously interpreted as a region of concentrated tidal heating (Tobie et al., 2005;Bierson & Nimmo, 2016), or as a decompacting boundary layer . Mode 1 of our model shows another manifestation of this decompaction hypothesis; a high melt fraction layer can arise due to freezing of deep refractory melts into the upper mantle. This is a result of the viscous resistance of the mantle to decompaction, and does not occur if the viscosity of the mantle is small, as shown by the solutions to the reduced model in figure 3, in which this viscous resistance is effectively ignored. A decompacting layer, whether caused by freezing or the strength of the crust , provides a means of generating high melt fractions in the upper mantle without requiring concentrated tidal heating in this layer.
Mode 1 predicts that no eruptions of refractory material take place. This could be considered consistent with the lack of observed olivine on the surface of Io, although this apparent absence may simply reflect an observational limitation (Geissler et al., 1999).
The key deviation of mode 1 from observations is that it does not predict any high temperature eruptions. For mode 1 to produce high temperature eruptions would require invoking processes like viscous heating on ascent (Keszthelyi et al., 2007).
In mode 2, refractory melts formed in the lower mantle rise to the base of the crust and are ultimately erupted. The relative lack of upper mantle emplacement in mode 2 means that melting throughout the upper mantle is caused predominantly by tidal heating (Moore, 2001;. A key strength of mode 2 in relation to observations lies in its prediction of high eruption temperatures. This provides a means of reconciling heat flow arguments that require heat transport by magmatic segregation (Moore, 2003;Breuer & Moore, 2015), with observations of high temperature eruptions (McEwen et al., 1998;de Kleer et al., 2014). Mode 2 supports the hypothesis of Keszthelyi and McEwen (1997) that eruptions of deep, refractory melts formed within a differentiated Io could produce very high temperature lavas. This study expands on that suggestion, demonstrating the dynamical conditions necessary for such eruptions. The rise of deep refractory melts to the surface is a means of recycling deep material to the crust, and so the upper mantle is never fully depleted in refractory material.
The eruption rate predicted by mode 2 is lower than that in mode 1. At steady state and given the negligible surface conduction, the heat lost through eruptions must equal that input by tidal heating . Increasing the temperature of erupted material means therefore that a lower eruption flux is needed (figure 6c). Despite this decreased eruption rate, in our model there is very little change in total melting. The combination of the decreased eruption rate and the approximately constant to-tal melt production means that more emplacement of intrusions takes place in bodies operating in mode 2. This effect was explained by , where it was shown that the emplaced fraction is given by C(T erupt − T s )/(L + C(T erupt − T surf )), where C is the specific heat capacity, L is the latent heat, T erupt is the eruption temperature, and T surf is the surface temperature. The increased emplacement yields a thinner crust than mode 1 for the same value of the crustal emplacement constant h C . However we note that the appropriate value of h C is not known, so larger crustal thickness could also be produced in mode 2, with emplacement spread over a larger region.
A conclusive determination of the presence or absence of olivine from Io's surface would provide a simple test for whether Io is likely to be in mode 1 or mode 2. Further, additional observations to constrain the globally averaged volcanic eruption rate and eruption temperature would also test whether refractory melts are migrating out of the deep mantle. On the basis of its ability to explain high eruption temperatures originating from a mantle governed by magmatic segregation, we propose that mode 2 is the more likely state for Io.

Mechanism of Ascent for Deep Refractory Magmas
A fundamental assumption of our model is that deep refractory melts are able to migrate out of the lower mantle without equilibration as they rise. From a modelling perspective, we assume that this occurs due to the accumulation of magmatic overpressure in the lower mantle, which enables melt to leave the lower mantle through some arbitrary 'magmatic plumbing system'. In the model, this plumbing system is treated in the same way as the plumbing system in the crust, which we envision as a network of 'heat pipes'. However, its physical manifestation in mantle may well be different. In this section we first discuss the assumption that refractory magmas can leave the lower mantle, and then discuss possible physical interpretations of the plumbing system.
If Io is indeed in a thermal steady state (Lainey et al., 2009), heat supplied to the lower mantle must be able to leave to the upper mantle. The heat being transported is primarily in the form of latent heat (Moore, 2001), which can only be lost by the freezing of lower mantle melts. If lower mantle melt was not extracted to a plumbing system, it would have to freeze at the top of the lower mantle where the temperature drops, passing its latent heat to fusible material at the base of the upper mantle, which would melt and continue heat transport upward. We consider such a perfect exchange of mass and energy unlikely due to the extreme liquid overpressures it would generate. We would expect these large liquid overpressures to cause melt to penetrate the overlying upper mantle, which is at its solidus and so is unlikely to have significant strength. Our mantle magmatic plumbing system is intended to capture the range of possible fates of this lower mantle melt. The ultimate freezing and heat transfer could take place at the very base of the upper mantle (large h M ); in a distributed region of the upper mantle (intermediate h M ); or in the crust and on the surface (small or zero h M ).
Assuming then that magma does leave the lower mantle, its rise could be accomplished in a number of ways. The lower mantle is hotter and, at the top, has a higher porosity than the overlying mid-and upper mantle. Together these create a lower bulk density that gives the potential for a Rayleigh-Taylor overturn. In our model, the entire mantle is on its solidus, so we would not expect significant resistance to such an overturn on long timescales. In this interpretation, h M parameterises the equilibration of rising refractory plumes with their surroundings. If the plumes are large and rise rapidly, the degree of equilibration may be very low, representing mode 2 of our model. Such an overturn represents a mode of convective heat transport. Another possibility is that lower mantle melts rise through a system of dikes. High magma pressure in the decompacting boundary layer may localise and nucleate fractures that are driven by magmatic buoy-ancy. It is possible that such conduits become semi-permanent features, although this would require large amounts of lateral melt transport in the decompacting boundary layer. Interpretations of our deep magmatic plumbing system as a system of dikes would presumably imply a higher value of h M than large convective plumes. Related to the concept of lower mantle melts rising through dikes is the formation of reactive channels. If rising refractory melts are corrosive to more fusible compositions, they can localise into high-flux channels (Kelemen et al., 1995;Rees Jones & Katz, 2018). Rising lower mantle melts are undersaturated in SiO 2 and so may disolve pyroxene and precipitate olivine. This could create high permeability, pure-olivine conduits that allow for the rapid upward rise of refractory melts.
We emphasise that our model makes no explicit assumption about the nature of this plumbing system, other than that it provides some mechanism for upward transport with an efficiency determined by the parameter h M . Further work might pursue a more detailed mechanistic interpretation, but that is beyond the current scope.

Model Limitations and Future Work
This work represents an initial step toward a full coupling of geodynamics and thermochemistry in heat-pipe bodies. We have used a simplified phase diagram that, whilst providing useful insight into the general processes of differentiation, could be significantly extended. Revisiting previous thermochemical modeling (Keszthelyi & McEwen, 1997;Keszthelyi et al., 2007) in light of the dynamics presented here could give a more realistic picture of the compositional structure of Io. This work also ignores the pressure dependence of melting temperature, the different latent heats of refractory and fusible material, solid-state phase changes, and the different densities of different compositions and their melts. Whilst we justified these simplifications, a more complete model would aim to incorporate their effects. In particular, Keszthelyi and McEwen (1997) noted that the 'mid-mantle' may be more dense than the refractory lower mantle, particularly if fusible melts there become Fe-rich. More detailed chemical modelling would enable such potential effects to be investigated.
In this work we have also neglected the radial distribution of tidal heating. In  it was demonstrated that the crustal balances of eruption, emplacement and crustal thickness depend only on the integrated heating from below, not its distribution. In the present case, the thicknesses and melt fractions of the different layers in the model would change with variable tidal heating with radius, but the general principles of differentiation and melt migration will hold. Future work may aim to couple dynamic models like that presented here with evolving tidal dissipation models.
Another significant simplification in our model is the assumption of spherical symmetry. Tidal heating is a function of not just radius but also latitude and longitude (Segatz et al., 1988;Ross et al., 1990), and may lead to lateral temperature differences on the order of ∼ 100 K (Steinke et al., 2020). Such considerations will be key to deciphering the links between interior dissipation and heat transport, and the surface expression of volcanism. If, as speculated above, convective overturn is a mechanism of upward migration of buoyant refractory melts, then future work should include this inherently symmetrybreaking process. The model here is developed to describe leading-order dynamics and compositional evolution; more detailed three-dimensional models are probably needed to facilitate close comparisons to specific surface observations. Such models would be best constrained by more detailed observations of eruptive heat fluxes, temperatures, and petrology.

Conclusions
In this work we have demonstrated that magmatic segregation and volcanic eruptions can rapidly lead to significant compositional differentiation of Io's mantle. This differentiation produces a refractory lower mantle and a fusible upper mantle and crust. Melting of the refractory lower mantle produces high-temperature melts that must leave the lower mantle in order to facilitate heat loss. The fate of these refractory melts controls the degree of differentiation of the mantle and the composition and temperature of erupted lavas. If high-temperature, refractory melts reach the surface, they can provide an explanation of the highest temperature observed eruption, but if they stall in the upper mantle, high temperature eruptions are not predicted. We hypothesise that Io's highest temperature eruptions originate from a deep, differentiated lower mantle, and that their eruption limits the differentiation of the upper mantle. Future observations of the petrology and temperature of eruptions will directly test this hypothesis.
The tidal heating scale ψ 0 is imposed, which gives the velocity scale q 0 which in turn gives the porosity scale φ 0 .

Appendix B Numerical implementation
Equations (A3), (A4), (A5b), (A6), (A2), and (A7) are solved for phase averaged composition c, plumbing system composition c p , compaction pressure P , enthalpy H, plumbing system flux q p , and plumbing system temperature T p respectively, using the finite volume method. Other variables are obtained from these six primary variables. In particular enthalpy and phase-averaged composition uniquely define temperature, porosity, solid composition, and liquid composition through the solidus and liquidus equations (1)-(2), the scaled definition of bulk enthalpy H = T + St φ 0 φ, and the definition of phase averaged composition c = φ 0 φc l + (1 − φ 0 φ)c s . This local (cell-wise) problem is solved with a Newton method.
For the numerical solution, we introduce a small amount of artificial diffusion of phase-averaged composition into the system as it helps to avoid discontinuous gradients in composition. The modified composition equation including this artificial diffusion is where D c is a constant that controls the size of the artificial diffusion. A value of D c ∼ 5 × 10 −4 is generally required for robust convergence, and can be decreased with grid refinement. The effect of this diffusion can be seen in figure 3d,h where the solid composition of the full model deviates slightly from that of the reduced model. Figure 3 (along with other tests not shown here) shows that the introduction of this diffusion does not affect the model results.
The monolithic system (equations (A3)-(A7)) is highly non-linear and tightly coupled (Katz et al., 2007). Robust convergence is obtained by splitting the system into three non-linear sub-system solvers shown schematically in figure B1. The first sub-system solves equation (A3) for phase averaged composition c, and equation (A6) for enthalpy H. Time integration is performed using the theta method. When θ = 0 the system is fully explicit, and is fully implicit when θ = 1. Initially θ = 0.5 is used, but if convergence fails an explicit timestep is taken. Sub-system 1 employs Newton's method (with globalization). As part of the residual evaluation for this sub-system, a local non-linear solve for porosity, temperature, and solid and liquid compositions (described above) is required.
Once a solution is found for sub-system 1, the result is passed to solver 2, which solves equation (A5b) for compaction pressure P using Newton's method (with globalization). This separates the non-linearity of permeability in equation (A5b) from the compositionenthalpy system in sub-system 1, which also computes porosity. Solver 2 also calculates the Darcy flux q and solid velocity u.
Upon convergence, the solutions to the previous two sub-systems are passed to solver 3, which contains the plumbing system equations (A4), (A2), and (A7). Placing the plumbing system equations in a separate non-linear solver separates them from the pressure dependence of extraction, and the temperature/plumbing system flux dependence of emplacement. Even so convergence can be poor when new regions of extraction emerge, which causes rapid changes to the solutions between timesteps. As per the previous two subsystems solvers, solver 3 also employs Newton's method (with globalization). If Newton fails to converge, we use a pseudo transient continuation method with implict (backward Euler) time integration. The pseudo transient problem is evolved to steady-state to yield the solutions to equations (A2), (A4), and (A7).
An adaptive time step is used. At the beginning of each time step k, a trial value for the step size ∆t k = 1.005 ∆t k−1 is selected. The time step is aborted if any of the solvers for the three sub-systems fail to converge, and the step size is reduced by 50%. In the event of multiple sub-system solve failures, when ∆t k < 1 × 10 −12 , an explicit timestep is taken using ∆t k−1 , and the process of step size reduction is repeated. The simulation is terminated if an explicit step with ∆t k < 1 × 10 −12 fails to converge.
After the convergence of all three non-linear sub-systems, a unified residual to the monolithic non-linear problem (A3)-(A7) is computed. Successive solution of the three sub-systems are continued until the 2 -norm of the residual of each discrete PDE is < 1×10 −7 . Once satisfied, the time step is accepted and the state of the time-dependent PDE is advanced in time from t k to t k+1 = t k + ∆t k . The discretisation and system of non-linear equations is solved using the Portable, Extensible, Toolkit for Scientific computation (PETSc) (Balay et al., 1997(Balay et al., , 2019(Balay et al., , 2020.

Appendix C Reduced model
Illuminating simplifications can be made to the full model by assuming small porosity and zero compaction length -this involves neglecting terms in φ 0 and δ within the scaled equations in Appendix A. Conservation of composition in the crust-mantle system becomes ∂c ∂t + 1 r 2 ∂ ∂r r 2 qc l + 1 r 2 ∂ ∂r r 2 uc s = −Ec l + M c p .
We assume that extraction E is zero outside of boundary layers at the base of any solid regions, where it acts to transfer any liquid flux q to the plumbing flux q p . E can therefore be thought of as a delta function on the boundaries between partial melt and solid (as boundary layers go to zero thickness in the zero-compaction-length approximation).
Darcy's law and the compaction relation become The reduced energy equation (A6) splits naturally into two cases: 'solid', in which case q = 0 and we have u ∂T ∂r = 1 Pe r 2 ∂ ∂r r 2 ∂T ∂r + St ψ + M (T p − T + St); (C3) Figure B1. Schematic of the solver used for the full model. The system is split into three non-linear solvers for enthalpy and composition, pressure, and the plumbing system. The solutions to each solver are iterated until all solvers agree to within some small tolerance. A pseudotransient solver is used for the pipe equations when convergence is poor.
The second term here is the melting due to the heat released when material is emplaced from the plumbing system. The final term comes from balancing energy at the interface r = r a ; since there is a temperature gradient below, the Stefan condition (jump condition for the enthalpy equation) gives a sudden melt flux where the temperature gradient here is known from the solution of (C7)-(C9). From these solutions we know the plumbing flux q p,c and liquid flux q c arriving at the crust mantle boundary r c (which is to be determined). Since the flux q c is then transferred to the plumbing system, the plumbing system in the crust subsequently has constant temperature and composition given by Note that if all refractory material has been emplaced beneath the crust, then q p,c = 0 and this simply says that the crustal plumbing system has c p = 1, and T p = T A . Within the region r c < r < R, we have to solve the system This system has the boundary conditions T = T A , ∂T ∂r = 0, q p = q c + q p,c , at r = r c , T = T s at r = R. (C17) This system is solved the same way as the mid-mantle solid region: a shooting method is used to find the position r c , as well as the crustal temperature distribution and the plumbing flux. Seeking a particular bulk composition for silicate Io, a guess can be made of r a , and a Newton method used on the resultant bulk composition to find the position of r a that gives the desired bulk composition. Figure 3 shows solutions to the reduced model as dashed lines, showing good agreement with the full model. There are slight differences in the position of the mid-mantle that arise in the full model due to the smoothed solidus (equation (1)).