Thermal evolution of Kuiper belt objects, with implications for cryovolcanism

Icarus 202 (2009) 694–714

Contents lists available at ScienceDirect

Icarus www.elsevier.com/locate/icarus

Thermal evolution of Kuiper belt objects, with implications for cryovolcanism Steven J. Desch a,∗ , Jason C. Cook b , T.C. Doggett a , Simon B. Porter a a b

School of Earth and Space Exploration, Arizona State University, P.O. Box 871404, Tempe, AZ 85287-1404, USA Southwest Research Institute, 1050 Walnut St., Suite 300, Boulder, CO 80302, USA

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 25 July 2008 Revised 4 February 2009 Accepted 8 March 2009 Available online 14 March 2009

We investigate the internal thermal evolution of Kuiper belt objects (KBOs), small (radii <1000 km), icy (mean densities <2500 kg m−3 ) bodies orbiting beyond Neptune, focusing on Pluto’s moon Charon in particular. Our calculations are time-dependent and account for differentiation. We review evidence for ammonia hydrates in the ices of KBOs, and include their effects on the thermal evolution. A key finding is that the production of the first melt, at the melting point of ammonia dihydrate, ≈176 K, triggers differentiation of rock and ice. The resulting structure comprises a rocky core surrounded by liquids and ice, enclosed within a >100-km thick undifferentiated crust of rock and ice. This structure is especially conducive to the retention of subsurface liquid, and bodies the size of Charon or larger (radii >600 km) are predicted to retain some subsurface liquid to the present day. We discuss the possibility that this liquid can be brought to the surface rapidly via self-propagating cracks. We conclude that cryovolcanism is a viable process expected to affect the surfaces of large KBOs including Charon. © 2009 Elsevier Inc. All rights reserved.

Keywords: Charon Ices Spectroscopy Thermal histories Trans-neptunian objects Volcanism

1. Introduction Volcanism on Earth, where silicates are the dominant rockforming mineral, plays an important role in reshaping the surface of this planet. Cryovolcanism is the process by which subsurface water, which as ice is a major “rock”-forming “mineral” in KBOs, is brought to the surface. Strong evidence for cryovolcanism comes from spacecraft images of the icy satellites of the giant planets. Moving outward in the Solar System, Europa has been recognized since the Voyager mission to have an extremely young surface (Smith et al., 1979), ∼50 Myr (Zahnle et al., 2003). Support for resurfacing by cryovolcanism exists (Figueredo and Greeley, 2004), although its role on large scales, and the extent of diapirism in resurfacing Europa, are currently debated (Fagents, 2003). Measurements by the Galileo magnetometer strongly indicate the existence of a subsurface ocean on Europa (Kivelson et al., 2000), supporting a role for cryovolcanism. Evidence for heat fluxes in excess of what could be produced by radioactive decay suggest tidal dissipation drives this surface activity (Ruiz and Tejero, 2000). The bright terrain of Ganymede has been interpreted as cryovolcanic in origin, with liquid filling low-lying grabens (Showman et al., 2004). High inferred heat fluxes in the past suggest tidal dissipation as the driver of this activity as well (Nimmo et al., 2002; Nimmo and Pappalardo, 2001; Dombard and McKinnon, 2001). Some of the surface of Enceladus appears extremely uncratered and young, <108 yr, and shows tectonic features highly indicative of cryovol-

*

Corresponding author. Fax: +1 (480) 965 8102. E-mail address: [email protected] (S.J. Desch).

0019-1035/$ – see front matter doi:10.1016/j.icarus.2009.03.009

©

2009 Elsevier Inc. All rights reserved.

canism (Kargel and Pozio, 1996). The Cassini mission has revealed the south polar region is especially young and active, with a plume of icy material escaping the moon from this region (Porco et al., 2006; Hansen et al., 2006). Whether the plume is attributable to the boiling of liquid water near the surface (Porco et al., 2006) or to decomposition of clathrates (Kieffer et al., 2006) is unknown at this time. The heat fluxes inferred to exist within Enceladus (Giese et al., 2008) also exceed the fluxes driven by radioactive decay alone. The Cassini Titan Radar Mapper observed a large dome or shield volcano, three calderas and, several extensive flows, consistent with cryovolcanism (Lopes et al., 2007). The inferred rheological properties of the flows are consistent with water–ammonia slurries (Lopes et al., 2007). Titan is large enough to maintain subsurface liquid by radiogenic heating alone, although tidal heating probably also contributed (Grasset et al., 2000; Tobie et al., 2005). Miranda and Ariel have relatively young surfaces in places, no older than 0.1–1 Gyr (Plescia, 1987, 1988; Zahnle et al., 2003), as well as grabens and flow features that have been interpreted as cryovolcanic (Kargel, 1994). There are no obvious sources of tidal heating for Miranda and Ariel, although orbital resonances in the past cannot be excluded (Peale, 1988). Images of Triton (itself likely a captured KBO: Agnor and Hamilton, 2006) show large uncratered terrains that are among the youngest in the Solar System, <10 Myr (Schenk and Zahnle, 2007). These regions contain large depressions that have apparently been extensively modified by flooding, melting, faulting and collapse, and several episodes of partial filling and removal of material, features that have been interpreted as cryovolcanic in origin (Croft, 1990; Kargel and Strom, 1990; Kargel, 1994). Triton could have been substantially tidally heated

Cryovolcanism on Kuiper belt objects

as a consequence of its capture (McKinnon, 1984), and there is evidence for high heat fluxes due to tidal dissipation in Triton’s past (Ruiz, 2003). In short, there is abundant evidence for cryovolcanism on the surfaces of icy satellites. It is likely that in each of these cases, tidal dissipation enhanced the heating in these satellites. It is therefore not known whether KBOs of similar size should have also experienced cryovolcanism or not. Cryovolcanism on KBOs has also been argued for, based on observations of crystalline water ice or ammonia/ammonia hydrates. Crystalline water ice is converted to amorphous ice by doses of Galactic cosmic rays (GCRs) ∼3 eV per molecule (Mastrapa and Brown, 2006; Cook et al., 2007, and references therein), and the spectral signatures of ammonia hydrates are likewise removed by doses ∼100 eV per molecule (Strazzulla and Palumbo, 1998). Using the GCR fluxes of Cooper et al. (2003), Cook et al. (2007) calculated that the NIR spectral signatures of crystalline water ice and ammonia hydrates (which probe to depths ∼350 μm) should be erased in times 1.5 and 50 Myr, respectively. When the ability of solar ultraviolet photons to amorphize ice was considered, these times were reduced to 0.03 and 1 Myr (Cook et al., 2007). The presence of these ices on the surface of a KBO therefore demands a physical explanation, possibly cryovolcanism. At the same time, ammonia or ammonia hydrates are especially suggestive of cryovolcanism because ammonia depresses the melting point of ice, from 273 to 176 K, enabling ice to be melted and mobilized (Hogenboom et al., 1997). Presumably any liquid brought to the surface would be ammonia-rich. Jewitt and Luu (2004) observed crystalline water ice on the KBO 50000 Quaoar, based on absorption at 1.65 μm, as well as absorption at 2.21 μm consistent with ammonia hydrates. More recently, Schaller and Brown (2007) have argued that the absorption at 2.21 μm is predominantly due to adsorbed methane, rendering it unclear whether ammonia hydrates exist on Quaoar’s surface or not. Cook et al. (2007) also observed crystalline water ice and ammonia hydrates on Pluto’s moon, Charon, and argued for cryovolcanism as a resurfacing process on this KBO. More recently, Barucci et al. (2008) have claimed detection of crystalline water ice and ammonia ices on the KBO 90482 Orcus. These KBOs may be experiencing cryovolcanism. The use of crystalline water ice as a diagnostic of surface renewal is strengthened by the apparent lack of physical mechanisms capable of converting amorphous ice back to crystalline form at the cold temperatures of KBO surfaces. Cook et al. (2007) quantitatively analyzed many possible annealing mechanisms and found none of them to be competitive against amorphization. The most effective mechanism they considered was thermal annealing by deposition of kinetic energy by micrometeorite impacts, which they calculated could anneal ice on Charon over timescales ∼5 Myr. This is to be compared to the amorphization by GCRs, which takes ∼1.5 Myr, and by solar ultraviolet (UV) photons, which Cook et al. (2007) calculate could take as little as 3 × 104 yr. It is possible that the micrometeorite dust fluxes within the Kuiper Belt are larger than those measured between 1 and 18 AU by Pioneer 10 (Humes, 1980) by orders of magnitude. These fluxes will be measured by the Student Dust Counter instrument on the New Horizons mission when it reaches the Pluto–Charon system in 2015 (Horányi et al., 2008). If they are very large, then localized heating by micrometeorite impacts may be why water ice exists in the crystalline state. Otherwise, according to the analysis of Cook et al. (2007), cryovolcanism may be the only explanation for the crystalline water ice and ammonia hydrates on the surfaces of KBOs. It is worth clarifying that while Cook et al. (2007) interpreted the spectral signature of crystalline water ice to be diagnostic of cryovolcanism on KBOs, their work should not be construed to mean that icy satellites with crystalline water ice necessarily experience cryovolcanism. Crystalline water ice is observed on essentially all icy satellites with water ice on their surfaces (Grundy et

695

al., 1999), including some with heavily cratered, ancient surfaces (e.g., Oberon), and others too small to have experienced internal heating (e.g., Nereid: Brown, 2000). However, due to the enhanced orbital velocities of satellites around their planets, as well as the increased dust fluxes around the giant planets, it is likely that micrometeorites could anneal ice on satellite surfaces, even if they are not quite able to anneal the ice on KBOs (Porter et al., 2009, in preparation). Observational diagnostics such as crystalline water ice being inconclusive, then, the question remains, do KBOs experience cryovolcanism? For now, the question must be answered approximately, with theoretical models. Direct imaging of Pluto and Charon will be possible in 2015 with the arrival of the New Horizons probe (Stern, 2008), and should largely settle the question of whether Pluto or Charon experience cryovolcanism. But even then, proper interpretation of the data will depend on a theoretical understanding of cryovolcanism on KBOs. We are strongly motivated, therefore, to study the possibility. This article represents a theoretical modeling effort to assess whether KBOs can experience cryovolcanism. A mechanism for delivering liquid to the surface is described in Section 6, and the focus in this paper is on whether subsurface liquid can persist on small, icy bodies, and Charon in particular. A prerequisite for maintenance of liquid is ammonia, and we begin in Section 2 by reviewing ammonia abundances on icy satellites, comets, and KBOs. In Section 3 we describe a numerical code written to calculate the thermal evolution of a KBO, including the effects of ammonia, differentiation, and heating by decay of long-lived radionuclides. Results for a Charon-like KBO are described in detail in Section 4, and sensitivity of the results to key parameters is explored in Section 5. We then draw conclusions in Section 6 concerning the persistence of subsurface liquid in KBOs, and the likelihood of cryovolcanism on them. We explore the implications of subsurface liquid in KBOs in Section 7. 2. Ammonia in the early Solar System The potential role of ammonia in depressing the melting point of water and affecting the rheology of ice has long been appreciated (Stevenson, 1982; Consolmagno, 1985; Croft et al., 1988; Kargel et al., 1991; Kargel, 1992; Durham et al., 1993; Hogenboom et al., 1997). A key parameter in this discussion is X , the weight fraction of ammonia relative to ammonia plus water. Based on the cosmochemical abundances of elements in the solar nebula, this primordial ratio could have been as high as X = 14.5%, assuming all N condensed into NH3 (Lodders, 2003). In fact, a significant fraction of the N must have condensed into organics, and an unknown fraction could have existed in the form of N2 , which was not obviously accreted. It is therefore difficult to predict from first principles the value of X in icy bodies in the early Solar System. Nonetheless, a variety of lines of evidence suggest that X is at least a few percent. Theoretical models (Charnley and Rodgers, 2002) predict that N2 is efficiently dissociated into atomic form in molecular clouds, and that the pathways to return the atomic N back to N2 are inhibited by a lack of oxygen radicals in the gas phase. Instead, a substantial amount of N is reduced to NH3 which then freezes onto icy grain mantles. This NH3 may then undergo photolysis by ultraviolet radiation in the icy mantles, which converts it to N-bearing organics. Once in the Solar System, atomic N may be converted to N2 or NH3 . This scenario is also consistent with the fractionation of N isotopes in various Solar System nitrogen reservoirs, as well. Observations of molecular clouds confirm a lack of gas-phase N2 (Maret et al., 2006), as predicted by the chemistry models. These same models predict that as the Solar System forms, roughly 25% of the N is in the form of NH3 ice (Charnley and Rodgers, 2002).

696

S.J. Desch et al. / Icarus 202 (2009) 694–714

Likewise, Maret et al. (2006) model the chemistry in the Barnard 68 core and infer X ≈ 7%. So, even if no NH3 formed within the Solar System from the remaining atomic N, a large value of X (≈ 4–7%) might be expected. On the other hand, observations of Solar System objects do not support such high values of X . Perhaps the best observational constraints come from observations of comets. The abundance of ammonia can be fixed by observing it directly in the radio spectrum, or by observing UV and optical lines of its photodestruction products, NH2 and NH. From NH2 observations, a value X ≈ 0.5% was found to be typical among a large database of comets (Kawakita and Watanabe, 2002). Note, though, that these same authors derived a range of values, from 0.25 to 1.0%. Many comets are depleted in ammonia. For example, 21P/Giacobini–Zinner is considered depleted (Beaver et al., 1990) and a value X < 0.1% is inferred from the NH2 observations (Kawakita and Watanabe, 2002). C/1996 B2 (Hyakutake) is inferred to have X ≈ 0.3% from radio lines of NH3 (Palmer et al., 1996) and X = 0.35% from lines of NH (Meier et al., 1998). Comet 1P/Halley has slightly higher ammonia abundance. From lines of NH a value X = 0.44–0.94% is inferred (Feldman et al., 1993), and from optical lines of NH2 , X = 0.75 ± 0.18% is inferred (Kawakita and Watanabe, 2002). Direct measurements from mass spectrometers on the Giotto probe yield values closer to 1.5%: 0.8–2.0% using the neutral mass spectrometer (Meier et al., 1994), and X = 1–2% using the ion mass spectrometer (Allen et al., 1987). Finally, radio observations of NH3 in C/1995 O1 (Hale–Bopp) yield X = 1.0–1.8% (Bird et al., 1997). The combined evidence seems to point towards short-period comets that have made repeated passages around the Sun being more depleted in ammonia. This implies that the primordial abundance of ammonia in comets must exceed the highest abundances observed in comets, X ≈ 2%. Observational searches for ammonia or ammonia hydrates on KBOs and icy satellites are less conclusive. Strong evidence exists for the presence of ammonia hydrates on the surface of Charon (Brown and Calvin, 2000; Dumas et al., 2001; Cook et al., 2007), based on absorption at 2.21 μm. Ammonia hydrates also appear present on Quaoar (Jewitt and Luu, 2004), although this absorption has also been interpreted as absorption by methane (Schaller and Brown, 2007). Tentative hints of ammonia hydrate absorption have also been reported for Miranda (Bauer et al., 2002), and evidence for ammonia on Enceladus is accumulating (Emery et al., 2005; Verbiscer et al., 2006). Intriguingly, observations of cryovolcanic plains on Miranda and Ariel yield estimates of the viscosity of the cryolavas; these are most consistent with a water–ammonia slurry (Schenk, 1991), with X large enough to reduce the viscosity. This implies X on the order of a few percent or more. Spectral measurements unfortunately probe only the top millimeter of each body’s surface, and do not provide a reliable estimate of the global values of X . Radiation damage at the surface, or cryovolcanism itself, may drastically alter the inferred values of X . Nonetheless, the spectral fits of Cook et al. (2007) (their model 2) do suggest that the surface ices are roughly 20–30% “ammonia hydrate”. The “ammonia hydrate” component is modeled based on optical constants obtained from experimental data by Brown et al. (1988), who froze a 3% molar solution of NH4 OH in water, or X ≈ 3%. The (globally averaged) value of X on the Charon’s surface is therefore probably ∼1% within the observed ices. To summarize, theoretical models of solar nebula chemistry indicate that X ≈ 7% or more are possible. Direct evidence for ammonia in the outer Solar System comes from observations of comets, although at levels closer to X ≈ 1–2%. Morphological features on icy satellites and the presence of ammonia hydrates on Charon support X ≈ a few percent on these bodies as well. In what follows, we consider X = 1% and X = 5% to bracket the reasonable range of ammonia abundances with which KBOs form.

3. Internal thermal evolution model 3.1. Background We have written a numerical code to calculate the thermal evolution of KBO interiors. This code, like many others, essentially solves the 1-D (spherical) heat conduction equation. As discussed below, a time-dependent solution is absolutely necessary. The role of differentiation is also of critical importance. Commonly solutions consider a completely homogeneous body or a completely differentiated body. As we show below, however, the maintenance of subsurface liquid is optimized when a rocky core is in contact with an icy layer, both of them surrounded by a crust of an undifferentiated rock/ice mixture. Also, as is well known, the effects of ammonia are important in lowering the thermal conductivity and viscosity of water ice, and especially in suppressing the melting point. To study the internal thermal evolution of KBOs, it is therefore essential to construct a time-dependent code, including the role of ammonia, and allowing for partial differentiation. In this section we describe such a code we have written. To our knowledge, ours is the first code including all the relevant effects (time-dependence, ammonia, partial differentiation) to be applied to KBOs. Previous codes used to study KBOs and icy satellites have been incomplete in one or more areas. Ruiz (2003) used a steady-state model to estimate the depth to the bottom of the ice layer on Triton’s surface. Similarly, Hussmann et al. (2006) considered steady-state, homogeneous or fully differentiated models of icy satellites and KBOs. De Sanctis et al. (2001) considered the thermal evolution of porous KBOs, and Shchuko et al. (2006) modeled the thermal evolution of Varuna, a porous, relatively small KBO. These models did not include differentiation (appropriate considering the small size of the objects they considered). Thermal evolution models of icy satellites were developed by Ellsworth and Schubert (1983), but did not apparently follow the change in structure due to differentiation. Prialnik and BarNun (1990) studied the thermal evolution of icy satellites subject to heating by short-lived and long-lived radionuclides but did not include differentiation or ammonia. McKinnon (2002) considered the thermal evolution of KBOs, but also did not include differentiation. The code used by Schubert et al. (2007), used to study Enceladus, incorporates time evolution and differentiation but does not include ammonia. To our knowledge, only one code described in the literature includes all of the relevant effects. The code used by Castillo-Rogez et al. (2007a) is a time-evolution model that includes ammonia and differentiation. This code has been applied to Iapetus (Castillo-Rogez et al., 2007a) and Rhea (Castillo-Rogez, 2006) and Ceres (Castillo-Rogez et al., 2007b), but has not been applied to KBOs. 3.2. Description of model Our code assumes 1-D spherical geometry and divides the KBO into N radial zones (typically N ≈ 200). The zones are evenly distributed in r so that the outer radius of each zone obeys r i = R p (i / N ), where R p is the KBO radius. The code tracks the timevarying internal energy E i (t ) within each zone. These are set at t = 0 to the energies corresponding to the temperature T surf ; that is, we assume the KBOs accreted “cold” (did not retain accretional heat). Given the energy at some time t, the energy at a later time t + t is updated, including contributions from radiogenic heating and the fluxes of heat energy into and out of each zone: E i (t + t ) − E i (t )

t

= Q i (t ) + 4π r i2−1 F i −1 − 4π r i2 F i .

(1)

Here F i refers to the heat flux out of zone i, into zone i + 1. At the innermost zone, F i −1 vanishes, since by symmetry F = 0 at

Cryovolcanism on Kuiper belt objects

697

Table 1 Adopted long-lived radionuclide heating parameters. Isotope 40

K Th 235 U 238 U 232

t 1/2 (Gyr)

Initial # per 106 Si

wt% (ppm)

E (MeV)

H0 (mW/kg elem)

1.265 14.0 0.704 4.47

5.244 0.04399 0.00592 0.01871

1.389 0.0676 0.0092 0.0295

0.6087 39.0 42.7 46.1

0.0254 0.0253 0.5444 0.0913

H0 (pW/kg rock) 35.2 1.70 5.02 2.69

r = 0. An “equation of state” is used to convert the energy E i in each zone to a temperature T i . At the outermost zone, the energy in that zone is not calculated; rather, the temperature is fixed: T N = T surf . This equation evolves the temperature in a physically reasonable way so long as the timestep satisfies a Courant stability criterion: t < min[(r )2 /(2κ )], where κ is the thermal diffusivity. In practice we find that rather than readjusting the timestep it is simpler to fix the timestep at t = 1000 yr, as this satisfies the Courant condition for all of the cases we consider and still allows the full thermal evolution to be computed in a few hours on a desktop computer (Apple G5, 2 GHz processor). The thermal evolution is then completely determined once the rate of radiogenic heating and the heat fluxes are determined, and once an equation of state for the matter, T ( E ), can be established. These quantities are in turn affected by the amounts of rock, ices and liquids in each zone, which are affected by differentiation. We now discuss each of these effects. 3.3. Radiogenic heating We consider radiogenic heating by four long-lived radionuclides only: 40 K, 232 Th, 238 U and 235 U. On the premise that formation of several-100-km diameter KBOs will take longer than ∼10 Myr, we do not include decay of 26 Al or 60 Fe here, although we intend to consider their effects in future work. The rate of radiogenic heating due to an isotope x with half-life t 1/2 , per mass of that isotope, is easily found to be ( E )x (ln 2/t 1/2 )/m X , where m X is the mass of an atom of x, and  E X represents the heat energy per decay. These  E X include heating due to emission of α and β particles and γ rays, but not the emission of neutrinos, which escape planets altogether. Because of the loss of neutrino energies, the radiogenic heating rate cannot be determined from the mass deficit between the parent and daughter nuclei alone. The uncertainties in the neutrino energy during radioactive decay are often on the order of 10%, and lead to uncertainties almost as large in  E X . We list our assumed values of t 1/2 and  E X for each radionuclide in Table 1. For each radionuclide, we also list in Table 1 each one’s assumed initial mass. The masses of each isotope decrease exponentially with time according to their half lives. Details on how we estimated neutrino energy losses and computed the initial abundances of each radionuclide are provided in Appendix A. We calculate an initial radiogenic heating rate of S 0 = 44.6 pW/kg of rock, dominated by decay of 40 K. For comparison, Ellsworth and Schubert (1983) used S = 52 pW/kg, and Multhaup and Spohn (2007) used S 0 = 45 pW/kg. For the present day, Ellsworth and Schubert (1983) used S = 6 pW/kg while Multhaup and Spohn (2007) find S = 5.2 pW/kg; we find, in comparison, S = 5.6 pW/kg. The variation of heating rate with time is depicted in Fig. 1. 3.4. Fluxes We generally assume that the heat fluxes are carried by conduction, so that F = −k∂ T /∂ r. We finite-difference this as Fi = −

ki + ki +1 2

T i +1 − T i . (r i +1 − r i −1 )/2

(2)

Fig. 1. Adopted heating rate per mass of rock due to the radioactive decay of 40 K, 232 Th, 235 U, and 238 U, as a function of time relative to the present day (t = 0 on this graph only).

As the energy in the outermost zone (E N ) is fixed, F N is not calculated. By symmetry F 0 = 0. Below (Section 3.6) we describe how we calculate the total thermal conductivity of a mixture of solids. There will arise in our calculations many instances where the heat fluxes are not carried by conduction alone but by convection as well, especially in the layer of pure water ice that we predict will develop during the thermal evolution of a KBO. A common prescription for determining whether a fluid heated from below will convect is to determine the Rayleigh number, comparing it to a critical value. In this context, which assumes a uniform viscosity η across the fluid layer, the fluid will convect if the Rayleigh number between the top of the layer (here, at r = R diff , the boundary between the ice layer and the undifferentiated crust, discussed below) and the base of the ice layer (at r = R base ) exceeds the critical value, ≈1100 (Spohn and Schubert, 2003). The Rayleigh number is calculated as Ra =

αρice g ( T )(r )3 , κη

(3)

where r = R diff − R base and  T = T ( R base ) − T ( R diff ). We calculate the other quantities—the expansivity α , the gravitational acceleration g, the thermal diffusivity κ and especially the viscosity η —at the midpoint of the ice layer. The thermal diffusivity κ = k/ρice /cP , where k is the thermal conductivity and c P is the heat capacity. Evaluation of these quantities is discussed below (Sections 3.6 and 3.7). This simple model for the onset of convection may not apply to icy satellites, however, where the viscosity contrasts across the fluid layer are great. Instead, satellites are understood to convect in the “stagnant lid” regime, where much of the ice layer will be high-temperature and low-viscosity, and will participate in vigorous convection; this is capped on top (and in some treatments, at the bottom as well) by a stagnant lid of ice that transports heat by conducting fluxes alone (Solomatov, 1995; McKinnon, 1998, 1999; Multhaup and Spohn, 2007). The onset of convection in the stagnant lid regime may require even higher Rayleigh numbers, ∼106 (McKinnon, 1999), so the onset of convection is not ensured in our simulations without modeling stagnant lid convection. Unfortunately, a true simulation of stagnant lid convection lies outside the scope of this paper. Inventing an algorithm for parameterizing fluxes in the stagnant lid regime is one obstacle; numerically resolving the temperature structure in the relatively thin stagnant lid is another.

698

S.J. Desch et al. / Icarus 202 (2009) 694–714

Despite the difficulties in modeling the stagnant lid regime, we can still approximate the thermal structure within a KBO by the following prescription. First, we are not particularly interested in the thermal structure within the ice layer per se, just the total temperature drop across it. Two regimes are relevant. The first is the case where convection is presumed to be inhibited, and fluxes are carried by thermal conduction alone. This case is easily accommodated in our models. The second case is where stagnant lid convection is initiated. In that case, we assume there is a nearly isothermal profile within the ice layer (at T = T c ), capped at the top and bottom by stagnant boundary layers across which the flux is carried by thermal conduction (e.g., Ellsworth and Schubert, 1983; Multhaup and Spohn, 2007). The thickness of the boundary layers is δ = (r )/Nu, where Nu is the Nusselt number of the convecting layer. The temperature drop across the ice layer is essentially the temperature drop across the stagnant lids, which is  T = F δ/k = F (r )/(k × Nu), where k is the thermal conductivity. Here it is understood that k represents the average of the thermal conductivities in the two stagnant lids on top and bottom, which will be equivalent to the thermal conductivity in the convecting layer. Our code will correctly compute the same temperature drop across the ice layer if the effective thermal conductivity keff is chosen appropriately; i.e., if  T = F (r )/keff = F (r )/(k × Nu), or if keff = k × Nu. This is how we implement convection in our code: within the convecting ice layer we increase the thermal conductivity by a factor of the Nusselt number. To calculate the Nusselt number, we adopt the parameterization Nu = (Ra/Rac )0.25 for Ra > Rac (Spohn and Schubert, 2003), where the Rayleigh number is calculated at the midpoint of the ice layer. While this treatment does not accurately reflect the likely temperature profile within the ice layer (nearly isothermal, except in the stagnant lids), it should well approximate the temperature drop across the ice layer. Our treatment of convection is crude, but does provide a means for calculating the effects of convection on the thermal evolution within a KBO, assuming convection is initiated. We have also performed runs where convection was shut off. As we will discuss below (Section 5.4), the mass of liquid within a KBO, and the time for which it can resist freezing, are not very sensitive to how convection is treated. Typically the mass of liquid varies by several tens of percent, and the time for which it stays liquid varies by a few ×108 years, depending on whether or not convection is shut off. So despite the crudeness of our treatment of convection, we expect it to describe the temperature drops across the ice layer adequately enough for our purposes. Within layers containing liquid, the Rayleigh number almost always greatly exceeds Rac because the viscosity of water is so low, yielding typical values of Nu ∼ 104 . The actual temperature gradient in the convecting liquid layer could be found by using a correspondingly high effective thermal conductivity ∼104 W m−1 K−1 ; however, this approach is inadvisable in the liquid layer, as the code is limited by a Courant condition from taking timesteps t < (r )2 /(2κ ), and high thermal conductivities would lead to proportionally small timesteps (103 yr). On the other hand, the physical significance of the large effective thermal conductivity is simply that the temperature drop across the liquid layer is small. It is likely, in fact, to be <10−2 K, but none of the main conclusions of this paper is altered if the temperature drop is allowed to be as large as 1 K. Within any layer with at least 2% liquid (by weight), we therefore impose an effective thermal conductivity keff = 400 W m−1 K−1 . (Low values of the liquid mass fraction presumably refer to heat transport by fluid flow through the rock, rather than liquid convection per se.) For a maximum heat flux <10 mW m−2 , the temperature gradient therefore will not exceed 0.025 K km−1 , and across a liquid layer 40 km thick the temperature drop will not exceed 1 K. This is a small discrepancy compared to the temperatures typically achieved in the core, which are

at least several ×102 K, and typically the heat fluxes and the temperature discrepancies are even smaller. We therefore assume that k = keff in any layers containing at least 2% liquid. We note that even the small temperature gradient this imposes is basically eliminated in the next step of our algorithm, which homogenizes the internal energy of the slush layer after each timestep (Section 3.9). 3.5. Ice viscosity In assessing the extent of differentiation or the effectiveness of convection, we must calculate the viscosity, especially of the water ice layer. At the low stress levels to be found in KBO interiors it is expected that the dominant deformation mechanism for ice is dislocation-accommodated grain boundary sliding (McKinnon, 1999; Durham and Stern, 2001). The effective viscosity is ηeff = σ /(3 ˙ ), where ˙ is the strain rate, σ is the stress, and d the grain size of the ice, and where the other constants are parameters dependent on the phase and temperature of the ice (Durham and Stern, 2001). We adopt the strain rate associated with Nabarro– Herring creep, in which case the strain rate is given by

˙ = 9.3

σ b3 D ν d2 kT

(4)

,

where b = 4.52 × 10−10 m is the Burgers vector, D ν = D 0ν exp(− Q / R T ) is the lattice diffusion frequency, D 0ν = 9 × 10−4 m2 s−1 , Q = 59.4 kJ mol−1 is the activation energy for selfdiffusion, and R is the gas constant (Song et al., 2006). Based on these parameters, we calculate



η ≈ 3.8 × 1014

2

d 1 mm

 exp 26.2

Tm − T T

 Pa s,

(5)

where T m = 273 K. For grain sizes slightly smaller than a millimeter, this formula matches well the commonly used viscosity law of Thomas et al. (1987):





η( T ) = η0 25

273 K T

 −1 ,

(6)

where η0 = 1.0 × 1014 Pa s. We choose to adopt the viscosity law of Thomas et al. (1987) in our runs. Many of the parameters entering into the more exact viscosity law are uncertain, and this viscosity law is commonly used and simple to implement. The difference will be slight in any case, because the viscosity is only ever used to change the effective thermal conductivity in the convecting ice layer. The conductivity is multiplied by the Nusselt number, which scales as the fourth root of the viscosity, so different treatments will not change the viscosity much. Moreover, relevant quantities such as the mass of liquid and the time at which it all freezes are fairly insensitive to whether the ice layer is convecting or is not, so minor changes in the viscosity law will not be expected to matter much anyway. Using the viscosity law of Thomas et al. (1987), we estimate that at temperatures just below 176 K, the viscosity in pure water ice is ≈1 × 1020 Pa s. If this ice contained some proportion of ADH, the viscosities would be lowered slightly. In an ice with X = 4.7%, just below 176 K the viscosity ≈1018 Pa s, while in an ice with X = 9.5%, the viscosity ≈ 1017 Pa s (Durham et al., 1993; Arakawa and Maeno, 1994). At temperatures just above 176 K, the viscosity of pure water ice is still ≈1 × 1020 Pa s. In an ice bearing ADH or ammonia, however, the melting of the ADH just above 176 K dramatically lowers the viscosity, to <1012 Pa s for a mixture with X = 4.0%, and to <1010 Pa s for a mixture with X = 8.4% (Arakawa and Maeno, 1994). It is not clear what the viscosity should be for an ice mixtures with abundance X ≈ 1%, as is inferred from observations of comets. Clearly there is some threshold between 0% and 4%, below which ammonia does not significantly

Cryovolcanism on Kuiper belt objects

699

alter the viscosity of the ice; probably it is on the order of 1–2%, but this threshold is unknown. 3.6. Thermal conductivities Besides the heat fluxes, no quantity is as important for determining the temperature gradient within an icy body as its thermal conductivity. We are somewhat surprised, then, that most of the models cited above have paid so little attention to this fundamental input. The conductivity of crystalline ice is well known to increase with decreasing temperature (e.g., Klinger, 1980), and the conductivities of rock can be quite variable. Nevertheless, with little or no justification in many cases, particular values of the thermal conductivities are merely imposed. For example, Hussmann et al. (2006) assume a temperature-independent kice = 3.3 W m−1 K−1 in their models of KBOs. Schubert et al. (2007) in their modeling of Enceladus assume a temperature-independent kice = 2.1 W m−1 K−1 and krock = 3 W m−1 K−1 . In their modeling of the KBO Varuna, Shchuko et al. (2006) assume kice = 5.67( T /100 K) W m−1 K−1 (Klinger, 1980), but a temperatureindependent value krock = 2.4 W m−1 K−1 . Castillo-Rogez et al. (2007a) in their modeling of Iapetus took kice = [0.47 + 4.88( T / 100 K)] W m−1 K−1 and krock = 4.2 W m−1 K−1 . Finally, McKinnon (2002) similarly adopted krock = 2 W m−1 K−1 , but recognized that the conductivities of amorphous quartz and feldspar near 100 K can be much lower, ≈1 W m−1 K−1 (Cahill et al., 1992). To summarize, modelers sometimes adopt a temperature-varying formula for the conductivity of (crystalline water) ice based on laboratory measurements, but usually adopt an arbitrary value for the conductivity of rock. The thermal conductivity of rock is an essential input, and laboratory measurements of suitable analogs are demanded. The typical values adopted in the literature, krock ≈ 2–4 W m−1 K−1 , are appropriate for terrestrial rocks at room temperature, but are not obviously applicable to the type of rock from which KBOs formed, at the typical temperatures of KBOs. Instead, we adopt the thermal conductivities measured for ordinary chondrites by Yomogida and Matsui (1983). Perhaps because of the aggregate nature of chondrites, these conductivities typically do not exceed 1 W m−1 K−1 ; also, they exhibit no variation with temperature from 100 K to over 500 K. We therefore adopt krock = 1.0 W m−1 K−1 , with no variation with temperature. Note that this means that rock in our model is effectively an insulator in comparison with ice. At higher temperatures, we do not expect the conductivities to increase significantly. We note that practically all terrestrial rocks tend to exhibit decreased thermal conductivities at higher temperatures (Clauser and Huenges, 1995). In fact, a variety of sedimentary, volcanic, plutonic and metamorphic terrestrial rocks have thermal conductivities that are nearly independent of temperature above about 600 K, generally not exceeding ≈2 W m−1 K−1 . Sedimentary rocks (limestone, dolomite, shale and quartz sandstone), which may be the closest terrestrial analogs to the rocks in KBOs, generally are lower in conductivity, ≈1–1.5 W m−1 K−1 (Clauser and Huenges, 1995). We accordingly consider the likely range of thermal conductivities to be krock = 1–2 W m−1 K−1 , and we favor the lower end of this range. For the thermal conductivity of water ice we adopt kH2 O(s) = 5.67( T /100 K)−1 W m−1 K−1 (Klinger, 1980). For the thermal conductivity of ammonia dihydrate, we use the data of Lorenz and Shandera (2001) and Kargel (1992) to infer a temperatureindependent conductivity kADH = 1.2 W m−1 K−1 . (See also the discussion in Ruiz, 2003, who came to the same conclusion.) For the sake of simplicity, we combine the conductivities in a water ice–ADH mixture by taking their geometric mean, using g H2 Os , the fraction of the H2 O(s)–ADH volume that is H2 O(s):

Fig. 2. Thermal conductivity of an ice–rock mixture as a function of temperature. The rock conductivity krock (lower solid line) conforms to measurements of ordinary chondrites (Yomogida and Matsui, 1983). The ice conductivity kice ( T ) (upper solid curve) assumes an ice that is mostly H2 O(s), with 1% NH3 by weight, or 3.1% ADH. Assuming a rock mass fraction of 0.63, we plot the total thermal conductivity (thick solid curve), calculated using the formulation of Sirono and Yamamoto (1997). This is well approximated as 3.49( T /100 K)−0.816 W m−1 K−1 (dash–dot curve). gH Os 1− gH2 Os

kice = kH O2 (s) kADH 2

(7)

.

With the thermal conductivities of the ices and rock found separately, we combine them to find the conductivity of an intimate mixture of ice and rock, using the formulation (motivated by percolation theory) of Sirono and Yamamoto (1997). Specifically, if V f rock is the volume fraction of rock in the combined ice–rock mixture, then ktot is found from:











V V 2k2tot − krock 3 f rock − 1 + kice 2 − 3 f rock ktot − krock kice = 0.

(8)

Fig. 2 shows the variation of the total conductivity of a mixture of 63% rock and 37% ice (by mass), with the ice being 1% NH3 , or 3.1% ADH by weight (which is potentially representative of Charon). As it happens, this conductivity is well described by a power law, ktot ≈ 3.49( T /100 K)−0.816 W m−1 K−1 . We calculate a thermal conductivity at 100 K (for a Charon-like composition with f rock = 0.63 and X = 0.01) of only 3.49 W m−1 K−1 in its surface layers. For a heat flux of 1 mW m−2 , this corresponds to a temperature change over 200 km of 57 K. Using the parameters and formulas of Castillo-Rogez et al. (2007a) for the same mixture of rock and ice, one would calculate a conductivity of 4.97 W m−1 K−1 ; over 200 km, the temperature change would be only 40 K. Finally, assuming a fully differentiated body with a pure water ice shell, the thermal conductivity in the outer layers would be 5.67 W m−1 K−1 , leading to a temperature jump of only 35 K over the same distance. Significant differences in temperature gradient can therefore arise depending on how the thermal conductivity is modeled. In particular, as this exercise illustrates, the inclusion of rock in an undifferentiated crust is insulating. 3.7. Other physical properties Besides thermal conductivities, there are other physical properties of materials to which attention should be paid. First are the densities. It is common when calculating the fraction of an icy body that is rock to adopt an ice density of 1000 kg m−3 (Gulbis et al., 2006; Hussmann et al., 2006; Schubert et al., 2007). If the ice contains contaminants such as CO2 , as Schubert et al. (2007) assume, the density may indeed be this high, but if not assuming this, a more appropriate ice density would be 918 kg m−3 , at least at the freezing point of ice. The density increases slightly with

700

S.J. Desch et al. / Icarus 202 (2009) 694–714

colder temperatures, and at 100 K, appropriate for KBO interiors, the density is 935 kg m−3 (Croft et al., 1988), which is the value we adopt here. For ammonia dihydrate (ADH), we adopt a density of 965 kg m−3 (Croft et al., 1988). For the density of rock, we consider chondrites to be a reasonable analog. Consolmagno and Britt (1998) and Wilkison and Robinson (2000) measured the densities of ordinary chondrites and found them to vary with porosity, ranging from ≈3600 kg m−3 for chondrites with near zero porosity, to ≈3200 kg m−3 for chondrites with ≈ 10% porosity. Yomogida and Matsui (1983) find the same trends; from their tabulations of the properties of L chondrites, we adopt a typical porosity of 10%, an intrinsic density 3600 kg m−3 and an average bulk density ρrock = 3250 kg m−3 . The latter value is more consistent with the thermal conductivities used above, since these are based on the average conductivities measured by Yomogida and Matsui (1983). If the ammonia content X and average density ρ¯ of a KBO are known, it is straightforward to use that information and the densities above to determine its rock (mass) fraction:

ρrock ρ¯ − ρice f rock = . ρ¯ ρrock − ρice

(9)

For example, taking Charon’s mean density, ρ¯ = 1710 ± 80 kg m−3 (Sicardy et al., 2006), and densities of rock ρrock = 3250 kg m−3 and ice ρice = 936 kg m−3 (appropriate for X = 0.01), we estimate Charon’s rock fraction to be f rock = 0.636 ± 0.036 (ignoring self-compression). For comparison, Person et al. (2006) find a rock fraction 0.58 ± 0.04 assuming ice and rock densities of 1000 and 3000 kg m−3 . For the heat capacity of H2 O(s), we use the approximation



c H2 O(s) ( T ) = 773.0



T 100 K

J kg−1 K−1 .

(10)

Over the range 40 < T < 273 K, this approximation matches the exact expression of Shulman (2004) to within a maximum deviation of 15% (at 75 K), and to within <1% above 200 K. For the heat capacity of ADH, we use the approximation

 c ADH(s) ( T ) = 1120

T 100 K



J kg−1 K−1 .

(11)

This matches the exact results of Fortes et al. (2003) to within 4% throughout the temperature range 40 < T < 176 K. For liquid water we set c P = 4188.5 J kg−1 K−1 , and adopt c P = 4700 J kg−1 K−1 for ammonia. Finally, the heat capacity of rock is a significant input, since the core will store significant amounts of heat over geological history. These heat capacities can be inferred from the data presented by Yomogida and Matsui (1983) for ordinary chondrites. We approximate the heat capacity of rock as

 T  ⎧ −1 −1 0  T < 275 K, ⎨ 770 275 K  J kg  K , −1 −1 c rock ( T ) = 607 + 163 T J kg K , 275  T < 1000 K, (12) 275 K ⎩ −1 −1 1200 J kg K , 1000 K  T . We note that the implied heat capacity is of the same order as common minerals; for example, our assumed value at 1000 K, 1200 J kg−1 K−1 , is comparable to the values c (1000 K) = 1250 J kg−1 K−1 for forsterite and 1210 J kg−1 K−1 for enstatite (Navrotsky, 1995). Heat capacities in mixtures combine arithmetically, weighted by the masses of the components. The expansivity of ice, α , is important for determining whether it can convect. It is common to assume a temperature-independent value for α that is appropriate to 273 K, e.g., α = 1.56 × 10−4 K−1 (Hussmann et al., 2006). In fact, α drops nearly proportionally with temperature, falling below 10−5 K−1 below 100 K, and actually changing sign below 60 K (Röttger et al., 1994). This means that at the frigid surface temperatures <60 K of many KBOs, ice cannot undergo solid-state convection, regardless of its viscosity, because

Fig. 3. Our assumed approximation to the phase diagram of the ammonia–water system. Region 1 corresponds to solid H2 O ice and ammonia dihydrate (ADH). ADH is assumed to melt in region 2 over a temperature range 174–178 K, producing H2 O and NH3 liquids. In region 3, solid H2 O ice gradually melts, in equilibrium with liquid H2 O and NH3 . At low values of the mass fraction of NH3 ( X < 0.0466), H2 O ice is assumed to melt in region 4 over a temperature range 271–275 K. In region 5, only liquid H2 O and NH3 are stable. Our assumed liquidus temperature T liq ( X ) is displayed in bold, and compares favorably with an exact function from Croft et al. (1988), displayed as a thin curve.

ice at these temperatures actually contracts upon heating. We approximate this behavior as α /(10−5 K−1 ) = 1.5( T /50 K) − 2.0 for T > 50 K, and = −0.5( T /50 K) for T < 50 K. We also adopt the expansivity of ADH from Fortes et al. (2003), which we approximate as α = 2.81 × 10−5 ( T /176 K) K−1 . We combine expansivities arithmetically, weighted by the volume of each component. 3.8. Equation of state Within each zone, we track 5 phases: rock; water ice, H2 O(s); solid ammonia dihydrate, ADH(s); liquid water, H2 O(l); and liquid ammonia, NH3 (l). At the beginning of each timestep, the mass of rock, the total mass of H2 O (in water ice, liquid water and bound in ADH), and the total mass of NH3 (in liquid ammonia and bound in ADH) are specified in each zone, as well as the total combined internal energy E i of all phases. The fraction of the total mass that is rock is denoted f rock,i . Also specified in each cell is the fraction of the non-rock mass that is in ammonia [either as NH3 (l) or the ammonia portion of ADH(s)], X i . Regardless of phase transitions (e.g., ADH melting) within each zone, the internal energy E i and total mass M i are conserved, as are X i and f rock . Given M i , X i , f rock,i and E i at a particular timestep, we define the temperature T i in that zone to be that temperature such that E i is the energy required to raise the mix of rock and ices from 0 K to the temperature T i , including phase transitions. Computing the energy needed to raise the mixture to a given temperature is done using an equation of state. To make the problem tractable, we have adopted a simplified version of the phase diagram of the ammonia–water system (Fig. 3). We do not consider ammonia monohydrate and so we ignore the distinction between the peritectic and eutectic points. We only consider 0  X  X c , where X c = 0.3210 (the fraction of the weight of ADH that is bound NH3 ). We also ignore the effects of pressure on the melting point of water ice. While the melting point of pure H2 O ice can be reduced to only 251 K at pressures of 200 MPa (Hogenboom et al., 1997), the peak pressures inside a KBO the size of Charon are only about (3/8π )G M 2 / R ≈ 135 MPa (assuming an homogeneous body); at the base of the ice shell in a partially differentiated Charon we predict even lower pressures, ≈80 MPa. In the bodies we consider, the melting point of H2 O

Cryovolcanism on Kuiper belt objects

ice is depressed by less than 10 K, and we neglect this effect. Up to pressures of 300 MPa, there appears to be no variation in the melting point of ammonia dihydrate (Hogenboom et al., 1997; ´ Leliwa-Kopystynski et al., 2002). The phase transition from Ice I to Ice II likewise requires pressures of about 220 MPa (Leliwa´ Kopystynski et al., 2002), that are unlikely to be achieved on the bodies in consideration here. A description of the phase diagram can be found in Appendix B. This phase diagram is used to determine the fractions of nonrock mass that is in each phase [H2 O(s), ADH(s), H2 O(l) and NH3 (l)], and the total specific energy e ( T ) required to raise a water–ammonia mixture (characterized by X i ) from 0 K to a temperature T. Once this is determined, the energy required to raise the rock and ices mixture in zone i from 0 K to a temperature T is readily found:

T E ( T ) = f rock,i M i

c rock ( T  ) dT  + (1 − f rock,i ) M i e ( T ).

(13)

0 K

In practice, however, we are faced with the inverse problem: given f rock , X and E, we must find the temperature T i in shell, by finding the temperature at which E ( T ) = E i . We accomplish this rapidly using simple bracketing techniques (Press et al., 1992). Thus, for a given internal heat energy E, we find the temperature T ( E ), as well as the apportionment of non-rock material into ices and liquids, that are consistent with that internal energy and temperature. Volume changes are explicitly neglected in this treatment. 3.9. Differentiation In our code we allow for differentiation of rock, liquids and ices. At the start of each simulation, we assume that each icy body is an homogeneous mix of rock and ice; but we allow that these phases may separate as the temperature is increased. If the conditions for differentiation are met, our code allows the process to proceed practically instantaneously (i.e., within the 1000 year timestep). Two types of differentiation must be considered. The first is the separation of the rock from the other components, forming a rocky core. The second is the separation of liquid water/liquid ammonia/ammonia dihydrate from water ice. Each depends critically but in different ways on the viscosity. We assume that separation of the rock from the ice occurs via solid-state creep of ice around the denser rock, and only via this mechanism. We consider the rock to exhibit Stokes flow, with a velocity V =

1 (ρrock − ρice ) g D 2rock 18

η

,

(14)

where D rock and ρrock are the diameter and density of the rock, and g (∼0.30 m s−2 inside a typical KBO) is the gravitational acceleration. We find that differentiation usually takes several ×107 yr, so we demand that rocks slip one grid zone (≈2 km) on timescales <106 yr to be consistent with our numerical code. For rocks with diameters of 1 m, differentiation over 1 Myr therefore demands η < 6 × 1011 Pa s, and differentiation over 1 Gyr demands η < 6 × 1014 Pa s, both thresholds increasing as the square of the rock diameter. The viscosities required for differentiation are lower than can be typically achieved in pure water ice, which are ∼1014 Pa s just below the melting point, and significantly higher at colder temperatures. In a KBO made of pure water ice, differentiation probably would not occur until the ice melted. But, if the ice contains just a few percent NH3 by weight, differentiation can be expected to occur when the ADH within it melts, at about 176 K. Just below this melting point, the measured viscosities are about 1017 Pa s

701

for X ≈ 5–10% (Durham et al., 1993), and differentiation would take a time greatly exceeding the age of the Solar System. But as soon as ADH begins to melt, producing some H2 O and NH3 liquid, the viscosity drops to <1012 Pa s for X = 4.0% or more (Arakawa and Maeno, 1994). It is not known what the viscosity of water– ice mixtures is for lower concentrations of ammonia ( X ), but a casual glance at the data of Arakawa and Maeno (1994) suggests that η ≈ 1013 Pa s for X ∼ 1%. Thus differentiation would proceed even with X ∼ 1%, provided the rock component is predominantly in objects at least 10 m in diameter. We assume within our code that differentiation occurs instantaneously wherever the temperature exceeds 176 K (actually, 174 K, the temperature at which melt is first produced in our approximate phase diagram). The maximum radius within a KBO at which the temperature ever has exceeded 176 K we denote R diff . (In general, R diff starts small and grows with time as the KBO heats up, to an asymptotic value.) All the rock inside R diff is assumed to have differentiated into a rocky core. Of the remaining phases, water ice is less dense (≈935 kg m−3 at 100 K; Croft et al., 1988) than ADH (≈965 kg m−3 ; Croft et al., 1988) and the liquid water–ammonia mixture (>950 kg m−3 at the same temperatures; Croft et al., 1988). At temperatures above 176 K, the slushy ADH and liquid water–ammonia should fall relative to the less dense water ice, and we assume that all the ADH and liquid phases end up sandwiched between a rocky core and a pure water ice layer inside R diff . Beyond R diff , material is undifferentiated. After the rock separates from other components, it is expected that ice and liquid in the remaining material can quickly separate via convection. The Rayleigh criterion for convection is that the Rayleigh parameter (see Eq. (3)) exceed the critical value Rac ∼ 103 . Assuming typical parameters of α ∼ 3 × 10−5 K−1 , ρice ∼ 1000 kg m−3 , g ∼ 0.10 m s−2 , κ ∼ 10−6 m2 s−1 ,  T ∼ 100 K, and r ∼ 100 km, we find Ra ∼ 3 × 1020 Pa s/η . Since the viscosity of the remaining ice is by assumption <1012 Pa s at the time of separation, vigorous convection is easily initiated in the remaining ice. This convection is expected to overturn the remaining phases in a time ∼ D 2 /κ /Nu, where D ∼ 100 km is the thickness of the ice layer, κ ∼ 2 × 10−6 m2 s−1 is the thermal diffusivity, and Nu ∼ 4 is the initial Nusselt number. Thus, within 40 Myr, the remaining phases should separate according to density. In particular, water ice will rise to the top and the ammonia will be concentrated in the liquid layer at the base. Although ADH and the liquid phases also have different densities, we do not assume that they differentiate. ADH is denser than its melt, and could be expected to sink to the bottom of the liquid layer upon freezing, but then it would warm up, melt, and rise again. This keeps the ADH and liquids well mixed, so we retain them within a single system, a “slush” layer that is well mixed chemically and energetically. This also may lead to effective heat transport, which we account for by assuming a high effective thermal conductivity (see Section 3.6). In our code, differentiation is implemented as follows. At each timestep, R diff is determined. The code searches for zones inside R diff containing rock, starting at the innermost zone. The mass and internal energy of rock from each zone are then used to fill up a new grid, starting with the innermost zones. After all the rock inside R diff has been moved, the code searches for all slush components (liquids and ADH) inside R diff and moves their mass and internal energy to the new grid, building on the rocky core there. Last, the code searches for pure water ice inside R diff and moves it to the new grid, filling up the new grid up to a radius of R diff . The densities of each component are chosen so that the total volume of the system does not depend on chemical phases; in this way the total volume of the system is conserved during differentiation, and the new grid is filled to exactly R diff . Mass outside R diff is moved to the new grid without change. As a last step, the chemical com-

702

S.J. Desch et al. / Icarus 202 (2009) 694–714

position and internal energies in the slush layer are homogenized, to account for the vigorous mixing expected in these layers. More specifically, we assign each cell associated with the slush layers the same ammonia concentration X equal to the average value in the slush layers, and the same energy per mass. This last step strongly imposes isothermality within the slush layers. Finally, we also calculate the gravitational potential energy U g before and after each timestep:

R p 4π ξ 2 ρ (ξ )

U g = −G

M (ξ )

ξ

dξ,

(15)

0

where M (ξ ) is the mass enclosed within a sphere of radius ξ . If differentiation has redistributed mass and changed U g during a timestep, we deposit that energy difference uniformly throughout all shells interior to R diff . We find that this additional heating is generally unimportant in KBOs, amounting to an increase of only a few K within the core during differentiation. Differentiation is a potentially complex phenomenon. For example, if rock is much smaller than the size we have assumed (meter-sized), or if liquid can somehow be lost from the ice before rock can differentiate, then differentiation would not proceed as we have described. We assume differentiation occurs in the manner described above; it remains to be seen how robust our assumptions are. 4. Thermal evolution of Charon As a canonical case, we now apply this code to the thermal evolution of a body comparable to Pluto’s largest moon, Charon. Charon’s radius has been reported as R c = 606 ± 8 km (Gulbis et al., 2006) and R c = 603.6 ± 1.4 km (Sicardy et al., 2006); combined observations yield R c = 606.0 ± 1.5 km (Person et al., 2006). Charon’s mean density has been reported as ρ¯ = 1720 ± 150 kg m−3 (Gulbis et al., 2006) assuming its mass is 1.60 ± 0.12 × 1021 kg (Olkin et al., 2003), or ρ¯ = 1710 ± 80 kg m−3 , assuming a mass 1.58 ± 0.07 × 1021 kg (Sicardy et al., 2006). Person et al. (2006) revised this downward to ρ¯ = 1630 ± 70 kg m−3 by assuming a mass 1.52 ± 0.064 × 1021 kg (Buie et al., 2006). Rather than fix particular values of mass and radius, here we consider a body with radius of 600 km and a mean density 1700 kg m−3 , yielding a mass 1.538 × 1021 kg. To be conservative, we consider a low ammonia fraction, X = 0.01. From these choices we derive a rock mass fraction f rock = 0.63. The mass of rock inside our Charon analog is therefore inferred to be 0.97 × 1021 kg. Assuming a value S ≈ 5.66 pW kg−1 in rock, this yields a steady-state surface heat flux due to the decay of long-lived radionuclides of 1.216 mW m−2 at the present day. The surface temperature of Charon is another input that is not well known. Fits to the shape of the absorption feature at 1.65 μm due to crystalline water ice suggest temperatures of 43 and 53 K on the sub-Pluto and anti-Pluto hemispheres of Charon (Cook et al., 2007). We regard these as lower limits to the surface temperature, because the absorption occurs in ice that is of higher albedo (and therefore lower temperature) than the rest of the surface. The presence of a poorly conducting regolith can also lead to a steep temperature gradient just below the surface. The thermal inertias of icy satellite regoliths have been determined to be typically (kρ c )1/2 ≈ 30–70 J m−2 s−1/2 K−1 (Spencer, 1989; Spencer and Moore, 1992; Spencer et al., 1999). From this we infer a thermal conductivity k ≈ 10−3 W m−1 K−1 within the regolith. Coupled to a heat flux ∼1 mW m−2 , this implies that a temperature jump of 1 K can occur across a regolith only 1 m deep. We therefore adopt a surface temperature that is a few degrees higher than the 53 K observed in ice on the anti-Pluto hemisphere of Charon; specifically, we adopt 60 K.

Fig. 4. Temperatures within a Charon analog, at times t = 0 (dotted line), t = 1 Gyr (dotted curve), and (solid curves, from top to bottom) t = 2, 3, 4 Gyr and the present day, t = 4.56 Gyr.

4.1. Numerical results The evolution of our Charon analog proceeds as follows. After accreting cold, the planet warms steadily as long-lived radionuclides decay. After roughly 75 Myr of evolution, temperatures near the planet center increase to the critical value of 176 K at which the ADH within the ice begins to melt, dramatically lowering the viscosity and allowing for differentiation. A rocky core quickly separates out; by t = 80 Myr it incorporates all rock within a radius R diff = 420 km and has a radius of about 292 km (radii are rounded off to the nearest grid point; here the grid resolution is 2 km). Roughly half of the core’s eventual mass differentiates during this rapid initial stage. For the next roughly 700 Myr, the differentiated region and the rocky core continue to grow, but much more slowly. By t = 0.75 Gyr the core reaches its maximum size of 356 km, and R diff reaches its maximum value of about 516 km. At this time a slush layer extends from about 356 km to the base of the water ice layer at about 368 km, and the pure water ice layer extends from about 368 to 516 km. The temperatures within the body are illustrated in Fig. 4 for several different times. Key features to note from this figure are as follows. First, temperatures in the core rise steadily for the first 2 Gyr; energy from the decay of long-lived radionuclides drives a heat flux outward, but at a rate reduced from the steady-state value, as energy also is consumed in heating the rocky core. The maximum temperature in the core is 1418 K, reached at 2.07 Gyr. It is worth noting that the solidus temperature of chondritic compositions is 1375 K (Agee et al., 1995), so a small fraction of the core should actually melt. After 2.07 Gyr, the temperatures in the core drop and the heat flux is enhanced by the release of stored heat. The heat flux emanating from the core is carried efficiently through the “slush” layer between 356 and 396 km because this layer is assumed to be partially liquid and convective, and is assigned a very large effective thermal conductivity. Convection also carries heat effectively through the ice layer, but is still associated with a non-negligible temperature gradient. The Nusselt number in the ice layer exceeds unity, but not by much, rising to its maximum of 4.4 at t = 1 Gyr, and dropping to 2.9 at t = 2 Gyr, 1.7 at t = 3 Gyr, and 1.0 at t = 4 Gyr. Thereafter convection shuts off. The distribution of phases within each layer, at t = 2 Gyr, is seen in Fig. 5. The distributions at other times look essentially the same as soon as differentiation is initiated, except that the values of R core and R diff change slightly, and the extent of the slushy layer vs. the ice layer will change. At the present day, R core and R diff are

Cryovolcanism on Kuiper belt objects

Fig. 5. Distribution of phases within a Charon analog at time t = 2 Gyr. The lines represent rock (solid line), solid or liquid water (dashed lines), and ammonia dihydrate or liquid ammonia (dash–dot lines). A core of pure rock exists within 355 km. Between 355 and 362 km there is a layer of water/ammonia liquid, and from 362 to 515 km there is a layer of pure water ice. From 515 km to the surface at 600 km there is an undifferentiated crust of rock, water ice and ADH.

703

Fig. 7. Mass fraction X of ammonia in the liquid inside a Charon analog, as a function of time. The eutectic concentration, 0.321, is plotted for comparison.

Fig. 8. Pressures inside a Charon analog, as a function of radius within the body, at a time t = +2 Gyr, when differentiation is complete. The vertical dotted lines represent (from left to right) the boundaries between the rocky core, the liquid layer, the water ice layer, and the undifferentiated ice/rock crust. Fig. 6. Total mass of liquids within a Charon analog (solid curve), and that mass of liquids derived from melting of ADH (dashed curve), as a function of time. The first liquid arises from ADH suddenly melted at t ≈ 75 Myr during initial differentiation. Continued heating melts water, but it eventually refreezes. By t = +4.5 Gyr, the remaining liquid refreezes as ADH, until by t = 4.61 Gyr it is completely frozen.

the same as in Fig. 5, except that much of the water liquid in the slush layer has frozen to water ice. At the present day, the liquid layer is predicted to extend from 356 to about 362 km. In Fig. 6 we display the total mass of liquids as a function of time. The total mass of ice on this Charon analog is 0.57 × 1021 kg, and the total mass of ADH is 0.0177 × 1021 kg. After about 80 Myr has elapsed, the temperatures rise to 176 K and the ADH within about 420 km all melts, producing 0.0061 × 1021 kg of liquid. (Eventually, as R diff increases, the mass of ADH that will be melted is 0.011 × 1021 kg.) At first, this liquid has an ammonia concentration near the eutectic, as seen in Fig. 7. Then, over the next several 100 Myr, water ice melts, diluting the ammonia. The liquid content peaks at about 1 Gyr, at a mass 1.8 × 1019 kg (about 1% the mass of the KBO), and an ammonia fraction X ≈ 0.2. Thereafter the core continues to heat but the heat fluxes escaping the core decrease. The water liquid begins to refreeze, a process that continues for several Gyr. The ammonia concentration increases as the water freezes onto the bottom of the ice mantle. As it happens, in these simulations the ammonia concentration reaches the eu-

tectic and the last of the liquid begins to refreeze right at around t = 4.5 Gyr, which is the present day. The last of the liquid takes about 200 Myr to refreeze. This duration partially reflects the temperature range of 4 K over which we presumed refreezing to take place, but is probably of the right magnitude. The latent heat associated with the freezing of 1.1 × 1019 kg of ADH over 200 Myr is equivalent to a surface heat flux of 0.05 mW m−2 , or 3% of the total heat flux. If the refreezing took place over 10 Myr, it would have an appreciable impact on the internal thermal structure. Finally, we plot in Fig. 8 the pressure as a function of depth within the KBO at t = 4.56 Gyr. The pressure was calculated simply by integrating dP dr

=−

G M (r )ρ (r ) r2

(16)

from the surface, where P = 0. Here M (r ) is the mass enclosed within a shell at radius r. For this modeled KBO, the pressure at the center of the body is 269 MPa, falling to 78 MPa at the edge of the rocky core. The central pressure is nearly doubled over its value before the differentiation, 135 MPa, due to the concentration of mass at the center. The pressures in the liquid layer are roughly 78 MPa, and those in the ice layer vary from 78 to about 38 MPa. The pressures in the undifferentiated crust are lower still, <38 MPa. These low pressures justify the neglect of higher-order

704

S.J. Desch et al. / Icarus 202 (2009) 694–714

phases of ice and other pressure effects. These effects could become important, however, on larger icy bodies. Taking 200 MPa as the limit for which our water–ammonia phase diagram can be considered important, we estimate these pressures would be reached in the liquid layer for bodies larger than about 900 km in radius. Our models therefore cannot capture the complex physics in the large satellites (Europa, Ganymede, Callisto, Titan), but might be applied with some care to large KBOs (Pluto, Triton, Eris). All other KBOs and icy satellites, with radii <800 km, should be described decently by these models. 4.2. Analysis The smooth variation of temperature with radius seen in our numerical results suggests that an analytical solution may be obtained. We are motivated to find such a solution, both as a check on the numerical results, and also because it will provide physical insight into the behavior of the system. To do so, we assume the same compositional stratification as found by the numerical simulations, then integrate the heat conduction equation. We assume a heat flux that is the sum of two components: a flux driven by the decay of long-lived radionuclides, F rad ( R p ); and an additional flux arising from the release of stored heat from the rocky core, ≡ F rad ( R p ). Because stored heat could be released from other layers besides the core, this approach is necessarily approximate; nonetheless, it does provide useful insights. We first examine the undifferentiated crust layer. Although this layer contains radioactive rock that generates heat, we assume no flux arises from release of heat stored within this layer, so refers to heat released from the rocky core and the ice and liquid layers. The flux at radius r < R p is straightforwardly found to be F (r ) = F

rad

 2   Rp r . (Rp) + Rp

T ( R diff ) T0

 1− p

 =

(17)

r

T (Rp)

 1− p

T0

  ×

1 2

1−

+ (1 − p )

R 2diff



 +

R 2p

F rad ( R p ) R p k0 T 0 Rp R diff

 −1 .

(18)

We assume T ( R p ) = 60 K, a thermal conductivity in the undifferentiated crust of k( T ) = 3.49( T /100 K)−0.816 W m−1 K−1 (see Fig. 2 and Section 3.6), and a present-day radiogenic heat flux at the surface F rad ( R p ) = 1.216 mW m−2 . These quantities are basically parameters of the problem; the last two are taken from the numerical simulations, from which we find R diff = 515.3 km = 0.859R p , and observe a heat flux at the surface of 1.724 mW m−2 at t = 4.56 Gyr, implying = 0.418. With these parameters in hand, our analytical solution predicts T ( R diff ) = 100(0.961 + 0.0632 )5.435 K = 93.2, which compares extremely well with the observed value, T ( R diff ) = 93.3 K. We would also predict a heat flux at r = 517 km (near the bottom of the undifferentiated crust) equal to F = 1.732 mW m−2 , which is within 1% of the value we observe in the simulations, F = 1.747 mW m−2 . We next integrate inward from R diff to the base of the ice layer at R base . The flux in the ice layer is

 F (r ) = F rad ( R p )

R diff Rp

3

 +

 rad   3  R diff F (Rp)Rp Rp Rp × exp . (20) − +

3 k0 T 0

Rp r

2 .

(19)

Combining this with the conductivity of pure water ice, k( T ) = k0 T 0 / T , where k0 = 5.67 W m−1 K−1 and T 0 = 100 K (Section 3.6), we derive the temperature at the base of the ice layer:

R base

R diff

Rp

Here it is assumed that the temperatures are low enough in the ice layer to render it non-convective. Our numerical simulations show that at t = 4.56 Gyr, R base = 362.2 km, so this analytical solution would predict T ( R base ) = T ( R diff ) exp(0.401 + 0.633 ) = 181.5 K, which is within 2.5% of the value observed in the simulations, T ( R base ) = 177.1 K. Likewise, we would predict a heat flux at r = 367 km (near the base of the ice layer) of F = 3.417 mW m−2 , which matches to within 2% of the observed flux, F = 3.362 mW m−2 . Our simulations show that the heat flux through the undifferentiated crust is augmented by 42% over the value arising from instantaneous radioactive decay alone. The additional heat flux is completely consistent with that expected from the release of stored heat from the rocky core, as we now demonstrate. Within the rocky core, assuming a spatially uniform rate of radiogenic heating and a constant thermal conductivity, the temperature in the core is easily shown to conform to





T (r ) = T (0) − T (0) − T ( R core ) x2 ,

(21)

where R core is the outer boundary of the core and x = r / R core . The temperature at the outer boundary ≈200 K, but the central temperature is much higher. The heat energy of the core is then easily found by integrating the internal energy (Eq. (13)) over the volume of the core: R

core





1

4π r 2 ρrock e rock T (r ) dr = 3M core

E core = 0

Assuming a power-law thermal conductivity in this rock–ice mixture, k( T ) = k0 ( T / T 0 )− p , the heat conduction equation is easily integrated to find the temperature at r = R diff :



T ( R base ) = T ( R diff )





x2 e rock T (x) dx. 0

(22) To proceed we must integrate the heat capacity of rock (Section 3.7), and use temperatures computed in the simulations. At t = 4 Gyr, T (0) = 1077 K and T ( R core ) = 197 K; at t = 4.56 Gyr, T (0) = 932 K and T ( R core ) = 176 K. Between these two times, it is tedious but straightforward to integrate the heat capacity of rock over the core to show that the core loses an amount of heat energy ≈6.85 × 104 J kg−1 . The core mass is 6.15 × 1020 kg, so the average heat loss over the last 0.56 Gyr of the simulation is equivalent to ≈ 0.43 times the present-day radiogenic heating. It is possible to quantify other sources of stored heat. From about t = 2.0Gyr to 4.5 Gyr, about 4 × 1018 kg of water freezes, releasing latent heat of freezing equivalent to = 0.003. The heat flux from freezing of water is therefore negligible. More significant is the freezing of the last remaining liquid from t = 4.51Gyr to 4.68 Gyr, which releases a burst of latent heat equivalent to = 0.049. The other sources of stored heat, cooling of the ice, liquid, and outer layers, are even less important. So it is the release of stored heat from essentially just the rocky core that is increasing the heat flux through the KBO. The heat stored in the core was generated during the earliest phases of the KBO’s evolution. While it was being stored, the heat escaping to the surface necessarily was lower than it would otherwise be; later, when the heat from the core was released, the heat flux is enhanced over the amount due to instantaneous radioactive decay alone. This is illustrated in Fig. 9, which shows that the point at which the core begins to release heat it has already stored is reached at about t = 1.6 Gyr. At the present day (t = +4.56 Gyr), the heat flux at the surface is enhanced over the instantaneous radiogenic heat flux by a factor of 1.42, or = 0.42. This is a significant increase in the heat flux, and demonstrates the importance of the KBO’s past history in determining its current thermal structure. In fact, without these stores of heat (setting

Cryovolcanism on Kuiper belt objects

Fig. 9. Heat flux at the surface of a Charon analog, as a function of time (solid curve). For comparison we plot the flux that would arise from the instantaneous rate of radioactive decay (dashed line). During the first 1.6 Gyr, the rocky core absorbs more heat than it emits, but thereafter the core releases heat and enhances the heat flux.

= 0), the temperature at R diff would be only 81 K, and the temperature at R base would be only 120 K, and no liquid layer would be possible. The value of is a function of the prior heating of the core. Thus, it is impossible to judge whether or not subsurface liquid can persist within KBOs without performing a time-dependent calculation. 5. Sensitivity to various parameters 5.1. Ammonia Ammonia abundance has a large bearing on the thermal evolution of a KBO, so we consider two departures from our canonical abundance X = 1%. The first case involves an increase in ammonia, to X = 5%, near the upper limit of what we deemed plausible in the ices accreted by KBOs (Section 2). Surprisingly, there is very little difference between this case and the canonical case. The additional ADH in the undifferentiated crust reduces the thermal conductivity somewhat, and allows differentiation to proceed out to roughly 528 km instead of 516 km, but this is a relatively minor effect. The amount of liquid that exists within R diff is, at all times, larger than in the canonical case, because there is more ammonia available as antifreeze; but the profile of M liq versus time looks identical to Fig. 6, merely scaled upward by a factor of exactly 5. Significantly, the time at which the liquid freezes is not changed much. The extra latent heat released upon freezing of the liquid buffers the freezing slightly more effectively than in the canonical case. Total freezing does not occur until 4.57 Gyr, instead of 4.54 Gyr, an insignificant difference (both simulations using 100 zones). The amount of liquid varies proportionally to X , but the behavior of the system is very insensitive to the value of X , provided it is large enough (a few percent) to allow differentiation at 176 K. It is especially instructive to compare the canonical results obtained above to the case where ammonia is absent. The presence of ammonia with X exceeding a few percent was found sufficient to trigger a change in viscosity at the melting point of ammonia dihydrate, 176 K. Here we consider a KBO with R p = 600 km, ρ = 1700 kg m−3 , but with X = 0. In this scenario without ammonia, the viscosity of the ice remains high at all temperatures. In ice just below the melting point, the assumed viscosity, 1014 Pa s, is high enough to prevent the Stokes flow of rocks <1 m in diameter from flowing faster than ∼10 km Gyr−1 , an inconsequential amount. We therefore consider differentiation to take place only

705

Fig. 10. Same as Fig. 6, but for the case X = 0. By t = 2.81 Gyr all the liquid is completely frozen.

when temperatures exceed 271 K, the point at which ice begins to melt in our code. Because of the lack of ammonia, no melt is produced until about 200 Myr, much later than the time taken to produce melt in our canonical case, roughly 70 Myr. When melt is produced and differentiation occurs, it ultimately is limited to R diff = 455 km, compared to 515.3 km in the canonical case, and R core = 313.2 km, compared to 354.5 km in the canonical case. The mass of the rocky core that forms in the absence of ammonia is only 70% as large as the core in the canonical case. The smaller size of the core leads to smaller temperatures within it. Temperatures at r = 0 peak at t = 1.94 Gyr at 1210.3 K, rather than at 1418.0 K at t = 2.06 Gyr. Instead of storing an average 4.90 × 105 J kg−1 in the core, the rock core in the X = 0 case stores only about 4.52 × 105 J kg−1 . The rocky core in the X = 0 case stores only 2/3 the heat of the X = 1% case (at t = +2 Gyr), about 1.9 × 1026 J instead of 3.0 × 1026 J. The lack of a substantial heat reservoir allows considerably more heat to be released into the outer layers early on, but considerably less heat later. In contrast to the X = 1% case, the rocky core absorbs more heat than it emits only until t ≈ 0.5 Gyr (as opposed to 1.6 Gyr). The difference in heat energy that is released early on, an additional 1.1 × 1026 J, is sufficient to melt roughly 3.3 × 1020 kg of water ice. producing an order of magnitude more liquid than in the X = 1% case. As illustrated in Fig. 10, the mass of liquid (all water) peaks between 1.0 and 1.5 Gyr, at about 2.3 × 1020 kg. Thereafter the liquid gradually refreezes. The transition to a completely frozen state is buffered by the immense amount of liquid that must refreeze, and the latent heat of fusion released as it does. The heat released by the refreezing of 2 × 1020 kg of water over 1 Gyr is equivalent to a surface heat flux of 0.5 mW m−2 , which alone is a significant fraction (≈10%) of the total heat flux escaping the body at this time. Essentially, the liquid freezes only as fast as heat can be released; this is in contrast to the X = 1% case, where the freezing of the liquid had less bearing on the overall heat flux. The liquid is completely frozen by t = +2.81 Gyr, and thereafter the heat flux drops considerably. 5.2. Thermal conductivity We have performed an additional simulation of our canonical case, with R p = 600 km, ρ = 1700 kg m−3 and X = 1%, but with the rock conductivity doubled to krock = 2 W m−1 K−1 . This slightly increases the conductivity in the undifferentiated crust, and greatly increases it in the core. The increased conductivity thus has two effects, both of which decrease the internal temperatures. For a fixed

706

S.J. Desch et al. / Icarus 202 (2009) 694–714

flux, the increased conductivity reduces the temperature gradient through the ice crust shell, and lowers the temperatures throughout the body. But an increased rock thermal conductivity reduces the fluxes, too, because less heat can be stored in the core while its temperature is increasing in the first 2 Gyr. This manifests itself at the present day via a reduced . Both these effects have the potential to depress the temperatures and limit the time for which liquid can exist. A quick comparison with the k = 1 W m−1 K−1 serves to illustrate the differences. In the increased conductivity case, at t = +4.56 Gyr, R diff = 515.0 km, and R core = 354.5 km, so these quantities are unchanged. On the other hand, the last liquid on this KBO freezes after 4.05 Gyr, instead of after 4.67 Gyr, as in the canonical case. The temperatures are generally lower: at the base of the undifferentiated crust, T ( R diff ) = 84.1 K, instead of 93.3 K, and at the base of the ice layer T ( R base ) ≈ 143 K, instead of 177 K. The reason is the lower heat stored in the core. At t = +4.56 Gyr, again, the temperature at r = 0 is 434 K, rather than 942 K. Less heat energy in the core translates into less heat flux escaping the core. Our results are consistent with a flux enhanced over the steadystate radiogenic heating rate, but with = 0.32 rather than 0.42. That is, the heat being released from the core at t = +4.56 Gyr is reduced below the canonical case by about 25%. The case described above shows that increasing the conductivity of rock by a factor of 2 has a significant quantitative effect, accelerating the freezing of the body by about 0.6 Gyr. We have also run a simulation with the rock conductivity reduced by a factor of 2, to 0.5 W m−1 K−1 . As might be expected, the lower conductivity leads to higher internal temperatures and delays the freezing of the body, by about 0.3 Gyr. The thermal conductivity of rock is not well known and variations of a factor of 2 may be considered reasonable; while we favor the measured value for cold ordinary chondrites, krock ≈ 1 W m−1 K−1 (Yomogida and Matsui, 1983), we do note that terrestrial rocks tend to have higher conductivity, krock ≈ 2 W m−1 K−1 . The uncertainty in the conductivity leads to an uncertainty of about 0.4 Gyr in the time at which the liquid in the body completely freezes. The evolution of a KBO with different thermal conductivity is not, however, qualitatively different from the canonical case we have considered. 5.3. Size and density Increasing the size or density of a KBO allows it to retain subsurface liquid longer. A higher mean density implies a higher rockto-ice fraction, which enhances the radiogenic heating at all times, and also provides a larger core that can store heat; both effects prolong the time for which some ice can remain liquid. Likewise, a larger object will contain a greater mass of rock. By inspection of Eqs. (18) and (20), the planet radius enters into the temperature most strongly as the product F rad ( R p ) R p . All other things being equal, increasing the planet radius by a factor of 2 should allow liquid to be retained with only half the surface flux, which is achieved on this larger body after an additional 4 Gyr of decay. Equivalently, an increase in radius by 10 km might be expected to extend the duration for which liquid exists by about 100 Myr, not accounting for changes in R diff , R base and due to heat stored in the rocky core. In Fig. 11, we display the results of 40 simulations we have performed, with ρ = 1300, 1500, 1700, 1900 and 2100 kg m−3 , and R p = 300, 400, 500, 600, 700, 800, 900 and 1000 km. An ammonia abundance of X = 1% has been assumed. Plotted are contours of the time at which all liquid in the KBO freezes. We note that these simulations were carried out at a lower numerical resolution (100 zones) than our canonical simulation (300 zones). Direct comparison between the cases with R p = 600 km and ρ = 1700 kg m−3 shows that the poorer numerical resolution accelerates the freez-

Fig. 11. Times after formation at which all liquid freezes within a KBO with initial ammonia content X = 1%. KBOs with mean densities and radii that place them to the right of the rightmost curve potentially can retain subsurface liquid to the present day.

ing by 0.13 Gyr (from t = 4.67 Gyr to t = 4.54 Gyr). The freezeout times are probably underestimated uniformly by ∼0.1 Gyr in this plot. The numerical simulations largely confirm that liquid can persist for many Gyr, to the present day in fact, in objects the size and density of Charon. The boundary between objects having liquid at the present day and those that do not is given roughly by ρ¯ R p ≈ 1.0 × 109 kg m−2 . This suggests that the heat flux at the surface (∝ f rock M p / R 2p ∝ f rock ρ¯ R p ) is the most important quantity controlling liquid content. The simulations demonstrate that for bodies with smaller density or radius, the time at which complete freezing occurred is roughly proportional to ρ¯ R p :



t freeze ≈ 4.6

ρ¯ 1700 kg m−3



Rp 600 km



Gyr.

(23)

The simulations also verify that every 10 km increase (decrease) in radius can prolong (shorten) the duration for which liquid can persist by about 100 Myr. We caution that this time of complete freezing is uncertain by many ×108 yr, as described below. These conclusions do not change substantially when the ammonia abundance is increased to X = 5%. We have run additional simulations with identical radii and densities as in Fig. 11, but with X = 5%. There is a slight tendency for the more ammonia-rich bodies to retain liquid for a longer time, but the difference amounts to only about 0.1 Gyr. The later freezeout time is due in part to a slight increase in R diff , by about 10 km, a consequence of a slightly lower thermal conductivity in the ice, now richer in ADH. An increase to X = 5% will increase the total amount of water–ammonia liquid there is in a KBO at any time, roughly as X ; it will also slightly increase the amount of water liquid. The later freezeout time is also due in part to the greater release of latent heat of freezing of this greater mass of liquid. Finally, we note that Charon’s position in these plots (R = 603.6 km, ρ = 1710 kg m−3 ) implies that it is just now reaching the point where its last liquid is about to freeze; however, it is difficult to assign too much certainty to this result. Several sources of systematic errors (e.g., in how convection is treated) and other uncertainties (e.g., in thermal conductivities) are all seen to shift the freeze-out times in Fig. 11 by several ×108 yr (even if the KBO radius is known accurately). Also, to generate results rapidly, the models used to create Fig. 11 used only 100 radial zones; however, convergence is found to require a numerical resolution of at least 300 zones, which would increase the freeze-out times by roughly 0.1 Gyr. The crudeness of the calculations preclude us from saying

Cryovolcanism on Kuiper belt objects

Fig. 12. Times after formation at which all liquid freezes within a KBO with initial ammonia content X = 5%, in the cases where convection is treated as described in the text (solid curves) and where convection is shut off (dashed curves). The differences in freezeout times are minor, and are comparable to the variations in freezeout time due to uncertainties in the inputs.

when the fluid within Charon should freeze entirely, other than sometime between 4 and 5 Gyr after formation (today, plus or minus 0.5 Gyr). 5.4. Treatment of convection Our treatment of convection is crude: we do not treat convection in the stagnant lid regime, and there are uncertainties in ice rheology. To quantify the effects of uncertainties due to how convection is treated, we have also performed a series of simulations with X = 5% (and 100 grid zones), but with convection turned off. Turning convection off effectively decreases the conductivity early on, allowing the heat to be stored more effectively in the core. KBOs that do not convect should therefore retain liquid at later times than the KBOs considered here. The results of our simulations with and without convection are shown in Fig. 12, and confirm our expectations. For example, a KBO with radius 600 km and density 1700 kg m−3 is seen to freeze at t ≈ 4.5 Gyr in the case with convection, and at a later time, t ≈ 4.7 Gyr, when convection is shut off. In all cases, however, the differences in freeze-out time are small, ≈0.1 Gyr, and are comparable to the variations in freeze-out time due to uncertainties in the inputs to the calculation. This result indicates that the details of how convection is treated, and the exact viscosity used, are relatively unimportant for KBOs. 6. Discussion 6.1. Summary of findings We have presented models of how KBOs thermally evolve over time due to the heating by long-lived radionuclides. A key feature of our models is the inclusion of ammonia. Ammonia plays three important roles, all of which aid in the retention of liquid. Obviously, ammonia is an antifreeze, depressing the melting point of water (at least at the eutectic) by 97 K. Ammonia is the reason liquid can exist even when the temperatures are near 176 K. More subtly, ammonia also has the effect of lowering the thermal conductivity of its host ice to a small degree. In the models presented here, however, the major role of ammonia is its effect on the viscosity. The presence of a few percent ammonia in water ice drastically lowers the viscosity by orders of magnitude once

707

the first melt is produced at 176 K. The drastically lowered viscosity enables differentiation of rock and ice, allowing the KBO to evolve to a new density configuration. This configuration allows liquid to exist just outside a rocky core but under a thermally insulating crust, and allows ammonia to be concentrated as the liquid refreezes. As demonstrated by the simulations without ammonia (Section 5.1), wherein differentiation takes place only at the melting point of water, retention of liquid on KBOs today is essentially not possible without at least ∼1% ammonia (by weight, relative to water). Another key feature of our simulations is the allowance for partial differentiation. Because differentiation in our model requires the production of some water–ammonia melt, at 176 K, differentiation does not reach the surface of the KBO, which is always much colder (≈30–60 K). For the typical temperature gradients in the crusts of KBOs (≈1 K km−1 ), the maximum radius of differentiation is likely to lie under about 100 km of rock/ice crust. The role of this crust is significant, because a rock/ice mixture is much less conducting than a pure ice layer (see Section 3.6). At depths of a few 100 km, the insulating properties of this crust raise the temperature by many tens of K typically. The existence of a high-density (∼1700 kg m−3 ) rock/ice crust overlying an ice layer of density <1000 kg m−3 is, of course, gravitationally unstable. To the extent that the crust can flow as a fluid, it seeks to sink under the ocean. The key to the stability of this configuration is that the crust is rigid and viscous and unable to flow. Throughout the undifferentiated crust, temperatures do not ever exceed 176 K (by definition), and the viscosity of ice is no less than 1017 –1018 Pa s. Stokes flow of rock through this ice occurs at velocities great enough to be relevant (e.g., a few km per Myr) only if the rock diameters exceed at least 300 m (cf. Section 3.9), Especially in light of the fact that ice and rock are likely to co-accrete into KBOs, it is difficult to see how kilometer-sized regions of pure rock would arise. Likewise, following the logic discussed in Section 3.9, this layer cannot initiate convection. Despite the technically unstable nature of the 100-km-thick rock–ice crusts at the surfaces of our model KBOs, we accept that they might remain in this state because we see no way for the layers to overturn themselves at the cold temperatures of KBOs. Finally, the third key feature of our simulations is that they are time-dependent, and can account for the storage of heat at early times when radiogenic heating is greater, and release of that heat later as radiogenic heating wanes. This effect is often neglected, but we find that for Charon-sized KBOs the release of stored heat is a major factor. The surface heat flux is increased by 42% in our canonical model, mostly by heat stored in the core and released at t > 2 Gyr, after the core has reached its maximum temperature. The combination of these effects is that KBOs and icy satellites can retain liquid water to the present day if they are as small as 600 km in radius or roughly 1.5 × 1021 kg in mass, provided they accreted with some ammonia ( X ∼ 1% or greater). The exact time for which liquid can be retained depends most strongly on the treatment of convection and on the thermal conductivity of rock, and is uncertain to within a few ×108 yr. 6.2. Implications for cryovolcanism A significant implication of our findings on the thermal evolution of KBOs is that they offer the possibility of bringing subsurface liquid to the surfaces of KBOs. On bodies the size of Charon and larger, subsurface liquid can exist, provided it is mixed with ammonia. We hypothesize that this liquid can be brought to the surface via self-propagating cracks that open in the brittle ice lying between the surface and the subsurface ocean. Crawford and Stevenson (1988) developed this model in the context of Europa. Using linear elastic fracture mechanics, they showed that for cracks

708

S.J. Desch et al. / Icarus 202 (2009) 694–714

of length l propagating from the bottom of the ice layer (from the top of the ocean), the stress at the tip of the crack is





K = (π l)1/2 T − 2g (ρ )l/π ,

(24)

where T is the tidal stress that they considered to initiate the crack, and ρ is the density of the liquid, minus the density of the ice or medium through which it is propagating. The liquid in the crack transmits pressure to the tip of the crack, which tends to widen it. If the stress exceeds the critical value K c = 6 × 106 N m−3/2 for water ice, Crawford and Stevenson (1988) showed, the crack will continue to propagate. On Europa, the appropriate densities are 920 kg m−3 for the ice and 1000 kg m−3 for the liquid (essentially just water), leading to ρ = +80 kg m−3 . Positive K therefore requires T to be large and positive; tidal stresses are required for cracks to propagate in the europan ice shell. On KBOs the situation is very different. The density of the ice is ≈930 kg m−3 , and the density of the liquid, being ammonia-rich, is likely to be much lower than the density of water alone. The density depends on temperature: a liquid with X = 0.32 has a density 880 kg m−3 at T = 288 K, rising with decreasing temperature to the density of ice at 230 K, and about 945 kg m−3 at T = 176 K (Croft et al., 1988). The density of the liquid is therefore less than the density of the ice if the liquid is warmer than about 230 K, and is never more than 15 kg m−3 denser. In our Charon analog, the liquid reaches temperatures ≈230 K at about t = +2 Gyr, during which time ρ < 0. Thus, even without tidal stresses (T = 0), K > K c and cracks would self-propagate, provided l starts with a large enough value. For the typical parameters being discussed here, that minimum length is l ∼ 1 km. If a crack is initiated in the ice at the bottom of the ice layer, and the length is at least 1 km, it potentially will propagate all the way through the ice layer to the undifferentiated crust. There the density contrast is even greater and the crack potentially can reach the surface. An additional difference between KBOs and Europa is that the base of the ice layer in a KBO, being at roughly 200 K, is far less ductile than the base of an ice layer on Europa. Initiation of a crack is likely due to stresses associated with the expansion of ice upon freezing. It is worth noting that despite the low density of an ammonia/water liquid, the liquid must always expand as the water within it freezes into water ice (this no longer holds once the eutectic composition is reached and the liquid freezes directly to ADH). A trapped pocket of liquid surrounded by ice will, upon freezing, generate enormous overpressures as it freezes, because of the incompressibility of the liquid; as Fagents (2003) have demonstrated, these overpressures are alone sufficient to drive liquid to the surface given a conduit. We view such pockets of trapped liquid likely locales for cracks to be initiated, if they exist. The volumes of liquid involved, and the needed size of the initial crack, demand that pockets of trapped liquid would need to be many kilometers in size. This mode of initiating cracks would serve to counteract the slight overdensity of the liquid that occurs for t > 2 Gyr, when the liquid is colder than 230 K. Once a crack is initiated, it will propagate to the surface at several meters per second (Crawford and Stevenson, 1988). Cracks could therefore traverse the 240 km to the surface, in our canonical case, in about 105 s, or about 1 day. It is worth noting that during this time, relatively little cooling of the liquid is expected. Taking a temperature of 200 K, we estimate the thermal diffusivity of a water–ammonia mixture to be no higher than 10−6 m2 s−1 , so that the timescale for losing heat out of a crack with radius R ∼ 1 m would be ∼106 s, about 10 times longer than the time taken to travel along the crack to the surface. While some of the water may freeze out of the water–ammonia liquid during its passage (raising the ammonia content of the liquid), for reasonable parameters this effect would seem to be limited. In fact, it

seems more probable that the fluid would entrain xenoliths of native “rock”, mostly water ice crystals. Such a possibility has been considered in the case of cryovolcanism on Titan and Enceladus, and xenolith weight fractions of several tens of percent are considered plausible, by analogy to terrestrial volcanos (Fortes, 2007; Fortes et al., 2007). This would add water ice to the surface. In any case, it seems likely that the liquid reaching the surface would be ammonia-rich, but with some water ice entrained in it. This material will then spill over the surface; or, if it contains dissolved gases, might more violently explode upon reaching the surface. The globally averaged rate at which liquid is brought to the surface is likely to be the rate at which the volume available to the liquid decreases, due to the expansion of the ice as it freezes. As discussed by Cook et al. (2007), the volume of ammonia–water mixtures with X < 0.321 must increase as they freeze, by an amount comparable to 7%, with little dependence on X . In our canonical simulation, between about t = 1 Gyr and the present day, about 7 × 1018 kg of liquid freezes, implying the displacement of 7% of that liquid, or 5 × 1017 kg, over 3.5 Gyr. Provided selfpropagating cracks provide a conduit for this liquid to reach the surface, it is energetically favored to do so and will resurface our Charon analog to a depth of 120 m, at a fairly uniform (global) rate of 0.033 km Gyr−1 . Freezing of ADH does not involve an expansion, and it seems unlikely that cryovolcanism could occur during this stage, nor, of course, once all the liquid is frozen. The rate at which liquid is predicted to be brought to the surface is consistent with both the presence of crystalline water ice and trace amounts of ammonia hydrates observed on Charon’s surface. First we consider the ammonia hydrates. Since the spectral signatures of ammonia hydrates on the surface are destroyed by GCRs (Section 1), ammonia hydrates are destroyed to depths ≈350 μm on 20 Myr timescales. For ammonia hydrates to be seen on the surface, they must be emplaced at a globally averaged rate of at least 2 × 10−5 km Gyr−1 . As the liquid brought to the surface is near the eutectic, the rate at which ammonia dihydrate melt is brought to the surface is less than but comparable to 0.033 km Gyr−1 . Ammonia hydrates can be emplaced on the surface about three orders of magnitude faster than they can be destroyed by GCRs, while cryovolcanism is occurring. Cook et al. (2007) likewise estimated that to preserve the spectral signatures of crystalline water ice, new crystalline water ice would have to be emplaced at a globally averaged rate equal to the penetration depth of the near-infrared (NIR) observations, ≈350 μm, divided by the amorphization timescale of 3 × 104 yr (if due to UV amorphization) to 1.5 Myr (by GCRs). Emplacement of crystallization water ice would outpace amorphization, then, if the globally averaged deposition rate of water exceeded 0.012 km Gyr−1 (UV) or 2.3 × 10−4 km Gyr−1 (GCRs). Just before cryovolcanism ends, the fraction of liquid brought to the surface that will freeze as water ice is small. For example, in our canonical run, 0.5 Gyr before cryovolcanism ends, X ≈ 0.28 in the liquid. Thus the fraction of liquid that can freeze as water ice is 1 − X /0.321 = 13%. To this several tens of percent more of the liquid mass may be present as water ice crystal xenoliths. At the present day on Charon, the rate at which water ice is emplaced on the surface is therefore on the order of several tens of percent of the total mass of cryomagma. We therefore estimate the rate of deposition of water ice on the surface to be about ∼0.01 km Gyr−1 near the time cryovolcanism ends, if xenoliths are significantly entrained in the flows. These rates are just adequate to compete with the rate of UV amorphization estimated by Cook et al. (2007), and are more than adequate to compete with GCR amorphization. The rates of emplacement of ADH and crystalline water ice on the surface will scale upward if X exceeds the conservatively low value of 1% assumed in the canonical case. Thus sufficient quantities of crystalline water ice and ammonia hydrates can be emplaced on

Cryovolcanism on Kuiper belt objects

the surface to account for the spectroscopic observations, even in the face of destruction by GCRs and UV. The cryomagma brought to the surface is necessarily ammoniarich. It cannot be at the eutectic composition or else it would not expand upon freezing and be extruded; but it does have X > 0.2 throughout our canonical simulation. Indeed, if the liquid becomes too water-rich it will lack the buoyancy possibly needed to drive a self-propagating crack through the ice layer. If all of the water ice and all of the ammonia dihydrate on the surface of Charon were emplaced by cryovolcanism, we should predict a high abundance of ADH in the spectrum. Using the estimates above, if water ice is emplaced at ∼0.01 km Gyr−1 and ADH at a rate ∼0.03 km Gyr−1 , then we would predict X ≈ 0.24, significantly higher than the value X ≈ 1% inferred from spectroscopic measurements (Cook et al., 2007). This discrepancy implies that water ice is more widespread than ADH on the surface. One explanation is that the spectroscopic measurements are including significant amounts of ancient, amorphized water ice as crystalline. Laboratory analysis of the spectra of ices shows that ice that is a mixture of 20% crystalline water ice and 80% amorphous water ice is difficult to distinguish from 100% crystalline water ice (R. Mastrapa, 2008, personal communication,). So if cryovolcanism only accounted for 20% of the ice observed on the surface, with ancient, amorphous water ice accounting for the other 80%. This would reduce our prediction of X on the surface by a factor of 5, to X ≈ 5%. Another factor which might ameliorate the discrepancy is if water ice simply spreads out more than ammonia-rich liquid upon reaching the surface. If the cryovolcanism is explosive, entrained water ice xenoliths may be launched in a plume around a vent, while ADH might ooze out. Of course, a third interpretation is that crystalline water ice is not exclusively produced by cryovolcanism (e.g., it is due to micrometeorite impacts), in which case only ADH on the surface would need to be explained by cryovolcanism. Cryovolcanism covering only a few percent of Charon’s surface would then be sufficient to produce the spectral signature of ADH on Charon. Finally, as a consistency check, we note that the rates at which fluid is transferred to the surface are not sufficient to violate the implicit assumption made here, that they do not transport significant heat to the surface. The fluids move from a region of about 176 K to the surface, at about 60 K, a drop of ∼102 K; using the heat capacities of water and ammonia, we estimate the excess heat stored to be ∼5 × 105 J kg−1 , plus the latent heats of fusion. No more than 106 J kg−1 can be stored in the liquid. At the globally averaged rate at which it is emplaced, the maximum heat flux associated with cryovolcanism is <10−3 mW m−2 , which is small compared to the large-scale heat flux. Likewise, the loss of liquid is necessarily limited to about 7% of the water fraction of the liquid, The largest effect this loss of liquid may have is in accelerating the moment at which the liquid reaches the eutectic and cryovolcanism stops. We did not account for this effect in our code, but we estimate it would accelerate that moment by about 0.07 × 3.5 = 0.2 Gyr, an amount within the range of other uncertainties. 7. Final speculations Our thermal evolution models have great significance toward understanding the evolution of KBOs generally, but Charon in particular. It may seem at first that the models cannot apply in the context of the hypothesized formation of the Pluto–Charon system via giant impact on Pluto (Canup, 2005). However, these models tend to favor scenarios in which Charon strikes Pluto and escapes more or less intact, but gravitationally bound. During these collisions, Charon is predicted to experience only minor (∼30 K) temperature increases; later tidal dissipation would only heat Charon an additional 20 K. This pulse of heat would then dissipate on the

709

thermal diffusion timescale R 2p ρ¯ c p /k which, for T = 100 K, is approximately 1.5 Gyr. There is hope, then, that after 4.5 Gyr the temperature difference between the real Charon and our models would approach 50 exp(−4.5/1.5) K, or only a few degrees. Our models are therefore probably not completely inapplicable. The largest discrepancy probably arises from the assumption that Charon accreted cold, thereby underestimating the effects of heating early in Charon’s evolution, and the extent of differentiation, A Charon that collided with Pluto would have a rocky core larger than we predict. This would increase the temperatures today, due to increased stored heat. Charon in particular, then, is somewhat favored over an identical but isolated KBO to experience cryovolcanism, but could still be cryovolcanic if it was isolated and heated only by long-lived radionuclides. Our models point out an interesting, relevant fact: Charon is likely to have partially differentiated at the time of the Pluto impact. The models of Canup (2005) favor the impact of an undifferentiated Charon in order to more likely produce a planet–Moon system. If the impact did indeed take place before Charon differentiated, our models imply that the Pluto–Charon binary was formed in the first 70 Myr of the Solar System. As for observational evidence of cryovolcanism today on Charon, we find it intriguing to note that the albedo maps of Charon by Buie et al. (1992) include not just apparent polar caps of highalbedo (>0.7) material, presumably water ice; they also reveal a region of high-albedo material centered on the sub-Pluto point. We estimate this region occupies about 10% of Charon’s surface area, in contrast to the two polar caps, which we estimate comprise about 15% of Charon’s total surface area. If cryovolcanism is restricted to one location on Charon, we speculate that it is this equatorial region. An additional complication to the above analysis is the possible presence of methanol in the ice. Methanol is also an effective antifreeze, and a eutectic mixture of water and methanol freezes at about the same temperature as a water–ammonia mixture, 180 K (Kargel, 1994). A ternary composition of water, ammonia and methanol (in the ratio 47:23:30) freezes at an even lower temperature, 153 K (Kargel, 1994). Methanol is routinely detected in comets, at levels relative to water that range from <0.15 to 6%; in most comets the ratio is ≈2% (Bockelee-Morovan et al., 2004). Methanol is in all likelihood at least as important as ammonia in reducing the melting point of ice in KBOs, and their effects will probably reinforce each other. This would manifest itself in KBOs in two ways. First, the lower melting point would probably trigger differentiation at a larger radius R diff , which would lead to the formation of a larger rocky core that could retain heat for longer. Second, the last liquid, being driven towards the ternary eutectic (just as ammonia was concentrated in the liquid in our models), the last freezing would not occur until temperatures of 153 K were reached. In our canonical model, the liquid layer cooled at a rate 12 K Gyr−1 until right before it froze (i.e., at t < 4 Gyr), then cooled at a rate 38 K Gyr−1 during the 0.5 Gyr before it froze. At these rates, it would take and additional 0.6–2 Gyr to cool from 176 K to 153 K. Inclusion of methanol is likely to reduce the minimum radius of a KBO with liquid today, probably to somewhere in the 400–500 km range. Clearly, future studies should include the effects of methanol as well as ammonia. Unfortunately, determination of whether other specific KBOs beside Charon could be experiencing cryovolcanism must await more precise determinations of their radii and densities. Based on tentative observations of ammonia hydrates on their surfaces, Quaoar (Jewitt and Luu, 2004) and Orcus (Barucci et al., 2008) are obvious candidates for cryovolcanic KBOs. Both KBOs are now known to possess satellites (Brown and Suer, 2007), so a determination of their masses should not be long coming. Unfortunately, radius determinations are notoriously difficult. For example,

710

S.J. Desch et al. / Icarus 202 (2009) 694–714

Stansberry et al. (2008) estimate Quaoar’s radius to be 422 ± 99 km based on its thermal emission, whereas Brown and Trujillo (2004), using the Hubble Space Telescope to image it, estimate its radius to be 615 ± 95 km. In the first case, Quaoar would have frozen solid 2 Gyr ago, but in the latter case, it would still be experiencing cryovolcanism. Likewise, Stansberry et al. (2008) estimate the radius of Orcus to be 473 ± 36 km, based on its thermal emission, but if thermal emission modeling underestimates the radius to the same degree implied by the data for Quaoar, then Orcus, too, may retain subsurface liquid and be experiencing cryovolcanism. Alternatively, these bodies are large enough that they may possibly experience methanol–ammonia–water cryovolcanism. We close with the provocative thought that KBOs may be astrobiologically significant. The proximity of subsurface oceans to a hot rocky core implies the possibility of hydrothermal circulation. The highly ammonia-rich nature of the liquid ( X > 0.2, typically), coupled with the extreme cold temperatures (≈200 K) certainly make these subsurface oceans a challenge for life, but not impossible in any obvious way. Certainly, there is abundant liquid in the outer Solar System. Assuming X = 2%, each KBO the size of Charon or larger would contain ≈2 × 1019 kg of liquid at the present time; KBOs the size of Pluto would contain ≈2 × 1020 kg of liquid. There are three KBOs (Triton, Pluto, Eris) with radii larger than 1200 km, and together these might contain about 6 × 1020 kg of liquid today. Adopting a conservative power law size distribution dN /dD ∝ D −3 (e.g., Donnison, 2006), we would estimate about 3 times as many KBOs in the 600–1200 km size range, or about 9. The total amount of liquid on objects in this size range would be about 2 × 1020 kg. Altogether, we estimate a total mass of liquid in the Kuiper Belt, today, ∼1021 kg, which is comparable to the mass of Earth’s oceans, 1.6 × 1021 kg (Mottl et al., 2007). It is probable that there exists as much liquid in the outer Solar System, now, as exists on Earth. Acknowledgments We thank Southwest Research Institute for hosting SJD during the final writing of this paper. We thank Wendy Hawley for assistance with running the calculations. We thank Bill McKinnon and Michael Brown for stimulating discussions. Appendix A. Radiogenic heating rates As described in the text, the radioactive heat released per kilogram of rock is only approximately known, because of uncertainties in the energy carried off by the neutrinos. Few descriptions of this calculation exist in the literature. For completeness, we provide details here. Consider the case of 40 K decay. In 89.14% of decays, this isotope emits an electron and a neutrino, decaying to 40 Ca. Only the kinetic energy of the electron, which is 0.5083 MeV, can heat the material. In 10.66% of decays, 40 K captures an electron and emits a neutrino, decaying to an excited state of 40 Ar. The 40 Ar then decays to the ground state with the emission of a 1.4609 MeV gamma ray, which can heat the material. The 0.20% of remaining decays are by electron capture directly to the ground state of 40 Ar, creating a neutrino but not heating the material. The average heat energy generated per decay of 40 K is therefore  E = 0.6088 MeV, which is considerably smaller than the average mass deficit of the decays, 1.332 MeV. Put another way, the neutrinos carry off an average of 0.72 MeV per 40 K decay, which is consistent with their predicted spectrum of energies (Araki et al., 2005). We note that the predicted heat energy is most sensitive to the uncertainty in the kinetic energy of the emitted β particles, and is estimated at about 10%.

The decay chains of U and Th also involve beta decays and there will be energy losses by neutrinos as well, although these are not as severe. The decay chains of 232 Th, 238 U and 235 U involve the emission of 4, 6 and 4 neutrinos, as they decay to 208 Pb, 206 Pb and 207 Pb, respectively. The energies of these neutrinos do not ever exceed 3.3 MeV, and their average energies are comparable to those emitted during beta decay of 40 K, about 1 MeV. For the sake of simplicity, we calculate the heat energy released in each decay chain as the mass deficit between the radionuclide and their stable daughters, minus 1 MeV per emitted neutrino. Our assumed values of the heat energies per decay of each isotope, and their half-lives, are given in Table 1. To converting the energy per decay to a heating rate per kilogram of rock, we use the initial abundances of each isotope as computed by Lodders (2003), extrapolated back to the origin of the Solar System. Their abundances (relative to 106 Si atoms) are provided in Table 1. As these elements should condense easily in the Kuiper Belt they are assumed to be present in KBO rock in the same proportions. To convert these abundances into weights per kilogram of rock, we assume that the rock is formed by complete condensation of the (Lodders, 2003) solar abundances of Ca, Al, Fe, Mg, Si and S (plus O), into gehlenite, spinel, olivine, pyroxene, and troilite. By appropriately averaging the densities of these components we calculate an average density of this mixture is 3590 kg m−3 , corresponding very closely to the average density of low-porosity chondrites (Wilkison and Robinson, 2000). Each gram of this combination of rock contains 1/151 moles of Si atoms. We calculate initial weight percents in rock of each radionuclide, listed in Table 1. The masses of each isotope decrease exponentially with time according to their half lives. Appendix B. Equation of state of water–ammonia system Our thermal evolution model requires knowledge of e ( T ), the energy per mass needed to raise the temperature of an ammonia– water mixture from 0 K to a temperature T, given the mass fraction that is ammonia, X . This appendix describes how we compute this quantity. We first require knowledge of which part of the phase diagram we are in, which is a function of X and T . Our approximate phase diagram is broken up into 5 regimes, as depicted in Fig. 3. Regime 1 refers to temperatures below the solidus, which we take to be 176 K, the melting point of ADH. Actually, we model the melting of ADH over a small range in temperatures from 174 to 178 K, so regime 1 applies to any X and T < 174 K. In this regime, the fractions of the total mass of water plus ammonia in each phase are given by g H2 O(s) = 1 − X / X c , g ADH(s) = X / X c , g H2 O(l) = 0, g NH3 (l) = 0.

(B.1)

Regime 2 refers to the small range of temperatures over which ADH is modeled to melt, 174 to 178 K. In this regime, the mass fractions are g H2 O(s) = 1 − X / X c , g ADH(s) = ( X / X c )(178 K − T )/(4 K), g H2 O(l) = X [1/ X c − 1]( T − 174 K)/(4 K), g NH3 (l) = X ( T − 174 K)/(4 K).

(B.2)

Cryovolcanism on Kuiper belt objects

Regime 3 refers to temperatures between the solidus (actually, above 178 K) and the liquidus. To simplify matters, we use a crude approximation to the liquidus temperature:



T liq ( X ) =

273 − 95( X / X c )2 K,

Xb  X  Xc ,

271 K,

0  X < Xb ,

X liq ( T ) = X c

273 K − T

T e 2 ( T ) = e 1 (174 K) +

,

178 K  T  271 K.

(B.4)

This quantity is used via the “lever rule” to find the mass fractions of the various phases in regime 3:

g NH3 (l) ( T  )c NH3 (l) dT 

174 K

T 174 K

g ADH(s) ( T  )c ADH(s) ( T  ) dT 

+

g ADH(s) = 0,

c H2 O(s) ( T  ) dT 

+ g H2 O(s)

T

g H2 O(s) = 1 − X / X liq ( T ),

174 K

g H2 O(l) = X / X liq ( T ) − X ,

T 

g NH3 (l) = X .

(B.5)

Regime 4 applies only if X < X b , so that some solid water ice remains unmelted even at 271 K. The melting of the remaining water ice is assumed to take place over the temperature range 271 to 275 K, and the mass fractions in this regime are given by

+

−

 ∂ g ADH(s) lADH dT  . ∂T

174 K



+

g ADH(s) = 0,

X

g NH3 (l) = X .

(B.6)

Finally, regime 5 refers to temperatures above the liquidus; that is, T > T liq ( X ) if X > X b , or T > 275 K for X  X b . In this regime, the mass fractions are simply g H2 O(s) = 0,

+



4K

T − 174 K 2K





lADH +

 X c c NH3 (l) +



182 K − T

c ADH(s) (174 K)

2K



T − 174 K 2K

 (1 − X c )c H2 O(l) .

g H2 O(l) ( T  )c H2 O(l) dT 

178 K

g NH3 (l) = X .

(B.7)

The next problem is determining the energy required to raise a mass of water and ammonia (with given X ) to a temperature T . We will find it convenient to define c H2 O(s) ( T  ) dT  = 3.87 × 104



T

2

100 K

J kg−1

(B.8)

0

T

g NH3 (l) c NH3 (l) dT 

+ 178 K

T

g H2 O(s) ( T  )c H2 O(s) ( T  ) dT 

+ 178 K

T  +

and c H2 O(s) ( T  ) dT  = 56.0 × 104



T 100 K

2

J kg−1 . (B.9)





T

c H2 O(s) ( T ) dT + g ADH(s) 0

c ADH(s) ( T  ) dT  , (B.10)

X  e ADH(s) ( T ) − e H2 O(s) ( T ) . Xc

(B.14)

After much tedious algebra we derive





 X + X c NH3 (l) − c H2 O(l) ( T − 178 K) + (1 − r ) Xc   lH2 O × + 2(95 K) c H2 O(l) − c H2 O(s) r

+

0

or e 1 ( T ) = e H2 O(s) ( T ) +

 ∂ g H2 O(s) lH2 O dT  . ∂T

e 3 ( T ) = e 2 (178 K) + e H2 O(s) ( T ) − e H2 O(s) (178 K)

Within regime 1 (T < 174 K), the specific energy e (per mass of water and ammonia) required to reach a temperature T is simple:

T

−

178 K

0

e 1 ( T ) = g H2 O(s)

e H2 O(s) ( T ) − e H2 O(s) (174 K)

Within regime 3 (178 K < T < T liq ( X )), one must again account for melting, this time of water ice, including its latent heat of melting, lH2 O = 3.335 × 105 J kg−1 , so

g H2 O(l) = 1 − X ,

e ADH(s) ( T ) =





(B.13)

e 3 ( T ) = e 2 (178 K) +

T

Xc





T

g ADH(s) = 0,

e H2 O(s) ( T ) =

X

T − 174 K

Xc



g H2 O(l) = [1 − X / X b ]( T − 271 K)/(4 K) + X / X b − X ,

(B.12)

Only because we have simplified the form of X liq ( T ) does this permit an analytic solution: e 2 ( T ) = e 1 (174 K) + 1 −

g H2 O(s) = [1 − X / X b ](275 K − T )/(4 K),

T

g H2 O(l) ( T  )c H2 O(l) dT 

174 K

T +

1/ 2

95 K

As the temperature is increased into regime 2 (174 K < T < 178 K), one must account for the fact that the abundances of products change as the ADH melts, and one must include the latent heat of melting of ADH, lADH = 1.319 × 105 J kg−1 . Thus,

(B.3)

where X b = 0.0466. The advantage of this approximation is that it can be algebraically inverted to find X liq ( T ), the ammonia concentration of the liquid at a temperature T :



711

(B.11)

2(95 K)2 3(273 K)



c H2 O(s) 1 + r + r 2

 

,

(B.15)

where r = [(273 K − T )/(95 K)]1/2 , and c H2 O(s) is evaluated at T = 273 K.

712

S.J. Desch et al. / Icarus 202 (2009) 694–714

If X > X b , then further heating takes the system into regime 5 (T > T liq ( X )), so

T





c H2 O(l) dT 

e 5 ( T ) = e 4 T liq ( X ) + g H2 O(l) T liq ( X )

T

c NH3 (l) dT  ,

+ g NH3 (l)

(B.16)

T liq ( X )

or











e 5 ( T ) = e 3 T liq ( X ) + Xc NH3 (l) + (1 − X )c H2 O(l) T − T liq ( X )

(B.17) Otherwise, if X < X b , the system must first pass through regime 4 (271 K < T < 275 K), in which the last of the solid water ice melts, so

T e 4 ( T ) = e 3 (271 K) +

g H2 O(l) ( T  )c H2 O(l) dT 

271 K

T + g NH3 (l)

c NH3 (l) dT 

271 K

T +

g H2 O(s) ( T  )c H2 O(s) ( T  ) dT 

271 K

T  −

+

 ∂ g H2 O(s) lH2 O dT  , ∂T

(B.18)

271 K

or





e 4 ( T ) = e 3 (271 K) + Xc NH3 (l) +

 + 1− 



X Xb

X Xb

T − 271 K

  − X c H2 O(l) ( T − 271 K)



4K 1

× lH2 O + c H2 O(l) ( T − 271 K) 2

 1 + c H2 O(s) (279 K − T ) ,

(B.19)

2

where c H2 O(s) is evaluated at 273 K. Finally, in regime 5 (T > 275 K), if X  X b ,

T e 5 ( T ) = e 4 (275 K) + g H2 O(l)

c H2 O(l) dT 

275 K

T + g NH3 (l)

c NH3 (l) dT  ,

(B.20)

275 K

or





e 5 ( T ) = e 4 (275 K) + Xc NH3 (l) + (1 − X )c H2 O(l) ( T − 275 K). (B.21) For a given temperature and X , these formulas allow us to determine what fraction of the water–ammonia mixture is in each phase, and the total specific energy required to heat that mixture to a temperature T .

References Agee, C.B., Li, J., Shannon, M.C., Circone, S., 1995. Pressure–temperature phase diagram for the Allende meteorite. J. Geophys. Res. 100, 17725–17740. Agnor, C.B., Hamilton, D.P., 2006. Neptune’s capture of its moon Triton in a binaryplanet gravitational encounter. Nature 441, 192–194. Allen, M., Delitsky, M., Huntress, W., Yung, Y., Ip, W.-H., 1987. Evidence for methane and ammonia in the coma of Comet P/Halley. Astron. Astrophys. 187, 502–512. Arakawa, M., Maeno, N., 1994. Effective viscosity of partially melted ice in the ammonia–water system. Geophys. Res. Lett. 21, 1515–1518. Araki, T., and 92 colleagues, 2005. Measurement of neutrino oscillation with KamLAND: Evidence of spectral distortion. Phys. Rev. Lett. 94. 081801. Barucci, M.A., Merlin, F., Guilbert, A., de Bergh, C., Alvarez-Candal, A., Hainaut, O., Doressoundiram, A., Dumas, C., Owen, T., Coradini, A., 2008. Surface composition and temperature of the TNO Orcus. Astron. Astrophys. 479, L13–L16. Bauer, J.M., Roush, T.L., Geballe, T.R., Meech, K.J., Owen, T.C., Vacca, W.D., Rayner, J.T., Jim, K.T.C., 2002. The Near infrared spectrum of Miranda: Evidence of crystalline water ice. Icarus 158, 178–190. Beaver, J.E., Wagner, R.M., Schleicher, D.G., Lutz, B.L., 1990. Anomalous molecular abundances and the depletion of NH2 in Comet P/Giacobini–Zinner. Astrophys. J. 360, 696–701. Bird, M.K., Huchtmeier, W.K., Gensheimer, P., Wilson, T.L., Janardhan, P., Lemme, C., 1997. Radio detection of ammonia in Comet Hale-Bopp. Astron. Astrophys. 325, L5–L8. Bockelée-Morvan, D., Crovisier, J., Mumma, M.J., Weaver, H.A., 2004. The composition of cometary volatiles. In: Festou, M.C., Keller, H.U., Weaver, H.A. (Eds.), Comets II. Univ. of Arizona, Tucson, pp. 391–423. Brown, M.E., 2000. Near-Infrared spectroscopy of Centaurs and irregular satellites. Astron. J. 119, 977–983. Brown, M.E., Calvin, W.M., 2000. Evidence for crystalline water and ammonia ices on Pluto’s satellite Charon. Science 287, 107–109. Brown, M.E., Suer, T.-A., 2007. Satellites of 2003 AZ_84, (50000), (55637), and (90482). IAU Circ. 8812, 1. Brown, M.E., Trujillo, C.A., 2004. Direct measurement of the size of the large Kuiper belt object (50000) Quaoar. Astron. J. 127, 2413–2417. Brown, R.H., Cruikshank, D.P., Tokunaga, A.T., Smith, R.G., Clark, R.N., 1988. Search for volatiles on icy satellites. I. Europa. Icarus 74, 262–271. Buie, M.W., Tholen, D.J., Horne, K., 1992. Albedo maps of Pluto and Charon — Initial mutual event results. Icarus 97, 211–227. Buie, M.W., Grundy, W.M., Young, E.F., Young, L.A., Stern, S.A., 2006. Orbits and photometry of Pluto’s satellites: Charon, S/2005 P1, and S/2005 P2. Astron. J. 132, 290–298. Cahill, D.G., Watson, S.K., Pohl, R.O., 1992. Lower limit to the thermal conductivity of disordered crystals. Phys. Rev. B 46, 6131–6140. Canup, R.M., 2005. A giant impact origin of Pluto–Charon. Science 307, 546–550. Castillo-Rogez, J., 2006. Internal structure of Rhea. J. Geophys. Res. 111. E11005. Castillo-Rogez, J.C., Matson, D.L., Sotin, C., Johnson, T.V., Lunine, J.I., Thomas, P.C., 2007a. Iapetus’ geophysics: Rotation rate, shape, and equatorial ridge. Icarus 190, 179–202. Castillo-Rogez, J.C., McCord, T.B., Davies, A.G., 2007b. Ceres: Evolution and present state. Lunar Planet. Sci. 38. Abstract 2006. Charnley, S.B., Rodgers, S.D., 2002. The end of interstellar chemistry as the origin of nitrogen in comets and meteorites. Astrophys. J. 569, L133–L137. Clauser, C., Huenges, E., 1995. Rock physics and phase relations. In: Ahrens, T.J. (Ed.), A Handbook of Physical Constants. In: AGU Reference Shelf, vol. 3. American Geophysical Union, Washington, pp. 105–126. Consolmagno, G.J., 1985. Resurfacing Saturn’s satellites — Models of partial differentiation and expansion. Icarus 64, 401–413. Consolmagno, G.J., Britt, D.T., 1998. The density and porosity of meteorites from the Vatican collection. Meteorit. Planet. Sci. 33, 1231–1241. Cook, J.C., Desch, S.J., Roush, T.L., Trujillo, C.A., Geballe, T.R., 2007. Near-infrared spectroscopy of Charon: Possible evidence for cryovolcanism on Kuiper belt objects. Astrophys. J. 663, 1406–1419. Cooper, J.F., Christian, E.R., Richardson, J.D., Wang, C., 2003. Proton irradiation of Centaur, Kuiper Belt, and Oort cloud objects at plasma to cosmic ray energy. Earth Moon Planets 92, 261–277. Crawford, G.D., Stevenson, D.J., 1988. Gas-driven water volcanism in the resurfacing of Europa. Icarus 73, 66–79. Croft, S.K., 1990. Physical cryovolcanism on Triton. Lunar Planet. Sci. 21. Abstract 246. Croft, S.K., Lunine, J.I., Kargel, J., 1988. Equation of state of ammonia–water liquid — Derivation and planetological applications. Icarus 73, 279–293. De Sanctis, M.C., Capria, M.T., Coradini, A., 2001. Thermal evolution and differentiation of Edgeworth–Kuiper belt objects. Astron. J. 121, 2792–2799. Dombard, A.J., McKinnon, W.B., 2001. Formation of grooved terrain on Ganymede: Extensional instability mediated by cold, superplastic creep. Icarus 154, 321– 336. Donnison, J.R., 2006. The size distribution of trans-neptunian objects. Planet. Space. Sci. 54, 243–250.

Cryovolcanism on Kuiper belt objects

Dumas, C., Terrile, R.J., Brown, R.H., Schneider, G., Smith, B.A., 2001. Hubble space telescope NICMOS spectroscopy of Charon’s leading and trailing hemispheres. Astron. J. 121, 1163–1170. Durham, W.B., Stern, L.A., 2001. Rheological properties of water ice-applications to satellites of the outer planets. Annu. Rev. Earth Planet. Sci. 29, 295–330. Durham, W.B., Kirby, S.H., Stern, L.A., 1993. Flow of ices in the ammonia–water system. J. Geophys. Res. 98, 17667–17682. Ellsworth, K., Schubert, G., 1983. Saturn’s icy satellites — Thermal and structural models. Icarus 54, 490–510. Emery, J.P., Burr, D.M., Cruikshank, D.P., Brown, R.H., Dalton, J.B., 2005. Near-infrared (0.8–4.0 μm) spectroscopy of Mimas, Enceladus, Tethys, and Rhea. Astron. Astrophys. 435, 353–362. Fagents, S.A., 2003. Considerations for effusive cryovolcanism on Europa: The postGalileo perspective. J. Geophys. Res. (Planets) 108, 5139. Feldman, P.D., Fournier, K.B., Grinin, V.P., Zvereva, A.M., 1993. The abundance of ammonia in Comet P/Halley derived from ultraviolet spectrophotometry of NH by ASTRON and IUE. Astrophys. J. 404, 348–355. Figueredo, P.H., Greeley, R., 2004. Resurfacing history of Europa from pole-to-pole geological mapping. Icarus 167, 287–312. Fortes, A.D., 2007. Metasomatic clathrate xenoliths as a possible source for the south polar plumes of Enceladus. Icarus 191, 743–748. Fortes, A.D., Wood, I.G., Brodholt, J.P., Voˇcadlo, L., 2003. The structure, ordering and equation of state of ammonia dihydrate (NH3 ·2H2 O). Icarus 162, 59–73. Fortes, A.D., Grindrod, P.M., Trickett, S.K., Voˇcadlo, L., 2007. Ammonium sulfate on Titan: Possible origin and role in cryovolcanism. Icarus 188, 139–153. Giese, B., Wagner, R., Hussmann, H., Neukum, G., Perry, J., Helfenstein, P., Thomas, P.C., 2008. Enceladus: An estimate of heat flux and lithospheric thickness from flexurally supported topography. Geophys. Res. Lett. 35. L24204. Grasset, O., Sotin, C., Deschamps, F., 2000. On the internal structure and dynamics of Titan. Planet. Space Sci. 48, 617–636. Grundy, W.M., Buie, M.W., Stansberry, J.A., Spencer, J.R., Schmitt, B., 1999. Nearinfrared spectra of icy outer Solar System surfaces: Remote determination of H2 O ice temperatures. Icarus 142, 536–549. Gulbis, A.A.S., and 12 colleagues, 2006. Charon’s radius and atmospheric constraints from observations of a stellar occultation. Nature 439, 48–51. Hansen, C.J., Esposito, L., Stewart, A.I.F., Colwell, J., Hendrix, A., Pryor, W., Shemansky, D., West, R., 2006. Enceladus’ water vapor plume. Science 311, 1422–1425. 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. Horányi, M., and 30 colleagues, 2008. The Student dust counter on the New Horizons Mission. Space Sci. Rev. 140, 387–402. Humes, D.H., 1980. Results of Pioneer 10 and 11 meteoroid experiments — Interplanetary and near-Saturn. J. Geophys. Res. 85, 5841–5852. Hussmann, H., Sohl, F., Spohn, T., 2006. Subsurface oceans and deep interiors of medium-sized outer planet satellites and large trans-neptunian objects. Icarus 185, 258–273. Jewitt, D.C., Luu, J., 2004. Crystalline water ice on the Kuiper belt object (50000) Quaoar. Nature 432, 731–733. Kargel, J.S., 1992. Ammonia-water volcanism on icy satellites — Phase relations at 1 atmosphere. Icarus 100, 556–574. Kargel, J.S., 1994. Cryovolcanism on the icy satellites. Earth Moon Planets 67, 101– 113. Kargel, J.S., Pozio, S., 1996. The volcanic and tectonic history of Enceladus. Icarus 119, 385–404. Kargel, J.S., Strom, R.G., 1990. Cryovolcanism on Triton. Lunar Planet. Sci. 21. Abstract 599. Kargel, J.S., Croft, S.K., Lunine, J.I., Lewis, J.S., 1991. Rheological properties of ammonia–water liquids and crystal-liquid slurries — Planetological applications. Icarus 89, 93–112. Kawakita, H., Watanabe, J.-i., 2002. Revised fluorescence efficiencies of cometary NH2 : Ammonia abundance in comets. Astrophys. J. 572, L177–L180. Kieffer, S.W., Lu, X., Bethke, C.M., Spencer, J.R., Marshak, S., Navrotsky, A., 2006. A clathrate reservoir hypothesis for Enceladus’ south polar plume. Science 314, 1764. Kivelson, M.G., Khurana, K.K., Russell, C.T., Volwerk, M., Walker, R.J., Zimmer, C., 2000. Galileo magnetometer measurements: A stronger case for a subsurface ocean at Europa. Science 289, 1340–1343. Klinger, J., 1980. Influence of a phase transition of ice on the heat and mass balance of comets. Science 209, 271. ´ Leliwa-Kopystynski, J., Maruyama, M., Nakajima, T., 2002. The water–ammonia phase diagram up to 300 MPa: Application to icy satellites. Icarus 159, 518–528. Lodders, K., 2003. Solar System abundances and condensation temperatures of the elements. Astrophys. J. 591, 1220–1247. Lopes, R.M.C., and 43 colleagues, 2007. Cryovolcanic features on Titan’s surface as revealed by the Cassini Titan Radar Mapper. Icarus 186, 395–412. Lorenz, R.D., Shandera, S.E., 2001. Physical properties of ammonia-rich ice: Application to Titan. Geophys. Res. Lett. 28, 215–218. Maret, S., Bergin, E.A., Lada, C.J., 2006. A low fraction of nitrogen in molecular form in a dark cloud. Nature 442, 425–427.

713

Mastrapa, R.M.E., Brown, R.H., 2006. Ion irradiation of crystalline H2 O ice: Effect on the 1.65-μm band. Icarus 183, 207–214. McKinnon, W.B., 1984. On the origin of Triton and Pluto. Nature 311, 355–358. McKinnon, W., 1998. Geodynamics of icy satellites. In: Solar System Ices. In: ASSL Series, vol. 227. Kluwer, Dordrecht, p. 525. McKinnon, W.B., 1999. Convective instability in Europa’s floating ice shell. Geophys. Res. Lett. 26, 951–954. McKinnon, W.B., 2002. On the initial thermal evolution of Kuiper Belt objects. In: Asteroids, Comets, and Meteors: ACM 2002. ESA SP-500, pp. 29–38. Meier, R., Eberhardt, P., Krankowsky, D., Hodges, R.R., 1994. Ammonia in Comet P/Halley. Astron. Astrophys. 287, 268–278. Meier, R., Wellnitz, D., Kim, S.J., A’Hearn, M.F., 1998. The NH and CH bands of Comet C/1996 B2 (Hyakutake). Icarus 136, 268–279. Mottl, M.J., Glazer, B.T., Kaiser, R.I., Meech, K.J., 2007. Water and astrobiology. Chem. Erde 67, 253–282. Multhaup, K., Spohn, T., 2007. Stagnant lid convection in the mid-sized icy satellites of Saturn. Icarus 186, 420–435. Navrotsky, A., 1995. Thermodynamic properties of minerals. In: Ahrens, T.J. (Ed.), Mineral Physics and Crystallography, A Handbook of Physical Constants. Amer. Geophys. Union, pp. 18–28. Nimmo, F., Pappalardo, R.T., 2001. Furrow flexure and ancient heat flux on Ganymede. Geophys. Res. Lett. 31. L19701. Nimmo, F., Pappalardo, R.T., Giese, B., 2002. Effective elastic thickness and heat flux estimates on Ganymede. Geophys. Res. Lett. 29. 070000-1. Olkin, C.B., Wasserman, L.H., Franz, O.G., 2003. The mass ratio of Charon to Pluto from Hubble Space Telescope astrometry with the fine guidance sensors. Icarus 164, 254–259. Palmer, P., Wootten, A., Butler, B., Bockelee-Morvan, D., Crovisier, J., Despois, D., Yeomans, D.K., 1996. Comet Hyakutake: First secure detection of ammonia in a comet. Bull. Am. Astron. Soc. 28, 927. Peale, S.J., 1988. Speculative histories of the uranian satellite system. Icarus 74, 153– 171. Person, M.J., Elliot, J.L., Gulbis, A.A.S., Pasachoff, J.M., Babcock, B.A., Souza, S.P., Gangestad, J., 2006. Charon’s radius and density from the combined data sets of the 2005 July 11 occultation. Astron. J. 132, 1575–1580. Plescia, J.B., 1987. Geological terrains and crater frequencies on Ariel. Nature 327, 201–204. Plescia, J.B., 1988. Cratering history of Miranda — Implications for geologic processes. Icarus 73, 442–461. Porco, C.C., and 24 colleagues, 2006. Cassini observes the active south pole of Enceladus. Science 311, 1393–1401. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in FORTRAN: The Art of Scientific Computing, second ed. Cambridge University Press, Cambridge. Prialnik, D., Bar-Nun, A., 1990. Heating and melting of small icy satellites by the decay of Al-26. Astrophys. J. 355, 281–286. Röttger, K., Endriss, A., Ihringer, J., Doyle, S., Kuhs, W.F., 1994. Lattice constants and thermal expansion of H2 O and D2 O ice Ih between 10 and 265 K. Acta Crystallogr. B 50, 644–648. Ruiz, J., 2003. Heat flow and depth to a possible internal ocean on Triton. Icarus 166, 436–439. Ruiz, J., Tejero, R., 2000. Heat flows through the ice lithosphere of Europa. J. Geophys. Res. 105, 29283–29290. Schaller, E.L., Brown, M.E., 2007. Detection of methane on Kuiper belt object (50000) Quaoar. Astrophys. J. 670, L49–L51. Schenk, P.M., 1991. Fluid volcanism on Miranda and Ariel — Flow morphology and composition. J. Geophys. Res. 96, 1887–1906. Schenk, P.M., Zahnle, K., 2007. On the negligible surface age of Triton. Icarus 192, 135–149. Schubert, G., Anderson, J.D., Travis, B.J., Palguta, J., 2007. Enceladus: Present internal structure and differentiation by early and long-term radiogenic heating. Icarus 188, 345–355. Showman, A.P., Mosqueira, I., Head, J.W., 2004. On the resurfacing of Ganymede by liquid-water volcanism. Icarus 172, 625–640. Shchuko, O.B., Coradini, A., Orosei, R., Shchuko, S.D., 2006. Varuna: Thermal evolution. Adv. Space Res. 38, 1946–1951. Shulman, L.M., 2004. The heat capacity of water ice in interstellar or interplanetary conditions. Astron. Astrophys. 416, 187–190. Sicardy, B., and 44 colleagues, 2006. Charon’s size and an upper limit on its atmosphere from a stellar occultation. Nature 439, 52–54. Sirono, S.-I., Yamamoto, T., 1997. Thermal conductivity of granular materials relevant to the thermal evolution of cometary nuclei. Planet. Space Sci. 45, 827–834. Smith, B.A., and 10 colleagues, 1979. The Galilean satellites and Jupiter — Voyager 2 imaging science results. Science 206, 927–950. Solomatov, V.S., 1995. Scaling of temperature- and stress-dependent viscosity convection. Phys. Fluids 7, 266–274. Song, M., Cole, D.M., Baker, I., 2006. Investigation of Newtonian creep in polycrystalline ice. Philos. Mag. Lett. 86, 763–771.

714

S.J. Desch et al. / Icarus 202 (2009) 694–714

Spencer, J.R., 1989. The thermal inertias of Rhea, Dione, and Tethys from groundbased radiometry. Bull. Am. Astron. Soc. 21, 983. Spencer, J.R., Moore, J.M., 1992. The influence of thermal inertia on temperatures and frost stability on Triton. Icarus 99, 261–272. Spencer, J.R., Tamppari, L.K., Martin, T.Z., Travis, L.D., 1999. Temperatures on Europa from Galileo PPR: Nighttime thermal anomalies. Science 284, 1514–1516. Spohn, T., Schubert, G., 2003. Oceans in the icy Galilean satellites of Jupiter? Icarus 161, 456–467. Stansberry, J., Grundy, W., Brown, M., Cruikshank, D., Spencer, J., Trilling, D., Margot, J.-L., 2008. Physical properties of Kuiper Belt and Centaur objects: Constraints from the Spitzer space telescope. In: The Solar System Beyond Neptune, pp. 161–179. Stern, S.A., 2008. The New Horizons Pluto Kuiper belt mission: An overview with historical context. Space Sci. Rev. 140, 3–21. Stevenson, D.J., 1982. Volcanism and igneous processes in small icy satellites. Nature 298, 142–144.

Strazzulla, G., Palumbo, M.E., 1998. Evolution of icy surfaces: An experimental approach. Planet. Space Sci. 46, 1339–1348. Thomas, P.J., Reynolds, R.T., Squyres, S.W., Cassen, P.M., 1987. The viscosity of Miranda. Lunar Planet. Sci. XVIII, 1016–1017 (abstract). Tobie, G., Grasset, O., Lunine, J.I., Mocquet, A., Sotin, C., 2005. Titan’s internal structure inferred from a coupled thermal-orbital model. Icarus 175, 496–502. Verbiscer, A.J., Peterson, D.E., Skrutskie, M.F., Cushing, M., Helfenstein, P., Nelson, M.J., Smith, J.D., Wilson, J.C., 2006. Near-infrared spectra of the leading and trailing hemispheres of Enceladus. Icarus 182, 211–223. Wilkison, S.L., Robinson, M.S., 2000. Bulk density of ordinary chondrite meteorites and implications for asteroidal internal structure. Meteorit. Planet. Sci. 35, 1203– 1213. Yomogida, K., Matsui, T., 1983. Physical properties of ordinary chondrites. J. Geophys. Res. 88, 9513–9533. Zahnle, K., Schenk, P., Levison, H., Dones, L., 2003. Cratering rates in the outer Solar System. Icarus 163, 263–289.