Tidal controls on the lithospheric thickness and topography of Io from magmatic segregation and volcanism modelling
TTidal controls on the lithospheric thickness and topography of Iofrom magmatic segregation and volcanism modelling
Dan C Spencer, Richard F Katz, Ian J Hewitt
Highlights • Io is modelled with a spatially-variable tidal heating model coupled to a magma-segregation and vol-canism model. • We predict long-wavelength lithospheric thickness variations that arise from spatially variable tidalheating. • Lithospheric thickness could either be correlated with surface heat flux (if intrusions form at a constantrate), or be weakly anti-correlated (if intrusions form at a rate proportional to magma flux). • A Pratt-like isostasy model predicts long-wavelength topography, showing that topography is expectedto be anti-correlated with lithospheric thickness.
Keywords
Io, tidal heating, heat piping, planetary volcanism, isostasy1 a r X i v : . [ a s t r o - ph . E P ] A ug bstract Tidal heating is expected to impart significant, non-spherically-symmetric structure to Jupiter’s vol-canic moon Io. A signature of spatially variable tidal heating is generally sought in observations ofsurface heat fluxes or volcanic activity, an exploration complicated by the transient nature of volcanicevents. The thickness of the lithosphere is expected to change over much longer timescales, and so mayprovide a robust link between surface observations and the tidal heating distribution. To predict long-wavelength lithospheric thickness variations, we couple three-dimensional tidal heating calculations to asuite of column models of magmatic segregation and volcanic eruption. We find that the lithosphericthickness could either be correlated with the radially integrated heating rate, or weakly anti-correlated.Lithospheric thickness is correlated with radially integrated heating rate if magmatic intrusions format a constant rate in the lithosphere, but is weakly anti-correlated if intrusions form at a rate propor-tional to the flux through volcanic conduits. We use a simple, Pratt-like isostasy calculation to predictlong-wavelength topography and find that long-wavelength topography anti-correlates with lithosphericthickness. These results will allow future observations to critically evaluate models for Io’s lithosphericstructure, and enable their use in constraining the distribution of tidal heating.
Io, the most volcanic body in the solar system, operates in a different tectonic regime to the terrestrialplanets. The eruption and burial of lava, combined with the low surface temperature leads to the growth ofa thick and cold lithosphere in spite of the high surface heat flux. This high heat flux is primarily exportedfrom the interior by magmatic segregation in the mantle (Moore, 2001), and through volcanic ‘heat-pipes’ inthe lithosphere (O’Reilly and Davies, 1981). Heat is supplied to the interior by tidal dissipation — a processof great importance in the Solar System (de Kleer et al., 2019) — and the distribution of input tidal heatingis expected to control the surface heat flux distribution (Ross et al., 1990; Tackley, 2001; Kirchoff et al.,2011; Beuthe, 2013; Rathbun et al., 2018; Steinke et al., 2020). Whilst theoretical links between the tidalheating distribution and the surface heat flux are well established, the implications for interior structure,magma dynamics, tectonics and topography are not well known.The spatial distribution of tidal heating is a longstanding and still largely unresolved problem in planetaryscience (Segatz et al., 1988; Roberts and Nimmo, 2008; Beuthe, 2013; Bierson and Nimmo, 2016; Renaudand Henning, 2018). The end-members generally considered for Io are that of lower-mantle heating or as-thenosphere heating (Segatz et al., 1988; Ross et al., 1990; Tackley et al., 2001; Hamilton et al., 2013),though magma-ocean dissipation has also been proposed (Tyler et al., 2015; Hay et al., 2020). Lower-mantledissipation predicts high polar heat fluxes, whereas asthenospheric dissipation predicts high equatorial heatfluxes. A number of works have sought to identify tidal dissipation patterns from surface heat fluxes (Veederet al., 2012), volcanic activity (Rathbun et al., 2018), and volcano distributions (Ross et al., 1990; Kirchoffet al., 2011; Hamilton et al., 2013). The primary hindrance to these works is the poor polar coverage ofobservations, so whilst a number of these works favour an asthenosphere heating model (e.g., Ross et al.(1990); Kirchoff et al. (2011)), the general consensus is that more polar observations are needed to fully ad-dress this question (Rathbun et al., 2018; de Kleer et al., 2019). Further, long-timescale, averaged heat fluxesare difficult to estimate, and it is unclear to what extent short-timescale observations of volcanic activityreflect the global dissipation structure. Tectonic features, which vary on vary on much longer timescales,may provide a more robust link between surface observations and the distribution of tidal heating.2n important tectonic feature that is expected to relate to the surface heat flux, and thus the tidal heatingdistribution, is the long-wavelength lithospheric thickness (Ross et al., 1990; Steinke et al., 2020). Recentstudies have proposed two hypotheses for the primary controls on the thickness of the lithosphere, whichwe define as the upper-most, fully solid layer of Io. Steinke et al. (2020) proposed a stagnant lid convectionmodel where a portion of mantle heat transport occurs by convection, and this convectively-transportedheat is transported though the lithosphere by conduction. This predicts that the lithosphere is thinnestwhere heat flux is highest. Alternatively, Spencer et al. (2020) proposed that the eruption and burial of lavaresults in the growth of a cold lithosphere, with a steady-state thickness that is controlled by the balance ofdownward advection and heat delivered by magmatic intrusions. Conduction plays a minor role in this modelbecause the rate of burial is so large. In such a system, the lithospheric thickness is primarily controlledby the rate of melt production and the rate of intrusive heating. It should be noted that both of theseprevious studies referred to the surface layer as the ‘crust’; here we use ‘lithosphere’ instead because thecrust is usually considered to be a petrologically distinct layer, a distinction that becomes important in theisostatic calculations below. If the thickness of the lithosphere can be related to topography and heat flow,then long-wavelength variations in lithospheric thickness could be used to infer the tidal heating distribution.Together with estimates of libration amplitude (VanHoolst et al., 2020), this can provide a powerful meansof investigating Io’s interior structure.Each of the models described above propose different controls on the lithospheric thickness and so may beexpected to predict different relationships between the tidal heating distribution and thickness. Steinke et al.(2020) used radially integrated tidal heating profiles to predict the effect of spatially variable tidal heating onlithospheric thickness, finding that the thickness anti-correlates with surface heat flux. In the present workwe extend the model of Spencer et al. (2020) to consider the effect of variable tidal heating on the eruptionand intrusion model for lithospheric thickness such that comparisons can be made between the models ofSpencer et al. (2020) and Steinke et al. (2020).We generalise the simplified, steady-state model of Spencer et al. (2020) to allow variable tidal heating. Iois divided into a set of laterally contiguous columns that are coupled to a viscoelastic tidal heating model.The tidal heating model calculates a three-dimensional heating rate from a spherically symmetric rheologicalstructure. This leads to a recognised limitation of these type of tidal heating models; the three-dimensionalheating rate that they produce generates a non-spherically symmetric structure that cannot be used torecalculate the heating distribution without averaging over spherical shells (Roberts and Nimmo, 2008;Bierson and Nimmo, 2016). Thus, models coupling such tidal heating calculations to dynamics cannot be fullythree-dimensional. We use this coupled, pseudo-three-dimensional model to investigate the links betweentidal heating, lithospheric thicknesses, and long-timescale eruption rates/heat fluxes. Our results show thatthe relationship between lithospheric thicknesses and heat flux depends on how magmatic intrusions formwithin the lithosphere. If the rate of formation of permanent magmatic intrusions is independent of the(non-zero) magma flux through volcanic conduits, as may be expected if the volcanic system exploits pre-existing fractures, we predict the lithosphere to be thickest where radially integrated heating rate (and thuseruption rate and heat flux) is highest. If, however, magmatic intrusions form at a rate proportional to themagma flux in volcanic conduits, as may be expected if volcanic conduits form due to basal magma pressurethat generates new pathways for magma to propagate into the lithosphere, the lithospheric thickness shouldbe weakly anti-correlated with radially integrated heating rate. We also use a simple, Pratt-like isostasycalculation to convert dissipation-derived long-wavelength lithospheric thickness to topography, predicting3hat topography anti-correlates with thickness. This relates a feature that generally has to be indirectlyinferred (lithospheric thickness), to an observation that is more readily obtained (topography). Improvedobservations of surface heat fluxes and their relationships to lithospheric thickness and topography will testdifferent models for the controls on Io’s lithospheric thickness (Steinke et al., 2020; Spencer et al., 2020).With a means of critically evaluating these models, the structure of tidal heating can feasibly be constrainedby future estimates of lithospheric thickness.
Our model consists of three parts: a theory for magmatic segregation and volcanism, another for tidaldissipation, and isostasy calculations. The one-dimensional magmatic segregation and volcanism model is ageneralisation of the asymptotic approximation in Spencer et al. (2020). In it, melting is driven by the tidaldissipation model, which most closely follows the approach of Beuthe (2013), utilising a Maxwell viscoelasticrhology. Rheological parameters required by the tidal calculation are predicted by the segregation andvolcanism model, completing the coupling of the two systems. The isostasy calculations utilise the equal-pressure formulation of Hemingway and Matsuyama (2017).The dynamics are described by the magmatic segregation and volcanism model. Spencer et al. (2020) derivea system where tidal heating causes the formation of magma in the mantle that rises buoyantly toward thesolid lithosphere (termed crust in that work). High magma overpressure just below the base of the lithospherefacilitates a transfer of magma from the pore space into a lithospheric magmatic plumbing system, which canbe thought of as a system of dikes. Magma rising in this plumbing system can freeze into the cold, surroundinglithosphere, forming permanent magmatic intrusions, delivering both mass and energy to the surroundings.The rest of the magma in the plumbing system rises to the surface and erupts, imparting a compensatingdownward flux of the (now cold) erupted products. Spencer et al. (2020) found that the delivery of heatfrom the freezing of magmatic intrusions is required to raise the temperature of cold, downwelling lithospheresuch that a lithospheric thickness within observational constraints can be maintained. This concept of theemplacement of permanent magmatic intrusions is an important one in the present work and is discussedbelow.The dynamic model is coupled to tidal dissipation to yield a consistent, three-dimensional structure. Weuse a spherically symmetric structure to calculate a three-dimensional heating rate. This heating rate distri-bution (which importantly is not spherically symmetric) is applied to a suite of column models, producinga three-dimensional structure. This structure is averaged over spherical shells and used to re-calculatethe heating distribution. This processes is iterated until the heating-distribution converges, yielding thethree-dimensional structures presented in this work. We utilise a Maxwell viscoelastic law despite the well-documented inability of such a rheological law to produce observed dissipation rates at realistic mantleviscosities (Bierson and Nimmo, 2016; Renaud and Henning, 2018). We also neglect all lateral flow, justifiedby the long-wavelength of the tidal forcing; the one-dimensional columns are considered isolated. This isa significant simplification that we discuss below, and we note that future work should aim to analyse thepropensity for lateral flow. We also inherit some of the assumption of Spencer et al. (2020), namely thatwe ignore the chemical composition and as a consequence neglect the possibility of compositional convectionin the mantle. Parameter values, where not explicitly stated, are given in table 1 of Spencer et al. (2020),4hough the shear viscosity used in the tidal heating calculation is different to that used in the column models,as explained below.
The model of Spencer et al. (2020) is based on conservation equations for mass, momentum, and energy ina compacting two-phase medium, together with conservation of mass in a magmatic plumbing system thattransports magma through the solid lithosphere. Here, we make use of the simplified model described inappendix B of that paper. In the mantle, which is at the melting temperature T m , tidal heating producesmelt, and mass conservation of the melt phase reads1 r ∂∂r (cid:0) r q (cid:1) = ψρL (1)where q = K φ n ∆ ρg/η l is the Darcy segregation flux related to the porosity φ , where ψ is the local vol-umetric heating rate (see section 2.2) and L is the latent heat. Here K is a permeability constant, n isthe permeability exponent, ∆ ρ is the density difference between the solid and liquid, and η l is the magmaviscosity (numerical values for these and other parameters are as in table 1 of Spencer et al. (2020)). Themagma flux is therefore q = 1 ρL r (cid:90) rr m ψr d r, (2)where r m is the base of the mantle. At the base of the lithosphere this flux is transferred to the plumbingsystem, in which the flux is denoted q p . Conservation of mass and energy in the lithosphere are describedby 1 r ∂∂r (cid:0) r ( u + q p ) (cid:1) = 0 , (3a)1 r ∂∂r (cid:0) r q p (cid:1) = − M, (3b)and 1 r ∂∂r (cid:0) r uT (cid:1) = 1 r ∂∂r (cid:18) r κ ∂T∂r (cid:19) + ψρC + M (cid:18) T m + LC (cid:19) , (4)where u is the solid velocity, T is the temperature, M is the emplacement rate (the rate at which magmaticintrusions remove material from the plumbing system), and C is the specific heat capacity. The final term inequation (4) represents the heating that emplacement provides to the downwelling lithosphere. The solutionof equations (3)–(4) together determines the temperature profile in the lithosphere as well as the lithosphericthickness (see Spencer et al. (2020) for details).In Spencer et al. (2020), we assumed a temperature-dependent parametrisation of the emplacement rate M ,but this assumption is problematic in the present case. With a coupled calculation of the tidal heating rate ψ ,there is very little tidal heating of the cold lithosphere, and hence more heat is required from emplacementto limit the growth of the lithosphere. Using the temperature-dependent form for emplacement, we findthat there is too little heating of the lithosphere to avoid it becoming unreasonably thick. We thereforeconsider two alternative parametrisations of the emplacement rate. First, that magmatic emplacement isat a constant rate, and second, that magmatic emplacement is a function of the amount of material in the5lumbing system. To allow for both possibilities we take the emplacement rate to be M = λ c + λ q q p , (5)and explore cases where only one of λ c or λ q is non-zero at a time. Taking a constant emplacement rate( λ c (cid:54) = 0 and λ q = 0) can be interpreted as modelling a system of dikes where the number of dikes is fixed butthe flux through them varies. If emplacement is a function of contact area with the host rock, such a systemcould result in emplacement rate being independent of the magma flux. This is similar to, but more simplethan the temperature dependence taken in Spencer et al. (2020). Taking emplacement to be proportionalto the amount of melt in the plumbing system ( λ c = 0 and λ q (cid:54) = 0) can also be interpreted as a systemof dikes, but where the dikes have equal fluxes and the number of dikes varies. As the flux (and thus thenumber of dikes) increases, the contact area with the host rock also increases, and so the total emplacementrate increases. In summary, we consider cases where emplacement is positively related to, or independentof the magma flux. We do not consider the possibility of a negative relationship between emplacement rateand magma flux as we cannot conceive of a realistic physical system that this would represent. For the calculation of tidal heating we most closely follow the methodology of Beuthe (2013). Volumetrictidal dissipation averaged over an orbit is given by (Tobie et al., 2005) ψ ( r, θ, ϕ ) = ω f σ ij )Re(˜ (cid:15) ij ) − Re(˜ σ ij )Im(˜ (cid:15) ij )] , (6)where ω f = 4 . × − s − is the orbital frequency, ˜ σ ij and ˜ (cid:15) ij are the components of the complex stressand strain tensors, and summation over components i and j is implied. We calculate the complex stressand strain tensors using the propagator matrix approach detailed in Sabadini and Vermeersen (2004), andexplained in appendix A of Roberts and Nimmo (2008).The tidal potential that forces the system arises from consideration of a synchronous eccentric orbit, to firstorder in eccentricity. It is given by (Kaula, 1964; Tobie et al., 2005)Ω = r ω f e (cid:20) − P (cos θ )cos( ω f t ) + 14 P (cos θ ) [3cos( ω f t )cos(2 ϕ ) + 4sin( ω f t )sin(2 ϕ )] (cid:21) , (7)where e = 4 . × − is the orbital eccentricity, θ and ϕ are the colatitude and longitude (the latter beingzero at the sub-Jovian point), t is the time, and P and P are associated Legrendre polynomials.To couple the tidal heating model to the dynamical model we follow the approach of Bierson and Nimmo(2016). We take the shear viscosity to be a function of temperature and porosity through the relationship(Katz, 2010) η = η exp (cid:20) E A R g (cid:18) T − T (cid:19) − Λ φ (cid:21) , (8)where E A = 3 × J/mol is the activation energy, R g is the gas constant, η is a reference viscosity at thereference temperature T (taken to be the melting point), and Λ = 27 is a positive constant. Temperature T and porosity φ are extracted from the model in section 2.1 and averaged over spherical shells, so η dependsonly on radius r . The value of η used is chosen so that the total global rate of tidal dissipation approximately6atches the observed dissipation rate of ∼ × W (Lainey et al., 2009). It is well documented that aMaxwell viscoelastic constitutive law requires a very low viscosity to produce the amount of tidal heatingobserved in Io (Segatz et al., 1988; Tackley, 2001; Bierson and Nimmo, 2016; Steinke et al., 2020). We assumethat this is a failure in the present understanding of the rheology that affects dissipative processes (Biersonand Nimmo, 2016; Renaud and Henning, 2018), rather than a reasonable assesment of Io’s long-timescalemantle viscosity. Bierson and Nimmo (2016) also take a porosity dependence of the elastic shear modulus,but we neglect this small effect in line with our simplified approach. We refer to the first coupled model,using (8), as the ‘mantle heating’ model.Numerous previous works have considered the possibility that tidal dissipation is concentrated within alower-viscosity asthenosphere (e.g., Segatz et al. (1988); Tackley (2001); Hamilton et al. (2013); Davies et al.(2015)). Such a dissipative layer does not arise in the above formulation, even when a large decompactingboundary layer is included in the dynamic model (Spencer et al., 2020), because the porosity dependence in aMaxwell viscoelastic model is too weak. In order to investigate the lithospheric thickness and long-wavelengthtopography implications of such a dissipation structure, we calculate an alternative ‘asthenospheric heating’model, where the shear viscosity in a 300 km layer beneath the lithosphere is set to be a factor of 1000 lowerthan the rest of the mantle. In the asthenospheric heating model we do not include the temperature andporosity dependence of shear viscosity; in this case the heating model is decoupled from the dynamical model.We do, however, set the shear viscosity in the cold lithosphere is be effectively infinite so no dissipation occursthere, consistent with the calculated dissipation structure in the coupled mantle heating model.The tidal heating code has been benchmarked against the radial functions in figure 2 of Tobie et al. (2005),against the TiRADE software used in Roberts and Nimmo (2008), and by reproducing figures 8 and 10 ofSegatz et al. (1988).
For our isostasy calculations we follow Hemingway and Matsuyama (2017) in using an equal-pressure for-mulation of isostasy in spherical coordinates. This assumes that compensated columns have equal pressuresat their bases (the compensation depth, r ccd ). Equal pressure isostasy assumes that we have (Hemingwayand Matsuyama, 2017) P = (cid:90) Rr ccd ρg d r, (9)where P is a constant (independent of latitude and longitude), ρ is the density profile, and R is local planetaryradius. We take gravity g to be uniform for simplicity, a reasonable assumption given the likely heavy core.We assume that density is a function of temperature only ρ = ρ [1 − α ( T − T m )] , (10)where α = 3 × − K − is the coefficient of thermal expansion, and ρ is the reference density at themelting temperature T m . It is at this point that the distinction between the crust and lithosphere becomesimportant for Io. In terrestrial systems, the base of the crust represents a petrological boundary between thelow-density crust and the high-density mantle. Here, however, the boundary is simply an isotherm. Effectiverecycling of lithospheric material into the partially molten mantle is expected to remove any significant7ompositional variation across this boundary (Spencer et al., In Review). As such, we expect that the coldlithosphere is more dense than the underlying hot, partially molten mantle. This inverted density structureis assumed to be stable because of the high viscosity of the rigid upper lithosphere. We also neglect densitydifferences associated with the presence of low-density melt in the upper mantle; the density of the mantle is ρ m = ρ . By considering density variations within the lithosphere, this calculation is closer to Pratt isostasythan Airy isostasy, though we are clear that this is not true Pratt isostasy as the base of the lithosphere isn’tconsidered to be flat (Fowler, 2004).The integral in equation (9) can be split at the base of the lithosphere, which has a thickness l to write P = ρ g ( R − l − r ccd ) + (cid:90) l ρg d z, (11)where z = R − l is the distance downward from the surface. Both l and ρ (in terms of temperature T )are known from the magmatic segregation and volcanism model, so this expression can be re-arranged todetermine the variable radius R relative to its spatial average R . Since P and r ccd are constant, we obtainthis topography h = R − R as h = l + (cid:90) l ρρ d z + constant , (12)where the constant is chosen to make the spatial average of h zero. Figure 1 shows model solutions for the lithospheric temperature distribution, mantle porosity, and tidalheating distribution at Io’s north pole and three points around the equator, for the (coupled) mantle heatingmodel and the (de-coupled) asthenosphere heating model. In the mantle-heating case (figure 1a–c), heatingrate is highest at the poles, and lowest at the sub- and anti-Jovian points, whereas in the asthenosphere-heating case (figure 1d–f), heating rate is highest at the sub- and anti-Jovian points, and lowest at the poles.A higher heating rate leads to increased melt production, though for the permeabilities used here meltfractions only vary by ∼ λ c (cid:54) = 0 and λ q = 0), integration of equation (3b) in the lithosphereyields q p = q e R r + λ c (cid:18) R r − r (cid:19) , R − l ≤ r ≤ R, (13)where q e = q p ( r = R ) is the eruption rate. Assuming negligible surface conduction, the eruption rate mustbe given by a column-wise energy balance as (Spencer et al., 2020) q e = 1 R ( ρL + ρC ( T m − T s )) (cid:90) Rr m ψr d r, (14)where the integral is the total tidal heating delivered to the column. From equation (2), the plumbing flux8t the base of the lithosphere is q p ( r = R − l ) = 1( R − l ) ρL (cid:90) R − lr m ψr d r. (15)Since negligible tidal heating takes place in the lithosphere (figure 1), the integrals in (14) and (15) areessentially identical. Thus, equating (15) with (13) at the base of the lithosphere yields an analyticalexpression for the lithospheric thickness in terms of the local eruption rate l = R − R (cid:18) − q e Rλ c C ( T m − T s ) L (cid:19) / . (16)A Taylor expansion of the term in brackets provides some intuition into this expression. Expanding to thefirst non-trivial term yields l ≈ C ( T m − T s ) L q e λ c . (17)The thickness of the lithosphere is controlled by the balance between latent heat release in the lithosphereand sensible heat loss at the surface. The greater the temperature difference between erupting lava andthe surface, the more heat that must be provided to downwelling material to raise it to its melting point.As the eruption rate increases, material downwells more quickly, and with no corresponding increase inemplacement rate, the thickness of the lithosphere grows. This effect can be seen in the main panels of figure1. A higher rate of emplacement means that downwelling material is heated more rapidly, reducing thelithospheric thickness. We note that an average lithospheric thickness can be estimated using the modelledglobal average eruption rate of Spencer et al. (2020).The insets in panels a and d of figure 1 show the lithospheric temperature profiles when emplacement rate isproportional to the plumbing system flux ( λ c = 0 and λ q (cid:54) = 0). In this case equation (3b) can be integratedto give q p = R q e r e λ q ( R − r ) . (18)Again assuming negligible surface conduction and equating equation (18) to the total melt production in theinterior (equation (15)) gives an expression for the lithospheric thickness, l = 1 λ q ln (cid:18) C ( T m − T s ) L (cid:19) . (19)Interestingly, this is independent of the melting rate, so lithospheric thickness is expected to be virtuallyconstant when emplacement rate is proportional to the plumbing system flux. A Taylor expansion of (19) tofirst order yields equation (17) but with q e /λ c replaced by 1 /λ q , illustrating that in this case, the relationshipbetween eruption flux and emplacement is fixed. The small variations in lithospheric thickness seen in theinsets in panels a and d of figure 1 are due to conduction (which is neglected in arriving at the estimate,equation (19)), with higher heating rates producing thinner lithospheres.Figure 2 shows lithospheric thickness, eruption rate, and topography as a function of latitude and longitudein the coupled mantle-heating model. The top row of figure 2 shows the case where emplacement rate isa constant and the bottom row shows the case where emplacement rate is proportional to the plumbingsystem flux. A constant emplacement rate means that lithospheric thickness correlates with the eruptionrate, as specified by equation (16). Lithospheric thickness varies by about 25 km, with the most pronounced9 igure 1: Lithospheric temperature profiles, mantle porosities, and tidal heating distributions at the poles andthree points around the equator, for the mantle heating (a–c) and asthenosphere heating (d–f) models, with constantemplacement rate, λ c = 1 .
66 Myr − . Panels a and d show the temperature variation in the upper-most mantle andlithosphere, indicated by the green region in the other panels. Where radially integrated heating rate is highest, meltproduction and porosity is highest. This results in an increased eruptive flux and the growth of a thicker lithosphere.The insets in panels a and d show the case when emplacement rate is proportional to plumbing system flux, with λ q = 0 .
05 km − . In this case, lithospheric thicknesses vary weakly (porosity and tidal heating profiles in the mantleare almost exactly the same as constant emplacement rate). Dots indicate estimates of lithospheric thickness usingequations (17) (panels a and d) and (19) (insets). Differences between the analytical estimates and the model arecaused by conduction, which is neglected in the analytical estimates. variation being between the thick polar lithosphere and the thin equatorial lithosphere. Our isostatic modelassumes that density is only a function of temperature, and so the cold lithosphere must be more densethan the underlying, partially molten mantle. This results in topographic highs where the lithosphere isthinnest. The coupled, mantle-heating model with constant emplacement rate predicts long-wavelengthtopography with an amplitude of about 250 m. In the case where emplacement rate is proportional to theamount of material in the plumbing system, shown in the bottom row of figure 2, the lithospheric thicknessonly varies by a couple of kilometres and the amplitude of long-wavelength topography is <
40 m. Thiscan be understood through equation (19); increased heating and the resultant increased eruption rate isbalanced by increased emplacement, resulting in an almost uniform lithospheric thickness. In this case, thelong-wavelength lithospheric thickness and topography variations are a result of different conductive heat10uxes and so lithospheric thickness is anti-correlated with eruption rate (Ross et al., 1990; Steinke et al.,2020).
Figure 2:
Solutions for lithospheric thickness, eruption rate, and topography in the case of coupled dynamics andtidal heating. Tidal heating is concentrated in the lower mantle in the coupled model, producing maximum eruptionrates at the poles (see figure 1). Panels a–c show the case where emplacement rate is constant, and panels d–f showthe case where emplacement rate is proportional to the plumbing system flux. Constant emplacement rate predictsa correlation of lithospheric thickness with eruption rate (or heat loss), and topographic lows where heat flux is high.An emplacement rate proportional to plumbing system flux predicts a relatively uniform lithospheric thickness andlittle long-wavelength topography.
Figure 3 shows the same plots as figure 2, but for the case of asthenospheric heating. All of the relationshipsbetween heating rate, eruption rate, lithospheric thickness, and topography are the same in this case, but thepattern of dissipation and so the pattern of the plotted solutions is different. Asthenospheric heating predictshigher eruption rates at the equator. If emplacement rate is constant, this predicts a thicker lithosphere atthe equator (amplitude ∼
30 km), and topographic highs at the poles (amplitude ∼
300 m). If emplacementrate is proportional to the amount of material in the plumbing system, lithospheric thickness is much moreuniform (amplitude ∼ ∼
90 m), with lithospheric thicknessvariations being controlled by variation in conductive heat fluxes.Assuming dominantly vertical flow — a significant assumption that we discuss below — the global patternof heat flow should be reflective of the tidal heating distribution, as has been noted elsewhere (e.g., Segatzet al. (1988); Tackley (2001); Veeder et al. (2012); Davies et al. (2015)). The primary means to distinguishbetween lower mantle and asthenospheric heating models is on the basis of heat flux. Lower mantle heatingpredicts higher polar heat flux, whereas asthenosphere heating predicts higher equatorial heat flux. Withthe present dearth of polar observations, this is a difficult distinction to make. Rigorous observation of Io’spoles is required to understand which mode of heating is more likely to be occurring. However, if the mode11 igure 3:
Solutions for lithospheric thickness, eruption rate, and topography in the case of asthenosphere heating.Panels a–c show the case where emplacement rate is constant, and panels d–f show the case where emplacementrate is proportional to the plumbing system flux. Relationships between lithospheric thickness, eruption rate, andtopography are the same as in figure 2, but patterns and amplitudes are different due to the different heating mode. of emplacement can be established, lithospheric thickness and topography can serve as a useful proxy forlong-timescale heat flux.This work predicts that the long-wavelength variations in lithospheric thickness should either correlate withthe long-timescale eruption rate/heat flux, or be weakly anti-correlated, as summarised schematically infigure 4. In the constant emplacement rate model, we predict that lithospheric thickness correlates witheruption rate. An explanation for why emplacement would be independent of magma flux is that volcanicconduits are not formed by magma pressure at depth, but rather tectonic processes in the lithosphere. Io’seruption and burial tectonics are thought to form mountains by thrust faulting (McKinnon et al., 2001;Kirchoff and McKinnon, 2009). If, for example, such faults can act as conduits for magma ascent, freezing ofascending magma on their walls may be largely independent of the flux through the conduit. Alternatively, inthe flux-proportional emplacement rate model, we predict that long-wavelength lithospheric thickness variesby only a few kilometers, and is weakly anti-correlated with heat flux. A rationale for why emplacement ratewould be proportional to volcanic plumbing flux may be that volcanic conduits are created by overpressuredmagma at the base of the lithosphere. It is plausible that higher melt production in the interior wouldlead to a larger number of conduits. If magma in each of these conduits has a chance of stalling withinthe lithosphere, this would imply a positive relationship between lithospheric magma flux and emplacementrate.The flux-proportional emplacement rate model makes predictions for variations in lithospheric thickness thatare similar to the results of Steinke et al. (2020). When comparing this work to Steinke et al. (2020), it is12 igure 4:
Schematic illustrating the primary results of this work. The lithospheric thickness is correlated withradially integrated heating rate if magmatic intrusions form at a constant rate (panel a), but is approximatelyuniform (or weakly anti-correlated) if intrusions form at a rate proportional to the flux through volcanic conduits(panel b). important to note that whilst both can predict a conductive control on lithospheric thickness variations, thecontrols on the absolute values of lithospheric thickness are different. In this work the lithospheric thicknessis primarily controlled by the rate of magmatic emplacement, whereas in Steinke et al. (2020) the lithosphericthickness is controlled entirely by conduction through a stagnant lid. To address the relative importance ofconvective heat transport in the mantle likely requires a model that couples two-phase flow and convection,a significant challenge due to the different timescales on which these processes operate.The proposed link between lithospheric thickness and topography provides a means of relating more readily-obtainable observations (topography) to the predictions of lithospheric thickness in works like this oneand Steinke et al. (2020). Topography can be compared to eruption rates and volcanic heat fluxes toclarify the heat-transfer and emplacement mechanisms in the lithosphere. Alongside with recent work thatdemonstrates a way to constrain interior structure from libration amplitudes (VanHoolst et al., 2020), thisprovides a powerful means to investigate Io’s interior structure and heating distribution. We predict long-wavelength topographic highs where the lithosphere is thin. However, our isostatic calculations assume thatdensity variations are due only to temperature differences, and that there is no compositional boundary atthe base of the lithosphere. It is likely that the compositional profile in the lithosphere is complex, reflectingshallow magma fractionation, sulfur cycling, and other processes. If the vertical structure of the lithosphereis approximately uniform with latitude and longitude, and simply scaled to lithospheric thickness, the resultsof this work should be largely unchanged. If, however, there is significant variation in lithosphere compositionwith latitude and longitude, the applicability of the isostatic model presented here would be reduced. It isnot clear, however, that any such variation would mirror the degree-two tidal forcing, and so may averageout on the long wavelengths considered here.White et al. (2014) created a partial stereo-topographic DEM of Io that found a system of longitudinally ar-ranged alternating basins and swells near the equator, with amplitudes ∼ ∼
400 km.13his amplitude is significantly greater than the largest topography predicted here (order hundreds of me-ters). Potential explanations for this discrepancy are that there are lateral compositional differences thatresult in significant additional isostastically compensated topography, or that dynamic topography causedby upwelling mantle plumes is significant (Tackley et al., 2001). It is important to note however that thereare considerable discrepancies between stereo-derived and limb-profile-derived long-wavelength topography(White et al., 2014), and hence that long-wavelength topography is not well constrained. Improved ob-servations of long-wavelength topography, particularly in the polar regions, are required to make robustcomparisons between modelled topography and data.A primary limitation of this work is the neglect of lateral flow in either the lithosphere or mantle. Differencesin lithospheric thickness are expected to be counteracted by deformation of the lithosphere. Such calculationsare common in studies of the ice-shells of icy satellites (Stevenson, 2000; Nimmo and Stevenson, 2001; Nimmo,2004), where there is a clear rheological and density transition at the base of the shell. The application ofsuch a model to Io is not straightforward because rheological and density transitions are expected to be moregradual (Spencer et al., In Review). It is not clear whether there is an easily defined petrological ‘crust’ ofIo. Nonetheless such lateral flow is possible, and would be best investigated by a two-dimensional model ofupper Io. Lateral flow is also possible in the partially molten mantle. Pressure gradients would be expectedto drive flow of the mobile magma phase. Pressure gradients could be produced by processes such as differentmelting rates or spatially variable extraction rates to the lithosphere. An investigation of lateral melt flowwould likely require a two-dimensional model of the partially molten mantle. Here we simply note that therelationships proposed in this model are expected to hold if vertical motion is much greater than lateralmotion, as generally expected in Io’s eruption and burial tectonics at long wavelengths.
We have demonstrated how spatially variable tidal heating leads to long-wavelength variations in lithosphericthickness in a model of magmatic segregation and volcanic eruptions. Our models predict that such variationsare controlled by how magma intrudes into the lithosphere. If permanent magmatic intrusions form at arate independent of the magma flux through volcanic conduits, the lithosphere should be thickest wheretidal heating is greatest. In this case the lithopshere thickness can vary by 10s of km. If however magmaticintrusions form at a rate proportional to the magma flux through volcanic conduits, lithospheric thicknesswill only vary by a few km, and will be anti-correlated with eruption rates. We also predict that if densitydifferences are predominantly derived from temperature differences, then areas of thin lithosphere will sit ontopographic highs. Improved observational constraints on eruption rates, heat fluxes, and long-wavelengthtopography, particularly at Io’s poles, will help distinguish between different models for the controls onlithospheric thickness.
Acknowledgements
This work was funded by the
Science and Technologies Facilities Council , the University of Oxford’s
OxfordRadcliffe Scholarship , and
University College, Oxford . This research received funding from the EuropeanResearch Council under the European Unions Horizon 2020 research and innovation programme grant agree-14ent number 772255. We thank F. Nimmo and P. England for their comments and suggestions. Source codeand data used in the production of figures can be found at https://zenodo.org/record/3975657 (Spencer,2020).
References
Beuthe, M. (2013). Spatial patterns of tidal heating.
Icarus , 223(1):308–329.Bierson, C. J. and Nimmo, F. (2016). A test for Io’s magma ocean: Modeling tidal dissipation with apartially molten mantle.
Journal of Geophysical Research: Planets , 121(11):2211–2224.Davies, A. G., Veeder, G. J., Matson, D. L., and Johnson, T. V. (2015). Map of Io’s volcanic heat flow.
Icarus , 262:67–78.de Kleer, K., McEwen, A. S., and Park, R. (2019). Tidal Heating: Lessons from Io and the Jovian System.In
Final Report for the Keck Institute for Space Studies .Fowler, C. M. R. (2004).
The Solid Earth: An Introduction to Global Geophysics . Cambridge UniversityPress, Cambridge, 2 edition.Hamilton, C. W., Beggan, C. D., Still, S., Beuthe, M., Lopes, R. M. C., Williams, D. A., Radebaugh, J.,and Wright, W. (2013). Spatial distribution of volcanoes on Io: Implications for tidal heating and magmaascent.
Earth and Planetary Science Letters , 361:272–286.Hay, H. C. F. C., Trinh, A., and Matsuyama, I. (2020). Powering the Galilean Satel-lites with Moon-moon Tides.
Geophysical Research Letters , n/a(n/a):e2020GL088317. eprint:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020GL088317.Hemingway, D. J. and Matsuyama, I. (2017). Isostatic equilibrium in spherical coordinates and implica-tions for crustal thickness on the Moon, Mars, Enceladus, and elsewhere.
Geophysical Research Letters ,44(15):7695–7705. eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2017GL073334.Katz, R. F. (2010). Porosity-driven convection and asymmetry beneath mid-ocean ridges.
Geochemistry,Geophysics, Geosystems , 11(11).Kaula, W. M. (1964). Tidal dissipation by solid friction and the re-sulting orbital evolution.
Reviews of Geophysics , 2(4):661–685. eprint:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/RG002i004p00661.Kirchoff, M. R. and McKinnon, W. B. (2009). Formation of mountains on Io: Variable volcanism and thermalstresses.
Icarus , 201(2):598–614.Kirchoff, M. R., McKinnon, W. B., and Schenk, P. M. (2011). Global distribution of volcanic centers andmountains on Io: Control by asthenospheric heating and implications for mountain formation.
Earth andPlanetary Science Letters , 301(1):22–30.Lainey, V., Arlot, J.-E., Karatekin, \ ., and Van Hoolst, T. (2009). Strong tidal dissipation in Io and Jupiterfrom astrometric observations. Nature , 459(7249):957–959.15cKinnon, W. B., Schenk, P. M., and Dombard, A. J. (2001). Chaos on Io: A model for formation ofmountain blocks by crustal heating, melting, and tilting.
Geology , 29(2):103–106.Moore, W. B. (2001). The Thermal State of Io.
Icarus , 154(2):548–550.Nimmo, F. (2004). Non-Newtonian topographic relaxation on Europa.
Icarus , 168(1):205–208.Nimmo, F. and Stevenson, D. J. (2001). Estimates of Martian crustal thickness from viscous re-laxation of topography.
Journal of Geophysical Research: Planets , 106(E3):5085–5098. eprint:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2000JE001331.O’Reilly, T. C. and Davies, G. F. (1981). Magma transport of heat on Io: A mechanism allowing a thicklithosphere.
Geophysical Research Letters , 8(4):313–316.Rathbun, J. A., Lopes, R. M. C., and Spencer, J. R. (2018). The Global Distribution of Active IonianVolcanoes and Implications for Tidal Heating Models.
The Astronomical Journal , 156(5):207. Publisher:American Astronomical Society.Renaud, J. P. and Henning, W. G. (2018). Increased Tidal Dissipation Using Advanced Rheological Models:Implications for Io and Tidally Active Exoplanets.
The Astrophysical Journal , 857(2):98.Roberts, J. H. and Nimmo, F. (2008). Tidal heating and the long-term stability of a subsurface ocean onEnceladus.
Icarus , 194(2):675–689.Ross, M. N., Schubert, G., Spohn, T., and Gaskell, R. W. (1990). Internal structure of Io and the globaldistribution of its topography.
Icarus , 85(2):309–325.Sabadini, R. and Vermeersen, B. (2004).
Global Dynamics of the Earth: Applications of Normal ModelRelaxation Theory to Solid-Earth Geophysics . Kluwer Academic Publishers, Dordrecht.Segatz, M., Spohn, T., Ross, M. N., and Schubert, G. (1988). Tidal dissipation, surface heat flow, and figureof viscoelastic models of Io.
Icarus , 75(2):187–206.Spencer, D. C. (2020). Spencer-space/Io tidal: IoTidal v1.1.0.Spencer, D. C., Katz, R. F., and Hewitt, I. J. (2020). Magmatic Intrusions Control Io’s Crustal Thickness.
Journal of Geophysical Research: Planets , 125(6):e2020JE006443.Spencer, D. C., Katz, R. F., Hewitt, I. J., May, D. A., and Keszthelyi, L. P. (In Review). Compositionallayering in Io driven by magmatic segregation and volcanism.Steinke, T., Hu, H., Hning, D., van der Wal, W., and Vermeersen, B. (2020). Tidally induced lateralvariations of Io’s interior.
Icarus , 335:113299.Stevenson, D. J. (2000). Limits on the variation of thickness of Europa’s ice shell. In , volume 1506.Tackley, P. J. (2001). Convection in Io’s asthenosphere: Redistribution of nonuniform tidal heating by meanflows.
Journal of Geophysical Research: Planets , 106(E12):32971–32981.Tackley, P. J., Schubert, G., Glatzmaier, G. A., Schenk, P., Ratcliff, J. T., and Matas, J.-P. (2001). Three-Dimensional Simulations of Mantle Convection in Io.
Icarus , 149(1):79–93.16obie, G., Mocquet, A., and Sotin, C. (2005). Tidal dissipation within large icy satellites: Applications toEuropa and Titan.
Icarus , 177(2):534–549.Tyler, R. H., Henning, W. G., and Hamilton, C. W. (2015). Tidal Heating in a Magma Ocean withinJupiter’s Moon Io.
The Astrophysical Journal Supplement Series , 218(2):22.VanHoolst, T., Baland, R.-M., Trinh, A., Yseboodt, M., and Nimmo, F. (2020). The Librations, Tides,and Interior Structure of Io.
Journal of Geophysical Research: Planets , 125(8):e2020JE006473. eprint:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020JE006473.Veeder, G. J., Davies, A. G., Matson, D. L., Johnson, T. V., Williams, D. A., and Radebaugh, J. (2012). Io:Volcanic thermal sources and global heat flow.
Icarus , 219(2):701–722.White, O. L., Schenk, P. M., Nimmo, F., and Hoogenboom, T. (2014). A new stereo topographic mapof Io: Implications for geology from global to local scales.