A model for the interior structure, evolution, and differentiation of Callisto

Icarus 169 (2004) 402–412 www.elsevier.com/locate/icarus

A model for the interior structure, evolution, and differentiation of Callisto K. Nagel, D. Breuer, and T. Spohn ∗ Institut für Planetologie, Westfälische Wilhelms Universität, Wilhelm-Klemm-Str.10, D-48149 Münster, Germany Received 20 December 2002; revised 25 September 2003 Available online 6 March 2004

Abstract The recently measured dimensionless moment of inertia (MoI) factor for Callisto of 0.3549 ± 0.0042 (Anderson et al., 2001, Icarus, 153, 157–161) poses a problem: its value cannot be explained by a model in which Callisto is completely differentiated into an ice shell above a rock shell and an iron core such as its neighboring satellite Ganymede nor can it be explained by a model of a homogeneous, undifferentiated ice–rock satellite. We show that Callisto may be incompletely differentiated into an outer ice–rock shell in which the volumetric rock concentration is close to the primordial one at the surface and decreases approximately linearly with depth, an ice mantle mostly depleted of rock, and an about 1800 km rock–ice core in which the rock concentration is close to the close-packing limit. The ice–rock shell thickness depends on uncertain rheology parameters and the heat flow and can be roughly 50 to 150 km thick. We show that if Callisto accreted from a mix of metal bearing rock and ice and if the average size of the rocks was of the order of meters to tens of meters, then Callisto may have experienced a gradual, but still incomplete unmixing of the two components. An ocean in Callisto at a depth of 100–200 km is difficult to obtain if the ice is pure H2 O and if the ice–rock lithosphere is 100 km or more thick; a water ocean is more plausible for ice contaminated by ammonia, methane or salts; or for pure H2 O at a depth of 400–600 km.  2004 Elsevier Inc. All rights reserved. Keywords: Callisto; Interiors; Satellites of Jupiter; Thermal histories

1. Introduction The gravity data returned by the Galileo mission from recent Callisto fly-bys allowed improved measurements of the mass and, for the first time, accurate measurements of the coefficient C22 of the gravity field (Anderson et al., 2001). From these, the dimensionless moment of inertia about the rotation axis C/MS RS2 , where MS is the satellite mass and RS its radius, was estimated to be 0.3549±0.0042. C/MS RS2 can be calculated from C22 assuming a satellite in hydrostatic equilibrium and is usually taken to be approximately equal to the MoI-factor defined as the dimensionless average moment of inertia
1 I = ρ(r)r 4 sin3 θ dr dθ dϕ, (1) MS RS2 MS RS2 where ρ(r) is the density, r is the radial distance from the center and θ and ϕ are the co-latitude and longitude, respec* Corresponding author.

E-mail address: [email protected] (T. Spohn). 0019-1035/$ – see front matter  2004 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2003.12.019

tively. The MoI-factor is an important quantity to constrain models of the interior structure of a planet or satellite because it measures the first moment of the interior density distribution. In the following, we will assume that C/MS RS2 equals the MoI-factor. For a homogeneous satellite the MoI factor is 0.4. If the density increase due to self-compression under the action of gravity and due to ice phase transitions is taken into account, the MoI-factor for a homogeneous Callisto can be as small as 0.38 (McKinnon, 1997). The observationally constrained MoI-factor is thus smaller than the value for a homogeneous Callisto but it is considerably larger than the value of 0.3105 ± 0.0018 found for Ganymede (Anderson et al., 1996). It was shown (Anderson et al., 2001; Sohl et al., 2002) that it is impossible to construct 2- or 3-layer models of Callisto that will fit the MoI-factor and have a pure ice shell above a rock–iron core or shells of the pure phases ice, silicate, and iron, as it is possible and widely accepted for Ganymede. This is generally taken to suggest that differentiation in Callisto has not run to completion. Stevenson (1998), for example, has proposed a struc-

Interior structure and differentiation of Callisto

ture with pure ice I overlying an ocean and a rock/ice core. While incomplete differentiation of Callisto appears to be a reasonable explanation, it is not clear how an icy satellite can evolve into such a state. It has been proposed that (partial) melting in near surface layers of undifferentiated Ganymede or Callisto sized satellites will lead to separation between the ice and the rock + metal components (Friedson and Stevenson, 1983; McKinnon, 1997; Anderson et al., 2001). Mueller and McKinnon (1988) and McKinnon (1997) have argued that convective layering due to the endothermic ice II to V phase transition and compositional layering may cause melting. These processes have not been modeled in detail, though, and in part rely on uncertain parameterizations of convective heat transport and compositionally driven flow. It is certainly true that the difference in gravitational energy between a homogeneous satellite and a satellite with an ice mantle and a rock + metal core is more than sufficient to provide the latent heat of ice melting. The gravitational energy stored in the rock sediment at the bottom of a 200 to 300 km deep ocean above an undifferentiated deeper interior, however, is not sufficient to melt the rest of the ice. A self-accelerating mechanism like the one proposed by Friedson and Stevenson (1983) has to be invoked. In any case, the consensus is that Callisto must have either formed cold such that melting can be avoided or must have differentiated so slowly as to remove the heat liberated upon differentiation. Recent models of the formation of the galilean satellites (Canup and Ward, 2002; Mosqueira and Estrada, 2003a, 2003b), therefore, form Callisto cold. But these models do not explain why and how it is incompletely differentiated. In this paper, we present a model for the incomplete differentiation of a solid Callisto by the gradual and slow unmixing of rock and ice.

2. A simple model of incomplete differentiation The mean density of Callisto of 1860 kg m−3 (Anderson et al., 2001) suggests that the satellite’s volume fraction of rock (and metal), ε0 , must be between roughly 20 and 40 vol%, mostly depending on the rock and the ice densities. If proto-Callisto consisted of a homogeneous mixture of metal bearing, perhaps hydrated rock and ice, then the density difference between the two phases will cause the rock particles to gradually settle towards the center if the ice deforms readily enough. The rate at which the settling will occur depends on the size of the rock particles and the viscosity of the ice in addition to its dependence on the density difference between ice and rock. An ice mantle depleted of rock will form above a rock–ice core in which the volume fraction of rock will increase with time. The settling rate will gradually approach zero as the volume fraction of rock in the core reaches the close-packing limit (CPL). The approximate minimum value of the fractional radius of the rock–ice core RCPL /RS will then be given by

403

Fig. 1. Ice mantle shell density and rock density in the rock–ice core as a function of core radius and core composition. The models satisfy the mass and moment of inertia factor calculated for Callisto from the Galileo gravity data. Cores larger than 2000 km in radius are unlikely because the ice shell density of less than 1175 kg m−3 will be too small to account for the high pressure ice phase transitions and rock remaining in the ice.

the cubic root of ε0 /εCPL , where εCPL is the rock volume fraction in the CPL. The latter depends on the size distribution and the shape of the particles and varies between 60 vol% for a homogeneous distribution of spheres of equal radius to about 90 vol% implying core radii between 1700 and 2000 km. It is possible to construct two-layer models of Callisto that satisfy the mass and the MoI-factor constraints and that have the core consist of ice and rock in the CPL (Fig. 1). The density of the rock component varies between 3400 kg m−3 for a CPL of 60 vol% and 2500 kg m−3 for an assumed maximum CPL of 90 vol% in these models. The CPL for randomly sized rock particles is likely to be larger than the random CPL for equally sized spheres recently determined to be as large as 0.74 (Torquato et al., 2000). The density of 3400 kg m−3 is close to the density of metal bearing anhydrous chondritic matter (Endress et al., 1996; Hutcheon et al., 1998), while the value of 2500 kg m−3 is consistent with hydrated rock and strongly oxidized iron phases such as magnetite (e.g., Lodders and Fegley, 1998). Given enough time, water will hydrate rock and irreversibly oxidize free metal (e.g., Wasson, 1974). Since the rock will be releasing radiogenic heat, thermal convection may interact with the settling. Rudman (1992, 1997) has developed a one-field theory of flow for particlesettling in the presence of thermally driven convection and has solved the equations of flow numerically in Cartesian geometry. In his theory, the matrix is a fluid whose viscosity is a function of the concentration of the (rigid) particles with µ(ε) = µ/f (ε) and f (ε) = (1 − ε/εCPL )β . The exponent β is about equal to 1.55. It was shown that a layer depleted of particles forms and grows in thickness at the top of the mixture while the concentration of particles increases with time in the bottom layer. Convection keeps the concentration in the bottom layer approximately homoge-

404

K. Nagel et al. / Icarus 169 (2004) 402–412

neous. The separation front between the two layers moves with the particle speed, which is given by v = f (ε) · S0 , where S0 = ρd 2 g/18µi is the well-known Stokes speed. ρ is the difference in density between the particles of diameter d and the fluid, g is the acceleration due to gravity, and µi is the (constant) viscosity of the pure fluid, in this case ice. The model is consistent with experimental data reported by Koyaguchi et al. (1990). Figure 4 further below shows the particle concentration (panel (b)) and temperature distribution (panel (d)) in the two layers for our own calculation in spherical, axisymmetric geometry. The separation front between the two layers is clearly developed and sharp. The temperature field shows the characteristics of thermally driven convection. These results can be used to construct a simple analytical model for the evolution of Callisto’s core that we will discuss before we introduce a numerical model in which we have modified and extended the Rudman model by allowing for temperature dependence of the viscosity in addition to its dependence on the particle concentration, spherical geometry, ice phase changes, and depth dependence of gravity. Let R denote the radial distance of the unmixing front to the center of the satellite. Then dR = −f (ε) · S0 dt implies that     RCPL 3 β dR =− 1− · S0 dt R

(a)

(2)

(3)

and β   S0 dε ε · . =− 1− dt εCPL RS

(4)

Equations (3) and (4) show that the unmixing front will asymptotically approach R = RCPL as the rock concentration gradually approaches εCPL . Figure 2a shows the general solution of Eqs. (3) and (4) and the moment of inertia calculated using Eq. (1) as functions of time in units of RS /S0 . Because the value of β is uncertain and experimental data (Durham et al., 1992) suggest values of 2 to 2.5, we have included results of a calculation with β = 2.5. The most rapid evolution occurs in the first unit of time. Thereafter, the distance to the unmixing front, the MoIfactor, and the rock concentration in the core slowly evolve towards their final equilibrium values. For β = 2.5 the evolution proceeds at smaller rates; taking three time units to proceed to similar values as in one time unit for β = 1.55. The question then is for what parameter values is RS /S0 equal to about one billion years? Reasonable values for the viscosity of ice are between 1014 and 1018 Pa s (see Durham and Stern (2001) for a recent review of ice rheology). With g equal to 1.25 m s−2 and ρ a few 103 kg m−3 , the size of the rock particles must be meters to tens of meters or larger (compare Fig. 2b). It is clear from these considerations that the model will not work if Callisto accreted as a more or less homogeneous mixture of snow and dust flakes. Rather, the model requires a mix of ice and at least boulder to satellitesimal size rock particles.

(b)

Fig. 2. (a) Evolution of the rock–ice core of Callisto in a simple model with temperature independent viscosity. The results are cast in terms of (R − RCPL )/(RS − RCPL ) (solid lines), (ε0 − ε)/(ε0 − εCPL ) (dash-dotted lines), and (I − ICPL )/(I0 − ICPL ) (dashed lines), where I0 = 0.4 MRS2 . R is the radial distance to the unmixing front, or the core radius and ε is the volumetric rock concentration in the core. The subscript CPL refers to the close packing limit and the subscript 0 to initial values for a homogeneous Callisto. RS is the surface radius of the satellite. The models have been calculated by integrating Eqs. (3)–(4) and by calculating the moment of inertia I from Eq. (1). The black lines represent models with β = 1.55, the red lines those with β = 2.5. (b) Rock particle size versus ice viscosity for constant values of RS /S0 between 0.5 and 3 Ga.

Interior structure and differentiation of Callisto

3. A more complete model 3.1. Model description In this section we describe the results of more detailed calculations where we solve the one-field equations of Rudman (1997) for temperature and rock concentration dependent rheology in spherical geometry while taking into account high-pressure ice phase transformations as well as the dependence of gravity on depth and of radiogenic heat production on time (see also Nagel, 2001, for more details of the model calculations). Solid-state creep of ice at the stress level expected in large icy satellites occurs by non-Newtonian grain boundary sliding and/or (at high stress levels) dislocation creep (e.g., McKinnon, 1998; Durham and Stern, 2001); the exponent to the stress is between 2 and 6. To reduce computational cost, we construct a Newtonian analogue by reducing the activation energy by a factor of 5/3 (Christensen, 1984; Forni et al., 1991). The viscosity of the ice µi is then given by    Tm (P ) 3 −1 , µi = µm · exp A (5) 5 T where µm is the viscosity of the ice extrapolated to the melting temperature, Tm (P ) is the pressure dependent melting temperature (e.g., Chizhov, 1993), and A a dimensionless parameter. Equation (5) is a simplified but commonly used rheological model (Durham and Stern, 2001). The applicability of this equation (using the homologous temperature Tm /T ) even for ice I at pressures below the triple point has been demonstrated experimentally by Durham et al. (1997) and for ice VI by Sotin et al. (1985). The activation parameter A is usually taken to be between 18 and 36 which

405

correspond to an activation energy between about 40 and 80 kJ mol−1 . Although the rheological parameters vary considerably between the various ice phases, we use a single rheological parameter set as given in Table 1, for simplicity. The presence of ammonia dihydrate and salts may affect the rheology. While the effects of salt are little understood (Durham and Stern, 2001), there are experimental data on the rheology of ammonia dihydrate–water mixtures (Durham et al., 1993). However, for reasonable concentrations of ammonia dihydrate of less than 45 wt% (equivalent to NH3 concentrations of less than 15 wt%) the effect on the rheology parameters is small and within the above range of parameter values. The presence of NH3 will cause a substantial reduction of the melting temperature for pressures up to 800 MPa, or 400 km depth (Hogenboom et al., 1997; Sotin et al., 1998). Therefore, temperature at the same homologous temperature will differ from the case of pure water ice and there may even be a substantial ocean (Spohn and Schubert, 2003). A detailed modeling of these effects is beyond the scope of the present paper. The viscosity of the rock–ice mixture is µi . µ(ε) = (6) f (ε) High-pressure ice phase transitions may have significant effects on the flow depending on their Clapeyron slopes. The endogenic ice II to V transition may cause layering of the flow (Bercovici et al., 1986; Mueller and McKinnon, 1988). However, the II to V transition is not likely to occur in Callisto when the ice mantle is convecting. We have checked the validity of this statement a posteriori. Exogenic phase transitions at shallow depths are likely to have negligible effects on the sluggish flow in the ice mantle. The ice VI to VII transition, finally, is temperature independent. Therefore, we

Table 1 Parameter values thought to be characteristic of Callisto and used for the modeling Radius Surface gravity Ice density

Heat capacity of ice Rock density Average rock volume fraction Rock volume fraction in the CPL Thermal expansion coefficient Thermal conductivity Thermal diffusivity Ice viscosity constant Ice viscosity activation parameter Radiogenic heat production rate (t = 0) Decay constant Scaling temperature difference

RS g0 ρI ρII ρIII ρV ρVI ρVII cP ρR ε0 εCPL α k κ µm A H0 λ T0

2400 km 1.2 m s−2 920 kg m−3 1193 kg m−3 1157 kg m−3 1235 kg m−3 1322 kg m−3 1680 kg m−3 1.5 × 103 J kg−1 K−1 3500 kg m−3 0.26 0.60 10−4 K 3 W m−1 K−1 2 × 10−6 m2 s−1 1014 Pa s 20 25 pW kg−1 1.43 × 10−17 s−1 150 K

Thermodynamic properties are from the collection of Sotin et al. (1998), rheological parameters have been estimated using the data in Durham and Stern (2001), and κ using the data in Ross and Kargel (1998).

406

K. Nagel et al. / Icarus 169 (2004) 402–412

have neglected phase transitions in the thermal buoyancy terms of the momentum conservation equation to be introduced below. The phase transitions have been taken into account in the compositional buoyancy terms and for the calculation of the particle velocity by making the density difference between rock and ice a function of r ρ(r) = ρR − ρi ,

ri
r < ri+1 ,

(7)

where ρR is the density of the rock component and ρi the density of the relevant ice phase in the appropriate depth interval ri to ri+1 (compare Table 1). Gravity depends on radius according to GM(r) (8) , r2 where G is the universal gravity constant, and M(r) is calculated using the ice densities ρi , ρR , and the radially varying concentration of the rock component. The radiogenic heat production rate in rock per unit mass, finally, is the usual function of time

g(r) =

H = H0 · exp(−λt).

(9)

We proceed by introducing the dimensionless and scaled one-field equations for mass, momentum and energy conservation. We use the surface radius RS as length scale and RS2 /κ, where κ is the thermal diffusivity, as time scale; µm scales the viscosity and the assumed initial temperature difference between the surface and the deep interior T0 = 150 K scales the temperature. The density difference ρ(r) is scaled by (ρR − ρVI ), where ρVI is the density of ice VI. Gravity is scaled by its surface value g0 and H is scaled by H0 RS2 /kT0 . ∇ u = 0,





(10)

∇P + ∇µ∇ u = g(r)RaT T − Fρ(r)ε , (11)    ∂T  (1 − ε)ρi St + u + ε(1 − ε)ρ(r) uS · ∇T 1+ ρ(r) ∂t = ∇ 2 T + εH, ∂ε = −∇ · εuR , ∂t where

(12)

u = εuR + (1 − ε) ui

(14)

and

(13)



uS = ui − uR =

0

Σ − ρ(r)g(r) µ

 .

(15)

In addition to already defined quantities, u is the mixtureaveraged flow velocity, uR the phase averaged rock particle velocity, and ui is the phase averaged ice velocity (averaged over a neighborhood of the size of the numerical grid spacing); uS is the relative velocity between ice and rock and equals the Stokes velocity. RaT is the thermal Raleigh number RaT ≡

αg0 ρVI T0 RS3 , κµm

(16)

where α is the thermal expansion coefficient, F is the ratio between the thermal and chemical buoyancy forces F≡

ρR − ρVI , αρVI T0

(17)

Σ is the dimensionless and scaled Stokes velocity, Σ≡

(ρR − ρVI )g0 d 2 RS 18µm κ

(18)

and St is the Stefan number, St ≡

L . cP Tm

(19)

The Stefan number accounts for thermal effects of ice melting. We cannot model large scale melting and the flow of liquid water. However, we model the effects of incipient melting by accounting for the latent heat. We have tested in additional calculations in Cartesian geometry (Nagel, 2001) the importance of accounting for small melt fractions in the buoyancy term in Eq. (11) and found it to be negligible. The term including the Stefan number in Eq. (12) has been derived as follows: Where the temperature attains the melting temperature, a term (1 − ε)ρi L∂Vm /∂t is added to the energy equation to account for the latent heat of melting. L denotes the latent heat of melting per unit mass and Vm is the volumetric concentration of melt. The latent heat of ice I is 3.3 × 105 J kg−1 (e.g., Lodders and Fegley, 1998). The latent heats of other ice polymorphs are believed to be within ±20% of this value (Friedson and Stevenson, 1983). Univariant phase transitions and sharp transition boundaries are difficult to model in numerical calculations. Therefore, we assume that melting occurs in a narrow temperature interval Tm centered at the melting temperature Tm . For Tm we take 3 K. We further assume that the melt volume grows linearly through Tm and obtain ∂Vm /∂t = 1/Tm ∂T /∂t. The term ρi L/Tm acts as an additional heat capacity in the energy balance equation and will prevent a rapid temperature rise above Tm . Division of ρi L/Tm by ρcP gives the ratio between the latent heat released per unit volume and unit temperature change and the heat capacity per unit volume (1 − ε)ρi /ρSt. Values are between 75 for pure ice and 15 for a rock–ice mixture in the CPL. The results of our model calculations, however, suggest a more representative number of 45. Note that for the time scales R2 RS = Σ −1 S . S0 κ

(20)

With the parameter values of Tables 1 and 2, RS /S0 is approximately 1 Ga. (RCPL is approximately 1800 km.) The flow is in the so-called supercritical regime of two-phase convection (Hansen and Yuen, 1987) since F  1 (Table 2). To solve these equations, we developed a code based on the spherical, axisymmetric version of ConMan by S.D. King. The coupled mass and momentum equations (10) and (11) are solved using the penalty method to eliminate the

Interior structure and differentiation of Callisto

407

Table 2 Dimensionless parameter values used for the model calculations in Section 3 Rayleigh number Buoyancy ratio Stokes velocity parameter Stefan number Scaled density difference

Ice I Ice II Ice III Ice V Ice VI Ice VII

Scaled viscosity Scaled radiogenic heat production rate Average rock volume fraction Rock volume fraction in the CPL

pressure dependence. The diffusion–advection energy equation (12) and the rock advection equation (13) are solved using a Streamline Upwind/Petrov Galerkin Method and a Predictor–Corrector time-stepping algorithm with a Courant time step criterion. The numerical method introduces some artificial numerical diffusion (Eq. (13) has no diffusion). Care had to be taken to use a grid size that kept the effects of the diffusion small enough. Figures 4(b) and 4(d) show that the algorithm can resolve the particularly sharp differentiation front for temperature independent viscosity. The grid has 200 lateral and 100 radial grid points. To avoid the singularity at the center of the sphere we exclude a small rocky central core of radius RI = 100 km from the grid. The boundary conditions are constant temperature Ts = 0.8 at r = RS and zero heat flow at r = RI , zero particle flux at both r = RS and r = RI , free slip at r = RS , and zero radial and lateral velocity at r = RI . Because of limitations in computer memory and speed we use a Rayleigh number of 109 and an initial scaled heat production rate H of 320 as listed together with the other dimensionless model parameters in Table 2. (We have also calculated models with Rayleigh numbers of 107 and 108 .) These values are too small in comparison with realistic values for Callisto and therefore our results have to be extrapolated. With the parameter values listed in Table 1 and a scaling temperature difference of T0 = 150 K the thermal Rayleigh number for Callisto is approximately 1012 , about three orders of magnitude larger than used in the present calculations. Since the heat flow in internally heated thermal convection scales as approximately Ra 0.2 , the heat production rate in equilibrium with the heat flow extrapolated from the scaled and dimensionless value of H (at t = 0) of 320 (Table 2) will be 1274. This value is equivalent to the chondritic heating rate of 25 pW kg−1 . The assumed radiogenic heat production rate accounts for long-lived isotopes only. Short lived isotopes such as 26 Al are neglected. The heat production rate of the latter exceeds that of long lived isotopes at the time of formation of the calcium–aluminum-rich inclusions CAIs by about a factor of 104 (e.g., Merk et al., 2002). After about 10 half-lives or

RaT F Σ St ρI ρII ρIII ρV ρVI ρVII µ H ε0 εCPL

109 50 100 45 1.178 1.053 1.070 1.034 0.995 0.831   exp 12(Tm /T − 1)   320 × exp −1.43 × 10−17 × t 0.26 0.60

7 Ma the heat production rate of 26 Al has become unimportant in comparison with the heat production rate of the long lived isotopes. Whether or not Al was effective will therefore depend on the accretion time scale and its relation to the formation of the CAIs, which preceded the formation of the Solar System. Since the recent models for the formation of Callisto form the satellite slowly, 26 Al was neglected. 3.2. Results Figure 3 shows profiles of the temperature, the ice melting temperature, and the rock concentration at times of t = 0, 0.0066, 0.013, and 0.02 dimensionless units equivalent to 0, 0.66Rp /S0 , 1.3Rp /S0 , and 2Rp /S0 , or, with the parameter values of Table 1, 0, 0.66, 1.33, and 2.0 Ga. Figures 4(a) and 4(c) show the rock concentration and the temperature at t = 0.02 (2Rp /S0 or 2.0 Ga). Figures. 4(b) and 4(d), finally, show results of a run with temperature independent viscosity for comparison. This model has parameter values as the model with temperature dependent viscosity but RaT = 105 to save computer time and is taken at t = 0.01, equivalent to 1Rp /S0 or 1 Ga. The concentration profiles in Fig. 3 show how a depleted ice mantle forms and develops. There are several differences with respect to the case where the viscosity is temperature independent: First, an undifferentiated layer remains immediately underneath the surface in which the rock concentration decreases towards the value in the depleted layer and second, the depleted mantle layer is not completely cleared of rock particles. The concentration gradients together with the temperature gradients will cause thermo-chemical convection in the layer above the differentiation front. In the rock/ice core, major concentration gradients are absent and convection is driven by temperature gradients as in the constant viscosity case. The low surface temperature and the high viscosity at shallow depth cause the undifferentiated near-surface layer. This layer is equivalent to a rheological lithosphere; it is eroded from below with time and as it thins and becomes further depleted of rock, it feeds the depleted mantle below.

408

K. Nagel et al. / Icarus 169 (2004) 402–412

(a)

(b)

(c)

(d)

Fig. 3. Profiles of temperature (solid line), volumetric rock concentration (dashed), and ice melting temperature (dotted) versus radial distance from the center of Callisto at times of t = 0 (a), 0.0066 (b), 0.013 (c), and 0.02 (d) dimensionless units equivalent to 0, 0.66Rp /S0 , 1.3Rp /S0 , and 2Rp /S0 , or, with the parameter values of Table 1, 0, 0.66, 1.33, and 2.0 Ga.

It is not likely that endogenic processes have ruptured the roughly 150 km thick lithosphere. This is consistent with the observation of an ancient, heavily cratered surface on Callisto. The unmixing front is initially not very well developed but the concentration gradient increases in steepness as the depleted mantle becomes more depleted. Within about 1 Ga, the unmixing front has almost reached a depth of 600 km, the final depth in the CPL. The front then moves only insignificantly; instead, the concentration further increases in the core and decreases in the ice mantle. After 2 Ga, the rock concentration in the core has reached a value of 0.46, still significantly smaller then the CPL of 0.6, and convection in the core is still dominated by the ice rheology. The concentration in the core thus increases considerably more slowly then in the simple temperature-independent viscosity theory. This is a consequence of the temperature dependence of the viscosity. However, the difference in terms of the separation front movement is not too large: Had we assumed a

viscosity of 2 × 1014 Pa s in calculating RS /S0 , then the latter would have been 2 Ga and the result of the simple theory would have been recovered. Heat is transferred through the depleted ice mantle by convection. The temperature profile through the lithosphere and the depleted mantle over time develops a boundary layer structure: The top thermal boundary layer comprises the lithosphere and the top of the depleted mantle. The bottom boundary layer extends from the depleted mantle to the unmixing front where it is joined by the top thermal boundary layer of the core. The Nusselt number, the dimensionless ratio between the surface heat flow and the surface heat flow in case of pure thermal conduction, is small, about 1.5. Figure 4 shows the rock concentration and the temperature at t = 2.0 Ga. For comparison, we also show (panels (b) and (d)) results for a similar model with temperature independent viscosity. In the latter case, the concentration distribution is simple. The ice layer is clear of rock particles and the core has a homogeneous concentration of 0.46.

Interior structure and differentiation of Callisto

409

(a)

(b)

(c)

(d)

Fig. 4. Rock concentration (a) and temperature (c) for the model with temperature dependent viscosity at time t = 0.02 or 2 Ga. Parameter values of this model are given in Table 2. For comparison, panels (b) and (d) show the results for a model with temperature-independent viscosity at a Rayleigh number of 105 and at time t = 0.01 or 1 Ga. Ice phase transitions and melting are neglected in the latter model.

The front between the two regions is extremely sharp. Convection in the absence of lateral concentration gradients is driven by thermal buoyancy. The thermal plumes in the ice mantle are clearly developed and drive the convection pattern in the core underneath. With temperature dependent viscosity (panels (a) and (c)), rock concentration gradients below the unmixing front in the core are still negligibly small. Across the unmixing front, the rock concentration decreases strongly in the radial direction but lateral variations are small. In the depleted mantle we observe a downwelling plume with a substantial rock concentration. This plume erodes the lithosphere from below and transfers rock to the core. The downwelling plume is also cold (panel (c)) but its main driving force is the negative buoyancy of the rock. The corresponding upwelling plumes are found underneath the poles, as it is typical for rotationally symmetric spherical shell convection. The plumes are elements of a two-cell convection pattern (l = 2) in the mantle. Temperature in the upwelling plumes reaches the melting temperature and melting occurs locally. The flow in the boundary layers and plumes can be discussed with the help of Fig. 4 and Eqs. (14) and (15). In the downwelling plume negative buoyancy is provided both by the relatively large rock concentration and the cold temperature. The rock flows past the ice in the downwelling plume with a relative velocity given by the Stokes velocity (15). In

the boundary layers the relative azimuthal velocity is zero. The radial relative velocity is again given by Eq. (15). This, together with the large lateral extent of the boundary layer, causes the separation front between the two layers to uniformly sink at the Stokes speed. Conservation of mass requires an upflow which occurs in the two upwelling plumes. The (positive) buoyancy here is entirely thermal and the plumes are mostly ice with the ice flowing past the rock at approximately negative Stokes speed. Combining Eqs. (14) and (15), we have u = εus + ui .

(21)

This shows that a positive (upward) mixture average velocity requires the ice velocity to be larger then −εus .

4. Discussion and conclusions We presented a model capable of explaining an incompletely differentiated Callisto slowly evolving from a homogeneous initial state to one with an ice–rock lithosphere in which the rock concentration decreases with increasing depth, an ice mantle depleted to a large extent of rock, and a rock–ice core that can satisfy the mass and moment of inertia constraints. Our conclusions do not depend significantly on the value of the CPL. Internal consistency requires

410

K. Nagel et al. / Icarus 169 (2004) 402–412

that the latter limit increase with decreasing rock density as Fig. 1 shows. The model predicts that the CPL will never be reached in the core since the particle speed will approach a value of zero in this limit. However, it is conceivable that the ideal model breaks down, for instance, because of melting in the top layers of the core or other non-ideal situations. We can, therefore, ask the question of what will happen if the CPL is reached. Although this is beyond our present modeling, we can present some simple calculations. In the CPL, convection as we have modeled it will stop altogether because the viscosity will be equal to the rock viscosity which is too large for convection at the relevant temperatures. Heat will be transferred by conduction and, possibly, by porous media convection, the details of the latter process depending on the porosity in the CPL and the permeability. It is unlikely however, that Callisto will then experience a rapid differentiation event. Warming the ice in the entire core to the melting temperature takes at least 1 Ga and melting the ice requires at least another 1 Ga, assuming that the heat production rate is constant and the core adiabatic. The times for warming the core to the melting temperature and melting the ice has been calculated as follows: According to Fig. 3d the core is on the average about 70 K away from the melting temperature. Assuming an adiabatic core, the time to warm by a given T is dtw = ρcP T /εH ≈ 107 T years, ρ ≈ 2000 kg m−3 is the average density in the core, and H is the present day chondritic heat production rate per unit volume of 10−8 W m−3 . The time to melt the ice is estimated from dtm = (1 − ε)ρi L/εH ≈ 1 Ga, where L is the latent heat of ice equal to 3.3 × 105 J kg−1 . Since the radiogenic heating rate will continue to decrease with time, large scale melting may never occur. Rather, the water may continue to react with silicates and metal to form more hydrous rock and oxides. This may result in some compaction. Further compaction will require that the silicates warm to temperatures where they can deform over geological time scales which will take at least as long as melting the ice in the core. Thus whatever differentiation may follow after the CPL has been reached will occur slowly. As a consequence we speculate that the slow differentiation we describe in this paper may continue beyond the CPL. Callisto may then end as a satellite with a mostly hydrated core underneath an ice layer. Because metal oxidation is irreversible and largely removes the density difference between free iron and silicate rock, it will then likely never form an iron-rich core. The time for differentiation in our model depends on particle size and ice rheology. There is considerable leeway in both parameters to make the model credible. Among the most important simplifications, though, is the assumption of a single size of the rock particles. Instead, there should be a range of particle sizes extending from perhaps micrometers to kilometers or more. Because of the differences in the Stokes velocities, some sorting may be expected in the core with rock particle size decreasing with increasing radial distance from the center. The less than meter size particles will

sink very slowly but increase the viscosity and by that slow the whole differentiation process. We have run one model in which we added dust particle of a size of 10−3 times that of the rock particles diameter. The effect was an increase of the time to reach the state shown in Fig. 3d by 500 Ma. From Fig. 2b, we can conclude that a reasonable timescale for running to the close packing limit is possible for particle sizes between 1 and 100 m if the viscosity is between 1013 and 1017 Pa s. For particle sizes of about 10 m the viscosity should be between about 1014 and 1015 Pa s. The size of rock and ice particles in the jovian nebula and in proto-Callisto is basically unknown. Solid matter can be transferred to the circum-jovian nebula effectively by gas drag for particle sizes  1 m (e.g., Canup and Ward, 2002), the exact size certainly depending on the density of the nebula. The satellitesimals, already consisting of rock/ice mixtures, can decouple from the gas and grow to satellites for sizes > 100 m (Canup and Ward, 2002; Mosqueira and Estrada, 2003a). The satellitesimals grow to this size by “sweeping up dust and rubble” (Mosqueira and Estrada, 2003a). The rock part of the rubble may end up as rock particles in proto-Callisto. The size range required by our model is consistent with the size range of sub-satellitesimals suggested by the recent accretion models. The little depleted ice–rock lithosphere is gravitationally unstable. The density difference between the lithosphere and the depleted ice mantle below depends on the density of the rock component and its concentration and should be of the order of a few 100 kg m−3 , roughly an order of magnitude larger than the thermal density difference between the two layers. The lid forms because the viscosity at near surface temperatures is orders of magnitude larger than the viscosity in the warmer and depleted layer. It is thus analogous to the well-known stagnant lid that forms in strongly temperature dependent viscosity convection (e.g., Davaille and Jaupart, 1993; Solomatov, 1995; Grasset and Parmentier, 1998). The thickness of the lid depends on the rheology, in particular on the parameter A which measures the rate of change of viscosity with temperature, and on the heat flow from the convecting layer underneath the lid. Since we are using a radiogenic heating rate that is four times smaller than the chondritic value, the lithosphere may actually be about four times thinner than the roughly 150 kilometers that we find in the present calculations. It may be argued that the little depleted ice–rock lithosphere should be subject to a Rayleigh–Taylor instability. However, one can view the chemically driven convection in our model as a Rayleigh–Taylor instability with temperature and concentration dependent rheology in the additional presence of thermal buoyancy. That the time scales are reasonable can be rationalized by applying the simple growth time for Rayleigh–Taylor instabilities given in, e.g., Turcotte and Schubert (1982, p. 257). The growth time is τ=

13.4µ , ρgl0

Interior structure and differentiation of Callisto

where l0 is a characteristic length. The latter cannot be the thickness of the lithosphere because τ will vary by many orders of magnitude through the lithosphere. Instead we take l0 to be the length across which the viscosity changes by a factor of e (if we use, instead, an order of magnitude viscosity variation, our conclusion will not be much different). The temperature at the base of the lithosphere is between 140 and 150 K. Calculating the viscosity at 150 K with the parameter values of Table 1 we find a growth time of hundreds of million years. For 140 K we already get growth times of Gyr. This shows that a lithosphere can survive for the age of the Solar System. The Galileo magnetometer data suggest that Callisto has a subsurface ocean (Khurana et al., 1998; Neubauer, 1998; Kivelson et al., 1999). Interpretation of the data is not unique with respect to depth and electrical conductivity of the ocean but a depth around 200 km is most plausible (Zimmer et al., 2000). Our results do not offer a simple explanation for an ocean at that or even shallower depth. The melt patches in the upwelling plumes are of local not global extent and may not explain the observed signal. Our results confirm an observation by Spohn and Schubert (2003) that a shallow ocean is unlikely in Callisto if ice at shallow depth contains a significant amount of rock. Melting is more conceivable in our model (Fig. 3d) in the transition region between the depleted mantle and the rock–ice core at about 600 km depth. Is it possible to have models with the rock–ice core extending to 200–300 km depth? According to Fig. 1, 400 km is a minimum depth to the core. But even at 2000 km core radius the density above the core is at 1175 kg m−3 unrealistically small considering ice phase transformations and the rock still in the ice layer. However, if the ice is not pure H2 O but contains other components such as ammonia, methane and/or salts (e.g., Kargel, 1998) then a near-surface ocean is possible. With 5 wt% ammonia, for instance, the ocean would extend from around 80 to 200 km depth (Spohn and Schubert, 2003). Ganymede and Callisto are similar in mass, radius and, therefore, most likely in composition. A model for Callisto may therefore be challenged to explain Ganymede. The Ganymede–Callisto dichotomy has long been recognized, however, and speculated to be attributable to differences in accretion, accretional heating rates, composition and, perhaps, differences in tidal heating (e.g., Squyres 1983; Friedson and Stevenson, 1983; McKinnon and Parmentier, 1986; Kirk and Stevenson, 1987; Mueller and McKinnon, 1988; Tittemore, 1990; Malhotra, 1991; McKinnon, 1997). Recent work on the accretion of the galilean satellites have further elucidated these differences (Canup and Ward, 2002, Mosqueira and Estrada 2003a, 2003b). Thus, the two satellites may have had initial conditions too different to allow the same evolutionary path. In any case, Ganymede must have experienced large scale ice melting while Callisto may have not. The formation of the iron core in Ganymede after the separation of the rock + iron from the ice can be attributed to heating by radioactives (e.g., Spohn and Breuer, 1997;

411

Schubert et al., 2003). More modeling is required to test evolutionary scenarios such as the runaway melting proposal of Friedson and Stevenson (1983) to get a better understanding of how and why the two satellites evolved differently.

Acknowledgments We have profited from discussions with G. Schubert, J.T. Wasson, W. McKinnon, and D.J. Stevenson. C. Sotin and D.J. Stevenson provided thoughtful reviews. This research was supported by the Deutsche Forschungsgemeinschaft.

References Anderson, J.D., Lau, E.L., Sjogren, W.L., Schubert, G., Moore, W.B., 1996. Gravitational constraints on the internal structure of Ganymede. Nature 384, 541–543. Anderson, J.D., Jacobson, R.A., McElrath, T.P., Moore, W.B., Schubert, G., Thomas, P.C., 2001. Shape, mean radius, gravity field and interior structure of Callisto. Icarus 153, 157–161. Bercovici, D., Schubert, G., Reynolds, R.T., 1986. Phase transitions and convection in icy satellites. Geophys. Res. Let. 13, 448–451. Canup, R.M., Ward, W.R., 2002. Formation of the galilean satellites: conditions of accretion. Astrophys. J. 124, 3404–3423. Chizhov, V.E., 1993. Thermodynamic properties and equation of state of high-pressure ice phases. Translated from. Prikl. Mek. i Tekhn. Fiz. 2, 113–123. Christensen, U., 1984. Convection with pressure- and temperaturedependent non-Newtonian rheology. Geophys. J. R. Astron. Soc. 77, 343–384. Davaille, A., Jaupart, C., 1993. Transient high-Rayleigh-number thermal convection with large viscosity variations. J. Fluid Mech. 253, 141–166. Durham, W.B., Stern, L.A., 2001. Rheological properties of water ice— application to satellites of the outer planets. Annu. Rev. Earth Planet. Sci. 29, 295–330. Durham, W.B., Kirby, S.H., Stern, L.A., 1992. Effects of dispersed particulates on the rheology of water ice at planetary conditions. J. Geophys. Res. 97, 20883–20897. Durham, W.B., Kirby, S.H., Stern, L.A., 1993. Glow of ice in the ammonia– water system. J. Geophys. Res. 98, 17667–17682. Durham, W.B., Kirby, S.H., Stern, L.A., 1997. Creep of water ices at planetary conditions: a compilation. J. Geophys. Res. 102, 16293–16302. Endress, M., Zinner, E., Bischoff, A., 1996. Early aqueous activity on primitive meteorite parent bodies. Nature 379, 701–703. Forni, O., Coradini, A., Federico, C., 1991. Convection and lithospheric strength in Dione, an icy satellite of Saturn. Icarus 94, 232–245. Friedson, A.J., Stevenson, D.J., 1983. Viscosity of rock–ice mixtures and applications to the evolution of icy satellites. Icarus 56, 1–14. Grasset, O., Parmentier, E.M., 1998. Thermal convection in a volumetrically heated, infinite Prandtl number fluid with strongly temperaturedependent viscosity: implications for planetary thermal evolution. J. Geophys. Res. 103, 18171–18181. Hansen, U., Yuen, D., 1987. Evolutionary structures in double-diffusive convection in magma chambers. Geophys. Res. Lett. 14, 1099–1102. Hogenboom, D.L., Kargel, J.S., Consolmagno, G.J., Holden, T.C., Lee, L., Buyyounouski, M., 1997. The ammonia–water system and the chemical differentiation of icy satellites. Icarus 128, 171–180. Hutcheon, I.D., Krot, A.N., Keil, K., Phinney, D.L., Scott, E.R.D., 1998. Manganese-53/chromium-53 dating of Fayalite in Mokoia: evidence for asteroidal alteration of CV chondrites. Meteorit. Planet. Sci. 32, A72– A73.

412

K. Nagel et al. / Icarus 169 (2004) 402–412

Kargel, J.S., 1998. Physical chemistry of ices in the outer Solar System. In: Schmitt, B., de Bergh, C., Festou, M. (Eds.), Solar System Ices. Kluwer, Dordrecht, pp. 3–32. Khurana, K.K., Kivelson, M.G., Stevenson, D.J., Schubert, G., Russel, C.T., Walker, R.J., Polanskey, C., 1998. Induced magnetic fields as evidence for subsurface oceans in Europa and Callisto. Nature 395 (6704), 777– 780. Kirk, R.L., Stevenson, D.J., 1987. Thermal evolution of a differentiated Ganymede and implications for surface features. Icarus 69, 91–134. Kivelson, M.G., Khurana, K.K., Stevenson, D.J., Bennett, L., Joy, S., Russell, C.T., Walker, R.J., Polanskey, C., 1999. Europa and Callisto: induced or intrinsic fields in a periodically varying plasma environment. J. Geophys. Res. 104, 4609–4625. Koyaguchi, T., Hallworth, M.A., Huppert, H.E., Sparks, R.S.J., 1990. Sedimentation of particles from a convecting fluid. Nature 343, 447–450. Lodders, K., Fegley, B., 1998. The Planetary Scientist’s Companion. Oxford Univ. Press, New York, 371 p. Malhotra, R., 1991. Tidal origin of the Laplace resonance and the resurfacing of Ganymede. Icarus 94, 399–412. McKinnon, W.B., 1997. Mystery of Callisto: is it undifferentiated? Icarus 130, 540–543. McKinnon, W.B., 1998. Geologic landforms and processes on icy satellites. In: Schmitt, B., de Bergh, C., Festou, M. (Eds.), Solar System Ices. Kluwer, Dordrecht, pp. 525–549. McKinnon, W.B., Parmentier, E.M., 1986. Ganymede and Callisto. In: Burns, J.A., Matthews, M. (Eds.), Satellites. Univ. of Arizona Press, Tucson, pp. 718–763. Merk, R., Breuer, D., Spohn, T., 2002. Numerical modelling of 26 Alinduced radioactive melting of planetesimals considering accretion. Icarus 159, 183–191. Mueller, S., McKinnon, W.B., 1988. Three-layered models of Ganymede and Callisto, compositions, structures, and aspects of evolution. Icarus 76, 437–464. Mosqueira, I., Estrada, P.R., 2003a. Formation of the regular satellites of giant planets in an extended gaseous nebula I: subnebula model and accretion of satellites. Icarus 163, 198–231. Mosqueira, I., Estrada, P.R., 2003b. Formation of the regular satellites of giant planets in an extended gaseous nebula II: satellite migration and survival. Icarus 163, 232–255. Nagel, K., 2001. Thermal-compositional convection and the differentiation process in Callisto. PhD thesis. Westfälische Wilhelms-Universität Münster.

Neubauer, F., 1998. Planetary science—oceans inside Jupiter’s moons. Nature 395, 749–751. Ross, R.G., Kargel, J.S., 1998. Thermal conductivity of Solar System ices, with special reference to martian polar caps. In: Schmitt, B., de Bergh, C., Festou, M. (Eds.), Solar System Ices. Kluwer, Dordrecht, pp. 33–62. Rudman, M., 1992. Two-phase natural convection: implications for crystal settling in magma chambers. Phys. Earth Planet. Inter. 72, 153–172. Rudman, M., 1997. One-field equations for two-phase flows. J. Austral. Math. Soc. Ser. B 39, 149–170. Schubert, G., Anderson, J.D., Spohn, T., McKinnon, W.B., 2003. Interior composition, structure and dynamics of the galilean satellites. In: Bagenal, F., Dowling, T.E., McKinnon, W.B. (Eds.), Jupiter. Univ. of Arizona Press, Tucson. Sohl, F., Spohn, T., Breuer, D., Nagel, K., 2002. Implications from Galileo observations on the interior structure and chemistry of the galilean satellites. Icarus 157, 104–119. Solomatov, V.S., 1995. Scaling of temperature- and stress-dependent viscosity. Phys. Fluids 7, 266–274. Sotin, C., Poirier, J.P., Gillet, P., 1985. Creep of high-pressure ice VI. In: Klinger, J., Benest, D., Dollfus, A., Smoluchowski, R. (Eds.), Ices in the Solar. Reidel, Dordrecht, pp. 109–118. Sotin, C., Grasset, O., Beauchesne, S., 1998. Thermodynamic properties of high pressure ices: implications for the dynamics and internal structure of large icy satellites. In: Schmitt, B., de Bergh, C., Festou, M. (Eds.), Solar System Ices. Kluwer, Dordrecht, pp. 79–96. Spohn, T., Breuer, D., 1997. Interior structure and evolution of the galilean satellites. In: Celnikier, L.M., Than Van, J.T. (Eds.), Planetary System: The Long View, Frontieres ed., pp. 135–144. Spohn, T., Schubert, G., 2003. Oceans in the galilean satellites of Jupiter? Icarus 161, 456–467. Stevenson, D.J., 1998. An ocean within Callisto? AGU Fall Meeting. Abstract P12 B-10, on CD-ROM. Tittemore, W.C., 1990. Chaotic motion of Europa and Ganymede and the Ganymede–Callisto dichotomy. Science 250, 263–267. Torquato, S., Truskett, T.M., Debenedetti, P.G., 2000. Is random close packing of spheres well defined? Phys. Rev. Lett. 84, 2064–2067. Turcotte, D.L., Schubert, G., 1982. Geodynamics. Wiley, New York, 450 p. Wasson, J.T., 1974. Meteorites: Classification and Properties. SpringerVerlag, New York, 316 p. Zimmer, C., Khurana, K.K., Kivelson, M.G., 2000. Subsurface oceans on Europa and Callisto: constraints from Galileo magnetometer observations. Icarus 147, 329–347.