UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Hydrologic modelling with variable fluid properties : a liquid sulfur aquifer beneath the surface of… Wijns, Christopher P. 1995

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1996-0110.pdf [ 3.44MB ]
JSON: 831-1.0052884.json
JSON-LD: 831-1.0052884-ld.json
RDF/XML (Pretty): 831-1.0052884-rdf.xml
RDF/JSON: 831-1.0052884-rdf.json
Turtle: 831-1.0052884-turtle.txt
N-Triples: 831-1.0052884-rdf-ntriples.txt
Original Record: 831-1.0052884-source.json
Full Text

Full Text

H Y D R O L O G I C M O D E L L I N G W I T H V A R I A B L E F L U I D P R O P E R T I E S : A LIQUID S U L F U R A Q U I F E R B E N E A T H T H E S U R F A C E O F IO By Christopher P. Wijns B. Sc. (Geophysics) McGill University A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E STUDIES G E O P H Y S I C S A N D A S T R O N O M Y We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A December 1995 © Christopher P. Wijns, 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Geophysics and Astronomy The University of British Columbia 129-2219 Main Mall Vancouver, Canada V6T 1Z4 Date: Tan.SI : 1 9 % Abstract The return of the first Voyager images of Io in 1979 provoked a controversy over the importance of sulfur versus silicate volcanism in resurfacing the satellite. The debate involves both the strength of crustal material as well as the mode of mass and heat transport to the surface of Io. The dynamics of sulfur magma transport will depend upon the properties of liquid sulfur. Unlike the Newtonian behaviour of silicate magmas, molten sulfur has anomalies in the heat capacity (a A-like transistion at 432 K) and in the viscosity (a maximum at 460 K) due to a polymerization reaction. Numerical modelling of a pure liquid sulfur supply region at depth in a porous crust shows that the fluid supply rate from such a magma reservoir is about seven orders of magnitude less than the observed mass and heat fluxes at the surface of Io. This requires accumulation of magma in a reservoir prior to eruption at the surface, but the accumulation time on the order of tens of thousands of years is much longer than would be consistent with the level of volcanic activity observed on Io. Low sulfur magma supply rates suggest that sulfur is not the dominant resurfacing material for the satellite. This is consistent with calculations performed for the viscous relaxation of topography on Io. Surface features of a largely sulfur crust would disappear within months because of material weakness. Significant topographic highs and steep slopes on Io provide evidence for a crust with stronger mechanical properties. Sulfur convection in the crust can still provide a large fraction of the global heat flow, within certain conditions on the permeability. Previous estimates of resurfacing rates have equated the cooling of erupted material to the entire global energy budget. Significant convective transfer of heat by liquid sulfur may call for a downward revision of resurfacing rates. u Table of Contents Abstract ii List of Tables vi List of Figures vii Acknowledgement ix Preface 1 1 Sulfur Volcanism on Io 2 1.1 Introduction 2 1.2 Observations and Constraints 4 1.3 Model 10 1.3.1 Initial and boundary conditions 10 1.3.2 Fluid thermodynamic equations 12 1.3.3 Parameters 14 1.3.4 Numerical model 18 1.4 Model Results 20 1.5 Discussion 23 1.5.1 Thermal stratification of the aquifer 23 1.5.2 Mass and heat fluxes 26 1.5.3 Crustal recycling 30 1.6 Conclusions 30 iii 2 Model Construction and Benchmarking 32 2.1 Detailed Mathematical Model 32 2.1.1 Basis Functions 33 2.1.2 Finite Element Equations 35 2.2 Validation Tests 38 2.2.1 One-Dimensional Thermal Conduction 38 2.2.2 One-Dimensional Thermal Advection 40 2.2.3 Two-Dimensional Thermal Convection 42 2.3 Summary 43 3 Numerical Experiments with Variable Fluid Properties 45 3.1 Heat Capacity 47 3.2 Viscosity 50 3.3 Summary and Application to Io 52 4 Viscous Relaxation of Topography on Io 55 4.1 Introduction 55 4.2 Model 56 4.2.1 Solution to the viscous flow equations 56 4.2.2 Boundary conditions 60 4.3 Model Parameters 65 4.4 Results and Discussion 67 Overview and Future Directions 72 References 74 Appendices 78 iv A Liquid Sulfur Thermodynamic Data and Equations of State 78 B Boussinesq Approximation 82 C A Note on the Derivation of the Heat Conservation Equation 84 v List of Tables 1.1 Sulfur Aquifer Parameters 16 3.1 Effect of Variable Heat Capacity on Fluid Velocity 48 4.1 Viscous Relaxation Parameters 66 A . l Density 79 A.2 Specific Heat Capacity at Constant Pressure 80 A.3 Dynamic Viscosity 81 vi List of Figures 1.1 Geotherm and Sulfur Phases 9 1.2 Crustal Structure and Aquifer Model 11 1.3 Liquid Sulfur Properties 15 1.4 Model A: Convection Patterns in Aquifer 21 1.5 Model A: Convected Heat Flux 22 1.6 Influence of Permeability on Convected Heat Flux 23 1.7 Model B: High-Permeability "Pipe" 24 1.8 Model C: Magma Withdrawal . 25 1.9 Power Radiated from Pele 27 2.1 Finite Element Grid 33 2.2 One Dimensional Thermal Conduction 39 2.3 One Dimensional Thermal Advection 41 2.4 Critical Rayleigh Number 43 3.1 Constant Property Fluid at Racr 46 3.2 Effect of Variable Heat Capacity on Fluid Velocities 49 3.3 Effect of Variable Heat Capacity Away From Racr 50 3.4 Effect of Variable Viscosity on Fluid Velocities 51 3.5 Effect of Variable Viscosity Away From RaCT 52 4.1 Mathematical Domain 56 4.2 Relaxation Flow Field 67 vii 4.3 Dependence of Relaxation Time on Viscosity 69 4.4 Dependence of Relaxation Time on Depth to Basement 70 4.5 Dependence of Relaxation Time on Wavelength 71 A . l Fit to Density Data 79 A.2 Fit to Heat Capacity Data 80 A.3 Fit to Kinematic Viscosity Data 81 C l Effect of Erroneous Heat Equation on Thermal Convection 88 vm Acknowledgement Thanks foremost to Susan Kieffer for all kinds of guidance and encouragement and pa-tience. She was always ready to devote full attention to any of my results, both feasible and not. Bruce Buffett and Roger Beckie were kind enough to review my thesis in the rush before Christmas. Thanks to Bruce for pointing out an error in my original heat conservation equation, which I have discussed in appendix 3. ix Preface The advent of powerful computers and workstations has opened up the study of non-linear fluid dynamics. Numerical simulations can be valuable aids in understanding sometimes complex flow patterns that arise due to variations in fluid properties and domain bound-aries. Models of mass and heat transport by fluids in porous media are increasingly able to take into account the many parameters involved in characterizing both the fluid and the porous matrix. The finite element model developed in the following chapters has been tested against problems for which analytic solutions exist, and subsequently used to investigate the effects of variable fluid viscosity and heat capacity on thermal convection in a porous medium with spatially heterogeneous permeability and porosity. The numerical code has been applied to the problem of sulfur volcanism on Io as a quantification of a hypothesis introduced by McEwen and Soderblom (1983). Chapter 1 tests their suggestion, within observed and calculated constraints for Io, that sulfur may feed volcanic eruptions from a liquid zone at depth in a porous silicate crust. The second chapter develops in detail the mathematical solution of the conservation equations by the finite element method and presents the benchmarking against known analytic solutions. The isolated effects of a temperature dependence in the heat capacity and viscosity of a fluid are treated in chapter 3. The final chapter focusses again on Io. In considering the merits of sulfur versus silicate volcanism for resurfacing the satellite, constraints imposed by surface topography are examined in chapter 4. An analytic treatment of the viscous flow equations for topographic relaxation provides insight into the likelyhood of an upper crust on Io which is rich in sulfur. 1 Chapter 1 Sulfur Volcanism on Io 1.1 Int roduct ion Active volcanism on Io was predicted by Peale et al. (1979) shortly before the Voyager spacecraft rendezvous with this satellite of Jupiter. They suggested that tidal heating would be dissipated by widespread volcanic activity. This prediction was borne out the same year by the two Voyager observations of volcanic plumes up to 300 km high. A debate emerged and continues today over the nature of the erupting material. The controversy between silicate and sulfur volcanism is based upon surface features of Io, infrared temperature estimates at hot spots, and spectroscopic evidence. Either side must address the nature of a magma source and the transport of fluid to the surface. The premise of a sulfur supply region lends itself to models of crustal structure and magma generation which differ markedly from those for silicate volcanism because the fluid dynamics and transport of a such a magma will depend upon the properties of liquid sulfur, which are quite different from those of molten silicate. Perhaps the most important fluid property that may allow us to distinguish sulfur from silicate transport is the variation of viscosity with temperature. As sulfur is heated through its liquid temperature range (393 K to 717 K at 1 atmosphere), it undergoes a viscosity increase of four orders of magnitude at 432 K, from less than 0.01 Pa-s to almost 100 Pa-s, due to a polymerization reaction. The viscosity decreases with temperature from 460 K onwards as the polymer chains are shortened by increasing thermal energy 2 Chapter 1. Sulfur Volcanism on Io 3 (Touro and Wiewiorowski,1966). The viscosity variation of silicate liquid is much simpler. A typical viscosity for a basaltic magma is about 1 Pa-s (Turcotte and Schubert, 1982), less than the sulfur maximum. Although a silicate melting temperature about 900 K higher than the melting point for sulfur will dictate a deeper source region for a silicate magma, the extreme temperature dependence of the viscosity of sulfur may exert greater control on a sulfur magma flow regime. McEwen and Soderblom (1983) considered the implications of this viscosity variation for volcanism on Io and speculated on the subsurface dynamics of a liquid sulfur supply region. They proposed that a reservoir of liquid sulfur within a silicate crustal matrix, at temperatures spanning the sulfur viscosity maximum, would contain two zones of more mobile fluid separated by the high-viscosity zone. They hypothesized that these two low-viscosity zones formed a sulfur aquifer associated with two different types of volcanic eruptions differing in magma temperature and composition - low temperature (?a400 K) S-SO2 events, and 650 K or higher eruptions of sulfur. McEwen and Soderblom (1983) related the colour of the surface material around a volcanic feature to the eruption temperature of the magma (obtained from infrared heat flow measurements). Elemental sulfur undergoes distinct colour changes at different temperatures in the liquid range. Sulfur dioxide and liquid sulfur quenched from below 432 K are white and pale yellow under laboratory conditions. These colours were recorded around low temperature calderas on Io. At 600 K and higher, sulfur is dark brown and black. Such dark deposits were observed around high-temperature calderas, and McEwen and Soderblom (1983) associated these deposits with hotter sulfur. If the magma ascent is isothermal, as suggested by Kieffer (1982), then the temperature of the magma source region can be inferred directly from the colour of surface deposits. McEwen and Soderblom (1983) used this argument to propose that the low-viscosity zone at the top of a liquid sulfur reservoir heats SO2 and drives the low-temperature eruptions which have Chapter 1. Sulfur Volcanism on Io 4 light-coloured deposits, and that the low-viscosity zone at greater depth supplies hotter sulfur magma to the high-temperature eruptions and is responsible for the dark surface deposits. The qualitative model of McEwen and Soderblom (1983) is developed in this thesis into a numerical model of the liquid sulfur aquifer. This is used to examine the feasibility of such a reservoir as a source region for sulfur magma on Io, within the context of crustal structure and the observational constraints of mass and heat flux at the surface. 1.2 Observations and Constraints Mass and heat flux estimates at the surface of Io provide some constraints on volcanic activity. Infrared spectroscopy has provided values for the global average heat flow and for temperatures and heat fluxes at localized hot spots coinciding with the locations of volcanic calderas (e.g. Pearl and Sinton,1982 ; McEwen et al.,1985 ; Veeder et al.,1994). Tidal forces are assumed to produce the bulk of the escaping heat energy, which has been estimated at up to 1014 W (Veeder et al.,1994), for a global average of 2.5 W/m 2 . Reseachers agree that most of this heat flow is radiated from the hot spots (Veeder et al.,1994) and is therefore not representative of the conductive heat flux through the lithos-phere. Reasonable values for this conductive flux are as low as 0.1 W/m 2 (Carr,1986) or as high as 0.5 W/m 2 (McEwen et al.,1989). Associated near-surface temperature gradi-ents range anywhere from 25 K/km to 160 K/km for a silicate lithosphere, depending on the thermal conductivity. A number of researchers have suggested that the upper crust of Io is rich in sulfur. This had already been proposed before the Voyager encounters (e.g Nash and Fanale, 1977). After the return of the first photographic images of the surface, sulfur was quickly championed based on the match between the colours of Io and those of various allotropes Chapter 1. Sulfur Volcanism on Io 5 of elemental sulfur (e.g. Smith et al., 1979). This colour-matching scheme has since been called into question by Young (1984) on the basis of the instability of those coloured allotropes. He argues that the conditions on the surface of Io should cause the reversion of red S 3 and S 4 to the stable Ss form of sulfur. A combination of low temperatures and bombardment of the exposed surface by x-rays would prevent S 3 and S 4 from existing for more than a few hours. Since Sg is colourless below 150 K, and the average surface temperature on Io is 130 K (Nash et al.,1986), sulfur should not be responsible for any colouration. Spectroscopic data remains, however, a source of inspiration for those who argue the case of abundant sulfur on the surface of Io. The similarity of Io's spectral reflectance in the UV through near infrared wavelength range with laboratory spectra of powdered sulfur (e.g. Fanale et al.,1974) led to crustal models such as that proposed by Smith et al. (1979). Their model consisted of layers of sulfur dioxide and sulfur near the surface, with possibly a molten sulfur "ocean" below this and overlying a silicate basement. However, topographical features on Io with relief of 2 km preclude a dominantly sulfur crust. Clow and Carr (1980) argue that structural failure would occur for the observed caldera depths and scarp heights. Chapter 4 details independent calculations for viscous relaxation of sulfur topography on Io which confirm this conclusion. Brief relaxation times of less than one week argue against the likelihood of a crust composed mainly of sulfur. Modifications of the Smith et al. (1979) layered crustal model must therefore be sought. McEwen and Soderblom (1983) propose a crust that is largely non-sulfur for a different reason. If sulfur exists as the molten ocean described by Smith et al. (1979), the Rayleigh number, although ill-defined for a fluid with such variable properties, will be very large (> 1015), and therefore the fluid will be freely convecting and thermally mixed. However, the two classes of volcanic eruptions require the existence of two distinct low-viscosity zones within the liquid sulfur layer. Invoking flow in a porous medium reduces the Chapter 1. Sulfur Volcanism on Io 6 Rayleigh number of the system by many orders of magnitude and may enable the sulfur layer to remain thermally stratified. For constant properties, the ratio of the Rayleigh numbers for non-porous and porous flow in a layer of depth b is Ra b2 where k is the permeability of the porous medium. The depth of a porous layer containing liquid sulfur will be greater than the depth of an entirely fluid sulfur layer because the thermal gradient through a porous silicate lithosphere will be lower. However, with the parameters adopted herein, the Rayleigh number is still more than thirteen orders of magnitude less, whereas the critical Rayleigh number for convection in a porous medium is little more than one order of magnitude less than that for convection in a purely fluid layer. This should prevent thorough mixing of the sulfur layer and ensure that the two low-viscosity temperature ranges remain at separate depths. For this reason McEwen and Soderblom (1983) propose a model where liquid sulfur is mobile within a matrix of silicates or impure sulfur compounds. Arguments for sulfur eruptions have been widely based on Voyager infrared data, which returned estimates of surface temperatures at volcanic calderas from 400 K to 650 K (Pearl and Sinton, 1982), spanning the liquid sulfur range. Early ground-based measurements similarly pointed to temperatures slightly greater than 600 K (Sinton, 1981). The type example of what McEwen and Soderblom (1983) classified as a high temperature sulfur eruption is associated with the volcano Pele. Pearl and Sinton (1982) used infrared data to model Pele as a 654 K circular area of radius 6 km superimposed on a larger, low temperature hot spot. McEwen and Soderblom (1983) match this result, and the dark surface deposits around Pele, to the extrusion and quenching of 650 K molten sulfur. More recent observations and calculations point to higher temperature infrared outbursts >1000 K. Blaney et al. (1994) modelled a 1990 outburst on Io as Chapter 1. Sulfur Volcanism on Io 7 being a large, erupting lava flow cooling from 1225 K. Two analyses of an earlier event in 1986 conferred temperatures of 900 K and 1550 K on that outburst. Such temperatures are consistent with the presence of silicate melts (Veeder et al.,1994), but are too high to account for liquid sulfur at surface pressures. Average mass flux estimates for eruptions can be obtained from resurfacing rates. An order of magnitude method for determining a resurfacing rate through energetics equates the heat energy lost by cooling erupted magmas to the global heat flux A i £ = 1 0 1 4 W (e.g. Blaney et al.,1994). This neglects the contribution of the kinetic energy of plumes, which Johnson and Soderblom (1982) have estimated at less than one percent of the total. Taking into account cooling through both liquid and solid phases, AE = PmcmATm + PsLs + p3c3ATa (1.1) where the subscripts m and s refer to the molten and solid states, respectively, of the density p, heat capacity c, temperature T, and latent heat of fusion (solidification) L. Using values from Turcotte and Schubert (1982) for silicate properties, with an initial melt temperature of 1300 K and a final temperature of 130 K, equal to the surface temperature of Io, a silicate resurfacing rate of 2.2 x l O 4 m3/s (1.66 cm/year) will provide the required energy budget. This is similar to the value 1.33 cm/year obtained by Blaney et al. (1994). Sulfur has a much lower heat of fusion ( « 4 x 104 J/kg), and so 1.4xl05 m3/s, or 10 cm/year, must be extruded to supply the same power upon cooling from 650 K. Veeder et al. (1994) report intense infrared emissions during 4% of their observation time. If this is taken as a measure of eruptive frequency, then during any one event the eruption rate is 25 times the average resurfacing rate. The volume flux for a sulfur eruption would be 3.5 x 106 m3/s, which is equivalent to about 7.2 x 109 kg/s. This translates into an advected heat flux of 4.7 x 1015 W during the course of a 650 K Chapter 1. Sulfur Volcanism on Io 8 eruption. The same analysis for resurfacing by 400 K sulfur yields mass and heat fluxes per eruption of 1.1 x IO10 kg/s and 5 x 1015 W, respectively. The silicate resurfacing rate only demands an effusion rate one tenth of that for hot sulfur, which again agrees with a calculation by Blaney et al. (1994) based on their observation of a diffusing infrared outburst. The geotherm will constrain the depth at which sulfur melts. Carr (1986) calculated various geothermal profiles based on the maximum heat flux that could be sustained through the lithosphere without causing thinning. A model geotherm from Carr (1986) for a 30 km thick lithosphere with a thermal conductivity of 4.5 W/m is redrawn in figure 1.1, along with the sulfur liquidus and vapourization curves. Also included are the steepest geotherm advocated by Smith et al. (1979) and a conservative geotherm for a pure sulfur upper crust based on a conductive heat flux of 0.1 W/m 2 . Although the melting temperature is relatively unaffected by depth, the boiling point is a strong function of pressure. The vapourization curve and the geotherm never intersect, even for the Smith et al. (1979) gradient. Sulfur will remain stable as a liquid deep into the lithosphere, and therefore it is warranted to model a single phase fluid in the aquifer. Carr's (1986) geotherm is adopted as a constraint on temperatures and associated fluid properties throughout the aquifer, providing conditions with which to test the hypothesis by McEwen and Soderblom (1983) that the two low-viscosity zones contribute to separate styles of volcanic eruptions. The derived estimates of mass and heat fluxes during volcanic events are compared with those of the theoretical modelling. Chapter 1. Sulfur Volcanism on Io 9 Figure 1.1: Geotherm and Sulfur Phases 0 10 =£15 20 25 30 0 i r i 1 "^""—v.-~ ~ ~" ~ f ^ ^ ^ ^ ^ C . p . ^ 15 MPa ^ ••v " — <£U Mra e bu Mra a \ \ b \ i i i i i i II 200 400 600 800 1000 1200 1400 T(K) (a) Model geotherm redrawn from Carr (1986). (b) Sulfur liquidus (from Vezzoli et al.,1969). (c) Vapourization curve of sulfur terminat-ing in critical point c.p. (from West,1950, c.p. from Tuller,1954) (d) Top of sulfur aquifer (depth of melting), (e) Upper limit on geotherm suggested by Smith et al.(1979). (f) Pure sulfur geotherm. Chapter 1. Sulfur Volcanism on Io 10 1.3 Model 1.3.1 Initial and boundary conditions A hydrological approach is adopted to model the magma source regions for sulfur volca-noes on Io. The region of the lithosphere below the melting depth of sulfur, illustrated in figure 1.2, constitutes the field of interest. Since liquid sulfur may exist deep into the crust, there is no well-defined lower boundary, and free inflow-outflow conditions are assumed at the bottom. The domain is modelled as a two-dimensional, single-phase aquifer of pure liquid sulfur. Saturated conditions are assumed, and the porous crust is non-deforming, although it may be heterogeneous in permeability and porosity. The assumption of pure sulfur is justified by an interest in the role of fluid thermodynamic properties in controlling flow in the layer. Attention is confined to the horizontally ex-tended region of liquid sulfur, and phenomena associated with a possible magma chamber or fluid conduit to the surface, which may be very complex and involve phase changes, are not dealt with. There are no internal heat sources within the layer, and convection is driven by thermal buoyancy forces alone. Three basic cases (A, B, and C) are examined. All models are run in a box of depth 10 km and width 6 km. The top of the aquifer remains at a constant temperature of 393 K, and the bottom at 1000 K. There is no mass or heat flux normal to the axis of symmetry at x=0. A superhydrostatic pressure of zero at the base of the domain permits fluid inflow or outflow as needed, in order to provide the least restrictive boundary condition for the ill-constrained lower region of the aquifer. Model A is designed to establish the convection patterns in the aquifer. A sinusoidal temperature perturbation is imposed and allowed to evolve in time. There is no flow through the vertical sides, in order to simulate a convective, pattern repeating infinitely in the x direction. Chapter 1. Sulfur Volca.nism on Io 11 130 K 393 K-Figure 1.2: Crustal Structure and Aquifer Model planetary surface ' v y x , \ \ ' v y y y v \ \ v v X A X x x pbrdUS • • ' &itfurVag'ma' chamber silicate liquid sulfur region / :1000 ^ " • ^ ^ > O ^ T ^ k e r n e l . <x.x outflow T=393 K | E = 0 5x i i 8x =0 P=Po 8y |E=0 8y 0 ) 9-CL 15 ctt £ C L 1c Y p=0 8x or p=0 8x or T=T0(y) T=1000 K A saturated aquifer of pure liquid sulfur exists within a permeable sil-icate crust. Partially molten silicate intrusions may supply a source of heating from below. The numerical model does not include a possible sulfur magma chamber and conduit to the surface. Some models in-clude a region of greater permeability representing a fracture network leading to the magma chamber. Boundary conditions are shown in the enlarged illustration of the mathematical domain. Chapter 1. Sulfur Volcanism on Io 12 As a test of the influence of the crustal permeability, a variation of model A is run at permeabilities of 2 x 10 - 1 1 m 2 to 8 x 10 - 1 1 m 2 and the associated average heat flow through the aquifer is calculated. Initial conditions for this suite of models are identical to those for model A with the exception that the bottom boundary is 1.5 km below the top of the aquifer, at a temperature of 484 K. As well, the domain is only 1.5 km long, corresponding to the most stable wavelength of perturbation. Model B has the same initial and boundary conditions as model A, but includes a vertically extended, higher-permeability zone along the reflection axis x=0. This het-erogeneity extends out to x=100 m, and is illustrated in the enlarged crustal section in figure 1.2. The permeability is increased by an order of magnitude inside the "pipe" to represent possible greater fracturing below, and leading to, the magma chamber. Model C allows fluid to escape from the top of the aquifer to feed a magma chamber or conduit. The initial conditions are the steady-state results of model A. The zero mass flux condition is removed for two outflow nodes feeding the magma chamber from the top boundary, from x=0 to £=100 m. Fluid is allowed to enter from the boundary at x=Q km to replenish the aquifer. This is accomplished by specifying the superhydrostatic pressure to be zero at all of these nodes. There is effectively an infinite supply of liquid sulfur from which to recharge the aquifer. The rise and escape of magma is driven solely by fluid density contrasts and not via a pressure drop. 1.3.2 Fluid thermodynamic equations The governing equations are those of mass, momentum, and energy balance. The Boussi-nesq approximation for an incompressible fluid in a non-deforming porous matrix leaves only the velocity v in the mass conservation equation. PVi,i = 0 (1.2) Chapter 1. Sulfur Volcanism on Io 13 The momentum equation is given by Darcy's Law k pvi = (pti ~{p- po)9i) (1.3) where p is the superhydrostatic pressure only. The velocity is pre-multiplied by the density in order to use the kinematic viscosity v as a variable, rather than the dynamic viscosity. This renders the equations more tractable because the density does not appear explicitly in each term. The density is also implicitly allowed to vary with temperature by virtue of the definition of v, so that the Boussinesq approximation is not strictly followed. The permeability is denoted by k, p0 is a reference density, and g is the gravitational force, which acts in the —y direction. In the heat equation dT (1.4) c is the specific heat capacity at constant pressure, A is the thermal conductivity, T is temperature, and t is time. Unsubscripted properties are those of the fluid and the subscript m denotes material averaged properties. These averages are weighted according to the porosity <j> as Xm = <f>Xfluid + (l — (f>) Xsolid Allowing for spatial variations in permeability and porosity, equation (1.4) is expanded in cartesian coordinates as dT puc + (A - As) dx dT dx + pvc+ (A - As) dy dT dy (1.5) "A m -r - A m „ _ — 0 dx* dy' Here X3 indicates the thermal conductivity of the solid porous medium, and u and v are the x and y components of velocity, respectively. Equation (1.6) is solved simultaneously Chapter 1. Sulfur Volcanism on Io 14 with the following flow equation obtained by the substitution of equation (1.3) into the continuity equation (1.2): (kdv_dT_ dk\ dp (kfodT dk\ dp k^p \v dT dx dx) dx \v dT dy dy j dy dx2 dy2 \ dp (dT dTo\ (dk kdv_dT_ + [ dT { dy ~ dy ) + ( P " po)[dy ~ v dT dy Fluid properties vary with temperature only. The physical constraints demanded by the Boussinesq approximation (including the neglect of pressure effects) are outlined in Appendix B. 1.3.3 Parameters The temperature dependence of the density, heat capacity at constant pressure, and kinematic viscosity of liquid sulfur are presented graphically in figure 1.3 for the liquid temperature range at atmospheric pressure. A comparison with the original data is illus-trated in Appendix A. The variation of density with temperature is constant everywhere, and the heat capacity remains constant at 1020 J/kg/K after 717 K. The viscosity in-creases linearly, beyond 717 K, from about 7 x 10~5 Pa-s to twice that value at 1000 K. Numerical experiments in chapter 3 reveal that the small variation in viscosity from 393 K to 432 K is not likely to affect fluid convection, and so the viscosity is approximated as constant in that temperature region. Fluid and crustal properties are listed in table 1.1 along with other physical parameters that are identical for all models presented below. Assumptions about fluid properties and boundary conditions far below the viscosity maximum are difficult to make. The depolymerization temperature for sulfur does not change appreciably with pressure (Vezzoli et al.,1969), and so the viscosity is likely to remain on the same order as it is at 717 K and atmospheric pressure. Experimental data from Heath (1995) at 1 GPa and 723 K shows a viscosity for sulfur of about 2 Chapter 1. Sulfur Volcanism on Io 15 Figure 1.3: Liquid Sulfur Properties 1 8 0 0 1 7 0 0 1 6 0 0 1 6 0 0 * 1 4 0 0 3 1 2 0 0 CL o 1 0 0 0 L 0 . 0 6 r -ST 0 . 0 4 -CM > 0 . 0 2 -01 M > 1 0 1 0 1 0 1 0 -6 4 0 0 (a) Density (Tuller,1954). (b) Specific heat capacity at constant pres-sure (West, 1959). (c) Kinematic viscosity (based on dynamic viscosity in Tuller,1954). (d) Log-linear plot of kinematic viscosity. The dotted line is the 432 K polymerization temperature. Chapter 1. Sulfur Volcanism on Io Table 1.1: Sulfur Aquifer Parameters gravitational acceleration g 1.7 m/s2 top boundary temperature 393 K silicate matrix thermal conductivity A 4.5 W/m/K density p 2700 kg/m3 specific heat capacity cp 1000 J/kg/K permeability k 6 X 1 0 _ u m 2 porosity (j> 0.10 liquid sulfur ~A 0.2 W/m/K p (kg/m3) 2032 - 0.59r cp (J/kg/K) T < 432 K 220 .1e( T - 4 2 5 ) / 8 + 1000 T > 432 K 5 0 8 e - ( r - 4 3 2 ) / 5 9 + 1020 kinematic viscosity v (m2/s) T < 432 K 5 x IO" 6 T > 432 K [ ( ^ p ) s e c h ( ^ ) ] 2 + ^ + 5 x l 0 -Chapter 1. Sulfur Volcanism on Io 17 Pa-s. Although this pressure is more than ten times the pressures encountered in the lithosphere of Io, the results lend value to the assumption of a relatively constant viscosity above 717 K. If the viscosity remains constant, the heat capacity is less likely to show extreme variation until critical point temperatures are reached. Thus the heat capacity of liquid sulfur is also assumed to remain constant below the high-viscosity barrier. The bottom of the modelling domain is kept no higher than 1000 K to avoid speculating on the fluid properties and dynamics at elevated pressures and temperatures. The thermal conductivity of the silicate crust is the same as that chosen by Carr (1986), in order to remain consistent with his geotherm. The density reflects that of an iron-poor, or silicic, crust. Smith et al. (1979) argued that if molten silicate intrusions came into contact with sulfur, iron in the silicate would tend to react with the sulfur to form sulfides. The denser sulfides would sink and permanently remove this sulfur from the supply region. Thus a differentiated, silicic crust is assumed. The natural permeability of unfractured silicates is much lower than 10 - 1 1 m 2. A typical value for granite, a silicic rock, is 10 - 1 9 m 2 (Turcotte and Schubert, 1982), which is basically impermeable. There would be almost no fluid circulation in such material. Thus fractures must play an important role in the flow. The adopted value represents a large-scale average which includes the contribution of fractures, and possibly fragmented pyroclastic silicate material. The mechanics of fracture flow are likely to better describe the fluid movement, and the use of Darcy's law for the momentum equation is a limitation of the numerical model. A permeability value becomes somewhat arbitrary, and is chosen to accomodate the upper limit on the velocity-grid cell size constraint in the model. The relationship between the porosity and the permeability depends upon the as-sumed geometry of the pores and their connections. It affects the calculations insofar as material-averaged properties are concerned. A porosity much smaller than 0.10 would Chapter 1. Sulfur Volcanism on Io 18 not allow the thermal properties of sulfur to contribute to the fluid dynamics and would change the effective Rayleigh number of the system. 1.3.4 Numer i ca l model Solutions to the pressure and temperature fields are approximated by the finite element method on a two-dimensional rectangular grid with triangular elements. Equations (1.6) and (1.6) are both of the form d£ dt, dt d2i d2i Cot, +C1-K7 + C2"S~ + C 3 7 T " + C 4 ^ - r + C 5 — + C 6 = 0 (1.7) ot ox dy ox2 dy2 where the a are constant during solution. The linear equation (1.7) is written L(0 = 0 and the finite element approach requires that the approximate solution £, when multiplied by an appropriate weighting coefficient u>i over the whole domain, satisfies JJ u>i L(i) dx dy = 0 (1.8) for i — 1,... 7i, where n is the number of grid nodes. The value of £ within elements is expressed as a linear combination of the values at the node points. i(x, y, t) = Nj(x, y)£j{t) , j = 1,... n Taking the weighting functions u>i in equation (1.8) equal to the basis functions Nj is referred to as the Galerkin formulation. These basis or interpolation functions are linear, of the form Ni = dix + biy + a Ni has a value of unity at node i and is zero at all other nodes. Substitution of the summation expression for £ into equation (1.8) results in a matrix equation for each Chapter 1. Sulfur Volcanism on Io 19 element. These element equations are solved simultaneously for all elements in the grid in the global matrix equation Tij + Sij (j + qi - 0 Tij is the coefficient matrix of transient terms, Sij contains the spatial coefficients, and qi is a vector of constant terms including known fluxes. The time derivative is dealt with by an implicit finite difference scheme, so that the above equation becomes Grid spacing is constrained according to typical expected velocities and the thermal diffusivity through a parameter Pe known as the Peclet number. Pe = < 2 K This experimentally determined criterion ensures that the distance AX between grid nodes is of the same order as the characteristic thermal diffusion distance y/ni. For a highly non-linear problem such as the one under investigation, it was found that Pe must be smaller than suggested above, in order to accomodate sudden changes in fluid properties. The Courant number Ci is a similarly derived parameter which controls the size of the time step by allowing fluid to move no further than the distance across one element (at velocity V) during one time step. « - ^ < ' Iteration within one time step is carried out between equations (1.3), (1.6) and (1.6) until convergence in the temperature field is satisfied. Boundary conditions allow the specification of a constant temperature or heat flux, and constant pressure or mass flux. Chapter 1. Sulfur Volcanism on Io 20 1.4 Model Results The nature of the convective regime of the sulfur aquifer can be established from the results of model A. Convection will develop fully only in the upper region (figure 1.4). There is no interaction between the upper and lower sections of the aquifer because they are separated by the high-viscosity zone. Velocities in the deeper low-viscosity zone are small since the viscosity there is more than an order of magnitude higher than in the upper region. The length and time scales of convection are different above and below the high-viscosity barrier. Maximum velocities near the top of the aquifer are about 3.5 mm/day, almost ten times those in the lower part. The Rayleigh number in the lower part of the aquifer is close to the critical value. Increasing the horizontal extent of the domain might allow a more strongly convecting system to develop in the lower aquifer at the most stable wavelength. The convected heat flux through the upper portion of the sulfur aquifer is presented in figure 1.5. The average value, integrated over the length of the domain, is about 2 W/m 2 . This is of the same order of magnitude as the global surface heat flow of 2.5 W/m 2 . This convected heat flux is highly dependent upon the permability of the crust. Figure 1.6 illustrates this, showing the maximum integrated heat flux through the upper zone of the aquifer for various permeability values. The system falls below the critical Rayleigh number for permeabilities less than 2 x 10 - 1 1 m 2, and the convected heat flux becomes zero. The high-permeability pipe in model B permits greater vertical velocities in the aquifer, but there is still no flow across the high-viscosity barrier, and the aquifer re-mains divided. Figure 1.7 shows only the top of the domain for model B. A narrow, higher velocity plume develops where the permeability is greater. Chapter 1. Sulfur Volcanism on Io 21 Figure 1.4: Model A: Convection Patterns in Aquifer x (km) Temperature contours and scaled fluid velocity vectors in the aquifer. Contours are every 10 degrees from 400 K to 460 K and every 100 degrees from 500 K to 900 K. The size of the arrow represents the magnitude of the velocity. Maximum velocities are about 3.5 mm/day in the upper region. The depth is measured from the top of the aquifer, which is 3.5 km below the surface of Io. Chapter 1. Sulfur Volcanism on Io 22 Figure 1.5: Model A: Convected Heat Flux Heat flux convected through the upper layer of the aquifer. The in-tegrated value is about 2 W/m 2 . The remaining model C simulates the withdrawal of fluid from the aquifer in order to feed eruptions at the surface. In this case fluid is allowed to escape out of the top boundary from x=0 to £=100 m. Thermal buoyancy forces power the withdrawal of magma. An infinite supply of fluid is available from the boundary at £=4 km. The mass and heat fluxes at the outflow nodes are calculated for each time step and plotted in figure 1.8. During a start-up period of about 15 years, the exit temperature increases from 393 K to a fairly steady 429 K as heat is advected towards an initially cold outflow region. Exit velocities increase by an amount equal to the increased buoyancy force due to the new temperature. The response of the heat flux is magnified because of the changing temperature. Mass and heat fluxes decrease in tandem as a steady state exit velocity is approached due to the establishment of a steady driving pressure field. The mass flux levels off at about 2.5 x 10~4 kg/s/m2, with an advected heat flux of 125 W/m 2 . Chapter 1. Sulfur Volcanism on Io 23 Figure 1.6: Influence of Permeability on Convected Heat Flux 1 2 3 4 5 k(m 2) x10"11 Maximum average heat flux convected through the upper layer of the aquifer for different values of the crustal permeability. 1.5 Discussion 1.5.1 Thermal stratification of the aquifer Within the confines of the chosen parameters, the result of model A supports the assump-tion by McEwen and Soderblom (1983) that two zones of low-viscosity sulfur would exist at different depths in a porous crust and would not be thermally mixed. This separation then requires that some mechanism exist to tap the deep region of hot, low-viscosity fluid and feed Pele-type eruptions according to the McEwen and Soderblom (1983) hypothesis. The role of crustal heterogeneity in this regard is examined using model B. A high-permeability pipe extending from the top of the aquifer down below the 460 K viscosity barrier does not succeed in bringing hot sulfur to the upper boundary of the domain. If fluid were to be allowed to escape from the top of the aquifer in model B, the higher permeability would simply enhance outflow velocities for sulfur at 400 K from the upper zone. Hotter fluid would not be drawn from the deeper zone. This conclusion would probably not change qualitatively for a large range of permeability contrasts, as long as Chapter 1. Sulfur Volcanism on Io 24 Figure 1.7: Model B: High-Permeability "Pipe" x(km) Top kilometre of the aquifer showing the effect of high permeability from x=0 to a;=100 m. an infinite supply of liquid sulfur is assumed at the right-hand boundary x=6 km. An alternative to crustal heterogeneity for bringing hot sulfur to the top of the reser-voir is the suppression of the viscosity maximum. This has been documented to occur with the addition of both halogens and H2S to liquid sulfur. Chang and Jhon (1982) de-veloped theoretical expressions for the viscosity of sulfur with halogen impurities which agree well with experimental results by Fanelli (Fanelli,1946, in Chang and Jhon,1982). The viscosity increase after 432 K is greatly damped, rising to less than 0.02 Pa-s at 460 K, compared to the pure sulfur viscosity of almost 100 Pa-s. The viscosity maximum for halogen-doped sulfur was not reached within the temperature range investigated by Chang and Jhon (1982), but at 520 K the viscosity was still under 0.06 Pa-s for all the halogens and concentrations involved. Although there is no reason to suspect that halogens are present on Io, Nash and Howell (1989) have proposed a similarity between the infrared reflection spectrum of Io and laboratory spectra of hydrogen sulfide frost at 77 K and 10 - 5 atm. Wiewiorowski and Touro (1966) calculated, in agreement with experimental data, the viscosity of liquid sulfur containing H2S, based on the average chain length in the polymerized state. As Chapter 1. Sulfur Volcanism on Io 25 0 i • 1 ' 0 5 0 100 1 5 0 2 0 0 t (years) Model (C). (a) Mass flux at upper boundary outflow, (b) Correspond-ing advected heat flux. with halogens, the chain termination reactions which take place drastically reduce the viscosity. At 510 K, sulfur saturated with H2S at 1 atm has a viscosity of only 0.125 Pa-s. Since the polymerization temperature of pure sulfur does not change by more than 20 K from atmospheric pressure to 60 MPa (10 km depth) (Vezzoli et al.,1969), the sup-pression of the viscosity maximum should remain stable with depth as well. This would most likely enable whole-domain convection to develop. In this case the aquifer would no longer be thermally stratified, defeating the hypothesis of McEwen and Soderblom (1983). Hot sulfur eruptions may be fed from thermal plumes in an entire, well-mixed reservoir. Lower-temperature eruptions, which may still be heated to nearly 400 K by the sulfur aquifer, might be fed solely by S0 2 circulating in a similar manner at a shallower depth. Chapter 1. Sulfur Volcanism on Io 26 1.5.2 Mass and heat fluxes The model mass and heat fluxes at the outflow regions are many orders of magnitude less than the observed and calculated fluxes at the surface of Io. The first conclusion to be drawn is that the sulfur aquifer, as modelled, cannot directly supply the surface eruptions. Rather, sulfur magma must be collected and stored at some point between the reservoir and the surface. Eruption intervals and durations may be roughly calculated based on the supply rate of magma from the reservoir. Before comparing the model results with surface data, a reinterpretation of the resur-facing rate as calculated in section 1.2 may be in order, considering the heat flux con-vected through the upper aquifer in model A. The calculations of section 1.2, and those of Blaney et al. (1994), rely on volcanism to supply the entire measured heat flux from Io, with the result that resurfacing rates are much larger than the minimum 0.1 cm/year needed to erase the cratering record (Johnson and Soderblom,1982). This need not be the case. Although conduction through the lithosphere probably does not contribute much to global heat flow, convecting fluids in the crust may transport heat much more efficiently to the surface. The heat flux convected through the top part of the aquifer is of the same order of magnitude as the globally averaged surface heat flow (2.5 W/m 2) observed by Veeder et al. (1994), if the permeability of the crust is high enough. This indicates that if the liquid sulfur aquifer (and possibly a SO2 aquifer) is present more or less continuously around Io, velocities may be sufficient to deliver a large part of this global heat flow by fluid circulation in the crust. The energetics of resurfacing may have to be equated with a fraction of the total heat flow. Veeder et al. (1994) have reported that more than four fifths of the measured heat flux is radiated from relatively low-temperature areas under 200 K. The 20% contribution from regions with an apparent surface temperature greater than 200 K includes small, Chapter 1. Sulfur Volcanism on Io 27 10000 5000 Figure 1.9: Power Radiated from Pele 1983 84 85 86 87 88-89 89-90 92 93 Year Observed infrared power radiation (from Veeder et al.,1994). (a) 650 K model source, (b) Large <200 K source. localized hot spots associated with volcanic centers. For example, as did Pearl and Sinton (1982), Veeder et al. (1994) model Pele as a 650 K circle of average radius 6 km, superimposed on a warm area (<200 K) of radius 300 km to 600 km. The areally averaged heat fluxes over these two complimentary regions are displayed in figure 1.9. The average 20% contribution from hot spots greater than 200 K could be considered a maximum supplied by cooling magma and volcanic plumes, with the remaining 80% derived from conduction and convection within the lithosphere. The heat flux from the large, relatively low-temperature Pele component is an order of magnitude greater than that convected through the upper part of the aquifer. However, such uncertainties exist in the model parameters that it is conceivable that fluid convection in the crust of Io could deliver a large part of that power. If resurfacing accounts for only 20% of the global heat flux, then the estimates of the Chapter 1. Sulfur Volcanism on Io 28 resurfacing rates may be revised downwards. The silicate resurfacing rate approaches the lower limit based on the obliteration of the cratering record on Io (Johnson and Soderblom,1982). A sulfur eruptive mass flux would be close to 1.5 x 10 9 kg/s, with an advected heat flux of 1 x 10 1 5 W. These flux values may be greater than the actual magma mass and heat fluxes. Besides omitting any convective heat transport through the crust, the resurfacing rate calculated in section 1.2 entirely neglects solid pyroclastic material in the volcanic plumes. Equation (1.1) assumes that all material is cooled from its liquidus temperature, whereas there may be significant entrained solids which are not supplied from a magma chamber. It would be difficult to estimate the relative proportions of initially solid and liquid material ejected from a volcano, but including solid near-surface ejecta in the resurfacing calculations would diminish the required magma flow rate. For example, Johnson and Soderblom (1982) point out that an analysis of the volcanic plume Loki, from Voyager images, yielded a fine particulate mass flux of up to 10 8 kg/s, equivalent to 1 cm/year of global resurfacing during the period of eruption. Much of this may have originated from above the magma supply region. A final point in the estimation of individual eruptive fluxes is that some volcanoes may be continuously active for years at lower temperatures (e.g. 400 K S-SO2 volcanoes in McEwen and Soderblom, 1983), while others may have short-lived eruptions correspond-ing to the observed high-temperature outbursts of Veeder et al. (1994). Individual flow rates will be highly variable, and it is simplistic to apply an average value to a particular volcano or type of eruption. Factoring in thermal convection in the crust and solid pyroclastic plume material as contributors to heat flow, magmatic eruption rates are assumed to be responsible for one tenth of the global 10 1 4 W. This yields a mass flux per eruption of approximately 5 x 10 8 kg/s and a corresponding advected heat flux of 5 x 10 1 4 W. Chapter 1. Sulfur Volcanism on Io 29 The rate of magma supply from the model reservoir C is about 2.5 x 10 - 4 kg/s/m2, and the heat flux is 125 W/m 2 (figure 1.8). The area of the reservoir outflow region will determine the total mass and heat flux supplied to a magma chamber. Since the domain is symmetric, the total width of the outflow region is 200 m. Assuming a 1 km length perpendicular to this, the area is 2 x 105 m2. The total mass flux would be 50 kg/s, which is 10~7 times the surface flux. At this rate, a one day eruption, bringing 4.3 x 1013 kg of liquid sulfur to the surface, would necessitate a supply time of almost 28 000 years. This is clearly in discord with the high rate of volcanic activity observed on Io. It would take 56 000 years to resurface Io completely with the minimum 0.1 cm/year rate required by the cratering record argument. The heat flux of 2.5 x 107 W leaving the reservoir is similarly only 5 x 10 - 8 times the eruptive surface value. If all the heat energy supplied to the magma chamber is transported to the surface in the above one-day eruption, the heat flux would be 2.5 x 1014 W. Averaged over 28 000 years, this is one millionth of the global heat flow assumed to be due to magma transport. Model results and surface observations are difficult to reconcile even within orders of magnitude. This calls into doubt the importance of a liquid sulfur aquifer in terms of a global contribution to volcanic output on Io. The assumption of Darcy flow for liquid sulfur, and the porosity and permeability distributions that accompany this assumption, do not provide the required mass and heat fluxes demanded by the observed rate of activity on the surface of Io. Although this conclusion lends support to resurfacing by silicate volcanism, convection of liquid sulfur in the crust could still supply a significant fraction of the observed global heat flow. Chapter 1. Sulfur Volcanism on Io 30 1.5.3 Crustal recycling If sulfur volcanism is important on Io, the continued existence of a liquid sulfur aquifer depends upon recycling what is erupted. Burial will result in eventual remelting and reincorporation into the supply region. The same is true of any equivalent SO2 layer, which would be at a shallower depth because S0 2 melts at a lower temperature. However, the assumption that reservoirs in a porous silicate crust supply uniquely sulfur and SO2 eruptions poses problems with the fluid recycling system. If the crust of Io is mainly silicate, as suggested by McEwen and Soderblom (1983) and from the magnitude of topographic variations on the surface (Clow and Carr,1980;chapter 4), then erupted sulfur must migrate downwards, through a relatively impermeable crust, to rejoin the aquifer. In the absence of silicate volcanism, continued burial by other sulfur and SO2 eruptions would tend to cause a sulfur-rich layer to accumulate over a silicate subcrust. This is the scenario envisaged by Smith et al. (1979) but questioned on the basis of inconsistency with the surface topography. If both sulfur (and SO2) and silicate volcanism coexist, the entire lithosphere would be extensively recycled. Burial of sulfur on the surface of Io by silicate lava flows and pyroclastic ejecta would increase the likelihood of returning sulfur to its melting depth in the lithosphere. Moreover, a young crust of fractured silicate lava flows and pyroclastics would be more permeable than an older, passive silicate crust. 1.6 Conclusions Although parameter uncertainties are significant, the result of modelling a sulfur magma reservoir according to the hypothesis of McEwen and Soderblom (1983) - as an aquifer of pure liquid sulfur within a porous silicate crust - suggests that mass and heat fluxes are insufficient to power the volcanic activity observed on the surface of Io. Chapter 1. Sulfur Volcanism on Io 31 The high-viscosity zone centered around 460 K impedes the transport of sulfur from the lower regions of the aquifer to the upper boundary. This remains the case when a high-permeability pipe extends to depth below the outflow area. Thus some mechanism other than purely thermal convection would need to exist in order to tap the deeper region of the reservoir and supply hot sulfur to a volcanic eruption. If impurities in the sulfur such as H2S suppress the viscosity maximum sufficiently, whole-aquifer convection would develop, and this would enable hot sulfur to rise and feed a magma chamber. While results indicate that the mass and heat fluxes from the reservoir are incompat-ible with those for eruptions at the surface, convection of liquid sulfur in the crust of Io may transport a significant proportion of the globally averaged heat flux. This convected heat flow is highly dependent upon the assumed permeability of the crust, and falls to zero below a critical permeability of about 2 x 10 _ n m 2. Whereas previous high resur-facing rates are based upon equating the cooling of magma alone to the global thermal budget, a lower rate of extrusion is required if crustal fluid circulation is responsible for part of the heat flow. Chapter 2 Model Construction and Benchmarking 2.1 Detailed Mathematical Model Finite element modelling constitutes a method of numerical solution which offers more flexibility, in terms of grid set-up and subsequent alteration to the equation being solved, than the more popular finite difference approach. Integral to the mathematical formu-lation is choosing the shape of the elements and the manner in which an approximate solution wi l l vary across these elements. Solutions to the equation are determined exactly at the node points in the grid, but require interpolation between the node points. The equations were solved on a rectangular grid with triangular elements, as illus-trated in figure 2.1. The n nodes are numbered along the y axis as shown, beginning in the top left hand corner. Numbering along the side of the domain with fewer nodes wi l l result in the solution of matrices with a smaller bandwidth. The advantage in using a rectangular grid is the ease with which the grid may be set up. The drawback lies in the fact that spacing is constant between two adjacent coordinates. Thus i f small elements are required in the vicinity of some (#;,?/;), this narrow spacing is propagated along all values of x and y throughout the grid. This may result in an unnecessarily large number of nodes and computation time. 32 Chapter 2. Model Construction and Benchmarking 33 Figure 2.1: Finite Element Grid The finite element grid has if node points along the vertical. The shaded region illustrates the bandwidth of the matrices to be solved. A node i is attached to 6 other nodes between i — if and i + if. This results in a bandwidth of 2if + 1. 2.1.1 Basis Functions The use of triangular elements provides the simplest approach to the construction of a finite element model. Linear interpolation within triangular elements can be combined with a dense grid for any desired accuracy, and results in the solution of matrices with the smallest possible bandwidth. An arbitrary triangular element has node coordinates (ajj,j/j). The solution of an unknown function £ is sought at the node points, using the approximation that it varies linearly across the element. Within the element, this approximate function is denoted by £, to distinguish it from the true function. It is a linear combination of the values of Chapter 2. Model Construction and Benchmarking 34 £ at trie node points. £{x,y,t) = Ni(x,y)ti(t) On a node i, Nj and Nk are zero, N{ is unity. In a global coordinate system £ = ax + by + c Since the approximate solution is equal to the true solution at the node points, C(xi, Vi) = 6 = axi + byi + c This provides the following three equations: xi yi 1 X2 2/2 1 \ x3 y3 1 J a b \ c ) 6 6 Solution of the above yields the cooefficients a, b, c of the linear interpolation function £. The computation was carried out symbolically using Mathematica 1 software. 1 a 2A 1 [(2/3 - 2/2)6 + (2/1 - 2/3)6 + (2/2 - 2/1)6] C = ^ X 3 V 2 ~ X 2 j / 3 ^ 1 + ( X l V 3 ~ a ; 32 / i )6 + (^22/1 - ^12/2)6] where A is the area of the element. Defining the following differences, « i = 2/3 - 2/2 «2 = 2/i - 2/3 « 3 = 2/2 - 2/1 B\ = x2 - x3 32 = x3 - xi 03 = X l - x2 7i = 3^2/2 - Z22/3 72 = a^iys - 332/1 73 = E2J/1 - xxy2 and substituting the expressions for a, fc, c back into £ = as + by + c, the basis functions iV,- are ^ = 2A (aiX + ^ + ^ (2.1) 1 Registered trademark of Wolfram Research, Inc. Chapter 2. Model Construction and Benchmarking 35 2.1.2 Finite Element Equations This section casts into finite element form a linear, second order, partial differential equation with constant coefficients. d£ d£ d£ d2i d2£ cot + c i ^ r + c 2— + c3— + c 4 — + c 5 — + c6 = 0 (2.2) ot ox dy dx2 dy2 Non-hnearity (e.g. Ci = Q(£) ) is dealt with by iterative solution, assuming each time a constant, but updated, coefficient Cj. Equation (2.2) is a linear operator on £ which will be denoted L(£). The approximate solution does not in general satisfy this equation, and the finite element approach seeks to minimize the residual L(i) - L{t) = L(i) ± 0 by requiring the integral over the entire domain JJWi L(£) dx dy = 0 , i = l,...n (2.3) The u>i are arbitrary weighting functions for each of the n node points. If the u>i are all equal, the system of equations is underdetermined, with fewer equations than un-knowns. Point collocation assumes u>i = 8(x — Xi)(y — yi) and so L(ti)—0 at every node point. However, with this approach the solution is sensitive to the positions of the nodes. Choosing these weighting functions to be equal to the basis functions Ni derived in the previous section is known as the Galerkin method. The residual is then made zero, in a piecewise sense, over each small domain defined by the basis function Ni. A full expansion is carried out using equation (2.2) with the summation expression £ = Nj£j and the basis functions from equation (2.1). The c0 term becomes JJ NiCo £ dx dy = c0 & JJ Ni Nj dx dy = C o ^ r r ^ ( 2 - 4 ) Chapter 2. Model Construction and Benchmarking 36 where A is again the area. The time derivative is dealt with in a finite difference manner and is considered constant during this spatial finite element expansion. JJ Ntd^dxdy = JJ NiNjdxdy = Cl^T—d7 (2-5) The spatial first derivatives are 9i j_ j.. _ _ t ff ,r W JJNiC2-^dxdy = c2 & JJ Ni-^-dxdy = c*ti^JJ Nidxdy cti A f c 2 a C2^j2A3 and similarly Hi (2-6) Green's identity must be invoked to ehminate the second derivative in the spatial terms. Otherwise, the choice of linear interpolation between node points produces a second derivative equal to zero. JJNiCtUixdy = CiJ N<%iy-CiJJltfx dxiy In both equations (2.8) and (2.9) the integral around the element boundary represents the normal flux of the quantity £ across that boundary. When these integrals are added Chapter 2. Model Construction and Benchmarking 37 together, all fluxes across shared element boundaries cancel out, leaving only the flux into or out of the outside boundary of the entire domain. If this is a known quantity it becomes a specified boundary condition. The final constant term is Nic6dxdy = ^ (2.10) The above expressions (2.4) through (2.10) are grouped into temporal, spatial, and con-stant terms. cx (A + ASjj) djj Ic0 (A + ASjj) c2 ctj c30j _ c 4 a,- ctj ^ c 5 Bj Bj \ 12 dt \ 12 6 6 4 A ~~ 4 A ) ^ + / ( e i W | | * + „ M | < t a ) + £ L d = 0 (2.11) Equation (2.11) is known as the element stiffness matrix equation. The entries i n each of these 3 x 3 matrices must be correctly placed into the global equation according to the numbering scheme introduced in figure 2.1. If the nodes 1,2,3 of the triangular element are numbered k in the global scheme, then an entry a\2 in an element stiffness matrix, for example, is placed at the position in the corresponding global matrix. The final entry i n the global stiffness matr ix is the sum of all contributions from each element stiffness equation. The global equation contains temporal, spatial, and constant coefficient matrices T^, Sij, and qi, respectively. T^ ^ + Sij ^j + qi = 0 The vector qi contains the boundary integrals from equations (2.8) and (2.9) as specified fluxes if they are known, along with the terms associated with the constant CQ. The finite difference treatment of the time derivative, in an implicit sense, allows all constant terms and those known from the previous time step to be grouped together, wi th the following rearrangement: ^j - $ + S i j ^ A t + * i + A t = o Chapter 2. Model Construction and Benchmarking 38 {hT-+Si>) &At+^+At - hTi& = 0 ^•^r A i = /*fc- ) ( 2 - 1 2 ) The right hand side of equation (2.12) is known, and the coefficient matrix Aij is updated at each time step. 2.2 Validation Tests In order to test the accuracy of the numerical model, results were benchmarked against analytic solutions for one-dimensional heat conduction and advection, and the critical Rayleigh number for two-dimensional thermal convection. All tests were carried out with constant property fluids in homogeneous porous media. 2.2.1 One-Dimensional Thermal Conduction The finite element code was tested against the analytic solution of the one-dimensional heat conduction equation dT d2T ~dt~Kdx^ = 0 ' *e[o>°°) (2-13) with the boundary conditions T(a;,0) = 100 T(0,t) = 150 T(oo,t) = 100 The value of K is 10~6. The analytic solution is derived, as in Turcotte and Schubert (1982), by means of a similarity variable. Using x V = Kt Chapter 2. Model Construction and Benchmarking 39 Figure 2.2: One Dimensional Thermal Conduction 100 ana ly t ic + + + n u m e r i c a l o 10 20 30 x(m) Analytic and numerical solutions after 100 time steps. I l l I I I I l l l 40 50 equation (2.13) becomes the ordinary differential equation d2T dT _ dry2 dn with the solution T = T(0, t) - (r(0, t) - T(x, 0)) erf x The solution was advanced one hundred time steps from the initial conditions on a 100 m grid of constant 1 m spacing. Each time step corresponds to 106 seconds, or about 11.6 days. This is equal to the characteristic diffusion time r of heat across one grid cell, r ~ L21K. The simulation time is such that the temperature disturbance never reaches the boundary at x = 100 m, satisfying the condition at x — oo. Figure 2.2 shows Chapter 2. Model Construction and Benchmarking 40 the simulation results and comparisons with the analytic solution. Agreement between numerical and analytic solutions is close at all times, and may be improved by employing a smaller time step. 2.2.2 One-Dimensional Thermal Advection The equation describing the transport of heat by both advection and conduction in a constant property flow field in one dimension is —- + U — - K — - = 0, xeO.oo 2.14) ot ox ox1 The same boundary conditions were chosen as for the conduction case. T(x,0) = 100 T(0,t) = 150 T(oo,<) = 100 Ogata and Banks (1961) derived the standard analytic solution to equation (2.14). The numerical solution was again calculated up to one hundred time steps from the initial conditions on a 100 m grid of constant 1 m spacing. Grid spacing is subject to the condition that the Peclet number Pe < 2, and the time step of 1 x 106 seconds is chosen so that the Courant number Ct < 1. The boundary at a;=100 m is sufficiently far to accomodate the condition at x = oo without influencing the solution. Results are shown in figure 2.3 for an advective velocity of 5 x 10 - 7 m/s. Numerical dispersion along the advective front occurs in the model calculations. This is a ubiqui-tous phenomenon due to the truncated Taylor series employed in defining the temporal derivative in a finite difference sense. Reduction of the time step results in improved Chapter 2. Model Construction and Benchmarking 41 Figure 2.3: One Dimensional Thermal Advect ion x(m) (a) Comparison of numerical and analytic solutions after 100 time steps. (b) Numerical solution at same time as (a) using one tenth the time step size. Chapter 2. Model Construction and Benchmarking 42 agreement between the numerical and analytic solutions. Figure 2.3b shows that a suffi-cient reduction in the time step can nearly eliminate numerical dispersion. 2.2.3 Two-Dimens ional The rma l Convect ion No analytic solutions exist to describe time-dependent convective flow. However, for a small perturbation to the conductive state of a fluid layer,a linear stability analysis can determine whether the perturbation grows or decays exponentially in time. From this the well-known relation between the critical Rayleigh number for the onset of convection and the aspect ratio of the disturbance is constructed. For a fluid with constant properties in an infinite layer of a homogeneous, porous medium heated from below, the Rayleigh number is defined as R a = §gckbAT for a layer of height b with a temperature difference of AT between the top and bottom boundaries. The numerical model was tested by comparison against the theoretical plot of critical Rayleigh number versus non-dimensional wavelength of the perturbation. Figure 2.4 shows that results are in agreement with theory. The theoretical curve is given by * ^ 2 + - ) 3 2-TTfe2 A The minimum critical Rayleigh number for convection in a porous medium is 39.48, at a wavelength corresponding to \=2b. The numerical calculations were carried out with a cosine temperature perturbation in a closed box one wavelength long, at selected values of Ra for various box aspect ratios. No heat flux was allowed out of the ends of the box, in order to simulate identical conditions repeated along an infinite layer, and the top and bottom were kept at constant temperatures. The initial response of the perturbation is plotted. Critical Rayleigh number theory does not predict the time evolution towards Chapter 2. Model Construction and Benchmarking 43 Figure 2.4: Critical Rayleigh Number 1 6 0 1 4 0 1 2 0 1 0 0 CO 4 0 2 0 0 -r 1 1 1 1 1 1 + pe r tu rbat ion g r o w t h o pe r tu rba t ion d e c a y -\ + + + / / \+ + o + O + o - O — -— -t o 1 1 1 1 1 4 5 6 2?rb 8 Comparison of numerical results against the theoretical curve of Racr. finite amplitude convective cells, and in general perturbation wavelengths revert to A=26, the most stable wavelength at the minimum critical Rayleigh number. 2.3 Summary The finite element method is used to provide a numerical solution to the equations of thermal convection in a porous medium. Included in the non-linear equations are varia-tions with temperature of the fluid density, heat capacity, and viscosity. Heterogeneous permeability and porosity are allowed in the porous medium. When tested with constant Chapter 2. Model Construction and Benchmarking 44 properties, the models agree with established analytic solutions for one-dimensional heat conduction and advection, and two-dimensional convection in an infinite layer. Chapter 3 Numerical Experiments with Variable Fluid Properties The viscosity and heat capacity of l iquid sulfur both display a strong temperature de-pendence (figure 1.3). Numerical experiments were performed in order to isolate the effect on convective behaviour of these two fluid properties and provide a better under-standing of the contribution of each to the dynamics of the l iquid sulfur aquifer. A l l runs were carried out in a closed box of aspect ratio 2, which corresponds to the wave-length parameter of the minimum critical Rayleigh number Racr for a fluid wi th constant properties save density (figure 2.4). This case is chosen as the benchmark against which to compare the following models and thus determine i f the variation in a fluid property increases or decreases the effective Rayleigh number Raeff- A measure of the evolution of a perturbation is made in terms of a normalized summation V of velocity magnitudes over the entire domain.The top and bottom boundaries of the box remain at a constant temperature and there is no heat flux out of the ends. The mass flux at all boundaries is zero. Figure 3.1 illustrates the time evolution of a cosine temperature perturbation for a constant property fluid at the theoretical value RaCT=¥). Real physical values of temperature, length, and fluid properties have been normalized for all models according to T' = where T0 = T(b/2) 45 Chapter 3. Numerical Experiments with Variable Fluid Properties 46 Figure 3.1: Constant Property F lu id at Racr (a) T(t=0) (b) T(t=400) 1i 1 1, (c) Evolution of Perturbation 1.011 > 1 (a) Init ial temperature contours, (b) Temperature contours after 400 time steps, (c) Time evolution of velocity. (T,t,x, and y refer to their primed counterparts in the text.) f = i St where b is the height of the box and St is the time step. The parameters for each model in the following sections are such that, based on v(T0) and c(T0), Ra=A0. The porosity of the medium is 10%. This does not enter into the calculation of Ra, but determines the weighting of fluid and solid properties in a material average. Chapter 3. Numerical Experiments with Variable Fluid Properties 47 3.1 Heat Capacity The high viscosity barrier at 432 K effectively isolates an upper low viscosity zone in the sulfur aquifer. This layer of fluid is fairly accurately modelled as having a constant viscosity and a lower impermeable boundary at a constant temperature of 432 K, the polymerization temperature. Figure 1.3 illustrates the relative positions of the maxima in viscosity and heat capacity. The steep positive slope in the heat capacity hes entirely within the low viscosity temperature range, and it is therefore instructive to examine the effect that such a variation in heat capacity produces in a fluid with otherwise constant properties. Heat capacity is a measure of the ability of a substance to contain thermal energy. An incremental change in this thermal energy content, or enthalpy, per unit mass is expressed as dh = cdT Comparing an equal enthalpy change Ah for a fluid particle with a constant heat capacity c0 to one with varying heat capacity c(T), which at initial conditions have equal values of c and T, j7"1 c0dT = fc(T)dT (3.1) JTo JTO A first order Taylor series expansion for the heat capacity (assuming a linear change over a small A T ) is c(T) = C o + J£(T - T0) Integrating (3.1) produces C o A r c s t = c0(T - To) + l-^(T - To)2 The ratio of the temperature response of the two fluids is Chapter 3. Numerical Experiments with Variable Fluid Properties 48 Table 3.1: Effect of Variable Heat Capacity on Fluid Velocity dc/dT < 0 dc/dT > 0 A T > 0 fluid sinks . . . A T > ATcst more slowly A T < A T c s t more quickly A T < 0 fluid rises . . . A T < ATcst more quickly A T > A T c s t more slowly for AT=T — TQ. If this ratio is greater than one, i.e. AT < ATcat, the fluid particle with a variable heat capacity equilibrates less quickly with the surrounding temperature field. This translates into a greater density difference and a consequently greater buoyancy force (positive or negative) which drives higher velocities. If the ratio (3.2) is less than one, temperature equilibration proceeds more quickly for a variable heat capacity fluid, and the buoyancy force is decreased. Whether this ratio is greater or less than one depends upon both the sign of A T and that of the slope of the heat capacity. Table 3.1 illustrates the competing effects of the buoyancy forces for rising and sinking fluid, as a function of the sign of the heat capacity slope. The overall effect on velocities is equivalent, whatever the direction of the heat capacity variation. Thus it seems that Raeff does not depend on the sign of the heat capacity slope. The ratio (3.2) can be rewritten in terms of the normalized variables T' and c'. AT,,,. . 1 = 1 + - A T ' ( — A T 2 I dT' (3.3) Since A T ' is small for this analysis, the effect of a variable heat capacity relies upon the magnitude of its normalized slope. If changes in the heat capacity of the fluid are sizable compared to the magnitude of the average value cQ, this will be reflected to a larger extent in the thermal buoyancy. Figure 3.2 shows the numerical results of a positive slope in the heat capacity function. Chapter 3. Numerical Experiments with Variable Fluid Properties 49 Figure 3.2: Effect of Variable Heat Capacity on Fluid Velocities Heat Capacity Evolution of Perturbation 0.8 1 1.2 0 400 T t The normalized heat capacity varies by 30%, 60%, and 120% respect-ively in a, b, and c. (T,c, and t refer to their primed counterparts in the text.) Results are identical for the equivalent negative slope of each figure, true to the prediction of table 3.1. The non-normalized slope of the heat capacity is the same for each case. A different reference value c0, however, results in different normalized slopes. Although Raeff is increased in each case, the effect of a variable heat capacity remains a second order contribution to fluid convection, and the final temperature contours are still quite similar to those of figure 3.1. For a variation in c' of 30% around some mean, Raeff is negligibly increased. When the variation is about 120%, fluid velocities increase by a Chapter 3. Numerical Experiments with Variable Fluid Properties 50 Figure 3.3: Effect of Variable Heat Capacity Away From Ra, (a) Constant Heat Capacity (b) Variable Heat Capacity 0 400 0 400 t t (a) Time evolution of perturbation, (b) Time evolution of perturbation for c varying by 120%. Ra/Racr=1.06 for both cases, (t refers to t' in the text.) factor of two at the steady state. This result is true whether the temperature dependence of the heat capacity is positive or negative. The effect of a variable heat capacity is most pronounced if the system is near Racr. At Ra only 6% higher than the critical value, steady state velocities are nearly equal for the cases of constant and variable heat capacities. Figure 3.3 illustrates this for the largest variation in c (120%). The difference in velocities is less yet for a smaller temperature dependence in the heat capacity. 3.2 Viscosity The fourfold viscosity increase due to polymerization dominates overall flow patterns in the liquid sulfur aquifer. As opposed to the effect of a variable heat capacity, changes in viscosity are straightforwardly reflected in the fluid velocity through Darcy's law k Vi = p,i where is the total pressure force including the buoyancy force. Relative to the heat capacity, a temperature dependent viscosity has more effect on Ra>eff- Figure 3.4 shows Chapter 3. Numerical Experiments with Variable Fluid Properties 51 Figure 3.4: Effect of Variable Viscosity on Fluid Velocities Kinematic Viscosity Evolution of Perturbation 0.8 1 1.2 0 400 T t The normalized viscosity varies by 30%, 60%, and 120% respectively in a, b, and c. (T, u, and t refer to their primed counterparts in the text.) the result of variations in the viscosity from 30% to 120%. Numerical results are identical for the cases of positive and negative slopes in the temperature dependence, as with the heat capacity. A variation of 30% in u causes a 50% increase in velocities at Racr, and a 120% variation in v leads to an increase of 4.5 times in the velocity. Differences between constant and variable viscosity cases become less pronounced away from the critical Rayleigh number. Figure 3.5 illustrates this for a 6% departure from Racr and a 120% variation in v. As v varies less, this difference becomes smaller also. Chapter 3. Numerical Experiments with Variable Fluid Properties 52 6 Figure 3.5: Effect of Variable Viscosity Away From Racr (a) Constant Viscosity (b) Variable Viscosity 1 6i : > 3 > 3 0 0 0 400 0 400 (a) Time evolution of perturbation for constant kinematic viscosity. (b) Time evolution of perturbation for kinematic viscosity varying by 120%. Ra/Racr —1.06 for both cases, (t refers to t' in the text.) 3.3 S u m m a r y a n d A p p l i c a t i o n t o Io A temperature dependence in the heat capacity of a fluid may appreciably influence the effective Rayleigh number of a system if variations in the heat capacity are large compared to the average value. When referenced to the case of a constant property fluid at Racr in a medium with 10% porosity: 1. A variation with temperature of 30% in the heat capacity increases fluid velocities by only about 5%. 2. A variation with temperature of 120% in the heat capacity increases fluid velocities by a factor of two. 3. Positive or negative slopes of equal magnitude have an identical effect and result in exactly the same velocity increase. Away from Racr of the system, a variable heat capacity has much less effect on velocities. A temperature dependence in viscosity of the same order of magnitude, whether positive or negative, has a greater effect on Raea. Chapter 3. Numerical Experiments with Variable Fluid Properties 53 1. A variation with temperature of 30% in the viscosity increases fluid velocities by about 50%. 2. A variation with temperature of 120% in the viscosity increases fluid velocities by a factor of 4.5. 3. Positive or negative slopes of equal magnitude have an identical effect and result in exactly the same velocity increase. The above results illustrate that the sudden increase in the heat capacity of sulfur, from about 1000 J/kg/K to 1500 J/kg/K between 393 K and 432 K, will promote con-vection in the upper low viscosity zone of the sulfur aquifer if that layer is near its critical Rayleigh number. However, the variation in c of about 40% relative to the average value of 1250 J/kg/K will result in velocity increases of at most 10% relative to a constant heat capacity system at the average value. The system could probably be modelled as having a constant heat capacity of 1250 J/kg/K. It would be inaccurate to model the system with a constant heat capacity of 1000 J/kg/K, since Ra depends strongly upon the magnitude of the heat capacity. After the onset of the polymerization reaction, the heat capacity of sulfur decreases with temperature. At this point, however, fluid motion is overwhelmed by the effect of the high viscosity in this region, from 432 K to about 600 K. Assumed sulfur properties at temperatures above 717 K are fairly constant, and the Rayleigh number will vary linearly with the depth to the bottom of the aquifer. Such an assumption is likely to be invalid at some point as higher temperatures and pressures are encountered deep in the lithosphere. The present modelling does not seek to speculate on the fluid dynamics far below the high viscosity barrier. The peak in heat capacity and the onset of polymerization in liquid sulfur occur in tandem. If impurities in the sulfur act to retard polymerization until higher temperatures Chapter 3. Numerical Experiments with Variable Fluid Properties 54 (Doi,1965 ; Chang and Jhon,1982), the upper convecting region of the aquifer wi l l be deeper. Whi le this region wi l l stil l cover the same heat capacity range, c wi l l remain at 1000 J / k g / K over a greater temperature range. This should still allow accurate modelling of the upper low-viscosity layer with a constant, but lower, average value of c. Experimental data from Kennedy and Wheeler (1984) indicates that the heat capacity peak of doped sulfur is broadened, but reaches the same maximum value wi th the addition of increasing amounts of impurities. This again would probably not alter the positive effect that the rise in heat capacity has on the effective Rayleigh number of the top low-viscosity zone. Thus it is likely that the conclusions arising in Chapter 1 from modelling a pure sulfur aquifer would not differ qualitatively i f impurities are included which change the polymerization temperature. Results would be altered only i f impurities lower the maximum viscosity values to a point where there is no longer a barrier to whole-aquifer convection. Chapter 4 Viscous Relaxation of Topography on Io 4.1 Introduction In the debate about the relative importance of sulfur versus silicates as crustal and volcanic material on Io, proponents of a sihcate crust often argue that the strength of sulfur is insufficient to maintain the observed topographic slopes and relief reported by Carr et al. (1979) of up to 10 km. Clow and Carr (1980) calculated the conditions for mechanical failure of scarps and caldera walls and concluded that it is improbable that sulfur is the main upper crustal material. Recent estimates of global surface heat flow are around 2.5 W/m 2 (e.g Veeder et al.,1993), and more than ten times higher around volcanic hotspots. But even for values less than 1 W/m 2 , Clow and Carr (1980) calculate that ductile behaviour will occur in a sulfur crust at depths much less than the kilometer scale depths of calderas on Io, causing the disappearance of such topographical variations within the span of a few months. They favour a sihcate crust with less than ten percent sulfur in order to achieve sufficient strength. Another approach to studying the feasibility of an upper crust of sulfur over the sur-face of Io, as opposed to a porous sihcate crust with interstitial sulfur, is to calculate the characteristic time for viscous relaxation of sulfur topography. This neglects mechanical failure, and so is suitable for small slope angles. 55 Chapter 4. Viscous Relaxation of Topography on Io 56 4.2 Model The crust of Io is modelled as a layer of sulfur, assumed to behave as a viscous fluid, overlying a rigid silicate lithosphere (figure 4.1). The viscosity contrast between sulfur and silicates is large enough that the two materials wi l l not deform on the same time scale. For this reason also, isostatic compensation is neglected. Topographic variations are treated as perturbations to the hydrostatic equilibrium state. A Figure 4.1: Mathematical Domain h=h 0 e ikx sulfur z=-H A n infinite, viscous layer of sulfur overlies a rigid silicate half-space at a depth H below the mean surface. Topographic variations are approximated by a sinusoidal function. 4.2.1 Solution to the viscous flow equations The governing fluid dynamic equations in two dimensions are those of mass and mo-mentum conservation. dp ^7 + (pvi),i = 0 Chapter 4. Viscous Relaxation of Topography on Io 57 and Dvi P^rjj- = P9 + TjiJ The fluid density is p, v is velocity, r is the viscous stress tensor, g is the gravitational acceleration, and t is time. For an incompressible fluid with a small Reynolds number, these become, respectively, Vi.i = 0 (4.1) and P9 + Tjij = 0 (4.2) To these is added the constitutive relation for the stress tensor Tij = -pSij + p (vitj + vjti) (4.3) wi th p being the pressure and p the viscosity. The hydrostatic equilibrium state, where velocities are zero, is given by Vi = 0 pg -p = 0 Tij Pb~ij Equations (4.1) to (4.3) are solved for the perturbations to this static state, yielding the following six equations for the perturbation terms: TxzyX ~r" rZZtZ 0 Vx,x + Vz,z = 0 l~xx = — P + 2pVx,x Txz = P>{vx,z + Vz,x) TZZ = ~P + 2pvZiZ (4.4) Chapter 4. Viscous Relaxation of Topography on Io 58 For irrotational flow, Tij=Tj{. The coordinate of horizontal distance is x, and z is the vertical coordinate. An analytical solution to this system of equations is facilitated by modelling topo-graphic variations as continuous sinusoids. The height h above a reference plane z=0 is ikx h = hQel Solutions of a similar form are sought for the perturbation quantities in (4.4). These are approximated by truncated sinusoidal series up to degree 2 in the wavenumber k. For a perturbation quantity £, the x and z dependencies are separated as £ = l^)(z)eikx + £W(z)e2ikx where £(m) is a complex function of z. Only real solutions will be kept in the end. These forms of the perturbation quantities are substituted into equations (4.4) and provide the relationships iWW + £ f £ > = 0 + Tz*« =0 ikv^ + = 0 dz dz f£> = - p W + 2p£vU (4.5) for the coefficients of the elkx terms, and a similar set of relationships for the e2lkx terms, with k replaced everywhere by 2k. Algebraic manipulation of the above results in a system of four equations for the four unknowns V ^ \ V ^ \ T ^ , and . Defining a vector Chapter 4. Viscous Relaxation of Topography on Io 59 u to be u = IT. (1) (4.6) this system can be written d_ dz U2 \ U* J or f(!) 0 k 1/p 0 —A; 0 0 0 Apk2 0 0 k 0 0 -fe 0 ™2 V y (4.7) The variables f ^ and can be determined from a combination of the elements of u. The coefficient matrix A is independent of z, and so equation (4.7) suggests a solution of the form u=ceXz where c is constant. This becomes the eigenvalue problem The characteristic equation (A-XI)c = 0 A 4 - 2k2X2 + k4 = 0 yields eigenvalues A=±Ai. Associated eigenvectors are I 1 X ( D c = - 1 2pk \ -2pk J ( 2V = 1 1 -2pk \ -2pk ) Chapter 4. Viscous Relaxation of Topography on Io 60 The remaining two eigenvectors satisfy the solutions u — ^c + ^ckz^j ekz and u = ^(4)c + ^ckzj e _ f c z , and with appropriate choices of constants to ensure orthogonal-ity, are (3), 1/2 (4) C = ' -1 /2 \ 1/2 jik \ -/** j These four solutions make up the general solution u = Y^=\ aj A similar procedure for the terms in e2lhx wi th the definition / aw ^  w = IT. (4.8) produces another set of four solutions with k replaced everywhere by Ik. Thus the entire general solution is the superposition of the two general solutions u and w. 4 s = u + w = (a,j + bj ^wj (4.9) 4.2.2 Boundary conditions The mathematical problem has thus far been linear. Non-linearity is introduced i n the formulation of the boundary conditions, which are simply the continuity of velocity and stress vectors across any interface. The problem being considered involves a stress free boundary at the surface and a no-slip velocity condition at the sulfur-silicate contact. For the deformed surface h=h0elkx, instead of applying them at the interface proper, a first-order Taylor series is used to bring the boundary conditions to a constant reference plane at the topographic mean z=0. The components of the stress vector t on some Chapter 4. Viscous Relaxation of Topography on Io 61 surface are related to the stress tensor by t{ — TjiTlj where n is the outward normal to the surface. Thus, with r 2 X = T X Z , tx(z + Az ) .= rxx(z + Az)nx + T X Z ( Z + Az)nz = (rxx(z) + A (z))nx + (TXZ(Z) + AzrXZiZ(z))nz with terms on the right hand side evaluated at z. The primed variables represent the perturbation quantities! Removing the contribution of the hydrostatic equil ibrium state, — ^x pt^x = R L n x + (p,z + TL,Z) A z n x + TX2nz + rxzzAznz Dropping the primes and eliminating products of three or more perturbation terms, tx(z + Az) — rxxnx + pgAznx + rxznz (4-10) The same treatment for tz yields tz(z + A z ) = rxznx + Tzznz + pgAznz (4-H) The sinusoidal expressions for the stress perturbations and surface normal are substituted into the above Taylor series expansions. For the surface in question, the normal is ft — — ikh0elkxx + z Care must be taken to incorporate the real part only. The real part of the quantity £ = £eikx Chapter 4. Viscous Relaxation of Topography on Io 62 can be written 2ft{£} = (eikx + (*e~ikx where £* denotes the complex conjugate. Using this expansion, equation (4.10) applied across an interface becomes l- (Af£V f e * + Arxx>e-ik* + Afffe 2 *- + Ar£>e-ik*) (-ih0keihx + ih0ke~ikx) + ^Apg [h0eikx + ih0e~ikx) [-ih0keikx + ih0ke'ikx) +Ari1Jeikx + AfHK-^ + Af£>e2ikx + Af£>*e-ikx = 0 Choosing the topography to vary symmetrically about x=0 ,i.e. h=hQ cos kx, means that h0—hl in the above equation. Collecting terms in powers of e%kx results in the following relations: A r « - A f « * = 0 [\ih0kAt£ + Af£)} eikx = 0 { - i iMAf£> - l-ih0kApgh0 + A f i 2 ) } e2ikx = 0 { - ^ o f c A f ( 2 ) * + A f W * } e - ^ = 0 {^MAfJi>* + l-ih0kApgh0 + A f i 2 > } e~2ikx = 0 (4.12) The same analysis applied to equation (4.11) for A £ 2 = 0 produces five more relations. AfJJ> - A f = 0 | ^ 0 f c A f ( 2 ) + A f « + Apgh.} eikx = 0 { - ^ M A f £ ) + A f 2 ) } e « t e = 0 {-itMAfW* + Af£> + Apyfco} e-ifc* = 0 |^ 0fcAfW* + Af£)*} e - 2 i f e a : = 0 (4.13) Chapter 4. Viscous Relaxation of Topography on Io 63 The first equations of (4.12) and (4.13) demonstrate that both A f ^ and A f £ ^ are real, and this in turn leads to A f ^ being real and the remaining coefficients Af^p al l imaginary. The ten relations (4.12) and (4.13) reduce to four boundary conditions. Af£> + i i M A f J ? = 0 (4.14) *Af£> + ^ M A f £ > + \h0kApgh0 = 0 (4.15) A r i 2 1 } + \ih0kAri2} + ApgK = 0 (4.16) .-Afg) + ^ M A f W = 0 (4.17) These equations illustrate that allowing solutions up to degree 2 in A; is equivalent to considering the product h0k to be relatively large. The analysis is therefore applicable to topography wi th a steeper slope, where (hQk)2 <C 1 but h0k~l. If terms of the order 1 0 _ 1 can be discarded, this implies that h0k= V O . l , and solutions are valid for amplitude to wavelength ratios as large as 1/20. In order to recast equations (4.14) to (4.17) into conditions for the solutions u and w, the relations (4.5) are employed. For instance, the first two boundary conditions can be writ ten iAi£ + UK—Af^J + hi0kApgh0 = 0 Taking the derivative of the second and substituting into the first, A U 3 + lh°d^AU3 = ° ( 4 - 1 8 ) an ordinary differential equation for u3. The other boundary conditions transform into 1 d2 Au4 + -hl—Aui + Apgh0 = 0 (4.19) Chapter 4. Viscous Relaxation of Topography on Io 64 Aw3 + -h20—Aw3 + -h0kApgh0 = 0 1 d2 Aw4 + -h2—Aw4 = 0 (4.20) (4.21) For the present problem, the no-slip and impermeabihty boundary conditions at z=—H are given by and the stress-free surface conditions at z = 0 by equations (4.18) to (4.21). Since the velocities in the sihcate basement and the stresses in the atmosphere are all zero, Ait; and Aw{ are equivalent to the values of U{ and Wi in the viscous sulfur layer. When the boundary conditions are evaluated at their respective values of z, they provide two sets of four equations for the four constant coefficients and bi in the general solution (4.9). Thus the solutions u and w are obtained. These are transformed into the real solutions for velocity and stress, recalling the definitions (4.6) and (4.8). Atti = Awi = Au2 = Au>2 = 0 v. X Ui(z) sin(A;a;) + Wi(z) sin(2A;x) (4.22) 1 x 2 ( 2 ) cos(A;a;) + w 2 ( z ) cos(2A;a;) (4.23) r. xz Uz(z) sin(Aia;) + w3(z) sin(2fca;) (4.24) T z z = u^(z) cos(kx) + Wi(z) cos(2kx) (4.25) Chapter 4. Viscous Relaxation of Topography on Io 65 The final step is to derive an expression for the relaxation time of the topography. This is defined as the time required for the amplitude of the perturbation to reach 1/e of its original value. The rate of change of the height h is dh dh dh —— = ——|- Vih i « —- = vz dt dt ' dt Consistent with the low Reynolds number approximation, vz is considered constant, and thus Ah = vzT = h — h0 — —h0 — h0 e where T is the characteristic relaxation time. The expression for T becomes T. = ^-(--l) (4.26) vz \e J The assumption of constant vertical velocity at x=0 was tested by time-stepping through the entire calculations and recalculating vz and h until 1/e of the original amplitude was reached. Agreement between the two methods was extremely close. The validity of equation (4.26) will deteriorate progressively below the relaxed profile as velocities begin to decrease to zero at equilibrium. 4.3 Model Parameters In order to carry out the calculations, material and topographic parameters are required. These are listed in table 4.3. The density is that used by Clow and Carr (1980). The yield stresses they use for brittle failure and plastic flow for solid sulfur are converted into equivalent viscosities through the simple relation cr where cr is the yield stress and e is the strain rate. The yield strength for brittle sulfur gives the highest viscosity (3 x 10n Pa-s) and is used in order to place an upper limit on relaxation times. Chapter 4. Viscous Relaxation of Topography on Io 66 Table 4.1: Viscous Relaxation Parameters 9 1.7 m/s2 P 2070 kg/m3 M- 3 x l 0n Pa-s hQ (mountain) 1000 m hQ (crater) -100 m H (mountain) 2000 m H (crater) 1000 m \ K I \ \ 1/20 A slope of 10° approaches the upper limit of the amplitude to wavelength range of this model. Heights of 2 km (amplitude 1 km) are typical of mountains on the surface of Io, although 10 km high relief has also been reported (Carr et al.,1979). The model is applied to small impact craters as well, using a crater depth of 100 m. The maximum depth to diameter ratio of 1:20 is imposed by the limits of the mathematics. Average depth to diameter ratios for fresh appearing bowl-shaped craters of a similar size on Ganymede are 1:6 to 1:12 (Passey and Shoemaker,1982). It may be misleading to compare impact craters on Io with those on the icy satellites: a sulfurous crust will be much softer and craters may be even deeper on Io. The depth to the rigid silicate basement is chosen as 1 km for the following reason. Since the thermal conductivity of sulfur is low (about 0.3 W/m/K, or one tenth that of silicates) the temperature gradient through a sulfur crustal layer will be very high in order to coincide with the heat flux through the silicate lithosphere. For a conductive heat flow of 0.1 W/m 2 , less than one tenth the overall global average from the surface of Io, sulfur will reach its melting point at little over 1 km depth. If-this occurs, the zero velocity boundary condition becomes a zero stress condition and flow will be greatly enhanced, diminishing relaxation times considerably. The chosen basement depth acts to provide conservative results, erring on the side of long relaxation times. This still Chapter 4. Viscous Relaxation of Topography on Io 67 Figure 4.2: Relaxation Flow Field N -4 1 i . t \ 1 \ ! \ s - * - - * ^ " " ^ "/ ~ ~ i ~ ~ r -/ 1 \ / 1 \ - ^ / I . t~ — i / i 1 N — 1 v. ~-— / 1 \ — — s 1 \ — " — / 1 s — — / t x i 1 rigid basement i i i -10 -5 0 5 x (km) 10 — Initial topography relaxed profile Scaled velocity vectors represent the magnitude and direction of viscous flow. neglects the probable ductile flow of solid sulfur at depths less than 1 km, which would also partially lubricate the bottom boundary and facilitate horizontal flow. Furthermore, the heat flow beneath volcanic structures would be still greater than that away from the hotspots, decreasing the probable depth to melting or ductile behaviour. 4.4 Resul ts and Discussion A typical viscous flow field for the calculated models is illustrated in figure 4.2 for one wavelength of topographic variation. The basement depth has been increased to allow better visualization of the velocity vectors. The sihcate subcrust is 1 km deep for the actual calculations. Relaxation profiles are reached within less than a week using the parameters from table 4.3, for both the mountains and craters. This suggests immediately that any enduring features on the surface of Io can not be composed of sulfur. Observation time during the Voyager missions was limited, and in general the same terrain was not Chapter 4. Viscous Relaxation of Topography on Io 68 photographed by the two spacecraft at a resolution which would allow comparison of landforms over time (Morrison,1982). However, in order for shield volcanoes composed of sulfur to be 2 km high, they must build up over multiple eruptions. Relaxation times on the order of a week or less would not allow successive eruptions to contribute to a single volcanic edifice. Such short relaxation times might, however, be invoked to explain the lack of impact craters on the surface, although large volcanic resurfacing rates have been the accepted explanation for this phenomenon. Global volcanic mass flux estimates based upon the ob-literation of the cratering record require a resurfacing rate of at least 0.1 cm/yr (Johnson and Soderblom, 1982). If the surface of Io is dominantly sulfur, then the resurfacing rate can be revised downwards to accomodate rapid relaxation of impact craters. The same rapid decay for elevated topography, however, makes it unlikely that silicate crust does not form a large fraction of the surface, and this fraction should preserve the cratering record. In the absence of vastly inhomogeneous global resurfacing and crustal composi-tion, extensive recycling of the crust seems the most likely candidate for explaining the lack of craters. The effect of varying physical parameters is shown in figures 4.3 to 4.5. Relaxation time is most sensitive to the viscosity of the crust. Reducing the viscosity by a factor of ten also cuts the relaxation time to one tenth. If a viscosity 7 x I O 1 0 Pa-s is calculated using the yield stresses for ductile sulfur given by Clow and Carr (1980), then craters and mountains both relax in one day. The depth to the rigid silicate basement is another important factor. Since velocities are zero along the sulfur-silicate boundary, a shallow basement inhibits rapid vertical velocities as well as spreading. Figure 4.4 illustrates the fact that if the sulfur crust is only 200 m thick, the relaxation time is increased by a factor of about eight. Thus 2 km high mountains would relax in 40 days, and the model craters in 33 days. This is still Chapter 4. Viscous Relaxation of Topography on Io 69 Figure 4.3: Dependence of Relaxation Time on Viscosity Viscosity (Pa s) x10 1 1 rapid and would explain the lack of craters in any sulfur-rich crust, but not the continued existence of raised relief. Although viscous relaxation theory predicts longer relaxation times for shorter wave-length features, this is not generally the case when a rigid basement is involved. Figure 4.5 illustrates this for three different basements depths H. For the infinite half-space problem, relaxation times vary inversely with topographic wavelength. As the silicate sub crust is moved closer to the surface, there appears a critical wavelength to amplitude ratio with a minimum relaxation time. This critical ratio is a function of the basement depth and decreases as H decreases. At shallow enough depths it disappears within the parameter range of the model, and relaxation time increases monotonically with wave-length. This can be explained by the increasing effect of the no-slip boundary condition along z=—H. As the wavelength increases in a thin viscous layer, horizontal flow be-comes relatively more important than vertical flow for moving mass. The zero velocity condition propagated upwards from the sulfur-silicate interface reduces the lateral mass flux to a relatively greater extent as the layer is made thinner. Thus wide, shallow craters Chapter 4. Viscous Relaxation of Topography on Io 70 Figure 4.4: Dependence of Relaxation Time on Depth to Basement Depth to basement/initial amplitude in a thin sulfur crust wiU last longer, but likely still less than 4 months for reasonable depth to diameter ratios. The conclusions to be made from viscous relaxation calculations do not uphold the presence of a dominantly sulfur crust on the surface of Io. This agrees with the findings of Clow and Carr (1980). Using a depth of 1 km to the sihcate subcrust, 2 km high relief will reach its relaxed profile within one week. Such rapid decay would prohibit the formation of shield volcanoes by accumulation of material from successive eruptions. The relaxation time for small craters is about the same, and increases to one month for a shallow sulfur layer 200 m thick. This can explain the apparent lack of impact craters on the surface, although a large volcanic resurfacing rate remains a better candidate for hiding the cratering record. The existence of the topographical relief observed on Io is better served by largely sihcate crustal material. This does not disallow a significant proportion of sulfur on and near the surface, but would seem to indicate that the bulk of planetary resurfacing is due to silicates. Chapter 4. Viscous Relaxation of Topography on Io 71 Figure 4.5: Dependence of Relaxation Time on Wavelength 2.5i 1 1 1 wavelength/initial amplitude Variat ion of the relaxation time with the ratio of wavelength to ini-t ial amphtude of topography for three different depths H to the rigid basement. Overview and Future Direct ions Hydrological modelling of a pure sulfur aquifer on Io as a source region for volcanoes leads to the conclusion that such a reservoir is not able to supply a sufficient volume of magma to resurface the satellite at a rate consistent with observation. Combined with the evidence from topographic features, this suggests a greater role for sihcate volcanism in resurfacing the satellite. However, fluid circulation in the crust of Io may carry a substantial fraction of the heat flow through the hthosphere. Future research should be directed towards determining what fluids are present in the crust and to what extent they contribute to the global heat flux. The adoption of common assumptions from hydrology for simulating flow through the crust of Io permits a focus on the effects of the variable properties of liquid sulfur. However, model parameters are poorly constrained and depend to a large extent upon a crust of unknown composition and structure. The use of Darcy's law to describe fluid motion may be found wanting in a quantification of flow through very heterogeneous or fractured material. An assumption of laminar velocity profiles and continual thermal equilibrium between fluid and porous matrix prohibits modelling velocities of more than a few centimeters per day. Turbulent fluid motion will not follow the hnear pressure-velocity Darcy law. The development of realistic descriptions of fracture and non-laminar flow beneath the surface of Io is crucial to any furthur investigation of the role of fluid circulation and heat transport. Since sulfur may be stable as a hquid until near the base of the hthosphere, a better understanding is needed of the behaviour of its thermodynamic properties at elevated pressures and temperatures. The very real possibility of substantial impurities in the 72 Overview and Future Directions 73 sulfur deserves attention, as this may drastically alter fluid responses to temperature variations. Studies of the composition of Io could constrain the type and amount of impurities that may be involved. Sulfur is most likely an important element on Io, and may be involved in volcanism in ways other than uniquely sulfur eruptions from sulfur reservoirs. The imminent return of new information from the Galileo spacecraft encounter with Jupiter and its satellites may offer insights into this complex planetary problem. References [1] Blaney, D.L., T.V. Johnson, D.L. Matson, and G.J. Veeder. Volcanic eruptions on Io: heat flow, resurfacing and lava composition. Icarus, 113, 220-225, 1995. [2] Carr, M.H. Silicate volcanism on Io. J. Geophys. Res., 91, 3521-3532, 1986. [3] Carr, M.H., H. Masursky, R.G. Strom, and R.J. Terrile. Volcanic features of Io. Nature, 280, 729-733, 1979. [4] Chang, M.C., and M.S. Jhon. Viscosity and thermodynamic properties of liquid sulfur. BuU. Korean Chem. Soc, 3, 133-139, 1982. [5] Clow, G.D., and M.H. Carr. Stability of sulfur slopes on Io. Icarus, 44, 268-279, 1980. [6] Doi, T. Physico-chemical properties of sulfur. Rev. Phys. Chem. Japan, 35, 1-10, 1965. [7] Fanale, F.P., T.V. Johnson, and D.L. Matson. Io: a surface evaporite deposit? Sci-ence, 186, 922-925, 1974. [8] Heath, D.R. Viscosity and density of liquid sulfur at high pressures and temperatures: implications for Earth's outer core and volcanism on Io. Bachelor's thesis, Department of Earth Sciences, University of Western Ontario, Canada, 1995. [9] Johnson, T.V., and L.A. Soderblom. Volcanic eruptions on Io: implications for surface evolution and mass loss, in Satellites of Jupiter, ed. D. Morrison. Univ. Arizona Press, Tucson, pp. 634-646, 1982. 74 References 75 [10] Jacob, C.E. Flow of groundwater, in Engineering Hydraulics, ed. H. Rouse. Wiley, New York, pp. 321-386, 1950. [11] Kennedy, S.J., and J.C. Wheeler. Heat capacity of doped liquid sulfur at the poly-merization transition: a non-classical model. J. Phys. Chem., 88, 1040-1044, 1984. [12] Kieffer, S.W. Dynamics and thermodynamics of volcanic eruptions: implications for the plumes on Io, in Satellites of Jupiter, ed. D. Morrison. Univ. Arizona Press, Tucson, pp. 647-723, 1982. [13] McEwen, A.S., and L.A. Soderblom. Two classes of volcanic plumes on Io. Icarus, 55, 191-217, 1983. [14] McEwen, A.S., D.L. Matson, T.V. Johnson, and L.A. Soderblom. Volcanic hot spots on Io: correlation with low-albedo calderas. J. Geophys. Res., 90, 12345-12379, 1985. [15] Morrison, D. Introduction to the satellites of Jupiter, in Satellites of Jupiter, ed. D. Morrison. Univ. Arizona Press, Tucson, pp. 3-43, 1982. [16] Nash, D.B., and R.R. Howell. Hydrogen sulfide on Io: evidence from telescopic and laboratory infrared spectra. Science, 244, 454-457. [17] Nash, D.B., M.H. Carr, J. Gradie, D.M. Hunten, and C.F. Yoder. Io, in Satellites, ed. J.A. Burns and M.S. Matthews. Univ. Arizona Press, Tucson, pp. 629-688, 1986. [18] Ogata, A., and R.B. Banks. A solution of the differential equation of longitudinal dispersion in porous media. USGS Professional Paper 411-A, 1961. [19] Passey, Q.R., and E.M. Shoemaker. Craters and basins on Ganymede and Callisto: morphological indicators of crustal evolution, in Satellites of Jupiter, ed. D. Morrison. Univ. Arizona Press, Tucson, pp. 379-434, 1982. References 76 [20] Peale, S.J., P. Cassen, and R.T. Reynolds. Melting of Io by tidal dissipation. Science, 203, 892-894, 1979. [21] Pearl, J . C , and W.M. Sinton. Hot spots of Io, in Satellites of Jupiter, ed. D. Mor-rison. Univ. Arizona Press, Tucson, pp. 724-755, 1982. [22] Rouse, H. Fundamental principles of flow, in Engineering Hydraulics, ed. H. Rouse. Wiley, New York, pp. 1-135, 1950. [23] Smith, B.A., E.M. Shoemaker, S.W. Kieffer, and A.F. Cook II. The role of S02 in volcanism on Io. Nature, 280, 738-743, 1979. [24] Sinton, W.M. The thermal emission spectrum of Io and a determination of the heat flux from its hot spots. J. Geophys. Res., 86, 3122-3128, 1981. [25] Touro, F.J., and T.K. Wiewiorowski. Viscosity-chain length relationship in molten sulfur systems. J. Phys. Chem., 70, 239-241, 1966. [26] Tritton, D.J. Physical Fluid Dynamics. 2nd ed, Oxford Univ. Press, New York, 519 pp, 1988. [27] TuUer, W.N. The Sulfur Data Book. McGraw-Hill, New York, 1954. [28] Turcotte, D.L., and G. Schubert. Geodynamics. Wiley, New York, 450 pp, 1982. [29] Veeder, G.J., D.L. Matson, T.V. Johnson, D.L. Blaney, and J.D. Goguen. Io's heat flow from infrared radiometry: 1983-1993. J. Geophys. Res., 99, 17095-17162, 1994. [30] West, E.D. The heat capacity of sulfur from 25 to 450°, the heats and temperatures of transition and fusion. J. Amer. Chem. Soc, 81, 29-37, 1959. [31] West, J.R. Thermodynamic properties of sulfur. Indust. Eng. Chem., 42, 713-718, 1950. References 77 [32] Wiewiorowski, T.K., and F.J. Touro. The sulfur-hydrogen sulfide system. J. Phys. Chem., 70, 234-238, 1966. [33] Young, A.T. No sulfur flows on Io. Icarus, 58, 197-226, 1986. Appendix A Liquid Sulfur Thermodynamic Data and Equations of State Thermodynamic data for the temperature dependence of the properties of liquid sulfur were fit with a minimum of splined functions to ensure maximum continuity in the derivatives of the functions. The data for the density of liquid sulfur are from Tuller (1954), listed in table A . l and plotted together with the model fit in figure A . l . The values for the heat capacity at constant pressure in table A.2 are a least squares fit by West (1959) to his experimental data. The fitted curve is plotted in figure A.2. The dynamic viscosity /z is divided by the extrapolated density at each corresponding temperature point and transformed into the kinematic viscosity v. Table A.3 contains the original viscosity data (Tuller,1954), and v is plotted alongside its mathematical representation in figure A.3. 78 Appendix A. Liquid Sulfur Thermodynamic Data and Equations of State 79 Ta ble A. l : Density T(K) P (kg/m3) T P T P 394.26 1803.7 416.48 1784.2 438.15 1771.4 397.04 1800.7 419.26 1781.8 444.43 1770.5 399.82 1798.1 422.04 1779.5 451.43 1767.1 402.59 1795.7 424.82 1777.3 457.15 1764.4 405.37 1793.5 427.59 1775.2 483.15 1750.9 408.15 1791.2 430.04 1773.9 512.65 1732.9 410.93 1788.8 431.65 1772.9 551.65 1709.6 413.71 1786.4 434.15 1772.3 630.15 1658.3 718.15 1606.0 p(r)(kg /m 3 ) = 2032 - 0.59T Figure A. l : Fit to Density Data I L_I i i i i : i i 400 450 500 550 600 650 700 T(K) o Data from TuUer (1954) — Fit Appendix A. Liquid Sulfur Thermodynamic Data and Equations of State Table A.2: Sped ic Heat Capacity at Constan' Pressure T(C) c (J/mol/K) T c T c T c 115.207 31.710 162 48.03 250 36.774 360 33.712 120 31.983 163 47.34 260 36.349 370 33.562 130 32.534 164 46.81 270 35.956 380 33.417 140 33.203 170 44.603 280 35.596 390 33.272 150 34.454 180 42.453 290 35.267 400 33.121 156 38.28 190 41.097 300 34.969 410 32.957 157 40.97 200 40.025 310 34.700 420 32.774 158 44.45 210 39.144 320 34.459 430 32.564 159 48.50 220 38.415 330 34.242 440 32.319 160 49.24 230 37.804 340 34.048 444.60 32.193 161 48.68 240 37.272 350 33.873 cp(T < 432)(J/kg/K) = 220.1e(r-425) /8 + 1000 cp{T > 432)(J/kg/K) = 508e-( r- 4 3 2) / 5 9 + 1020 o Data from West (1959) — Fit Appendix A. Liquid Sulfur Thermodynamic Data and Equations of State Table A.3: Dynamic Viscosity T(K) /x (Pa-s) T T 399.82 0.00997087 438.71 20.8346 510.93 24.5551 405.37 0.00892913 444.26 46.1339 522.04 16.3701 410.93 0.00803622 455.37 89.2913 533.15 11.1614 416.48 0.00744094 460.93 93.1606 544.26 7.44094 422.04 0.00669685 466.48 90.7795 555.37 5.20866 427.59 0.00654803 472.04 81.8504 566.48 3.42283 430.37 0.00669685 477.59 72.9213 577.59 2.23228 432.04 0.0119055 488.71 56.5512 588.71 1.78583 435.93 11.9055 499.82 38.6929 v(T < 432)(m2/s) = 5 x 10 - 6 / m „ „ « w 9 / x | Y T - 432\ , / T - 4 3 0 M 2 (T-432) „ u(T > 432)(m2/s) = sech + ^ - 1 + 5 x v ~ LV 76 J 28 4 x 10" Appendix B Boussinesq Approximation Certain conditions must be fulfilled in order to approximate a fluid as incompressible in fluid thermodynamic equations. The requirement that density changes remain small with respect to the scale of the problem is reflected in the following two criteria (Tritton,1988): 1. Either the temperature dependence of the density or the magnitude of typical density perturbations must be small. where Ap is the maximum variation of density across the length scale L of the problem. This ratio is about 0.013 in the upper low-viscosity zone of the sulfur aquifer and 0.19 for the deeper zone of low-viscosity fluid. 2. The length scale L of the problem must be small compared with the length scale over which pressure effects on density are significant. LgpB < 1 (B.l) for an isothermal compressibility B. The expression (gpB)-1 is referred to as the isothermal pressure scale height. The density of liquid sulfur at 1 GPa and 623 K is about 1500 kg/m3 (Heath,1995). At the same temperature and atmospheric pressure it is 1665 kg/m3, giving an approximate value for B of IO - 1 0 P a - 1 . For the problem examined in chapter 1, LgpB=0.003. Equivalent to equation (B.l) is the 82 Appendix B. Boussinesq Approximation 83 constraint that the driving pressure forces due to buoyancy must be much smaller than the total pressure variations across the domain. Ape? < 1 The neglect of any viscous heating terms requires that heat loss due to viscous energy dissipation must be much less than that due to advection. The dissipation of viscous energy is proportional to the viscosity and to the square of the velocity gradient Darcy's equation assumes laminar flow in pore capillaries, which minimizes viscous effects. In free convection there is a balance between the viscous force and the buoyancy force ^ ~ gpaAT where 8 represents the length scale for the importance of viscous effects. Replacing the viscous heating expression by gpocATv and comparing to the advection term, gpctATv gaL ~ <C l pCpvVT cp The sulfur aquifer model satisfies this with a ratio of 0.006. The quantity ga/cp is often called the adiabatic temperature scale height. Appendix C A Note on the Derivation of the Heat Conservation Equation Fluid convection problems involve the simultaneous solution of the equations of mass, momentum, and energy conservation. The advent of powerful computers has made pos-sible the numerical simulation of complex flow problems, including fluids with variable thermodynamic properties. In treating fluids with a variable density and heat capac-ity, it has come to my attention that a basic inconsistency between various scientific publications exists in the derivation of the heat conservation equation. The Boussinesq approximation is a widely used approach to the problem of convec-tion, one which assumes a constant fluid density except in the buoyancy force. This approximation also neglects variations in the heat capacity, which is reasonable for many fluids. The use of the Boussinesq approximation sidesteps the erroneous derivation of the heat equation that is present in some publications. The problem manifests itself only when a variable heat capacity is included. Heat conservation equation The equation describing the conservation of energy can be derived in either an inertial or Langrangian reference frame. Both methods of course lead to the same result. The following derivation is performed in the context of a control volume V of fluid in an inertial frame. The rate of change of energy of the fluid in the volume must be balanced by the rate at which work is done both on the boundary of V by surface forces and on the fluid in V by body forces, and by the rate at which heat is conducted and convected across 84 Appendix C. A Note on the Derivation of the Heat Conservation Equation 85 V. For the purpose of illustrating the error discovered in certain texts, the work terms are neglected. This is akin to neglecting viscous dissipation and fluid compressibility, and will focus attention on the targeted terms in the heat equation and leave a simpler equation to manipulate. This simplified energy equation, in conventional summation notation, is d Jvpedr- Js kTjdSj + J^evjdSj = 0 (C. 1) dt where p is the density, e is the energy per unit mass of fluid, k is the thermal conductivity, T is the temperature, and v is the fluid velocity. The bounding surface of the volume V is denoted by S, with unit normal Sj, and dr is a unit element of volume. Green's identity transforms surface integrals to corresponding volume integrals, and equation (Cl) becomes Wt Iv PedT ~ Iv (kT'j)'i d T + Iv ip£Vj)'i d T = ° Since this expression must hold for any volume V, we may discard the integrals and write jt{pe) - (kTtj). + (peVj). = 0 (C.2) Making use of the continuity equation equation (C.2) reduces to dp P ^ T - ( K T ^ J + P v ^ = ° (C-3) The energy of the fluid per unit mass may be divided into a kinetic part \v2 and the internal energy u, without including any chemical or other energy. The kinetic energy terms are cancelled by contributions from the neglected work terms in a full expansion of the heat equation, so we are left with only the internal energy part. The problem noted Appendix C. A Note on the Derivation of the Heat Conservation Equation 86 arises in the definition of this internal energy. By the first law of thermodynamics, du = dq — dw where dq is the differential quantity of heat absorbed by the system (in this case the volume of fluid), and dw = pdV is the work performed by the system. In keeping with previous simplifications, the differential quantity of work is discarded and thus du — dq — Tds For £ equal to either i or a spatial coordinate Xj, de du ds ds dT dZ = d t = d£= dT~d£ Substituting this relation into equation (C.3), and making use of the definition cp = T (ds/'dT) , conservation of heat energy is described by dT pep — - (kTA. + pcrVjTj = 0 (C.4) This derivation has made no assumption about the constancy of the heat capacity. Certain texts (e.g. Chandrasekhar,1961;Turcotte and Schubert, 1982) erroneously de-fine an absolute heat content per unit mass as u = cpT rather than fT u = c„dT JT, IT0 where T0 is a reference temperature which may be absolute zero. The first expression is true only in the case of constant heat capacity. Regardless of any assumptions about cp, the first definition leads to du = cpdT + Tdcp Appendix C. A Note on the Derivation of the Heat Conservation Equation 87 and cp is included in the time and spatial derivatives. Thus gradients of cp would exist in a full expansion of the heat equation. This is false and leads to dramatically different results for a convecting fluid with a moderate variation in the heat capacity. Numerical results The difference in fluid convection when an erroneous heat equation is employed is illus-trated by a series of numerical experiments. The reference model is the case of a constant property system at the critical Rayleigh number (figure 3.1 in chapter 3). Velocities do not evolve in time and the final and initial conditions are equal. Figure C l shows the time evolution of the same perturbation in the case of the correct and incorrect heat equations. All variables are those normalized and defined in chapter 3. The heat capac-ity varies by 120% around a mean value (figure 3.2c), and the time steps are one tenth those used in chapter 3. Velocities are relatively unchanged by a variable heat capacity (figure C.lb). However, when the internal energy is incorrectly defined, the derivative of the heat capacity enters directly into the heat equation, and a steep positive gradi-ent in the temperature dependence of c greatly enhances the convective vigour of the system (figures C.lc and C.ld). The opposite effect occurs with a negative temperature dependence for c — Raefj is less, and velocities decrease. Conclusion The assumption of a constant heat capacity is often valid in modelling fluid dynamic problems. In this case, the error in the heat equation discussed above does not affect results. When a variable heat capacity is modelled with the incorrect heat equation, how-ever, fluid velocities may be greater or less than they should be. A positive temperature dependence in the heat capacity leads to an increase in the effective Rayleigh number of the system, whereas a negative slope results in a reduced value of Raeff. The difference Appendix C. A Note on the Derivation of the Heat Conservation Equation 88 Figure C. 1: Effect of Erroneous Heat Equation on Thermal Convection (a) T(t=400) (b) Evolution of Perturbation 1| 1 1 0 2 0 400 x t (a) Temperature contours after 400 time steps with correct heat equa-tion, (b) Corresponding evolution of fluid velocities, (c) Temperature contours after 400 time steps with erroneous heat equation, (d) Cor-responding evolution of fluid velocities. (x,y,T, and t refer to their primed counterparts in the text.) between the correct and incorrect cases becomes more pronounced as the slope of dc/dT increases. If the heat capacity varies on the order of 100%, fluid velocities may be off by an order of magnitude. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items