Convection and lithospheric strength in Dione, an icy satellite of Saturn

ICARUS ~ ,

2 3 2 - 2 4 5 (1991)

Convection and Lithospheric Strength in Dione, an Icy Satellite of Saturn OLIVIER FORNI lnstitut d'Astrophysique Spatiale hat.120. Uniuersitt~ Paris-Sud, 91405 Orsay Cedex, France

ANGIOLETTA CORADINI I.A.S. Reparto di Planetologia tCNR), viale dell'Unicersitt'~ I1, 00185 Rome, Italy AND

COSTANZO FEDERICO Dipartimento di Scienze della Terra. Universita degli Studi, Piazza dell'Uniw, rsita, 06100 Perugia, Italy

Received January 3, 1991; revised June 5, 1991

ridges, and troughs (Consolmagno 1985). The faulting may Finite element numerical simulations for the structure of convec- be related either to external events, such as large impacts, tion cells within a self-gravitating, Boussinesq fluid sphere with or to internal stresses caused by the convective flow inside variable viscosity are used to determine temperature distribution, the satellite, and the overall contraction or extension due stream functions, velocities, and intensities of the stresses within to thermal evolution. Dione, an icy satellite of Saturn. A non-Newtonian, i.e., stressIn this p a p e r we discuss thermal convection inside an dependent, viscosity is considered, but the numerical simulation, icy body like Dione in order to assess the characteristics of according to U. R. Christensen (1984a, Geophys. J. R. Astron. the convective motions and the related stress distribution. Soc. 77, 343-384), is carried out knowing that the properties of Parameterized convection has been used to date to invespower law convection can be closely simulated by Newtonian flow tigate the thermal and structural evolution of icy m o o n s with a reduced value of the activation enthalpy. The initial temperature profile has been determined taking into account the amount in order to find a reasonable explanation for the o b s e r v e d of energy trapped as heat during the accumulation process and the tectonic features (Ellsworth and Schubert 1983, Federico and Lanciano 1983). The use o f p a r a m e t e r i z e d convection other energy source is the radioactivity of the long-lived isotopes. The obtained results seem to show a long duration endogenic was justified by the very limited amount of experimental activity in agreement with the geological observations. The intensi- data on rheology of the involved materials and by its ties of the stress are able to fracture the icy crust both in extensional simplicity. Recent d e v e l o p m e n t s in numerical studies of and in compressional regimes. The thickness of the mechanical convective systems (Christensen 1984a,b, 1985a,b, 1989) lithosphere evolves in a way determined by the complex interplay have drawn attention to the role played in thermal convecof the lithospheric temperature gradient and the strain rate. , t991 tion by a variable viscosity. T h e y have also shown that it Academic Pre~, Inc. is not possible to correctly describe this p h e n o m e n o n by the use of only two parameters, that is the Rayleigh (Ra) and Nusselt (Nu) numbers. M o r e o v e r , the effects of a 1. INTRODUCTION spherical g e o m e t r y and internal heating are not well unVoyager images of icy satellites of the outer Solar Sys- derstood in this modeling. In this p a p e r a finite element method is used to model tem have shown that on the surface of these bodies tecconvection in a spherical object of infinite Prandtl number, tonic and volcanic features are present. These features Boussinesq fluid having variable viscosity, and internal are generally believed to have originated through the comheating. We have studied the global flow structure and bination of several geologic processes (Murchie 1990). heat transfer in order to estimate the temporal evolution Three types of terrain have been recognized by the Voyof the mechanical lithosphere thickness. It is assumed that ager images: heavily cratered areas, s m o o t h e r plains, 232 0019-1035/91 $3.00 Copyright © 1991 by Academic Press, Inc. All rights of" reproduction in any form reserved.

CONVECTION AND LITHOSPHERIC STRENGTH

Dione is a homogeneous mixture of ice and silicate and that heat sources include gravitational energy r e l e a s e d . during satellite formation (accretional heat)-and radioactive decay of long-lived radionuclides (U, Th, K). In the numerical calculations of the thermal evolution some simplifying assumptions have been used and are discussed in Section 2. The most recent experimental results on the rheology of ice have been included in the model. Dione was selected as the object of our study, because its mass and density imply a reasonable amount of radioactive sources and because Dione appears to be a complex satellite from a geological point of view. In fact several geological units have been identified: heavily cratered terrains, plain units, lobate deposits, bright wispy material, and crater rim deposits (Piescia 1983). Within the cratered terrains, the crater floors appear either fiat or slightly convex depending on their diameter. This indicates that the floors have evolved as the result of viscous relaxation of the ice crust. Wide expanses of plains have been attributed to processes leading to liquid eruptions along fissures or to explosive events producing deposits on the surface (Plesica 1983). Tectonic features like pit chains, faults, rifts, troughs, and ridges have also been observed (Plescia 1983, Moore 1984). One of the hemispheres of Dione is covered with bright wispy markings, associated with extensive systems of ridges and valleys. On the northern hemisphere a large and complex valley system is present and the northern polar region is crossed by a 500-km graben-like structure. To conclude, on the satellite geological processes were still active after the cessation of the heavy bombardment. They lasted approximately 1 Gyr. (Plescia 1983) with predominant extensional tectonics. 2. MODEL AND ASSUMPTIONS The determination of the thermal evolution of a body requires the knowledge of different physical parameters. Certain assumptions and simplifications are also necessary. These are: I. It is assumed that the icy satellite was formed by accumulation from a circumplanetary swarm which results in a characteristic initial temperature profile inside the object. This profile has been obtained by Federico and Lanciano (1983) and Squyres et al. (1988). Thus the earliest source of energy is the percentage of heat irreversibly trapped in the satellite during the accumulation. The obtained initial temperatures are such that melting of water ices is generally not possible, but a peritectic temperature of 175 K of the NH3-H20 system is reached. This point is of interest because cosmochemistry seems to indicate that temperatures in the primitive solar nebula were probably low enough to condense ammonia-water mixtures in the region of Saturn and beyond. Ammonia hydrates may be significant constituents of icy satellites,

233

comprising as much as 20% of the total mass (Lewis 1972). However, the melting of even this maximum mass fraction does not necessarily imply a differentiation of the satellite with a possible formation of a silicatic core. Thus we can assume that Dione is a homogeneous mixture of ices and silicatic particles. 2. Heat sources used in the calculations are the long half-life radioactive elements with chondritic concentrations (U = 0.013 ppm, Th = 0.04 ppm, and K = 0.078%, Turcotte and Schubert 1982). 3. In Dione no higher density ice phases are present and it has been assumed that thermal expansion, thermal conductivity, and specific heat are not dependent on the temperature. 4. The most important, but difficult to estimate, parameter is the viscosity and its dependence on temperature and pressure. In our case the increase of pressure with depth is too small to have an influence on the viscosity. Rheological laws of pure water ice are becoming better understood as a result of recent experimental works. These experiments, however, deal with pure ice and we are interested in the behavior of mixtures, whose flow properties do not greatly change with respect to those of pure ices (Durham et al. 1989). Moreover laboratory strain rates are orders of magnitude different from the ones needed. Different viscosity models, based on theoretical and experimental results, have been critically and extensively reviewed. For more details refer to Schubert et al. (1986). It is important to notice that several mechanisms which lead to water ice deformation have been envisaged, and the dynamic viscosity/x is independent of effective stress if effective shear stress is less than 0. I MPa. During deformation by thermally activated creep under the above presented conditions, viscosity can be expressed as /.t = /ZmeA¢-Tm/7"" l )

(1)

where/z m is the viscosity at the melting temperature Tm and A is a constant related to the activation enthalpy of creep. Values for A between - 18 and - 26 are reported from laboratory work, and for pure ice it can be assumed that ,o,m 1.0 × 1013 Pa sec (Goodman et al. 1981). The activation enthalpy, which corresponds to the widely used value of A = - 2 5 (Reynolds and Cassen 1979), is H* = 57 kJ/mole and it results in a significant temperature dependence. Using Eq. (I) with A = - 25, surface viscosity values do not allow viscous relaxation of any topographic features in a reasonable time interval. Experimental data (Durham et al. 1983, Kirby et al. 1985) on the deformation of polycrystailine water ice at high pressure and low temperature have shown that at a confining pressure of 50 MPa and T < 195 K the flow law =

234

FORNI. C O R A D I N I . A N D F E D E R I C O

is non-Newtonian. The fitting of the experimental data gives a law of the type b = CCr"e (H'/R7)

(2)

where b is the strain rate, o- is the effective stress. C and n are material parameters, and R is the gas constant. In the range of t e m p e r a t u r e in which we are interested, these are

n

=

3.7 +- 0.4

C = 3.94 × 10 -4 MPa "see

tt* = 27.4 k J/mole. The viscosity in this case is obtained a s / z -- cr/(3b). In a satellite h o m o g e n e o u s l y mixed, the presence of silicate in the ice would increase the viscosity by an order of magnitude. Friedson and Stevenson (1983) have found that a 50% by volume ice rock mixture has a viscosity 6 - 1 6 times greater than that of pure ice, depending on the size distribution of the silicate particles. If a m m o n i a is also present in a satellite, the mixture is likely to be dominated by water ice, with N H 3 in the form of a hydrate at a concentration of only 20% by mass. As already mentioned the eutectic point of a m m o n i a hydrate is reached at a t e m p e r a t u r e of 175 K where partial melting also occurs. We have not included the effect of melting in our numerical simulation. Our objective was to avoid the complexities which arise from the strong nonlinear effects of the melting of the material, through the acquisition or the release of latent heat. Moreover, the rheology of a solid mixture of water and a m m o n i a ice is not sufficiently known, although recently Durham et al. (1988) carried out experimental studies to determine the rheoiogical properties of this mixture under pressure and t e m p e r a t u r e conditions that probably exist in the icy moons. The preliminary results also in this case show a nonlinear flow with a stress exponent on the order of 6 and an activation enthalpy of lI* = 65 kJ/mole. In particular, at T = 152 K and at a strain rate of 4 × 10 -5 sec ~, the flow stress is approximately 65 MPa. These results are consistent with the idea that a satellite with a m m o n i a is weaker than a pure water ice one. H o w e v e r , the effects of different concentrations of N H 3 as well as those of partial melting, which are of interest in our simulations, are still unknown. From a theoretical and experimental point of view the rheology of ice and of an a m m o n i a - i c e mixture seems to be non-Newtonian. M o r e o v e r two competing effects act on the viscosity: the presence of a m m o n i a , which m a k e s the mixture weaker, and the presence of silicatic grains. which makes it harder. In order to model the thermal evolution with stress-, temperature-, and pressure-dependent viscosity, we can follow Christensen (1984a), who

has shown that the properties of p o w e r law convection can be closely simulated by Newtonian flow with a reduced value of the activation enthalpy. As reported belore, water ice in the T - P range of interest has a nonNewtonian behavior with a stress exponent equal to 3.7 and an activation enthalpy of 27.4 k J/mole. According to Christensen (1984a), for a time evolution numerical simulation the activation enthalpy is to be reduced by a factor of 0.6. In this way a Newtonian analogue to the case of a non-Newtonian rheology can be constructed (Christensen and Yuen 1989). In our simulation it is possible to determine the preexponential factor of this N e w t o n ian analogue law, requiring that at T = 150 K the viscosity assumes the same value as the one given by the viscosity law (1). Thus we obtain /z = 2.42 x lOt6e ¢)/RI,

(3)

where Q = 16460 J/mole. The viscosity law in Eq. (3) hardens the material with respect to the viscosity value given by Eq. (!). Note that in Eq. (1) a value of A = 25 is used to take into account the presence of silicatic particles. Besides, in the subsurface layers the values of the viscosity are such that relaxation of crater basins in reasonable time scales is possible. An estimation of this relaxation interval for a crater with a diameter D on a body with density p and gravity g is given by -

r~ - 3.21/z~ff

(4)

pgD (Passey and S h o e m a k e r 1982). The effective viscosity takes into account the variation of the subsurface viscosity with depth in the lithosphere of the satellite. For D = 10, 50, and 100 km and for Dione at different times of its evolution: 0. I, I, and 2 Gyr., we obtain the results shown in Fig. 1. in order to obtain the effective viscosity from Eq. (3), a mean value for the t e m p e r a t u r e was used. This value c o r r e s p o n d s to the average of the t e m p e r a t u r e s to a depth of D in the lithosphere of Dione. In fact, conditions in a layer of size c o m p a r a b l e to the topographic feature considered must be taken into account (Thomas and Squyres 1988). In Fig. I the vertical bars represent r~ at different times (0. I, i, 2 Gyr.) during the history of the satellite for craters with diameters I0, 50, and I00 km. In determining the viscous time scale we have taken into account the time evolution of the satellite, using the corresponding thermal profiles. F r o m the inspection of Fig. I, it follows that for craters D > 50 km it is possible to have viscous relaxation of the floors on reasonable time scales, during the early evolution of the satellite. This result is in agreement with Passey (1983), who found values for the viscosity of between 1023 and 10 24 Pa sec at the top of the

CONVECTION

I0 e

,

,

.

.

.

.

.

.

.

.

.

.

.

.

,

.

.

.

.

AND

,

.

LITHOSPHERIC

, , , ,

10 /

lO s j/

/. j

//""

//'l

/

50 //

/"

:

/

10 ( .< ¢..) r.,i 10 a y-

10 ? V

r10 l 1 0 ~'z

..." /"

..! "

1023

VlSCOSl'rY ( P a

l O'-',t

lO~S

s~l

FIG. 1. Time scale for viscous relaxation of the crater floor on Dione. The upper, the middle, and the lower lines refer to craters with diameters equal to 10, 50, and 100 kin, respectively. The vertical bars represent the time scales at t = 0. I, 1, and 2 Gyr. after Dione formation for craters whose diameters range from 10 to 100 kin.

lithosphere by studying the forms of the flattened craters on Enceladus. Moreover, the rheologicai law in Eq. (3), following Quareni et al. (1985), can be called " h a r d " since the increase in the value of Q causes an increase in the viscosity at the surface of the satellite, consequently making the lithosphere harder. This rheology is also relevant for the study of convection beneath an immobile lithosphere, as in the case of planetary bodies in which socalled "one-plate tectonics" are operative.

tion number Di = a g D / % , which determines the influence of viscous dissipation and adiabatic compression on convection, is much less than unity (about 0.005), and adiabatic compression and viscous dissipation have negligible influence on convection. Moreover since the Prandtl number, Pr = v / k , is at least equal to 1023, the inertial forces are unimportant. Specifying the viscosity function appropriate to the problem, as we have done in the previous section, the solution of the governing equations (equation of continuity, Navier-Stokes equations with variable viscosity in spherical geometry, energy balance, or heat equation) furnishes the values of the dependent variables:radial velocity (ur), tangential velocity (Uo), components of the stress tensor (o-0), stream function (tO), temperature (T). We adopt a spherical coordinate system with origin in the center of the satellite that is axially symmetric about the 0 = 0 axis so that the longitudinal velocity (u6) is everywhere zero and all derivatives with respect to 4) are likewise zero. 3.1. I n i t i a l a n d B o u n d a r y C o n d i t i o n s

The calculation is carried out as a transient problem by retaining the time derivative term in the energy equation. Since the problem has been posed as a transient problem, initial conditions are required. We have chosen to start with a spherical domain of no flow and with an initi',d temperature distribution similar to that already obtained by Federico and Lanciano (1983) for Dione at the end of its accumulation around Saturn. The mechanical equations are solved under the following boundary conditions, due to symmetry of the problem: at the axis of symmetry: u0=0 OUr _ OUo _ 0 O0 Or

3. MATHEMATICALFORMULATION Utilizing the Boussinesq approximation, many investigators have determined the governing equations for the thermal convection (see, for example, Normand and Pomeau 1977, Oxburgh and Turcotte 1978), and they will not be reported here. We consider an infinite Prandtl number and Boussinesq fluid in a sphere; one diameter is an axis of symmetry for the system. The fluid is heated entirely from within with a rate H(t), which is the energy rate for unit volume. The local acceleration of gravity is directly proportional to the radial position, as is appropriate for a homogeneous fluid sphere of constant density p. The viscous dissipation term is neglected and it is assumed that the satellite interior is too viscous for inertia and the Coriolis force to be effective. In fact, the dissipa-

235

STRENGTH

and at the outer surface for shear stress-free conditions

I~ r --

l Ou r _ r O0

r

O( us/r) - o. Or

The boundary conditions at the outer surface imply that this surface is free and is a mechanical boundary of an isolated sphere in free space (Chandrasekhar 1961). The energy equation requires a thermal boundary condition at all boundaries. At the outer surface, temperature is constant and equal to To = 80. At the axis OT/O0 = 0 because of the symmetry. At the center of the sphere the heat flux must be zero.

236

FORNI. CORADINI. AND FEDERICO field, and then solve the energy equation with the velocity field obtained from the first integration. We have specialized the code incorporating the dependence of viscosity on t e m p e r a t u r e and pressure, the dependence of the gravitational acceleration on the radial distance in spherical s y m m e t r y for a constant density body, and the d e p e n d e n c e of the radioactive energy source on time. In order to follow the time evolution, different criteria can be used. In all the results presented here the time increment At always satisfies the C o u r a n t - L e v i - F r e d richs condition, based on the distance between two radial points and the m a x i m u m radial velocity. 5. RESULTS

FIG. 2. Mesh used in this study, consisting of a finite element grid of 1350 elements, with 3(1 radial points and 90 tangential points.

4. NUMERICAL APPROACH We use the method of finite elements to obtain numerical solutions of the governing equations. For the sake of brevity we will not go into the details of the finite element method. The reader may refer to Zienkiewicz (1977). The numerical code we use is M A N T L E and it is vastly referenced in the literature ( T h o m p s o n and H a q u e 1973. T h o m p s o n 1975, Dawson and T h o m p s o n 1978) and is especially written to solve coupled equations for creeping viscous flow and heat transfer (Schubert and Anderson 1985). The equations describing the flow are expressed in variational form and are written directly as a statement of virtual work with the constraint of incompressibility directly incorporated. The spherical mesh used has 1350 elements and is illustrated in Fig. 2. The mesh has N r = 30 radial points and N~ = 90 tangential points. This c o r r e s p o n d s to a mean radial distance between two neighboring points of 18.5 km. Machetel and Yuen (1987) havc shown that for a mean radial distance of less than 36 km, in the case of the spherical convecting Earth, the numerical solution of the coupled thermomechanicai equations is independent from the mesh grid. being N , = 3 N r , as in our simulation. In each element, quadratic interpolation of the velocities and linear interpolation of the pressure to assure incompressibility are used. An efficient method for use of this code, from the computational point of view, is to solve first the flow equation with a given t e m p e r a t u r e

We have studied the evolution of the satellite Dione, whose reference p a r a m e t e r s are given in Table 1. The model we run has the Newtonian analogue viscosity given by Eq. (3). As the temporal evolution of the system is lbllowed, the stream lines at six successive times are reported in Fig. 3. Although the actual orientation of the convection within the satellite is not known and the equator of the model need have no relation to the geographic equator of the satellite, we o b s e r v e a high degree of symmetry in the stream lines a b o v e and below the equator. This pattern is obviously e x p e c t e d and it is obtained because of the a c c u r a c y of the numerical solution. A closer look at a single pattern of circulation reveals an a s y m m e try, which depends on the internal heating. Since velocity in the direction of the stream lines is inversely proportional to the distance between adjacent streamlines, the closer spacing indicates a higher velocity. Thus, the close spacing of the descending streamlines at the equator and the wide spacing of the upwelling streamlines indicate a relatively rapid narrow descending flow and a broad region of slow upwelling. Although these results are not new (see, for example, Jarvis and Peltier, 1989) they seem to confirm the accuracy of the numerical solution. In Fig. 4 the fundamental p a r a m e t e r in the theory of convection, the Rayleigh n u m b e r as function of time, is presented. The Rayleigh n u m b e r for variable viscosity

TABLE I Physical Parameter Values Utilized in the Simulation

Parameter Radius Density Expansivity Thermal conductivity Specific heat generation rate at t = 0 Mean decay constant

Value 5.6 × 105 m 1430 kg/m~ 10 4 K-I 3.2 W m : K 2 × 10 s W m 1.5 x 10-~ sec "

CONVECTION AND LITHOSPHERIC STRENGTH

237

a

I I I I l l l l l ~ l

I I I I I I l l l l l l l l l

Illll

I IIII

/i

I I L I I I I I I I l ] l l l l l l l l l l l l l l l

I llllll

II

I l l l l t l l l

II

11

I II

II

II

I111

I II

I111111

t~

I I I I I I l l l l l l l l l l l l l l l l l l l l

FIG. 3. Two-dimensional stream function fields for a time-dependent convection at six subsequent times. Contours of stream function are plotted with a constant interval of 0.04 at t = 100 Myr, of 0.4 at t = 500 and 9(10 Myr, of 0.2 at t = 1.3 Gyr., of 0.2 at t = 1.9 Gyr., and of 0.02 at t = 3. I Gyr. The arrows on the streamlines indicate the direction of the flow. This direction has been determined by knowing the velocity vectors at each grid point. The outer solid line represents the surface of the satellite.

238

FORN1.

CORADINI,

AND

FEDERICO

b

......

o~I

"\.

P.

~ .

.............

" J / .....

I 1_..±./

~

/----.-

//././~-- ...... ~.~',."\til!/--'0 ~.'.\\,,

I/~:

I

-

" i 1_1. I I I li_J

J ..1

._J

o

\

..." ........

! " .....

I ~ I i I I . . .L_I._;,

I I . I I 1.,

m.

L_

J 1 ! I I : 1 ±_L.I

.

. ;~_ I 1 1~'..,

I l~_i

I I I 1 ~ I I 1 LS._II Flti.

3--('ontim,'d

.oh'

\

./

1

I I I ~ I I J..

.

.

.

1 1._1. I , I I I i LL

.

: I I t I '. ] I I 1 1 I I I I_.

]

i 1 1 I

CONVECTION AND LITHOSPHERIC STRENGTH

°,I I J

104~

i /

;

"~ .\\

!

/

10 3

J

I0: 500

I000

1500

2000

" 2500

3000

3500

4000

4500

5000

T I M E [Ma)

FIG. 4. Rayleigh number as function of time during the evolution of the satellite. The convective vigor first increases, during the delivering of the accumulation energy, and then slowly decreases. The critical Rayleigh number is assumed to be 2000.

convection has been defined following Christensen (1985a). The " i n t e r n a l " Rayleigh number was chosen by using a viscosity calculated for a globally averaged temperature. The mean value of the temperature, /'mean, is obtained by averaging over the modulus of the velocity

Vu~ + u0. In spherical symmetry and for an internally heated system we have Ra-

4rrGapa6H 3 Kktx( Tmean)'

(5)

where a , / 9 , k, and K are thermal expansion coefficient, density, thermal expansivity, and conductivity, respectively, and a is the radius of the satellite. The vigor of convection increases during the early phases because there is a nonnegligible energy contribution by the initial thermal profile. In fact, during the first 500 Myr the energy delivered by the radioactive elements is 60% of the energy associated with the initial thermal profile. Successively the convection is sustained by the radioactive energy and then slowly decreases, reaching values lower than the critical one at t = 3300 Myr. Following Busse (1989), we assume that 2000 is the critical value for the Rayleigh number. In Figs. 5a-5c radial variations of the horizontally averaged temperature profiles are reported. The general behavior of these profiles is analogous to the corresponding profiles published by Peltier et al. (1989). In making this comparison we must keep in mind that our Ra is on the order of 103-104 (Fig. 4) and our convective system is heated only from within. As time goes on, the temperature rises and a larger quantity of the internal material reaches

239

a temperature higher than the eutectic temperature of the N H 3 - H 2 0 system. As it is obtained from experimental results of Croft et al. (1988), the density of the ammon i a - w a t e r peritectic liquid at 175 K is 946 kg/m 3. In an undifferentiated icy satellite with mean density of 1430 kg/m 3 this liquid should have no problem rapidly reaching the surface. We can try to estimate the needed time using a firstorder result of the theory of magma migration through a porous medium, as given in Turcotte and Schubert (1982). Assuming that the pressure gradient drives the magma upward, we can obtain the relative velocity between the magma (liquid ammonia) and the matrix (solid ice-silicate mixture). Using a magma viscosity of I Pa sec (Kargel 1988), a density contrast of At) = 484 kg/m 3, a grain size of the order of 10 -3 m, and a mass fraction of 20%, a relative velocity of 5 × 10 7 m/sec is obtained. This velocity corresponds to a time of 10 4 years to move a parcel of N H 3 liquid magma from the center of the satellite to its surface. Thus it can be assumed that this process does not significantly alter the mechanics of overall convection. The temporal behavior of the surface heat flow (Fig. 6a) shows a rapid early decrease, due to the initial high thermal lithospheric gradient connected with the initial profile. A more gradual decrease throughout most of the lifetime of the satellite follows the first phase. A slight increment is observed at t = 800 Myr, when the internal temperature reaches its maximum. The effect of the initial temperature profile is made clearer if the two curves, shown in Fig. 6a are compared. In fact the higher curve is simply the convective cooling of a satellite whose initial heat flux and other characteristics are equal to that of Dione. In the presence of decaying heat sources, it is important to compare the behavior of the heat flowing outward and the heat produced within the convective system. For this purpose the Urey ratio can be defined, which in our notation is given by

Ur - H(t) q,(t)'

(6)

where q~(t) is the averaged heat flux at the surface of the satellite. We find (Fig. 6b) that the initial slope of the Ur(t) depends on the initial thermal condition, and being that Ur > 1, the satellite is heated up. Then Ur decreases toward the equilibrium value I. The turning point of this curve is at t = 500 Myr, while at t = 2500 M y r a value comparable to the mean life of the radioactive source, a further slight increase of Ur is found. At this point it is useful to introduce the Nusselt number in order to obtain the N u - R a relationship. Following Christensen (1985a), we define Nu as

240

FORNI, CORADINI, AND F E D E R I C O

a

:

TIME

~00

r

180

•

0

-

20

-

,50 - 100 Ma

b

•

200 f . .

160 L

•

140

~

.

.

--

-~-. g---i._I

180 !..

.

v '

: 200 - 1000 Ma

TIME

220,"

160F ~---. •

"

,..... ...

-

140,

,

120,

',iil~,,

lO0 i 80

10o, .

G

.

.

.

.

.

100

200

300 RADIUS

soL__

600

500

400

.

.

O

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

200

1 O0

(Ki'l~)

.

.

.

.

.

300 RADIUS

.

.

.

.

.

.

4(1(I

.

.

.

500

6(1(I

(Ella)

TIME : 1.1 - 4.7 Gs ~:~0

11.o ..... "

i

.......

I~O,

1401

""

120"

lO0,

!

00

. . . . . . . . . 0

FIG. 5.

1

O0

:200

30(1

4 O0

500

t100

Averaged temperature profiles as a function of the radius at different times during all the history of the satellite.

Nu-

"/'~ond

(7)

T~ .... '

where

3q41

7"~,,°d = -~-

+ "/~,.

Figure 6c shows the time behavior of Nu. After the initial phases, Nu is larger than i between 500 and 2000 Myr. For this time interval it is possible to compute the coefficient and the power law exponent that define the N u - R a relationship. As is shown in Fig. 7, a coefficient equal to 0.06 and an exponent equal to 0.31 are obtained. Thc viscosity contrast between the upper boundary and the center of the satellite is constant and equal to 6 x I0 (', as can be argued by looking at the temperature profiles,

reported in Fig. 5. According to Christensen (1985a, 1989), for variable viscosity convection at least two parameters are needed to correctly describe the convective system. In our case, one of these parameters, the viscosity contrast, remains constant, and the value of the exponent seems to be in agreement with the values obtained by Quareni and Yuen (19841 for a constant viscosity fluid, heated from within and with a stress-free boundary condition, and by Quareni et al. (1985) for a temperature dependentviscosity fluid, heated at the base and with a stress-free boundary condition. We now focus our attention on the successive thermal phases encountered by Dione during its evolution. (i) The first 100 Myr (Figs. 5a, 8a, and 9a): During the first 100 Myr after the formation of the satellite, the evolution of the temperature profile (Fig. 5a) is similar to that obtained by Squyres et al. (1988). The pcak temperature decreases and moves deeper into the interior, while the temperature in the center of the satellite rises. However. the stress be-

241

CONVECTION AND LITHOSPHERIC STRENGTH

b

4.5

3.8 36

/ \

3.4

3.5

3.2. m

2.5

~

-

~ ' ~

/

\

3

""

]

2.8

~

2.6

1.5 2.4 1

2.2

0.5 0 -

2 -

500

tO00

2000

1500

2500

3000

3500

4000

4500

I

I.A

5000

500

1000

1500

2000

TIME (Mi)

C

2500

3000

3r)O0

4000

4500

5000

TIME (Nit)

2.5

2

1.5 ~

I ....................

- ~ ~

~........ ~ - 5 -

~ ~. . . .

m 5~ ~

i

0.5L

oi

1

0

500

1000

1500

I

2000" 2500

. . . . . . .

3000

3500

4000

4500

5000

TIME (Ma)

FIG. 6. Surface heat flux (a), Urey ratio (b), and Nusselt number (c) as a function of time during the lifetime of Dione. The upper curve in (a) represents a simple model of convective cooling neglecting the effects due to the initial temperature profile.

0.12

. . . . . . .

/i

!

////

0.1

J

0.08

3/ / /

"~i 0.06 3

J

.(/

004 ~

{

/ 0.02

////

0,'" 395

~ 4

405

. 4.1

.

.

.

.

4 15

.

. 42

.

.

.

.

4,25

I

. 4,3

4.:)5

I.og(Rll)

FIG. 7. N u - R a relationship for a developed convection in our simulation of the history of Dione. in this case the exponential factor is equal to 0.3 I.

havior is very different because solid state convection sets in early (Fig. 3). Thus, convective motions of material are mainly responsible for the induced sublithospheric stresses. In Fig. 8a, the state of convective stress in the subsurface of the satellite at a depth of about 6 km is shown as a function of the latitude, and at t = 100 Myr. This state of stress is given by the vertical stress (o'~), the horizontal stress (or00), and the shear stress (o'r0), and we consider compressive stresses negative and tensile stresses positive. o-~ is not reported here because in our computations it is always the intermediate principal stress. The 6 km depth can be interpreted as a representative thickness of an elastic lithosphere, as will be shown later. Under the action of this stress field that the convective mantle imposes on the overlying shell, the lithosphere itself deforms. We have not computed the elastic displacements because we are interested only in determining where the stress field is compressive or extensive. In fact this pattern will be used to define the lithospheric strength envelopes (Golombek and

242

FORN[. CORAI)INI. TIME:IO0 t;0 .

AND FEDERICO

b

Ma ...,_

80

TIME : 900 ................ •.........................

M~

t30

V 10"

40

A c

;~O,-

o.

"~ 0 i.

-F

20

~

t

o s

~

- ;30,

-~0 ~4(1

-40

t~c

..-

.60

60 -HO HI)

I O0.

-IOO; HO

60

--10

-L~0

0

40

:_~0

t;O

H{I

-HO

-60

40

0

20

4H

BO

HO

I~TITUI)E

I.ATI I ' I l I ) E

C

-20

TIME

3.9

(;.

40

:2O • l(}

'

O"

"~

"10' -;20 • - ;|(I •

-4(},

-60

-40

- 20

(l

20

40

~(}

HO

LATITI'D[e~

F I G . 8. Subsurface vertical, horizontal, and shear stress as a function o f the latitude and at different times at a d e p t h oi"6 k m . N o t e that at the equator the shear stress is roughly zero so that the other [w,o c o m p o n e n t s of the stress can be considered as principal.

Banerdt 1986). At 100 Myr (Fig. 8a) and at the equator o% is roughly zero and the remaining components of the stress tensor can be interpreted as principal stresses. In the equatorial region (rrr > o'00and the state of stress is compressive. At midlatitude -+ 40 ° the state of stress is tensional with (rr, < (r~o. From the inspection of Fig. 8a, a "'quasi"-compressional state of stress is also located at ±70 ° of latitude. It is to be noted that this pattern remains unchanged during the evolution, and there are only changes in the intensities of the deviatoric stresses, as it results from the inspection of Figs. 8b and 8c. (ii) The heating phase (200 Myr-I Gyr.) (Figs. 5b, 8b. and 9b): The cooling from the surface competes with the warming of the interior by radiogenic heating. In the interior of the satellite the temperature increases and reaches a maximum, because the heat produced by the decaying radionuclides overcomes the heat removed by convection. This maximum temperature of about 210 K is reached around 800 Myr. The motion of the convective

cells is relatively rapid and the values of the stream function and of the Rayleigh number at this time reach their maximum, i.e.. 2.5 (Fig. 3a) and 12,500, respectively (Fig. 4). During this period, the two convective cells that were located near the surface at the beginning of the evolution sink into the interior. Their final location is reached after 500 Myr as is shown by the stream function (Fig. 3a) and by the averaged temperature profiles (Fig. 5b). (iii) The cooling phase (I. I-4.7 Gyr.) (Figs. 5c, 8c, and 9c): During this period there is cooling and a decrease in the convective activity: in tact, Ra goes from 10,0(R) to less than 2000 (Fig. 4) and the maximum value of the stream function decreases from 2.3 to less than 1 (Fig. 3b). Moreover, as results from the inspection of Fig. 5c the cooling approaches a conductive behavior. In order to obtain a more representative view of the rheological stratification of the mechanical lithosphere of this icy satellite, we have studied the transition between brittle and ductile regimes. The importance of this transi-

243

CONVECTION AND LITHOSPHERIC STRENGTH R

TIME : I 0 0 M a

b

/" /

0

-5

!i!, _,5~

/

-I0 i

I -I"

-I,5

_ ot 1,5~"-.~ ....

,

f

I

-25 -30 -35 -40 - -50

I

i

.J

-20

TIME : 9 0 0 Ma

-40

-30

-20

-I0

,

-~SL

"

-30 I

I

- 35 [

.-J

0

I0

20

•.

_40 I_. . . . . . . . . . .

30

-50

-30

-40

-20

STRESS DIFFERENCE (MPL)

TIME

o

10

~0

:~0

3 . 9 Ga

7- - - -

,/

-I,5

d

"~-.

-2oI

\

30[

i

-2s[

i

J

.:

-3,5 -40 -50

- 10

S T R E S S DIFFERENCE (Mf'a~

or

v

. . .

:.. - -

\

'\

/

i

'\

- -

-40

-30

-20 ~RESS

-I0

0

10

20

30

DIFFERENCE (MPn}

FIG. 9. Frictional (linear branches) and viscous (exponential branches) yield stress versus depth for compression (to the left) and extension (to the right) at three successive times during the life of Dione.

tion in the development of extensional and tensional instabilities in the lithosphere of icy satellites has been recently shown by Herrick and Stevenson (1990). The depth of this transition separates an upper zone where material deforms by brittle failure from a more internal zone where deformation occurs under ductile conditions. The value of the stress at the brittle-ductile transition gives an estimate of the strength of the lithosphere. The frictional stress difference is obtained by applying rock mechanical models. The maximum stress difference supported by friction can be computed for two limiting cases, horizontal compression and horizontal extension, and it is given in Beeman et al. (1988); we consider compressive stresses negative and tensile stresses positive. The ductile stress difference is a result of our numerical simulation and it also depends on the depth; we have taken into account only the maximum extensional and compressional values. From our simulation, as we have already shown, it results

that the maximum compressional stress difference is located at the equator, while the maximum extensional stress is located at a latitude of about 40 °. The "strength" of the lithosphere at any given depth is the lower of the brittle and the ductile stress differences. In Figs 9a-9c the stress difference is reported as a function of depth and at three times. The strength under tension is found to be smaller than or equal to 10 MPa, in agreement with Golombek and Banerdt (1987). This value is greater than that required to initiate frictional sliding on preexisting fractures and faults (2 MPa). This suggests that tensile fracture can constitute a mode of fracture throughout the "active" lifetime of this satellite. The maximum tensile thermal stress obtained by Hillier and Squyres (1990) for a viscoelastic satellite in the case of Dione is about 6.6 MPa. This value is roughly equal to the one we obtained at t = 100 Myr in the extensive region (Fig. 9a). In fact, the deviatoric stress of such a value

FORNI, CORADINI, AND FEDERICO

244 IIRITI'I.E/I)UCTII.E

TI~ANNII'IDN

~'~

I'IMI'~

I I

Li ItJ k t,,

g

(; t;

:~(: K t~L

tl f ,

,J I

/

o :_"

11 I j

:

-

II)" IM

tt)

J

FIG. 10. Thickness of the mechanical lithosphere of Dione as functions of strain rate and temperature gradient. Note that at the beginning of its evolution, Dione is characterized by the higher value of the lithospheric thermal gradient.

in a brittle shell is able to reactivate preexisting faults (Beeman et al. 1988). Our modeling, in any case, is more accurate because it is capable of predicting latitudinal and depth variations in the stress field behavior. Brittle compressional deformation probably occurs by shear failure; the shear strength of unfractured ices is about 10 M Pa and this value is also reached by the strength in compressional regime. In both conditions the thickness of the brittle lithosphere varies in time, but in the extensional regime there is a larger variation. In fact, in the early phases the mean thickness is about 15 km. Thickening continues until a quasi-constant value of 28 km is reached, while in later stages, due to the complex interplay among damping of the convection, lithospheric thermal gradient, and viscosity, there is a further thinning of the brittle lithosphere. This process is represented in Fig. 10 where the thermal gradient as a function of the strain rate for different values of the brittle lithospheric thickness is shown. In this plot the trajectory, described by Dione, goes from relatively high thermal gradients to lower ones. It is clearly shown that when the vigor of convection decreases there is a slight thinning of the lithosphere. Although uncertainties in the rheoiogical parameters exist, as previously discussed, we consider that the main features of these results seem to be established with a reasonable degree of confidence.

profile that allows the early onset of convection. The resuiting stress distribution reflects both the presence of an early convection and the peculiar rheology chosen. Nevertheless, in this early stage, it is assumed that the stresses are able to fracture the primordial crust of the satellite. Moreover, the magnitude of these stresses is higher than those induced only by the thermal cooling and recently studied by Squyres et al. (1988) and Hillier and Squyres (1990). In the subsequent history of the satellite, the stresses increase during further heating and the contemporary sinking of the convective cells down to the interior of the satellite. It is during this period that the rise of liquid mixtures of N H 3 - H 2 0 seems to be possible. This is in agreement with the presence of extensional features, recognized by Plescia (1983) and Moore (1984). These structures are formed during the first evolutionary phases of the life of Dione together with a large amount of resurfacing. This period of internal heating produces an increase in the convective activity that reaches its maximum after about 800 Myr. In the same time interval the stresses reach their maxima, and the values are large enough to allow fracturing both in extensional and in compressionai regimes. This "high stress period" lasts for more or less than 3 Gyr., notwithstanding the slow satellite cooling and the sinking of the convective cells. This behavior is the consequence of two competing effects: (I) the increase in the viscosity as the temperature decreases, and (2) the decrease in the strain rate due to the slow turning off of the convective activity. During a long period of time these two competing effects balance each other, maintaining a more or less constant stress level large enough to reactivate already formed fractures. This could be a possible explanation for the relatively long tectonic activity recognized on a small satellite like Dione. It must be pointed out that the stress regime that induces this tectonic activity in our model is due to the convective motions. The stress induced by overall contraction or extension due to thermal and tidal effects and by differentiation should be superimposed to obtain a complete "stress history." In conclusion, our modeling is an approximation of the real non-Newtonian viscosity behavior. H o w e v e r , the obtained results seem to indicate that it would be interesting to simulate the satellite evolution using nonlinear rheology like that proposed for water ice by Durham et al. (1983) and for a w a t e r - a m m o n i a ice mixture by Durham et al. (1988).

6. CONCLUDING REMARKS

ACKNOWLEDGMENTS

Based on the above presented results, a tentative history for the satellite Dione can be suggested. During the first hundred millions years, the thermal history is controlled by the effects of the primordial accretion heating and consequently, by the presence of a nonuniform initial thermal

This work was partially supported by MURST (C. Federico) and by INSU-ATP Plan~tologie (O. Forni). The authors acknowledge the computing center of Rome, Universit~ "La Sapienza," and the computing center (CIRCE) of the CNRS. The authors thank Dr. J. V. A. Keller for his constructive comments and the two reviewers for their helpful remarks.

CONVECTION AND LITHOSPHER1C STRENGTH REFERENCES BEEMAN, M., W. B. DURHAM, AND S. H. KIRBY 1988. Friction of ice J. Geophys. Res. 93, 7625-7633. BUSSE, F. H. 1989. Fundamentals of thermal convection. In Mantle Convection: Plate Tectonics and Global Dynamics (W. R. Peltier, Ed.), pp. 23-95. Gordon & Breach, New York.

CHANDRASEKHAR,S. 1961. Hydrodynamic and Hydromagnetic Stability. Oxford Univ. Press, New York.

CHRISTENSEN,U. R. 1984a. Convection with pressure- and temperaturedependent non-Newtonian theology. Geophys. J. R. Astron. Soc. 77, 343-384. CHRISTErqSEN, U. R. 1984b. Heat transfer by variable viscosity convection and implications for the Earth's thermal evolution. Phys. Earth Planet. Inter. 35, 264-282. CHRISTENSEN, U. R. 1985a. Heat transfer by variable viscosity convection. 11. Pressure influence, non-Newtonian rheology and decaying heat sources. Phys. Earth Planet. Inter. 37, 183-205. CHRISTENSEN, U. R. 1985b. Thermal evolution models for the Earth. J. Geophys. Res. 90, 2995-3007. CHRISTENSEN, U. R. 1989. Mantle rheology, constitution, and convection. In Mantle Convection: Tectonics and Global Dynamics (W. R. Peltier, Ed.), pp. 594-655. Gordon & Breach, New York. CHRISTENSEN, U. R., AND D. A. YUEN 1989. Time-dependent convection with non-Newtonian viscosity. J. Geophys. Res. 94, 814-820. CONSOLMAGNO, G. J. 1985. Resuffacing Saturn's satellites: Models of partial differentiation and expansion. Icarus 64, 401-413. CROFT, S. K., J. I. LUNINE, AND J. KARGEL 1988. Equation of state of ammonia water liquid: Derivation and planetological applications. Icarus 73, 279-293.

245

KARGEL, J. 1988. Liquidus phase relations and liquid properties in the system H20-NH3-CO2-H2CO. Lunar Planet. Sci. XIX, 583-584. KIRBY, S. H., W. D. DURHAM, AND H. C. HEARD 1985. Rheologies of H20 ices l, II and Ill at high pressures: A progress report. In Ices in the Solar System (J. Klinger and D. BeHest, Eds.), pp. 89-I07. Reidel, Dordrecht. LEwis, J. S. 1972. Low temperature condensation from the solar nebula. Icarus 16, 241-252. MACHETEL, P., AND D. A. YUEN 1987. Chaotic axisymmetrical spherical convection and large-scale mantle circulation. Earth Planet. Sci. Lett. 86, 93-104. MOORE, J. M. 1984. The tectonic and volcanic history of Dione. Icarus 59, 205-220. MURCHIE, S. L. 1990. Tectonics of icy satellites. Adv. Space Res. 10, 173-182. NORMAND, C., AND Y. POMEAU 1977. Convective instability: A physicist's approach. Rev. Mod. Phys. 49, 581-624. OXBURGH, E. R., AND D. L. TURCOTTE 1978. Mechanism of continental drift. Rep. Prog. Phys. 41, 1249-1312. PASSEY, Q. R. 1983. Viscosity of the lithosphere of Enceladus. Icarus 53, 105-120. PASSEY, Q. R., AND E. M. SHOEMAKER 1982. Craters and basins on Ganymede and Europa: Morphological indicators of crustal evolution. In Satellites o f Jupiter (D. Morrison, Ed.), pp. 379-434. Univ. of Arizonia Press, Tucson. PELTIER, W. R., G. T. JARVIS, A. M. FORTE, AND L. P. SOLHEIM 1989. The radial structure of the mantle general circulation. In Mantle Convection: Plate Tectonics and Global Dynamics (W. R. Peltier, Ed.), pp. 765-815. Gordon & Breach, New York. PLESCIA, J. B. 1983. The geology of Dione. h'arus 56, 255-277.

DAWSON, P. R., AND E. G. THOMPSON 1978. Finite element analysis of steady state elasto-visco-plastic flow by the initial stress rate method. Int. J. Num. Methods Eng. g, 3-16.

QUARENI, F., AND D. A. YUEN 1984. Time-dependent solutions of mean field equations with applications for mantle convection. Phys. Earth Planet. Inter. 36, 337-353.

DURHAM, W. B., H. C. HEARD, AND S. H. KIRBY 1983. Experimental deformation of polycrystalline H20 ice at high pressure and low temperature: Preliminary results. J. Geophys. Res. 88, 377-392.

QUARENI, F., D. A. YUEN, J. SEWELL, AND O. R. CHRISTENSEN 1985. High Rayleigh number convection with strongly variable viscosity: A comparison between mean-field and two-dimensional solutions. J. Geophys. Res. 90, 12,633-12,644.

DURHAM, W. B., S. H. KIRBY, AND L. A. STERN 1988. Rheology of the water-ammonia ice: First results. Lunar Planet. Sci. XIX, 285-286. ELLSWORTH, K., AND G. SCHUBERT 1983. Saturn's icy satellites: Thermal and structural models. Icarus 54, 490-510. FEDERICO, C., AND P. LANCIANO1983. Thermal and structural evolution of four satellites of Saturn. Annal. Geophy. 1,469-476. FRIEDSON, A. J., AND D. J. STEVENSON 1983. Viscosity of rick-ice mixture and application to the evolution of icy satellites. Icarus 56, 1-14. GOLOMBEK, M. P., AND W. B. BANERDT 1986. Early thermal profiles and lithosphere strength of Ganymede from extensional tectonic features. Icarus 68, 252-265. GOLOMBEK, M. P., AND W. B. BANERDT 1987. Failure strength of icy lithospheres. Lunar Planet. Sci. XVIII, 337-338. GOODMAN, D. J., H. J. FROST, AND M. F. ASHBY 1981. The plasticity of polycristalline ice. Philos. Mag. A 43, 665-695. HERRICK, D. L., AND D. J. STEVENSON 1990. Extensional and compressional instabilities in icy satellite lithosphere. Icarus 85, 191-204. HILLIER, J. K., AND S. W. SQUYRES 1990. Differential thermal stresses: A source for early extensional tectonism on the Saturnian and Ur/mian satellites. Lunar Planet. Sci. XXI, 514-515. JARVIS, G. T., AND W. R. PELTIER 1989. Convection models and geophysical observations. In Mantle Convection: Plate Tectonics and Global Dynamics (W. R. Peltier, Ed.), pp. 479-593. Gordon & Breach, New York.

REYNOLDS, R. T., AND P. M. CASSEN 1979. On the internal structure of t he major satellite of the outer planets. Geophys. Res. Lett. 6, ! 21 - 129. SCHUBERT, G., AND C. A. ANDERSON 1985. Finite element calculation of very high Rayleigh number thermal convection. Geophys. J. R. Astron. Soc. 80, 289-31 I.

SCHUBERT,G., T. SPOHN, AND R. T. REYNOLDS 1986. Thermal histories, compositions, and internal structures of the moons of the Solar System. In Satellites (J. A. Burns and M. S. Matthews, Eds.), pp. 224-292. Univ. of Arizona Press, Tucson. SQUYRES, S. W., R. T. REYNOLDS, A. L. SUMMERS, AND F. SHUNG 1988. Accretional heating of the satellites of Saturn and Uranus. J. Geophys. Res. 93, 8779-8794. THOMAS, P. J., AND S. W. SQUYRES 1988. Relaxation of impact basins on icy satellites. J. Geophys. Res. 93, 14,919-14,932. THOMPSON, E. G. 1975. Average and complete incompressibility in the finite element method. Int. J. Num. Methods Eng. 9, 925-932. THOMPSON, E. G., AND M. I. HAQUE 1973. A high order finite element for completely incompressible creeping flow. Int. J. Num. Methods Eng. 6, 315-321. TURCOTTE, D. L., AND G. SCHUBERT 1982. Geodynamics. Wiley, New York. ZIENKIEWICZ, O. C. 1977. The Finite Element Method. McGraw-Hill, London.