Thermophysical evolution of planetesimals in the Primordial Disk
Abstract
The Primordial Disk of small icy planetesimals, once located at – from the Sun, was disrupted by giant planet migration in the early Solar System. The Primordial Disk thereby became the source region of objects in the current–day Kuiper Belt, Scattered Disk, and Oort Cloud. I present the thermophysics code “Numerical Icy Minor Body evolUtion Simulator”, or NIMBUS, and use it to study the thermophysical evolution of planetesimals in the Primordial Disk prior to its disruption. Such modelling is mandatory in order to understand the behaviour of dynamically new comets from the Oort Cloud, as well as the activity of Centaurs and short–period comets from the Scattered Disk, that return pre–processed to the vicinity of the Sun. I find that bodies in the midst of the Primordial Disk with diameters ranging – lost all their CO ice on time-scales of order 0.1– depending on size, through a combination of protosolar and long–lived radionuclide heating. CO and other hypervolatiles therefore require a less volatile host for their storage. I consider two possible hosts: amorphous water ice and ice. Because of the high luminosity of the protosun, some Primordial Disk bodies may have sustained significant crystallisation, segregation, and sublimation in the uppermost few tens of meters. I discuss how this may affect coma abundance ratios and distant activity in dynamically new comets.
keywords:
methods: numerical – comets: general – Kuiper belt: general – Oort Cloud – protoplanetary discs1 Introduction
The Nice model (Tsiganis et al., 2005; Morbidelli et al., 2005; Gomes et al., 2005) postulates that the giant planets once were located in a compact orbital configuration roughly – from the Sun. Such a configuration appears to be caused by migratory behaviour of the giant planets within the gas disk from which they formed (Masset & Snellgrove, 2001; Morbidelli & Crida, 2007; Morbidelli et al., 2007; Walsh et al., 2011). Beyond lay a vast population of icy minor bodies with a sharp truncation near known as the Primordial Disk (Gomes et al., 2004). Beyond lay a substantially smaller population that still is present (Parker & Kavelaars, 2010; Parker et al., 2011) in the form of the dynamically cold Kuiper Belt. The Nice model demonstrates that a gravitational instability among the giant planets, followed by planetesimal–driven migration to their current locations, disrupted the Primordial Disk. A fraction of its members was distributed between four large populations that formed because of this event and remain to this day: 1) the Trojan swarms in 1:1 mean motion resonance with Jupiter (Morbidelli et al., 2005); 2) the dynamically hot Kuiper Belt that is superimposed on the previously mentioned dynamically cold population (Fernández & Ip, 1984; Malhotra, 1993; Levison et al., 2008); 3) the Scattered Disk with objects having perihelia near the orbit of Neptune and aphelia far beyond the Kuiper Belt (Duncan & Levison, 1997); 4) the Oort Cloud (Brasser & Morbidelli, 2013).
Numerous authors have investigated the thermophysical evolution of bodies in the Kuiper Belt (e.g. De Sanctis et al., 2001, 2007; Choi et al., 2002; Merk & Prialnik, 2003, 2006; Guilbert-Lepoutre et al., 2011; Prialnik et al., 2008; Shchuko et al., 2014). Substantial efforts have also been made to understand the thermophysical evolution of Centaurs (e.g. De Sanctis et al., 2000; Capria et al., 2000, 2009; Guilbert-Lepoutre, 2011; Sarid & Prialnik, 2009). Centaurs evolve dynamically from the Scattered Disk into the – region and a fraction eventually become Jupiter Family Comets, henceforth abbreviated JFCs (Fernández, 1980; Duncan & Levison, 1997; Levison & Duncan, 1997). The majority of current–day Centaurs therefore populate the region that once hosted the Primordial Disk. However, the solar–driven thermophysical evolution endured by Primordial Disk objects is expected to have been stronger than that of most Centaurs, and much stronger than that of current Kuiper Belt objects, for two reasons. Firstly, the luminosity of the protosun at the birthline111The birthline marks the transition from a Class I protostar to a visually detectable T Tauri star (a Class II protostar), see e. g., Larson (2003). This coincides with the end of strong accretion because the envelope and thick disk have been consumed. The thin disk remaining is referred to as the “Solar Nebula” in Solar System science terminology. The Solar Nebula ceases to exist when most gas has been removed by being heated to escape velocity (photoevaporation), consumed by the giant planets, or subjected to slow protosolar accretion. was times the current solar luminosity, and it was still times higher after , only falling below the current value into Solar System history (Palla & Stahler, 1993). Therefore, the protosun was substantially more luminous than the current Sun throughout the lifetime of the Solar Nebula (Zuckerman et al., 1995; Haisch et al., 2001; Sicilia-Aguilar et al., 2006), during which most planetesimals formed and giant planet growth was completed. Consequently, Primordial Disk objects located at could initially have evolved as if being placed at today, and the conditions at could have resembled those currently found at . Opaqueness of the Solar Nebula may have shielded planetesimals for some time period, that remains to be determined. Secondly, the estimated mean dynamical half–life of Centaurs with semi–major axes and aphelia below ranges from (Horner et al., 2004) to (Di Sisto & Brunini, 2007). The lifetime of the Primordial Disk is a matter of debate. The original version of the Nice model suggested that the gravitational instability that disrupted the disk coincided with the Late Heavy Bombardment (Gomes et al., 2005), which would imply a lifetime of at least – and perhaps as much as (Morbidelli et al., 2012; Bottke et al., 2012; Marchi et al., 2013). In their “Nice II” model scenario, Morbidelli et al. (2007) and Levison et al. (2011) strengthened the case for a late gravitational instability. However, Nesvorný & Morbidellli (2012) demonstrated that models with five giant planets and an early gravitational instability successfully reproduced a number of observable criteria, suggesting that the Primordial Disk lifetime could have been as short as . The shorter estimate of the Primordial Disk lifetime is still longer than the typical Centaur dynamical lifetime by a factor of a few (although the lifetime of some resonance–hopping Centaurs may exceed ; Bailey & Malhotra, 2009). In conclusion, Primordial Disk objects were thermally processed at a higher radiative flux and for a longer time than most modern–day Centaurs.
For these reasons it is crucial to study the thermophysical evolution of icy minor bodies in the Primordial Disk. This is particularly important because: 1) this thermal processing is not only expected to be stronger, but also pre–dates the ones taking place among Kuiper Belt objects and Centaurs, thereby providing appropriate initial conditions for thermal studies of those populations; 2) the Primordial Disk thermophysical processing would also have affected objects currently placed in the Oort Cloud, thereby challenging the concept of dynamically new comets as being largely unevolved; 3) the predicted level of thermal processing in the current Nice model paradigm, if found to be at odds with observed properties, may provide novel constraints on the conditions and durations of early Solar System key processes and events.
This paper is therefore devoted to the study of thermophysical processing of minor icy bodies in the Primordial Disk. It focuses on a number of specific questions that are inspired by previous work and recent observations, as motivated in the following. Theoretical studies of Kuiper Belt objects have shown that hypervolatiles such as , when stored within the bodies as separate and pure condensed ices, readily sublimate at such distances. De Sanctis et al. (2001) considered a diameter body orbiting at –, having a dust/ice mass ratio of 1 or 5, an porosity, containing condensed ice relative water, where the dust component contained the long–lived radionuclides , , , and with chondritic abundances. They solved the coupled heat and gas diffusion problem for a spherical–symmetric body. They modelled the body for and found that the sublimation front had withdrawn – below the surface at that time, depending on the assumed dust/ice mass ratio. Based on the asymptotically declining loss rate, De Sanctis et al. (2001) speculated that a steady–state might be reached where sublimation is balanced by recondensation and net loss of eventually ceases. If so, an inner core of ice could remain. Among several different model bodies, Choi et al. (2002) considered a case at , having a dust/ice mass ratio of 1, a porosity, containing condensed ice relative water, with long–lived radionuclides included. Their model considered heat conduction for a spherical–symmetric body, but assumed immediate escape of vapour because a proper treatment of gas diffusion was considered computationally prohibitive. They modelled the body for the Solar System lifetime and found that all was lost in little over . Their other models additionally considered the short–lived radionuclide that, if present at the time of planetesimal formation, could speed up the loss substantially. These works, that apparently are the only ones that have attempted to quantify the loss time–scale, inspire to the following related questions.
-
Q#1. Under what conditions is loss of condensed ice complete or partial in the Primordial Disk? If loss of condensed ice is complete, what are the time–scales for that loss in bodies of different size in the Primordial Disk?
Next, consider the problem of the activation distance of most comets and Centaurs, as well as extreme exceptions to the general rule. There is mounting evidence that a fundamental change in the level of activity of icy minor bodies takes place as they pass a transition region from the Sun inbound. Jewitt (2009) demonstrated that all Centaurs known to be active at the time were confined to the – region and that the lack of active Centaurs at larger distances was not an observational bias. Radio telescope searches for gaseous in a combined sample of 18 Centaurs and Kuiper Belt objects located beyond by Bockelée-Morvan et al. (2001) and Jewitt et al. (2008) yielded no detections but upper limits only. A Hubble Space Telescope survey of 53 Centaurs and Kuiper Belt objects with perihelia at did not detect dust activity in any of those targets (Li et al., 2020). Cabral et al. (2019) observed 20 Centaurs at –, of which 15 were beyond , and none had a coma. Statistics on the discovery distances of 2,096 dynamically new and long–period comets assembled by Meech et al. (2017b) shows that all distances within the orbit of Saturn are well represented. However, there is a sharp transition near beyond which merely a dozen discoveries have been made. This suggests a sudden systematic onset of activity and brightening near . Among the Centaurs that do display activity, gas production has been verified in 29P/Schwassmann–Wachmann 1 (Senay & Jewitt, 1994), 95P/Chiron (Womack & Stern, 1999), and 174P/Echeclus (Wierzchos et al., 2017). Carbon monoxide is also one of the most common and abundant gases observed in JFCs (e.g. Bockelée-Morvan et al., 2004; A’Hearn et al., 2012). As was pointed out by Jewitt (2009), the ubiquitous presence of outgassing in icy minor bodies near the Sun, but its absence in the majority of such objects beyond the orbit of Saturn, is inconsistent with the concept of large deposits of condensed ice. If condensed ice was abundant, the icy minor bodies would display measurable activity beyond the orbit of Neptune. The conclusion must be that most currently available is trapped (occluded) within a host medium that does not release the unless a boundary near is crossed. Jewitt (2009) proposed that amorphous water ice is the most likely host, particularly because the sunward emission of inferred from the blue–shifted lines in 29P/Schwassmann–Wachmann 1 requires a host located sufficiently close to the surface to be sensitive to diurnal temperature variations. The same sunward emission is seen in 174P/Echeclus (Wierzchos et al., 2017). This excludes a deep source of condensed ice in combination with thermal lag effects as a viable explanation for the delayed onset of outgassing. The capability of amorphous water ice to trap hypervolatiles like , and the detailed conditions governing their release, have been well documented in laboratory experiments (e.g. Bar-Nun et al., 1985, 1987). The inclusion of amorphous water ice, crystallisation, and release of occluded species are standard features of advanced comet nucleus thermophysical codes since the seminal papers by Prialnik & Bar-Nun (1987, 1990) and used, e. g., to explain cometary outbursts. A systematic investigation of the level of amorphous water ice crystallisation at different distances and body latitudes, for a range of spin axis obliquities, was made by Guilbert-Lepoutre (2012). She finds that crystallisation is an efficient source of comet activity at –, and that crystallisation–driven activity beyond that distance is very unlikely. This reinforces the suggestion by Jewitt (2009) that most activity in Centaurs and distant comets may be governed by the crystallisation of amorphous water ice. This, in fact, has been the working assumption in comet science for decades (e. g., reviews by Prialnik et al., 2004, 2008), though empirical evidence has been lacking.
The surveys of activity in distant dynamically new or long–period comets by Sárneczky et al. (2016) and Kulyk et al. (2018) confirm that development of comae and even tails are common at – and that observationally confirmed activity beyond such distances is very rare. There are, however, notable exceptions. Comet C/2010 U3 (Boattini) was active at a record distance of – inbound with a dust production of , apparently keeping a similar level of activity all the way down to from the Sun (Hui et al., 2019). Comet C/2017 K2 (PANSTARRS) was active at inbound, with an estimated dust production rate of at (Jewitt et al., 2017). Hui et al. (2018) estimate the dust production rate for the – segment of the orbit to . Evidently, Comet C/2017 K2 is producing at least two orders of magnitude more dust than Comet Boattini at comparable heliocentric distances. Thermophysical modelling by Meech et al. (2017b), including an estimate of the dust production rate and the associated coma brightness, shows that the observed magnitude of C/2017 K2 from down to can be reproduced by the model if the hypervolatile is driving activity. If the substantially more stable supervolatile is considered instead, Meech et al. (2017b) find that the predicted brightness at would be fainter than the observed one, which corresponds to a two orders of magnitude smaller dust production rate. Therefore, –driven distant activity is potentially relevant for Comet Boattini, but not for Comet C/2017 K2.
Although the activity of Comet C/2017 K2 may indeed be driven by the sublimation of condensed , there are alternative options. One possibility is that is partially stored by an additional host medium that releases the at lower temperatures than possible for amorphous water ice. Based on comet abundances (Bockelée-Morvan et al., 2004), the most likely host besides amorphous water ice would be condensed ice. Observations of protostars show that roughly of the is intimately mixed with water ice, but that the remainder forms a separate system (Pontoppidan et al., 2008). Considering that the ESA Rosetta mission demonstrated that a fraction of ices in Comet 67P/Churyumov–Gerasimenko appear to have presolar origin (Calmonte et al., 2016; Marty et al., 2017), icy minor Solar System bodies could have inherited mixtures from the interstellar medium. Laboratory experiments show that readily traps hypervolatiles such as (Luna et al., 2008), (Satorre et al., 2009), and (Simon et al., 2019). Such substances are released from the ice above their individual sublimation temperatures (the process is called segregation or distillation), primarily at a temperature interval of –. As pointed out by Jewitt et al. (2017), the activity of Comet C/2017 K2 took place at a heliocentric distance where a nucleus temperature of – is expected. This makes segregation of out from a host a relevant process. Furthermore, Luspay-Kuti et al. (2015) demonstrated that the production rates of and were correlated with that of in Comet 67P/Churyumov–Gerasimenko, and distinctively different from the production of water. Building on that study, Gasc et al. (2017) showed that the production rates of , , and also correlated with , while was correlated with water. This led Gasc et al. (2017) to propose the presence of two main phases of separate condensed ices – and – each of them trapping a specific set of more volatile species.
The possibility that segregation of mixtures is responsible for activity at extreme distances, as well as its potential role played in the activation distance of most comets, needs to be explored further. However, to address this issue it is first necessary to understand the level of segregation that may have taken place already in the Primordial Disk. This inspires to the second question.
-
Q#2. What degree of release from entrapment is expected to take place in the Primordial Disk? Specifically, at what depth is the segregation front located?
As previously mentioned, the objects in the Primordial Disk at – could temporarily have been subjected to protosolar radiation equivalent to being located at – in the current Solar System. As seen from the works by Jewitt (2009) and Guilbert-Lepoutre (2012) this coincides with the conditions where crystallisation of amorphous water ice is expected to be an important driver of cometary activity. The observations of Comet C/2015 ER61 (PANSTARRS) at down to by Meech et al. (2017a), and their thermophysical modelling efforts to reproduce the observed coma brightness, indicate that the activity of this particular object is driven by sublimation of condensed ice. This inspires to the last question considered in the current work.
-
Q#3. At what depths are the sublimation front and the crystallisation front of amorphous water ice expected to be located after thermal processing in the Primordial Disk?
An additional key goal of the current paper is to make the first in–depth presentation and application of the thermophysics modelling code “Numerical Icy Minor Body evolUtion Simulator”, or NIMBUS, written by the author. It has previously only been briefly mentioned and applied by Davidsson et al. (2021b) and Masiero et al. (2021).
The rest of the paper is organised as follows. The thermophysical code is described in Sec. 2. Readers who wish to fast–forward to the results are encouraged to first read the Sec. 2 introduction, Sec. 2.1, and Sec. 2.12, for a minimum summary of the model. The model parameters used for the simulations in this paper are motivated in Sec. 3, the results are presented in Sec. 4, and they are discussed in Sec. 5. Finally, the conclusions are summarised in Sec. 6.
2 The thermophysical model
The in situ exploration of comet nuclei, initiated when spacecraft from Europe, the Soviet Union, and Japan explored Comet 1P/Halley in 1986, sparked a strong interest in the thermophysics of porous and sublimating ice/dust mixtures. It intensified during the years preparation for the ESA cornerstone mission Rosetta to Comet 67P/Churyumov–Gerasimenko in 2014–2016. A large number of thermophysical models were developed and were used to systematically explore the importance and behaviour of various physical processes in comets (e.g. Weissman & Kieffer, 1981; Fanale & Salvail, 1984; Kührt, 1984; Rickman & Fernández, 1986; Prialnik et al., 1987; Prialnik & Bar-Nun, 1987, 1988, 1990; Mekler et al., 1990; Rickman et al., 1990; Espinasse et al., 1991; Kömle & Dettleff, 1991; Kömle & Steiner, 1992; Prialnik, 1992; Kossacki et al., 1994; Kührt & Keller, 1994; Tancredi et al., 1994; Benkhoff & Huebner, 1995; Orosei et al., 1995; Prialnik & Podolak, 1995; Skorov & Rickman, 1995; Capria et al., 1996, 2001; Benkhoff & Boice, 1996; Podolak & Prialnik, 1996; Enzian et al., 1997, 1999; Shoshani et al., 1997; Markiewicz et al., 1998; De Sanctis et al., 1999; Skorov et al., 1999). In parallel to model development, numerous laboratory experiments were carried out to further the understanding of thermophysical evolution and to validate models (e.g. Bar-Nun et al., 1985, 1987; Benkhoff & Spohn, 1991; Grün et al., 1991, 1993; Benkhoff et al., 1995; Kossacki et al., 1997). The appearances of bright Comets C/1996 B2 (Hyakutake) and C/1995 O1 (Hale–Bopp), the discovery of the Kuiper belt (Jewitt et al., 1992; Jewitt & Luu, 1993), and NASA flybys of Comets 19P/Borrelly, 81P/Wild 2, 9P/Tempel 1, and 103P/Hartley 2 in 2001–2010 triggered further development of models, particularly considering larger body sizes, where gravitational compression and radiogenic heating become important (e.g Gutiérrez et al., 2000, 2001; De Sanctis et al., 2001, 2007; Capria et al., 2002; Capria et al., 2009; Choi et al., 2002; Davidsson & Skorov, 2002; Davidsson & Gutiérrez, 2004, 2005, 2006; Merk & Prialnik, 2003, 2006; Podolak & Prialnik, 2006; González et al., 2008, 2014; Rosenberg & Prialnik, 2009, 2010; Marboeuf et al., 2010; Marboeuf et al., 2011; Prialnik et al., 2008; Sarid & Prialnik, 2009; Gortsas et al., 2011; Guilbert-Lepoutre et al., 2011). Work to understand the thermophysics of the Rosetta target comet is ongoing (e.g. De Sanctis et al., 2015; Keller et al., 2015; Kossacki et al., 2015; Hu et al., 2017; Höfner et al., 2017; Mousis et al., 2017; Hoang et al., 2020; Davidsson et al., 2021b; Davidsson et al., 2021a), and has partially challenged previous ideas about comet physics. This primarily concerns ongoing discussions about dust mantle fragmentation and dust coma formation (Gundlach et al., 2015; Jewitt et al., 2019), as well as the storage medium for hypervolatiles (Gasc et al., 2017). The inspiration to construct NIMBUS came from this vast body of work, and the code should be considered an independent implementation of what may be considered a “standard comet nucleus thermophysical model”. It contains most, though perhaps not all, physical processes that are common to state–of–the–art thermophysical models. NIMBUS relies on and benefits from decades of efforts in the scientific community to better understand the thermophysics of porous media consisting of refractories and various types of ice. The only feature that may be considered new is the inclusion of as a host of more volatile substances (in this first implementation, exclusively CO) and the associated segregation process. The current paper attempts to describe NIMBUS in quite some technical detail, with the hope that this will benefit future users222It is the intention of the author to make the NIMBUS source code publicly available at GitHub..
NIMBUS considers a spherical body of radius consisting of a porous mixture of refractories and different types of ice. Currently, NIMBUS includes , , and but more species can be added in the future. The water ice may be amorphous and during heating it first transforms to cubic ice and then to hexagonal ice. Amorphous water ice may trap other species that are released gradually in three stages – during crystallisation, the cubic–hexagonal transition, and sublimation. Also may trap other species except water, that are released through segregation prior to the onset of sublimation. All other volatiles may be condensed333The term “condensed” is here used for pure ice substances. Because they are not trapped in another medium, they are free to sublimate when the temperature is sufficiently high. NIMBUS considers all vapour–to–ice transitions to form clean ices, which is why the term “condensed” is used interchangeably with “pure” or “free” ice. or trapped in or in some proportion.
NIMBUS makes no explicit assumptions about the geometric arrangements of matter on small local scales, e. g., whether it is monolithic (with or without holes), consists of evenly distributed grains, or forms a lumpy medium where grains locally are concentrated into larger pebbles. However, a particular geometric model may be implicitly applied by the model physical relations chosen in order to deal with, e. g., the porosity–correction of heat conductivity, the calculation of gas diffusivity, or evaluation of volume mass production rates. Such model physical relations are retrieved through function calls, and those functions may contain several different versions of such relations, published by different authors (and more can be added over time). It is up to the user to select the exact model physical relation to be applied, via an identification number that is submitted to the function. Therefore, the implied geometric arrangement of matter may change depending on the choice of model physical relation. This paper describes the relations chosen for the current simulations, but they are not hard–wired in NIMBUS.
NIMBUS tracks the conduction of heat and the diffusion of vapour both radially and latitudinally. Heat is transported through solid–state conduction, by radiative transfer, and during gas diffusion (advection). Sources of heat include solar irradiation, short– and long–lived radionuclei, energy release during crystallisation of amorphous water ice, and latent energy released when vapour is recondensing. Sinks of heat include thermal infrared radiative loss of energy from the surface into space, and consumption of energy during ice sublimation, segregation (here, escaping entrapment in ice), and crystallisation (here, and released from amorphous water ice). The model body can be treated as a slow rotator (to fully account for daily variations in illumination conditions), or as a fast rotator (where each considered latitude slab receives the diurnally averaged radiative flux valid for that latitude). Any elliptic, parabolic, or hyperbolic orbit can be considered, and any spin axis orientation may be applied. Orbital elements and spin angles can change with time, allowing for, e. g., consideration of spin axis precession and nutation.
The following sub–sections describe NIMBUS in detail. Specifically, they deal with governing equations and boundary conditions (Sec. 2.1), composition and porosity (Sec. 2.2), heat capacity and heat conductivity (Sec. 2.3), volume mass production rates and latent heats (Sec. 2.4), crystallisation and the cubic–hexagonal transition (Sec. 2.5), release of and from sublimating (Sec. 2.6), release of from (Sec. 2.7), gas diffusion (Sec. 2.8), radiogenic heating (Sec. 2.9), erosion (Sec. 2.10), the spatial grid, temporal resolution, orbit, and spin (Sec. 2.11), illumination conditions (Sec. 2.12), and implementation (Sec. 2.13). Finally, verification of the code is described in Sec. 2.14.
2.1 Governing equations and boundary conditions
The coupled differential equations for conservation of energy (Eq. 1), vapour mass (Eq. 2) and ice mass (Eq. 3) govern the thermophysical evolution of the body. For the symbols used throughout this work, see Tables 1–3 (functions), Table 4 (parameters), and Tables 5–7 (constants). The current version of NIMBUS considers number of species for which the following labelling convention is applied: refractories (), amorphous water ice (), cubic water ice (), hexagonal/crystalline water ice (), carbon monoxide (), and carbon dioxide (). Future versions of NIMBUS will incorporate a larger number of species, but the current focus is on the top three in terms of abundance. Note that the index nominally is used to refer to species, but in the governing equations it is necessary to distinguish between species acting as hosts (denoted by ) and species that are guests (denoted by ) within hosts.
The energy conservation equation, formulated in a spherical geometry, reads
| (1) |
Dependencies on for all functions are implicit, and other dependencies on functions have been written explicitly (primarily to highlight coupling between Eqs. 1–2, omitting parameters and constants for clarity). From left to right the terms in Eq. (1) denote: 1) the rate of change of internal energy (in units ); 2) radial heat conduction; 3) latitudinal heat conduction; 4) consumption of energy by net volume sublimation or release of energy by net condensation; 5) energy release during phase transitions and energy consumption during escape of an occluded species from a host medium; 6) radial and latitudinal convection; 7) energy release by radioactive decay.
The mass conservation equation of vapour reads,
| (2) |
where . In Eq. (2) the terms denote: 1) the rate of change of vapour density within pore spaces; 2) radial gas diffusion; 3) latitudinal gas diffusion; 4) vapour production during sublimation or vapour removal during condensation; 5) release of vapour from a host ice.
The mass conservation equation for solid ices reads,
| (3) |
where . In Eq. (3) the terms are: 1) the rate of change of ice bulk density; 2) the rate of removal by sublimation or formation by condensation; 3) the rate by which the ice is formed through phase transition of a precursor, or is destroyed by production of a successor.
| Symbol | Description | Unit |
|---|---|---|
| Matrix, analytic gas diffusion method | ||
| Area of cell and interface | ||
| Area of cell interface, generic | ||
| , | Area if northern/southern cell interfaces | |
| Total grain surface area within a volume | ||
| Initial molar fraction of in uranium | ||
| Initial molar fraction of in uranium | ||
| Current number of potassium | ||
| atoms in sample | ||
| Initial–to–current molar ratio of | ||
| atoms in sample | ||
| Diffusivity of in ice | ||
| Change of cell internal energy | ||
| Mass ratio: released guest | ||
| to transformed host | ||
| Initial effect of radionuclide decay | Table 12 | |
| Energy release during phase transition | ||
| Latent heat of sublimation of species | ||
| Luminosity factor | ||
| Water ice mass per cell | kg | |
| Mass of species per cell | kg | |
| External pressure | ||
| Gravitational pressure | ||
| Compressive strength | ||
| Mass loss rate of species to space | ||
| Radial cell boundary | ||
| Total radiogenic energy release | ||
| Solar flux at | ||
| Temperature | K | |
| Surface temperature | ||
| , | Temperature change | |
| Total mass within a grain | ||
| Cell volume | ||
| Initial mass fraction of radioactive | ||
| isotope in chondrite | Table 12 | |
| mass flux in | ||
| , | Zenith function |
| Symbol | Description | Unit |
|---|---|---|
| Final core cell thickness | ||
| Vector, analytic gas diffusion method | ||
| Specific heat capacity of mixture | ||
| Specific heat capacity of species | ||
| , , | Initial, intermediate, final common ratio | |
| Distance between cells and | ||
| Volumetric fraction of refractories | ||
| Volumetric fraction of volatiles | ||
| Specific heat capacity of vapour | ||
| Hertz factor | ||
| Mixing rate coefficient | ||
| Vapour number density of species | ||
| Number of radial cells | ||
| Pressure array for mantle cells | ||
| Fractal porosity | ||
| Partial gas pressure of species | Pa | |
| Pressure in cell | Pa | |
| Saturation pressure of species | Pa | |
| Volume vapour mass production | ||
| rate of species | ||
| Host mass transformation rate | ||
| of species | ||
| bulk density within ice | ||
| density within ice | ||
| density outside ice | ||
| Newton–Raphson function (Eq. 55) |
| Symbol | Description | Unit |
|---|---|---|
| Thermal inertia based on | ||
| Total heat conductivity | ||
| Solid–state heat conductivity | ||
| Refractories–to– ice mass ratio | ||
| Local bulk density | ||
| Global bulk density | ||
| Bulk density of condensed ice | ||
| Bulk density of hexagonal ice | ||
| Compact density of all ices | ||
| Bulk density of species | ||
| Density of metals | ||
| , | Density of minerals | |
| , | Density of refractories | |
| Bulk density of trapped in | ||
| within a large volume compared | ||
| to grain size | ||
| Transformation rate to/from | ||
| Radial mass flux rate of species | ||
| Porosity | ||
| Latitudinal mass flux rate of species |
| Symbol | Description | Unit |
|---|---|---|
| Targeted core cell thickness | ||
| Biot number | ||
| Co–declination | ∘ | |
| Surface cell thickness | ||
| Index of host species | ||
| Index of guest species | ||
| Newton–Rapson iteration index | ||
| Generic index of species | ||
| Latitudinal coordinate | rad | |
| Latitude cell index | ||
| 67P non–water vapour loss | ||
| 67P water loss | ||
| Radial cell index | ||
| Top cell label, analytic | ||
| gas diffusion method | ||
| Number of species | ||
| Rotation period | s | |
| Radial coordinate | m | |
| Body radius | m | |
| Index of radionuclide | ||
| Heliocentric distance | ||
| Time | s | |
| Disk clearing time | ||
| Time of local midday | s | |
| Radial coordinate within a grain | ||
| Polar angle within a grain | ||
| 67P dust–to–all–ice mass ratio | ||
| Coma dust–to– water vapour mass ratio | ||
| Azimuth coordinate within a grain |
| Symbol | Description | Value | Unit |
|---|---|---|---|
| Bolometric | 0.04 | ||
| Bond albedo | |||
| Crystallisation pre– | |||
| exponential factor | |||
| Crystallisation | |||
| activation energy | |||
| Cubic–to–hexagonal | |||
| transformation rate | |||
| Disk clearing speed | |||
| Segregation | |||
| activation energy | |||
| Energy release of | Table 11 | ||
| single isotope | |||
| decay | |||
| Newtonian gravity | |||
| constant | |||
| Crystallisation | |||
| energy release | |||
| Cubic–to–hexagonal | 0 | ||
| energy release | |||
| mass transfer | |||
| coefficient in | |||
| Pore length | |||
| Solar luminosity | |||
| 67P mass loss | |||
| molar mass | |||
| molar mass | |||
| molar mass | |||
| molar mass | |||
| Initial mass | 0.524 | ||
| fraction of | |||
| that is condensed | |||
| Initial mass | 0.717 | ||
| fraction of | |||
| that is condensed | |||
| Initial mass | 0.158 | ||
| fraction of | |||
| trapped in | |||
| Universal gas | 8.314510 | – | |
| constant | |||
| Solar constant | 1367 | ||
| Sky temperature | 0 | ||
| Time of perihelion | |||
| passage | |||
| Accretion impactor | 15 | ||
| velocity |
| Symbol | Description | Value | Unit |
|---|---|---|---|
| Semi–major axis | 23 | ||
| Grain surface area | |||
| Eccentricity | 0 | ||
| Initial mass | 1 | ||
| fraction of | |||
| that is amorphous | |||
| Current chondritic | Table 11 | ||
| mass fraction | |||
| of element | |||
| Chondritic mass | Table 11 | ||
| fraction of element | |||
| in dry rock prior | |||
| to aqueous alteration | |||
| Current molar fraction | Table 11 | ||
| of radionuclide | |||
| in element | |||
| Fraction of released | 1 | ||
| during crystallisation | |||
| Fraction of released | |||
| during the cubic–to– | 0 | ||
| hexagonal transition | |||
| Fraction of released | 1 | ||
| during crystallisation | |||
| Fraction of released | |||
| during the cubic–to– | 0 | ||
| hexagonal transition | |||
| Mass fraction of element | Table 11 | ||
| in refractories | |||
| Chondritic hydrogen | 0.021015 | ||
| mass fraction | |||
| Initial molar fraction | Table 11 | ||
| of radionuclide | |||
| in element | |||
| Chondrite water | 18.8 | ||
| mass fraction | |||
| Inclination | 0 | ∘ | |
| Boltzmann constant | |||
| molecule mass | kg | ||
| molecule mass | kg | ||
| molecule mass | kg | ||
| Isotopic mass of | Table 11 | ||
| radionuclide | |||
| Atomic mass unit | |||
| Perihelion distance | 23 | AU | |
| Grain radius | |||
| Pore radius | |||
| Ca–Al–rich | |||
| Inclusion | |||
| formation time | |||
| Half–life of | Table 11 | ||
| radionuclide | |||
| Solar System lifetime | |||
| NIMBUS time zero | |||
| Grain volume |
| Symbol | Description | Value | Unit |
|---|---|---|---|
| Spin axis right ascension | 0 | ∘ | |
| coefficient | Table 9 | ||
| Spin axis ecliptic latitude | ∘ | ||
| and coefficient | Table 9 | ||
| and coefficient | Table 9 | ||
| Spin axis declination | |||
| and coefficient | Table 9 | ||
| Emissivity | 0.9 | ||
| degrees of freedom | 6 | ||
| degrees of freedom | 5 | ||
| degrees of freedom | 5 | ||
| Spin axis ecliptic longitude | 0 | ∘ | |
| Initial refractories–to– | 4 | ||
| mass ratio | |||
| Initial ratio | 0.155 | ||
| by number | |||
| Initial ratio | 0.17 | ||
| by number | |||
| Argument of perihelion | 0 | ∘ | |
| Longitude of the | 0 | ∘ | |
| ascending node | |||
| Tortuosity | 1 | ||
| Iron density | |||
| Troilite density | |||
| Forsterite density | |||
| Accretion impactor | 300 | ||
| bulk density | |||
| Organics density | |||
| Refractories compact | 3000 | ||
| density | |||
| Amorphous | 500 | ||
| compact density | |||
| Cubic | 917 | ||
| compact density | |||
| Hexagonal | 917 | ||
| compact density | |||
| compact density | 1500 | ||
| compact density | 1500 | ||
| Stefan–Boltzmann | |||
| constant | |||
| Decay –folding scale | Table 11 | ||
| of radionuclide |
The surface boundary condition of Eq. (1),
| (4) |
envisions an infinitely thin layer responsible for absorbing solar radiation (and possibly cosmic background or Solar Nebula radiation) and emitting thermal radiation into space. The net difference between these fluxes equals a conductive energy flux that causes energy to either be added to, or removed from, an underlying layer of finite thickness (in practice the top grid cell). The common temperature of these two layers, , is not uniquely provided by Eq. (4), but is obtained from the combination of Eqs. (1) and (4) applied to the finite layer. is therefore determined through the joint effects of all processes captured in Eq. (1), such as sublimation of ices, radiogenic heating et cetera, in addition to solar heating and radiative cooling. Most thermophysical models of comets explicitly include a sublimation term in the boundary condition, thereby mixing fluxes over the boundary surface with one specific volumetric energy sink while ignoring other processes accounted for in those models (such as heating by crystallisation or radioactive decay). Equation (4) is a more stringently formulated boundary condition, that reinstates the expression found in older literature (e.g. Mekler et al., 1990; Prialnik, 1992).
2.2 Composition and porosity
The initial composition is specified by the mass ratio between refractories and water ice, the fraction of the water ice mass that is amorphous (assuming no initial cubic ice), the abundance of species relative that of water by number, and the densities at full compaction . The total mass of water ice in a cell of volume and porosity is given by
| (6) |
and the initial cell masses of species are therefore given by
| (7) |
Note that these values change with time because of the physical processes described in Eq. (1)–(3), and as a consequence, changes to are calculated as well. The bulk density of the cell (see Eq. 1) is consequently
| (8) |
Hyper– and super–volatiles () are initially distributed among different storage modes. In terms of mass fractions, is condensed ice, is trapped in ice, and the remainder is trapped in amorphous water ice. Note that ( is not considered trapping itself). If there is trapped in amorphous water ice, NIMBUS allows for a fraction to be released during crystallisation (i.e., a fraction is being transferred to the newly formed cubic water ice). During the cubic–to–hexagonal transition, a fraction of the original amount may be released, so that a fraction is released once the water starts sublimating. Similar mass fractions and are defined for that do not necessarily have to be equal to their counterparts. Water vapour that condenses is assumed to form pure crystalline water ice, even in the presence of and/or vapour, and regardless of temperature. During sublimation of water ice, clean crystalline water ice and “dirty” hexagonal ice (containing and/or remaining from the crystallisation and cubic–to–hexagonal processes) are assumed to be consumed in equal amounts, as far as possible.
Although the decay of radioactive species in reality has some effect on the overall mass budget, it is ignored in the context of calculating, i. e., the bulk density and heat capacity of the body. The amount of radioactive nuclei are only accounted for in terms of the heating resulting from their decay, see Sec. 2.9.
The porosity of the body can either be set to a constant initial value, or a radial porosity distribution can be calculated by solving the hydrostatic equilibrium equation (e.g. Henke et al., 2012) where the sum of the gravitational pressure
| (9) |
and a potential external pressure is balanced by the pressure at which a granular material just resists compression below porosity . At the relation is taken as an average of the one measured for silica particles by Güttler et al. (2009) according to the omnidirectional version of their Eq. (10), and the one measured for water ice particles by Lorek et al. (2016). The weight factors are taken as the volumetric fractions of ices,
| (10) |
and refractories, . At the function is taken as that measured by Yasui & Arakawa (2009) for a mixture at (see their Fig. 2).
The external pressure is set to the dynamic pressure
| (11) |
in order to mimic the compaction effect during accretion of impactors with density colliding at velocity during the formation of the body.
2.3 Heat capacity and heat conductivity
The specific heat capacities applied in the nominal version of NIMBUS are measured in the laboratory as functions of temperature. for refractories is represented by tabulated data for forsterite () as given by Robie et al. (1982). All water ice phases are represented by crystalline water ice with
| (12) |
as fitted by Klinger (1981) to data from Giauque & Stout (1936). For carbon monoxide ice I apply
| (13) |
from Tancredi et al. (1994) and for carbon dioxide ice () I apply the tabulated values from Giauque & Egan (1937) in their Table IV.
The specific heat capacity of vapour at constant volume is given by
| (14) |
applying the (translational and rotational) degrees of freedom and (vibrational modes are not excited at the relevant temperatures).
The total specific heat capacity (see Eq. 1) is given by the mass–weighted average,
| (15) |
The heat conductivities of pure and compact species are taken as follows: refractories have equal to the values tabulated for the H5 ordinary chondrite Wellman by Yomogida & Matsui (1983); amorphous water ice has
| (16) |
according to Kührt (1984) based on a model by Klinger (1980); cubic and hexagonal water ice has
| (17) |
according to Klinger (1980). Currently lacking better alternatives I here apply (see, e.g., Capria et al., 1996; Enzian et al., 1997) for similar approximations).
The heat conductivity of a granular medium is strongly dependent on its porosity. The current simulations use the “generation 2” model of Shoshany et al. (2002) to calculate the Hertz factor (the factor by which conductivity is reduced from that of compacted material due to porosity)
| (18) |
where the fractal porosity is given by
| (19) |
Note that this formula only is valid for . For the current simulations a ceiling was applied so that would default to a Hertz factor . The solid state component of the heat conductivity is defined as
| (20) |
i. e., a mass–weighted average conductivity is applied for compacted material (e.g. Prialnik, 1992), although there are alternative formulations (e.g. Enzian et al., 1997). The conductivity is modified continuously as the temperature, porosity, and composition change during simulations. For reference, Table 8 shows the thermal inertia corresponding to the solid–state component, for distinct values of for a mixture of dust and crystalline water ice.
| 720 | 790 | 990 | 1300 | 1500 | 1700 | 1900 | 2000 | 2100 | 2200 | |
| 440 | 490 | 610 | 770 | 950 | 1100 | 1200 | 1200 | 1300 | 1400 | |
| 160 | 170 | 220 | 270 | 330 | 370 | 410 | 440 | 460 | 480 | |
| 57 | 62 | 78 | 100 | 120 | 140 | 150 | 160 | 170 | 180 | |
| 8.3 | 9.1 | 11 | 14 | 18 | 20 | 22 | 23 | 24 | 26 |
NIMBUS also considers a contribution to heat conduction from radiative transfer (see, e.g. Kömle & Steiner, 1992; Enzian et al., 1997), thus the heat conductivity in Eq. (1) is given by
| (21) |
where is the pore radius. For , the radiative heat transfer component adds to the values in Table 8, peaking at the highest and . The radiative heat transfer component only becomes important at substantially higher temperatures that are relevant for dust mantles.
2.4 Volume mass production rates and latent heats
Condensed , , and will experience net sublimation whenever the local partial gas pressure is below the saturation pressure at the temperature in question. For a granular medium with grain radii the volume vapour mass production rate is given by
| (22) |
(Mekler et al., 1990; Prialnik, 1992; Tancredi et al., 1994) for . These production rates do not depend on the abundances of the species in question, i.e., multi–layer sublimation is assumed (a zeroth–order process in context of the Polanyi–Wigner equation, e.g., Suhasaria et al., 2017). Formally, Eq. (22) breaks down when only (sub)monolayers of species– molecules remain on the grain surfaces of more refractory elements. However, the errors introduced by not switching to first–order sublimation at extremely low concentrations have no practical implications on the solutions, neither from a mass nor an energy point of view. If the partial vapour pressure
| (23) |
is higher than the saturation pressure then and the vapour will experience net condensation. It is assumed that amorphous water ice and cubic water ice do not sublimate (). The saturation pressure of water vapour is certainly not zero below the temperatures at which crystallisation and the cubic–to–hexagonal transformation process are rapid ( and , respectively, see Sec. 2.5). Therefore, those species in principle do sublimate according to Eq. (22). However, the time scales of significant net sublimation are long when deviations from saturation pressure are small, as is the case except at the very surface. Therefore, the errors introduced by these simplifications are small.
The saturation pressures are functions of temperature according to laboratory measurements fitted by analytical functions on the form
| (24) |
with coefficients in Table 9 provided by Huebner et al. (2006) for and .
Data from Azreg-Aïnou (2005) are used for instead of Eq. (24).
2.5 Crystallisation and the cubic–hexagonal transition
I here define the functions and in Eqs. (1)–(2), as well as , , , and for in Eq. (1). See Table 10 for a summary of all and functions and values. Also and in Eq. (3) are defined.
Amorphous water ice crystallise into cubic water ice at a rate
| (26) |
where and according to Schmitt et al. (1989). The energy release during crystallisation is taken as (Ghormley, 1968). Part of this energy is being consumed by occluded species upon their release (see González et al. (2008) for a discussion on this topic).
Following Enzian et al. (1997), Eq. (1) assumes that this energy equals the latent heat during sublimation for the trapped species in question, explaining how the fifth term from the left in Eq. (1) is defined. At , where the crystallisation time–scale is , the applied latent heats are such that transition from an exothermic to an endothermic overall process takes place if the amorphous water ice contains or by number and the release is complete ( or ). However, crystallisation may still be exothermic at higher concentrations of trapped species as long as or .
The –values in Table 10 are obtained by calculating the mass ratio of guest to host, and are based on the requested abundances , with appropriate corrections for the partitioning of a species among hosts (, ) as well as the fraction of hosted species that are released during a given transition (, , , and ). Note that yields the volume production of a released guest species during the phase transition of its host (), so that yields the energy consumption with appropriate units () in Eq. (1).
| Host | guest and guest | |
|---|---|---|
| Amorphous | ||
| () | ||
| Cubic | ||
| () | ||
| Hexagonal | ||
| () | ||
| Carbon monoxide | ||
| () | ||
| Carbon dioxide | ||
No energy release is assumed to take place when cubic ice transforms to hexagonal ice (). For simplicity the transformation is assumed to proceed at a fixed rate when the temperature exceeds , thus
| (27) |
assuming , that corresponds to a half–life of cubic ice of . This process is endothermic because of the consumption of latent heat.
Finally, the definitions of and in Eq. (3) reflects the fact that amorphous water ice crystallises irreversibly at a rate , while equation for cubic ice has both a source term (formation at ) and a sink term (transformation to hexagonal water ice at ),
| (28) |
2.6 Release of and from sublimating
Crystalline water consists of two phases that are tracked separately by NIMBUS – pure condensed water ice, and (potentially) “dirty” hexagonal water ice. The former may be present from the start and/or forms when water vapour recondenses. It is pure in the sense that it does not host foreign species. The latter forms exclusively through the phase transition of cubic water ice, and may contain and/or . The combined sublimation or condensation rate of the two equals , and that contribution to the energy and mass conservation is fully accounted for by the fourth terms in Eqs. (1) and (2). For that reason, , i. e., the fifth term in Eq. (1) does not include energy consumption caused by the sublimation of the hexagonal water ice host (it has already been included in the fourth term). However, it is necessary to account for the energy consumption of its entrapped guests.
Let and be the local bulk densities of condensed and hexagonal water, respectively. Note, that . The transformation rate (equivalently, sublimation rate) of the hexagonal water ice host equals that of water itself,
| (29) |
because sublimation is here considered a zeroth–order process. The expressions for and in Table 10 specifies the guest/host mass ratio within hexagonal water ice, but also account for the fact that only a fraction of the available water contains trapped species. Therefore, provides the mass release () of () and ().
The ice mass conservation equations for and read,
| (30) |
The first two rows in Eq. (30) show that sublimation () reduces the amount of condensed water ice at a rate that is a fraction of , but that condensation () allocates the full water ice formation to the condensed phase. The two last rows in Eq. (30) show that the amount of hexagonal ice only increases during the cubic–to–hexagonal transformation (at rate ) because condensation () is not creating “dirty” hexagonal ice, while sublimation () proceeds at a rate that is a fraction of . Adding rows 1 and 3 (or 2 and 4), yields
| (31) |
Comparison with Eq. (3) shows that Eq. (31) is an identical statement and that
| (32) |
i. e., that water ice (the sum of condensed and hexagonal phases) form as a combination of recondensation and cubic–to–hexagonal transformation, and is being destroyed through sublimation.
2.7 Release of from
Because sublimation and condensation of are fully covered by the fourth terms in Eqs. (2) and (3), and because is not considered a host neither of itself nor of any other species, then . Because is assumed not to change structure and release energy during segregation (), and is not considered a host of itself (), the only remaining problem in the fifth terms of Eq. (1) and (2) is to define the segregation rate and .
In order to derive an expression for , consider a spherical grain of ice that contains at local mass density , where are spherical coordinates within the grain. The evolution of with time is governed by Fick’s second law,
| (33) |
where is the diffusivity (here referring to moving through ice). Integrating both sides of Eq. (33) over the grain volume, assuming a spherical–symmetric distribution of within the grain, defining the mass flux
| (34) |
and applying the Gauß theorem, yields
| (35) |
where is the total mass inside the grain, , and is the grain surface area. The mass flux across the grain surface can also be written
| (36) |
where is the mass transfer coefficient and is the partial density in the gas that surrounds the grain. The gas density is here assumed to be negligible in comparison with the concentration of molecules near the surface of the grain. If the Biot number is much smaller than unity,
| (37) |
then the mass loss is sufficiently slow to remove internal abundance gradients, so that everywhere, including at the surface (). Combining Eqs. (35)–(37) therefore yields,
| (38) |
Equation (38) can be generalised to a medium of volume containing of a large number of grains with combined surface area with a total mass of trapped ,
| (39) |
Denoting the bulk density of trapped within the entire volume by , and recognising that (the mass of released within unit volume and time), it yields
| (40) |
In order to proceed, must be evaluated. Applying Einstein’s relationship (Cooke et al., 2018),
| (41) |
and recognising that the mixing rate coefficient in the segregation process follows the Arrhenius equation according to laboratory experiments (e.g. Öberg et al., 2009),
| (42) |
Eq. (40) transforms to
| (43) |
where is the proportionality constant (pre–exponential factor) in Eq. (42) times , and is the activation energy of segregation. Recognising that the mass ratio between the trapped and the host is
| (44) |
the expression for the segregation rate used by NIMBUS consequently is,
| (45) |
Evaluating and is a delicate problem. Laboratory studies of the behaviour of mixtures have been initiated fairly recently. Cooke et al. (2018) studied the segregation of out of at low temperatures (–) and found a low activation energy , suggesting efficient diffusion along pore surfaces rather than bulk diffusion. However, the conditions within icy minor Solar System bodies are different from the interstellar conditions studied by Cooke et al. (2018). Firstly, the relevant temperatures are higher and the porosity of ice decreases drastically with increasing temperature (Satorre et al., 2009). The smaller presence of pores at higher temperatures may strongly reduce the diffusivity and result in a much higher –value. Secondly, the experiments are performed at ultrahigh vacuum in order to simulate interstellar conditions. However, the gas pressure in the interior of Solar System bodies may be significant, and the loss rate at a grain surface depends on the density difference between the grain and the surroundings (see Eq. 36). Although that difference was ignored in the derivation above, it could still have the effect of increasing the activation energy. Simon et al. (2019) studied the entrapment efficiency of in as function of temperature, ice layer thickness, and concentration during deposition. They found that entraps more efficiently than below the sublimation temperature. They also studied the desorption rate as function of temperature and found that released some near (corresponding to the sublimation temperature and likely representing desorption of deposited on the surface of the ice rather than authentic segregation). They found that low–level emission of took place at –, and that most was released during sublimation that peaked at . It is interesting to compare these results with similar experiments for (Luna et al., 2008) and (Satorre et al., 2009) entrapped in . Both species have a first desorption peak near (corresponding to and sublimation), a second broad desorption peak at – (attributed to the removal of porosity in the host), and a third peak centred at – during sublimation. The experiments of Luna et al. (2008) and Satorre et al. (2009) were performed at a comparably high pressure of which makes sublimate strongly at a temperature that was higher than in the experiment of Simon et al. (2019) that took place at . By doing so, Simon et al. (2019) potentially may have superimposed the two release peaks corresponding to densification and sublimation. It remains to be seen if displays three–peak segregation behaviour as do and in segregation experiments at higher pressure. The pressure within Solar System minor icy bodies may become substantially higher than in these experiments. For example, the average pressure in the inner of a Hale–Bopp–sized model body during of steady–state in the early Solar System is (see Sec. 4.2), while the average temperature is (i. e., at such a pressure the temperature must reach even higher, to roughly – to cause strong sublimation). In order to provide –parameters relevant for modelling segregation within Solar System bodies, such experiments should be performed at substantially higher pressures than currently available.
In the meanwhile, a simple and admittedly arbitrary method of estimating values was employed to enable the test simulations of this paper. Defining as the timescale of complete segregation (), inserting that into Eq. (43), and logarithmising that equation, yields a linear relation with slope , having on the x–axis and on the y–axis. Fixing two pairs of is therefore sufficient to determine the corresponding values. If currently is present in icy minor Solar System bodies, such ice must have survived the intense radiation of the protosun. Postulating a clearing–time of (at which the solar luminosity was times the current value, corresponding to a radiative equilibrium temperature of ), and demanding a segregation time–scale of in such conditions provided the first pair . Considering that the second segregation peak is near , it is reasonable to attempt defining a time–scale of segregation at that temperature. Simon et al. (2019) find that segregation propagates at monolayers per minute, corresponding to full segregation in for a thick ice layer ( monolayers). However, that time–scale was measured at ultrahigh vacuum, and may very well correspond to sublimation rather than segregation. The interior of a Solar System body would be close to saturation, having a pressure at that exceeds that at by a factor according to Eq. (24) and Table 9. Assuming that segregation is slower by a corresponding factor, it yields the second pair . The corresponding Eq. (45) parameters become .
2.8 Gas diffusion
The mass flux of vapour flowing through a porous medium depends on local temperature and pressure gradients, as well as on geometric parameters used to calculate the diffusivity (here referring to gas moving through the porous mixture of refractories and ices in the body interior). NIMBUS applies the Clausing formula (e.g. Skorov & Rickman, 1995; Davidsson & Skorov, 2002)
| (46) |
It is here envisioned that gas flows in cylindrical tubes of radius and length . The formulation in Eq. (46) is accurate for any ratio (Skorov & Rickman, 1995), whereas expressions used in some thermophysical models only hold for infinitely long tubes. The tortuosity is the ratio of the path length through a series of tubes with different orientation, and the net path passed along the vertical. A similar expression is used for the latitudinal mass flux , except that the differential of has to be replaced by .
2.9 Radiogenic heating
It is possible to include the effect of radiogenic heating in NIMBUS, with the caveat that it does not allow for melting of volatiles or dust. The code currently includes the long–lived radionuclides , , , and , as well as the short–lived radionuclides and . Although the simulations in the current work only includes the heating of long–lived radionuclides, the nominal abundances of short–lived radionuclides are here discussed as well for future reference.
NIMBUS uses four different parameters to define the abundances and physical properties of the radionuclides. These are: 1) the fraction of the refractories mass consisting of the chemical element in question, ; 2) the fraction by number of a given element that consisted of the radioactive isotope at after CAI, ; 3) the half–life of the radioactive isotope, ; 4) the energy released per atom during decay . The index identifies the radionuclide as follows: (), (), (), (), (), and ().
The default abundances of the radionuclides are obtained as follows. The current elemental mass fractions in CI carbonaceous chondrites are taken from Table 3 in Lodders (2003), see Table 11. It is assumed that such material, prior to water melting and aqueous alteration of the CI chondrite parent body (that integrated water into the rock phase), consisted of a mass fraction of water ice and a fraction of dry rock. In order to calculate the elemental mass fractions in the dry rock component of the CI carbonaceous chondrite parent body (that is the quantity relevant for NIMBUS) it is recognised that the current CI chondrite hydrogen mass fraction is (Lodders, 2003) and it is assumed that all this hydrogen once was locked up in the water ice. The water mass fraction is then
| (47) |
Garenne et al. (2014) performed thermogravimetric analysis of several carbonaceous chondrites, where the samples are heated and the amount of water loss from the rock phase is measured. They find for CI chondrite Orgueil, an average for three CR chondrites, and an average for 16 CM chondrites. This indicates that Eq. (47) is a reasonable estimate. The dry–rock mass fractions are given in Table 11.
The current molar fractions of the isotope to that of their chemical elements in Table 11 are taken from Table 6 in Lodders (2003). Note that the entries for and are zero because those radioisotopes are currently extinct. At after CAI those isotopes had fractional abundances by number of for and for (Henke et al., 2012). The other –values are calculated from in the following. However, first it is necessary to specify the half–life that is related to the –folding scale ,
| (48) |
that are taken from Henke et al. (2012) except the value that is from Norris et al. (1983). The values are given in Table 11.
| Quantity | ||||||
|---|---|---|---|---|---|---|
| 0 | 0 | 1 | ||||
| 0.705 | 2.6 | 1200 | 14000 | 690 | 4500 | |
| 1 | 0.2558 | 0.7442 | ||||
| 25.986891692 | 59.934071683 | 39.963998475 | 232.038055325 | 235.043929918 | 238.050788247 |
Let be the current number of potassium atoms in a sample. Then the current number of isotopes is and the current number of and isotopes is . The latter number of isotopes has remained constant over time, while the number of isotopes at was where is the age of the Solar System. Then the fraction at is obtained as
| (49) |
For uranium I proceed in a similar way, and ignore isotopes other than and (the third most common isotope, , constitutes merely by weight in uranium; Goldin et al., 1949). Defining and as the original number of the isotopes in what currently is of uranium, it is found that the primordial molar fractions were,
| (50) |
The larger abundances of chemical elements in the past, owing to the presence of radioactive isotopes that now are partially or fully gone, are calculated by equating the current mass of stable isotopes in chondritic rock, with that in the past, . This assumes that the total chondritic rock mass has not changed, which is not entirely true (e. g., a nucleus is somewhat heavier than its daughter), but such minor modifications are here ignored. It also equates isotope molar and mass fractions, which is acceptable given other uncertainties, because the molecular mass differences among isotopes of the same element are small. If so,
| (51) |
This works for aluminium, iron, and potassium (), that do have stable isotopes. As seen in Table 11, the differences between and are very small, simply because the radioactive isotope fractions among those elements were always small. For thorium and uranium that do not have stable nuclei, the elemental fractions are calculated as
| (52) |
That explains the values shown in Table 11. Finally, Table 11 also lists the energy being released when one isotope decays, with values taken from Henke et al. (2012) except for taken from Castillo-Rogez et al. (2009).
The literature often provides the applied mass fraction of the radionuclide in the whole rock mass at (e.g. De Sanctis et al., 2001; Choi et al., 2002; Guilbert-Lepoutre et al., 2011). Those values are calculated for the default NIMBUS composition and presented in Table 12 to ease a comparison. Table 12 also provides the effect produced at , where is the energy provided if worth of the radioactive nuclide would decay at once ( is here in seconds). This quantity is often reported in the literature as well (e.g. Henke et al., 2012).
| Radionuclide | ||
|---|---|---|
The default NIMBUS values for , , and are intermediate between those of De Sanctis et al. (2001) and Choi et al. (2002) and closer to the latter. The Guilbert-Lepoutre et al. (2011) values are somewhat higher still. For the NIMBUS are higher than both De Sanctis et al. (2001) and Choi et al. (2002), but not as high as for Guilbert-Lepoutre et al. (2011). The default NIMBUS -values are somewhat higher for , , and but somewhat lower for and compared to Henke et al. (2012).
2.10 Erosion
If the near–surface gas pressure becomes sufficiently strong to overcome the local tensile strength, solid material will be torn off and ejected into the coma. The nominal version of NIMBUS does not take such erosion into account. However, an alternative version called NIMBUSD has been developed (“D” for “dust”, although this version of the code allows for erosion that ejects ices as well as refractories, that may contribute to an extended source of vapour in the coma). Because the different latitudes will evolve at different rates, the initially spherical nucleus will start to deform into a non–spherical shape. NIMBUSD keeps the original radial cell boundaries on a normalised axis (with zero at the centre and unity at the surface), but re–labels the cell boundary depths in response to the erosion. Physical quantities are interpolated to comply with the new grid after each “erosion campaign” (that typically is implemented when a certain fraction, e. g., 10% or 50% of the top cell has been eroded). Because the erosion progresses at different speeds at each latitude, neighbouring slabs eventually get out of phase with each other. For that reason, the latitudinal diffusions of heat and vapour are switched off in NIMBUSD. Heavily eroding objects (e. g., comets passing close to the Sun) are typically modelled on time–scales that are sufficiently short (a single parabolic passage, or a handful of apparitions for a short–period comet), that latitudinal heat and gas flows are negligible. NIMBUSD allows the user to define types of erosion rates, including: 1) observed rates stored on file as function of heliocentric distance; 2) a certain amount of solids is eroded in some proportion to the current outgassing rate of vapours; 3) a certain amount of solids is eroded in some proportion to the peak sub–surface net force caused by the gas pressure.
2.11 Spatial grid, temporal resolution, orbit, and spin
The volume of the body is divided into grid cells. The radial distribution of grid cell boundaries is generated with various methods depending on the application, including geometric progression that allows for high spatial resolution near the surface, and equi–distant or equal–volume cells when coarser resolution is desirable. The latitudinal distribution of grid cell boundaries is uniform in or depending on application. Once all cell boundaries have been set up, all cell volumes and the surface areas of the four walls associated with each cell are calculated analytically.
As an example, the grid–generation is here described for the geometric progression (radially) and equal–angle (latitudinally) case. The properties of the radial discretisation is determined by specifying the body radius , the thickness of the outermost cell and a targeted innermost cell thickness . The number of cells are calculated as
| (54) |
where is an initial guess of the common ratio. The Newton–Raphson method is then used in order to find a more exact value for the common ratio, and the final inner most cell thickness ,
| (55) |
which is iterated until convergence is reached, upon which is evaluated. The locations of the outer boundaries of the cells are obtained as
| (56) |
For example, , , and yields , and .
If the desired number of latitudinal slabs is , the equal–angle variant trivially has northern boundaries at with and , in terms of the angular distance from the north pole. A cell with label radially and latitudinally has lower and upper surface areas given by
| (57) |
with and . Because NIMBUS is not geometrically three–dimensional, an arbitrary longitudinal width of a cell can be assigned. Here, that width is the full , which means that the cell shape resembles a torus. The northern and southern surface areas are given by
| (58) |
The volume of the cell is given by
| (59) |
As described in more detail in Sec. 2.13, NIMBUS initially assigns species masses and internal energies to each cell, and then essentially operates by injecting new energy due to radioactivity and absorption of solar energy, and by moving mass and energy between cells (including escape of mass and radiative energy loss across the outer surface). Because the body is not modelled in three dimensions, special conditions prevail in cells that contain the polar axis (including all core cells). A NIMBUS core cell (with ) only trades energy and mass with the core cells just north and south of it, and with the cell just exterior to it. All these cells share the same local hour. A real core would be simultaneously influenced by exchange with exterior cells having the full variety of local hours. However, as long as the time–scale for conducting energy from the surface to the core is much longer than the rotation period, the error in the core temperature will be negligible. Cells with and trade energy and mass with cells immediately under and above, as well as with the cell just south of it (or, if , just north of it). Again, all these cells have the same local hour, but on a real body they have physical proximity to regions with very different local hours. In such regions, the quality of NIMBUS calculations would only be compromised if longitudinal lateral flows of energy and mass cannot be ignored. For these reasons, extreme slow–rotators should not be considered with NIMBUS, because they require a full three–dimensional treatment.
NIMBUS uses a dynamical time step that adjusts to the prevailing conditions. The dynamical time steps are controlled by forcing the change of internal energy at each step to fall within a certain percentile range, by making sure that the changes in gas partial pressures remain within certain limits, and that net vapour mass removal is not so large that saturation conditions cannot be re–established (because of insufficient amount of ice).
The orbit of the body is described by the standard elements: semi–major axis , eccentricity , inclination , argument of perihelion , the longitude of the ascending node , and the time of the perihelion passage (in case , then is replaced by the perihelion distance ). The orientation of the rotational axis of the body is given by the right ascension and declination in the equatorial system, or by the longitude and latitude in the ecliptic system. All these parameters can be considered functions of time in order to account for changes to the orbit and spin state during long–term simulations.
2.12 Illumination conditions
The solar flux at equals the solar constant when modelling the contemporary Solar System. However, when modelling the early Solar System NIMBUS accounts for changes of the protosolar luminosity according to the Hayashi– and Henyey–tracks presented by Palla & Stahler (1993), , where the luminosity factor is a dimensionless number that should be multiplied to the current solar luminosity to obtain the luminosity at a given moment. Specifically, on the “birth–line”, that marks the end of main accretion and the start of the Solar Nebula. This is considered or “time zero” in NIMBUS. It is likely that the formation of Calcium–Aluminium–rich Inclusions (CAI) at , generally considered the first Solar System solids that define the age of the Solar System (MacPherson et al., 1995), took place during the protostellar Class 0 stage (Tscharnuter et al., 2009) because of its unique capacity of generating temperatures sufficiently high for CAI formation. If so, there is a time difference of order (André & Montmerle, 1994). For simplicity, time is here measured with respect to instead of .
The luminosity factor falls to in the first and to in the first . At the end of the Solar Nebula lifetime (taken as based on the mean lifetime of gas disks in foreign solar–mass systems (Zuckerman et al., 1995; Haisch et al., 2001; Sicilia-Aguilar et al., 2006) the luminosity factor falls to . Therefore, the time period during which giant planet formation and significant planetesimal growth in the outer Solar System takes place is characterised by substantially higher luminosities than in the current Solar System. The luminosity factor continues to drop and reaches a minimum of at . At that point, marking the transition from the Hayashi–track to the Henyey–track, the luminosity factor starts to increase anew, reaching at the onset of hydrogen fusion at . It is assumed that the luminosity factor then increases linearly to the current value in the following .
Prior to substantial dust agglomeration and planetesimal formation the disk midplane is likely to be shielded by the opaque disk. Thus, an option exists to reduce the illumination to a low level at and increase it exponentially to the unshielded value at a disk clearing time ,
| (60) |
where the coefficient determines the rate of change. If planetesimal formation is as rapid as suggested by models of the gravitational collapse of pebble swarms (Nesvorný et al., 2010) formed in streaming instabilities (Youdin & Goodman, 2005; Johansen et al., 2007) the clearing time could be as short as a few times . However, if the metallicity needs to increase to a certain value (through gas removal from the disk) in order for streaming instabilities to become efficient (Johansen & Mac Low, 2009), could be substantially larger. Slower planetesimal growth and a longer clearing time are also predicted in other planetesimal formation scenarios (Weidenschilling, 1997; Davidsson et al., 2016).
NIMBUS can operate with time–resolved rotation, where the cosine of the solar zenith angle is
| (61) |
where
| (62) |
Alternatively, NIMBUS operates in fast–rotator mode were insolation is rotationally averaged for each latitude,
| (63) |
Equation (61) is suitable for studying the detailed thermophysics of bodies on time–scales corresponding to a few orbits, while Eq. (63) enables simulations extending over – timescales. Note that the co–declination (the angle between the positive spin vector and the direction to the Sun) is a function of time because of orbital motion and/or spin orientation changes. Note that , given by or Eq. (60), and that , given by Eq. (61) or Eq. (63), enter into Eq. (4). The same geometrical grid can be applied regardless if Eq. (61) or Eq. (63) is considered. Both provide a certain radiative flux at any given moment. The only difference is that Eq. (61) contains a (potentially strong) diurnal variation, while Eq. (63) provides a constant flux (if integrated over an rotational period the total fluxes are the same).
2.13 Implementation
The NIMBUS numerical scheme is explicit, i.e., the properties at a current time step are calculated based on knowledge of the properties at the previous time step. Cells are assumed to be isothermal, i. e., a single temperature is applied to the entire interior of a cell. It means that NIMBUS does not rely on temperature specification at grid points (cell interception points) and it does not attempt to interpolate temperatures between such grid points. Each grid cell carries a certain amount of mass for each species (solids and vapours are stored separately) and a certain internal energy, and those masses and energies are updated over time by considering local sources and sinks as well as exchange with neighbouring cells through diffusion and transport processes. Temperatures are calculated from the internal energies, considering the local heat capacity, i. e., if the energy change is , and the prevailing temperature is , the temperature change is obtained by solving the equation
| (64) |
Equation (64) can be written as a simple function that can be solved analytically for if approximating by a locally linear function, which is an acceptable approximation because is small compared to temperature ranges over which changes significantly and non–linearly. Exchanges of masses and energies between cells are calculated based on the temperature and pressure gradients prevailing between neighbouring pairs of cells. For example, solid–state heat conduction is calculated by applying the Fourier law: 1) using the temperature difference between the cells to obtain ; 2) using the distance between the cell centres to obtain ; 3) using the average heat conductivity values for the two cells to obtain ; 4) evaluating the energy exchange between the cells during a time step as , where is the surface area that the two cells have in common, through which the heat transfer takes place. The exact amount of energy removed from one cell is added to the other, leading to global energy conservation. The governing equations are fulfilled by modelling the underlying physical processes that give rise to these equations, rather than subjecting the equations themselves and their terms to mathematical manipulation (as would have been done in, e. g., the finite element method).
The code first evaluates all changes to masses of solids and vapour consistent with the current time step, using the latest known set of temperatures and pressures. The code moves on to evaluate energy changes only if the mass changes are consistent with the following requirements: 1) no cell should lose all its vapour when diffusion net losses exceed gains of vapour from sublimation, segregation, or crystallisation; 2) local changes in gas pressure exceeds a certain user–defined threshold. If the mass changes are unacceptable, the time step is reduced and the evaluation starts from the beginning. Reducing the time step primarily has the effect of not allowing net mass loss (at a given rate) to proceed for so long that a cell runs out of gas. If the mass changes are approved, the energy changes are evaluated and compared with user–defined criteria. If the percentile changes are above a certain maximum threshold, the time step is reduced and the entire process starts from the beginning. Repeating the mass–change evaluation for the new smaller time–step is necessary, because many mass–changing mechanisms (e. .g, net sublimation) lead to specific energy changes that need to be properly evaluated. If both the mass and energy changes have been approved, those changes are implemented and the code advances to the next time step. At this point, there is also a possibility to increase the time step, if the largest percentile energy change has fallen below a minimum threshold.
The amount of vapour in a grid cell is typically extremely small compared to the amount of the frozen volatile (often on the ppm–ppb level). The net gas diffusion fluxes are often large (particularly near the surface where temperature and pressure gradients are steep). Net mass loss at such high rates (if considered constant) could remove all available cell vapour in a very short time. Imposing the requirement that only a small percentage of the cell vapour mass can be removed before re–evaluation of the local conditions, would lead to extremely small time steps. One would find that the time–scale for re–establishing saturation conditions in the cell due to ice sublimation typically is even shorter, by orders of magnitude. An excessive amount of computing time would be spent on calculating small deviations from saturation pressure caused by gas diffusion, followed by guaranteed and practically instantaneous re–establishment of saturation conditions through ice sublimation (when ) or vapour condensation (when ). NIMBUS therefore assumes that cells that do contain ice of species have vapour of species at saturation pressure.
The magnitude of the mass flux rate of vapour between neighbouring cell pairs that both contain condensed ice , is determined by the difference in saturation pressure between the cells according to Eq. (46). As long as the two cells do not have the same temperature, the saturation pressure difference is non–zero. The amount of mass that a given cell loses or gains across a given boundary during the duration of the time–step is determined. This includes both radial and latitudinal flows. That mass is then subtracted from the losing cell and added to the gaining cell. This is done for all cells and all their mutual boundaries. The simple procedure of subtracting a given amount of mass from one cell and adding it to another guarantees mass conservation. At that point, most cells will no longer have vapour densities corresponding to saturation. The volume mass production rate (Eq. 22) is evaluated for the difference between saturation pressure and actual current pressure, whatever that difference might be. That provides a time–scale for re–establishment of saturation conditions. If this time–scale is shorter than the duration of the time–step, re–establishment of saturation condition is considered guaranteed. Re–introduction of saturation conditions in an under–dense cell means that the an amount of vapour is added to the gas phase (to reach saturation), and the same amount is subtracted from the ice phase. For an over–dense cell, the excess vapour is removed and the amount of ice increases accordingly. Again, the simple procedure of shifting specific masses between accounts guarantees mass conservation.
That allows for the usage of larger time steps where the total vapour mass change in a cell by outflow or inflow are calculated (typically being much larger than the current content of vapour), and where re–establishment of saturation conditions are made through a one–time larger net sublimation or condensation of ice. Such assumptions of local saturation conditions are common in thermophysical codes (e.g. Espinasse et al., 1993; Orosei et al., 1995; Capria et al., 1996; De Sanctis et al., 1999). Doing so, there is only a need to use small time steps when: 1) a cell is about to run out of ice (making sure the demands for vapour does not exceed the availability of ice); 2) when it is so cold (i. e., the sublimation is so slow) that the ice is not capable of providing the requested amount of vapour in the time available; 3) or, in case condensed ice is lacking, injection of vapour from segregation and crystallisation are not sufficiently fast to prevent full net vapour removal by diffusion.
Saturation conditions for vapour are generally not prevailing in regions where ice is not present. That includes the dust mantle that is devoid of all ices, and regions above the and sublimation fronts. In order to calculate accurate vapour pressures in such regions, without necessarily having to resort to the time–consuming short time steps described above, NIMBUS calculate those pressures analytically in the radial direction according to the method described in the following. The lateral flows have substantially smaller temperature and pressure gradients, which means that those flows can be calculated with the nominal numerical procedure (without requiring too small time steps). That is to say, the latitudinal flows are calculated numerically, and the radial flows above the sublimation front are calculated analytically.
Consider the uppermost cell that contains plenty of ice (constituting the sublimation front of that element), ignoring cells closer to the surface that might contain very small amounts of ice deposited by temporary recondensation. Label the cells with being the previously mentioned icy cell and being the surface cell. The method is applied to all latitudes and species but the indices are not written explicitly below for brevity. The cell temperatures are considered known (those have just been updated during the previous time step), but the gas pressures for (and the corresponding cell vapour masses ) are unknown. Cell is assumed to have local saturation pressure, . Additionally, a “ghost cell” in the coma just above the body surface is assumed to have . Steady–state is required, which means that the cell vapour masses do not change during the time step, i. e., that the inflow of mass to a given cell is exactly balanced by the outflow. However, the flow rate changes from one time–step to the next. Using Eq. (46), that mass transfer balance can be stated as
| (65) |
where and are the common cell wall area and distance between cells, respectively, with the cell pairs in question specified by the sub–scripts. The average porosity, e. g., , is applied for each cell pair, and because are assumed identical for all cells, those factors cancel on both sides. Equation (66) can be re–arranged in three different ways,
| (66) |
Now define an vector containing the unknown pressures. Also define a sparse matrix :
| (67) |
| (68) |
| (69) |
All elements in that have not been explicitly assigned in Eqs. (67)–(69) are set to zero. With formulated like this, the product equals the left–hand sides of the expressions in Eq. (66). Specifically, the two left–side terms in the first expression of Eq. (66) are used for and in Eq. (67), the three left–side terms in the second expression of Eq. (66) are used for , , and with in Eq. (68), and the two left–side terms in the third expression of Eq. (66) are used for and in Eq. (69). In order to complete the equations in Eq. (66), an vector needs to be defined, so that . It has a single non–zero entry,
| (70) |
corresponding to the right side of the first expression of Eq. (66), because the right–hand side of the second expression in Eq. (66) is explicitly zero (applied for elements in ), and as well, because it was decided to apply a zero pressure for the ghost cell, . Evidently, both and only contain geometric factors, porosities, and temperatures that are all known. It is therefore possible to calculate the unknown pressures,
| (71) |
Applying a zero pressure for space at the body/coma interface is standard in many thermophysical models (e.g. Mekler et al., 1990; Prialnik, 1992; De Sanctis et al., 1999; Capria et al., 2009). However, in reality the coma does have a pressure that is non–trivial to evaluate. This is because the near–surface coma is strongly deviating from thermodynamic equilibrium, i. e., the molecular distribution function is not given by the Maxwell–Boltzmann distribution function within the so–called Knudsen layer (e.g. Cercignani, 2000; Davidsson, 2008). Davidsson & Skorov (2004) solved the Boltzmann equation numerically, using the Direct Simulation Monte Carlo (DSMC) method, in order to calculate Knudsen layer properties for a range of body surface temperatures and near–surface temperature gradients. NIMBUS does not apply the near–surface density solutions of Davidsson & Skorov (2004) because those assumed presence of water ice up to the very surface, which is inconsistent with the presence of the dust mantle considered in Eq. (71). Furthermore, the work of Davidsson & Skorov (2004) only concerned water sublimation, while NIMBUS also includes and . However, NIMBUS can easily be upgraded to include a non–zero coma pressure by defining as the right–hand side of the third expression in Eq. (66) and use a proper value for for species obtained from Knudsen layer simulations.
Once the pressures have been calculated, the corresponding total amount of vapour mass integrated over cells 1– may have changed since the last time step. If the mass has increased, this mass enhancement is assumed to be provided by sublimating the corresponding amount of ice in cell . If the new vapour mass is lower than that currently present, the excess is vented to space and contributes to . In certain conditions, the calculated steady–state pressures may exceed the local saturation pressures. This typically happens soon after sunset, when the interior is still warm and the saturation pressure of cell is relatively high, while the surface cells are cooling off and their saturation pressures drop drastically with temperature. If so, NIMBUS allows for the frost formation in those cells, by locally condensing the excess mass ( is reduced accordingly). At sunrise, the analytically calculated gas pressure typically is below the local saturation pressure, which triggers sublimation of the near–surface frost (increasing by adding to the vapour arriving from larger depths).
This analytic method is not suitable for gas release due to crystallisation or segregation because those zones of release are thick compared to the single–cell sublimation fronts. A suitable cell is difficult to define, and pressures may be far from saturation. If vapour of species is released by segregation and/or crystallisation within a region above the sublimation front of ice, those contributions are added to the vapour mass radial profile set up by the analytical method. As long as those extra contributions keep the cell pressures below the local saturation pressures, they will only give rise to an elevated pressure compared to the pure analytical solution. In case the saturation pressure is exceeded, the excess will condense. If vapour of species is released by segregation and/or crystallisation below the sublimation front of ice, or if all condensed ice has been exhausted, that release takes place in an environment where no analytical solution exists (remembering that the validity of Eq. (71) is limited to a particular near–surface region that lacks ice, and that the approach loses meaning when all ice in the body has been exhausted). In such cases, the accurate numerical evaluation takes over.
If the local temperature and pressure gradients are shallow, comparably long time steps may still be afforded, because cells with a net loss rate would empty slowly. For example, radiogenically heated bodies typically have a large and almost isothermal core, while the temperature falls rapidly near the surface. If the core temperature is high enough to allow for release by slow crystallisation and ice sublimation, such vapour will diffuse upwards at a gentle pace. Most will recondense near the surface where it is sufficiently cold for freezing, and before the temperature and pressure gradients have become so steep that diffusion calculations would require a substantial reduction of the time step. In such situations simulations will proceed relatively fast. However, if a substantially more volatile species like is released during crystallisation in the same scenario, it will not encounter any region cold enough for condensation, and NIMBUS would have to deal with its diffusion across the very surface, where temperature and pressure gradients might be very steep. In such cases, the time step might become very small. In one numerical experiment, a body was studied right after the last amounts of condensed ice at its core were exhausted. The time it took to evacuate the gas that filled all pores from the core to the surface was measured. That complete removal took 6.8 years in the simulated world, and costed of CPU time. Prior to the removal of ice, the model object was advanced at the same cost. Because such drastic slow–downs may be calculationally prohibitive, NIMBUS comes with the option to omit detailed calculations of diffusion of vapour when ice is not present. Specifically, when segregation and/or crystallisation are solely responsible for production of vapour, gas may be vented directly to space. Such a simplification is acceptable because: 1) during multi–kyr simulations, a diffusion time–scale of a few years can be considered instantaneous and there is little value in resolving it temporally (on the minute level); 2) the release mechanisms do not depend on the vapour pressure of element , thus gradual diffusion or immediate venting does not change the way vapour is being produced; 3) all energy release and consumption associated with the release of vapour are still calculated accurately, except that of advection that is very small compared to other sources and sinks of energy. However, the outgassing of vapour might be exaggerated by not accounting for the possibility of vapour diffusion towards the interior. Such inaccuracies are minimised when the temperature systematically increase with depth (as when radiogenic heating takes place), and when there is a background outward flow of vapour from a sublimation front at larger depths. If the user chooses not to activate the venting option, gas diffusion (including inward flow) is calculated according to Eq. (2), and energy transport due to advection is accounted for in Eq. (1).
NIMBUS stores snapshots of all the physical data it produces, in different ways depending on the type of body rotation: 1) for fast rotators, data is written to file at a fixed temporal resolution; 2) for rotationally resolved bodies, data is stored with a certain rotation–angle resolution for selected full rotational periods. Such sequences can be used to reconstruct the total production rate of the full body, including all latitudinal, solar zenith angle, and day/night effects. Simulations can be stopped and re–started seamlessly on demand, because all information required for the continuation of a simulation are stored. Special programs have been developed to change, e. g., the spatial resolution of a model and interpolate all physical quantities of a previous simulation onto a new grid. The files produced by that program can be read by NIMBUS, which allows for continued simulations under new conditions when necessary. There are also dedicated programs for plotting and analysing stored solutions. NIMBUS is implemented in MATLAB® and also runs on GNU Octave.
2.14 Code verification
The confirmation of the correctness of a computerised physics model is made by validation (comparison with measured data), and/or verification (comparison with other independently implemented computer codes that solve the same problem). Validation requires laboratory measurements of temperatures, gas densities, abundances et cetera as functions of depth and time, for ice/dust mixtures that are subjected to illumination cycles. As of now, no proper attempt to validate NIMBUS has been made.
NIMBUS has very simple basic numerical principles. Cells have a certain mass for each species and phase, and a particular internal energy. If a given mass is being removed from one cell, the exact same mass will be added to another phase of the same cell, added to a neighbouring cell, or ejected to space, depending on the considered process. Energy transfer processes is a similar zero–sum game, with the added component that new energy is injected to the cell network by solar energy absorption and radiogenic heating. Because of this, mass and energy conservation are guaranteed to within machine accuracy by design, as long as there are no outright implementation errors. However, because machine accuracy is finite, small errors will accumulate with time. To exemplify the degree of mass and energy conservation, NIMBUS considered a body in conditions where sublimation was active but segregation, sublimation, crystallisation, and were dormant. The model was run for in the simulated world, during which of the was lost to space. The simulation consisted of time steps, i. e., cycles of mass and energy swapping between cells. Comparing the difference in abundance of the body at the beginning and end, with the documented loss to space, the quantities differed by . The difference in internal energy at the end compared to the beginning (including absorbed solar energy and lost thermal radiation), compared with the documented removal of energy by escaping vapour, differed by . It is evident that coding errors are absent, and error accumulation is slow.
However, mass and energy conservation by themselves do not guarantee correctness. It is also important to verify that the transferred amounts have the correct magnitude. The basic correctness of the heat conductivity and radiogenic heat generation treatment was made by comparing NIMBUS simulations with an analytical solution (e.g. Hevey & Sanders, 2006) to the energy conservation equation for a spherical homogeneous body that is heated radiogenically, cooled radiatively at its surface, and that has temperature–independent heat conductivity and heat capacity. In one particular case, using spatial and temporal resolutions typical of simulations in this paper, the –driven heating of a body from at to a temperature peak at the centre of , followed by a drop to at , was carried out with maximum and mean errors in NIMBUS with respect to the analytical solution of and , respectively. At a depth of below the surface, where temperature gradients are steep and change relatively fast with time, a peak of and cooling to at was obtained with maximum and mean errors of and . Although errors can be reduced further by considering higher spatial resolution and shorter time–steps, this level of uncertainty would be acceptable in most scientific applications.
|
|
Huebner et al. (1999); Huebner et al. (2006), in the context of an International Space Science Institute (ISSI) Team in Bern, Switzerland, made detailed comparisons between five advanced thermophysical codes for a number of standard model cases. Their model ISSI 4a was selected for verification of NIMBUS. That model considers a comet on a typical JFC orbit, that consists of a porous mixture of dust, water ice, and CO ice, assuming sublimation front withdrawal into the nucleus as the result of not considering dust mantle erosion (for a detailed list of parameters, see Huebner et al., 2006). NIMBUS was set up to be as similar as possible to ISSI 4a (which included some functions in Secs. 2.3–2.8 being replaced by those used by Huebner et al., 2006). Despite being set up as similarly as possible, the five codes applied by Huebner et al. (2006) did not produce identical solutions. Figure 1 (left) shows the upper and lower enveloping boundaries of the equatorial daily peak surface temperatures for the five codes of Huebner et al. (2006) as red curves. It is not necessarily the same code that is responsible for the lowest or highest temperature throughout the orbit, but the red curves show the full range of the assemble of model solutions. The blue curve in the left panel of Fig. 1 shows the result of running NIMBUS, and as can be seen, it is consistent with the results of the other five models. The NIMBUS pre–perihelion temperatures most closely resemble the Huebner et al. (2006) algorithms A (Enzian et al., 1997) and C2 (Orosei et al., 1999), while the post–perihelion temperatures are closer to algorithm D (Prialnik, 1992).
The upper right panel of Fig. 1 shows the equatorial daily peak production rate, with the two red curves enveloping the set of Huebner et al. (2006) models. For unknown reasons, Huebner et al. (2006) chose not to plot the first of the orbit. Some of the Huebner et al. (2006) models display local fluctuations due to numerical instabilities that sometimes become as large as the difference between the five models. Those fluctuations were not reproduced in the upper right panel of Fig. 1, but the red curves trace the upper and lower production rates for the numerically stable curves in regions where some models had problems. The NIMBUS solution (blue solid curve) closely follows the upper envelope curve, and is most similar to algorithm C2. The lower right panel of Fig. 1 shows the equatorial daily peak production rate. Huebner et al. (2006) only showed the results of algorithms B (Benkhoff & Boice, 1996) and C2, being the upper and lower red curves in the lower right panel of Fig. 1, respectively. The NIMBUS pre–perihelion solution closely follows algorithm B, whereas the post–perihelion solution is intermediate between those of B and C2.
Obtaining identical numerical solutions to the same (highly non–linear) physical problem is not trivial because of differences in numerical schemes and certain variability in the treatment of physical processes, as explained in detail by Huebner et al. (2006). It is clear from Fig. 1 that NIMBUS performs as well (or as poorly) as comparable advanced thermophysical modelling codes.
3 Model parameters
All model bodies studied here are assumed to have formed without short–lived radionuclides. That is done in order to isolate the thermal evolution caused by the unavoidable heating by the early Sun and by long–lived radionuclides. Therefore, the presented cases should be considered lower limits on the degree of processing (for the particular combinations of rock abundances and porosities considered here). Whether the real Primordial Disk objects contained short–lived radionuclides at the time of their formation is currently unknown. This uncertainty is because we do not yet know if the Solar Nebula contained short–lived radionuclids at birth, or if they were injected locally during formation, e. g., as suggested in the supernova–triggered model by Boss & Vanhala (2000) and Boss et al. (2008). Until returned samples from the outer Solar System has settled this issue, it is motivated to consider models both with and without short–lived radionuclides: this paper focuses on the latter.
Three different body sizes have been considered, in order to better understand the effect of thermophysical evolution in the Primordial Disk for different types of bodies. The three body classes are named after similarly–sized objects in the current Solar System for easy reference. Although the simulations describe how objects of these sizes may have evolved in the early Solar System, the identifications with real objects should not be taken literally. Because of differences in formation distance, formation time, and subsequent evolution, the model bodies are not necessarily good representation of the current properties of their namesakes. The smallest body under study has a diameter of and is referred to as the “67P/Churyumov–Gerasimenko” case (or the 67P case for short). It is representative for the majority of JFCs, dynamically new comets, and small members of the Centaur and Kuiper Belt populations. Such bodies are not heated significantly by long–lived radionuclides, therefore, those models target the solar–driven evolution most directly. The intermediate body is an analogue to Comet C/1995 O1 (Hale–Bopp) with (Szabó et al., 2012). The largest model body has which is similar in size to one of the largest Centaurs, (2060) Chiron (Fornasier et al., 2013), as well as to the saturnian satellite Phoebe (Castillo-Rogez et al., 2012) that likely was captured (Johnson & Lunine, 2005) from the Primordial Disk.
The Hale–Bopp and Chiron/Phoebe cases consider bodies sufficiently large to make long–lived radionuclides important heat sources, although the exact amount of heating depends on the abundance of refractories (radionuclides are only present in the rocky material). The refractories–to– mass ratio is not well known for the bodies that once populated the Primordial Disk. It is easiest to measure for bodies so large that they cannot possibly have substantial macro porosity. Pluto and Charon, with well–determined sizes and masses thanks to the New Horizons flyby, have and (McKinnon et al., 2017), respectively. However, those may be upper limits to the original composition, considering that substantial amounts of volatiles may have been lost in the collision that formed the Pluto–Charon system. Several attempts to constrain for Comet 67P/Churyumov–Gerasimenko have been made. The total mass loss of the nucleus during the Rosetta mission is well known, (Pätzold et al., 2019), and the total loss of water vapour and other gases can be estimated from in situ measurements. This yields the mass ratio of dust to water vapour that is not necessarily identical to of the nucleus. Ice within large coma boulders may escape detection, which would lower with respect to . Dust that release water vapour before returning to the nucleus as airfall (Thomas et al., 2015a, b; Keller et al., 2015, 2017; Davidsson et al., 2021b) could increase above . Biver et al. (2019) obtain from gas measurements by MIRO, while Combi et al. (2020) obtain based on ROSINA data. Several estimates made with different methods and based on data from other Rosetta instruments are available (see review by Choukroun et al. (2020)). Based on COSIMA data, Choukroun et al. (2020) obtain the mass ratio of refractories to all ices, which translates to if using and from Combi et al. (2020), or if using the values from Biver et al. (2019). The lower limit is the ice–richest model of 67P published to date. The ice–poorest estimate is based on GIADA data (Fulle et al., 2017). Davidsson et al. (2021a) used NIMBUS to reproduce the water production curve of 67P and found that was necessary to make that possible. The albedo analysis of the interior of a boulder that was broken open by the Philae lander indicated (O’Rourke et al., 2020). As a compromise, is applied in the current work. It is assumed that all water ice initially is amorphous.
Gerakines et al. (1999) observed a diverse set of interstellar environments (massive protostars, sources near the Galactic Centre, and Taurus dark clouds) with the Infrared Space Observatory and found that the molar abundance of ice relative water had a narrow range of –. The average concentration was . Accordingly, it is here assumed that the amount of is relative by number. Gerakines et al. (1999) also find that , hence it is assumed that is relative by number. In terms of the total body mass, the composition of all model bodies is refractories, water, carbon dioxide, and carbon monoxide (i. e., the mass ratios of and relative water are 0.41 and 0.24, respectively).
The partitioning of and between entrapment in amorphous water ice and other storage modes aimed at making the crystallisation process energy–neutral. That is to say, the crystallisation energy should be entirely consumed when liberating the occluded and . This (admittedly arbitrary) assumption was made in order to minimise the evolutionary changes triggered by the energy release during crystallisation, in order to isolate the effects of solar and radiogenic heating. Furthermore, it was required that should absorb of energy during crystallisation, and that the abundance of in condensed should be by number. This was made by considering the latent heat of at and of at (the highest and lowest temperatures, respectively, for which the Huebner et al. (2006) formulae formally are valid), which is the temperature interval where the bulk of crystallisation takes place (on the long time–scales relevant for the considered problem). With those constraints, the partitions as follows: is condensed, is trapped in , and is trapped in amorphous . The is predominantly condensed () while the remainder () is trapped in amorphous . The large contribution of during absorption of energy released during crystallisation was selected in order to place roughly half the in the condensed phase. If a 50/50 division of energy absorption between and is considered, only of the is been placed in the condensed phase, potentially biasing the loss time–scales towards too low values.
The hydrostatic code that solves Eq. (9) requires that compacted densities and mass fractions are specified for ices and refractories. Ices have the compact density
| (72) |
where the mass fraction of water in volatiles is according to the composition defined above, (Weast, 1974), and the densities of both and are taken as (Satorre et al., 2009). For a chondritic composition, the abundances of the 12 most common elements, according to Lodders (2003), are applied. Specifically: all Fe, S, and Ni are used to form metal/sulphides; all Mg, Si, Al, Ca, Na, and three O per Si are used to form rock; all C, N, and three H per C are used to form organics. The combination of metal, sulphides, and rock are here called “minerals”, and the combination of minerals and organics is called “refractories”. If so, there is a ratio of minerals:organics, and the minerals have metal/sulphide by mass. The density of metal/sulphide is taken as , which is the average of the densities for iron and troilite (Tesfaye Firdu & Taskinen, 2010). The density of rock is taken as that of forsterite, (Horai, 1971). If so, the density of minerals is . Taking the density of organics as that of “spark tholins”, (Hörst & Tolbert, 2013), the density of refractories is .
However, the composition of bodies in the outer Solar System is not necessarily chondritic. The 67P dust analysed by COSIMA is iron–poor and substantially more carbon–rich than chondritic material, if normalising both to the silicon abundance (Bardyn et al., 2017). With metal/sulphide down to the mineral density is . With a minerals:organics ratio drastically changed to , the density of refractories drops to . The hydrostatic code was evaluated for refractories by mass, using , and ices with . The initial total mass was set such that a body should have a bulk density , as for 67P (Preusker et al., 2015). Under these conditions, and in the absence of an external dynamic pressure (Eq. 11), the self–gravity compresses the model body to a diameter , having , and an average porosity of . In order to obtain the desired and it was necessary to apply an accretion velocity of , resulting in . The body is nearly homogeneous. The same composition and –value was applied to the Hale–Bopp and Chiron/Phoebe cases. The porosity profiles of the three model bodies are shown in Fig. 2. Self–gravity effects are rather weak at , lowering the core porosity to 60% and yielding a bulk density and average porosity . Self–gravity effects are stronger at , lowering the core porosity to 34% and yielding a bulk density and average porosity . Meech et al. (1997) found that the observed coma of Chiron, if interpreted as a bound ballistic atmosphere, suggested . The radiometrically determined diameter of the primary in the binary Centaur (65489) Ceto/Phorcys system is based on Herschel Space Observatory observations according to Santos-Sanz et al. (2012), i.e., a similar size to the Chiron/Phoebe case. Using the mass derived by Grundy et al. (2007), the bulk density would be (Santos-Sanz et al., 2012). However, Grundy et al. (2007) obtained a smaller Ceto diameter of based on Spitzer Space Telescope observations, and consequently, a higher bulk density of . It is reassuring that the model of the Chiron/Phoebe case is consistent with the estimate of Meech et al. (1997), and is intermediate between the estimates of Grundy et al. (2007) and Santos-Sanz et al. (2012) for a similarly–sized body. However, it should be remembered that Phoebe itself has a density (Matson et al., 2009), which likely is caused by heating and compaction (Castillo-Rogez et al., 2012).
Comets are often modelled with albedo and emissivity and those values are applied here despite the possibility that near–surface ice may have made the Primordial Disk objects substantially brighter. The solar energy absorption in the current models could therefore be somewhat exaggerated, although solar wind sputtering and cosmic ray bombardment may have caused non–thermal removal of ice in the top few millimetres and a very low albedo. Constituent grains in comet material are typically of order in size (e.g. Brownlee et al., 2006; Bentley et al., 2016), hence is applied. Such grains are assumed to cluster into –sized pebbles (Blum et al., 2017), motivating the choice to use and (i. e., diffusivity is calculated based on the assumptions that the cylindrical tube radii are similar to the diameters of pebbles, and that channels may be reasonably straight over distances of tens times their radii, motivating the usage of and ). The orbits of all model bodies are assumed to be located in the midst of the Primordial Disk () and they are assumed to be circular and placed in the ecliptic plane. Fixed spin axes of were applied, using an intermediate obliquity to prevent the polar regions to be permanent cold–spots (spin axes are not likely stable on the time–scales considered here, but this may be considered a reasonable average orientation).
The current simulations always considered a detailed diffusion treatment for and . The same treatment was applied to as long as condensed ice existed within the model body. After condensed ice had been removed, vapour produced by segregation and crystallisation was vented directly to space.
4 Results
4.1 The “67P/Churyumov–Gerasimenko” case
Three versions of the 67P case were considered, having disk clearing times of , , and . The orbital skin depth (i. e., the depth over which the amplitude of temperature oscillations, caused by obliquity and orbital motion, are damped by an –folding scale) was for the warmest body () and for the coldest body (). In order to resolve the orbital skin depth and properly evaluate near–surface heat fluxes, near–surface radial grid cell thicknesses should be smaller than the skin depth by a factor of a few. In these simulations, the thickness of the uppermost grid cell was , and cells grew to at the core by geometric progression, resulting in 102 radial cells. With 18 latitudinal slabs, the 67P case model body was spatially resolved by a total of 1,836 cells.
The case exposed the model nucleus to the peak illumination of the young protosun ( was avoided to allow for a gradual thermal adjustment to the harsh conditions). It was modelled for . Under these conditions, half the condensed has been removed by , is remaining at , and all of it has been lost by . This illustrates the gradual slow–down of the (close to exponential) loss rate with time: losing the last 10% took almost as long as losing the first 90%. The ice is not a safe haven for the , because of the has segregated and left the model body by . Extremely small deposits ( of the original amount) survive near the poles, in a layer 80– below the surface, but even those layers have lost of their original CO abundance. The sublimation of is rather substantial near the surface. The front withdraws below the surface in the first and then stops there. The degree of crystallisation is substantial near the surface, but of the global amorphous water ice reservoir remains at the end of the simulation. The degree of crystallisation is at least 10% in the top , and more than in the top (all water ice at the surface is cubic). Therefore, if comets formed right after the accretion of the envelope and the massive disk, the only plausible storage medium of hypervolatiles like would be amorphous water ice. Model results are summarised in Table 13.
|
|
|
|
Figure 3 shows temperature, concentrations of condensed and –hosted , as well as the gas pressure for the 67P case at . The temperature in the top few hundred meters is , which is too warm for condensed or –hosted but sufficiently cool to preserve most and amorphous water ice. The sublimation front has withdrawn under the surface. Note that a substantial fraction of the core has a concentration that has been elevated by above the initial value because of downward diffusion and recondensation of vapour. The region where has segregated out from is significantly thicker at the equator than at the poles ( versus ). The pressure map reflects the location of the two fronts: a broader and lower pressure peak around the sublimation front, and a narrower and higher pressure peak around the segregation front.
The second 67P case was started with a fully opaque disk at that reached full transparency at . It was propagated to . Because of the lower peak protosolar luminosity ( compared to for the case) the loss rate is lower. Nevertheless, after onset the model body has lost all condensed . The peak temperature is sufficiently low to prevent significant segregation: 98% of the trapped in the remained at the end of the simulation. The mixture is still present at the surface, though the concentration has diminished by . The losses of and amorphous water ice are insignificant, with both species being present in the top cell.
The third 67P case started at , with the disk fully cleared up at , and it was first propagated until . This is a time period during which the protosolar luminosity is expected to have been very similar to the current one. The model object lost all condensed in , while segregation out from , loss of through sublimation, and loss of amorphous water ice through crystallisation all were insignificant.
All 67P cases were run without activating the long–lived radionuclides. This can be motivated by considering Eq. (11) of Prialnik & Podolak (1999) with parameters relevant for the current investigation, showing that the cooling rate exceeds the heating rate of long–lived radionuclides at diameters below .
4.2 The “Hale–Bopp” case
The Hale–Bopp case was started at and the disk cleared at . The body was spatially resolved by 109 radial cells (growing from thickness at the surface to at the core), and 18 latitudinal slabs. The model was propagated for . Radiogenic heating by long–lived radionuclides is an important heat source for a body with . The condensed was lost in . The bound to is reduced to half on a time–scale and the amount was down to at the end of the simulation. Crystallisation was substantial as well, with half the amorphous water ice and the trapped and gone after . At the end of the simulation, amorphous ice remained. The amount of condensed ice is 126% of the amount present at the beginning of the simulation. The reason for that enhancement is that much of the that is driven out of amorphous water ice during crystallisation, unlike , does not diffuse into space but recondenses within the nucleus. The , , and amorphous water abundances are final because a steady–state has been reach were the body is no longer evolving.
The upper panel of Fig. 4 shows the core temperature as a function of time. The curve is best understood if compared with the core abundance, shown in the middle panel of Fig. 4. The core is heating from its initial temperature during the first . At that point, sublimation of at the core is initiated, establishing a steady temperature just below for the next . During the following the core temperature falls somewhat and reaches a local minimum of at because of intensified sublimation. At this point all condensed is gone, and the temperature starts climbing rapidly to reached at . At this point, segregation out of is initiated at the core. Because of the cooling effect the temperature stabilises anew for until the core segregation is completed at . The core temperature rises to at . At this point, crystallisation is initiated, and the middle panel of Fig. 4 shows that the abundance of amorphous water starts to decline as the abundance of cubic water is rising. The abundances of and had been set so that the energy needed to liberate these species (taken as their latent heats of sublimation) should equal and neutralise the energy release of crystallisation (see Sec. 3). However, the saturation pressure of is extremely low at these temperatures, and much lower than the pressure of the released during crystallisation. As a consequence, the condenses and therefore releases its latent energy. This causes rapid heating and the temperature reaches at . There is some further radiogenic net heating and the core temperature peaks at at . However, at this point a balance is reached between the radiogenic energy production rate and the thermal radiation cooling rate at the surface of the body. A steady–state temperature is reached at around . There is no net sublimation, and it is also too cold for the cubic–to–hexagonal transition.
The lower panel of Fig. 4 provides further insight into the mass loss processes. It shows the depths of the condensed sublimation front, the segregation front, and the crystallisation front, as functions of time. Condensed is being removed through combined protosolar and radiogenic heating, but the former process is dominating. Although the core abundance reduction of condensed is initiated at , the complete removal is a wave that spreads from the surface and downwards. When the core loses the last of its condensed , the rest of the body is already free from pure ice. The segregation process, on the other hand, is a wave that spreads from the core and upwards. As previously mentioned, the core has lost all its at . The region of –free then expands outwards over the following until the wave grinds to a halt about under the surface. Likewise, crystallisation is a process that spreads from the core and outwards, albeit at a significantly faster rate – the entire process takes and stops below the surface. Note that there are no movements of the segregation or crystallisation fronts at , showing that the body has reached a thermal steady–state. A warm interior and a substantially cooler thin near–surface region (explaining the survival of amorphous water ice and –laden near the surface) is typical of radiogenically heated and radiatively cooled bodies (e.g. Hevey & Sanders, 2006).
Figure 5 shows the abundances of trapped in , ice, amorphous water ice, and cubic water ice, as functions of depth at the end of the simulation. The three panels show different levels of zoom on the surface region. Below a depth of the entire volume has crystallised and turned to cubic water ice. As already has been mentioned, the occluded driven out during this process has recondensed. However, Fig. 5 shows that this recondensation has not necessarily been immediate and local. Instead, there has been a substantial upward diffusion of vapour before it recondensed, forming a mass enhancement – below the surface. There has even been a small reduction of ice within a –radius core, showing that some net sublimation of has taken place as well. That vapour has diffused upwards and contributed to the enhancement at depth –. Although the temperature at – depth has not been sufficiently high to cause widespread crystallisation, it has been high enough to drive out from the . The outermost of the body shows the effects of the early short–duration protosolar heating. Sublimation of ice has emptied the upper . Note that a fraction of the liberated vapour has diffused downwards and created a somewhat elevated ice concentration down to a depth of . Additionally, the entire region down to has experienced various levels of crystallisation. At depth, the amount of amorphous water ice is about half of the original abundance, and it falls to zero at the surface. The only ice that would be visible at the surface of this Hale–Bopp analogue would be cubic water ice, at this stage of its existence. The protosolar heating has also driven out all from its host down to a depth of . Because of the protosolar heating from above, and the radioactive heating from below, –bearing ice is only present in a comparably thin layer located – below the surface.
Note, that a later clearing time than would have reduced and perhaps eliminated the crystallisation, segregation, and sublimation in the near–surface region. Also, if the internal crystallisation, segregation, and re–distribution would become smaller and perhaps negligible. However, in this context it is worth remembering that C/1995 O1 (Hale–Bopp) is considered a dusty object, with according to Jewitt & Matthews (1999) and according to Dello Russo et al. (2000). Therefore, although the current simulation is not necessarily representative of all –class bodies, it may describe the deep interior of Comet Hale–Bopp itself fairly well.
4.3 The “Chiron/Phoebe” case
Three versions of the Chiron/Phoebe case were considered, each having 91 spatial cells (growing from thickness at the surface to at the core) and 10 latitudinal slabs. The first considered the same time–dependent protosolar luminosity as the previously discussed models, starting at , having a disk clearing at , and included long–lived radionuclides. A second model also included radiogenic heating but applied a fixed current–sun luminosity to better understand the sensitivity of the results to the heating sequence. A third model also considered current illumination conditions but removed the radiogenic heating. Such a model would be more representative of bodies having a substantially lower fraction of refractories than the other models considered here.
The radiogenically heated body exposed to protosolar heating was simulated during . Compared to the Hale–Bopp case, the Chiron/Phoebe body has a times larger volume and a 2.7 times lower surface–to–volume ratio (i. e., it is less capable of dissipating the radiogenic heat it produces, and is therefore heated to a higher temperature). Because of its lower porosity, the Chiron/Phoebe case body has a mass that is almost 30 times higher than the Hale–Bopp case body. Despite having 30 times more to get rid of, the time–scale for doing that is just , marginally longer than the time–scale for the Hale–Bopp case. At the end of the simulation only of the –trapped remains, and of the amorphous ice. The condensed is 137% of the initial abundance because vapour released during crystallisation eventually recondenses with very small losses to space.
|
|
The left panel of Fig. 6 shows the core temperature versus time, and the final distribution of ices as function of depth at the end of the simulation. The temperature curve is reminiscent of the Hale–Bopp case (Fig. 4), but the values are higher. The plateau caused by ice sublimation is (compared to for the Hale–Bopp case). Because the vapour needs to diffuse over a longer distance, and because the lower core porosity reduces the diffusivity, the Chiron/Phoebe body evolve into having higher pressures and sublimation temperature than the Hale–Bopp case. The late steady–state temperature reaches , compared to for the Hale–Bopp case because of its smaller surface–to–volume ratio. Due to relatively high internal temperatures, all amorphous ice below a depth of has crystallised (it is too cold for the cubic–to–hexagonal transformation). The vapour release during crystallisation has diffused upwards and recondensed in a zone located 3– below the surface, where the ice abundance has increased by up to a factor 3.4 above the initial abundance. The radiogenic heating has caused segregation of all out of ice below a depth of .
Because of the intense protosolar heating in the first few Myr of the simulation the abundance of amorphous water ice falls from nearly the original value at a depth of to zero at the surface. The abundance of cubic water ice increases from zero to unity over the same near–surface region. The ice in the uppermost has sublimated. Furthermore, the has been segregated out of in the top . That means that the only remaining –bearing ice is found in a thick zone located roughly under the surface.
The second Chiron/Phoebe case, shown in the right panel of Fig. 6, was propagated for . By employing a fixed current–day solar luminosity, this model avoids the early scorching and is roughly equivalent of using . The time–scale for losing all condensed (), the amounts of amorphous water ice (), condensed (), and trapped in (1.9%) are very similar to that of the previously discussed Chiron/Phoebe case. That is because radiogenic heating is the major mechanism that drives the evolution. The only significant difference between the models are found at the surface: the entire upper layer is rich in amorphous water ice, the entire upper layer is additionally rich in –bearing , which means that amorphous water ice and ice are present at the very surface, both trapping .
Finally, a Chiron/Phoebe case without heating by long–lived radionuclides was considered, assuming a current–day solar luminosity. As expected, the loss of condensed is comparably slow in this case. By about of the remained, falling to at . At the end of the simulation () about of the condensed still remained (the time–scale of complete loss is estimated to be ). The model body only lost of the stored in , it lost of its , and the amount of crystallisation was even smaller.
| Model case | Diameter | Cond. | in | Cond. | Amorp. | ||
|---|---|---|---|---|---|---|---|
| 67P, | 0% | 96.9% | 94.6% | ||||
| 67P, | 98.3% | 99.9% | 100% | ||||
| 67P, | 99.9% | 100% | 100% | ||||
| Hale–Bopp | 13.2% | 126% | 32.5% | ||||
| Chiron/Phoebe | 1.7% | 136% | 7.2% | ||||
| Protosun, LLR | |||||||
| Chiron/Phoebe | 1.9% | 137% | 7.1% | ||||
| Current Sun, LLR | |||||||
| Chiron/Phoebe | 100% | 100% | 100% | ||||
| Current Sun, no LLR |
5 Discussion
The simulations presented in Sec. 4 rely on NIMBUS and not NIMBUSD, i. e., erosion of near–surface solids due to the gas flow is not accounted for. This can be motivated as follows. When the reduction from full to negligible opacity of the Solar Nebula takes a sufficient amount of time, the CO front withdraws underground so slowly that the outward force associated with the gas pressure is too weak to break up and eject the surface material. Once the surface is fully exposed to the protosolar early luminosity, the CO is already so far underground that pressure gradients are negligible compared to material cohesion. Regarding , its volatility is so low, that it is unable to cause significant erosion even when being exposed to at when present at the surface. To illustrate these points, a version of the was run with much finer spatial resolution ( surface cells) to evaluate outward forces more accurately. The model was run for (about half an orbit) from (when the effective luminosity was , corresponding to illumination at in the current Solar System). After , the CO front had withdrawn to – depending on latitude, and the strongest force recorded was merely . At that point, the model object was subjected to the full , causing a force spike of that rapidly decreased as the CO drew even deeper. The force is a generous upper limit, as the CO in reality would have withdrawn even deeper before Solar Nebula clearing at . The force due to sublimation from within the top cell peaked at , and declined rapidly as well. Considering that dust mantle analogues with a porosity of have a tensile strength of according to laboratory measurements (Güttler et al., 2009), vapour escape without erosion and ejection of solids is a reasonable assumption.
Furthermore, in their modelling of Centaur (2060) Chiron, Capria et al. (2000) found that the gas flow only was capable of removing loose dust on the surface, that sustained weak activity at perihelion and waned after a few orbits (as sublimation and crystallisation fronts withdrew underground). With the water ice remaining intact at the low temperatures of Chiron, there was no mechanism in the model to liberate the grains and to replenish the surface deposit of loose dust. These simulations included CO ice and were made for heliocentric distances that overlap and somewhat subceeds the (effectively) simulations in Sec. 4 (Chiron has a perihelion distance of ). One mechanism that possibly could remove water ice in a surface layer and create a weakly cohesive mantle of dust is solar wind sputtering. In laboratory measurements of that process, Muntean et al. (2016) found that a layer of water ice could be removed in objects located at during the Solar System lifetime, constituting an upper limit on the thickness of dust mantles formed accordingly, for objects in the Scattered Disk that later become Centaurs. A global layer of dust (with an assumed density of ) would be able to sustain the observed dust production rate of Chiron (; Luu & Jewitt, 1990) for . Other Centaurs have dust production rates of up to order – (Jewitt, 2009), that could be sustained for at most – (considering that most Centaurs are smaller than Chiron). A working hypothesis is therefore that Centaurs (the best analogues of Primordial Disk objects) do not experience large–scale erosion, but merely removal of minor unreplenishable deposits of surface dust. This further motivates the assumption in Sec. 4 that erosion can be ignored.
Based on the NIMBUS simulations it is possible to address the questions formulated in Sec. 1. The questions are here re–stated for convenience.
-
Q#1. Under what conditions is loss of condensed ice complete or partial in the Primordial Disk? If loss of condensed ice is complete, what are the time–scales for that loss in bodies of different size in the Primordial Disk?
According to the simulations in Sec. 4, the loss of condensed CO proceeds at an exponentially declining rate. If the bodies remain at sufficiently long, the process runs to completion. The current simulations therefore do not support the hypothesis of De Sanctis et al. (2001), who suggested that a core of condensed CO might be preserved if the body reaches a state where sublimation is balanced by condensation. Such a quasi–stationary condition does not arise, at least not for the combination of model parameters applied in the current work. One scenario that would allow for a long–lived reservoir of ice is the existence of an internal compacted shell or core with sufficiently low porosity to prevent diffusion from that region (at least until the temperature and gas pressure become high enough to forcefully create channels through which vapour can diffuse). However, considering that self–gravity is not capable of removing porosity even in –class bodies (see Fig. 2) it is hard to imagine that such compacted interiors would be common, particularly for JFC–sized objects. The time–scale for complete loss of ice is about 70– for a body, depending on the clearing time of the Solar Nebula (Table 13). The exact number depends on the actual CO abundance, conductivity, and diffusivity of the body, but a loss time–scale of order is realistic for JFC–sized bodies. Bodies with sizes ranging between those of Hale–Bopp and Chiron have CO ice loss time–scales of roughly 10– when being relatively dust–rich (). This time–scale is comparable to the shortest proposed lifetime of the Primordial Disk (e. g., according to Nesvorný & Morbidellli, 2012). If the Primordial Disk indeed was dispersed that early, it is likely that bodies in the – range could preserve a fraction of their original CO ice content, at least if they were less dusty than considered in the current simulations. Because the loss is partially driven by radiogenic heating, the long–term survival of CO ice in such bodies depends largely on what heliocentric distance they are scattered to. Bodies ending up in the Oort Cloud would obtain sufficiently low near–surface temperatures for CO vapour diffusing from the deep interior to condense near the surface (in a fashion akin to the near–surface concentrations seen in Figs. 5 and 6). Bodies ending up in the Kuiper Belt would likely lose all condensed CO, considering the work by De Sanctis et al. (2001) and Choi et al. (2002) in combination with the results in Sec. 4. Sufficiently large and dust–poor bodies ending up in the Scattered Disk may have been able to conserve some CO ice if their semi–major axes are sufficiently large. However, the apparent lack of activity beyond (Bockelée-Morvan et al., 2001; Jewitt et al., 2008; Jewitt, 2009; Li et al., 2020) suggests that such bodies are not common among Centaurs (and by inference, the Scattered Disk) or Kuiper Belt objects. The most likely scenario is therefore that bodies indeed lost their CO ice, which places constraints on the Primordial Disk lifetime: it could not have had a lifetime shorter than 10– (unless the initial CO abundance was significantly lower than considered here). If the bodies contained less radionuclides than assumed in this work, the Primordial Disk lifetime increases accordingly. If the actual dust–to–water–ice mass ratio of Centaurs could be determined accurately (being the former Primordial Disk large members that can be reached by spacecraft with less difficulty), thermophysical simulations of the type considered in this paper could be used to determine the shortest allowable Primordial Disk lifetime. If the Primordial Disk disruption was associated with the Late Heavy Bombardment, all bodies with would have lost all their CO ice prior to disruption (by , or after at most half of the Primordial Disk lifetime), even in the unrealistic scenario that they are dust–free and experienced no radiogenic heating at all.
Observations of debris disks around young stars sometimes reveal emission of CO gas from within those structures. Confirmed cases in systems of known age include NO Lup (; Lovell et al., 2020), TWA 7 (; Matrá et al., 2019), HD 138813 (; Lieman-Sifry et al., 2016; Hales et al., 2019), HD 129590 (10–; Kral et al., 2020), HD 131835 (; Moór et al., 2015; Lieman-Sifry et al., 2016; Hales et al., 2019), HD 181327 (; Marino et al., 2016), Pictoris (; Matrá et al., 2017a), and Fomalhaut (; Matrá et al., 2017b). This CO is not considered primordial by these authors (in the sense of having remained in gas form since the time of formation), but originates from exocomets that store the CO as ice. These CO disks often extend up to 100– from the central star, beyond the CO snow line for most systems. These authors therefore suggest that CO ice from the interior of exocomets is being exposed during collisional cascades in these debris disks, and is vaporised through UV desorption. That may certainly be the case, but the simulations in Sec. 4 suggest that an additional CO gas source could be radiogenically heated planetesimals, perhaps in combination with cratering that penetrates surface regions where CO might recondense, thus allowing for venting of CO gas from deeper and warmer regions. Such a source of is potentially preferable in systems where the gas disk extends beyond the debris disk (i.e., presence of CO vapour in regions without detectable evidence of high collisional activity), and could explain why gas–rich disks are not preferentially associated with dust–rich debris disks (Lieman-Sifry et al., 2016). Bodies in the 70– class with can sustain CO release for 10–, but that time–scale can easily be extended if the bodies are less dusty (i. e., the CO release is slowed down by a smaller radiogenic power) or simply have larger size.
It is clear that Primordial Disk bodies with sizes of a few kilometres stand no chance of still containing primordial CO ice at the time the disk was dispersed and some of them were placed in the Oort Cloud. The notion that dynamically new comets are pristine in the sense that they still carry all species they formed with and have not experienced any form of global–scale processing, is therefore not necessarily correct. If a dynamically new comet contains ice there are at least two possible reasons: 1) the comet experienced heating strong enough to cause segregation and/or crystallisation just prior to its ejection, and the released condensed near the surface when the object reached sufficiently large heliocentric distances; 2) the comet did not originate from the Primordial Disk. The first scenario requires that the object spent some time substantially closer to the Sun than prior to ejection, because both mixtures and amorphous ice are stable at that distance. The viability of that scenario must be tested in future simulations, because it is not clear whether the CO has time to escape before sufficiently large distances are reached, or if the accumulated amount of condensed CO is sufficiently large to drive observed levels of activity. Furthermore, if this was a common and efficient process, one might have expected a larger number of active comets at large distances. One version of the second scenario concerns planetesimals that formed within the giant planet zone and were moved to the Oort Cloud during the first after the Solar Nebula cleared up (even shorter time–scales may be necessary due to the higher level of protosolar illumination compared to ). The work of Brasser et al. (2007) shows that an object with a perihelion near Jupiter needs to reach to avoid orbit circularisation by Solar Nebula gas, and it needs – (depending on inclination) to experience changes to during close encounters with Jupiter that are of a size similar to itself (i. e., enabling a rapid increase of that allows stars in the birth cluster and galactic tides to raise the perihelion out from the planetary region and incorporate it into the Oort Cloud). In order to successfully create an Oort Cloud orbit, this process must necessarily take place before the Solar Nebula gas drag has time to move the perihelion away from Jupiter’s orbit. That time–scale is if but approaches if (Brasser et al., 2007). At Saturn, – is needed to avoid circularisation, – is needed to initiate large orbital perturbations, and Oort Cloud insertion must take place in less than –. For Neptune, the insertion must take place within –. The time–scale of Oort Cloud injection by Jupiter is shorter than the CO loss time for bodies, therefore constituting a mechanism that could create a subset of Oort Cloud comets that actually contain CO ice. The time–scales for Saturn, Uranus, and Neptune are longer, but it must be remembered that once has reached values high enough to initiate the Oort Cloud insertion sequence, a very small fraction of the time is spent at distances small enough to allow for efficient loss. It is therefore possible that all giant planets contributed to populating the Oort Cloud with some objects rich in CO ice. However, the fraction of such objects in the Oort Cloud could be small. The number of planetesimals in the giant planet region that avoided accretion may have been low, and of these merely 1– end up in the Oort Cloud according to the simulations of Brasser et al. (2007). Furthermore, objects with never reach the Oort Cloud because the Solar Nebula gas drag is too strong. Comet C/2017 K2 could be a relatively rare example of a body that formed among the giant planets and was inserted into the Oort Cloud long before the planetary gravitational instability took place and the Primordial Disk was disrupted. If that is the case, it also shows that the giant planets had grown sufficiently large to feed the Oort Cloud before the Solar Nebula in that region became transparent, otherwise Comet C/2017 K2 would not have contained ice. A second version of the second scenario is that some, perhaps a majority, of Oort Cloud objects were captured from foreign stellar system in the solar birth cluster (Levison et al., 2010). If that is the case, the current simulations offer little insight into their histories.
As mentioned in Sec. 1, segregation of CO out from could be another potential reason for unusually distant activity, in addition to, or instead of, ice. The second group of questions concern mixtures.
-
Q#2. What degree of release from entrapment is expected to take place in the Primordial Disk? Specifically, at what depth is the segregation front expected to be located?
According to Sec. 4 simulations, small JFC–like objects could have segregated fully if the Solar Nebula cleared very early, while almost complete preservation takes place if . –class objects preserve roughly while –class objects preserve about –trapped , respectively, when (Solar Nebula clearing time matters little at these sizes). If radiogenic heating is insignificant, preservation is essentially complete. Admittedly, exact numbers depend on the preliminary –values assigned in Sec. 2.7. However, the object reaches peak temperatures of and at the core and on the surface, respectively, and the laboratory measurements of Luna et al. (2008) and Satorre et al. (2009) show that no hypervolatile can survive within at such temperatures. The object never exceeds at the core, where entrapment is very stable. With surface temperatures reaching the level of segregation in the surface region could be substantial, and perhaps more so than the current simulations suggests. Yet, it seem clear that small objects exposed to the protosun within the first million years could display a very wide range of abundances, while later exposure likely allows for primordial abundance–levels. Small objects with sufficiently large would have mixtures throughout their volumes, including the surface. The 70––class models studied here would have mixtures confined to the upper – (because of radiogenic heating if as dusty as ), and potentially lack in the upper few hundred meters if exposed to the protosun very early.
The suggestion that Comet 67P contains intimate mixtures (Gasc et al., 2017) excludes an early Primordial Disk exposure to intense protosolar radiation. The need for an extended opaque era may imply a late formation time for comets as well, considering that planetesimal formation and Solar Nebula clearing necessarily are related. In this context, “late” means after CAI (potentially earlier), which still is consistent with early giant planet formation. If cometary mixtures are confirmed, and our knowledge of the segregation process is improved, the lack of ice and presence of –trapped may provide important means of dating the clearing of the Solar Nebula. It should also be remembered, that 67P/C–G may have formed in the outermost regions of the Primordial Disk (considering its unusually high fraction of deuterated water; Altwegg et al., 2015). Comets formed closer to the Sun could have experienced earlier formation and Solar Nebula clearing.
|
|
If intimate mixtures indeed are present in comets, they could play an important role in distant activity of comet nuclei and Centaurs. To illustrate this point, Fig. 7 compares two models of a dynamically new comet: one with clean where is traditionally trapped in amorphous water ice (left), and another where the amorphous is clean and the is trapped in (right). These models were run with NIMBUSD, i.e., dust mantle erosion was included. The latter object behaves qualitatively similar to observed comets. It reaches activation around (compare with the discussion in Sec. 1), where the and productions have reached 1– of the strong activity near . At a somewhat larger distance, e. g., at , the activity is 2–3 orders of magnitude lower than at . At , the production rate is times higher than that of . The combination of the two might drive the dust production that often is observed in the region, although telescopes yet do not have the sensitivity to routinely measure the and directly at such distances, to explore their absolute and relative abundances. Segregation should be capable of driving the production observed in Centaurs, with almost the same intensity throughout the – region. The equatorial segregation front in the model of Fig. 7 (right) is merely below the surface at , and such shallow depths could perhaps be maintained over long periods of time if the erosion rate is sufficiently high. If so, segregation could explain the sunward activity observed in Centaurs 29P/Schwassmann–Wachmann 1 and 174P/Echeclus (Jewitt, 2009; Wierzchos et al., 2017).
Figure 7 (right) also shows that significant and extremely remote activity (e. g., a dust production of for C/2017 K2 at ; Hui et al., 2018) may not be possible to maintain through segregation after all, unless such activity represents the “ emission peak” seen in segregation experiments (Luna et al., 2008; Satorre et al., 2009). At the moment, the best explanations for the activity of C/2017 K2 appear to be an origin in the giant planet region or a foreign stellar system, or alternatively, minor near–surface condensates formed during ejection of an actively segregating/crystallising object originating from the Primordial Disk. However, the latter alternative is less compelling because of the perceived difficulty of reconciling a probably common process with the unusual distantly active comets.
The more classical model in Fig. 7 (left) also becomes active around . However, somewhat surprisingly considering the discussion in Sec. 1, the driver is not released from crystallising amorphous water ice. Instead, sublimating is responsible for the activity. At the production from crystallisation is times weaker than that of segregation, and the production rate is times below that of . The production from crystallisation becomes comparable to production at and the model object needs to get within before exceeds . Guilbert-Lepoutre (2012) found that crystallisation becomes an important source of activity at . Her thermophysical model did not consider sublimation of any volatile, but solved the heat flow problem within a rotating spherical (and, in practice, asteroidal) body at different latitudes, for different spin axis orientations. The thermophysical model contained a module similar to Eq. (26), allowing her to calculate crystallisation rates as function of depth, latitude, and time. In the absence of sublimation cooling, the temperatures become comparably high, which explains why crystallisation reaches full force as distant as from the Sun. However, when ice is included in a thermophysical model, the NIMBUSD simulations in Fig. 7 show that the cooling from is substantial. With , the temperature does not become high enough to allow for substantial crystallisation until the object is to within of the Sun. When discussing the activation barrier for Centaurs, Jewitt (2009) dismissed with an argument similar to that for , e. g., his calculations showed that should be active throughout the planetary region while Centaurs are not. However, the simulations in Fig. 7 show that the outgassing drops at least two orders of magnitude between and , becoming completely insignificant further away. Potential reasons for the discrepancy is that NIMBUSD includes heat conduction, rotation (including night–time cooling), and a finite diffusivity for vapour from sub–surface ice, that Jewitt (2009) did not consider in his work. I therefore propose that the activation distance around , observed in numerous long–period comets and Centaurs, is caused by the onset of sublimation, potentially in association with release during segregation, and not necessarily by crystallisation.
The last question concerned ice and amorphous water ice:
-
Q#3. At what depths are the sublimation front and the crystallisation front of amorphous water ice expected to be located after thermal processing in the Primordial Disk?
Bodies that are exposed to the most intense protosolar heating lose all ice in the upper and the degree of crystallisation is at least in that layer. Complete preservation of amorphous water ice is only expected below a depth of . Although small bodies would be amorphous throughout the rest of their interiors, sufficiently large and dusty objects would only have amorphous ice in the upper because of radiogenic heating. Such bodies would also have substantially enhanced ice abundances at – depth because of upward displacement. Bodies with would have ice and amorphous water ice up to their very surfaces (while the internal crystallisation and displacements would still be present in the sufficiently radioactive ones). Whether the and crystallisation fronts are located at the surface, or at depths of, e. g., , , or would have a tremendous influence on their capacity to produce and when entering the inner Solar System for the first time as new arrivals from the Oort Cloud. The observed coma abundance ratio varies by orders of magnitude among comets passing close to the Sun, and the coma abundance ratio varies by orders of magnitude (A’Hearn et al., 2012). Although such differences routinely are interpreted as reflecting variability in the bulk composition of comet nuclei, another possibility is that the abundance ratio diversity is caused by differences in front depths (the larger the depth of the source of vapour, the smaller the mass flux that reaches the surface). That would imply that at least parts of the Solar Nebula (in what would become the Primordial Disk) cleared during the first , when the protosolar radiation was strong enough to cause near–surface sublimation and crystallisation. If so, cometary coma and abundance ratios could be powerful indicators of the birth distance of comets: the closer to the protosun, the earlier could the clearing have taken place and the stronger was the processing, resulting in significantly lowered present coma abundance ratios because of the withdrawn fronts. If coma abundance ratios truly reflected the nucleus bulk abundance, one might have expected a rather bimodal distribution: extremely low abundances in objects formed within the snow line of the species in question, and the full available abundance in objects formed beyond the snow line. The model of the chemical gradients of the Solar Nebula by Dodson-Robinson et al. (2009) shows ice abundances (e. g., of CO) raising from zero to full freeze–out over just . Because the radial distributions of and essentially are step–functions, it is difficult to comprehend why such narrow transition regions should be so well–represented in the observational record. The model by Dodson-Robinson et al. (2009) suggests that all planetesimals formed beyond should have roughly the same ratio. Similar bulk compositions among comet nuclei, combined with monotonically decreasing and amorphous water front depths with increasing heliocentric distance, could possibly be the best explanation for the wide range and lack of bimodality in the observed and distributions.
Comets that were exposed to solar radiation late and/or at large distance would have comparably shallow and amorphous water front depths (and, as may be the case with 67P, would contain mixtures). If placed in the Scattered Disk, and later evolving dynamically into Centaurs, such objects would readily produce observable amounts of dust, and . Other objects, formed closer to the Sun and/or exposed to its heat earlier, would have deeper fronts and perhaps be observed as inactive Centaurs. A Centaur that lacks observable activity may not necessarily have exhausted its near–surface supervolatiles during its time as a Centaur, but its – and –poor crust could have developed already in the Primordial Disk. If such Centaurs are evolving dynamically into JFCs, one can imagine a time–period during which water–driven activity and erosion gradually brings the and amorphous fronts closer to the surface. The and ratios of such bodies would increase with time. A’Hearn et al. (2012) pointed out that periodic “comets with the fewest perihelion passages are depleted in CO relative to the others, exactly contrary to what one would expect from simple evolutionary models.” If the apparent over–abundance of objects with low ratios among recently injected JFCs is statistically significant, that is certainly inconsistent with the notion that physically ageing comets gradually “run out” of supervolatiles. However, a comparably large fraction of low– objects among dynamically young JFCs would be consistent with the idea that some comets may need to erode significantly before their and amorphous fronts are sufficiently shallow to inject large amounts of and into the coma.
If variability in and amorphous front depths exist in dynamically new comets, Centaurs, and recently injected JFCs, and if that variability indeed was created because of Primordial Disk processing, it has an important corollary. The time period during which the variability needs to be established, during the first at the centre of the disk, is short compared to the Primordial Disk lifetime. If this near–surface stratification is destroyed at a later time during Primordial Disk evolution, it cannot be re–established. This includes a consideration of heating of Oort cloud comets by passing O– and B–stars, as well as by supernovae. Stern & Shull (1988) demonstrated that the former have heated of the Oort cloud comets to , while there is a 50% probability that the latter have heated the comets to (also see Stern, 2003). These temperatures are not sufficient to remove ice or to crystallise amorphous water ice near the surface (according to Table 13, a object is heated to without losing and amorphous water ice). Therefore, if one could demonstrate the existence of such primordial stratification, it means that comet nuclei did not experience significant catastrophic collisional processing in the Primordial Disk. Arguments for the perceived lack of heavy collisional processing in 67P and other JFCs, as well as a comet formation and early evolution scenario in which collisional cascades are avoided, were presented by Davidsson et al. (2016).
6 Conclusions
The thermophysical code NIMBUS (Numerical Icy Minor Body evolUtion Simulator) is described at length for the first time, and it is used to investigate the thermophysical evolution of small icy bodies in the Primordial Disk (here, ) during the earliest epoch of Solar System history. Specifically, porous bodies with diameters have been considered, consisting of amorphous water ice and ice (that both trap in different proportions), separate “condensed” ice, and refractories that contain long–lived radionuclides. Time–dependent protosolar illumination consistent with solar–mass Hayahsi– and Henyey–tracks was applied, but coupled with a simple parameterisation of Solar Nebula opacity that exposed the icy bodies to the protosun at different clearing times, , , or into Solar System history. The main conclusions are the following:
-
Small comets () lose all ice that is not trapped inside a less volatile medium in 70– depending on the Solar Nebula clearing time. Comets that were placed in the Oort Cloud, Kuiper Belt, or Scattered Disk at the time of the Primordial Disk dispersal ( after Solar System formation) cannot possibly contain large deposits of CO ice.
-
Large comets (–) lose all ice that is not trapped inside a less volatile medium in 10– if they are dusty (here, a dust–to–water ice mass ratio of 4). Dust–free –class objects lose all free CO ice in .
-
Small comets () do not preserve any mixtures if they are exposed to the protosun too early. Exposure after CAI or later could imply almost complete survival.
-
Large comets (–) may have a thin near–surface layer of mixtures at – that survives radiogenic heating (and if applicable, early exposure to protosolar radiation).
-
Primordial Disk bodies, of any size, that get exposed to protosolar radiation within the first of the Solar Nebula, lose ice in a surface layer being up to thick, and experience partial crystallisation in the upper that reaches completeness at the very surface.
Based one these conclusions and associated discussions, I propose that:
-
1.
Dynamically new comets are not necessarily pristine – the majority have experienced global loss of free CO ice, and some may have crystallised and lost ice in the top few tens of meters, because of thermal processing in the Primordial Disk prior to their emplacement in the Oort Cloud.
-
2.
The large ranges in and observed in dynamically new comets is not necessarily due to intrinsic bulk abundance variability, but because the previously described diversity of and amorphous front depths.
-
3.
The apparent increase of and in short–period comets with number of perihelion passages is because their and amorphous fronts are gradually brought closer to the comet surfaces by erosion driven by sublimation of crystalline water ice.
-
4.
If propositions i–iii are correct, the Primordial Disk cannot have experienced a collisional cascade, because that would destroy the near–surface stratification that only forms in the first .
-
5.
CO gas observed in the debris disks of young stars could originate from radiogenic heating of sufficiently large planetesimals, in addition to, or instead of, collisional cascades and UV desorption.
-
6.
Comet C/2017 K2 (PANSTARRS) probably did not form in the Primordial Disk, but could be a rare example of a body originating in the giant planet region and ejected very early, or was captured from a foreign system in the solar birth cluster.
-
7.
The activity switch–on observed in long–period comets and Centaurs is due to ice sublimation and/or mixture segregation. Presence of ice prevents crystallisation from being an important source of released hypervolatiles beyond .
-
8.
The lack of activity in large Centaurs beyond implies that they are sufficiently dusty, and/or the Primordial Disk was sufficiently long–lived, to deprive them of all free CO ice. The observed CO is likely released from in the 8– region, and from a combination of and amorphous water ice within .
Acknowledgements
Parts of the research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a
contract with the National Aeronautics and Space Administration. The author is indebted to Dr. Pedro Gutiérrez from Instituto de Astrofísica
de Andalucía, Granada, Spain, for invaluable discussions, suggestions, and support during the writing of this paper. The author
appreciates the detailed, constructive, and important recommendations of the reviewer, Dina Prialnik, that
substantially improved the paper.
COPYRIGHT. © 2021. California Institute of Technology. Government sponsorship acknowledged.
Data Availability
The data underlying this article will be shared on reasonable request to the corresponding author.
References
- A’Hearn et al. (2012) A’Hearn M. F., et al., 2012, Astrophys. J., 758, 29
- Altwegg et al. (2015) Altwegg K., et al., 2015, Science, 347, 1261952
- André & Montmerle (1994) André P., Montmerle T., 1994, Astrophys. J., 420, 837
- Azreg-Aïnou (2005) Azreg-Aïnou M., 2005, Monatshefte für Chemie, 136, 2017
- Bailey & Malhotra (2009) Bailey B. L., Malhotra R., 2009, Icarus, 203, 155
- Bar-Nun et al. (1985) Bar-Nun A., Herman G., Laufer D., Rappaport M. L., 1985, Icarus, 63, 317
- Bar-Nun et al. (1987) Bar-Nun A., Dror J., Kochavi E., Laufer D., 1987, Phys. Rev. B, 35, 2427
- Bardyn et al. (2017) Bardyn A., et al., 2017, Mon. Not. R. Ast. Soc., 469, S712
- Benkhoff & Boice (1996) Benkhoff J., Boice D. C., 1996, Planet. Space Sci., 44, 665
- Benkhoff & Huebner (1995) Benkhoff J., Huebner W. F., 1995, Icarus, 114, 348
- Benkhoff & Spohn (1991) Benkhoff J., Spohn T., 1991, in Kömle N. I., Bauer S. J., Spohn T., eds, , Theoretical modelling of comet simulation experiments. Verlag der Österreichischen Akademie der Wissenschaften, Wien, pp 31–47
- Benkhoff et al. (1995) Benkhoff J., Seidensticker K. J., Seiferlin K., Spohn T., 1995, Planet. Space Sci., 43, 353
- Bentley et al. (2016) Bentley M. S., et al., 2016, Nature, 537, 73
- Biver et al. (2019) Biver N., et al., 2019, Astron. Astrophys., 630, A19
- Blum et al. (2017) Blum J., et al., 2017, Mon. Not. R. Astron. Soc., 469, S755
- Bockelée-Morvan et al. (2001) Bockelée-Morvan D., Lellouch E., Biver N., Paubert G., Bauer J., Colom P., Lis D. C., 2001, Astron. Astrophys., 377, 343
- Bockelée-Morvan et al. (2004) Bockelée-Morvan D., Crovisier J., Mumma M. J., Weaver H. A., 2004, in Festou M. C., Keller H. U., Weaver H. A., eds, , Comets II. Univ. of Arizona Press, Tucson, pp 391–423
- Boss & Vanhala (2000) Boss A., Vanhala H. A. . T., 2000, Space Sci. Rev., 92, 13
- Boss et al. (2008) Boss A. P., Ipatov S. I., Keiser S. A., Myhill E. A., Vanhala H. A. T., 2008, Astrophys. J., 686, L119
- Bottke et al. (2012) Bottke W. F., Vokrouhlický D., Minton D., Nesvorný D., Morbidelli A., Brasser R., Simonson B., Levison H. F., 2012, Nature, 485, 78
- Brasser & Morbidelli (2013) Brasser R., Morbidelli A., 2013, Icarus, 225, 40
- Brasser et al. (2007) Brasser R., Duncan M. J., Levison H. F., 2007, Icarus, 191, 413
- Brownlee et al. (2006) Brownlee D., et al., 2006, Science, 314, 1711
- Cabral et al. (2019) Cabral N., et al., 2019, Astron. Astrophys., 621, A102
- Calmonte et al. (2016) Calmonte U., et al., 2016, Mon. Not. R. Astron. Soc., 462, S253
- Capria et al. (1996) Capria M. T., Capaccioni F., Coradini A., De Sanctis M. C., Espinasse S., Federico C., Orosei R., Salomone M., 1996, Planet. Space Sci., 44, 987
- Capria et al. (2000) Capria M. T., Coradini A., De Sanctis M. C., Orosei R., 2000, Astron. J., 119, 3112
- Capria et al. (2001) Capria M. T., Coradini A., De Sanctis M. C., Blecka M. I., 2001, Planet. Space Sci., 49, 907
- Capria et al. (2002) Capria M. T., Coradini A., De Sanctis M. C., 2002, Earth, Moon, Planets, 90, 217
- Capria et al. (2009) Capria M. T., Coradini A., De Sanctis M. C., Mazzotta Epifani E., Palumbo P., 2009, Astron. Astrophys., 504, 249
- Castillo-Rogez et al. (2009) Castillo-Rogez J. C., Johnson T. V., Lee M. H., Turner N. J., Matson D. L., 2009, Icarus, 204, 658
- Castillo-Rogez et al. (2012) Castillo-Rogez J. C., Johnson T. V., Thomas P. C., Choukroun M., Matson D. L., Lunine J. I., 2012, Icarus, 219, 86
- Cercignani (2000) Cercignani C., 2000, Rarefied gas dynamics. From basic concepts to actual calculations. Cambridge University Press, Cambridge
- Choi et al. (2002) Choi Y.-J., Cohen M., Merk R., Prialnik D., 2002, Icarus, 160, 300
- Choukroun et al. (2020) Choukroun M., et al., 2020, Space Sci. Rev., 216, 44
- Combi et al. (2020) Combi M., et al., 2020, Icarus, 335, 113421
- Cooke et al. (2018) Cooke I. R., Öberg K. I., Fayolle E. C., Peeler Z., Bergner J. B., 2018, Astrophys. J., 852, 75
- Davidsson (2008) Davidsson B. J. R., 2008, Space Sci. Rev., 138, 207
- Davidsson & Gutiérrez (2004) Davidsson B. J. R., Gutiérrez P. J., 2004, Icarus, 168, 392
- Davidsson & Gutiérrez (2005) Davidsson B. J. R., Gutiérrez P. J., 2005, Icarus, 176, 453
- Davidsson & Gutiérrez (2006) Davidsson B. J. R., Gutiérrez P. J., 2006, Icarus, 180, 224
- Davidsson & Skorov (2002) Davidsson B. J. R., Skorov Y. V., 2002, Icarus, 159, 239
- Davidsson & Skorov (2004) Davidsson B. J. R., Skorov Y. V., 2004, Icarus, 168, 163
- Davidsson et al. (2016) Davidsson B. J. R., et al., 2016, Astron. Astrophys., 592, A63
- Davidsson et al. (2021a) Davidsson B. J. R., Samarasinha N., Farnocchia D., Gutiérrez P., 2021a, Mon. Not. R. Astron. Soc. (In preparation)
- Davidsson et al. (2021b) Davidsson B. J. R., et al., 2021b, Icarus, 354, 114004
- De Sanctis et al. (1999) De Sanctis M. C., Capaccioni F., Capria M. T., Coradini A., Federico C., Orosei R., Salomone M., 1999, Planet. Space Sci., 47, 855
- De Sanctis et al. (2000) De Sanctis M. C., Capria M. T., Coradini A., Orosei R., 2000, Astron. J., 120, 1571
- De Sanctis et al. (2001) De Sanctis M. C., Capria M. T., Coradini A., 2001, Astron. J., 121, 2792
- De Sanctis et al. (2007) De Sanctis M. C., Capria M. T., Coradini A., 2007, Mem. S. A. It. Suppl, 11, 135
- De Sanctis et al. (2015) De Sanctis M. C., et al., 2015, Nature, 525, 500
- Dello Russo et al. (2000) Dello Russo N., Mumma M. J., DiSanti M. A., Magee-Sauer K., Novak R., Rettig T. W., 2000, Icarus, 143, 324
- Di Sisto & Brunini (2007) Di Sisto R. P., Brunini A., 2007, Icarus, 190, 224
- Dodson-Robinson et al. (2009) Dodson-Robinson S. E., Willacy K., Bodenheimer P., Turner N. J., Beichman C. A., 2009, Icarus, 200, 672
- Duncan & Levison (1997) Duncan M. J., Levison H. F., 1997, Science, 276, 1670
- Enzian et al. (1997) Enzian A., Cabot H., Klinger J., 1997, Astron. Astrophys., 319, 995
- Enzian et al. (1999) Enzian A., Klinger J., Schwehm G., Weissman P. R., 1999, Icarus, 138, 74
- Espinasse et al. (1991) Espinasse S., Klinger J., Ritz C., Schmitt B., 1991, Icarus, 92, 350
- Espinasse et al. (1993) Espinasse S., Coradini A., Capria M. T., Capaccioni F., Orosei R., Salomone M., Federico C., 1993, Planet. Space Sci., 41, 409
- Fanale & Salvail (1984) Fanale F. P., Salvail J. R., 1984, Icarus, 60, 476
- Fernández (1980) Fernández J. A., 1980, Mon. Not. R. Astr. Soc., 192, 481
- Fernández & Ip (1984) Fernández J. A., Ip W.-H., 1984, Icarus, 58, 109
- Fornasier et al. (2013) Fornasier S., et al., 2013, Astron. Astrophys., 555, A15
- Fulle et al. (2017) Fulle M., et al., 2017, Mon. Not. R. Astron. Soc., 469, S45
- Garenne et al. (2014) Garenne A., Beck P., Montes-Hernandez G., Chiriac R., Toche F., Quirico E., Bonal L., Schmitt B., 2014, Geochim. Cosmochim. Acta, 137, 93
- Gasc et al. (2017) Gasc S., et al., 2017, Mon. Not. R. Astron. Soc., 469, S108
- Gerakines et al. (1999) Gerakines P. A., et al., 1999, Astrophys. J., 522, 357
- Ghormley (1968) Ghormley J. A., 1968, J. Chem. Phys., 48, 503
- Giauque & Egan (1937) Giauque W. F., Egan C. J., 1937, J. Chem. Phys., 5, 45
- Giauque & Stout (1936) Giauque W. F., Stout J. W., 1936, J. Am. Chem. Soc., 58, 1144
- Goldin et al. (1949) Goldin A. S., Knight G. B., Macklin P. A., Macklin R. L., 1949, Phys. Rev., 76, 336
- Gomes et al. (2004) Gomes R. S., Morbidelli A., Levison H. F., 2004, Icarus, 170, 492
- Gomes et al. (2005) Gomes R. S., Levison H. F., Tsiganis K., Morbidelli A., 2005, Nature, 435, 466
- González et al. (2008) González M., Gutiérrez P. J., Lara L. M., Rodrigo R., 2008, Astron. Astrophys., 486, 331
- González et al. (2014) González M., Gutiérrez P. J., Lara L. M., 2014, Astron. Astrophys., 563, A98
- Gortsas et al. (2011) Gortsas N., Kührt E., Motschmann U., Keller H. U., 2011, Icarus, 212, 858
- Grün et al. (1991) Grün E., et al., 1991, in R. L. Newburn J., Neugebauer M., Rahe J., eds, , Vol. 1, Comets in the Post–Halley era. Kluwer Academic Publishers, pp 277–297
- Grün et al. (1993) Grün E., et al., 1993, J. Geophys. Res., 98, 15091
- Grundy et al. (2007) Grundy W. M., et al., 2007, Icarus, 191, 286
- Guilbert-Lepoutre (2011) Guilbert-Lepoutre A., 2011, Astron. J., 141, 103
- Guilbert-Lepoutre (2012) Guilbert-Lepoutre A., 2012, Astron. J., 144, 97
- Guilbert-Lepoutre et al. (2011) Guilbert-Lepoutre A., Lasue J., Federico C., Coradini A., Orosei R., Rosenberg E. D., 2011, Astron. Astrophys., 529, A71
- Gundlach et al. (2015) Gundlach B., Blum J., Keller H. U., Skorov Y. V., 2015, Astron. Astrophys., 583, A12
- Gutiérrez et al. (2000) Gutiérrez P. J., Ortiz J. L., Rodrigo R., López-Moreno J. J., 2000, Astron. Astrophys., 355, 809
- Gutiérrez et al. (2001) Gutiérrez P. J., Ortiz J. L., Rodrigo R., López-Moreno J. J., 2001, Astron. Astrophys., 374, 326
- Güttler et al. (2009) Güttler C., Krause M., Geretshauser R. J., Speith R., Blum J., 2009, Astron. J., 701, 130
- Haisch et al. (2001) Haisch K. E., Lada E. A., Lada C. J., 2001, Astrophys. J., 553, L153
- Hales et al. (2019) Hales A. S., Gorti U., Carpenter J. M., Hughes M., Flaherty K., 2019, Astrophys. J., 878, 113
- Henke et al. (2012) Henke S., Gail H.-P., Trieloff M., Schwarz W. H., Kleine T., 2012, Astron. Astrophys., 537, A45
- Hevey & Sanders (2006) Hevey P. J., Sanders I. S., 2006, Meteo. Planet. Sci., 41, 95
- Hoang et al. (2020) Hoang M., et al., 2020, Astron. Astrophys., 638, A106
- Höfner et al. (2017) Höfner S., et al., 2017, Astron. Astrophys., 608, A121
- Horai (1971) Horai K.-I., 1971, J. Geophys. Res., 76, 1278
- Horner et al. (2004) Horner J., Evans N. W., Bailey M. E., 2004, Mon. Not. R. Astron. Soc., 354, 798
- Hörst & Tolbert (2013) Hörst S. M., Tolbert M. A., 2013, Astrophys. J. Lett., 770, L10
- Hu et al. (2017) Hu X., et al., 2017, Astron. Astrophys., 604, A114
- Huebner et al. (1999) Huebner W. F., Benkhoff J., Capria M. T., Coradini A., De Sanctis M. C., Enzian A., Orosei R., Prialnik D., 1999, Adv. Space Res., 23, 1283
- Huebner et al. (2006) Huebner W. F., Benkhoff J., Capria M.-T., Coradini A., De Sanctis C., Orosei R., Prialnik D., 2006, Heat and gas diffusion in comet nuclei. ESA Publications Division, Noordwijk, The Netherlands
- Hui et al. (2018) Hui M.-T., Jewitt D., Clark D., 2018, Astron. J., 155, 25
- Hui et al. (2019) Hui M.-T., Farnocchia D., Micheli M., 2019, Astron. J., 157, 162
- Jewitt (2009) Jewitt D., 2009, Astron. J., 137, 4296
- Jewitt & Luu (1993) Jewitt D., Luu J., 1993, Nature, 362, 730
- Jewitt & Matthews (1999) Jewitt D., Matthews H., 1999, Astron. J., 117, 1056
- Jewitt et al. (2008) Jewitt D., Garland C. A., Aussel H., 2008, Astron. J., 135, 400
- Jewitt et al. (2017) Jewitt D., Hui M.-T., Mutchler M., Weaver H., Li J., Agarwal J., 2017, Astron. J. Lett., 847, L19
- Jewitt et al. (2019) Jewitt D., Agarwal J., Hui M.-T., Li J., Mutchler M., Weaver H., 2019, Astron. J., 157, 65
- Jewitt et al. (1992) Jewitt D., Luu J., Marsden B. G., 21992, IAU Circ., 5611, 1
- Johansen & Mac Low (2009) Johansen A., Mac Low A. N. Y. M.-M., 2009, Astrophys. J., 704, L75
- Johansen et al. (2007) Johansen A., Oishi J. S., Low M.-M. M., Klahr H., Henning T., Youdin A. N., 2007, Nature, 448, 1022
- Johnson & Lunine (2005) Johnson T. V., Lunine J. I., 2005, Nature, 435, 69
- Keller et al. (2015) Keller H. U., et al., 2015, Astron. Astrophys., 583, A34
- Keller et al. (2017) Keller H. U., et al., 2017, Mon. Not. R. Astron. Soc., 469, S357
- Klinger (1980) Klinger J., 1980, Science, 209, 271
- Klinger (1981) Klinger J., 1981, Icarus, 47, 320
- Kömle & Dettleff (1991) Kömle N. I., Dettleff G., 1991, Icarus, 89, 73
- Kömle & Steiner (1992) Kömle N. I., Steiner G., 1992, Icarus, 96, 204
- Kossacki et al. (1994) Kossacki K. J., Kömle N. I., Kargl G., Steiner G., 1994, Planet. Space Sci., 42, 383
- Kossacki et al. (1997) Kossacki K. J., Kömle N. I., Leliwa-Kopystyński J., Kargl G., 1997, Icarus, 128, 127
- Kossacki et al. (2015) Kossacki K. J., Spohn T., Hagermann A., Kaufmann E., Kührt E., 2015, Icarus, 260, 464
- Kral et al. (2020) Kral Q., Matrá L., Kennedy G. M., Marino S., Wyatt M. C., 2020, Mon. Not. R. Astron. Soc., 497, 2811
- Kührt (1984) Kührt E., 1984, Icarus, 60, 512
- Kührt & Keller (1994) Kührt E., Keller H. U., 1994, Icarus, 109, 121
- Kulyk et al. (2018) Kulyk I., Rousselot P., Korsun P. P., Afanasiev V. L., Sergeev A. V., Velichko S. F., 2018, Astron. Astrophys., 611, A32
- Larson (2003) Larson R. B., 2003, Rep. Prog. Phys., 66, 1651
- Levison & Duncan (1997) Levison H. F., Duncan M. J., 1997, Icarus, 127, 13
- Levison et al. (2008) Levison H. F., Morbidelli A., van Laerhoven C., Gomes R., Tsiganis K., 2008, Icarus, 196, 258
- Levison et al. (2010) Levison H. F., Duncan M. J., Brasser R., Kaufmann D. E., 2010, Science, 329, 187
- Levison et al. (2011) Levison H. F., Morbidelli A., Tsiganis K., Nesvorný D., Gomes R., 2011, Astron. J., 142, 152
- Li et al. (2020) Li J., Jewitt D., Mutchler M., Agarwal J., Weaver H., 2020, Astron. J., 159, 209
- Lieman-Sifry et al. (2016) Lieman-Sifry J., Hughes A. M., Carpenter J. M., Gorti U., Hales A., Flaherty K. M., 2016, Astrophys. J., 828, 25
- Lodders (2003) Lodders K., 2003, Astrophys. J., 591, 1220
- Lorek et al. (2016) Lorek S., Gundlach B., Lacerda P., Blum J., 2016, Astron. Astrophys., 587, A128
- Lovell et al. (2020) Lovell J. B., et al., 2020, Mon. Not. R. Astron. Soc., TBD, TBD
- Luna et al. (2008) Luna R., Millán C., Domingo M., Satorre M. Á., 2008, Astrophys. Space Sci., 314, 113
- Luspay-Kuti et al. (2015) Luspay-Kuti A., et al., 2015, Astron. Astrophys., 583, A4
- Luu & Jewitt (1990) Luu J., Jewitt D., 1990, Astron. J., 100, 913
- MacPherson et al. (1995) MacPherson G. J., Davis A. N., Zinner E. K., 1995, Meteoritics, 30, 365
- Malhotra (1993) Malhotra R., 1993, Nature, 365, 819
- Marboeuf et al. (2010) Marboeuf U., Mousis O., Petit J.-M., Schmitt B., 2010, Astrophys. J., 708, 812
- Marboeuf et al. (2011) Marboeuf U., Mousis O., Petit J.-M., Schmitt B., Cochran A. L., Weaver H. A., 2011, Astron. Astrophys., 525, A144
- Marchi et al. (2013) Marchi S., et al., 2013, Nature Geosci., 6, 303
- Marino et al. (2016) Marino S., et al., 2016, Mon. Not. R. Astron. Soc., 460, 2933
- Markiewicz et al. (1998) Markiewicz W. J., Skorov Y. V., Keller H. U., Kömle N. I., 1998, Planet. Space. Sci., 46, 357
- Marty et al. (2017) Marty B., et al., 2017, Science, 356, 1069
- Masiero et al. (2021) Masiero J. R., Davidsson B. J. R., Liu Y., Moore K., Tuite M., 2021, Planet. Sci. J.
- Masset & Snellgrove (2001) Masset F., Snellgrove M., 2001, Mon. Not. R. Astron. Soc., 320, L55
- Matrá et al. (2017a) Matrá L., et al., 2017a, Mon. Not. R. Astron. Soc., 464, 1415
- Matrá et al. (2017b) Matrá L., et al., 2017b, Astrophys. J., 842, 9
- Matrá et al. (2019) Matrá L., Öberg K. I., Wilner D. J., Olofsson J., Bayo A., 2019, Astron. J., 157, 117
- Matson et al. (2009) Matson D. L., Castillo-Rogez J. C., Schubert G., Sotin C., McKinnon W. B., 2009, in Dougherty M. K., Esposito L. W., Krimigis S. M., eds, , Saturn from Cassini–Huygens. Springer Science, pp 577–612
- McKinnon et al. (2017) McKinnon W. B., et al., 2017, Icarus, 287, 2
- Meech et al. (1997) Meech K. J., Buie M. W., Samarasinha N. H., Mueller B. E. A., Belton M. J. S., 1997, Astron. J., 113, 844
- Meech et al. (2017a) Meech K. J., et al., 2017a, Astrophys. J., 153, 206
- Meech et al. (2017b) Meech K. J., et al., 2017b, Astrophys. J. Lett., 849, L8
- Mekler et al. (1990) Mekler Y., Prialnik D., Podolak M., 1990, Astrophys. J., 356, 682
- Merk & Prialnik (2003) Merk R., Prialnik D., 2003, Earth Moon Planets, 92, 359
- Merk & Prialnik (2006) Merk R., Prialnik D., 2006, Icarus, 183, 283
- Moór et al. (2015) Moór A., et al., 2015, Astrophys. J., 814, 42
- Morbidelli & Crida (2007) Morbidelli A., Crida A., 2007, Icarus, 191, 158
- Morbidelli et al. (2005) Morbidelli A., Levison H. F., Tsiganis K., Gomes R., 2005, Nature, 435, 462
- Morbidelli et al. (2007) Morbidelli A., Tsiganis K., Crida A., Levison H. F., Gomes R., 2007, Astron. J., 134, 1790
- Morbidelli et al. (2012) Morbidelli A., Marchi S., Bottke W. F., Kring D., 2012, Earth Planet. Sci. Lett., 355, 144
- Mousis et al. (2017) Mousis O., et al., 2017, Astrophys. J. Lett., 839, L4
- Muntean et al. (2016) Muntean E. A., Lacerda P., Field T. A., Fitzsimmons A., Fraser W. C., Hunniford A. C., McCullough R. W., 2016, Mon. Not. R. Astron. Soc., 462, 3361
- Nesvorný & Morbidellli (2012) Nesvorný D., Morbidellli A., 2012, Astron. J., 144, 117
- Nesvorný et al. (2010) Nesvorný D., Youdin A. N., Richardson D. C., 2010, Astron. J., 140, 785
- Norris et al. (1983) Norris T. L., Gancarz A. J., Rokop D. J., Thomas K. W., 1983, J. Geophys. Res., 88, B331
- O’Rourke et al. (2020) O’Rourke L., et al., 2020, Nature, 586, 697
- Öberg et al. (2009) Öberg K. I., Fayolle E. C., Cuppen H. M., van Dishoeck E. F., Linnartz H., 2009, Astron. Astrophys., 505, 183
- Orosei et al. (1995) Orosei R., Capaccioni F., Capria M. T., Coradini A., Espinasse S., Federico C., Salomone M., Schwehm G. H., 1995, Astron. Astrophys., 301, 613
- Orosei et al. (1999) Orosei R., Capaccioni F., Capria M. T., Coradini A., De Sanctis M. C., Federico C., Salomone M., Huot J.-P., 1999, Planet. Space Sci., 47, 839
- Palla & Stahler (1993) Palla F., Stahler S. W., 1993, Astrophys. J., 418, 414
- Parker & Kavelaars (2010) Parker A. H., Kavelaars J. J., 2010, Astrophys. J. Lett., 722, L204
- Parker et al. (2011) Parker A. H., Kavelaars J. J., Petit J.-M., Jones L., Gladman B., Parker J., 2011, Astrophys. J., 743, 1
- Pätzold et al. (2019) Pätzold M., et al., 2019, Mon. Not. R. Astron. Soc., 483, 2337
- Podolak & Prialnik (1996) Podolak M., Prialnik D., 1996, Planet. Space Sci., 44, 655
- Podolak & Prialnik (2006) Podolak M., Prialnik D., 2006, in Thomas P. J., Chyba C., McKay C., eds, , Comets and the origin and evolution of life. Springer–Verlag, Berlin
- Pontoppidan et al. (2008) Pontoppidan K. M., et al., 2008, Astrophys. J., 678, 1005
- Preusker et al. (2015) Preusker F., et al., 2015, Astron. Astrophys., 583, A33
- Prialnik (1992) Prialnik D., 1992, Astrophys. J., 388, 196
- Prialnik & Bar-Nun (1987) Prialnik D., Bar-Nun A., 1987, Astrophys. J., 313, 893
- Prialnik & Bar-Nun (1988) Prialnik D., Bar-Nun A., 1988, Icarus, 74, 272
- Prialnik & Bar-Nun (1990) Prialnik D., Bar-Nun A., 1990, Astrophys. J., 363, 274
- Prialnik & Podolak (1995) Prialnik D., Podolak M., 1995, Icarus, 117, 420
- Prialnik & Podolak (1999) Prialnik D., Podolak M., 1999, Space Sci. Rev., 90, 169
- Prialnik et al. (1987) Prialnik D., Bar-Nun A., Podolak M., 1987, Astrophys. J., 319, 993
- Prialnik et al. (2004) Prialnik D., Benkhoff J., Podolak M., 2004, in Festou M. C., Keller H. U., Weaver H. A., eds, , Comets II. Univ. of Arizona Press, Tucson, pp 359–387
- Prialnik et al. (2008) Prialnik D., Sarid G., Rosenberg E. D., Merk R., 2008, Space Sci. Rev., 138, 147
- Rickman & Fernández (1986) Rickman H., Fernández J. A., 1986, in , The Comet Nucleus Sample Return Mission. ESA SP–249, pp 185–194
- Rickman et al. (1990) Rickman H., Fernández J. A., Gustafson B. Å. S., 1990, Astron. Astrophys., 237, 524
- Robie et al. (1982) Robie R. A., Hemingway B. S., Takei H., 1982, American Mineralogist, 67, 470
- Rosenberg & Prialnik (2009) Rosenberg E. D., Prialnik D., 2009, Icarus, 201, 740
- Rosenberg & Prialnik (2010) Rosenberg E. D., Prialnik D., 2010, Icarus, 209, 753
- Santos-Sanz et al. (2012) Santos-Sanz P., et al., 2012, Astron. Astrophys., 541, A92
- Sarid & Prialnik (2009) Sarid G., Prialnik D., 2009, Meteo. Planet. Sci, 44, 1905
- Sárneczky et al. (2016) Sárneczky K., et al., 2016, Astron. J., 152, 220
- Satorre et al. (2009) Satorre M. Á., Luna R., Millán C., Santonja C., Cantó J., 2009, Planet. Space Sci., 57, 250
- Schmitt et al. (1989) Schmitt B., Espinasse S., Grim R. J. A., Greenberg J. M., Klinger J., 1989, ESA SP, 302, 65
- Senay & Jewitt (1994) Senay M. C., Jewitt D., 1994, Nature, 371, 229
- Shchuko et al. (2014) Shchuko O. B., Shchuko S. D., Kartashov D. V., Orosei R., 2014, Planet. Space Sci., 104, 147
- Shoshani et al. (1997) Shoshani Y., Heifetz E., Prialnik D., Podolak M., 1997, Icarus, 126, 342
- Shoshany et al. (2002) Shoshany Y., Prialnik D., Podolak M., 2002, Icarus, 157, 219
- Sicilia-Aguilar et al. (2006) Sicilia-Aguilar A., et al., 2006, Astrophys. J., 638, 897
- Simon et al. (2019) Simon A., Öberg K. I., Rajappan M., Maksiutenko P., 2019, Astrophys. J., 883, 21
- Skorov & Rickman (1995) Skorov Y. V., Rickman H., 1995, Planet. Space. Sci., 43, 1587
- Skorov et al. (1999) Skorov Y. V., Kömle N. I., Markiewicz W. J., Keller H. U., 1999, Icarus, 140, 173
- Stern (2003) Stern S. A., 2003, Nature, 424, 639
- Stern & Shull (1988) Stern S. A., Shull J. M., 1988, Nature, 332, 407
- Suhasaria et al. (2017) Suhasaria T., Thrower J. D., Zacharias H., 2017, Mon. Not. R. Astron. Soc., 472, 389
- Szabó et al. (2012) Szabó G. M., Kiss L. L., Pál A., Kiss C., Sárneczky K., Juhász A., Hogerheijde M. R., 2012, Astrophys. J., 761, 8
- Tancredi et al. (1994) Tancredi G., Rickman H., Greenberg J. M., 1994, Astron. Astrophys., 286, 659
- Tesfaye Firdu & Taskinen (2010) Tesfaye Firdu F., Taskinen P., 2010, Densities of molten and solid alloys of (Fe, Cu, Ni, Co) – S at elevated temperatures – Literature review and analysis. Aalto University Publications in Materials Science and Engineering, Aalto University School of Science and Technology, Espoo
- Thomas et al. (2015a) Thomas N., et al., 2015a, Science, 347, aaa0440
- Thomas et al. (2015b) Thomas N., et al., 2015b, Astron. Astrophys., 583, A17
- Tscharnuter et al. (2009) Tscharnuter W. M., Schönke J., Gail H.-P., Trieloff M., Lüttjohan E., 2009, Astron. Astrophys., 504, 109
- Tsiganis et al. (2005) Tsiganis K., Gomes R., Morbidelli A., f. Levison H., 2005, Nature, 435, 459
- Walsh et al. (2011) Walsh K. J., Morbidelli A., Raymond S. N., O’Brien D. P., Mandell A. M., 2011, Nature, 475, 206
- Weast (1974) Weast R. C., 1974, Handbook of chemistry and physics, 55 edn. CRC Press, Cleveland
- Weidenschilling (1997) Weidenschilling S. J., 1997, Icarus, 127, 290
- Weissman & Kieffer (1981) Weissman P. R., Kieffer H. H., 1981, Icarus, 47, 302
- Wierzchos et al. (2017) Wierzchos K., Womack M., Sarid G., 2017, Astron. J., 153, 230
- Womack & Stern (1999) Womack M., Stern S. A., 1999, Sol. Sys. Res., 33, 187
- Yasui & Arakawa (2009) Yasui M., Arakawa M., 2009, J. Geophys. Res., 114, E09004
- Yomogida & Matsui (1983) Yomogida K., Matsui T., 1983, J. Geophys. Res., 88, 9513
- Youdin & Goodman (2005) Youdin A. N., Goodman J., 2005, Astrophys. J., 620, 459
- Zuckerman et al. (1995) Zuckerman B., Foreille T., Kastner J. H., 1995, Nature, 373, 494