© The Authors 2024

1 Introduction

The origin of water on Earth is still a subject of heated debate that involves many unknowns (e.g. O’Brien et al. 2018). Water is an essential ingredient for building life (e.g. McKay 2014) and understanding its origin would be of paramount relevance. This knowledge could then be extrapolated to extra-solar systems and help estimate the uniqueness of Earth compared to other exo-planets. In the currently accepted paradigm, it is thought that the inner Solar System (e.g. the Earth’s feeding zone) was too warm to hold water ice and that the water had come from the outer parts of the Solar System. The amount of water on terrestrial planets is relatively small; for instance, the Earth containing only 0.023% water by mass¹, whereas the giants of the outer Solar System contain up to 40% water (e.g. Helled & Lunine 2014).

The most widespread theory to date is that the water may have been brought in by impacts from icy bodies, mainly asteroids and possibly comets in smaller quantities (e.g. Raymond et al. 2006; Chen et al. 2019). One key constraint supporting this scenario is that Earth’s D/H ratio is very similar to that of carbonaceous chondrites (Alexander et al. 2012), which are thought to have originated from C-type asteroids, suggesting that they were Earth’s main source of water. There are several different scenarios to explain how outer icy bodies could have been brought to the inner Solar System, but they all share some common features, which can be summarised as follows. Some icy asteroid-like bodies from the outer regions of the Solar System (4–10 au for the most part) were perturbed (either by scattering induced by long term resonances or by more violent instability), allowing them to populate the asteroid belt and, at the same time, cause impacts on terrestrial planets. As these bodies are icy, they can then supply water to the planets via these impacts (e.g. Raymond et al. 2006; Walsh et al. 2011; Raymond & Izidoro 2017a).

However, some alternative scenarios suggest that the Earth’s water could be local in origin and may have been inherited very early in the Earth building blocks from a material similar to enstatite chondrites, which are very dry (but can contain enough hydrogen to provide sufficient water) and very similar to the material that formed the Earth (Piani et al. 2020). This may, however, be in disagreement with the composition of the Earth’s mantle and core, where models suggest that water is added to the Earth mainly at the end of its formation (Wood et al. 2008; Rubie et al. 2015). It should also be noted that there is evidence that water could have been present on Earth rather early, probably < 150 Myr after the Sun’s birth, based on observations of very old zircons (Mojzsis et al. 2001; Wilde et al. 2001), which leads to some wide interval for the water arrival.

However, most of the data suggest that the water was delivered before the late veneer, which corresponds to the latest stages of planetesimals impacts, occuring after the Moon-forming impact that is dated between 30 and 200 Myr (Kleine et al. 2009), but probably closer to 50-60 Myr (Barboni et al. 2017) after the birth of the Sun. The late veneer is expected to have contributed an additional 5 × 10⁻³ M_e of chondritic material based on measurements of highly siderophile (iron-loving) elements in the Earth’s mantle (e.g. Morbidelli & Wood 2015). However, at most 10% of the water may have been delivered to the Earth by impacts during late veneer according to geochemical and isotopic arguments (e.g. noble gas budget, Halliday 2013; Morbidelli & Wood 2015) and impact models showing that the number of impacts needed to explain the mass of terrestrial water would lead to an atmosphere that is too massive compared with the current atmospheric mass (Sinclair et al. 2020). It is therefore likely that most of the water had already been delivered before ~60 Myr.

Although impact scenarios are currently favoured to explain the delivery of water to planets, they are in fact highly contingent. Indeed, in most of these impact scenarios, it takes a complex dynamic history (e.g. resonances between planets) to move icy bodies from the outer regions of the Solar System towards the planets at the right time. Most of these scenarios require fine-tuning in one way or another. It would be interesting to find a more universal mechanism and this is the motivation for this work.

In this paper, we propose studying a new mechanism that could provide water to planets without impacts. The main idea is that the primordial asteroid belt must have contained icy asteroids, given that most C-type asteroids (which represent the largest fraction of asteroids in the main belt) are hydrated. This water ice would have sublimated in a few tens of Myr, creating a disk of gas that could then spread towards the inner planets and bring them water.

Therefore, our mechanism relies on the asteroid belt composition and mass and there are several possibilities. The young asteroid belt may have been born massive, with a mass up to a thousand times greater (e.g. Wetherill 1980) than its current mass of ~4 × 10⁻⁴ M_⊕ (e.g. Pitjeva & Pitjev 2018), according to the MMSN model (Hayashi 1981). It could also have been born with a low mass and subsequently became populated by bodies from the outer region of the Solar System (e.g. Raymond & Izidoro 2017a). The current asteroid belt is composed of two main classes of asteroids, with S-class objects being most common in the inner main belt and C-class objects dominating the outer main belt, with substantial mixing between the two populations (e.g. DeMeo & Carry 2013). C-types are spectroscopically related to carbonaceous chondrites, which generally contain ~10% water by mass (Alexander et al. 2018).

As we show in more detail in the next section (Sect. 2), young C-types may have been rich in water ice and the young asteroid belt would have had a fairly high sublimation rate, building up a disk of water gas.

It is even possible that disks of gaseous water exist outside of our Solar System, in exo-planetary systems. It would then build up from gas released in exo-asteroid belts. However, those water disks have not yet been detected but it may well be in the close future. Indeed, it would then become similar to exo-Kuiper belts, which are the equivalent of the Solar System’s Kuiper belt around other stars, that are now observed to be releasing significant quantities of gas, in particular CO and carbon (e.g. Moór et al. 2017; Cataldi et al. 2023). This phenomenon was not predicted by previous models of planetary formation and has only recently been discovered (e.g. Köspäl et al. 2013), but it shows that CO ices are present in the outer regions of extrasolar systems and that they can sublimate and produce disks of gas. These disks can propagate towards exo-planets and deliver them gaseous CO, carbon, or oxygen in a relatively efficient late accretion process (Kral et al. 2020a).

The closest observation to a disk of water in an exo-planetary system comes from Spitzer observations of HD 69830 (a star not very different from our Sun, with an asteroid belt at ~1 au), showing that water ice may be present in its exo-asteroid belt and that it should be sublimating, building up a disk of gas (Lisse et al. 2007). Gaseous water could also spread towards the inner regions of exo-planetary systems and provide water for exo-planets in the system’s habitable zone. This is a new mechanism that could be universal across all planetary systems and we explore it in this paper from a Solar System perspective. We show that current facilities, such as ALMA, may be able to detect these disks of gas and thus test this new theory using extra-solar system data.

2 Water ice in the young asteroid belt

The foundation of our model is that young asteroids were icy. This hypothesis, which is the basis of this paper may not have been examined in sufficient detail, even in the case of previous (impact) models working on a similar hypothesis.

When we look at asteroids today, they appear to be mostly free of water ice, with the exception of Ceres, the largest asteroid, where exposed water ice was detected locally (Combe et al. 2016; Platz et al. 2017; Raponi et al. 2018), and potentially Themis (Rivkin & Emery 2010; Campins et al. 2010) as well as (65) Cybele (Licandro et al. 2011), large ~200 km bodies in the outer asteroid belt. The dwarf planet Ceres contains ~20% water ice in the near-surface (Prettyman et al. 2017), and the bulk crustal average water content is estimated to be greater than 60 vol.% (Fu et al. 2017; Park et al. 2020). Moreover localised sources of water vapour have been detected (Küppers et al. 2014). In addition, Ceres appears to have an icy mantle and may even have an ocean of water beneath its surface (Castillo-Rogez et al. 2020). Therefore, from current observations it appears that water ice is only detected on a few of the larger asteroids in the outer parts of the asteroid belt. However, we note that ice on the surface is short-lived because the totality of the asteroid belt is located within the snow line (the distance at which water transitions from ice to gas; see Sect. 6.8 for a detailed discussion) so that the ice reservoir should come from further down. Moreover, it is now known that water ice is not easy to detect on surfaces of dark objects. Indeed, in situ missions such as Rosetta clearly show that even on comets, it is difficult to be observe ice on the surface. For 67P/Churyumov-Gerasimenko, Fornasier et al. (2023) showed that 67P/CG surface is dominated by refractory material and that the detected exposed water ice during the Rosetta mission represents only 0.1% of the nucleus surface. Additionally, the characteristic size of exposed chunks of water ice on 67P/CG is of the order of 1 m² or less, on average, implying that high spatial resolution is required to detect it. However, there are several compelling arguments that indicate that even if it is difficult to observe ice on the surface, water ice is abundant just below the top refractory-dominated layer, at a depth of 10 cm to 1 m, as can be seen after cliff collapses and scarps formation (e.g. Pajola et al. 2017; Fornasier et al. 2017, 2019) or by the imprint of Philae on a boulder (O’Rourke et al. 2020).

Even if ice is rare on the surface of asteroids, the observations of meteorites show that some of them are full of hydrated minerals. Given that meteorites observed today do not have water ice, their water content is defined as the amount of water that could be reformed after heating the solids; thus, it goes on to become a gas mixture, where H and O are assumed to be recom-bined to potentially form water. There is a clear segregation between the inner belt, dominated by water-poor S-type asteroids and the outer belt dominated by water-rich C-type asteroids (e.g. Kerridge 1985; Alexander et al. 2018). The water fraction in S-type asteroids is as low as 0.1%, whereas that of C-type asteroids is around 10%. It is therefore generally thought that most of the water delivered to Earth must come from C-type asteroids. We note that water in asteroids currently observed today is in the form of hydrated minerals and not in the form of ice. However, it is important to stress that the composition of today’s asteroids could be very different from their make-up in the early Solar System. It is indeed possible that most asteroids (at least the C-types rich in hydrated minerals) were initially composed of water ice that sublimated over the 4.6 Gyr of evolution of the Solar System. Several pieces of evidence point in this direction:

First, the study of meteorites show that carbonaceous chondrites (which are very likely to have C-type asteroids as parent bodies) present evidence of fluid-rock interactions because they host numerous aqueously formed minerals such as phyllosilicates, carbonates, and magnetite (Brearley 2006). Clays produced by the effect of water (Alexander et al. 1989) are also found on chondrites. It is thought that aqueous alteration may have been driven primarily by impact on water ice, rather than by radiogenic heating given the correlation of aqueous alteration with petrofabric strength (most likely arising due to shock deformation; e.g. Suttle et al. 2021).

The second argument comes from the study of asteroids using telescopes or in situ missions. From observations of a great number of C-type asteroids, Fornasier et al. (2014) show that aqueous alteration becomes important beyond 2.3 au in the asteroid belt. Moreover, Rivkin et al. (2015) show that 70% of C-type asteroids are hydrated. Recent in situ missions targeting two carbonaceous asteroids Bennu (OSIRIS-REx mission) and Ryugu (Hayabusa 2 mission) confirmed the presence of hydrated minerals on both (Hamilton et al. 2019; Kitazato et al. 2019) and of organic and carbonate material on Bennu (Simon et al. 2020). From the samples returned to Earth, measurements in the lab found even deeper lines, namely, even more hydrated minerals, than what could be observed from the spacecrafts (Hamilton et al. 2019; Praet et al. 2021; Pilorget et al. 2022; Hamilton et al. 2024). Bennu and Ryugu are near-Earth asteroids with semi-major axes close to 1 au, but originating from the inner main belt, around 2.1–2.5 au (Bottke et al. 2015; Wada et al. 2018). Ryugu’s isotopic analysis and paleomagnetic studies date the onset of fluid activity and hence aqueous alteration between <1.8 and 6.8 Myr after the formation of the CAI (McCain et al. 2023; Maurel et al. 2024).

Thirdly, there are numerous observations of so-called main belt comets or activated asteroids, which are asteroids showing signs of activity which may be comet-like (Hsieh & Jewitt 2006; Jewitt & Hsieh 2022). Several of these main belt comets come from the Themis family (Nesvorný et al. 2008; Hsieh et al. 2012) in the outer main belt. Fornasier et al. (2016) studied some members of the Themis family and found spectral features near 3 μm that could be associated with hydrated minerals or water ice. These features could be confirmed in the future using JWST.

Fourth, and perhaps more crucially, models show that the water snowline that currently lies beyond the asteroid belt was much closer in the past, allowing water ice to form throughout the asteroid belt (see Sect. 6.8 for a simplified model of the snowline). The snowline occurs at temperatures of around 145–170 K (Hayashi 1981; Lodders 2003) and recent simulations of proto-planetary disks show that the snowline moves rapidly inwards to stabilise around 2 au after a few million years (Baillié et al. 2015) and only changes when the disk dissipates to a new position beyond 3 au. It is even possible that the snowline moved inwards to less than 2 au as accretion slowed and the disk cooled and then got stuck (or fossilised) around 2.3–2.7 au (Morbidelli et al. 2016), which corresponds to where the separation between water-rich and water-poor asteroids occurs in the current belt (e.g. Abe et al. 2000; Morbidelli et al. 2000). Fossilisation would occur in the proto-planetary disk phase because of Jupiter blocking the flow of pebbles as it reaches the isolation mass, so that even if the inner main belt becomes cooler, it cannot be filled with icy material anymore and water ice never comes back on those asteroids (Morbidelli et al. 2016). Interestingly, this scenario predicts that extra-solar systems without giant planets should be much richer in water in their inner regions. Another possibility is that the magneto-rotational instability (MRI) was not active in the inner regions of the belt; furthermore, due to the resulting dead zone, matter can accumulate beyond the dead zone, becoming more massive and reaching higher temperatures (because of viscous heating), thereby pushing the snowline further out than in a fully MRI turbulent model (Martin & Livio 2012).

Fifth, most current scenarios for the transport of water to the planets are based on the injection into the asteroid belt of planetesimals from the outer regions (4–10 au) and their accretion onto the planets at the same time due to intense scattering (e.g. Raymond et al. 2006; Walsh et al. 2011; Raymond & Izidoro 2017a). These wet asteroids are considered to quantify the amount of water that can be delivered to Earth, but studies have never consider that wet asteroids implanted in the main belt will be icy and their subsequent evolution has been omitted. The aim of this paper is to rectify this gap and explore whether sublimated ice in the main belt is an effective means of transporting water to planets that could dominate over impacts.

Sixth, another way to learn more about ice in the asteroid belt is to study extra-solar systems with warm exo-asteroid belts. The water content can be studied at the beginning of the proto-planetary disk phase in the inner few au and the extent of the water snowline can be explored at high resolution (Facchini et al. 2024). Indeed, various observations of water vapour have found cold and warm water lines in some young disks by Herschel (van Dishoeck et al. 2021), Spitzer (Pontoppidan, et al. 2010), JWST (Banzatti et al. 2023), and also from the ground (Salyk et al. 2015). This is not surprising as water vapour is expected to be released near the snowline as dust or pebbles migrate into the disk interior, which is an important prediction of disk models during planet formation. What is even more interesting is that water ice is detected at a later stage, at the debris disk stage, when the star is on the main sequence and the disk resembles our asteroid belt. Indeed, in the system around the KOV-type star HD 69830, aged ~3–10 Gyr, there is an exo-asteroid belt close to 1 au (making it a close analogue of our Solar System), where water ice is detected on dust grains produced by collisions in the warm belt (Lisse et al. 2007). This is somewhat surprising because the system is old and the temperature at 1 au is higher than the sublimation temperature of water. However, it could be that water is sublimated in this belt, which creates a certain opacity to incident light and does not allow for much heat to reach the exo-asteroids, making them cold enough to keep their water ice. We note that it is likely that many systems with exo-asteroid belts should show the presence of water ice, particularly in the youngest systems.

3 The model

We developed a new model that is able to follow the sublimation of ice in the asteroid belt from its youth to Gyr timescales accounting for the varying luminosity of the Sun. The gas released from sublimation then makes up a gas disk composed of water, which evolves viscously within the frame of our model and spreads radially inwards and outwards towards the planets. We then follow the water gas mass that gets accreted onto the terrestrial planets for two different scenarios: 1) the asteroid belt was initially massive, as indicated by the MMSN model, and depletes early at ~50 Myr, 2) the asteroid belt started with its current mass. We describe those two scenarios in detail and what effects they correspond to in the next section.

3.1 Mass, spatial distribution and size distribution of the asteroid belt

In this paper, we explore two scenarios that lead to different initial belt masses. First, we modelled a belt of similar mass to today’s, which has only undergone very minor collisional evolution, which would be the case if it had started with the current mass of ~4 × 10⁻⁴ M_⊕ (e.g. Pitjeva & Pitjev 2018). This scenario is consistent with the likely end result of the “Grand Tack" (Walsh et al. 2011) or the “empty belt model" (which gets refilled by planet interactions, Raymond & Izidoro 2017b).

The second scenario is based on a more massive belt that gets depleted to current asteroid belt mass level after a time, fdep. The depletion is implemented as a step function reduction in the belt mass after t_dep. In our fiducial model, we used an initial rocky mass of 0.1 M_⊕ and t_dep = 50 Myr. This is a generic scenario assuming that the belt started with a mass similar to the MMSN model and it got emptied over time, either via an instability or dynamical depletion due to resonances with a planet, such as Jupiter; however, we note that in this latter case, the depletion would be smoother than assumed here, but our results would not change much (as discussed later in this paper).

The massive case scenario is consistent with the early instability model (that can deplete the asteroid belt at the 99–99.9% level, Clement et al. 2018, 2019), which is an updated and amended version of the “Nice model" (Gomes et al. 2005) for which the instability occurs much earlier between 10 and 100 Myr. Some newer studies have even asserted that the instability may have occurred between 60 and 100 Myr (Avdellidou et al. 2024), but this is still strongly debated, as explained in Izidoro et al. (2024). We note that our model is generic and not limited to an early instability as an instability may not be necessary to deplete the initially massive asteroid belt. Indeed, natural processes such as mean motion resonances with Jupiter (Raymond et al. 2006) or the ν₆ secular resonance with Saturn at 2.1 au (Morbidelli et al. 2000) would be able to clear the belt over long timescales. This clearing can operate over less than one Gyr and deplete more than 99% of the belt (Petit et al. 2001).

We emphasise that the primordial asteroid belt may have been even more massive; according to previous works, it could have reached ~0.5 M_⊕ (e.g. Morbidelli et al. 2000; Petit et al. 2001). Furthermore, simulations accounting for self-gravity or planet stirring find an upper limit of 2 M_⊕ (e.g. O’Brien et al. 2007; Deienno et al. 2018; Clement et al. 2019).

We base our fiducial model on the characteristics of the spatial and size distributions of the current asteroid belt, which we describe here. The main belt extends from ~2 to ~3.3 au (between the strong 4:1 and 2:1 resonances creating Kirkwood gaps). These are the values we use in our model to define the radial extent of the asteroid belt. The radial distribution of the early asteroid belt is defined using the MMSN model such that the surface density varies as r^−3/2. The size distribution we use is similar to that of the current asteroid belt with three power laws to define three regimes between different solid body diameters: D_min = 1 m, D_med =20 km, D_big =120 km, and D_max = 1000 km, with slopes q_low/ for D_min <D< D_med, q_med for D_med < D < D_big, and q_high for D_big < D < D_max. In our fiducial simulation, we use q_low = 3.6, q_med = 1.2, and q_high = 4.5 (Bottke et al. 2005). This is represented in Fig. 1, where we note N the number of solid bodies in each size bin. Later in this paper, we discuss the fact that given the very long collisional timescales of the largest bodies in the belt, the size distribution is not evolving much during the lifetime of the Solar System.

Based on observations, the separation between wet and dry asteroids is around 2.3 au (Fornasier et al. 2014, or up to 2.7 au, see the fossilised snow line in Morbidelli et al. 2016) and we therefore initially only place ice between 2.3 and 3.3 au in our model. Because there are no signs of hydration below 2.3 au, we consider that asteroids there may have never hosted water ice.

Fig. 1 Particle size distribution of bodies in the asteroid belt used in our model represented by the number of bodies in each diameter bin. There are three regimes with different slopes as described in the main text.

3.2 Composition of primordial asteroids

Many arguments for the presence of water ice in the young asteroid belt were already provided in Section 2. According to models, the snow line in the Solar System has fallen below 2 au (e.g. Baillié et al. 2015) in the proto-planetary disk phase. Thus, asteroids would have been cold enough in this young disk phase (<5 Myr) to contain water ice, which would be released as gas later when the young disk dissipated and the snowline moved outwards when entering the debris disk phase. Indeed, as the snow line is now beyond the asteroid belt (e.g. Morbidelli et al. 2016), this ice should have sublimated and should no longer be found there, if it has been sublimated within the age of the Solar System. In Section 6.8, we show a simple model of the temperature evolution from the proto-planetary disk phase to the debris disk phase, which explains rather straightforwardly how the snow line moves between the different phases.

Lodders (2003) gives a theoretical estimation of the fraction of water that may condense beyond the snowline in an MMSN model and find that we should expect ~50% water by mass, which is indeed similar to comets (Bockelée-Morvan et al. 2015). As for asteroids, it is clear that liquid water was present initially given the presence of hydrated asteroids as liquid water needs to interact with rock to create the hydrated minerals (Doyle et al. 2015; Vacher et al. 2017). Moreover, liquid water is not stable on the surface of asteroids due to the low gas pressures there and the original form of water would have been ice. Indeed, the most likely cause of the presence of liquid water is the melting of water ice following impacts (e.g. Suttle et al. 2021).

To investigate how much water might have interacted with minerals, oxygen isotopes can be used to obtain initial water/rock ratios. Using this technique, Marrocchi et al. (2018) find that the water/rock ratio can vary between 1% and 40% depending on the type of carbonaceous chondrites (increasing in that order CO, CV, CR, and CM), showing some heterogeneity in the liquid water content. This gives an idea of the potential initial ice content that must have been greater than the liquid content.

However, there are no direct connection between the liquid and icy content and it is very difficult to determine the initial amount of water ice (that has since sublimated), as it leaves no trace on the surface of asteroids or meteorites. As explained below, one way of observing the early sublimation of water ice would be to study extra-solar systems and check whether this phase of water gas release is common. Another way would be to study the effect of the released gas on planet formation, as we have done in this article. For our study, given the uncertainties, we leave the amount of water ice that can sublimate as a free parameter lower than 50%, but for our fiducial model, we assume that the asteroids initially contained f_ice = 20% water by mass.

This leads to an ice mass between 2.3 and 3.3 au of ~6.4 × 10⁻⁵ M_⊕ for scenario 1 (low mass) and ~1.6 × 10⁻² M_⊕ for scenario 2 (high mass). The ice mass is added in addition to the rocky mass in our model so that we get a total solid mass of ~4.6 × 10⁻⁴ M_e for the first scenario and ~0.12 M_⊕ for the second.

3.3 Sublimation of asteroid ices

The temperature of solid bodies after the dissipation of the proto-planetary disk is assumed to be a black body temperature so that $T_{\text{bb}}^4=L_{\star }(1-A)/\left(16\pi \sigma r^2\right)$ where L_* is the solar luminosity, cc is the Stefan-Boltzmann constant, A is the albedo, which we take to be equal to 0.06 (because C-type asteroids are very dark², Delbo et al. 2023), and r is the distance of the body from the Sun. We note that the young Sun after ~3 Myr was much less luminous than it is today. We follow the evolution of luminosity as a function of time L_*(t) with a state-of-the-art model in Section 3.7. Finally, we obtain: $$T_{\text{bb}}~278\text{\hspace{0.17em}K}{\left(\frac{L_{\star }(t)}{L_{\odot }}\right)}^{1/4}{\left(\frac{r}{1\text{au}}\right)}^{-1/2}{(1-A)}^{1/4}.$$(1)

This leads to a temperature range of 130–170 K for the main belt at the beginning of the simulation when the solar luminosity is ~0.6 L_⊙ increasing to 150–190 K for 1 L_⊙. This means that sublimation can be efficient from the very beginning given that the sublimation temperature of water ice is around 145–170 K. For simplicity, we use Eq. (1) for the temperature T of both solids and gas. We will discuss later how a different gas temperature would change our results.

We assume that the sublimation rate per unit surface area of water ice Z for an icy body in vacuum is given by $$Z(T)=\frac{a_1}{\sqrt{T}}\exp \left(-\frac{a_2}{T}\right)\left[\text{in}{\text{m}}^{-2}{\text{s}}^{-1}\right],$$(2)

where a₁ = 7.08 × 10³⁵ m⁻² s⁻¹ K^1/2 and a₂ = 6062 K (Lichtenegger & Komle 1991). Therefore, the solid mass loss rate for a solid body of radius a is given by $$\dot{m}=-4\pi a^2{\mu }_w{m}_{p}Z(T),$$(3)

where μ_w = 18 is the mean molecular weight of the water molecule and m_p is the mass of a proton. We note that we keep a constant over the whole simulation. We remain agnostic concerning the exact distribution of ice with respect to rock (e.g. mixed or in an upper layer). In fact, Eq. (2) works either for the case where the entire surface of the asteroid is covered in ice, or for a mixture of rock and ice, on the assumption that the ice is sublimated at depth and the vapor produced rises rapidly to the surface. As will be clear later on, we keep track of how much ice mass has been lost over time for each asteroid, and we stop producing water once the ice reservoir has been exhausted.

We can then estimate the survival time of a certain mass of ice M_ice on an icy body of size a and density ρ such that $$t_{\text{ice}}=\frac{M_{\text{ice}}}{\dot{m}}=\frac{a\rho }{3Z(T){\mu }_w{m}_p},$$(4)

which is close to a Gyr for the largest bodies (~500 km) in the asteroid belt and, for example, 10⁵ times smaller for a 5 m size body, that is 10 000 yr. This means that the lower end of the size distribution will sublimate its water content quite quickly (and even more so near the Sun) but that there will still be a long-lived water vapour input from the sublimation of the larger bodies that will diminish with time but which can last for hundreds of Myr, as we observe in our simulations.

In the formalism of Lichtenegger & Komle (1991), the equilibrium vapour pressure is given by $$P_{\text{eq}}=Z(T)\sqrt{2\pi {\mu }_w{m}_p{k}_{b}T},$$(5)

where k_b is Boltzmann’s constant and the gas vapour pressure can be calculated from the surface density of water gas $$P_{\text{vap}}=\frac{{\Sigma }_w{k}_{b}T}{2H{\mu }_w{m}_p}.$$(6)

When P_vap < P_eq, using Eq. (3), we convert the water ice to steam until we reach P_vap = P_eq. In all subsequent simulations, the ice mass is not sufficient to reach P_eq and we convert the lost ice mass into vapour every time step and for each size bin, thus increasing Σ_w. When P_vap > P_eq, the water vapour condenses to form solid ice. In this case, excess water vapour is removed from the gas phase to maintain equilibrium vapour pressure and ice is returned to the smaller dust by adding the contribution ${\text{ΔΣ}}_{\text{cond}}=\left(P_{\text{vap\hspace{0.17em}}}-P_{\text{eq\hspace{0.17em}}}\right)2H\frac{{\mu }_w{m}_p}{kT}$ and removing an equal amount from the vapour phase. However, the case where P_vap > P_eq does not occur in the simulations we have carried out for the Solar System.

3.4 Viscous evolution of gas released from asteroids

At each time step, we make the gas evolve viscously using the following equation (Lynden-Bell & Pringle 1974): $$\frac{\partial \text{Σ}}{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[\sqrt{r}\frac{\partial }{\partial r}(\nu \text{Σ}\sqrt{r})\right]+{\dot{\text{Σ}}}_{\text{io}},$$(7)

where ν is the kinematic viscosity, Σ(r, t) is the surface density of the gas, r is the radial variable, and ${\dot{\text{Σ}}}_{\text{io}}(r,t)$ is the rate of input/output to the surface density at radius r and time t which corresponds to the gas sublimated by asteroids or recondensed on dust.

In practice, we calculate how much mass m⁺(a) = N(a)ṁ(a)Δt is lost to vapour by sublimation for the different size bins of radius a containing N(a) solid bodies, over one time step and at a specific radius r in the belt, and how much is recondensed m⁻, to finally sum over all size bins to obtain: $${\dot{\text{Σ}}}_{\text{io}}(r)=\frac{{\displaystyle {\sum }_{a_{\min }}^{a_{\max }}{m}^{+}}(a)-m^{-}}{\pi \left[{(r+\delta r)}^2-r^2\right]\text{Δ}t}~{\displaystyle \sum _{a_{\min }}^{a_{\max }}\frac{N(a)\dot{m}(a)}{2\pi r\delta r}},$$(8)

where Δt is the time step of the simulation, δr is the width of the radial bin and a_min and a_max are the mininum and maximum bin sizes considered (see Sect. 3.1).

For the viscosity ν, we use the standard a-parameterisation of Shakura & Sunyaev (1973), such that $$\nu =\alpha c_{s}H,$$(9)

where $c_s=\sqrt{k_b{T}_g/\left(\mu m_p\right)}$ is the speed of sound, T_g is the gas temperature, μ is the molecular weight of the gas mixture (which may be different from μ_w), and H = c_s/Ω is the local scale height of the disk, with Ω the orbital frequency. We use a fiducial value of α = 10⁻² in the rest of the paper and discuss other values in Sect. 5.1.

To evolve Eq. (7) forward, we use a finite difference scheme with an explicit spatial integration and implicit temporal evolution. In practice, we numerically integrate the differential equation using the python function solve_ivp from scipy.integrate with the integrator LSODA and an adaptive timestep.

3.5 Shielding, ionisation, and dissociation of water

Initially, when water ice sublimates into the gaseous phase, it is released as a water molecule, H₂O. However, under the effect of UV photons (from the star and the interstellar medium), water can dissociate or ionise over fairly short periods, of the order of 1 to 10 days, at the location of the asteroid belt. However, if enough water molecules are present, Bethell & Bergin (2009) have demonstrated that water can protect itself against photons. This water shielding mechanism not only protects water molecules, but also shields a wide range of wavelengths, as ozone does in the Earth’s atmosphere, which can allow other molecules (if present) to survive. This is a process known as water UV shielding.

The absorption cross-section of a water molecule is σ_w = 5 × 10⁻¹⁸ cm² between 91.2 and 200 nm, which is similar to that of OH (Bethell & Bergin 2009). To reach a water column density N_w such that the optical depth τ_w is 1, we need N_w = τ_w/σ_w = 1/σ_w and then the water starts to protect itself from dissociation. This corresponds to N_w = 2 × 10¹⁶ cm⁻² or a critical surface density Σ_{crit_water} = N_wμ_wm_p = 6 × 10⁻⁶ kg/m² (given that the shielding operates vertically but we note that some radial shielding may also be at work). This is somewhat similar to carbon shielding in exo-Kuiper belts (except that C shields CO in this case, Kral et al. 2019) for which the cross-section is a factor of 3 larger with σ_c = 1.6 × 10⁻¹⁷ cm². The survival time scale of water can then be approximated as t_shield = t_w exp(σ_wN_w), where t_w is the standard survival time scale of unshielded water (which would be fixed by the stellar flux when radial shielding is not strong and by the interstellar medium flux otherwise). This is well known in the proto-planetary disk community where gaseous water is observed in such quantities that shielding is necessary (Najita et al. 2007).

For the case where there is no shielding, we use the photodissociation and photo-ionisation time scales given by Huebner & Mukherjee (2015) at 1 au where we give the values as a range to take into account periods where the Sun is faint or active³:

• For H₂O – the dominant photo-dissociation channel is H+OH (0.66–1.13 days): photo-ionisation leads to H₂O⁺ (14.0–35.0 days);

• For OH – the dominant photo-dissociation channel is O(3P)+H (1.62–1.77 days according to theory, or 0.85–0.93 days from experiment);

• For O(3P) – photo-ionisation leads to O⁺ (17.6–47.5 days);

• For H – photo-ionisation leads to H⁺ (67.4–160 days).

We note that in simulations starting with a high mass of solids (scenario 2), we have a case where water is self-protecting and can accumulate. This aspect becomes clear later on when we present our results. Otherwise, water will photo-dissociate to O and H atoms and some of it could photo-ionise and produce electrons. Therefore, here we consider a mean molecular weight μ = 18 for the high mass cases (scenario 2), and μ = 6 otherwise (because of the mixture of atomic O and H, scenario 1). We note that we fixed μ to its initial value for the whole simulation.

3.6 Gas accretion onto planets

As the gas flow radially spreads inwards to the Sun and outwards to the outer planets of the Solar System, some of the gaseous material interacts with the planets by entering their Hill spheres. To model this process, we add a sink cell for each planet to calculate the mass accreted M_acc onto them. We calculated the radial mass flux through each cell by computing⁴ ${\dot{M}}_r=6\pi \sqrt{r}\frac{\partial }{\partial r}(\sqrt{r}\text{Σ}\nu )$ at the locations r_p (0.38, 0.72, 1, and 1.52 au) of the four terrestrial planets (Mercury, Venus, Earth, Mars) and multiplying by the time step Δt. We then calculated the Hill radii Rh = r_p(M_p/(3M_⊙))^1/3 of the planets (using masses M_p of 0.0553, 0.815, 1, and 0.107 M_⊕) and compared the scale height of the gas disk H to R_h. For the cases where H > R_h (so part of the gas mass cannot be accreted by the planet⁵), we multiplied M_acc by a factor R_h/H (see details in Kral et al. 2020a). For instance, at the beginning of the simulations, R_h/H equals 0.24, 0.49, 0.49, and 0.20 for Mercury, Venus, Earth, and Mars, respectively. We note that H scales as $T_g^{1/2}$, which increases in our model with L_*, so that $H\propto L_{*}^{1/8}$ and R_h/H can at most be divided by a factor 1.05 for a solar luminosity.

We also take into account the fact that 3D hydro simulations show that some gas can pass through or circulate back to where it comes from without being able to be accreted onto the planets by multiplying M_acc by a factor f_hydro = 1/2 (e.g. Kley 1999; Robert et al. 2018; Bergez-Casalou et al. 2020). Finally, we use the work of Kral et al. (2020a) to calculate the efficiency of planets to accrete gas once this gas has reached the Hill sphere. We find that the planets can cool fast enough relative to M_r for all of the incoming gas. As a consequence, we assume that all incoming gas that can be accreted is accreted and fix the cooling efficiency parameter f_cool = 1. Finally, we can compute the mass accreted onto each planet at a given time step as $$M_{\text{acc\hspace{0.17em}}}={\dot{M}}_r\text{Δ}t{f}_{\text{cool}}{f}_{\text{hydro}}\min \left(\frac{R_h}{H};1\right).$$(10)

We checked that with this new recipe to include planets in our simulations, our code still conserves mass and angular momentum throughout all the simulations.

3.7 Model of the Sun’s evolution

In our fiducial model, we start the release of gas in the asteroid belt at 5 Myr, which is a good estimate for the end of the proto-planetary disk phase. It is confirmed by observations showing that aqueous alteration may have begun between <1.8 and 6.8 Myr after the formation of the CAI (McCain et al. 2023; Maurel et al. 2024), which can only happen when the snowline moves back outwards from ~2 au in the proto-planetary disk to >3 au in the debris disk phase. We will discuss later what a change in this starting time would imply to our results.

Figure 2 shows the output of a state-of-the-art code, Cesam2k20 (Morel 1997; Morel & Lebreton 2008; Marques et al. 2013), used to model the Sun and stellar evolution in general. Cesam2k20 is used to simulate the Sun’s evolution from its youth to its current luminosity as it is described in more detail in Appendix B. In Fig. 2 showing the Sun’s luminosity as a function of time, we can see that initially the luminosity decreases because the Sun keeps on collapsing, becomes denser and cools. The Sun is initially completely convective with a large radius, but a radiative zone appears around 1 Myr. Between 20 and 90 Myr, an additional transient convective zone appears next to the core, as it is now sufficiently hot and dense to initiate the first stage of the CNO cycle (nuclear fusion converting hydrogen into helium via nuclear reactions involving carbon, nitrogen and oxygen), i.e. the conversion of ¹²C into ¹³N, while releasing a large excess of energy. This sudden excess of energy destabilizes the medium and leads to the formation of this internal convective zone next to the core. The excess energy is then released from the surface, leading to a luminosity surge that is clearly observed between 20 and 40 Myr. The second stage of the CNO cycle (conversion of ¹³N to ¹³C) requires much higher temperatures, so that the CNO cycle remains stuck at the first stage. Once all the ¹²C has been converted, this energy source runs out, the inner convective zone disappears, and luminosity drops sharply.

But in the meantime, the temperature has risen sufficiently for the proton-proton chain (another set of nuclear fusion reactions by which hydrogen can be converted into helium) to get underway. There is a lot more fuel available, but the reaction is much weaker than the CNO cycle and does not generate enough energy to destabilize the medium, so luminosity grows slowly over time to reach the current solar luminosity (dashed line in Fig. 2).

Fig. 2 Solar luminosity (in L⊙) as a function of time in Myr (see details of the solar model in Sect. 3.7). The horizontal dashed line is showing the current solar luminosity. Our gas release model starts at 5 Myr (vertical dotted line) when the proto-planetary disk dissipates and the sun’s luminosity is around 0.6 L⊙.

4 Results

4.1 Evolution of the gas disk

Gas is injected into the disk as it sublimates from the initially icy bodies in the asteroid belt. The M rate at which it is injected (integrated over all sizes and radii so that it is equal to $1/\delta r{\displaystyle {\int }_{2.3\text{au}}^{3.3\text{au}}{\dot{\text{Σ}}}_{\text{io}}}(r){\mu }_w{m}_{p}dr$ depends mainly on the initial mass of the belt once we have fixed the solar luminosity, the location of the icy belt and its size distribution. We run our models for 1 Gyr using the two scenarios presented earlier, where we start with the current asteroid belt mass ~4 × 10⁻⁴ M_⊕ (scenario 1), and with a mass of 0.1 M_⊕ that gets depleted to current asteroid belt level after 50 Myr (scenario 2).

In Fig. 3, we show the rate of injection of gas into the belt for both scenarios⁶. We can see that the rate of injection into the belt is high at the beginning (but for a short period, beware of the logscale in time), as this is when the greatest mass is available for sublimation. Very quickly, the smallest bodies become completely rocky, and the ice sublimation can then only happen for the largest bodies. A peak appears around 20 Myr as the Sun’s luminosity increases rapidly, leading to an increase in the temperature T of solid bodies $\propto L_{\star }^{1/4}$ and an increase in Ṁ, which scales as exp(−a₂/T) (see Eq. (2)). The increase in Ṁ is a factor of ~100 compared to the values before and after the surge, over a considerable period of time (tens of Myr) and it is in fact the period when most of the gas mass is released from the asteroids. The injection rates of the high (dashed line) and low mass (solid line) scenarios are a factor 250 apart, which corresponds to the ratio of the initial masses used in these two scenarios. Ṁ can reach values of the order of 10⁻³ M_⊕/Myr in scenario 2 during the surge, which is also a typical value in colder exo-Kuiper belt releasing CO (e.g. Kral et al. 2020b) and could have been even higher in the initially massive Kuiper belt beyond Neptune (Kral et al. 2021).

The released gas will spread over time due to viscous evolution. In our fiducial model, we start with an a value of 10⁻² but we later explore how a change in this parameter could affect the results. In Fig. 4, we show the viscous evolution of the gas disk from t close to zero, where we can notice that initially the gas is injected between 2.3 and 3.3 au, i.e. at the location of the initially icy asteroids, and that as time increases gas spreads inwards and outwards due to viscosity (see Kral & Latter 2016 and Cui et al. 2024 to explore what mechanisms may explain the viscosity). The disk never reaches a steady state, as Ṁ varies over time and the viscosity varies with temperature (which depends on the varying solar luminosity). By observing the evolution of the radial profile of the surface density at different times, we can see that gas first accumulates in the vapour producing 2.3–3.3 au region and the disk then progressively spreads radially. Up to ~10 Myr, its surface density decreases because of the combination of the spread and the decrease of Ṁ. Between 10 and 25 Myr, there is an upturn in the gas surface density because of the sudden increase of Ṁ (Fig. 3). Then, after 25 Myr, Ṁ decreases again and so does the disk density.

Viscous evolution takes place on time scales t_Visc ~ r²/ν ~ 300/α years in the asteroid belt. This means that for α = 10⁻² the gas injected into the main belt will take around 30 000 years to spread in significant quantities into the inner region close to the Sun. This can be thought of as the time needed to fully react to changes in Ṁ, or to erase the memory of a given input rate. If α increases (decreases) by a factor of 10, Ṁ will remain the same⁷ (as it is mostly fixed by the asteroid belt mass) but the viscous evolution timescale will become smaller (greater) by a factor of 10, thus lowering (increasing) Σ by a factor of 10.

Estimating the surface density in the gas disk is of crucial importance because two essential processes depend on its value. If the density is high enough then the gas disk can become optically thick to stellar light, effectively shielding the water molecules from photo-dissociation. The critical density at which the shielding occurs can be estimated to be Σ_{crit_water} = 6 × 10⁻⁶ kg/m² (see Sect. 3). On the contrary, for very low gas densities, water molecules can become affected by stellar wind (see Kral et al. 2023). We estimate the density value below which this can happen with the formula $${\text{Σ}}_{\text{crit\_SW}}=\frac{2H\mu m_p}{\text{Δ}R{\sigma }_{\text{col\hspace{0.17em}}}}=7\times {10}^{-8}\text{\hspace{0.17em}kg}/{\text{m}}^2\left(\frac{r}{2\text{\hspace{0.17em}au}}\right){\left(\frac{\text{Δ}r}{1.5\text{\hspace{0.17em}au}}\right)}^{-1},$$(11)

where we assumed T_g = 197 K (blackbody temperature) at 2 au and a scaling of r^−0.5 (it is represented as a dotted line in Fig. 4). For the collisional cross-section, σ_col, we took a polarisability of ${\alpha }_{{\text{H}}_2\text{O}}=1.501\AA$ (Olney et al. 1997) and calculated the typical radius of the particle modeled by an imaginary hard sphere (Van der Waals radius) given by (3α_X/(4π))^1/3 so that we obtain $R_{{\text{H}}_2\text{O}}=0.71$ Å or σ_col = 1.6 × 10⁻²⁰ m². We note that the typical radii for H and O are 0.54 and 0.58Å, respectively, which would lead to similar cross-sections. This gives an upper limit for the critical density at which the solar wind can start pushing the gas away onto radial orbits, as taking into account collisions between protons at high velocities (>100 km/s) and H₂Owould lead to a higher cross-section due to the permanent dipole of water⁸.

We could naively assume that as scenario 2 starts with a mass that is 250 times greater, the gas disk would then evolve with a surface density that is 250 times greater, but this is not the case. Indeed, the mean molecular weight of the atomic gas (μ = 6 because the gas disk is made up of atomic H and O rather than the H₂Omolecule) is smaller in scenario 1 by a factor 3 compared to scenario 2, thus implying a greater viscosity (by a factor 3) and a lower viscous timescale by the same factor. It is similar to having a greater value of α so that the disk spreads faster and less gas can accumulate in scenario 1, compared to a gas disk with μ = 18 (scenario 2); thus the total surface density would drop and lead to an increase in the 250 factor.

In the case of high mass (scenario 2), the surface density is greater than that of the water shielding critical density, and the gas remains mainly in the form of water molecules. In the low-mass case, water photo-dissociates into H+OH, then OH dissociates into O+H and we obtain a disk composed of atomic gas with an average molecular weight of 6. However, we note that for low a values (<10⁻³), because the gas surface density gets higher, water can start shielding even in the low mass case scenario (scenario 1). For both scenarios, the gas density in the main belt may eventually (after typically a few 100 Myr) becomes low enough for the solar wind to begin blowing the gas away (depending on the exact collisional cross-section between high velocity protons and H or O), creating a weak atomic gas wind that outflows rather than a viscously spreading disk, similar to the predictions for the current Kuiper belt (Kral et al. 2021).

In Figure 5, we show the temporal evolution of the total gas mass in the disk for the low (solid) and high (dashed) mass scenarios (corresponding to scenarios 1 and 2, respectively). The mass located within 20 au is shown in blue and up to 400 au in black. The total gas mass is highest during the luminosity surge (at around 20 Myr) and can reach 2 × 10⁻⁶ M_⊕ for scenario 1 and 10⁻³ M_⊕ for scenario 2. We note that the difference in gas mass between the two scenarios is not exactly 250 for the same reasons as the surface density (i.e. different μ in the two scenarios).

Fig. 3 Injection rate of gas Ṁ from the asteroid belt into the gas disk as a function of time (since the start of sublimation, i.e. when the Sun was ~5 Myr) for scenario 1 (low mass) and scenario 2 (high mass).

Fig. 4 Evolution of the surface density profile as a function of time for scenario 1 (top) and scenario 2 (bottom). The different colours show the temporal evolution as indicated in the legend with brighter colours indicating later times. The dashed line is the density above which water is self-protected against dissociation and the dotted line is an upper limit of the critical density where gas will start to be blown away by the solar wind and create an outflow. The dash-dotted line shows the location of icy asteroids that release water in the belt.

Fig. 5 Total gas mass (in M⊕) in the disk for scenarios 1 (solid) and 2 (dashed) as a function of time (in yr). The black lines are for masses including gas up to 400 au and the blue lines for up to 20 au.

Fig. 6 Total ice mass left on solid bodies as a function of time for the scenario 1 (solid) and scenario 2 (dashed). The depletion at 50 Myr for the high mass case (scenario 2) mimics an early instability that can deplete the asteroid belt.

4.2 Evolution of the ice content

It is also instructive to examine the temporal evolution of ice left in solid bodies as a function of time and size. Fig. 6 shows that the total ice mass starts from ~6.4 × 10⁻⁵ M_⊕ (low mass, solid line) and 1.6 × 10⁻² M_⊕ (high mass, dashed line) for the scenarios 1 and 2, respectively. The qualitative evolution is then very similar in both cases with an almost constant level for the first 10 Myrs, followed by a sharp decrease during the luminosity surge around 20 Myr after the start of sublimation, and then a slow decrease up until 1 Gyr where only ~8% of the initial mass remains for scenario 1. The sharp decrease at 50 Myr for the high mass case (scenario 2) is artificial as this is the moment where we choose to deplete the belt to mimic, e.g. an early instability, and the final mass is then ~0.03% of the starting mass.

Figures 7 and 8 show in more detail the evolution of the ice content as a function of radial distance and solid sizes. For the sake clarity, we only present the results for scenario 1, the results for scenario 2 being qualitatively very similar and essentially offset by a scaling factor of 250. Let us first focus on a large size bin (because small solid bodies lose ice too rapidly) and follow the ice depletion as a function of time. Fig. 7 shows the remaining ice mass at different distances in the main belt for a ~100 km body as a function of time. We can indeed verify that for such large bodies, the inner belt gets depleted inside out and becomes devoid of ice after the surge in luminosity at ~20 Myr. After 100 Myr of evolution, only asteroids (for a size ~100 km) closest to 3.3 au remain icy but they eventually lose all their ice after 1 Gyr of evolution. This means that most primordial ice should indeed have disappeared in the current asteroid belt as we confirm with the next figure showing the evolution of even larger bodies. Fig. 8 shows the remaining ice mass for different size bins (colour) and different times (line style), demonstrating that indeed the smallest objects are being depleted rapidly and that ice only remains in the outer parts of the main belt on the largest objects ≳100 km for more than 100 Myr evolution and only bodies close to 3.3 au and greater than ~300 km still have ices after 1 Gyr of evolution.

Fig. 7 Ice mass left over total initial mass in the different radial bins in the main belt for a large 100 km body as a function of time (brighter colours are for longer times) for the scenario 1. Most of the mass gets depleted from the inner region after the surge in luminosity at ~20 Myr. There is no ice left on the 100 km body after 1 Gyr evolution.

Fig. 8 Ice mass left over total initial mass in different size bins (brighter colours are for larger objects) after 10 Myr (solid), 100 Myr (dotted) and lGyr (dashed) for the scenario 1. The smallest size bin (~40 m, black colour) never shows up on the plot as it gets depleted before 10 Myr.

4.3 Accretion onto planets

We now focus on the mass accreted onto the 4 terrestrial planets Mercury (blue), Venus (orange), the Earth (green), Mars (red), and show its temporal evolution in Fig. 9 for the low (solid) and high (dashed) mass scenarios (corresponding to scenarios 1 and 2, respectively). One notices that most of the water gets accreted around 15–25 Myr after the start of the simulation, which is when the surge in the luminosity of the Sun happens. Interestingly, the amount of water accreted on the different planets is not radically different in this mechanism (up to a factor 3 difference). In fact, it mainly depends on the radial flux Ṁ_r through each cell, which is relatively constant as a function of distance from the Sun, once the disk has had time to spread viscously. Note that Mercury and Mars will naturally accrete less gas, due to their smaller R_h/H values close to 0.2 in our model while the Earth and Venus have values closer to 0.5.

The Earth is the biggest accretor even though its Rh/H value is similar to that of Venus (and equal to 0.49; see Sect. 3.6) because it is located closer to the asteroid belt and it gets served first thus reducing the gas quantity available for Venus. After 1 Gyr evolution, the total mass accreted onto all planets is 3 × 10⁻⁵ M_⊕ for scenario 1 and 7 × 10⁻³ M_⊕ for scenario 2 (approximately 250 times more), which represents ~50% of the ice mass released after 1 Gyr. The rest of the gas is made up of two components, one is the gas disk mass as shown in Fig. 5 and the other the gas lost to the star at the inner edge and at the outer border of the simulation.

In Table 1, we list the total gas masses accreted by the different planets after 1 Gyr evolution for both scenarios. To get a comparison point, we note that the total water mass on the Earth is estimated to be around 1–10 oceans of water, including water within the mantle, or between 2.3 × 10⁻⁴ to 2.3 × 10⁻³ M_⊕. The high mass case (scenario 2) result could then explain the presence of the majority of water on the Earth while the low mass case scenario (scenario 1) would only explain a portion of it. We will discuss these crucial results concerning the Earth and other planets in more detail in Sect. 6.1.

Fig. 9 Cumulative mass accreted onto the different planets (Mercury, Venus, Earth, Mars - in blue, orange, green, red, respectively) as a function of time for the scenario 1 (solid) and 2 (dashed). We see that indeed most gas gets accreted during the surge in luminosity around 20 Myr.

5 Parameter exploration

Thanks to the simplicity of the model, we do not need to run more numerical simulations and can explore the parameter space using some analytical reasonings.

5.1 Change in the a-viscosity parameter

In debris disks, a is not very well constrained as it depends on molecular viscosity or complex hydro instabilities (e.g. the VSI) or even MHD instabilities (e.g. the MRI) as explained in more detail in Kral & Latter (2016); Cui et al. (2024). A change in a will have effects on the gas disk. Indeed, the gas disk viscosity is proportional to a and the viscous timescale scales as 1/α. Decreasing α will thus lead to an increase of the residence time of gas in the disk (before it gets accreted onto the star) and the total gas surface density (given enough time) will increase by a factor α. However, the input rate of gas into the belt Ṁ does not change nor is the radial mass flux going through each cell Ṁ_r⁹ and our results concerning the total mass accreted onto planets will not be affected.

To check the last point more firmly, we ran one simulation with a lower α value. We notice that the only difference is the time for the gas to reach planets and start accreting because it takes longer to spread as expected by the greater viscous timescale t_visc ~ 300/α years at the belt location (see Sect. 4.1). A side effect is that it means that if α becomes smaller than 10⁻⁵ so that t_visc > 30 Myr, the total accreted mass on the planets would not change but it may take longer than 30 Myr to get it accreted onto the planets. Thus, it may change some of our conclusions on the period where most accretion happens, which would become decorraleted to the surge in luminosity and would happen at, e.g. ~300 Myr for α = 10⁻⁶. However, α values lower than 10⁻⁵ may not be very realistic because in low density gas disks such as those studied here, the molecular viscosity can become important leading to α values greater than 10⁻⁵ (Cui et al. 2024).

5.2 Change in gas temperature

In our model, we assumed a gas temperature, T_g, which is the same as that of the solid bodies. It means that we assumed that T_q scales as r^−1/2, which may not be correct and that the absolute temperature is given by the black body approximation, which might also be false. For instance, one could imagine that in the less dense disks (not shielded and dominated by atoms), the gas temperature could go up due to its inability to cool. In fact, it is very complicated to calculate the temperature self-consistently because many heating and cooling mechanisms can intervene (e.g. see Kral et al. 2016, 2017a), some of which depending on the unknown composition or unknown quantities. In a simplified analytical model in Kral et al. (2017a), we show that accounting for carbon photo-ionisation and recombination with electrons, the temperature does not increase exponentially but remains rather close to a black body temperature. However, there could be some large variations and therefore, we simply analyze the impact that a change in the absolute temperature and profile would have without trying to further quantify the gas temperature. First, a change in the absolute temperature of the gas will have an effect on the viscosity because $\nu =\alpha c_s^2/\text{Ω}\propto T_g$ and the end effect is similar to a change in α for the surface density (see previous paragraph). For example, an increase in T_g by a factor of 3 would lead to a decrease in the surface density of the gas by a factor of 3. This will also change the scaleheight of the disk $H=c_s/\text{Ω}\propto T_g^{1/2}$, which could affect the total mass accreted on the planets as it depends on min($\frac{R_h}{H}$; 1) (see Eq. (10)). For example, an increase in T_g by a factor of 3 would lead to a decrease in mass accreted onto the planets by a factor of 1.7 if it is in the regime where $\frac{R_h}{H}$ < 1. We also note that the gas surface density profile will be affected by a change in the temperature slope ß. Indeed, the gas surface density profile at steady state is expected to vary as ß − x with x = 3/2 interior to the belt and x = 2 exterior to it (Kral et al. 2019). On the other hand, the radial mass flux at steady-state going through each cell Ṁ_r ∝ r^3/2−x is independent of β and thus this will not imply changes on the planet mass accretion rates.

5.3 Change in initial ice mass

The final results depend on the initial ice content as shown via our low and high mass cases (scenarios 1 and 2, respectively). For instance, an initial ice mass that is twice as large will lead to a release rate Ṁ that is twice higher and twice as much mass accreted on the different planets. In a way, our study shows that the initial ice mass could not have been much greater than that for scenario 2, i.e. ~2 × 10² M_⊕ because then too much water gets accreted on, e.g. the Earth. This fixes an upper limit on the initial asteroid belt mass of $~0.1\times \frac{f_{\text{ice\hspace{0.17em}}}}{0.2}{\text{M}}_{\oplus }$, where f_ice is the fraction of ice-to-solid on initial asteroids. This is in agreement with models showing that the initial mass of the asteroid belt could have reached ~0.5 M_⊕ (e.g. Morbidelli et al. 2000; Petit et al. 2001) and is much smaller than the upper limit of 2 M_⊕ (e.g. O’Brien et al. 2007; Deienno et al. 2018; Clement et al. 2019) fixed using numerical simulations.

5.4 Change in starting luminosity

The moment in time when we start the luminosity evolution could be important. Here, we decided to start the gas release after 5 Myr of the Sun’s birth when its luminosity is close to 0.6 L_⊙. Starting between 5 and 10 Myr would not change our results because most of the gas mass is released during the surge in luminosity at 20–30 Myr. On the other hand, starting at 2 Myr when the luminosity is close to 1 L_⊙ would increase the initial solid body temperature by a factor ~1.14, which would boost the evaporation rate by a factor ~2.7 over the first few Myr. We note that this boost will not change our results by a large factor because most of the mass will still be released over the longer period of time (> 10 Myr, mind the logscale in Fig. 2) when the luminosity surge happens.

5.5 Change in albedo

If we take the current Bond albedo of C-type asteroids that is roughly 0.02 (Helfenstein & Veverka 1989; Verbiscer et al. 2022) instead of the 0.06 assumed, we would have a gas input rate that would increase by a factor 1.5 (see Eq. (1)), thus speeding up the process. On the contrary, the initial albedo may have been even greater than the fiducial albedo assumed in this study. In this case, the temperature of solid bodies would have been smaller by a factor (1 − A)^1/4, thus changing by a factor 1.04 if going from an albedo of 0.06 to 0.2. This is not a great factor but because of the exponential factor on the evaporation rate, the latter could be reduced by a factor ~0.2 and the evolution would then be slower overall. However, we note that the surface of active comets like 67P/CG are dominated by refractory material and that the detected exposed water ice represents only 0.1% of the nucleus surface (Fornasier et al. 2023) giving an albedo very close to 0.06 (Davidsson et al. 2022). Primordial icy asteroids could actually be dark, similar to active comets, and the albedo used in our model may be reasonable.

5.6 Change in accretion efficiency

In the fiducial model, we assumed that the accretion efficiency due to cooling is maximum with f_cool = 1. This is based on the model by Kral et al. (2020a) where we have computed at each time step the theoretical accretion rate Ṁ_r due to cooling and compared to the actual potential accretion rate Ṁ_r. Using a dust-free atmospheric composition (e.g. Freedman et al. 2014; Lee & Chiang 2015), we find that Ṁ_t is always greater than Ṁ_r for all planets considered, even during the luminosity surge when Ṁ_r increases greatly. Therefore, the accretion rate onto the planets is controlled by Ṁ_r rather than by cooling. However, we note that the specifics of Ṁ_t depend on the early atmospheric mass and composition of the planets, which is very hard to predict. Even though it is expected from the analytical model of Kral et al. (2020a) that Ṁ_t is rather large, a more complex numerical model of planetary accretion would be useful to consider such complex young systems, which leaves some uncertainties on the f_cool parameter.

We have also considered that the amount of mass that gets in the Hill sphere without being accreted is f_hydro = 1/2. This recipe has been used widely in the proto-planetary disk community with varying values of f_hydro and some refinements such as considering that more mass gets removed farther away in the Hill sphere (e.g. Kley 1999; Robert et al. 2018; Bergez-Casalou et al. 2020). Those results have been solidified lately with the use of 3D hydro simulations showing that gas flow is complex and can, for instance, circulate back to the circumstellar disk or just circulate around the planet (e.g. in a circumplanetary disk) without leading to accretion (e.g. Lambrechts et al. 2019; Schulik et al. 2019; Moldenhauer et al. 2022). Considering those 3D simulations, f_hydro = 1/2 seems to be on the low side but specific simulations with terrestrial planets in a debris-disk like configuration would be required to confirm that.

6 Discussion

6.1 Delivery of water to planets and the Moon

6.1.1 Mercury

Some water ice hides in Mercury close to its poles in regions of permanent shadow (e.g. Harmon et al. 2011; Neumann et al. 2013). A current total mass of water ice of 3.4 × 10⁻¹²–3.4 × 10⁻¹⁰ M_⊕ is inferred from the MESSENGER spacecraft data (Lawrence et al. 2013). Thermal models show that water ice could be stable in cold traps near the poles in permanently shaded regions over geological timescales (e.g. Paige et al. 1992) given the long-term stability of Mercury’s obliquity.

It has been suggested that external sources may be necessary to explain the water found on Mercury because of the plausible dryness of Mercury’s crust and mantle meaning that interior sources such as volcanic outgassing may prove difficult to release large quantities of water, even over billions of years (Moses et al. 1999). Impacts by interplanetary dust particles, asteroids or comets may explain the low water content on Mercury but it is not very efficient due to Mercury’s low gravity, because most of vaporized debris after impacts have velocities higher than the escape velocity and are not retained by the planet (Moses et al. 1999).

Let us quantify the fraction of the delivered water to Mercury (3 × 10⁻⁶ M_⊕ for scenario 1 and 6 × 10⁻⁴ M_⊕ for scenario 2) that could survive at the poles over geological time scales with our gas delivery mechanism. According to Lawrence et al. (2013), the area of permanently shadowed regions at both poles amounts to ~5.6 × 10¹⁰ m², which is roughly 7.5 × 10⁻⁵ of the total surface area of Mercury (~7.5 × 10¹³ m²). As an order of magnitude and assuming that the water delivered via our mechanism arrives uniformally on the planet, we find that the water trapped in the permanently shadowed regions would be 2.3 × 10⁻¹⁰ M_⊕ for scenario 1 and 4.5 × 10⁻⁸ M_⊕ for scenario 2 if all the water is retained over Gyr timescales, which is higher than the currently observed mass of 3.4 × 10⁻¹²−3.4 × 10⁻¹⁰ Me.

However, we note that some external factors could reduce the delivered water content such as the erosion by micrometeoritic impact on exposed ice deposits (Moses et al. 1999). The erosion rate is uncertain and depends on many factors. For instance, once the ice is covered with fine layers of dust, it can be protected from sputtering, micrometeorite impacts, and sublimation (Moses et al. 1999). Given the uncertainties on this erosion rate, one needs a scenario where the delivery rate is high enough that it can explain the current water mass estimates to suggest a plausible way of delivering water to Mercury, which is what our mechanism provides.

6.1.2 Venus

The present-day water atmospheric mass of Venus is ~10⁻⁹ M_⊕ (Marcq et al. 2018), which corresponds to ~5 × 10⁻⁶ Earth ocean mass. This does not account for water that could be trapped in Venus' mantle though no estimates exist to quantify that amount (see discussion in Gillmann et al. 2009). Therefore, it appears that Venus is very dry today but it could have been wetter in the past. Indeed, the D/H ratio in Venus of 120 ± 40 the telluric ratio (de Bergh et al. 1991) is extremely high, which has been used abundantly in the literature to infer a large primordial content of water accounting for Jeans escape over several Gyr (e.g. Donahue et al. 1982).

However, more recent works put that large water content into question (e.g. Grinspoon 1993). Indeed, escape mechanisms are not very well known and they could have varied over time. For instance, Venus could have suffered hydrodynamic escape in its youth, Jeans escape around 100–500 Myr, and at later times pickup ion escape (Gillmann et al. 2009). Venus express mission even shows that water (or rather O and H) is still escaping today at a significant rate due to the solar wind (Futaana et al. 2017). On the other hand, some water may not be primordial because it could have been resupplied at a later stage, either via volcanism and/or exogenic sources such as comets in the most recent ~10⁹ years of planetary evolution (Bengtsson et al. 2013). Hence, the large D/H ratio results from the complicated history of escape over very long Gyr timescales and of possible recent water sources, which are very hard to properly quantify.

Therefore, the original water content in Venus is highly unconstrained. However, given the great similarity in bulk densities and its close proximity with the Earth, it may seem reasonable that Venus and Earth started with similar water abundances, which is the best assumption we can make at present. Because the mechanism we suggest in this paper could deliver the right amount of water to the Earth, and that it would bring roughly an equal amount to Venus, we can state that our new delivery scenario would in that sense satisfy this assumption.

6.1.3 Earth

The total amount of water on the Earth is estimated to be between 2.3 and 23 × 10⁻⁴ M_⊕, including water within the mantle. Today, the escape rate of hydrogen from the Earth is very low (~3 kg/s, Catling & Kasting 2017). It is limited by the low mixing ratio of water vapor (~1–5 ppm in the upper atmosphere), due to the so-called cold trap (i.e. cloud formation) at the tropopause. Therefore, it is not very likely that the total amount of water on the Earth, once accretion was over, has changed significantly over geological time scales due to the water cycle, the cold trap and hence the very low escape rate. However, during the Archean, the escape rate of hydrogen could have been ~2 orders of magnitude higher and was controlled by the mixing ratio of CH₄. CH₄ is indeed photolysed in the upper atmosphere and was more abundant (likely 100–1000 ppm, Catling et al. 2001) during the Archean than today.

The current deuterium to hydrogen ratio (D/H) for ocean water is known very precisely to be 1.5576 ± 0.0001 × 10⁻⁴ with respect to the Vienna Standard Mean Ocean Water (VSMOW). It is still a matter of debate whether this value is representative of the bulk of Earth’s water because of some measurements of a lower D/H ratio in deep mantle materials that may be more primordial (Hallis et al. 2015). The D/H ratio is a mixture of all sources that could have contributed to Earth’s reservoirs and is abundantly used to identify the sources of Earth’s water. As already stated in Sections 1 and 2, it appears that the best match to the current ocean water D/H ratio is that of carbonaceous chondrites, in particular CI and CM types (Alexander et al. 2012) suggesting that they were the principal source of Earth’s water.

In this paper, we find that 10⁻⁵ M_⊕ (scenario 1, low mass) and 3 × 10⁻³ Me (scenario 2, high mass) of water could have been accreted onto the Earth, mostly between 20 and 30 Myr 2after the Sun’s birth, via our viscous delivery mechanism. Therefore, the high mass scenario (scenario 2) could deliver all of Earth’s water and even in the low mass case (scenario 1) this disk-delivered component cannot be neglected (~4% of an Earth ocean).

It is also important to note that most of the water content of asteroids is sublimated rapidly in our model and the D/H in the gas state is then expected to be the same as that on the solids. Hence, within our mechanism, the D/H should be the same as that of initially icy bodies in the main belt, i.e. mainly C-type asteroids, whose D/H is indeed close to that of Earth’s water (Alexander et al. 2012). We note that our mechanism does not need impact to deliver water and it would be interesting to find ways to distinguish both the impact and viscous delivery mechanisms, which we attempt in Sects. 6.5 and 6.9.

6.1.4 Mars

Today, it is estimated that 34 m GEL¹⁰ (or 8.0 × 10⁻⁷ M_⊕) of water lies in the Martian polar-layered deposits and shallow ground ice (Carr & Head 2015). It is a lower limit because some additional water ice could hide deeper underground, notably within hydrated minerals depending on assumptions on the thickness of the hydrated crust (Mustard et al. 2012).

As for the quantity of water on the early Mars, there are quite large uncertainties. If oceans were present at the surface then a GEL of 100 m to 1 km liquid water would have been required near the surface (e.g. Baker et al. 1991; Clifford & Parker 2001). By measuring the current D/H ratio in different locations on Mars, Villanueva et al. (2015) find that water has a representative D/H value enriched by a factor of ~7 with respect to VSMOW. This is indicative of great water loss and they estimate that early Mars could have had about 137 m GEL (3.2 × 10⁻⁶ M_⊕), which is in agreement with subsurface observations with the MARSIS/Mars Express radar data (Mouginot et al. 2012). More recently, Wernicke & Jakosky (2021) estimated the amount of water required to form hydrated minerals observed on Mars and found 70–860 m GEL, which corresponds to a lower limit of ~2 × 10⁻⁶ − 2 × 10⁻⁵ M_⊕.

Our model predicts that up to 1.5 × 10⁻³ Ms of water could be accreted onto Mars in scenario 2, which is roughly larger by two orders of magnitude to these lower limits of early quantities of water on Mars. Within the framework described in this paper, this could be due to several reasons: 1) either Mars' history differs from other planets concerning its dynamics (e.g. it moved significantly due to an instability), or 2) its water inventory is underestimated, or 3) the accretion efficiency may be lower than for larger planets (see discussion in Sect. 5.6) leading to less water accreted, or 4) some unknown physics is missing in the model.

6.1.5 Moon

Water ice is now detected on the Moon but its origin is still strongly debated (Li et al. 2018). The quantity of water ice in the permanently shadowed regions on the Moon is very uncertain and estimated to be around ~10⁻¹¹ M_⊕ if all cold-traps hide a ~10-m-thick pure subsurface ice deposit (Rubanenko et al. 2019) but it could be around two orders of magnitude lower based on surface detections and the LCROSS impact results (Colaprete et al. 2010; Hayne et al. 2015; Fisher et al. 2017; Li et al. 2018). The Moon’s axis may have varied over time and rotated to its current axis ~3 Gyr ago (Siegler et al. 2016), which allowed for water ice to accumulate over billions of years. The permanently shadowed regions are expected to have slightly evolved over time with a current surface of ~36000 km² for both the North and South poles taken poleward 80·¹¹.

Because the Moon is within the Hill radius of the Earth, we can use our estimated accretion rate for the Earth to derive an accretion rate for the Moon after a few corrections. Using our model (see Fig. 9), we can estimate that the Earth accretes of order a few 10⁻⁶ Me after the moon-forming impact (i.e. after 30–100 Myr). The Moon being smaller than the Earth, we need to correct for the smaller Hill radius by a factor (M_moon/M_⊕)^1/3~0.23 leading to ~5 × 10⁻⁷ M_⊕ that may accrete on the Moon. Assuming that water is accreted uniformly on the Moon and that only the water in the shadowed regions survives, we calculate that the ratio of the shadowed area to the total surface area (3.79 × 10⁷ km²) is equal to 9.5 × 10⁻⁴ and therefore find that the surviving water mass is ~5 × 10⁻¹⁰ M_⊕, which is very close to the current estimates of water ice mass on the Moon.

In addition, we note that there is ~5 × 10⁻⁷ M_⊕ of water ice that has not yet sublimated from large bodies after 1 Gyr (see Fig. 6), which is expected to turn into vapour in the next 3.6 Gyr evolution up to now, of which roughly 50% will be accreted onto planets and 15% onto Earth according to our previous results. Accounting for the differences in Hill radii between the Earth and the Moon, we find ~2 × 10⁻⁸ M_⊕ could be accreted on the Moon and ~2 × 10⁻¹¹ M_⊕ could survive in the permanently shadowed regions. Because of the change in the Moon’s axis after about 1 Gyr, it is not clear if water previously accreted onto the Moon could survive but the latter calculation shows that even after 1 Gyr, our mechanism would allow for the accretion of the right amount to explain currently observed water ice amounts on the Moon. One caveat is that it is not certain that the gas released after 1 Gyr evolution would still viscously evolve inwards towards the planets and the Moon because it may well be blown out in a wind-like fashion if the gas density is too low (see Sect. 4.1).

6.2 Migration of planets

In our model, the planets end up embedded in a disk with substantial amounts of gas. The resulting planet/gas interactions could therefore lead to the formation of a gap or to the migration of the planet (for a recent review see Paardekooper et al. 2023). We conducted calculations described in Appendix A to compute the drift rate of planets embedded in a gas disk originating from an (exo)asteroid belt. We find that there is probably not enough gas to lead to subsequent drift rates in both type I or II migrations schemes in the Solar System. We find that in some very specific conditions (e.g. massive disk, planets at large distances), the drift rate may become more important and it will be necessary to redo the calculations for each specific extra-solar system.

6.3 Water delivery in the face of the moon forming impact

The Moon likely formed as a result of a collision between a young Earth and a Mars-sized object called Theia. The Moon-forming impact happened between -30 and 200 Myr after the birth of the Solar System, as constrained by the ¹⁸²Hf–¹⁸²W chronometer (Hafnium-Tungsten, Kleine et al. 2009) and reinforced by the U/Pb systematics (Uranium-Lead, Rudge et al. 2010). Using potentially more precise uranium-lead dating of Apollo 14 zircon fragments, Barboni et al. (2017) find that the moon should have formed within 50–60 Myr after the birth of the Solar System, so that most water delivered via our mechanism would have been accreted (at 20–30 Myr) before the Moon-forming impact.

Interestingly, it is now possible to assess whether accretion of water happened before or after the moon-forming impact by exploring the terrestrial basalt/lunar oxygen isotope difference (Greenwood et al. 2018). Assuming a nearly complete mixing and isotopic homogenisation in the aftermath of the impact, Greenwood et al. (2018) find that differences in Δ¹⁷O could possibly indicate a contribution of 5–30% of material that were accreted after the moon-forming impact. This means that the bulk of Earth’s water should have been accreted before the moon-forming impact. This is quite in agreement with our model where most water gets accreted within 30 Myr.

Another hint that water was present early comes from numerical simulations of giant impacts with the Earth. Indeed, simulations show that the presence of an ocean enhances the loss of atmosphere during the giant impact phase in the final steps of planetary formation, which may explain differences in the ³⁶Ar quantities between the Earth and Venus, being 50 times greater in the latter (Genda & Abe 2005). More interestingly, simulations show that most water in oceans can be retained despite the violent impacts of this terminal accretion phase or of the moon-forming impact (Genda & Abe 2003, 2005).

Contrary to our water delivery mechanism where most water is delivered before the moon-forming impact, most impact-delivery scenarios have more uncertain timings. For instance, for the early instability scenario, because of its unknown timing (between 10 and 100 Myr), it is not clear whether the moon-forming impact happened before or after it (Joiret et al. 2024). This may be important to decipher the amount that could have gotten in the mantle as its composition can only be changed before the moon-forming impact. Indeed, after magma ocean solidification there is a decoupling of internal (mantle and crust) and atmospheric reservoirs and water will be mainly retained in the atmosphere. Before the moon-forming impact, the Earth was often in a runaway greenhouse state interacting with a magma ocean surface (e.g. Young et al. 2023) and volatile exchanges between the mantle and the atmosphere were facilitated by liquid-gas interactions at the surface with possible partial volatile retention in the solid or melted portion of the mantle (e.g. Elkins-Tanton 2008). Shortly after the impact, most volatiles are expected to have partitioned between the atmosphere and the mantle according to their solubility in silicate melt with much of the water going in the mantle, whereas CO, CO₂ and most other gases likely go into the atmosphere (e.g. Elkins-Tanton 2008). Once the mantle has solidified, the mantle no longer convects easily and the heat flow becomes lower thus passing from a runaway greenhouse atmosphere to an atmosphere where water rains down, which facilitates the creation of surface oceans.

6.4 Current status of water-delivery models and comparison to our new mechanism

The origin of Earth’s water is still to this day an unsolved mystery (O’Brien et al. 2018). The first issue is whether or not Earth’s water has a local primordial origin or was later delivered from other regions of the Solar System, which we explore in the next two sub-sections. Currently, the most popular theory is that Earth formed dry and that water was delivered by a subsequent mechanism (e.g. Albarède 2009). There exists a wide variety of later-delivery scenarios, which all imply a delivery by impact with planetesimal (sometimes called left-over planetesimals), asteroidal or cometary material, but can significantly diverge with respect to the exact mechanism causing the delivery as well as its timing. Our aim here is not to present an exhaustive summary of all existing scenarios, but to focus on key features of the most popular ones, which will serve as a reference for comparison with our new disk-delivery mechanism.

6.4.1 Local models

Some studies suggest that water was created locally, or was already there when the local (~1 au) primordial bricks that formed the Earth assembled and that it was released later via, e.g. volcanism. For instance, water could have been incorporated into olivine grains at around 1 au through adsorption directly from the gaseous nebula (e.g. Asaduzzaman et al. 2015). Another study explores the creation of water from the oxidation of an early hydrogen-rich atmosphere by FeO in the magma ocean of the Earth (Genda & Ikoma 2008). However, both scenarios would lead to a D/H ratio close to solar, and some additional unlikely processes would be needed to explain current terrestrial D/H (e.g. hydrodynamic escape, fractionation processes).

Finally, Ciesla & Lauretta (2005) suggest that early in the proto-planetary disk, some water-bearing phyllosilicate dust could drift inwards and be incorporated into planetesimals around 1 au. However, many phyllosilicates found in meteorites appear to have been formed by processes in the parent body after the proto-planetary disk had dissipated. It seems that all of the in situ or very early addition of water scenarios have problems to match data (e.g. D/H ratio). In addition, the first materials that accreted to form the Earth are expected to have been greatly reduced, which would not be the case if they contained a significant amount of water (Wood et al. 2008). On the other hand, enstatite chondrite meteorites that have almost identical isotopic compositions to those of the Earth (suggestive of significant contributions to the latter) are drier than the Earth and highly reduced, and it would be difficult to explain if significant water-bearing material was present, unless significant quantities of H atoms could have been released and reacted with oxygen to form water¹² as suggested in Piani et al. (2020).

6.4.2 Main impact models

Here, we start by listing the main scenarios for water delivery via impacts.

The “extended feeding zone model" of water delivery proposes that the giant planets formed near their current locations, and that bodies located at 2.5−4 au, beyond the primordial snowline, and thus icy, can be excited and then accreted onto a forming Earth. Indeed, simulations by Raymond et al. (2006) show that after roughly 10 Myr, large bodies with sufficient mass form in the outer regions, which can start scattering smaller bodies by close encounters. It is the end of the oligarchic phase and the outer region is no longer dynamically isolated and can interact with forming planets in the terrestrial region. Bodies in the outer regions can also see their eccentricities pumped up via resonant excitation and secular forcing by Jupiter in addition to mutual close encounters. However, this model is currently disfavored because it fails to create large mass ratios between neighbouring planets like Mars/Earth or Venus/Mercury (e.g. Raymond et al. 2009). Moreover, it would lead to too much carbonaceous-chondrite-like material accretion in disagreement with data. Indeed, it would lead to volatile element abundances and oxygen isotope ratios of the final planets that are not consistent with Earth’s values (Drake & Righter 2002; Marty 2012). To be consistent, carbonaceous chondrites should have only contributed to about 2% of the Earth’s mass (Marty 2012) but they contribute to 15–20% in the extended feeding zone model (Raymond et al. 2006).

The “Grand Tack model" proposes that Jupiter migrated inwards in the young proto-planetary disk and then was caught up by Saturn when Jupiter was at 1.5–2 au (Walsh et al. 2011). This is when the migration of both planets were reversed and they migrated back out until the disk dispersed (Masset & Snellgrove 2001). In this scenario, the inner part of the asteroid belt is filled with S-type asteroids scattered outwards from the terrestrial planet zone during Jupiter’s inward migration, while C-type asteroids have been implanted from outer orbits (> 3 au) during the outward migration. In the process, some scattered asteroids would have impacted the planets and ice-rich asteroids may have been able to deliver water early to those forming planets. In the Grand Tack, water can be delivered early during the inwards and outwards migrations (Ogihara 2023) or later during the scattering by giant planets of planetesimals located in the outer regions that can last for more than 30 Myr (O’Brien et al. 2014). In this scenario, much less carbonaceous chondritic-like material gets accreted and it remains consistent with carbonaceous chondrites only contributing to about 2% of the Earth’s total mass (Drake & Righter 2002; Marty 2012). Similar considerations regarding the isotopic ratio of Zn find that only 5% of C-types are necessary for explaining the volatiles on Earth, while 95% would come from local enstatite-chondrite like material (Steller et al. 2022; Savage et al. 2022). We also note that recent hydrodynamical simulations at low viscosities¹³ show that migrating giant planets (like Jupiter and Saturn) go inwards without ever tacking outwards, which may invalidate the Grand Tack scenario (Pierens et al. 2014; Griveaud et al. 2023). One positive aspect of this model is that it explains the small Mars size naturally as its feeding zone gets depleted very early, not allowing Mars to grow significantly.

The Early Instability model proposes that a dynamical instability was triggered by the giant planets (after the protoplanetary disk has dissipated, i.e. after a potential Grand Tack for instance) and that it quite likely happened early, between 10 and 100 Myr (Clement et al. 2018, 2019). The instability led to a phase of violent events, which would have strongly depleted the young asteroid belt, leading to impacts on the planets and subsequent water delivery as well as leading to a reduced Mars' feeding zone. It is a very compelling scenario in that it can reproduce features both in the inner and outer regions of the Solar System.

The “destabilisation by growing planets model" from Raymond & Izidoro (2017a) proposes that when giant planets are growing, there is a natural mechanism that delivers water to planets. Indeed, they explain that the scattering of planetesimals by destabilizing growing gas giant planets naturally lead to both the implantation of a variety of planetesimals that condensed between 4 and 9 au in the asteroid belt (the precursors of C-types) and water delivery to growing terrestrial planets (because objects are scattered in all directions due to growing planets).

The “late veneer" stage is the addition of mass via impacts to the Earth in the late stage of its formation (i.e. the proto-planetary disk is no longer present), notably after the moon-forming impact. It is expected from measurements of highly siderophile elements in the Earth’s mantle that an addition of 5 × 10⁻³ M_⊕ of chondritic material was delivered after the moon formed, i.e. between 30 and 200 Myr (Kleine et al. 2009) but probably closer to after 50–60 Myr (Barboni et al. 2017). It is sometimes argued that the late veneer provided volatile-rich elements but there are many reasons to think that it is actually unlikely that a dominant fraction of Earth’s volatiles (including water) could be delivered by the Late Veneer (see more details in Morbidelli & Wood 2015). Modelling the effect of impacts on the Earth’s atmosphere¹⁴, Sinclair et al. (2020) show that because most planetesimals accreted during the late veneer are left-overs from planet formation, and hence very dry, even assuming that they contain 0.001% water by mass would only explain 0.25 Earth oceans of water. In addition, Sinclair et al. (2020) show that the number of impacts needed to explain the Earth’s water content would result in too massive an atmosphere compared to current atmospheric mass. They conclude that at most 10% of water may have been delivered to Earth from impacts (including comets, asteroids and left-over planetesimals) during the late veneer, which is in agreement with geochemical and isotopic arguments also ruling out that a majority of volatiles were delivered during the late veneer (Halliday 2013; Morbidelli & Wood 2015).

Independently of which impact-scenario is the most efficient, comets have also been suggested as a potential source to deliver water given their obvious ice-rich content (e.g. Owen & Bar-Nun 1995). However, their D/H ratio is diverse and it is often stated to be on average close to twice the terrestrial value (O’Brien et al. 2018). There has been some hope for Jupiter family comets when 103P/Hartley 2 and 45P/Honda-Mrkos-Pajduvsáková were found to have a terrestrial D/H ratio (Hartogh et al. 2011; Lis et al. 2013), but 67P/CG, also a Jupiter-family comet, was found to have a D/H ratio about 3 times terrestrial (Altwegg et al. 2015) and thus more data is needed to look for any potential trend. Another major issue with comets is that the probability of collisions with, e.g. the Earth, is very small and even delivering one Earth ocean of water is thought unlikely (Morbidelli et al. 2000). Moreover, their ¹⁵N/¹⁴N ratios are systematically higher than Earth’s (Junk & Svec 1958; Marty et al. 2016), which shows that their contribution should have been less than 10% overall (Chen et al. 2019).

We note that the mechanisms listed above are not mutually exclusive and may happen in parallel (e.g. comet impacts and extended feeding zone model) or one after the other (e.g. destabilisation by growing planets, Grand Tack, and late veneer). There are other models among the impact-like scenarios that are not described here, such as the “empty primordial asteroid belt model", “pebble accretion scenario", “convergent migration model", or “Jupiter-Saturn resonant chaotic excitation model" (see the review by Raymond 2024, and references therein). Most models are capable of reproducing the broad strokes of the inner Solar System formation, while also allowing for some water to be delivered.

6.4.3 Comparison between delivery of water via impacts and accretion from a gas disk

We note that our water delivery mechanism implies a relatively late arrival of water relative to the dissipation of the proto-planetary disk, with most of it arriving when the Sun is around 25 Myr old, which is consistent with the findings of Wood et al. (2008); Rubie et al. (2015) that early material accreting to form Earth must have been highly reduced and the oxidation state increased with time (the most likely oxidant being water). It also means that the Earth can accrete from dry planetesimals (e.g. we do not need the feeding zone to extend beyond 2.5 au) and get its water at a later stage to be consistent with the aforementioned requirement that carbonaceous chondrites only contributed to about 2% of the Earth’s total mass (Drake & Righter 2002; Marty 2012). Our mechanism would also lead to the correct D/H ratio because water sublimates from C-type asteroids (Alexander et al. 2012) and to the correct ¹⁵N/¹⁴N ratio because the latter is of chondritic origin and may be already present in the building blocks of Earth (Marty 2012). Therefore, we do not see any obvious compositional or isotopic reasons that would rule out our delivery mechanism and it even seems that it allows all previous constraints to be matched.

We would like to emphasise that we do not rule out the possibility that some of the water could have been transported by impacts before or after the late veneer. However, our scenario comes on top of this, and could be less contingent. Indeed, the impact scenarios described in the previous section (Sect. 6.4.2) are mostly based on complex dynamical configurations that may require fine-tuning, such as specific resonances between planets with specific timing to obtain the correct architecture of the solar system and the correct amount of water delivered (e.g. Grand Tack, Early Instability). Our mechanism is essentially agnostic with respect to the complex underlying dynamics that can occur in the system and relies solely on the existence of an early asteroid belt, which does not require fine-tuning.

Another advantage of our mechanism is that it seems very efficient because most water ice implanted in the asteroid belt will end up on the planets (half according to our model) with about 1/6^th on the Earth (for our fiducial model) whereas one needs many icy objects scattered to be able to have a few collide with the Earth and contribute their water. But we insist that the effectiveness of our mechanism does not rule out impacts. A rigorous comparison between the efficiency of our disk-delivery mechanism and scenarios invoking impacts would be complicated because it needs to be done on a case by case basis, which is beyond the scope of this paper.

A corollary of our work is that most of the primordial ice in the young asteroid belt will be released before 30 Myr and that, consequently, later impacts will not provide much water to the planets. This potential water depletion of asteroidal projectiles is never taken into account in impact scenarios and could be important in cases where the contributions are late and mostly come from the asteroids (e.g. contribution from asteroids in the late veneer).

Another aspect is that our disk delivery mechanism could relax certain isotopic constraints. For example, if water is mainly supplied by our mechanism, there is a decoupling between the delivery mechanism of volatiles such as N, C or Zn, which may come from chondrites, and the water (some H and O), which could be acquired by means of the disk delivery mechanism. This means that, for example, the D/H ratio constraint does not need to be adjusted within the same scenario as for other volatiles, which may allow for a re-investigation of constraints deduced from isotopic abundances.

6.5 Testing the abundance of noble gases as a plausible way to discriminate between impact and disk-driven scenarios

From an observational perspective, it may be hard to distinguish the impact-type scenarios from our mechanism because they both predict that the D/H ratio should be close to that of carbonaceous chondrites, in particular CI and CM types (Alexander et al. 2012). Could other isotopes, like those of noble gases, be used to discriminate between impacts and our disk-delivery mechanism?

Several sources have been proposed for the Earth’s noble gases origin: 1) the implanted solar wind in primordial building blocks of the Earth, 2) accretion of solar nebula gas in the young forming Earth, and 3) a cometary source and/or a chon-dritic component (for a review see Mukhopadhyay & Parai 2019). An interesting point is that some of the noble gases are thought to have been delivered to Earth via impacts of comets or asteroids, which may allow us to make some distinct predictions between the impact scenarios and our mechanism, because the disk-delivery scenario does not lead to any noble gases delivery. Let us explore what we know about each of the noble gases of importance Ne, Ar, Kr, and Xe concerning their origins and whether our mechanism is consistent with noble gases observations on Earth and Mars.

For the light noble gases Ne, and Ar, there are still large uncertainties concerning their origin on the Earth. For instance, Péron et al. (2018) note that some atmospheric Neon may have been brought via chondrites and/or comets but it seems more uncertain for Argon because it is hard to decipher whether there is a significant difference between the atmospheric and mantle compositions. Neon isotopes have a solar composition in the Earth’s mantle and the origin of this is still debated between either an atmosphere captured from the solar nebula that would get dissolved in a magma ocean, or stellar wind implantation on small grains that later formed the Earth (e.g. Moreira 2013).

On the other hand, Kr in the Earth’s mantle seems to have a chondritic origin mostly originating from early building blocks that assembled into the Earth as measured recently in Péron et al. (2021). As for the Krypton isotopic composition of the Earth’s atmosphere, it is intermediate between chondritic and solar, which may be explained as a mixture of degassed Kr from the mantle and a late input of cometary material but this is still debated (Péron et al. 2021). For Xe, it also seems that a cometary contribution is needed to explain the atmospheric isotopic composition on Earth (Avice et al. 2017) and, potentially, also for the Earth’s mantle (Péron et al. 2021).

Hence, it seems possible to explain the composition of noble gases in the Earth’s mantle without late impacts, for example, for Ne (as explained above), and if the building blocks of the Earth have a chondritic composition for Kr or Xe. On the other hand, the Ne and Kr isotopic composition of the Earth’s atmosphere cannot be explained by simple outgassing from the mantle, so we probably need a late input (carbonaceous chondrites or comets) after the moon-forming impact¹⁵ (Mukhopadhyay & Parai 2019).

The mass of impacting comets needed to explain the Ne and Kr atmospheric isotopic ratios is about 3 × 10⁻⁶ Me (Halliday 2013; Marty et al. 2016; Mukhopadhyay & Parai 2019), constituting a minor portion ~0.06% of the probable mass accreted during the late veneer (e.g. Morbidelli & Wood 2015). It means that the isotopic ratios of atmospheric noble gases can be explained quite easily for the Earth with a few impacting comets during the late veneer. This is because comets are expected to be much richer in noble gases than chondrites and this minor portion of the late veneer does not influence other major volatiles (like the H or N budget). If instead there are no cometary impacts and one assumes that only carbonaceous chondrites get accreted, the noble gases present in the atmosphere are not well fitted and it would lead to significant overabundances in the terrestrial budget of platinum group elements (Mukhopadhyay & Parai 2019).

Some recent results on noble gases on Mars have also been obtained by measuring the isotopic composition of the martian meteorite Chassigny, representative of the planet’s interior (Péron & Mukhopadhyay 2022). They find chondritic krypton isotope ratios, which implies early incorporation of chondritic volatiles. But the atmosphere of Mars has different (solar-like) krypton composition, indicating that it is not a byproduct of magma ocean outgassing and may come from early accretion of solar nebula gas after formation of the mantle. Hence, no late veneer impacts are needed to explain Mars’s Kr data.

As a consequence, noble gases cannot be used to discriminate between our disk-driven mechanism and impact-driven scenarios, because only a few late impacts are needed to explain noble gases atmospheric data, which should happen whether or not there is a disk-delivery of water in addition to it. It is thus required to look for more tracers that would be able to differentiate the impact-delivery from our mechanism. We suggest that the best test to know whether our scenario is valid and quantify it better may be to look for water gas disks in extrasolar systems as we explain in further detail in Sect. 6.9. It is well possible that after discovering the ubiquity of CO gas disks at the debris disk stage in exo-Kuiper belts (e.g. Moór et al. 2017), we discover the ubiquity of water gas disks at the same stage in exo-asteroid belts.

6.6 Collisional evolution of solid bodies in the solar system

If the asteroid belt starts with a low mass of 4 × 10⁻⁴ Me, then the largest bodies do not have time to collisionally evolve significantly; namely, the change in mass for the most massive asteroids is negligible. Indeed, using Eq. (31) in Löhne et al. (2008), we can compute that for solid bodies between 100 km and 1000 km with a steep size distribution in −4.5, the colli-sional lifetime of a 200 km asteroid in the middle of the main belt (~2.65 au) is ~50 Gyr, i.e. much larger than the age of the Solar System. For a belt of 0.1 Me similar to that of the second scenario, the collisional lifetime gets down to 0.2 Gyr but it is still much longer than the 50 Myr when we deplete most of the mass in scenario 2 to mimic an early instability. We conclude that the collisional evolution in the main belt should not affect our results significantly because most of the ice mass is in large bodies > 100 km.

6.7 Ejection of dust caused by water ice sublimation

Based on 67P/CG in situ analysis, it is clear that, for comets, the sublimation of gas releases dust along with it. The dust-to-water mass ratio for escaping material is around 1.5 (likely between 0.86 and 2.6, Choukroun et al. 2020). The size distribution index of cometary dust that comes off ranges from −3.5 to −4 (Engrand et al. 2023) with a minimum size around 10 microns (with smaller sizes existing but not dominating the cross-section, Moreno et al. 2018; Güttler et al. 2019) and a maximum size between 1 mm and 1 m that depends on the gas sublimation rate (Marschall et al. 2020) because more activity can allow to lift larger grains.

However, the physics of ejection is expected to be much different for the large (>100 km) asteroids that are probably the main ice reservoir for our disk-delivery mechanism, as compared to comets that are typically 1 to a few tens of km (e.g. ~3 km for 67P/CG) from which material escape more easily. For these large 100–1000 km asteroids, the escape velocity is 30–300 times greater than for 67P/CG, and the dust mass that can be dragged away by sublimating gas is thus necessarily greatly reduced. Indeed, using the cometary model of Zakharov et al. (2018) with standard values for large asteroids, we find that the largest grain size that can escape, dragged along with gas, is around 1 micron, while all larger dust falls back onto the asteroid after ejection (hence paradoxically no bright dusty tails are expected for those large bodies). These small micron-sized grains are not expected to be very abundant on comets (Crifo & Rodionov 1997) and may not be abundant on large asteroids (Gundlach & Blum 2013). Overall, we expect that it will not be a large mass of dust that is released from large asteroids because most mass is contained in larger fragments that do not escape.

If the dust ejection level is reduced for large asteroids, there may still be activity releasing dust for bodies smaller than 100 km. Due to the flatter size distribution of asteroids between 20 and 100 km (q_med = 1.2), the mass as a function of body size varies greatly in this range. The initial ice mass in bodies smaller than 30 km is ~10⁻⁶ for scenario 1 and ~2 × 10⁻⁴ M_⊕ for scenario 2. To gauge how much dust mass could be removed from asteroids during water sublimation from bodies smaller than 30 km, we assume a dust-to-water mass ratio of 1.5 so that we find that the activity could release a dust mass of ~ 1.5 × 10⁻⁶ M_⊕ for scenario 1 and ~3 × 10⁻⁴ M_⊕ for scenario 2. Most of it will be released between 20 and 30 Myr during the surge in luminosity and its subsequent evolution will depend on the grain sizes. Because of the steep size distribution of dust released during activity (based on cometary results, see above), most grains will have a size close to the minimum size around 10 microns for which the ratio of radiation pressure to gravitational forces β_rp ~ 10⁻² (Burns et al. 1979). This dust will thus spiral inwards from the main belt to the Sun due to Poynting Robertson drag over timescales of 2.5 × 10⁵ yr (Wyatt 2005). Assuming that the dust is released over 10 Myr (between 20 and 30 Myr after the birth of the Sun) the dust mass that will be present in the solar system at any moment during the surge in luminosity is thus a factor 0.025 smaller and closer to 4 × 10⁻⁸ M_⊕ for scenario 1 and ~8 × 10⁻⁶ M_⊕ for scenario 2. We note that this component is far from negligible as it is close to typical dust masses measured in exo-asteroid belts (Chen et al. 2014) and it should be accounted for in models trying to predict the dust quantity in warm belts but also when trying to compute the amount of warm and hot dust, called exozodis, in a specific system (Kral et al. 2017b).

6.8 Consequences for extrasolar systems: A universal mechanism

Our disk-based water delivery mechanism seems more universal than impact-type scenarios, because it is less contingent and can occur regardless of the complex dynamical history of the system. It seems universal across a large range of planetary systems (including a large variety of stellar types) because asteroid belts in younger proto-planetary disks are expected to be colder than in debris disks (see the model below), leading to a snowline that always moves outwards when the primordial gas dissipates and thus allowing degassing of water ice that can now sublimate. To better quantify this displacement of the snowline, we compare a simple model of temperature evolution of a young proto-planetary disk to a black body temperature (i.e. the temperature of solids after dissipation of the young disk).

For an optically thick proto-planetary disk, where the radiation from the star impinges at an angle ϕ(r), Dullemond et al. (2001) show that the midplane temperature profile due to stellar irradition can be simplified to $$T_{\text{mid},\text{irr}}={\left(\frac{\varphi }{2}\right)}^{1/4}{\left(\frac{R_{\star }}{r}\right)}^{1/2}{T}_{\star },$$(12)

where R_* is the radius of the Sun, T_* the effective temperature of the Sun and ϕ(r) = 0.4R_*/r + 0.05r^2/7 (e.g. Chiang & Goldreich 1997; Sasselov & Lecar 2000). The midplane temperature due to accretional heating can be derived as (e.g. Lecar et al. 2006): $$T_{\text{mid},\text{acc}}={\left[\frac{3}{4}\left({\tau }_R+\frac{2}{3}\right)\right]}^{1/4}{T}_{\text{eff},\text{acc}},$$(13)

where τ_R is the Rosseland optical depth and T_eff,acc is the effective temperature corresponding to the total flux released by accretional heating with (e.g. Lecar et al. 2006): $$\sigma T_{\text{eff},\text{acc}}^4=\frac{3}{8\pi }\frac{G{M}_{\star }{\dot{M}}_{\text{acc}}}{r^3}\left(1-\sqrt{\frac{R_{\star }}{r}}\right),$$(14)

where Ṁ_acc is the accretion rate onto the star of mass M_*. Accounting for both the irradiation from the central star and accretional heating, the final midplane temperature can be calculated as $$T_{\text{mid}}^4=T_{\text{mid},\text{irr}}^4+T_{\text{mid},\text{acc}}^4,$$(15)

which we plot in Fig. 10 (black). We note that the Rosseland opacity k_r near the sublimation temperature is around 1 m²/kg (Jang-Condell 2008; Baillié et al. 2016), which we use as a simplification in our model (because the dependence of the temperature scales as k^1/4) using: τ_r = k_RΣ_MMSN with the minimum mass solar nebula density equal to Σ_MMSN = 1.7 × 10⁴(r/lau)^−3/2 kg/m². However, viscous evolution quickly reduces the initial MMSN surface density in the inner 10 au and the radial slope of the surface density evolves to −1 (e.g. Baillié et al. 2019). To be more accurate, we use the MMSN beyond 10 au and reduce it to 2000(r/lau)⁻¹ kg/m² within 10 au, which is the MMSN value after 1−2 Myr. This is in line with an Ṁ_acc value of 10⁻⁸ M_⊙/yr (e.g. Baillié et al. 2019).

Figure 10 shows the results of our temperature evolution model. We can see that accretional heating dominates when the disk is still young with a strong accretion rate onto the star (dotted line where Ṁ_acc = 10⁻⁸ M_⊙/yr), and pushes the snowline further out at around 3 au (similar to the black body temperature in blue). But when viscous heating diminishes and that the gas disk gets to lower masses, the midplane temperature decreases and the snowline (the thin horizontal solid line shows a temperature of 170 K close to the sublimation temperature) moves inwards to ~1 au. When the disk dissipates the temperature increases to a value close to the black body temperature (in blue, Eq. (1)) which is higher and the frozen ice can now sublimate (the snowline moves back out). We used realistic values of R_*, T_*, and L_* from our solar evolution model (see Sect. 3.7) and tested for two different ages of the Sun, 2 and 10 Myr, which led to very similar results and we only show the former in Fig. 10. Note that this is an approximate model but it shows results that are not at odds with more complex self-consistent models (e.g. Baillié et al. 2016).

Comparing T_mid,irr with T_bb (Eq. (1)) is instructive, as it can be used to explore this behaviour for different stellar types. Indeed, after using that $4\pi R_{\star }^2\sigma T_{\star }^4=L_{\star }$ and rearranging, we obtain T_mid,_irr/T_bb = (2ϕ(r)/(1 − A))^1/4, with ϕ, defined above, being dominated by the term 0.05r^2/7 at exo-asteroid belt locations; hence, it has values close to 0.05−0.1, making the mid-plane temperature due to irradiation around twice lower than the black body temperature for all stellar types. We note that this explains why all asteroid belts located at the right position (i.e. with temperatures that can reach >150 K once the gas has dissipated to allow for sublimation to start) should initially sublimate but the sublimation rate will evolve over long timescales in different ways, depending on the temporal evolution of the stellar luminosity. For instance, in the case of the Sun, most of the gas is released when its luminosity increases around 20 Myr and this behaviour may become even more complex for more massive stars having more phases of luminosity surge (e.g. Iben 1967) that may lead to several phases of water delivery.

The conclusion is that because of the opacity of young massive proto-planetary disks, asteroid belt locations are colder before the disk dissipates and water ice can accumulate. But when the young disk dissipates, icy asteroids can suddenly find themselves interior to the snowline and start sublimating, which creates a gas disk that can viscously expand to the planets and potentially deliver large quantities of water onto them. This means that the prerequisites for our scenario to unfold should be generically met (at least qualitatively) in young asteroid belts. It seems therefore essential to start looking for those water disks in extrasolar systems, especially because warm belts seem to be located close to the primordial snowlines and the same mechanism as explained in this section should be at work in extrasolar systems (Morales et al. 2011; Ballering et al. 2017). In Sect. 6.9, we investigate whether we could detect these water gas disks with current technology.

Fig. 10 Temperature as a function of radial distance to the star in a proto-planetary disk (black) vs. black body temperature (blue). For the proto-planetary disk case, the different linestyles represent different accretion rates onto the star (0 for the solid line, 10−8 M⊙/yr for the dotted line, and 10−10 M⊙/yr for the dash-dotted line), which corresponds to different timescales (e.g. roughly 1−2 Myr for 10−8 M⊙/yr and later for lower accretion rates). The dashed horizontal line shows a temperature of 170 K, close to the sublimation temperature of water ice. We see that the snowline moves initially inwards when the disk accretion rate gets lower and then outwards when the disk fully dissipates and the temperature of asteroids gets close to a black body temperature.

6.9 Observations of water in extrasolar systems

6.9.1 ALMA

Water gas can be targeted with ALMA and explored at high resolution to look for the snowline position in proto-planetary disks as shown recently in Facchini et al. (2024). There are at least three transitions of water of interest in bands 5 and 7 that are p-H₂O 3₁₃ – 2₂₀ and 5₁₅ – 4₂₂, at 183.31 and 325.15 GHz, respectively, and o-H₂O 10₂₉ – 9₃₆ at 321.22 GHz, where two lines are of para-water and one line of ortho-water.

It is beyond the scope of this paper to fully explore the parameter space of plausible detections with ALMA. Here, we focus here on an illustrative example which is the promising system around the star HD 69830 that hosts a warm belt at ~1 au (Smith et al. 2009) in which water ice has been detected (Lisse et al. 2007).

We ran several simulations with the radiative transfer code RADMC3D (Dullemond et al. 2012) to estimate the integrated flux of water lines in bands 5 and 7 as a function of the total water content and gas temperature. Our simulations show that the water transition line o-H₂O 10₂₉ − 9₃₆ at 321.22 GHz (band 7) is the most sensitive to detect gaseous water released from the icy belt in HD 69830. We assume LTE in our simulations. Fig. 11 shows the integrated line flux for this transition in the water mass vs. gas temperature parameter space. This shows that if the gas temperature is greater than 100 K (which is likely given the dust temperature of ~275 K in the belt), we could readily detect a total mass of water gas above that needed for water shielding (i.e. above the black dashed line). This is promising because the dust mass in the warm belt of HD 69830 is about 10 times that of the current asteroid belt, and using Fig. 5, it indicates that the gas disk around HD 69830 may indeed be more massive than 10⁻⁶ M_⊕. The comparison is appropriate given that HD 69830 is a K0V star (Tanner et al. 2015) with a belt close to 1 au but dedicated simulations would be needed to make accurate predictions. The black dashed line shows that a gas mass above 2 × 10⁻⁷ M_⊕ ensures that water will be shielded rather than photo-dissociated (see Sect. 3.5). We can only use water lines for masses above this line because we need water molecules and not their photo-dissociation products (H and O) to probe the water lines. The estimated gas mass of 10⁻⁶ M_⊕ thus shows that the gas will be composed of water and the total water mass predicted in HD 69830 is high enough for detections with ALMA even in the case of higher a values.

We note that it is intriguing that this system still shows traces of water ice given its age (3–10 Gyr) but it may be explained by several factors: 1) the star is less massive than our Sun (~0.86 M_⊕, Lovis et al. 2006) and its luminosity is about 0.6 L_⊙, meaning that sublimation will happen on longer timescales for a given radial position to the star, or 2) the stellar evolution luminosity profile may differ and a surge in luminosity may happen at different times, or 3) the initial water ice-to-rock fraction or albedo may have been higher.

This system is all the more interesting because it hosts 3 known planets, the most distant of which, HD 69830 d (at ~0.68 au with a mass greater than 18 M_⊕), may be in the habitable zone with an orbit of ~200 days (Lovis et al. 2006). The possibility that there is enough water left to be efficiently accreted by these planets is a fascinating one. We therefore suggest that ALMA could open a new avenue as it can be used to target water gas down to low levels. ALMA would even have the capability to resolve the gas disk at ~1 au by going to high resolution (e.g. C-7 configuration or more extended) because the system is located at 12.6 pc, which could be an opportunity to refine the gas models and look for gaps due to planets in the profile.

Fig. 11 Alma observation predictions to target water in HD 69830 showing the integrated flux in Jy km s−1 for the o-H2O 1029 − 936 transition at 321.22 GHz as a function of the temperature and the quantity of water. The black dashed line indicates the critical mass level when water starts self-shielding. The level at 0.012 Jykm s−1 is for a signal-to-noise ratio of 5 represented by the solid white line labelled 5c. This line is rather sensitive to low amounts of water and may allow for the detection of the first water gas disk in the near future.

6.9.2 JWST (exploratory work)

Observations with the JWST, though at lower spectral and spatial resolutions, would also be insightful. Indeed, water lines are present everywhere in the mid-infrared. Xie et al. (2023) have for example shown that, for the young (1–3 Myr) class II disk Sz 114 the JWST can detect myriads of water lines from the hot bending rovibrational band at ~6.6 μm to the cooler pure rotational lines at ~25 μm with, e.g. MIRI. Spitzer has targeted warm belts in the mid-IR and though it may have detected water ice, it did not find any water gas lines (e.g. Lisse et al. 2007). However, the JWST has a collecting surface area 58 times larger than that of Spitzer, which may allow for those faint gas lines to be detected.

There are water lines that could be targeted at small wavelength with NIRSpec in addition to the longer wavelengths with MIRI (e.g. MRS) but the full study of the best observing strategy and comparison to ALMA are out of the scope of this paper. We leave the JWST predictions for future work but note that the HD 69830 system may be the most obvious target because of the potential presence of water ice in an asteroid belt. However, it may be interesting to make predictions for a wider variety of systems and target them with the JWST in the near future as well as to look into archival data.

6.10 Consequences on the habitability of exo-Earths

Liquid water is thought to be a necessary ingredient for the development of life on a planet (e.g. McKay 2014). Exo-Earths, which are planets with a size and temperature similar to Earth’s, have a similar problem to our blue planet in that they probably form dry and need exogenous sources to provide water. Impact-type scenarios often need a complex dynamical history requiring the right resonances or planetary migration rates to scatter a certain amount of icy bodies from outer planetary regions onto the forming exo-Earths. This could mean that it is difficult to deliver water to exo-Earths and far from generic. However, our disk-delivery mechanism is less contingent and can occur in all systems with exo-asteroid belts, meaning that it may not be so difficult to deliver water to young exo-Earths in formation.

If our disk-delivery mechanism is dominant over the other scenarios, this means that future missions targeting exo-Earths may need to search for systems with exo-asteroid belts that may contain exozodiacal dust (because of dust produced by collisions in the exo-asteroid belt and of activity due to sublimation) that could add noise to the signal and prevent the detection of exo-Earths using high contrast direct imaging (e.g. Roberge et al. 2012). It is therefore important to swiftly understand how water is transported and whether we should target exo-Earths close to dusty exo-asteroid belts to optimise our chances of having water on the surfaces of the exo-Earths in the quest for detecting life. This point has profond implications as it could change the way we want to approach exo-Earth detections from an observational standpoint in the coming decades.

7 Conclusions

There are many reasons to think that young asteroids were made up of a fraction of water ice (e.g. most C-type asteroids are hydrated, which could be the signature of an earlier icy component) that sublimated over time to lead to the current situation where most asteroids appear to be ice-free. In this article, we explore the effect of these young icy asteroids sublimating over time and show that this offers a new mechanism capable of transporting water to the inner planets in the Solar System. It may even be universal across all exo-planetary systems with respect to habitable-zone planets.

We have developed a new model that allows us to follow the sublimation of these young icy asteroids over time, taking into account the evolution of the luminosity of the young Sun. Here, we show that a disk of gas filled with water can be rapidly created in the inner region of the Solar System. We follow the viscous evolution of this disk and find that it propagates inwards and outwards towards the surrounding planets. The gas disk can last for hundreds of Myr because it is predominantly sustained by the sublimation of the largest objects in the belt, which can retain ices for considerable amounts of time. We find that the mass of the gas disk reaches a maximum when the Sun’s luminosity reaches ~1 L_⊙, around 20–30 Myr after its birth. This is when most of the gas is released from the icy asteroids, which could be an ideal spot to try and detect these new types of aqueous gas disks in extrasolar systems.

In our model, we quantify the amount of water that can be accreted onto the terrestrial planets (Mercury, Venus, Earth, and Mars) and the Moon with our viscous spreading transport mechanism for two distinct scenarios: 1) the asteroid belt was initially massive, as indicated by the MMSN model, and depleted early at ~50 Myr and 2) the asteroid belt started with its current mass. We find that in the case of an initially massive belt, our water transport mechanism can explain the water content data for the terrestrial planets and the Moon in our solar system without the need of any additional source. For example, for the Earth, we find that we can easily transport ten terrestrial oceans of water, which corresponds to the total amount of water in the hydrosphere and mantle measured in the Earth. We note that our model gives an upper limit on the initial asteroid belt mass of $~0.1\times \frac{f_{\text{ice\hspace{0.17em}}}}{0.2}{\text{M}}_{\oplus }$, where f_ice is the fraction of ice-to-solid on initial asteroids. This is because a more massive belt would lead to too much water accreted onto the Earth.

Our mechanism would also be in agreement with the Earth’s D/H ratio measurements, because most of the water would initially have come from the sublimation of C-type asteroids, which have a terrestrial D/H ratio. Other important constraints such as the ¹⁵N/¹⁴N ratio, the late veneer and noble gases would not be affected by our mechanism. Our new mechanism could have an influence on the results of works attempting to fit the isotopic ratios of the various volatiles together. In fact, we propose that there is a possible decorrelation of the D/H ratio (which would be set by disk delivery) in relation to the other volatiles that could be contributed by the building blocks or subsequent impacts, which could change the conclusions of many studies.

We find that this new mechanism of disk-driven water distribution is inevitable and may in fact be universal in systems with exo-asteroid belts around a wide range of stellar types. The main requirements to trigger it are indeed that (exo)asteroids form in an initially colder environment allowing them to retain water ice and that later on the snowline moves outwards, leading to the sublimation of these ices and the creation of a viscous gas disk. These requirements should be met in most young planetary systems, where the opacity of the proto-planetary disk should maintain low-enough temperatures for the first ~5 Myr, followed by a temperature increase after the dissipation of this primordial disk, allowing ices to sublimate. This sublimation process then naturally creates a water-vapour secondary disk that inevitably spreads towards the inner region and deliver water to exoplanets.

The best observational constraint for better quantifying the importance of our mechanism for the Solar System and exo-planetary systems would be able to detect disks of aqueous gas in extrasolar systems. We estimate that ALMA (and potentially JWST) could detect these disks, provided that they are targeted at the right moment when the water gas disk is not overly diluted. More precise predictions for different stellar types would be useful to increase the probability of detecting these aqueous gas disks in the coming years.

Acknowledgements

This paper is dedicated to Gregory. We thank the referee for their helpful report. We wish to thank Louis Manchon for providing a state-of-the-art model of the solar evolution and some description about it. We also thank Kévin Baillié, Bruno Bézard, Benjamin Charnay, Sandrine Péron, Laurette Piani, and Vladimir Zakharov for interesting discussions.

Appendix A Migration rates of planets in gas disks originating in asteroid belts

The speed at which a planet migrates depends on its ability to form a substantial gap: a planet producing a small gap (i.e. with a density within the gap which is less than half of the density outside) can rapidly migrate in a Type I migration scheme, whereas a gap-opening planet migrates more slowly in a Type II migration scheme. Kanagawa et al. (2015) show that the gap depth – the ratio of the minimum density in the gap Σ_min to the unperturbed density out of the gap Σ_un,p – can be written as $$\frac{{\text{Σ}}_{\min }}{{\text{Σ}}_{\text{un},\text{p}}}=\frac{1}{1+0.04K},$$(A.1)

where, $$K={\left(\frac{M_p}{M_{\star }}\right)}^2{\left(\frac{H_p}{r_p}\right)}^{-5}{\alpha }^{-1}={\left(\frac{M_p}{M_{\star }}\right)}^2{\left(\frac{r_p{k}_B{T}_p}{G{M}_{\star }\mu m_{\text{p}}}\right)}^{-5/2}{\alpha }^{-1}\text{,}$$(A.2)

with M_p and M_* being the mass of the planet and the star respectively, H_p = c_s(r_p)/Ω_p (with Ω_p = Ω(r_p)) is the disk scale height at the planet position r_p and a is the turbulent viscosity parameter.

It is possible to estimate the total torque and the resulting migration speed from this gap depth. In the small gap cases, the total torque that the planet feels can be approximated by (for K < 5) Γ = −cΓ₀(rp)/(1 + 0.04K) ≈ −cΓ₀(r_p) (Kanagawa et al. 2018) where ${\text{Γ}}_0(r)={\left(M_p/M_{\star }\right)}^2{\left(H_p/r\right)}^{-2}{\text{Σ}}_{\text{un}}{r}^4{\text{Ω}}_p^2$ is the normalizing torque at distance r of an unperturbed disk of surface density Σ_un(r) and c is a numerical constant varying between 1 and 3, depending on α (Kanagawa et al. 2018). This torque leads to the migration of the planet and the new planet position can be estimated by (Takeuchi et al. 1996) $$\frac{d{r}_p}{dt}=\frac{2\text{Γ}}{M_p{\text{Ω}}_p{r}_p},$$(A.3)

which leads to the following drift rate in the case of small gaps: $$\begin{array}{l}\hfill \frac{d{r}_p}{dt}\approx -3\times {10}^{-6}{\left(\frac{M_{\star }}{1{\text{M}}_{\odot }}\right)}^{-1/2}\left(\frac{M_p}{1{\text{M}}_{\oplus }}\right)\left(\frac{\mu }{18}\right)\left(\frac{{\text{Σ}}_{un}}{{10}^{-3}\text{\hspace{0.17em}kg}/{\text{m}}^2}\right)\\ \hfill {\left(\frac{r_p}{1\text{au}}\right)}^{1/2}{\left(\frac{T_p}{278.3\text{\hspace{0.17em}K}}\right)}^{-1}[\text{au}/\text{Myr}].\end{array}$$(A.4)

On the other hand, for gap-opening planets, the total torque applied on the planet can be written as Γ ≈ −cΓ₀(r_p)/(0.04K) (Kanagawa et al. 2018) leading to: $$\begin{array}{l}\frac{d{r}_p}{dt}\approx -2.6\times {10}^{-7}{\left(\frac{M_{\star }}{1{\text{M}}_{\odot }}\right)}^{-1}{\left(\frac{M_p}{1{\text{M}}_{\text{Jup}}}\right)}^{-1}{\left(\frac{\mu }{18}\right)}^{-3/2}\hfill \\ \left(\frac{{\text{Σ}}_{un}}{{10}^{-4}\text{\hspace{0.17em}kg}/{\text{m}}^2}\right){\left(\frac{r_p}{5.2\text{au}}\right)}^3{\left(\frac{T_p}{122\text{\hspace{0.17em}K}}\right)}^{3/2}\left(\frac{\alpha }{{10}^{-2}}\right)\left[\text{au}/\text{Myr}\right].\hfill \end{array}$$(A.5)

In our disk with μ = 18, α = 10⁻² in the most massive disk (i.e. with Σ_un,p ~ 10⁻³ kg/m², see Fig. 4 bottom), the terrestrial planets hardly open gaps (see Table A.1). However, for smaller values of a, for instance, 10⁻³, the gaps for the Earth and Venus would become more important with Σ_min/Σ_un,p ~ 0.4. For the largest planets, such as Jupiter, deeper gaps are opened (Table A.1). For these cases, the resulting type II migration is highly inefficient, barely pushing the giant. If Jupiter indeed creates such a deep gap, we can wonder whether the asteroid-belt produced gas would manage to cross its orbit outward towards Saturn, though it does not affect our conclusions for the accretion of water on terrestrial planets. Indeed, it is possible that due to the presence of Jupiter, the gas piles up at the inner edge of its gap, preventing this water-rich gas from reaching the outer planets. We note that the presence of a gap or a sharp density drop could then be used as a signature for the presence of a planet (Bergez-Casalou & Kral 2024).

As for the drift rates of planets, Eqs. (A.4) and (A.5) show that there is not enough gas in the disk interacting with planets for them to be affected and significantly migrate in this secondary water gas disk in the case of both small and large gaps. We note that for extrasolar systems having even more massive disks for extended periods of time, or less massive gap-opening planets (e.g. sub-Neptunes) or for planets at larger distances from their stars (the migration rate scales as $r_p^3$ for gap-opening planets), the drift rate will be more important and these calculations would need to be adjusted to consider the potential migration of planets in subsequent studies.

Appendix B Details on the Sun’s model used in this paper

The Solar evolution model used for this work was computed with the Cesam2k20¹⁶ stellar evolution code (Morel 1997; Morel & Lebreton 2008; Marques et al. 2013; Manchon et al. 2024). Stellar structure equations (for the pressure, radius, temperature and luminosity as functions of the mass coordinate) are solved using a collocation method where solutions are represented as a linear combination of piecewise polynomials, projected over a B-Spline basis. The present model uses Asplund et al. (2009)’s determination of the solar chemical composition, with opacity tables from the OPAL team (Rogers & Iglesias 1992; Iglesias & Rogers 1996), adapted to these abundances. The equation of state (EoS) uses tables from the OPAL5 EoS (Rogers & Nayfonov 2002). The nuclear reaction rates follow the compilations from NACRE (Aikawa et al. 2006) except for the ¹⁴N(p, γ)¹⁵O reaction where the LUNA facility has been used (Broggini et al. 2018). The convection is modelled with the mixing-length theory formalism under the formulation of Henyey et al. (1965), taking into account the optical thickness of the convective bubble. The atmosphere is retrieved using the Hopf function q(τ) given in Hubeny & Mihalas (2014). Effects of rotation and diffusion are not considered.

The solar model was calibrated using the program Optimal Stellar Model¹⁷ (OSM), interfaced with Cesam2k20, that implements a Levenberg-Marquardt algorithm. Initial helium content Y₀ and the mixing-length parameter α_MLT were tuned in order to retrieve, at solar age, the luminosity 1L_⊙ (L_⊙ = 3.828 × 10²⁶W; Mamajek et al. 2015) and the effective temperature 5772 K of the Sun, with a maximum error of respectively 10⁻⁵ and 1K. For this model, optimal values are found to be Y₀ = 0.2553 and α_MLT = 1.631.

References

All Tables

All Figures

	Fig. 1 Particle size distribution of bodies in the asteroid belt used in our model represented by the number of bodies in each diameter bin. There are three regimes with different slopes as described in the main text.
In the text

	Fig. 2 Solar luminosity (in L⊙) as a function of time in Myr (see details of the solar model in Sect. 3.7). The horizontal dashed line is showing the current solar luminosity. Our gas release model starts at 5 Myr (vertical dotted line) when the proto-planetary disk dissipates and the sun’s luminosity is around 0.6 L⊙.
In the text

	Fig. 3 Injection rate of gas Ṁ from the asteroid belt into the gas disk as a function of time (since the start of sublimation, i.e. when the Sun was ~5 Myr) for scenario 1 (low mass) and scenario 2 (high mass).
In the text

Fig. 4 Evolution of the surface density profile as a function of time for scenario 1 (top) and scenario 2 (bottom). The different colours show the temporal evolution as indicated in the legend with brighter colours indicating later times. The dashed line is the density above which water is self-protected against dissociation and the dotted line is an upper limit of the critical density where gas will start to be blown away by the solar wind and create an outflow. The dash-dotted line shows the location of icy asteroids that release water in the belt.

In the text

	Fig. 5 Total gas mass (in M⊕) in the disk for scenarios 1 (solid) and 2 (dashed) as a function of time (in yr). The black lines are for masses including gas up to 400 au and the blue lines for up to 20 au.
In the text

	Fig. 6 Total ice mass left on solid bodies as a function of time for the scenario 1 (solid) and scenario 2 (dashed). The depletion at 50 Myr for the high mass case (scenario 2) mimics an early instability that can deplete the asteroid belt.
In the text

	Fig. 7 Ice mass left over total initial mass in the different radial bins in the main belt for a large 100 km body as a function of time (brighter colours are for longer times) for the scenario 1. Most of the mass gets depleted from the inner region after the surge in luminosity at ~20 Myr. There is no ice left on the 100 km body after 1 Gyr evolution.
In the text

	Fig. 8 Ice mass left over total initial mass in different size bins (brighter colours are for larger objects) after 10 Myr (solid), 100 Myr (dotted) and lGyr (dashed) for the scenario 1. The smallest size bin (~40 m, black colour) never shows up on the plot as it gets depleted before 10 Myr.
In the text

	Fig. 9 Cumulative mass accreted onto the different planets (Mercury, Venus, Earth, Mars - in blue, orange, green, red, respectively) as a function of time for the scenario 1 (solid) and 2 (dashed). We see that indeed most gas gets accreted during the surge in luminosity around 20 Myr.
In the text

Fig. 10 Temperature as a function of radial distance to the star in a proto-planetary disk (black) vs. black body temperature (blue). For the proto-planetary disk case, the different linestyles represent different accretion rates onto the star (0 for the solid line, 10−8 M⊙/yr for the dotted line, and 10−10 M⊙/yr for the dash-dotted line), which corresponds to different timescales (e.g. roughly 1−2 Myr for 10−8 M⊙/yr and later for lower accretion rates). The dashed horizontal line shows a temperature of 170 K, close to the sublimation temperature of water ice. We see that the snowline moves initially inwards when the disk accretion rate gets lower and then outwards when the disk fully dissipates and the temperature of asteroids gets close to a black body temperature.

In the text

Fig. 11 Alma observation predictions to target water in HD 69830 showing the integrated flux in Jy km s−1 for the o-H2O 1029 − 936 transition at 321.22 GHz as a function of the temperature and the quantity of water. The black dashed line indicates the critical mass level when water starts self-shielding. The level at 0.012 Jykm s−1 is for a signal-to-noise ratio of 5 represented by the solid white line labelled 5c. This line is rather sensitive to low amounts of water and may allow for the detection of the first water gas disk in the near future.

In the text