© W. Neumann et al. 2020

1 Introduction

Ceres is a remarkable example and a unique case study of a water-rich dwarf planet located in the asteroid belt. From the Dawn mission and related research, an unprecedented amount of observational data and a variety of modeling results have become available, allowing for constraints on the origin of Ceres and a more holistic understanding of Ceres’ interior. A mean radius of 470 km classifies this body as a dwarf planet, a potential representative of those that contributed to the early Earth during its accretion (Peslier et al. 2017; Sarafian et al. 2014). A rather low density of ≈2162 kg m⁻³ indicates an ice-rich composition (15−55 vol.% for an average rock density of 2400−3700 kg m⁻³), while geomorphological analyses suggest the presence of liquids in the subsurface in geologically recent times (Ruesch et al. 2016; Fu et al. 2017), supported by a mean surface temperature of T_surf = 168 K (Tosi et al. 2015) with a maximum of T_surf = 235 K. Spectral observations show pervasive aqueous alteration of the surface itself with a variety of hydrated minerals (De Sanctis et al. 2015), indicating the mobility of liquids in the subsurface and in the deeper interior of Ceres in the past. Recent surface activity in the form of water vapor plumes (Küppers et al. 2014) and putative haze, as an interpretation of brightness variations specific to the crater Occator (Thangjam et al. 2016), indicates that the internal activity could be still ongoing. Both exogenic and endogenic sources have been proposed for the altered material and organics detected (Marchi et al. 2019; Pieters et al. 2018). However, since Ceres shows clear signatures of pervasive hydrothermal activity and mobility of liquids (the mountain Ahuna Mons that was proposed to have formed as a cryovolcanic dome), both altered rock and organic-rich areas are probably the products of internal processes (De Sanctis et al. 2017). Detection of surface water ice deposits (Platz et al. 2016; Prettyman et al. 2017; Combe et al. 2019) and the potential presence of ice in Ceres’ subsurface (Bland et al. 2016; Ermakov et al. 2017) suggest geological similarities to differentiated icy bodies (Castillo-Rogez et al. 2020), as predicted by pre-Dawn models (McCord & Sotin 2005; Castillo-Rogez & McCord 2010).

Several lines of evidence indicate that Ceres was at least partially differentiated by internal processes involving water mobility and involved tectonic and cryovolcanic processes. Dawn’s gravity and shape measurements show that Ceres is in hydrostatic equilibrium with a mean moment of inertia of 0.37, implying partial differentiation (Park et al. 2016). Admittance analysis points to an ≈ 1200−1400 kg m⁻³ dense and ≈40 km thick crust, assuming a two-layer model (Ermakov et al. 2017) that overlays an ≈ 2400 kg m⁻³ dense mantle. A decrease of the crust viscosity with depth from a mechanically strong maximum of ≈ 10²⁵ Pa s to < 10²¹ Pa s in the upper ≥60 km of the mantle were inferred by Fu et al. (2017) and attributed putatively to the presence of pore liquids trapped below the crust. Residual sodium and potassium chloride brines and ammonia in marginally small amounts from the freezing of a global ocean were suggested by recent geochemical modeling as candidates for these fluids (Castillo-Rogez et al. 2019).

The structure suggested by recent modeling is at least partially differentiated and has multiple layers: a thin crust of dirty ice (rock and salts), a frozen ocean, a liquid basal water layer, a hydrated rocky outer core and a dry rocky inner core (e.g., Neveu et al. 2015). This is, however, challenged by further interior models that match an only weakly differentiated structure with observations (Castillo-Rogez et al. 2016, 2019) or infer no differentiation at all (Zolotov 2009, 2020), and by density constraints derived from gravity data that argue against silicate dehydration in the interior (Ermakov et al. 2017).

Another question that is under debate pertains to the circumstances of Ceres’ formation. From the current state of the art, the following scenarios can be summarized (McKinnon 2012; Castillo-Rogez et al. 2013). Ceres could have accreted within ≈ 5 millions of years near or at its present orbit in the asteroid belt (AB; Dodson-Robinson et al. 2009), or on a similar timescale in the inter-planetary part of the accretion disk (IPD, inter-planet disk between the orbits of Jupiter and Neptune) and migrated into the asteroid belt later, according to recent scattering models (“grand tack” scenario, e.g., Raymond & Izidoro 2017, Walsh et al. 2012). Lastly, it has also been suggested that Ceres could have accreted on a much slower timescale in the Kuiper belt (KB) and migrated to the inner solar system during the late heavy bombardment phase in the Nice model (McKinnon 2012; Castillo-Rogez et al. 2013). Table 1 illustrates that these accretion scenarios result in distinct thermo-chemical evolution paths and influence, in particular, the timing of the water-rock differentiation and the retention of water (the main parameters here are the accretion duration, the ice content and the current temperature). It should be noted that current understanding of proto-Kuiper belt objects converges on the accretion within the first 10 Myr of solar systemevolution (e.g., Morbidelli & Nesvorny 2019; Nesvorny et al. 2019) in contrast to the extended accretion timescale of up to 100 Myr considered previously (Castillo-Rogez et al. 2013; McKinnon 2012).

Further implications arise for the alteration and metamorphism of the material of smaller planetesimals that accreted into Ceres. Assuming the asteroid belt or inter-planet disk scenario, it has been proposed that these small objects were not able to evolve far enough to experience aqueous alteration prior to their accretion onto Ceres (Castillo-Rogez et al. 2013). On the other hand, if accretion took place in the Kuiper Belt, the extended accretion duration of Ceres would enable these bodiesto grow to considerable sizes, but the lack of short-lived radionuclides and a higher ice fraction, as well as low ambient temperatures, would not allow for aqueous alteration either (Castillo-Rogez et al. 2017). This impliesthat aqueous alteration would have taken place on Ceres itself after accretion and not on its smaller precursors. However, the present study shows that aqueous alteration of planetesimals is possible for all three accretion scenarios.

Aqueous alteration is a series of chemical reactions that result in the production of hydrated minerals (e.g., serpentine) from liquid water in contact with a “dry” reactant (e.g., olivine). Experimental studies on the CV chondrite Allende showed that aqueous alteration of its minerals requires weeks or even days at the temperatures of T > 420 K (Jones & Brearley 2006). These experiments observed the very earliest stages of the alteration process under conditions close to those experienced by carbonaceous chondrites, producing alteration assemblage resembling most closely CI chondrites and establishing the general alteration behavior of chondritic materials in the presence of aqueous fluids. As the rate of serpentinization is a strong function of temperature, the time necessary for water to remain in contact with the reactant increases for lower temperatures. For instance, at T = 273 K, serpentinization completes over approximately 200 yr as calculated with a time reaction equation derived by Wakita & Sekiya (2011) from the results of Jones & Brearley (2006). Investigations of carbonate formation in CM chondrites indicate that aqueous alteration of the CM parent body occurred in the temperature range of 293–344 K (Guo & Eiler 2007). For the CI parent body, lower temperatures of 243–338 K were even derived (Guo 2009) from the experimental study of the Orgueil CI chondrite. For T < 273 K, the presence of liquid requires chemical compounds that act as antifreeze, such as ammonia. For comparison, a temperature interval of 320–520 K has been suggested by Allen & Seyfried (2004) for the alteration of oceanic peridotites.

Although an initially ice-rich object can experience the loss of a substantial fraction of water caused by intense heating by, for example, radioactive nuclides or by impacts, it is obvious that Ceres retained at least a part (if not nearly all) of its original water content. Among the variety of criteria that influence water retention, the most relevant ones considered here are the accretion time scale, the ice content (quantified in terms of dust/ice ratio) and the ambient nebula temperature T_amb (the initial temperatureT_init = T_amb, and the surface temperature T_surf calculated from T_amb). They are essential for the thermal evolution and drive the course and time scale of the differentiation and, eventually, the present internal structure. The accretion scenarios for Ceres determine these parameters and direct the outcome of a thermal model along a specific track. The parameters T_init and T_amb have higher values for the asteroid belt and for the inter-planet disk and lower by up to 150 K for the Kuiper belt. The accretion time scale determines the original ²⁶Al and ⁶⁰Fe content and the interplay of heating and cooling during accretion via the evolution of the surface-to-volume ratio. The accretion of planetesimals and protoplanets was very fast in the inner solar system and proceeded more slowly in the outer solar system. Therefore, ²⁶Al and ⁶⁰Fe have a substantial effect in the first case and a less important effect in the second. Furthermore, the formation distance influences heating via the accreted ice fraction and the regime of the differentiation (water percolation for a dust-dominated composition, Stokes settling of dust for an ice-dominated composition).

In the present study, the internal evolution of Ceres is calculated from one-dimensional differentiation and thermal modeling. We elaborate on results that were already presented at several conferences (Neumann et al. 2016b, 2017a,b) and account for the latest results from the Dawn mission. The calculations consider three sets of parameters that correspond to the accretion scenarios of Ceres described above. The simulations track the entire evolution, starting with the accretion of the small precursors from dust and ice particles, followed by Ceres’ accretion, through its differentiation and a prolonged cooling phase, producing, finally, the present structural and thermal state of this dwarf planet for comparison with the Dawn observations summarized above.

2 Modeling aspects

2.1 Accretion

Theoretical accretion models predict that big planetesimals with R ≫ 10 km formed within a few $\mathcal{O}({10}^5)$ yr, at most, from an initial swarm of much smaller planetesimals (e.g., Nagasawa et al. 2007; Weidenschilling 2011). This means that the duration of the accretion was much shorter than the initial heating time scale (one half-life of ²⁶Al). These planetesimals must have reached their final sizes before substantial temperature increase could have occurred (e.g., Vernazza et al. 2014). Therefore, for those bodies, it would suffice to assume an “instantaneous” formation at some instant t₀ as an approximation, as it is usually done in numerical models of meteorite parent bodies. Henke et al. (2013) showed that at least for the H chondrite parent body, the thermo-chronological data are best fitted by models with short accretion duration of ≤ 0.3 Myr, and a similar conclusion was drawn by Monnereau et al. (2013). However, the accretion duration could have been longer for water-rich planetesimals that accreted farther away from the Sun. Particularly for such bodies, thermal and structural evolution during an extended accretion can be quite important (Neumann et al. 2015; Dehant & et al. 2016). It should be noted that accretion in a gas-bearing disk that quantitatively matches observations of Kuiper belt object size distributions and obliquities has come forward as the preferred scenario (e.g., Morbidelli & Nesvorny 2019; Nesvorny et al. 2019). The need for nebular gas implies that accretion took place within the first 10 Myr of solar system evolution, that is, within a few half-lives of short-lived radionuclides. In contrast with most previous thermal evolution models of Ceres that approximated accretion within a short time period relative to the heating time scale by instantaneous formation, we explore the effect of continuous accretion (see Sect. 3.1 of this paper).

2.2 Heating

For the early solar system’s planetesimals, the short-lived isotopes ²⁶Al and ⁶⁰Fe are believed to be the dominating heat sources. Their presence in the parent bodies of all meteorite classes is ascertained by the pervasive presence of their decay products (e.g., Dauphas & Chaussidon 2011; Kita et al. 2013). The heat produced by these nuclides could produce strong thermal metamorphism observed in chondrites and primitive achondrites in bodies larger than 10 km (Ghosh et al. 2006). If only heat conduction is assumed, a sufficient net amount of heat is produced to melt larger planetesimals completely (e.g., Neumann et al. 2015). The long-lived radionuclides ⁴⁰K, ²³²Th, ²³⁵U and ²³⁸U dominate the energy production after almost complete decay of the short-lived ones, approximately 10− 20 Myr after CAIs. They remain active on timescales of the age of the solar system and provide a nearly constant heating source relative to that provided by short-lived radionuclides. However, in contrast with planets and largest asteroids, their total contribution to the heating of the relatively small chondrite parent bodies is much less significant than that provided by short-lived isotopes, and, therefore, of marginal importance. Furthermore, any substantial heating by impacts can be ruled out for sub-Vesta sized planetesimals (see, e.g., Merk et al. 2002, and references therein).

2.3 Composition

Since no space rocks originating from Ceres have yet been identified in the collections, constraints on the chemical, mineralogical, and mechanical properties of Ceres’ material are limited to and based on spectral observations, topological properties, etc. Pre-Dawn telescopic observations and models suggested some aqueous alteration on Ceres, including the formation of clay-like materials (Rivkin et al. 2011), and possibly salts covered by a regolith layer of only small-scale compositional variations (McCord & Sotin 2005; Castillo-Rogez & McCord 2010). Observations by Dawn confirmed aqueous alteration and presence of hydrated clay minerals and salts. Analyses of Ceres’ near-infrared spectra indicate that the bulk surface composition consists of ammoniated phyllosilicates (NH₄ -montmorillonite or NH₄-annite, Mg-phyllosilicates, carbonates, and a dark component) and does not show spectral signatures of water ice at a global level (De Sanctis et al. 2015). The mineral textures of Ceres’ surface are characterized by grains with a varying size of $\mathcal{O}(1)-\mathcal{O}({10}^2)$ μm. The mineralogical composition of the surface indicates a chondritic origin, where the closest match is provided by the volatile-rich CI and CM carbonaceous chondrites (De Sanctis et al. 2015). Therefore, the pristine material of the small planetesimals that accreted into Ceres must have been a chondritic precursor.

A density of ≈2160 kg m⁻³ indicates an ice-rich composition while pervasive aqueous alteration of the surface itself, with a variety of hydrated minerals, is clear evidence for the mobility of liquids in the subsurface and in the deeper interior of Ceres in the past. This prospect is supported by the detection of surface water ice deposits that are preserved in permanently shadowed areas in the northern hemisphere of Ceres by Dawn (Platz et al. 2016; Prettyman et al. 2017). Insufficient relaxation of Ceres’ largest craters challenges the pre-Dawn expectations for an ice-rich subsurface (Bland et al. 2016). Crater relaxation simulations suggest a subsurface that is by three orders of magnitude more viscous than water ice, and imply that the shallow subsurface is < 30−40% ice and >60−70% rock, salts and potentially clathrates by volume (Bland et al. 2016). Based on multiple surface sites enriched in carbonates and ammonium chlorides (De Sanctis et al. 2016), salt compounds are abundant in Ceres’ crust and are possibly associated with other chlorides (Castillo-Rogez et al. 2018). The presence of ice within Ceres’ crust (Bland et al. 2016) raises the prospect of geological processes similar to differentiated icy bodies (McCord & Sotin 2005). The potential for a relatively thin hydrosphere on Ceres could help enhance the prospect for material transfer from the deep interior and for surface expression of endogenic activity.

2.4 Initially porous micro-structure and compaction

It is widely accepted that the planetesimals that evolved to form small planetary objects observed in the solar system accreted as porous objects with porosity present as empty pores (i.e., “dust porosity”) from fine dust and ice grains in the protoplanetary accretion disk (Ormel et al. 2007). In the course of thermochemical evolution, some of them underwent processes that caused partial consolidation, as indicated by meteorites. The micro-structure of chondritic and primitive achondritic meteorites and in situ observations of near-Earth asteroids confirm the partial consolidation (Neumann et al. 2018a; Grott et al. 2019; Jaumann et al. 2019). Compaction occurred first by “cold” isostatic pressing at low temperatures and then by creep (“hot” pressing) at higher temperatures that enabled plastic flow of the material (Henke et al. 2012; Neumann et al. 2014b). A higher porosity occurs typically in smaller objects and in the surface layers of larger objects. Thermal evolution models infer that a porous surface layer with a thickness of several to several tens of kilometers was retained due to sub-critical temperatures (and pressures) for compaction (Neumann et al. 2013, 2014b). A chemically homogeneous but porous interior has been suggested for Ceres (Zolotov 2020). Regardless of the present porosity, compaction was one of key processes during the accretion and differentiation of Ceres. However, most thermal evolution models for Ceres do not consider porosity.

Retention of pores in a material with some specific composition is essentially controlled by the thermal evolution, that is, by its accretion time relative to the formation of CAIs. An early accretion within less than one half-life of ²⁶Al of 0.717 Myr would implysufficient heating for compaction to already have occurred to some extent on the smaller precursors (Henke et al. 2013). However, icy planetesimals, in particular, KBOs, likely accreted later (McKinnon et al. 2008) and experienced less heating. Therefore, weakly heated icy precursors of Ceres must have retained substantial porosity prior to their accretion onto Ceres (see also our calculations for such bodies in Sect. 4 of this paper), and Ceres most likely had a high degree of the initial dust porosity (the later the end of the accretion relative to CAIs, the higher it is expected to be). The typical temperatures required for compaction of olivine-dominated material by diffusion creep are ≈ 650−750 K (Yomogida & Matsui 1984; Neumann et al. 2014b), while “soft” materials are mobilized in the low-temperature regime (e.g., creep of serpentine is effective at T < 500 K, see, e.g., Hilairet et al. 2007). At elevated temperature and pressure, the volume fraction of pore space gradually decreases due to creep processes. A typical approach to modeling of creep-driven compaction of planetesimals was introduced by Yomogida & Matsui (1984) and adopted in several studies of small bodies. For more details, we refer to Neumann et al. (2014b, 2015) and Sect. 3 of this paper.

2.5 Constraints on the thermal evolution

Information available on the minimum and maximum temperatures experienced by Ceres can be used to constrain the models of its evolution and therefore the formation variables. The determinative process of these temperatures is the separation of dust and H₂ O. Inefficiency of solid-state separation (e.g., by a very slow settling of dust in an ice layer) and signatures of aqueously altered minerals on Ceres indicate that liquid water was present. Therefore, a very rough minimum requirement is T > 273 K. It is supported by the thermal evolution models that produce a variably high peak temperature for different formation times of Ceres (e.g., McCord & Sotin 2005; Castillo-Rogez & McCord 2010; Neumann et al. 2015; Neveu et al. 2015). On the other hand, an H₂ O-rich present-day interior indicates that conditions that would have led to the loss of water must be avoided. Such a condition would be ebullition and evaporation of the water layer after the differentiation. The peak temperature in the core can be regarded as yet another crude estimate. Combined with the lack of any indications of an inner iron core, the solidus temperature of the silicates of ≈ 1425 K or even the eutectic of the Fe,Ni-FeS system of ≈1261 K is not expected to have been exceeded. At least the temperature corresponding to aqueous alteration must have been kept up over a prolonged period of time, but this temperature varies substantially (see above). Differentiation models indicate that temperatures of up to 800 K are needed for the differentiation in the matrix compaction regime. Presently, warm interior and recent activity indicate by “backward engineering” argumentation that Ceres was warm in the past. The conditions for the upkeep of liquids in the shallow subsurface in the geologically recent past necessitate temperatures of > 273 K for liquid water or even around <200 K for various brines. The freezing point of water containing a mixture of hydrates, carbonates, and chlorides can decrease to 190 K (Kargel et al. 2000; Kargel 1991; Wynn-Williams et al. 2001), while the melting point of ammonia dihydrate of ≈ 176 K is even lower (Desch et al. 2009). A further constraint is the mean moment of inertia of ≈ 0.37 (Park et al. 2016), matching a differentiation degree that is established by certain thermal conditions. This constraint is somewhat weak as well, due to the fact that distinct internal structures and differentiation degrees can produce the same moment of inertia.

2.6 Differentiation

A minimum requirement for the separation of H₂O from the dust component is production of liquid water, typically by the melting of water ice. This is, however, not sufficient for the differentiation by itself. The CI and CM chondrites that are mineralogically closest to Ceres’ surface material are highly porous and have experienced pervasive aqueous alteration on their parent bodies. However, the flow of liquid water was very limited, indicating a closed system at a scale of $\mathcal{O}(100)$ μm and more (Bland et al. 2009). The reason for that is the low permeability of plausible fine-grained chondritic starting materials that is several orders of magnitude smaller than the permeability that was derived from terrestrial analogues and was typically used in numerical models of water-rich asteroids. Bland et al. (2009) predicted water flow to be limited to distances of $<\mathcal{O}(100)$ μm, even in a high porosity, water-saturated asteroid over millions of years if the average matrix grain sizes are $<\mathcal{O}(10)$ μm. An idealized structure with pores that are completely filled with liquid water becomes reasonably more complex if a three-phase system is considered where the pores are only partially filled with water. A fast differentiation and a clean stratification into a rocky core and a water mantle modeled in previous numerical studies of Ceres’ differentiation seem, therefore, to be somewhat naive and idealized.

The differentiation of planetesimals consisting initially of a homogeneous mixture of ice and dust grains could have proceeded along two fundamentally different paths (see Fig. 1 for a visualization and, e.g., Malamud & Prialnik 2013 for further information). Those are (1) settling of dust particles in a layer of liquid water, and (2) upward percolation of liquid water through a rock matrix to the top of a solid-like layer. The first scenario necessitates that the initial composition of an icy planetesimal is dominated by ice (a small dust/ice ratio), while in the second case, it is dominated by the dust (a high dust/ice ratio). The consideration of the initial microstructure of the material provides further details of these scenarios. Here, two conceptual models can be considered (a) the initial porosity is zero, and (b) the initial porosity is > 0. Prior to the melting of ice, the first model has one phase (solid) and two species (ice and dust), and it evolves to a two phase mixture when the ice melts (liquid water and solid dust). The second model initially has two phases (solid material and empty pores) with two species (ice and dust) contained in the solid phase. A third phase (liquid) forms when the water ice melts. The produced water does not fill the pores, because the ice already occupied a volume.

Assuming the model (a) and a small dust/ice ratio, scenario (1) would lead to a quasi-instantaneous formation of a rocky core by Stokes settling of dust particles in a liquid environment. For a high dust/ice ratio, core formation by percolation of water and matrix compaction would occur under high-temperature conditions, allowing the compaction of the matrix on a much larger time scale. If model (b) and a small dust/ice ratio are assumed, the composition is dominated by H₂ O. Therefore, the pores will be closed quasi-instantaneously when the ice melts and differentiation will occur, as in model (a). Finally, if the composition is dominated by the dust, and if the concept of model (b) applies, the differentiation will be the slowest. Here, no percolation will occur as long as there is some empty space (porosity). Therefore, the percolation will be delayed until the pore fraction that is not filled with water becomes reduced to zero. After that, the phase separation will be possible; however, it will be controlled by the matrix deformation under high-temperature conditions. First-order estimates of the timing of differentiation are provided by Neumann et al. (2016a).

The pores evolve to smaller volume due to heating and pressure. This process is frequently called “compaction” (or, more specifically, “sintering”), and it starts before the melting of ice by creep of ice and weak species (e.g., salts, antigorite), proceeding subsequently by creep of more rigid species (e.g., olivine).

Fig. 1 Overview of the main differentiation regimes that arise for different assumptions on the proportions of ice and rock (Neumann et al. 2016a). Left: composition is dominated by ice whose melting results in the formation of an ocean. The rock mantle grows from the settling of rock particles in the Stokes regime and their deposition at the bottom of the ocean. Right: composition is dominated by rock. Melting of ice creates a water network within a porous matrix. The ocean grows due to the percolation of water which is driven by the gravity-induced compaction of the matrix (e.g., McKenzie 1984). Differentiation in this regime requires a longer time scale. In general, planetesimals of equal size but with different ice fractions would differentiate on a different time scale and result in structures with different layer dimensions.

2.7 Physical processes in porous layers

The porosity can be present either as empty pores (usually termed as “dust porosity” in the literature) or as pores that are filled with a fluid (usually termed as “melt porosity” but being, in fact, the volume fraction of the fluid). The former is present from the beginning of the accretion on, and is sustained until the liquid and the solid separate. If the liquid solidifies, no further change of the pore volume occurs in the model. The latter is created when a “melt” is produced, that is, when the ice melts. Therefore, it will exist in an undifferentiated layer during the phase when the dust pores are not closed (i.e., before the differentiation for a dust-dominated material), but the temperature exceeds the melting temperature of water ice. The changes of the phases, that is, the reduction of the pore volume, and the reduction of the solid volume during the increase of the liquid (i.e., melting of water ice) are essential processes for the modeling of the differentiation. Furthermore, when the core is being formed during the phase separation either by the upward percolation of water or by downward settling of dust particles, its structure resembles an agglomerate of grains with interstitial water, that is, it has a melt porosity. These pores are closed in the porous core from the center towards the outer core boundary by creep because temperatures are always highest at the center and decrease toward the surface. The reduction of the melt pore volume, the associated dimensional change of the core, and the increase of the ocean volume are important modeling aspects as well.

2.8 Heat transport in the ocean and in largely solid layers

Under favorable thermal conditions, layers established in consequence of the differentiation undergo processes that transport heat (convection) and material (percolation, Stokes settling), and change locally the chemical and mineralogical composition (dissolution, precipitation, hydration and dehydration). Thereby, the heat fluxes are not carried by conduction alone but by convective processes as well. This can apply, among others, to a water ocean, a solid ice layer, and a porous water-filled outer core. We predict these layers to develop in the course of the model calculations for Ceres.

An ocean that is being heated from below and cooled from above can establish current patterns due to thermal gradients. This mechanism, known as liquid-state convection, provides an effective way to transport mass within the ocean and the energy to the upper boundary of the ocean. Its occurrence and effectiveness can be assessed by calculating the Rayleigh number Ra across the depth of the ocean (Ra > 10³ is necessary, see, e.g., Spohn & Schubert 2003). If the structure of outer part of the core in contact with a liquid ocean has interconnected void spaces, they are filled with the interstitial water. Due to a thermal gradient, cooler water from the ocean can circulate down into the voids, warm up, and return to the ocean. The effectiveness of this type of convective process (termed usually as the hydrothermal circulation) can be estimated as well by the Rayleigh number being super-critical (e.g., Cherkaoui & Wilcock 2001) and it results in mass and energy transfer. Another type of convection (solid-state convection) in the solidified top part of the ocean (or even in the entire ice-rock crust) can also transport more heat than is transported by conduction (see, e.g., Neumann 2019). The transport of the salts dissolved or mineral particles entrained in the water to the subsurface or to the surface can be further consequences.

3 Model features

Numerical calculations have been performed using a numerical model that builds upon the model by Neumann et al. (2015) and now includes, most importantly, the differentiation of water and dust. It solves a 1D heat conduction equation in spherical symmetry considering heating by short- and long-lived radioactive isotopes as well as by impacts, temperature-, porosity- and composition-dependent parameters, and compaction of initially porous material (as well as the compaction of the water-filled rocky core). In addition, it considers melting of water ice and differentiation of a rocky core and water mantle, as well as convection in the water ocean, convection in a solid ice shell, and hydrothermal circulation in a water-filled porous layer. The non-linear set of partial differential equations described further below is solved in both models by an implicit finite difference method that allows for rather big time-steps without running into stability problems. In the following, we provide an overview over the new aspects of the model.

3.1 Accretion

Previous models of Ceres consider instantaneous formation at some time instance relative to CAIs. Neumann et al. (2015) thermal evolution and compaction models, that included continuous accretion, were calculated. Here the increase of Ceres’ dimensions is calculated in the same way by increasing its radius between the time steps until the final radius is reached. The exponential accretion law $$R(t)=R_0\exp \left(\log \left(\frac{R_{\text{f}}}{R_0}\right)\frac{t-t_0}{t_{\text{ac}}}\right)$$(1)

describes the late runaway accumulation of material by impacts of planetesimals that are substantially smaller than the target object modeled. Here, R₀ is the initial radius, R_f is the final radius, t₀ is the time of the accretion onset (relative to CAIs) and t_ac is the duration of the accretion. Starting with the initial radius of R₀ ≈ 100 m, the radius increase is calculated for time t_ac with respect to the accretion law R(t) and to the change of the dust porosity due to compaction by creep. Consequently, the major part of the mass accretes late within the time t_ac : within the last 18% of the accretion duration, around 99% of the mass accretes (for the last 9% of the accretion duration this is 90% of the mass). Free parameters are t₀, t_ac, and the actual form of the accretion law R(t). The former two are varied here according to Table 1, while R(t) is exponential throughout all calculations (it is reasonable that dwarf planets accreted in a runaway scenario, but more details are not available).

3.2 Heating

Using typical data (see, e.g., Henke et al. 2012; Neumann et al. 2015, and Table 2), we include both short-lived radionuclides ²⁶Al and ⁶⁰ Fe and long-lived radionuclides ⁴⁰K, ²³²Th, ²³⁵U and ²³⁸U as heat sources.

We assume that the radionuclides are homogeneously distributed within the dust phase, therefore the amount of the heat available scales with the dust/ice ratio and with the porosity. When water and dust separate, the local heating rates at different depths change because the radionuclides migrate along with the dust, and upon core formation, they become concentrated in the core, while the water ocean and ice layer do not contain any heat sources.

In addition, the heating due to the release of the gravitational energy E_g and the thermal energy E_t by the impacting small planetesimals is included (Šrámek et al. 2012). This energy input is calculated according to the mass ingrowth due to the accretion law R(t). It changes both the surface temperature T_S and the temperature of the impacted body within the burial depth ΔR: $$\begin{array}{lll}\hfill & \hfill & 4\pi \overline{\rho }{H}_{\text{g}}{R}^2\Delta R\approx f\left(\frac{\text{d}{E}_{\text{g}}}{\text{d}t}+\frac{\text{d}{E}_{\text{t}}}{\text{d}t}\right)\hfill \\ \hfill & \hfill & {\sigma }_{\text{SB}}\left(T_{\text{S}}^4-T_{\text{amb}}^4\right)=-k\frac{\text{d}T}{\text{d}r}{|}_{r=R}+\frac{(1-f)}{4\pi R^2}\left(\frac{\text{d}{E}_{\text{g}}}{\text{d}t}+\frac{\text{d}{E}_{\text{t}}}{\text{d}t}\right)\hfill \\ \hfill & \hfill & \frac{\text{d}{E}_{\text{g}}}{\text{d}t}=\frac{16}{3}{\pi }^{2}G{\overline{\rho }}^2{R}^4\dot{R}\text{, and }\frac{\text{d}{E}_{\text{t}}}{\text{d}t}=4\pi \overline{\rho }{c}_p\left(T_{\text{Imp}}-T_{\text{amb}}\right)R^2\dot{R}\hfill \end{array}$$

with the heat source term H_g arising from the impact energy buried within the depth ΔR, the fraction of the impactor kinetic energy converted to thermal energy f, the Stefan-Boltzmann constant σ_SB, the ambient nebula temperature T_amb, the thermal conductivity k, the radial coordinate r, the gravitational constant G, the average density $\overline{\rho }$, the rate of the radius change Ṙ, the average heat capacity ${\overline{c}}_p$ and the average temperature of the impactor T_Imp. Free parameters for the calculation of the accretional heating are f, ΔR and T_Imp. For the former two, same values as in Šrámek et al. (2012) are used, and for the latter one, a temperature of T_Imp = T_amb + 5 K is assumed in all calculations. Accretional heating was not included in previous models of Ceres in such a way.

Further energy sinks or sources are the consumption and the release of latent heat during melting of water ice and the crystallization of water, respectively. These processes are considered in a temperature interval of ten degrees around the melting temperature of T = 273 K in order to avoid numerical issues with too sharp a phase transition at 273 K. The latent heat contributes to the energy balance via Stefan number (e.g., Neumann et al. 2012).

3.3 Composition

Based on the conclusions from studies of crater relaxation regarding the composition of the subsurface (Bland et al. 2016), we use an H₂ O volume fraction of 0.25 and a dust volume fraction of 0.75 as standard values for the majority of calculations. In addition, considering Kuiper belt objects, we also use an H₂ O volume fraction of 0.5 (Desch et al. 2009; Bland & Travis 2017). From a bulk density of 2162 kg m⁻³ and water density of 918 kg m⁻³, a core density of 2580 kg m⁻³ follows for 25 vol.% H₂ O. Assuming some small average dustporosity of Ceres of 2%, the corresponding intrinsic density would be 2200 kg m⁻³ and the core density 2630 kg m⁻³. Both core densities fit into the interval spanned by aqueously altered CI and CM chondrites (≈ 2400−3000 kg m⁻³). Studies of spectral properties of the surface material (De Sanctis et al. 2015) indicate that its composition is rich in ammonium-bearing annite, antigorite (Mg-serpentine) and, possibly, magnetite (suggested for the “dark component”) as well as different amounts of chlorine salt (Na,Ca-perchlorates or chlorides). The abundances suggested by De Sanctis et al. (2015) reflect the composition of the material at specific locations on the surface. It is essential to note that the abundances retrieved from analyses of spectra are optical contributions and are, thus, at best areal fractions, not necessarily mass or volume fractions. For example, a small areal dark component fraction can considerably darken a mixture and further factors, such as grain size effects, can tone down spectral features non-linearly. Therefore, the range of reported magnetite abundances from spectral analyses would lead to too high a bulk density if assumed to convert 1:1 into a mass or volume fraction. These abundances should be changed considering the bulk composition of the rock phase on Ceres. Therefore, without attempting to model the chemical evolution of the rock, we use a composition that reproduces the above core density and is reasonable for Ceres to calculate compaction by creep by combining various species identified in spectra. For the calculation of the water-dust separation, we distinguish between solid ice and liquid water on one hand and “dust” on the other (the latter assumed as nearly a mixture of salts, e.g., chlorides, as well as antigorite, annite, and magnetite, neglecting less abundant species) within the composition outlined above. For the calculation of compaction, individual creep laws are used for the species (thereby, ice and salts are considered as one species before the melting of ice). For the volume fractions of the species, see Table 3 and Fig. 2. This is a rough representation of the composition and does not contain small species, for example, clathrates. It is, however, reasonably representative for Ceres in terms of the compaction and differentiation behavior and not far from typical CI and CM compositions in the terms of mineralogy and chemistry (Howard et al. 2011; King et al. 2015). Notably, representatives of both CI and CM chondrites contain roughly around 25 wt.% of iron, which is present mostly in the form of the oxides. A possible outcome of the calculations for the highly energetic thermal conditions in the planetesimals accreting close to the formation of the CAIs would be the formation of a metallic core made of this iron fraction (cf. Neumann et al. 2018b).

3.4 Compaction

We use a typical approach to model the compaction of planetesimals adopted in several studies of small bodies. For more details, we refer to Neumann et al. (2014b, where compaction of rocky planetesimals was explored) and Neumann et al. (2015), where the hypothesis of a porous interior for Ceres was tested with such a model. Here, the variation of the bulk pore space volume fraction (termed as “porosity” ϕ) is described by a time-dependent differential equation which establishes the relation between the strain rate $\dot{\epsilon }$ and the applied (effective) stress σ: $$\frac{\partial \left(1-\varphi \right)}{\partial t}=\dot{\epsilon }=\mathcal{B}{b}^{-m}{\sigma }^n\exp \left(-\frac{\mathcal{E}}{\mathcal{R}T}\right)$$(5)

with a pre-factor $\mathcal{B}$, the grain size b, the grain size exponent m, the effective stress σ, the stress exponent n, the activation energy $\mathcal{E}$ and the gas constant $\mathcal{R}$. For the mineralogical composition considered, the volume fractions v_i, i = 1, ..., 4 (where v₁ corresponds to ice plus salts or only to salts if the ice has melted, v₂ to antigorite, v₃ to annite, and v₄ to magnetite relative to the fraction of the solid phase), and the associated porosities ϕ_i. The creep laws for the single mineral phases provide the strain rates ${\dot{\epsilon }}_i$.

The ice and the salts are considered together for simplicity. Deformation studies ofice-perchlorate mixtures indicate more efficient creep than for pure ice. The grain boundary sliding-accommodated basal slip regime used here can be described by the following equation (Goldsby & Kohlstedt 1997, 2001; Lenferink et al. 2012): $${\dot{\epsilon }}_1=\frac{\partial \log (1-{\varphi }_1)}{\partial t}=6.31\times {10}^{-14}{\sigma }^{1.8}{b}^{-1.4}\exp \left(-\frac{49}{\mathcal{R}T}\right)\text{.}$$(6)

A Peierl’s law equation (Katayama & Karato 2008) was derived for the deformation of antigorite at pressures of ≤ 200 MPa by Amiguet et al. (2012): $${\dot{\epsilon }}_2=\frac{\partial \log (1-{\varphi }_2)}{\partial t}=3.98\times {10}^{-22}{\sigma }^2\exp \left(-\frac{27}{\mathcal{R}T}\left(1-\frac{\sigma }{2.7\times {10}^9}\right)\right)\text{.}$$(7)

An exponential law for biotite is used for annite (Kronenberg et al. 1990; Amiguet et al. 2012) (thereby the pre-factor accounts for the deformation between the easiest and hardest slip direction): $${\dot{\epsilon }}_3=\frac{\partial \log (1-{\varphi }_3)}{\partial t}=4.6\times {10}^{-8}\exp \left(4\times {10}^{-7}\sigma \right)\exp \left(-\frac{82}{\mathcal{R}T}\right)\text{.}$$(8)

This creep law describes a more effective deformation than other fits described in the literature (e.g., the exponential and power laws derivedby Shea & Kronenberg 1992). Finally, grain boundary diffusion creep law is used for magnetite (Till & Moskowitz 2013): $${\dot{\epsilon }}_4=\frac{\partial \log (1-{\varphi }_4)}{\partial t}=5.63\times {10}^{-25}\sigma b^{-3}\exp \left(-\frac{124}{\mathcal{R}T}\right)\text{.}$$(9)

In the above equations, the stress σ is in Pa, the grain size b in m, the activation energy $\mathcal{E}$ in kJ mol⁻¹, the gas constant $\mathcal{R}$ in kJ and the temperature T in K.

The bulk strain rate of the material with the mineralogical composition described above can be calculated as the average of the single strain rates: $$\dot{\epsilon }=\frac{\partial \log (1-\varphi )}{\partial t}=v_1{\dot{\epsilon }}_1+v_2{\dot{\epsilon }}_2+v_3{\dot{\epsilon }}_3+v_4{\dot{\epsilon }}_4$$(10)

and therefore the bulk porosity ϕ.

In the model, we use this approach on one hand for the calculation of the evolution of the dust porosity (i.e., we define ϕ_dust,0 = ϕ₀, ϕ_dust and ϕ_dust,i as well as ${\dot{\epsilon }}_{\text{dust}}$ and ${\dot{\epsilon }}_{\text{dust,}i}$ for the calculation of creep while ϕ_dust > 0), and on the other hand for the evolution of the melt porosity (i.e., we define the water fractions ${\varphi }_{\text{melt,}0}=v_{{\text{H}}_2\text{O,}0}$, ϕ_melt and ϕ_melt,i as well as ${\dot{\epsilon }}_{\text{melt}}$ and ${\dot{\epsilon }}_{\text{melt,}i}$ for the calculation of creep in the water-filled matrix where ϕ_dust = 0; we note that v_i change here because the ice is missing, see Table 3). Important free parameters are grain size b and the initial porosity ϕ₀. A grain size of b = 50 μm is based on the surface grain sizes derived by De Sanctis et al. (2015). An initial porosity of 40− 50% is used typically in compaction models of planetesimals and is based on the porosities of the random loose and random close packings as well as on the porosity of 40−50% expected after “cold pressing” acted on a granular material (e.g., Henke et al. 2012; Neumann et al. 2014b).

Fig. 2 Left: initial composition used for the simulation of compaction and differentiation of Ceres. Presence of abundant aqueously altered species is based on the calculations of the thermal evolution of small precursors that clearly predict aqueous alteration for the building blocks of Ceres (see Sect. 4.1 of this paper). Right: sketch of the initial three-phase micro-structure with rock grains, ice grains and voids (dust porosity).

3.5 Differentiation

Differentiation was included in several previous models. However, these models considered a composition that was dominated by ice and would lead to a fast separation of water and dust from each other upon the melting of ice. Typically, dust particles are assumed to settle instantaneously in a water layer in these calculations. Non-instantaneous water-rock separation was considered in the Stokes differentiation model by Wakita & Sekiya (2011) and in sophisticated differentiation models by Malamud & Prialnik (2013, 2015, 2016) for icy planetesimals. These studies did not consider Ceres, though. A continuous course of the differentiation was only due to the outward movement of the ice melting front. In the present study, the composition is dominated by dust and, more importantly, after the melting of ice, a three-phase system with empty pores, water and a dust matrix emerges. The separation of water and dust is guided by the compaction of ice and dust that leads to the reduction of the dust porosity, by the compaction of the dust matrix that contains liquid water but no dust pores, and by the settling of dust particles in a differentiated water layer. This way a continuous deformation of matrix that creates conditions for the percolation of water (dust porosity equal to zero, further deformation of the matrix and reduction of the melt porosity) is taken into account.

First, the closure of the dust porosity due to the elevated temperature proceeds from the center of the planetesimal where the temperature maximum occurs toward the surface following a monotonously decreasing temperature profile. In the volume in which the pores shrank to close the dust porosity, if the deformation proceeds further, the melt pores are reduced. The displaced water volume moves to the top boundary of the volume and forms a water layer while the dust fraction is displaced to the center, according to the melt porosity profile, and forms a proto-core. When the zero dust porosity front moves upwards across a layer between two time steps, the dust grains from this layer (which is located above the water layer) settle quasi-instantaneously at the top of the proto-core as a lose agglomerate with the initial melt porosity ϕ_melt,0. Upon heating, they undergo deformation, squeezing the water upwards into the ocean. See Fig. 3 for comparison.

Fig. 3 Left: differentiation sequence implemented in the model on a length scale of a few tens of microns. (a) Initial three phase system with ice particles (dark blue), rock particles (brown) and empty voids (white). (b) Partial closure of the dust pores, melting of water (light blue). (c) Complete closure of the dust pores, onset of the compaction of the melt pores. (d) Formation of a rocky core with a compacted inner part and a porous outer part. Right: locations of the local micro-structures a–d in a globally partially differentiated structure that can occur at some time during the evolution of a planetesimal. The central area consists of a rock core with a compacted inner and a porous outer layer. It is overlain by a water ocean below an undifferentiated crust. Liquid water is present in the lower crust, while the colder upper crust contains water ice.

3.6 Hydrothermal circulation, convection in the ocean, convection in the ice layer

The hydrothermal circulation in a porous water-filled layer, the liquid-state convection in the water ocean, and the solid-state convection in the solid ice layer are processes that take place under very different thermo-mechanical conditions. However, they can be parameterized in a similar way. In each of these cases, a Rayleigh number can be calculated. If its value is super-critical, then these processes contribute to the energy balance (the thermal conductivity is substituted by k_eff to simulate cooling). This contribution can be modeled using the effective thermal conductivity k_eff, that is either calculated from Darcy’s law (hydrothermal circulation, see, e.g., Neveu et al. 2015), or from the thermal conductivityand the Rayleigh number itself (convection of ice, see Desch et al. 2009; and convection of water in the ocean, see Neumann et al. 2014a). We use these modeling approaches and refer to the publications mentioned for detail.

3.7 Moment of inertia factor I/MR²

The factor I/MR² is calculated during each time step of a model run. The final value of I/MR² at t = 4.567 Gyr should correspond to the mean moment of inertia 0.37, based on the observations by Dawn (Park et al. 2016). Assuming a three-layer model with a rocky core, a water or ice mantle and a porous dust-rock crust, the moment of inertia factor is calculated as $$\frac{I}{{\text{MR}}^2}=\frac{2}{5}\left(\frac{{\rho }_{\text{cr}}}{\overline{\rho }}+\frac{{\rho }_{\text{co}}-{\rho }_{\text{ma}}}{\overline{\rho }}{\left(\frac{R_{\text{co}}}{R}\right)}^5+\frac{{\rho }_{\text{ma}}-{\rho }_{\text{cr}}}{\overline{\rho }}{\left(\frac{R_{\text{ma}}}{R}\right)}^5\right)$$(11)

where the indices co, ma and cr correspond to core, mantle and crust, respectively, and radii to the upper boundaries of the respective layers, while M is the mass, R is the average radius and $\overline{\rho }$ is the average density of the body. Before the formation of core and mantle, the object has only the “crustal” layer and, therefore, I/MR${}^2=(2/5){\rho }_{\text{cr}}/\overline{\rho }=0.4$ (assuming a homogeneous density profile). Since, prior to the differentiation, the density profile is not homogeneous due to the reduction of the dust pores, a small error occurs for the undifferentiated structure (i.e., I∕MR² ≈ 0.4). This error can be disregarded because all endmembers for Ceres are differentiated (see Sect. 4 of this paper).

3.8 Material properties

The material properties, such as density, thermal conductivity and specific heat capacity, vary with the porosity, the composition and the temperature. The local bulk density scales with the local porosity as ρ_bulk = (1 − ϕ_dust)ρ_intrinsic. The intrinsic density ρ_intrinsic is the density of the compacted material with ϕ_dust = 0 and is a volume fraction weighted arithmetic mean of the densities of water and dust: ${\rho }_{\text{intrinsic}}={\varphi }_{\text{melt}}{\rho }_{{\text{H}}_2\text{O}}+(1-{\varphi }_{\text{melt}}){\rho }_{\text{dust}}$. For ρ_dust see Table 2.

The bulk heat conductivity k_b varies with porosity as k_bulk = k_intrinsicf(ϕ), where $k_{\text{intrinsic}}={\varphi }_{\text{melt}}{k}_{{\text{H}}_2\text{O}}+(1-{\varphi }_{\text{melt}})k_{\text{dust}}$ and f(ϕ) varies exponentially between $\mathcal{O}(0.01)$ and 1 if the dust porosity varies between its high initial value of ϕ₀ and 0 (see Neumann et al. 2014b for details). Thereby, $k_{\text{dust}}={\displaystyle \prod k_i^{v_i}}$ with i ∈{1 = annite, 2 = salts, 3 = antigorite, 4 = magnetite, k_i is the thermal conductivity and v_i is the volume fraction of the dust species i relative to the dust fraction. For $k_{{\text{H}}_2\text{O}}$, and k_i and v_i see Table 3.

The heat capacity c_p is a mass fraction weighted arithmetic mean of the heat capacities of water and dust: $c_p=x_{\text{melt}}{c}_{p{\text{, H}}_2\text{O}}+(1-x_{\text{melt}})c_{p\text{, dust}}$. Thereby, c_{p, dust} = ∑ x_ic_p,i, with the mass fractionx_i and the heat capacity c_p,i of the species i. For $c_{p{\text{, H}}_2\text{O}}$ and c_p,i see Table 3; x_i are easily derived as products of v_i and the ratio of the species density and the intrinsic density of the dust.

4 Results

Before calculating the thermal evolution of Ceres, the thermal evolution and the internal structure of potential building blocks of Ceres are addressed (Sect. 4.1 of this paper). Subsequently, Ceres is considered, and its differentiation is calculated for different accretion scenarios (Table 4). The evolution of the shallow subsurface and of the water layer is analyzed. The conclusions regarding Ceres’ present state and the results for the aqueous alteration, water retention, structure, hydrothermal circulation, etc. for Ceres-like protoplanets as well as for smaller icy planetesimals and KBOs are summarized.

4.1 Ice-rich precursors of Ceres and other ice-rich small bodies across the solar system

Applying the model to icy planetesimals in general, the question is addressed of how far they can evolve in terms of heating and aqueous alteration, that is, for which parameter spaces can liquid water be produced in the course of thermal evolution, what is the maximum temperature reached, what is the duration of the super-liquidus phase, and is the water being retained.

Here the composition is assumed to be a mixture of nearly chondritic dust and ice, such that the average intrinsic density of the material is 2200 kg m⁻³. Here we use the material properties from Table 3 for simplicity. Although this assumption does not reflect an initial mixture of dry silicates and water ice, the thermal diffusivities differ only by a few percent and a potentially resulting temperature deviation is not expected to have any noticeable impact on the thermal evolution. The evolution of the temperature at the center and its maximum value is an indicator of the changes the material may undergo. If the temperature stays below the threshold temperature interval of ≈ 180−240 K for creep of ice, the whole body retains its porous structure. If the central temperatures exceed the melting point of water ice of 273 K, aqueous alteration starts and differentiation is theoretically possible. The maximum temperature T_c,max depends mainly on the dimensions of the body (the radius R is the measure of the removal of the heat from the interior through the surface by heat conduction) and on the amount of the short-lived radionuclide ²⁶Al (the accretion time t₀ + t_ac is the measure of how much of the initial inventory of ²⁶Al and radioactive material already decayed).

Table 4 shows the parameters adopted for every formation scenario investigatedfor the potential building blocks of Ceres. The accretion onset time used is t₀ = 0 Myr after CAIs, and the accretion duration t_ac is varied within the appropriate interval for the formation scenario assumed. The dust/ice ratio is set to 3:1 in every case, corresponding to the composition of Ceres’ subsurface and in addition to 1:1 for the Kuiper belt as an upper bound on the amount of ice.

Figure 4 shows the dependence of T_c,max on R and t₀ +t_ac for models of planetesimals accreted in the asteroid belt (upper left panel) and, for comparison, in the inter-planet disk (lower left panel), and in the Kuiper belt (upper right panel: dust:ice = 3∕1, lower left panel: dust:ice = 1∕1). Obviously the thermal evolution of equal bodies differs strongly between the scenarios with a small t₀ +t_ac ≤ 5 Myr and those with t₀ + t_ac ≥ 7 Myr.

4.1.1 Asteroid belt

Accretion of small planetesimals in the asteroid belt within 5 Myr can have a wide range of outcomes in terms of the peak temperature and internal structure (Fig. 4, upper left panel). The governing factor is the ²⁶Al inventory and the representative parameter isthe accretion time t_ac. The peak temperatures range from <300 K for a slow accretion of bodies with R ≤ 10 km to those surpassing the liquidus of silicates of ≈1800 K for an accretion within 1 Myr and an arbitrary radius ≥10 km. The high-temperature endmembers (t_ac ≲ 1−1.8 Myr) will clearly experience ice melting, aqueous alteration of silicates and even melting of silicates (although an explicit treatment of the latter two processes and of the dehydration is not implemented here). In particular, they will differentiate into a rocky core and a water ocean. Dehydration of the rocky core material should occur at silicate sub-solidus temperatures of T < 1425 K. Melting of Fe-FeS metal and of silicates would cause stratification of the core into an inner metallic and an outer silicate core. Such an energetic environment would result in an eventual loss of the ocean due to boiling and evaporation, leaving only anhydrous, mantle-stripped dense cores behind. Compaction in low-temperature endmembers can lead to the retention of a porosity that is different from zero throughout the interior. Such bodies are not differentiated at all, but the silicates will be pervasively altered (accretion within 1.5−5 Myr, R = 10−50 km). Finally, parameters in the intermediate range lead to only partially differentiated bodies with a small rocky core having a water-filled porous outer core, a water mantle, and a porous undifferentiated crust with an altered lower part and an upper part that has never experienced melting of ice. Such bodies accrete between 1 and 1.5 Myr for very small planetesimals and between ≈1.8 and up to 5 Myr for larger bodies. A variety of structures results from a substantial variation of the heating rate within the first few Myr after CAIs and an analogously strong variation of the cooling rate for bodies with different sizes (Fig. 4, structures 2–6 shown at the bottom and their locations in t_ac -R-plots). Layers that could form comprise differentiated metallic cores and dry rock mantles in extremely strongly heated bodies (Fig. 4). In a less energetic interior, dry rock inner cores, hydrated rock intermediate core layers, hydrated rock porous outer cores, H₂ O layers, undifferentiated porous hydrated rock lower crusts with interstitial H₂ O, and porous upper crusts composed of mixtures of dry rock and water ice can form. Finally, the alteration in the interior of extremely weakly heated cases is only marginal.

4.1.2 Inter-planet disk

Crucial parameters and conditions in the inter-planet disk are not too different from those in the asteroid belt. Although the initial and the surface temperatures differ by around 60 degrees (110 vs. 170 K), the accretion times are still small (≤5 Myr after CAIs). Therefore, ²⁶Al (and ⁶⁰Fe) still produces a substantial amount of energy to compensate for these less favorable initial and boundary conditions. The results regarding maximum temperatures change only marginally, shifting to slightly earlier formation times and larger radii. An analogous shift occurs for structures 2–6, while an additional structure 1 emerges for the coldest objects composed of dry rock and water ice that never produced any liquids (Fig. 4, top right panel).

The results for small icy bodies that evolved early in the inner part of the solar system imply, that, in general, a “fast” accretion of Ceres within ≥ 5 Myr would involve both dry (water-stripped, with a possibly differentiated core) and wet (aqueously altered or such with unmelted ice) planetesimals that contributed to its composition. In particular, for both accretion scenarios, hydrated minerals (e.g., Mg-serpentine antigorite) can be produced in abundance on small precursors of Ceres on a short time scale of ≪ 1 Myr after CAIs.

Fig. 4 Maximum central temperature Tc,max for 10 km ≤ R ≤ 200 km, t0 = 0 Myr and tac varying according to the accretion scenarios from Table 4. Upper left panel: asteroid belt, dust/ice = 3∕1. Middle left panel: inter-planet disk, dust/ice = 3∕1. Upper right panel: Kuiper belt, dust/ice = 3∕1. Middle right panel: Kuiper belt, dust/ice = 1∕1. Lower panel:final structures of initially icy planetesimals, numbered 1 through 6 according to the degree of differentiation. Not all of the structures are obtained for the accretion at different heliocentric distances (see Table 5 for an overview). Hydration of silicates is indicated by the dashed blue lines, starting at the water ice melting temperature of 273 K (serpentinization at 273 K takes ≈200 days according to Wakita & Sekiya (2011) and Jones & Brearley (2006); the temperature interval of 320− 520 K has been suggested by Allen & Seyfried 2004). We note that consideration of endothermic processes, such as melting of metal and silicates, could reduce the temperature and shift the isolines with T > 1260 K slightly to earlier accretion times and larger radii.

4.1.3 Kuiper belt, dust/ice = 3∕1

The most severe consequence for the accretion of small bodies in the outer solar system is the potentially comparatively long span of accretion times. While accretion of proto-Kuiper belt objects within 10 Myr after CAIs is currently considered as more likely (e.g., Morbidelli & Nesvorny 2019; Nesvorny et al. 2019), we present calculations also for extended accretion times of up to 100 Myr suggested previously for Ceres and its precursors (Castillo-Rogez et al. 2013; McKinnon 2012). Considering t_ac = 10−100 Myr after CAIs, ²⁶Al will be active only during the initial growth stage during which only a minor fraction of the final mass accretes to a small body. A simple calculation shows that after 5 Myr after CAIs, the amount of ²⁶Al decreases by three orders of magnitude (relative to its initial abundance at the time of CAI formation) and by ≈ 9 Myr after CAIs the energy produced by ²⁶Al falls below the cumulative energy of the long-lived radionuclides. For ⁶⁰Fe this happens by ≈4 Myr after CAIs. The short-lived radionuclides have, therefore, a negligible contribution to the overall heat budget. Missing this energetic peak, lower maximum temperatures are achieved, compared to the bodies in the asteroid belt and in the inter-planet disk. The thermal evolution is governed by the long-lived radionuclides that provide a nearly constant heating rate relative to short-lived ones but also produce considerably less energy of decay mainly due to their long half-lives. Although the duration of the accretion is a free parameter within the interval considered, it has a marginal effect for a planetesimal with any given size (Fig. 4). The isolines for T_c,max of Kuiper Belt planetesimals shown on the lower left plot in Fig. 4 have very small positive slopes, indicating that the exact accretion timing is only of secondary relevance. Independently of t_ac, only in largest planetesimals with R = 170−200 km partial separation of water and rock can be achieved. A radius of >70 km is necessary to hydrate dry silicate minerals. In particular, all objects within the size range considered retain water (both water of hydration and ice that never melted) and do not nearly reach the melting temperature of metal or silicates.

4.1.4 Kuiper belt, dust/ice = 1∕1

The most extreme case in the terms of heat production considered with a dust/ice ratio of 1 : 1 implies a further reduction of the internal heating rate because the water molecules store no radionuclides. Another drop of the T_c,max values occurs (Fig. 4, lower right panel). More ice-rich planetesimals do not differentiate (not even partially) into a rocky core and a water layer, as they remain too cold for the matrix deformation even if the ice melts. Only the largest among them (R > 110 km) produce enough heat for aqueous alteration, and, naturally, they keep their initial bulk water content. Some uncertainty regarding water-dust separation arises from a high water volume fraction of 50%, that could indicate that differentiation of water occurs by settling of dust in a liquid water layer instead of percolation of water in a matrix of partially sintered dust. In this case, differentiation would occur catastrophically in the layers with T >273 K. On the other hand, convective stirring of a water-dust mush could prevent phase separation for a while (or indefinitely, see Bland & Travis 2017) until a potential re-freezing of the water. As in the case of dust/ice = 1∕1, no boiling and no loss of water from the interior occurs in the Kuiper belt objects with R ≤ 200 km investigated.

It should be noted that icy planetesimals considered above have initially an ice fraction of ≈ 50%. The differentiation would be initiated by settling, but it could change to percolation upon hydration given fast hydration (Wakita & Sekiya 2011). Therefore, the mechanism we used here does not fully apply to the entire course of differentiation of these bodiesand they could have somewhat stronger differentiated structures than in our calculations (this does not change the general occurrence of structures shown in Fig. 4, though). It applies, however, fully to their outer core regions (i.e., regions where dust has been already concentrated such that the dust/ice ratio is locally enhanced) and to Ceres itself which accretes hydrated planetesimals.

4.1.5 Summary for precursors of Ceres

Icy planetesimals considered above have radii of ≤200 km, and they cover the size range reasonable for the precursors of Ceres. An icy object with a radius of 100 km represents approximately 1% of Ceres’ mass, while an object with a radius of 200 km represents ≈8%. Larger building blocks are unlikely for the dwarf planet considered.

A simple calculation based on the accretion law R(t) shows that for the accretion duration of t_ac = 5 Myr, the material accreting onto Ceres within the last 1000 yr of t_ac corresponds to a planetesimalwith R = 80 km, assuming an equal density. For the last 10⁴ and 10⁵ yr of t_ac, this is R = 171 km and R = 343 km, respectively. For t_ac = 50 Myr, the corresponding radii are 37, 80 and 171 km. These are hypothetical, because it is more reasonable that the accreting mass per time unit of 10³ −10⁵ yr is delivered by several building blocks instead of a single one. Therefore, planetesimals with R = 10−200 km shown in Fig. 4 cover the likely size distribution range. Assuming some formation scenario and a corresponding dust/ice ratio, the calculations presented above can estimate the degree of alteration, the extent of differentiation, and the possibility of water retention in Ceres’ precursors. In particular, they provide insight into the water content of these building blocks – accretion from “wet” planetesimals only, or from a combination of “wet” and “dry” objects. Table 5 provides an overview of these results.

The porosity of the material that covers the surface of Ceres during accretion is another uncertainty. Models presented above calculate the average porosity of building blocks during their evolution and provide their final average porosity. Its value is determined by the efficiency of compaction and varies between < 1% if a thin outer porous layer atop a compacted interior is present and >20% for an incomplete reduction of the porosity throughout the interior, and a high porosity in the outer part of the body. Impacts can generate additional porosity due to the debris of the back-falling material. The connection between the low-density regions and impact sites on the Moon imply that the shockwaves generated by on-going impacting resulted in micro-fractures and additional porosity of the lunar crust (Huang & Wieczorek 2012). Most pores in impact-generated chondrites appear to result from shock-related microfractures as well (e.g., Consolmagno et al. 2008). The average porosity of the impactors is, therefore, a lower bound. Although it would be more reasonable to adopt a somewhat smaller porosity for t_ac ≤ 3 Myr and a higher porosity for a longer accretion duration, most calculations for the evolution of Ceres are carried out with an initial porosity of ϕ₀ = 0.5 in the present study and only some with ϕ₀ = 0.4.

As shown further above, aqueous alteration is highly probable for a wide range of R and t_ac. If dry material is also accreted (either such material that was never hydrated or such that formed in hot planetesimals upon dehydration), the early thermal evolution of Ceres will deviate somewhat due to the perturbations in the energy balance caused by different material properties of dry silicates, but only if dry material is accreted in a substantial amount. However, as soon as hydration of this material is achieved on Ceres, further evolution of Ceres that is much more relevant on the R and t scale will proceed as it does in the calculation presented in the following.

4.2 Thermal evolution and differentiation of Ceres

In a broader context of the accretion of planets in the early solar system, small icy planetesimals as those considered above should have contributed not only to Ceres, but also to the water budget of the Earth, and so should have Ceres-like icy protoplanets themselves. Considering Ceres in the current section not only produces a potential evolution path for the physical structure of the dwarf planet and for the chemical properties of its material, but also provides important clues regarding the retention of water in icy planetesimals and connects to the overarching question of water delivery to Earth and other rocky planets.

In the following, the thermo-chemical and geophysical evolution of Ceres is demonstrated for selected accretion times appropriate for the formation scenarios at different heliocentric distances with the associated initial and boundary conditions. Each calculation covers the evolution from the accretion onset contemporaneously with the CAIs, over a varying accretion duration between 2 and 100 Myr, until present (t = 4.567 Gyr). The heating associated with the impacts by pre-heated ice-rock planetesimals that were hydrated in the course of their evolution (see previous chapter) is included, but is secondary compared to the internal heating by ²⁶Al and other radiogenic nuclides. Processes, such as dehydration, melting of metal and silicates, and formation of a metallic core and a dry or hydrated silicate mantle are addressed based on the temperature evolution, since they are not implemented in the numerical model.

For the formation at small to middle heliocentric distances, that is, in the asteroid belt and in the inter-planet disk, three cases with t_ac = 2 Myr (Fig. 5, top row), t_ac = 4 Myr (Fig. 5, middle rows) and t_ac = 7 Myr (Fig. 5, bottom row) are shown in the following. Roughly 13% of the initial ²⁶Al inventory remains at the moment when Ceres has grown to its final size in the first case, less than 3% in the second, and less than 0.2% in the third. But even in the last case, Ceres begins with a heat production rate that is one order of magnitude higher than that resulting from the decay of the long-lived radionuclides. Unequal initial ²⁶Al inventory entails a faster (small t_ac) or slower (large t_ac) warming of the interior. This dynamics is inherited by all processes that are activated or influenced by the temperature.

For an estimate, using volume fractions for the composition shown in Fig. 2 and assuming subsequent and not overlapping compaction, the porosity can be reduced in four steps as the temperature increases. The ϕ values that correspond to these steps depend on ϕ₀. For ϕ₀ = 0.48, these are ϕ = 0.38, 0.29, 0.12, and 0; and for ϕ₀ = 0.39, the values of ϕ₀ = 0.31, 0.23, 0.09, and 0. A similar calculation can be done for the stepwise development of the melt porosity. In reality, different species will compact more or less simultaneously with different compaction rates, and the porosity will be reduced in a temperature range with varying lower and upper boundaries for different depths due to pressure variations.

In general, the weakest species are ice and salts, and they facilitate the initial stage of the compaction. Relatively fast local porosity is reduced at T ≈ 180−220 K to ≈ 0.39. Subsequently, creep of antigorite reduces the porosity to ≈0.25 at T ≈240−270 K. After some stagnation, creep of annite (bringing the local porosity to ≈0.09) and then ofmagnetite facilitate closure of the pores at the temperatures between 480 K and up to those close to the critical point of water of 647 K, if the temperature increases faster than the compaction occurs (otherwise slower creep rate at lower temperatures can facilitate the pore closure). These values vary to some extent with the pressure and with the duration of heating at some specific temperature. Globally, compaction follows the maximum values of temperature and pressure, starting, therefore, at the center of Ceres and proceeding from the inside out. At some point, it comes to a halt, producing a completely compacted interior below an only partially compacted crust that remained cold and, therefore, porous.

4.2.1 Formation and evolution in the asteroid belt

The initial reduction of the porosity to small values throughout the deeper interior is facilitated by the early temperature spike at the center due to the decay of ²⁶Al and by a slow temperature increase during the accretion. The following thermal evolution is similar for different accretion durations, but absolute temperature values and the timing of temperature changes are strongly influenced by the accretion end time. The central area cools until accretion is finished, and beyond that, until the overall warming interior facilitates a positive trend. Energy input in the subsurface by accretional heating can produce a local temperature maximum (if the interior is cold) that can be inherited for a time of up to several hundreds of millions of years, depending on the accretion end time. The interior becomes continuously warmer, reaching the maximum temperature at the center between 1 and 2 Gyr after CAIs. This is followed by a continuous cooling phase until the present day (Fig. 5).

After the completion of accretion, the interior is still porous. If non-negligible amounts of ammonia are contained in the ice, first liquids are produced during the accretion since the initial temperature of 170 K and the melting point of ammonia dihydrate of 176 K are separated by 6 degrees, and such a temperature gap is overcome easily. This liquid occupies a volume fraction of approximately the same value as its volume fraction in the form of ice before melting. Since dust pores are still present, the pore water pressurization needed for the differentiation of water by upward percolation is lacking. Liquid water from the melting of pure water ice can react with the rock when both are in contact with each other (this can occur, in principle, prior to the end of accretion, but only locally at the center for a late accretion). The hydration time scale decreases with an increasing elevation of the temperature above 273 K. The time scale for the serpentinization of olivine of ≈200 yr at T = 273 K (and even faster at higher T; see Wakita & Sekiya 2011) is very short compared to the compaction time scale.

Closure of dust pores from the center upwards at T > 500 K (an effect of prolonged heating at relatively low temperatures) marks the onset of the differentiation. It passes continuously into the reduction of the melt pores. A downward compaction of the dense rocky matrix due to creep creates a pressurization of the interstitial water that is squeezed upwards. The growth of a core and of an ocean that is injected atop the core by the percolating water is controlled by the time scale of the compaction and occurs mainly on the time scale of 100 millions of years. The growth of the last few percent of the core radius can take one billion years or more because compaction in the shallower and colder interior is sluggish and comes slowly to a halt. The core develops a non-flat melt porosity profile with porosity of 0% at the center (completely compacted rock) and a melt porosity of nearly 100% at the core-ocean boundary (COB) corresponding to the ocean. The porosity in the outer core decreases continuously to either zero (early accretion with t_ac ≤ 2 Myr) or to an exponentially shaped branch from 0 to 100% in an ≈ 1.5 km thick layer (late accretion with t_ac ≥ 3 Myr).

When compaction of the subsurface stops, the upper boundary of the ocean reaches its maximum radial position. This happens approximately when cooling of the subsurface advances and the water crystallization front starts moving through the crust towards the ocean. The overall duration of the differentiation is < 2 Myr to < 5 Myr (t_ac≤ 2 Myr), ≈ 0.3 to ≈ 1.3 Gyr (t_ac= 4 Myr) and ≈0.54 to ≈ 1.3 Gyr (t_ac= 7 Myr).

The differentiated structure established comprises a rocky core (out to a radial distance of r = 420, 385 and 382 km, respectively) with a porous outer core for t_ac ≥ 3 Myr, an H₂O mantle (out to r = 462, 424 and 420 km, respectively) and an undifferentiated porous crust (5, 48 and 52 km thick) that is composed of partially sintered and hydrated rock grains, trapped liquid water and empty pores in its lower part and rock, water ice, and pores in its upper part. The water-rock separation concentrates radionuclides in the core, correlating the timing of the differentiation with the efficiency of the radiogenic energy production and with the temperature structure. An initially flat temperature profile developed a maximum in the core (Fig. 5). Two cooling processes – hydrothermal circulation in the outer core and liquid-state convection in the water ocean – can take place in parallel. The former is mainly active during the growth of the core when the melt porosity and, therefore, the permeability are high. The latter is active during the entire life time of the ocean because necessary conditions (Rayleigh number of > 3000 and a super-adiabatic heat flow at the upper boundary of the ocean) are satisfied and levels the temperature profile in the ocean. A layer of water ice forms when the water crystallization front traverses the upper ocean boundary. By this time, hydrothermal circulation has long ceased, but a thin porous outer core is still present as a structural layer. At a supercritical Rayleigh number, solid-state convection of ice is initiated and provides further cooling, establishing together with the shrinking ocean a two-layered convection system. By t = 4.567 Gyr after CAIs, the ice layer encompasses most of the H₂O layer. A basal water ocean is retained for a late accretion of ≥6 Myr, encompassing the lower half of the H₂O layer (i.e., ≈25 km) for t_ac = 7 Myr. While the core is still growing, its inner part dehydrates in a temperature interval of 770−870 K. For t_ac = 2 Myr, the entire core is dehydrated because most of it is at T > 870 K for several tens of Myr to >2 Gyr, while the outer 20 km are at T > 770 K for >2 Myr. For t_ac = 4 Myr and 7 Myr, the inner 180 km and 140 km, resp., are at T > 870 K for around half Gyr (dehydrating, therefore, completely) while the outer 125 and 142 km, resp., stay below T = 770 K (remaining hydrated). A high-temperature phase at T > 1262 K in the innermost 110 km occurs for t_ac = 2 Myr and implies production of metallic melt while the time scale of ≈0.6 Gyr of this phase implies migration of this melt to the center and formation of a small metallic core below a thin mantle made of dry silicates.

Fig. 5 Evolution of the temperature (left column) and density (right column) assuming accretion of Ceres in situ (first, second and fourth rows) and accretion and subsequent evolution in the inter-planet disk (third row). The accretion time is tac = 2 Myr (top row), tac = 4 Myr (second and third rows) and tac = 7 Myr (bottom row). Due to the sluggish compaction of the rock, the water separates very slowly on a time scale of hundred millions of years. The differentiation is not complete because a porous undifferentiated crust is retained. The water layer can be frozen at present (accretion in AB, tac < 7 Myr) or remain liquid at the base (tac ≥ 7 Myr). Free H2O is present inthe ocean or ice layer, while in the undifferentiated crust and in the porous outer core, interaction between the liquid water and the rock results in aqueous alteration if the precursor material contained dry silicates. The interior is colder for the inter-planet disk scenario. A warmer present-day interior is obtained for a longer accretion, which is an effect of a weaker compaction and, as a consequence, a slower cooling of the interior.

4.2.2 Formation and evolution in the inter-planet disk

The thermal and structural evolution, the timing of the processes, the final structure and the degree of alteration of the materialassuming the inter-planet disk accretion scenario are very similar to those for the in situ formation scenario for equal accretion times. The differences arise mainly due to the initial condition for the temperature T_init = 110 K (lower by 60 degrees than for the in situ formation scenario) and by the boundary condition for the surface temperature imposed by the ambient temperature of 110 K. The initial condition creates a gap of 60 degrees thatmust be bridged in the course of the thermal evolution, while the boundary condition keeps the surface and, therefore, the subsurface colder than in the previous accretion scenario. While the first constraint is only important for the early evolution, the second one is effective permanently and influences the final state to a substantial degree.

Within the typical evolution of the temperature that consists of an initial heating phase, a long-lasting high-temperature phase with one maximum and a final phase of slow cooling, formation of Ceres in the inter-planet disk results in a slower temperatureincrease prior to the maximum, a smaller maximum and a faster subsequent cooling than for an accretion in situ. The smallest differences to the first scenario regarding thermal conditions (minor), timing of the differentiation (a shift by ≈ 0.5 Myr) and the differentiated structure (negligible) are observed for a fast accretion with t_ac = 2 Myr. A shorter high-temperature phase of ≈0.3 Gyr within a radius of ≈50 km in this scenario still allows for the formation of a very small metallic core and an associated dry silicate mantle. These effects are owed to a still sufficiently large fraction of ²⁶Al by the end of accretion at 2 Myr after CAIs that overpowers the effects of small T_init and T_amb. Much more substantial differences arise if Ceres accretes slowly. The single stages, such as the melting of the ammonia dihydrate at 176 K, melting of pure water ice at 273 K, dehydration of hydrated silicates in the temperature interval of 770−870 K, compaction and differentiation are shifted in time by ≈0.1 Myr to ≈0.1 Gyr (see Fig. 5, second row for the formation in situ with t_ac = 4 Myr and third row for the formation in the inter-planet disk with t_ac = 4 Myr). The duration of the differentiation is cut short by this shift relative to the first scenario. The cores that are smaller by ≈ 7 km and the crusts that are thicker by ≈8 km are observed, while the H₂O mantle thicknesses are the same. The dehydrated inner cores are smaller by ≈ 20 km. The present-day state of the interior (in particular, of the shallow sub-surface) differs substantially in this scenario independently of the accretion duration (for example, the H₂O mantle is completely solid and it also does not convect). These results are presented in Sect. 4.2.6 of this paper.

Fig. 6 Evolution of the temperature (left column) and density (right column) assuming accretion of Ceres in the Kuiper belt for a dust/ice ratio of 3 : 1 (top and middle rows) and 1 : 1 (bottom row). The accretion time is tac = 10 Myr (top row) and tac = 100 Myr (middle and bottom rows). Thicker crust and smaller cores are obtained compared to models shown in Fig. 5. A warmer present-day interior is obtained for models with a longer accretion, which is an effect of weaker compaction and slower cooling.

4.2.3 Formation and evolution in the Kuiper belt

In this rather extreme scenario, the evolution of Ceres was considered for a small ice volume fraction of 0.25 (a dust/ice ratio of 3 : 1) as for other scenarios (Fig. 6, top and middle rows) and a high ice volume fraction of 0.5 (a dust/ice ratio of 1 : 1, Fig. 6, bottom row). We explore the effect of possible long accretion times by assuming t_ac = 10 Myr (Fig. 6, top row) and a potential although less likely maximum of t_ac = 100 Myr (Fig. 6, middle and bottom rows). The results vary little with t_ac, but differences arise while varying the dust/ice ratio. The outcome of the calculations is also considerably different from the previous scenarios. By 10 Myr after CAIs, the decay of ²⁶Al advanced so far that the heat is generated almost only by the long-lived radionuclides. Lacking an initial energy boost, the temperature increases very slowly. A period of ≈0.1−0.3 Gyr is necessary to climb from T_init = 50 K to T = 176 K and ≈ 0.2−0.4 Gyr to arrive at T = 273 K for t_ac = 10−100 Myr. After reaching a maximum temperature of 820−870 K (dust/ice = 3∕1) or of 706−752 K (dust/ice = 1∕1) at the center by t ≈ 2 Gyr or t ≈ 2.1−2.4 Gyr, the planetesimal cools continuously. The interior is, on average, substantially colder than for the in situ formation scenario, by ≈50−500 degrees at equal depths. The differentiation is still possible, however, it occurs considerably later (between ≈ 0.7 Gyr and > 2 Gyr after CAIs for dust/ice = 3∕1 and between1.1 Gyr and >2 Gyr for dust/ice = 1∕1) and is less pronounced. For a small ice fraction, the radial extension of the core and the mantle is r = 370−365 km and r = 408−402 km, resp., the thickness of the crust is 67− 74 km, and that of the porous outer core is ≈ 1−2 km for t_ac = 10 − 100 Myr.

For a high ice fraction of 0.5, the core is ≈90 km smaller (but the outer porous core layer varies between ≈2 and ≈3 km), the H₂ O layer is doubly thick as for dust/ice = 3∕1, and the crust thickness is almost doubled. An only partially dehydrated inner core of r ≈ 140−200 km forms for dust/ice = 3∕1, but there is no dehydrated inner core for dust/ice = 1∕1. A typical outcome for the H₂O mantle is a complete solidification by ≈0.7−1 Gyr (dust/ice = 3∕1) and by ≈0.5−0.8 Gyr (dust/ice = 1∕1) before present. However, in a H₂O mantle obtained for dust/ice = 1∕1 that is twice as thick and located deeper within Ceres, the solid-state ice convection proceeds until present.

Fig. 7 Evolution of the temperature (left column) and density (right column) assuming accretion of Ceres in the inter-planet disk (top row) and in the Kuiper belt (bottom row) and subsequent migration starting at 10 Myr after CAIs and 600 Myr after CAIs, respectively, within 1 Myr into the outer asteroid belt. The accretion time is tac = 7 Myr (top row) and tac = 100 Myr (bottom row).

4.2.4 Migration scenarios: IPD to AB, KB to AB

Formation in the inter-planet disk (IPD) or in the Kuiper belt (KB) implies migration into the asteroid belt (AB) at some point. We considered such scenarios as well, assuming that Ceres moves from the inter-planet disk starting at 10 Myr after CAIs, and from the Kuiper belt starting at 600 Myr after CAIs within 1 Myr into itspresent orbit in the outer asteroid belt (in both cases after the end of accretion, and in the second due to a less likely albeit possible late dynamical instability of the outer solar system potentially related to the late heavy bombardment). This is implemented in the model simply by increasing the ambient temperature within 1 Myr from 110 to 170 K and 50 to 170 K, for the inter-planet disk or the Kuiper belt origin, respectively. Due to the colder initial conditions, Ceres retains more porosity, that is, a thicker crust than for the accretion in the asteroid belt, promoting a warmer interior (cf. Figs. 5 and 7). Relative to non-migration scenarios (Figs. 5, 6), through moving into a warmer orbit with an ambient temperatureof 170 K, the interior becomes warmer only in the subsequent evolution and produces essentially the same differentiated structures but with a warmer present-day interior. The inter-planet disk migration scenario produces an H₂ O layer that is almost completely liquid (shown for t_ac = 7 Myr in Figs. 7 and 9) and Kuiper belt migration scenario produces even liquid water in the lower part of the undifferentiated crust and a completely liquid H₂ O layer (shown for t_ac = 100 Myr, i.e., the “coldest” case for KB, in Figs. 7 and 9).

The accretion and evolution scenarios presented above show similar trends if the accretion duration t_ac is varied within a single scenario. A slow accretion results in a late and less pronounced differentiation with a smaller core and a thinner H₂ O layer than for a fast accretion. The undifferentiated porous crust is thicker, which, in turn, results in a slower cooling and a warmer present-day interior. A similar trend is also observed while moving the formation orbit of Ceres farther into the outer solar system (if AB scenarios are compared with IPD scenariosforan equal accretion duration, or even combined with a longer accretion when comparing AB and IPD scenarios with KB scenarios).

A partial differentiation with a thick undifferentiated crust is an ultimate outcome of the models where the H₂ O fraction is relatively small, as it is in our model. This represents a challenge in terms of crustal stability, which is discussed in the next section.

Fig. 8 Evolution of the critical viscosity for the occurrence of Rayleigh-Taylor instabilities (red) and evolution of the viscosity at the base of the crust (blue) for the cases shown in Figs. 5–7. Top row: models from left to right and top to bottom in order of their appearance in Figs. 5–7.

4.2.5 Crust stability

A dense crust (with a depth-dependent density varying between 1400 and 2100 kg m⁻³) atop a less dense water or ice layer (ρ ≈ 1000 kg m⁻³) is a potentially gravitationally unstable configuration. Such crust can develop further by subduction and overturn, mixing crustal material with the H₂O layer. Whether this happens depends critically on the viscosity of the crust. A stiff crust with a viscosity that is higher than a “critical” viscosity for the development of Rayleigh-Taylor instabilities will not overturn and will remain stable on a geological timescale (see, e.g., Desch et al. 2009).

The stability of Ceres’ crust was studied by Shoji & Kurita (2014) and by Formisano et al. (2016) for mixed ice-rock compositions by comparing the mixture viscosity with a critical viscosity for the development of Rayleigh-Taylor instabilities. Considering a composition that is dominated by phyllosilicates (see above) and calculating viscosities of each species involved, bulk viscosity can be calculated as their mean value (see, e.g., De Meer et al. 2002; Grocott et al. 2009; Dygert et al. 2016). The formulation of bulk properties of composite systems uses arithmetic, geometric or harmonic average values of each material property. For the viscosity, the volume fraction weighted arithmetic mean (Voigt or isostrain approximation) provides an upper bound, while the harmonic mean (Reuss or isostress approximation) provides a lower bound. The geometric mean (Tullis model) is also frequently used for calculating the average viscosity of multiphase aggregates. The Voigt model is appropriate for aggregates with a strong matrix and isolated inclusions of a weaker phase, while the Reuss model represents the rheology of an aggregate for which the strain is accommodated by interconnected layers of a weak phase (Handy 1990). Since we assume that the weaker phase (water ice or liquid water) is no more than 25% by volume, the Voigt model appears more appropriate here: $${\eta }_{\text{mix}}=v_{\text{ice}}{\eta }_{\text{ice}}+v_{\text{ant}}{\eta }_{\text{ant}}+v_{\text{ann}}{\eta }_{\text{ann}}+v_{\text{mag}}{\eta }_{\text{mag}}$$(12)

with the viscosity of ice η_ice calculated as shown above and the viscosities of antigorite, annite and magnetite calculated using the relation $\eta =\sigma /\dot{\epsilon }$ (cf. Eqs. (6)–(9)) between the viscosity, the differential stress and the strain rate. For example, values obtained for ice, antigorite, annite and magnetite vary within the temperature range of 50− 500 K between 10⁶² and 10⁻³, 10⁴² and 10¹⁷, 10⁹⁸ and 10²⁰, and 10¹⁴⁰ and 10²³ Pa s for a differential stress of 10 MPa.

Figure 8 shows the evolution of the viscosity at the base of the crust (the minimum of the crust viscosity with depth, i.e., a lower bound on the average crust viscosity) and of the critical viscosity for the stability of the crust for the cases presented in Figs. 5–7. The critical viscosity for the development of Rayleigh-Taylor instabilities at the base of the crust is calculated as in Shoji & Kurita (2014) (see also Rubin et al. 2014; Formisano et al. 2016). We can confirm crustal viscosities calculated by Shoji & Kurita (2014), since they match approximately those that we calculated using the harmonic mean. Since these authors assume a higher H₂ O fraction, the rheology of the crust will be dominated by ice, and the calculation of the mixture viscosity according to the Reuss model is appropriate in their case. For our calculations, using the Voigt model, the viscosity at the base of the crust is larger than the critical viscosity, except for the otherwise improbable case with a very fast accretion (t_ac < 3 Myr after CAIs). This implies that if Ceres retained a largely undifferentiated crust in the course of the water-rock separation in the interior, this crust could have remained stable as a separate layer until present.

Fig. 9 Present temperature profiles for the cases shown in Figs. 5–7 in the same order as in Fig. 8. Dotted horizontal lines label the extension of the pure H2 O layer, dashed vertical lines correspond to the melting point of ammonia dihydrate of 176 K and to that of pure H2O ice of 273 K.

4.2.6 Present conditions in shallow sub-surface

Thermal conditions at 4.567 Gyr after CAIs for the cases presented in Figs. 5–7 are compared in Fig. 9 with the melting window of 176−273 K for various salt solutions and with the location of the ocean layer. The thermal models show that pure liquid wateris confined to deep interior and that only for the migration scenarios a basal ocean in the pure H₂ O layer (migration from IPD) is retained, or the H₂O layer remains entirely liquid (migration from KB). The latter is also the only case with liquid H₂ O in the lower 20 km of the crust below a depth of ≈40 km. In any caseconsidered, presence of ammonia dihydrate melt at depths of ≲ 2 km appears likely due to the surface boundary condition T_surf = 170 km and the temperature profile.

Figure 10 shows temperature profiles in the top 50 km of the crust, comparing formation in situ with migration scenarios and resulting thermal conditions with the eutectic temperatures of various salt solutions (Hennings 2014), as well as with the depth range of craters where floor faculae were detected according to the shape model. This comparison indicates that ammonia dihydrate solutions could be mobile up to shallow depths of ≲ 2 km, and some salts, such as Ca(ClO₄)₂ and Mg(ClO₄)₂ up to depthsof ≲10 km, while water remains confined to a deeper interior, explaining the low relaxation of the crust. While the eutectic temperatures are given at a pressure of 1 atm, subsurface pressure increases to >100 atm. at 20 km depth, which may further reduce the melting temperature of the ice-salt mixture. The combination of the thermal models with the shape model indicates that cryovolcanism is a likely and expectable process on Ceres. Aqueous alteration of the surface material after its accretion on Ceres is unlikely due to cold surface conditions. Observations of altered material imply alteration on precursors that formed the surface and/or extrusion of altered material from the interior by cryovolcanism. Some craters on Ceres show post-impact resurfacing by deposition of extended plains with pits, multiple flows and widely dispersed deposits building a diffuse veneer on the preexisting surface. Spectral analysis of the material of these features can suggest a small amount of hydrated salts added to ammoniated phyllosilicates. The context of smooth plains and flow-like features as well as mountainous material with distinct spectral variations and surface temperature variations suggests an endogenic cryovolcanic rather than exogenic origin. However, the correlation with impacts point to an exogenic triggering of cryovolcanic processes (Buczkowski et al. 2016; De Sanctis et al. 2015).

It should be noted that the depths for liquid species shown in Fig. 10 could change if phase diagrams of antifreeze or their thermo-physical properties were accounted for. These aspects are not included in our model. However, we do not expect a dramatic change in the results since fractions of the species from Fig. 10 are likely very small. In addition, solidification of those species that have high eutectic temperatures (e.g., MgSO₄, NaCl) would release latent heat, therefore acting toward decreasing the depth below which liquids of other species could be present, or extending the time scale for the existence of these liquids at a given depth.

Fig. 10 Temperature profiles in the top 50 km of Ceres (blue, red, and pink) compared with the eutectic temperatures of various salt solutions and pure water ice (black) and with the depth of craters with floor faculae (e.g., Ikapati, Occator, Haulani, Dantu, Urvara) derived from the shape model. The thermal models show that ammonia dihydrate solutions can be mobile up to shallow depths of < 2 km, and some salt solutions (Ca(ClO4)2, Mg(ClO4)2) up to depths of <10 km while water remains confined to a deeper interior. This combination of the thermal models with the shape model indicates that cryovolcanism is a likely and expected process on Ceres. The eutectic temperatures are given at a pressure of 1 atm. We note that Ceres’ subsurface pressure of ≈100 atm. at 20 km depth may further reduce the melting temperature of various ice-salt mixtures.

4.3 Summary for Ceres and Ceres-like protoplanets

Bodies addressed above are approximately the size of Ceres and have accreted in the asteroid belt, or farther away from the sun and have stayed in their orbit, or have migrated into the asteroid belt. Small objects in the IPD, such as Jupiter trojans, could be remnants of a part of such a body that was destroyed by impacts (given that the estimated mass of Jupiter trojans is smaller that that of Ceres, e.g., Jewitt et al. 2000). By contrast, trans-Neptunian objects of similar size to Ceres are abundant, and our investigations cover their evolution as well.

Had a part of the material that accreted to such bodies not been hydrated in the first place, it could have hydrated early in the course of the differentiation that proceeds slowly and leads to a partially porous core with interstitial water overlain by a water ocean below a porous undifferentiated crust. While the core becomes transiently warm, short-term hydrothermal circulation in its outer part and long-term liquid-state convection in the ocean cool the interior. Later on, ice convection in the crystallized upper part of the ocean establishes a two-layer convection system along with the liquid layer. Chemical structure consists of a dehydrated inner core, a hydrated outer core, liquid and/or solid H₂ O layers, a hydrated lower crust and an upper crust made of dry rock and ice (unless it contains hydrated material brought by smaller precursors). Icy protoplanets in the IPD evolve only marginally differently, but would end up with a substantially colder interior at present. Objects considered in the KB heat up more slowly and to a smaller maximum temperature, and they have thinner differentiated layers and thicker undifferentiated crusts. Their cores are only partially dehydrated for a moderate ice fraction andnot dehydrated at all for equal dust and ice volume fractions. The early evolution and the differentiated structure for the migration scenarios are very similar to those of the bodies at their original orbits. After the insertion into the warmer asteroid belt, the evolution paths of bodies from IPD approach those that evolved in AB from the beginning. Bodies that accreted very late and migrated from KB develop, however, to a much warmer present-day interior than in all other scenarios (Figs. 9 and 10) having a stronger crustal thermal shielding. Thereby, partial differentiation with a thick undifferentiated crust occurs in all scenarios except for t_ac < 3 Myr. Calculationsof the crust viscosity show that icy bodies with a composition similar to that considered here could have retained such an undifferentiated crust until present. It is also expected that pure liquid water can be found mostly in the outer core and only for a few models in the H₂O layer, while ammonia dihydrate and some salt solutions can be liquid at shallow depths of down to 2 km on Ceres.

5 Summary and discussion

We investigated the development of the temperature and internal structure of small icy planetesimals and of Ceres using a physical thermal evolution and differentiation model. Consequences of formation at different heliocentric distances with the associated accretion duration and of migration scenarios were investigated. Our results agree with the formation of Ceres in the asteroid belt or inter-planet disk, if the accretion occurs not too early and with a late accretion in the Kuiper belt. Both evolved icy building blocks of Ceres and the dwarf planet itself differentiate only partially. In particular, the differentiation of Ceres occurs slowly on a time scale of $\mathcal{O}(100)$ Myr relative to the melting of ice. Thereby, aqueous alteration is possible on the building blocks of Ceres, regardless of their accretion location.

While models produce an undifferentiated crust, its stability cannot be taken for granted, as analyzed in Sect. 5.2 of this paper. Another extreme scenario would produce a completely differentiated body if the crust were assumed to have a low viscosity. An intermediate case is represented by the mixing of crust with the water mantle, resulting in a mixed H₂ O-rock mantle (Fig. 11, endmembers 1.1 and 3.1) that could experience differentiation processes at present or where differentiation could have been interrupted due to cooling in a geologically recent past. Evidence of recent surface activity detected by Dawn could reflect active small-scale differentiation processes on Ceres.

5.1 Icy precursors of Ceres

Small icy objects with radii of ≤200 km were shown to be able to produce aqueously altered material for a wide range of R and t_ac (Fig. 4). Bodies that accreted in the outer asteroid belt and in the inter-planet disk are at least aqueously altered, and some of them are partially differentiated. Small icy bodies that accreted in the Kuiper belt can be aqueously altered as well, if their sizes surpass a critical value, but only the largest of them can be partially differentiated. Those bodies that accreted very early within ≈1.8 Myr after CAIs lost most of their water, while all others retained their H₂ O fraction in the form of ice or aqueously altered minerals. Thus, aqueous alteration of Ceres’ material could have occurred on small precursors. In particular, Ceres could have accreted from both dry and wet planetesimals in the asteroid belt and in the inter-planet disk, but only from wet bodies (i.e., ice-rock or aqueously altered) in the Kuiper belt (Table 5). Recent studies suggesting that a fraction of Ceres’ accreted material may be CI- or CM-like (Marchi et al. 2019; McSween et al. 2018) support our result that alteration could have occurred on smaller precursors of the dwarf planet.

5.2 Ceres’ internal structure and crust evolution

The final global interior structure obtained consists of a differentiated rocky core with a thin porous outer layer filled with interstitial water or ice, a H₂O mantle (its state ranging from completely solidified, over a basal liquid layer, to a completely liquid ocean for different models), and an undifferentiated crust (a three-phase partially compacted layer) where H₂ O can still be liquid in some accretion scenarios, but liquid brines are likely to be present at shallow depths. It should be acknowledged that these melts were not explicitly computed by the numerical model. The thermal evolution of the differentiated rocky core indicates that it can consist either of a dehydrated dry rock inner layer and hydrated rock outer layer, or of a locally moderately dehydrated material.

Depending on the dynamics of the undifferentiated crust (Fig. 11), the evolution of the crust and mantle can, in principle, take several different paths, some of which are not covered by the present model. The most straightforward path is retention of a stiff crust with a high viscosity (cf. Figs. 8 and 11, scenario 2). While in the deeper area where no dust pores are present, differentiation by compaction occurs fast, in the area with dust porosity, no differentiation is possible since the permeability is low (cf. Bland et al. 2009), as pores are locally preserved due to small strain rates. The study by Ruesch et al. (2019) suggests, however, that convection was active in the crust in geologically recent past. This could contradict the inability of the crust to overturn and would imply a lower crust viscosity than in our calculations, for example due to a different composition. Should the viscosity at the base of the crust be sub-critical, such that Rayleigh-Taylor instabilities could occur globally, the opposite development would result in the subduction and overturn of the crust (Fig. 11, scenario 1). In this case, undifferentiated material would be mixed with the H₂ O layer and the insulating porous surface layer would be destroyed. Should the energy budget suffice to keep all the water liquid (its volume fraction now being >0.25), and, thus the permeability high, a second differentiation stage would occur (scenario 1.2). The rocky core would grow at the expense of the formerly crustal rock, while the ocean would stretch all the way to the surface. If, however, only a small fraction of the water is liquid, then the permeability would remain relatively low and the second differentiation stage would not result in a complete phase separation, but rather in a more or less homogeneous mixing of crustal and mantle material (scenario 1.1). Furthermore, crust subduction can be facilitated by strong and numerous impacts (scenario 3). In this case, the two paths described above are conceivable (endmembers 3.1 and 3.2).

Following scenario 1.1 or 3.1, prevention of a complete water-rock separation and, simultaneously, sufficient energy production along with the inhibition of cooling such that liquids would be still present on Ceres represent a challenge. Mixing the crust with a purely H₂ O water ocean would result in a super-critical water fraction (e.g., ≈ 55% estimated for the migration scenario from IPD to AB with t_ac = 7 Myr if the crust would mix with the ocean). Differentiation in the Stokes regime could be prevented only if the liquid fraction were sub-critical (i.e., <50%) throughout the mixed layer. This could be achieved if this layer had a temperature of < 273 K, that is, if water were in solid state. Lacking further potent heat sources, no pure liquid H₂ O would be present today, only brines. Paths 1.2 and 3.2 lead to bodies that have no rock material remaining at the surface. This contradicts Dawn’s observations, unless there was a larger proto-Ceres that was stripped of the H₂ O layer and left its rocky core which represents the dwarf planet presently observed behind (Castillo-Rogez et al. 2020). In any of these scenarios, except in the one where the crust never subducts, the porous thermal shield is destroyed. However, this porosity is essential for keeping the sub-surface warm enough to have brines at a very shallow depth and pure liquid H₂ O in the sub-surface (unless some surface insulation can be introduced in some way in the other scenarios).

The structures that arise from the crust-ocean mixing (endmembers 1.1 and 3.1) raise the prospect that, given favorable thermal conditions in the interior, Ceres could experience some small-scale differentiation processes at present or that differentiation of the mixed layer has been interrupted due to cooling in a geologically recent past. Observations of recent surface activity act in favor of this prospect.

Fig. 11 Crust evolution scenarios. Differentiation models produce an undifferentiated crust that could remain stable following an extreme scenario 2. Another extreme scenario (endmembers 1.2 and 3.2) would produce a body with free water completely differentiated from rock if the crust were assumed to have a low viscosity. An intermediate case is represented by the mixing of crust with the water mantle resulting in a mixed H2 O-rock mantle (endmembers 1.1 and 3.1) that could experience some weak differentiation processes at present or where differentiation could have been interrupted due to cooling in a geologically recent past.

5.3 Solid-state convection in the crust

A first-order estimate of convective strength is the Rayleigh number Ra. For a subcritical Rayleigh number Ra < Ra_c≈ 10³, the heat in a system is transferred primarily by conduction, and for a supercritical Ra by convection (Schubert et al. 2001). As shown in Sect. 4.2.5 and Fig. 8, crustal convection is unlikely, if the viscosity is approximated by the Voigt model. In fact, a simple calculation yields Ra ≪Ra_c, and shows that a crust viscosity of ≲10¹⁸ Pa s is required to obtain a supercritical Ra value. Such low values correspond to a description of the crust composite system intermediate between the geometric mean and harmonic mean, counteracting crust stability and the warm present-day interior that is facilitated by a stable crustal thermal shield. Low viscosity values required for a supercriticalRayleigh number could be achieved for crust evolution scenarios 1.1 and 3.1 in Fig. 11 if the interior of the corresponding endmembers were as warm as that of an endmember with an intact crust. Since thermal evolution modeling for these scenarios goes beyond the scope of the present study, we refer to the convection analysis by Ruesch et al. (2019) that is applicable to these scenarios and results in supercritical Ra in a temperature range of ≈200−270 K for different convection modes and parameters appropriate here.

5.4 Consistency with observations

Crater relaxation models suggest that an ice-rich subsurface should relax large craters on short geologic timescales (Bland et al. 2016). Observations reveal large craters that are several kilometers deep, contradicting a pure ice subsurface and implying a subsurface material that is >1000 times more viscous than water ice. A composition accounting for crater morphologies consists of more than ≈ 70% of rock and salts, and less than ≈30% of ice (Bland et al. 2016). This agrees very well with an undifferentiated crust that has 25% ice in the calculations and is thicker by roughly a factor of ten than Ceres’ maximum crater depth (Figs. 5–7).

The density of the undifferentiated crust (Fig. 11, scenario 2) produced by the differentiation calculations is higher than the estimate by Ermakov et al. (2017). Since this estimate that was produced assuming a two-layer structure reflects Ceres’ crust density, there are several ways to reconcile our results with it. A lower crust density could be obtained with our model by adopting a higher initial porosity or if compaction were less efficient. While we cannot exclude a higher initial porosity or a more rigid initial composition, we considered these aspects as less likely and did not explore them here. Alternatively, scenarios 1.1 and 3.1, which involve mixing of the crust with the H₂ O layer after overturn and weak or no differentiation in the mixed layer afterwards, would produce a density profile consistent with the current estimates. This may require assuming some minor enrichment of rock downward and of less dense material upward, having potentially crustal porosity upon overturn. Such a scenario is, in general, worth exploring. Since our crust viscosity results argue for a stable crust for the composition considered, we do not explore crust subduction with a numerical model here. Another possibility to obtain a lower crust density is introducing clathrate hydrates as low-density phases. This would have the advantage that while an ocean-crust mixture would have 40 to 50% H₂ O, a water fraction of <0.4 consistent with crater morphologies would be possible by considering low-density clathrates. In addition, their role as a thermal insulator has been shown to be important for maintaining a long-lived subsurface ocean (Kamata et al. 2019). These compounds have been proposed to be present in Ceres’ shallow subsurface (e.g., Bland et al. 2016; Ermakov et al. 2017; Castillo-Rogez et al. 2018, 2019). While our work shows that invoking clathrates may not be necessary to explain the thermal, mechanical, and compositional state of the crust, they could resolve the mismatch between crustal densities calculated and those estimated by Ermakov et al. (2017). Since clathrates were not observed by Dawn, their abundance is uncertain, and experimental studies on their deformation behavior when mixed with the species assumed here are lacking, we did not include them in the model.

Our conclusion regarding the potential existence of a dry inner part of the core (suggested by previous models, e.g., Castillo-Rogez & McCord 2010; Neveu et al. 2015, and supported by our results) has not been a part of the model by Ermakov et al. (2017) who considered a two-layer structure. While not supported by observational evidence, the volume of a dry inner core volume should be small enough not to affect the moment of inertia factor significantly and the data do not preclude models with more than two layers compatible with the present modeling work.

Spectral observations of Ceres’ surface reveal spectral color variations correlating with morphology and topography, geologically young features, bright spots, and domes. A correlation of spectral variations with features of different geological ages (Jaumann et al. 2015) points to possible mixing with material from the subsurface, such as hydrated salts (De Sanctis et al. 2015, 2016; Ruesch et al. 2018). This agrees with the warm subsurface obtained in the models (Figs. 9 and 10), that indicate mobile brines and explain potentially cryovolcanic features observed.

Gravity measurements point to a moment of inertia factor I∕ MR ² of ≈ 0.37 (Park et al. 2016). In our calculations, I∕MR² varies between 0.36 and 0.38, bracketing these estimates (Table 6). It should be noted, however, that different interior models can result in the same I∕MR² value, such as for the two-layer structure suggested by Park et al. (2016) and three-layer structures produced by our models.

5.5 Cryovolcanism

Potentially cryovolcanic features on Ceres and mechanisms that led to their formation are in striking contrast to the cryovolcanic activity on Saturn’s moon Enceladus, where gas plumes caused by local extrusion of water from long linear depressions up to 130 km known as “Tiger Sripes” were revealed by the Cassini mission (see, e.g., Neumann & Kruse 2019, for a summary).

Ceres’ cryovolcanic features that are associated with craters, that is, floor faculae, have been suggested to originate from crystallization of brines that reached the surface from below and were dated as relatively recent (≈ 4 Myr) in the case of Occator (De Sanctis et al. 2016). The presence of brines at a depth of only a few km appears likely based on the comparison of present temperature profiles with the eutectics of various salt solutions shown in Fig. 10. Ascension ofcrustal material enriched in brines, either as part of a deep subsurface layer or accumulated in smaller local reservoirs within the crust at different depths, could have been promoted by tectonic weakening. These reservoirs could also have been exposed by impacts (Jaumann et al. 2015). In the case, however, that salts were not present as liquid solutions, they can be mobilized by impact heat (Hesse & Castillo-Rogez 2019). Another possible emplacement scenario is the formation of subsurface reservoirs enriched by a salt-bearing component and driven by density or temperatureinhomogeneities. A possible example for such a mechanism is the tectonic structure at the southern part of Occator, where a pattern of cracks indicates breaking of the plain material by upwelling and the partly polygonal arrangements suggest formation by dehydration (e.g., Buczkowski et al. 2016, 2018; Scully et al. 2018, 2019). Another explanation for this feature is refreezing of impact melt that could force material up and out, building a dome. However, the plains in Occator are more likely the result of impact-triggered cryovolcanic activity (Bowling et al. 2019).

Figure 12 sketches a mechanism that could have brought liquid brines to the surface. Impacts hitting subsurface reservoirs enriched in hydrated salts will mobilize these compounds, increase the fraction of brines by the production of additional melt as a result of impact-produced heating, and pour them out as flows. Ejected material covers the flanks of craters as flows and as fine dispersed erupted ash-like deposits of possibly salty composition (e.g., Buczkowski et al. 2016, 2018; Scully et al. 2018, 2019). Our thermal modeling indicates that thermal conditions can, at depths of only a few km below Ceres’ surface, favor the presence of brines and, thus, their mobility starting at this depth (Figs. 9, 10). This would also explain the buoyancy of ice and salt-enriched crustal reservoirs. Impacts may have cut into layers of salt-rich crustal reservoirs, triggering the mobility that formed the faculae and sustained long after the impact energy decayed, as indicated by the significantly younger ageof the plains and flow materials (Jaumann et al. 2015; Hesse & Castillo-Rogez 2019; Raymond et al. 2020).

Another cryovolcanic mechanism was suggested by Ruesch et al. (2019) as a source of the dome Ahuna Mons. This study constrained a specific form of cryovolcanism producing mud volcanoes that is different from cryovolcanism on other icy bodies. Convective processes in a deeper subsurface can cause formation of a regional uplift and of a cryochamber containing solid–liquid slurry (i.e., a fluid entrained with ≈30−45 vol.% of non-soluble particles) of brine and solid particles (Ruesch et al. 2019). As a result of a lower specific weight compared to the surrounding materials, the cryochamber moves upwards pushing the ice-rich crust, and conduits form at structural weak points. The slurry is extruded to the surface, solidifies in the cold of space, and piles up to form a dome. Mechanisms that initiate the mantle uplift are either the ascent of a mantle plume formed due to ongoing convection, or inflow of liquids into the crust due to active hydrothermal circulation (Ruesch et al. 2019). Here we suggest an alternative in the sense of our model results. While Ruesch et al. (2019) proposed slurry extrusion from the mantle, we reconcile such mechanism with our model, where a slurry reservoir could potentially form in the lower crust and the liquids could be extruded to the surface. While the crust that contains liquids is not completely water-filled on global scale, since some dust porosity is retained, lateral inhomogeneity in the terms of crustal viscosity and mobility could be present, such that locally voids in the crust could be fully filled with liquids (potentially involving some liquids inflow from the ocean) that could ascend to the surface. With a density of up to 1300 kg m⁻³, such liquids (termed as slurry, brines or mud in Ruesch et al. 2019) are less dense than any porous rocky part of the crust (cf. Figs. 5–7). After the extrusion of the low-density liquids, the interior would be locally dominated by the rock component producing a mascon. It must be noted that an ongoing convection in the crust contradicts the viscosity shown in Fig. 8. In principle, our calculations could overestimate the crust viscosity and do not consider lateral inhomogeneity of crust properties, that would potentially allow for a sluggish convection in the crust (see also Sect. 5.3 of this paper).

Extrusion of slurry from a mud-bearing interior (where mud is a viscous mixture of water and suspended dust particles, see Bland & Travis 2017; Ruesch et al. 2019) on partially or weakly differentiated Ceres suggested is strikingly different from water-dominated cryovolcanism on icy satellites, such as Enceladus, that did not retain an undifferentiated crust, and is directly linked to the differentiation degree of the surface. Slurry extrusion does not necessarily require impacts and, in fact, Ahuna Mons is not associated with an impact crater. The scenario presented in Ruesch et al. (2019) requires, however, ongoing or recent solid-state convection in a mud-bearing interior.

Fig. 12 Scenario for the formation of cryovolcanic features associated with craters: spectrally bluish material originates from subsurface, cryoplutonic sources released or triggered by impact. Upwelling of salt rich hydrated (bluish) material forms subsurface reservoirs. Impacts mobilize such reservoirs, producing lava-like flows ponding in the crater. Ejected salt-rich hydrated material covers the flanks of craters as flows and as fine dispersed ejected or erupted salt deposits.

5.6 Comparison with previous models

Multiple global thermal evolution and differentiation models and only a few studies of the crust stability of Ceres were published in the past.

While we did not include all aspects of previous models, our model setup has many similarities to them, simply by the nature of a thermalevolution and differentiation code. However, we also considered key processes that were not modeled in previous studies. In general, our model complements these studies (e.g., McCord & Sotin 2005; Castillo-Rogez & McCord 2010; Castillo-Rogez 2011; Neveu et al. 2015; Neveu & Desch 2015; Bland & Travis 2017; Castillo-Rogez et al. 2019)by considering continuous accretion (included previously in Neumann et al. 2015, 2016a) and by including a complementary differentiation mechanism by water percolation which was not investigated in any previous model except in our own preliminary studies (Neumann et al. 2016a,b, 2017a,b). In particular, this mechanism utilized both dust and melt porosities and compaction for a mineralogical composition derived from spectral observations of Ceres. These model aspects produced a partially differentiated interior structure with a retained undifferentiated crust that was not obtained by previous workers.

Our results of the inevitable water melting and water-rock differentiation in the deeper interior regardless of the accretion time are consistent with predictions of previous models. We added consideration of the precursor thermal evolution and differentiation, of the accretion of Ceres in the inter-planet disk and in the Kuiper belt, as well as migration scenarios, that produced, in general, the same type of a partially differentiated structure. In particular, we demonstrated that small precursors of Ceres could have been hydrated regardless of their accretion location. Thereby, one of the new results of our model is the porosity of the crust, which has been considered as negligible in the past (Castillo-Rogez 2011) or constrained mainly to the outer ≈ 10 km (Neumann et al. 2015).

Further general results that we confirm is the wide-spread water activity during the evolution of Ceres, the critical role for the final structure and liquid water content played by the amount of ²⁶Al (i.e., the accretion time), potential present-day presence of liquid water in the sub-surface for a narrow parameter space, likely presence of brines, possibility of cryovolcanism, and a potentially dehydrated inner core.

As expected, our model produces a quite different general outcome than the “mudball” model by Bland & Travis (2017) that is based on the opposite assumption of no lithification of the material (i.e., no pore reduction) and does not result in partially or fully differentiated layers, but in a thoroughly mixed “muddy” ice-dust interior.

The current understanding of Ceres’ interior appears to converge on a brine layer at depth (Nathues et al. 2017; Castillo-Rogez et al. 2020; Quick et al. 2019). Assuming a structure with a volatile-dominated crust and rocky mantle with remnant brines, Castillo-Rogez et al. (2020) showed that Ceres should preserve a warm crust with a temperature of > 220 K at the base until present if the crust is rich in clathrate hydrates, allowing for a brine layer. The structure inferred is similar to that produced in our calculations, but the density of the brine layer is higher than in the upper crust, while in our models the situation is opposite. As pointed out further above, clathrates could resolve the mismatch between crustal densities calculated here and those estimated by Ermakov et al. (2017). An additional factor is a potential suspension of micron-sized fines in the H₂ O layer suggested by Neveu & Desch (2015), that would on one hand increase its density and keep this layer warm due to radioactive heating by the radionuclides contained in the fines on the other.

Ceres’ crust stability was studied by Shoji & Kurita (2014) and by Formisano et al. (2016) for mixed ice-rock compositions by comparing the mixture viscosity with a critical viscosity for the development of Rayleigh-Taylor instabilities. Shoji & Kurita (2014) assumed a higher H₂O fraction than in our study and a crust rheology dominated by ice. Therefore, they used the Reuss model for viscosity calculations and obtained considerably lower viscosity values than in our study. Formisano et al. (2016) adopted a mixed composition with 7/3 to 9/1 olivine/ice ratios and reported a crust viscosity of ≈ 10²¹ Pa s using Eqs. (13)–(15). Testing the viscosities of an ice, rock and ice-rock mixture with $$\begin{array}{lll}\hfill & \hfill & {\eta }_{\text{ice}}={10}^{14}\exp \left(25\left(\frac{273}{T}-1\right)\right)\hfill \\ \hfill & \hfill & {\eta }_{\text{rock}}={10}^{21}\exp \left(\frac{525}{0.00831\times 1425}\left(\frac{1425}{T}-1\right)\right)\hfill \\ \hfill & \hfill & {\eta }_{\text{mix}}={\eta }_{\text{rock}}^{0.7}{\eta }_{\text{ice}}^{0.3},\hfill \end{array}$$

we obtain values of 9 × 10¹⁷, 7.5 × 10¹³⁸ and 4 × 10¹⁰² Pa s, respectively for T = 200 K (the crust temperature calculated in Formisano et al. 2016). Under such conditions the crust is not expected to develop Rayleigh-Taylor instabilities.If crustal movement would occur at all, it would be potentially due to other processes such as brittle failure. For the temperatures of 273 and 500 K calculated in the undifferentiated crust in this paper, crust viscosities of 1.6 × 10¹⁰² and 4 × 10⁵⁶ Pa s are obtained assuming a water viscosity of 0.0016 Pa s at T > 273 K. This outcome was to be expected, since the upper equation for η_rock corresponds to olivine-dominated rock that has a viscosity of 10²¹ Pa s at 1425 K and higher viscosities at lower temperatures. Calculations of the critical viscosity for the stability of Ceres’ crust starting from differentiation to the present shown in Sect. 4.2.5 of this paper provide η_crit < 10²⁴ Pa s, indicating that an olivine-dominated crust such as that assumed by Formisano et al. (2016) would remain stable on a time scale of the age of the solar system. However, olivine rheology is hardly appropriate for the rock component of Ceres’ crust, as demonstrated by Dawn’s observations. Our calculations presented in Sect. 4.2.5 indicate crust stability also for a less viscous rock composition than olivine.

6 Conclusions

The presentation of a new differentiation model for icy bodies and its application to Ceres in this paper provides an explanation for the observation of a largely undifferentiated crust on Ceres. Our results indicate that Ceres’ material could have been aqueously altered on small precursors and that some of these precursors were dry if the dwarf planet accreted within Pluto’s orbit. The differentiation in the water percolation regime proceeds slowly on the time scale of several hundred million years and leads to a final structure of a partially differentiated body with a liquid layer that could still be present between the porous undifferentiated crust and the silicate core, which is consistent with Dawn observations. The evolution paths for Ceres that migrates from the inter-planet disk are close to those for the evolution in the asteroid belt, while migration from Kuiper belt leads to a considerably warmer present-day interior than in all other scenarios having a stronger thermal shielding by a thicker porous crust. Calculations of the crust viscosity indicate crust retention and presence of liquid brines at shallow depths and of liquid water in the deeper interior. This is consistent with the observations of geologically recent water activity and cryovolcanism on Ceres and further implies that the crust could still be evolving.

Acknowledgements

The initial work for this study was supported by the special research program 1833 “Building a Habitable Earth,” (grant number NE 2055/ 1-1), sponsored by the Deutsche Forschungsgemeinschaft (DFG). We thank an anonymous reviewer for constructive comments and helpful suggestions which led to a substantial improvement to the original manuscript.

References

All Tables

All Figures

Fig. 1 Overview of the main differentiation regimes that arise for different assumptions on the proportions of ice and rock (Neumann et al. 2016a). Left: composition is dominated by ice whose melting results in the formation of an ocean. The rock mantle grows from the settling of rock particles in the Stokes regime and their deposition at the bottom of the ocean. Right: composition is dominated by rock. Melting of ice creates a water network within a porous matrix. The ocean grows due to the percolation of water which is driven by the gravity-induced compaction of the matrix (e.g., McKenzie 1984). Differentiation in this regime requires a longer time scale. In general, planetesimals of equal size but with different ice fractions would differentiate on a different time scale and result in structures with different layer dimensions.

In the text

Fig. 2 Left: initial composition used for the simulation of compaction and differentiation of Ceres. Presence of abundant aqueously altered species is based on the calculations of the thermal evolution of small precursors that clearly predict aqueous alteration for the building blocks of Ceres (see Sect. 4.1 of this paper). Right: sketch of the initial three-phase micro-structure with rock grains, ice grains and voids (dust porosity).

In the text

Fig. 3 Left: differentiation sequence implemented in the model on a length scale of a few tens of microns. (a) Initial three phase system with ice particles (dark blue), rock particles (brown) and empty voids (white). (b) Partial closure of the dust pores, melting of water (light blue). (c) Complete closure of the dust pores, onset of the compaction of the melt pores. (d) Formation of a rocky core with a compacted inner part and a porous outer part. Right: locations of the local micro-structures a–d in a globally partially differentiated structure that can occur at some time during the evolution of a planetesimal. The central area consists of a rock core with a compacted inner and a porous outer layer. It is overlain by a water ocean below an undifferentiated crust. Liquid water is present in the lower crust, while the colder upper crust contains water ice.

In the text

Fig. 4 Maximum central temperature Tc,max for 10 km ≤ R ≤ 200 km, t0 = 0 Myr and tac varying according to the accretion scenarios from Table 4. Upper left panel: asteroid belt, dust/ice = 3∕1. Middle left panel: inter-planet disk, dust/ice = 3∕1. Upper right panel: Kuiper belt, dust/ice = 3∕1. Middle right panel: Kuiper belt, dust/ice = 1∕1. Lower panel:final structures of initially icy planetesimals, numbered 1 through 6 according to the degree of differentiation. Not all of the structures are obtained for the accretion at different heliocentric distances (see Table 5 for an overview). Hydration of silicates is indicated by the dashed blue lines, starting at the water ice melting temperature of 273 K (serpentinization at 273 K takes ≈200 days according to Wakita & Sekiya (2011) and Jones & Brearley (2006); the temperature interval of 320− 520 K has been suggested by Allen & Seyfried 2004). We note that consideration of endothermic processes, such as melting of metal and silicates, could reduce the temperature and shift the isolines with T > 1260 K slightly to earlier accretion times and larger radii.

In the text

Fig. 5 Evolution of the temperature (left column) and density (right column) assuming accretion of Ceres in situ (first, second and fourth rows) and accretion and subsequent evolution in the inter-planet disk (third row). The accretion time is tac = 2 Myr (top row), tac = 4 Myr (second and third rows) and tac = 7 Myr (bottom row). Due to the sluggish compaction of the rock, the water separates very slowly on a time scale of hundred millions of years. The differentiation is not complete because a porous undifferentiated crust is retained. The water layer can be frozen at present (accretion in AB, tac < 7 Myr) or remain liquid at the base (tac ≥ 7 Myr). Free H2O is present inthe ocean or ice layer, while in the undifferentiated crust and in the porous outer core, interaction between the liquid water and the rock results in aqueous alteration if the precursor material contained dry silicates. The interior is colder for the inter-planet disk scenario. A warmer present-day interior is obtained for a longer accretion, which is an effect of a weaker compaction and, as a consequence, a slower cooling of the interior.

In the text

Fig. 6 Evolution of the temperature (left column) and density (right column) assuming accretion of Ceres in the Kuiper belt for a dust/ice ratio of 3 : 1 (top and middle rows) and 1 : 1 (bottom row). The accretion time is tac = 10 Myr (top row) and tac = 100 Myr (middle and bottom rows). Thicker crust and smaller cores are obtained compared to models shown in Fig. 5. A warmer present-day interior is obtained for models with a longer accretion, which is an effect of weaker compaction and slower cooling.

In the text

	Fig. 7 Evolution of the temperature (left column) and density (right column) assuming accretion of Ceres in the inter-planet disk (top row) and in the Kuiper belt (bottom row) and subsequent migration starting at 10 Myr after CAIs and 600 Myr after CAIs, respectively, within 1 Myr into the outer asteroid belt. The accretion time is tac = 7 Myr (top row) and tac = 100 Myr (bottom row).
In the text

	Fig. 8 Evolution of the critical viscosity for the occurrence of Rayleigh-Taylor instabilities (red) and evolution of the viscosity at the base of the crust (blue) for the cases shown in Figs. 5–7. Top row: models from left to right and top to bottom in order of their appearance in Figs. 5–7.
In the text

	Fig. 9 Present temperature profiles for the cases shown in Figs. 5–7 in the same order as in Fig. 8. Dotted horizontal lines label the extension of the pure H2 O layer, dashed vertical lines correspond to the melting point of ammonia dihydrate of 176 K and to that of pure H2O ice of 273 K.
In the text

Fig. 10 Temperature profiles in the top 50 km of Ceres (blue, red, and pink) compared with the eutectic temperatures of various salt solutions and pure water ice (black) and with the depth of craters with floor faculae (e.g., Ikapati, Occator, Haulani, Dantu, Urvara) derived from the shape model. The thermal models show that ammonia dihydrate solutions can be mobile up to shallow depths of < 2 km, and some salt solutions (Ca(ClO4)2, Mg(ClO4)2) up to depths of <10 km while water remains confined to a deeper interior. This combination of the thermal models with the shape model indicates that cryovolcanism is a likely and expected process on Ceres. The eutectic temperatures are given at a pressure of 1 atm. We note that Ceres’ subsurface pressure of ≈100 atm. at 20 km depth may further reduce the melting temperature of various ice-salt mixtures.

In the text

Fig. 11 Crust evolution scenarios. Differentiation models produce an undifferentiated crust that could remain stable following an extreme scenario 2. Another extreme scenario (endmembers 1.2 and 3.2) would produce a body with free water completely differentiated from rock if the crust were assumed to have a low viscosity. An intermediate case is represented by the mixing of crust with the water mantle resulting in a mixed H2 O-rock mantle (endmembers 1.1 and 3.1) that could experience some weak differentiation processes at present or where differentiation could have been interrupted due to cooling in a geologically recent past.

In the text

Fig. 12 Scenario for the formation of cryovolcanic features associated with craters: spectrally bluish material originates from subsurface, cryoplutonic sources released or triggered by impact. Upwelling of salt rich hydrated (bluish) material forms subsurface reservoirs. Impacts mobilize such reservoirs, producing lava-like flows ponding in the crater. Ejected salt-rich hydrated material covers the flanks of craters as flows and as fine dispersed ejected or erupted salt deposits.

In the text