© The Authors 2023

1 Introduction

The light elements carbon, hydrogen, nitrogen, oxygen, and sulfur (CHNOS) play an important role in the evolution of rocky planets. For example, their absolute and relative surface and atmospheric abundances in the form of volatile molecules such as H₂O, CO₂, or N₂ play a key role in the surface conditions and, in extension, habitability of planets (e.g. Kasting et al. 1993; Kasting & Catling 2003; Kopparapu et al. 2013). The interior properties and evolution of a planet are also affected by the planetary budgets of CHNOS. CHNOS abundances have profound effects on the physical state (solid versus liquid) of planetary cores (e.g. Trønnes et al. 2019), on the melting temperatures and mineralogy of their silicate mantles (e.g. Kushiro 1969; Dasgupta & Hirschmann 2006; Hakim et al. 2019), and on volcanic outgassing speciation (e.g. Bower et al. 2022).

In order to understand the evolution of a planet, it is crucial to identify how much CHNOS a planet initially inherits from its building blocks. Planets are thought to form in a few megayears from micrometre-sized dust inferred to be present in planet-forming disks around young stars (e.g. Andrews 2020; Raymond & Morbidelli 2022). The first stage of planet formation involves the coagulation of these micron-sized dust grains into millimetre- to centimetre-sized particles through pairwise collisions (Dominik & Tielens 1997; Birnstiel et al. 2012; Krijt et al. 2016b). These may subsequently grow into planetesimals either through continuous coagulation or through gravitational collapse triggered by, for example, streaming instabilities (Youdin & Goodman 2005; Johansen et al. 2007, 2014; Okuzumi et al. 2012).

In the colder regions of the disk, a considerable fraction of the solid-phase CHNOS mass budget is likely incorporated as ices onto dust grains in molecular species such as H₂O, CO, CO₂, CH₄, NH₃, and H₂S (Boogert et al. 2015; Öberg & Bergin 2021; Krijt et al. 2022). On the one hand, evidence from the Solar System suggests that these ices present on dust grains are inherited from the interstellar medium at least to some degree (e.g. Altwegg et al. 2017; Drozdovskaya et al. 2019). On the other hand, the amount of ice is also affected by the local balance of molecule adsorption and desorption, which are highly sensitive to the local temperature, radiation field, and composition of the gas phase (e.g. Cuppen et al. 2017). As dust grains grow into larger aggregates through collisions, dynamical processes including vertical settling, radial drift and turbulent diffusion result in significant displacement of dust throughout the disk (Weidenschilling 1977; Armitage 2010; Ciesla 2010, 2011). This dynamical transport exposes individual dust grains to a wide range of local conditions, which could have profound consequences on the amount and composition of ice (Ciesla 2010, 2011).

Many modelling efforts have attempted to constrain the effects of dust transport and collisional processes, and suggest large scale transport of volatiles throughout the disk (Cuzzi & Zahnle 2004; Krijt et al. 2016a, 2018, 2020; Bosman et al. 2018; Bergner & Ciesla 2021). However, these studies usually focus on one molecular species, such as H₂O (Krijt et al. 2016a; Schoonenberg et al. 2018), CO (Kama et al. 2016; Krijt et al. 2018, Krijt et al. 2020; Van Clepper et al. 2022), or CO₂ (Bosman et al. 2018), or emphasize a limited subset of migration, collision and ice proccesses (Ciesla 2010, 2011; Krijt & Ciesla 2016; Bergner & Ciesla 2021), or disk chemistry (Krijt et al. 2020). Alternatively, the models are sometimes limited to either the radial (Schoonenberg et al. 2018; Booth & Ilee 2019) or vertical disk dimension (Krijt et al. 2016a; Krijt & Ciesla 2016).

In this work we investigate for the first time the full coupling of dynamical transport (vertical settling, radial drift, and turbulent diffusion), collisional processes (coagulation and fragmentation), and the adsorption and desorption of multiple ices. We develop a stochastic model where we track the behaviour of dust grains with an ice mantle containing multiple CHNOS -bearing molecules in response to these processes. All processes are treated in a fully coupled, 2D fashion. In our model, we follow the evolution of individual dust grains. This enables predictions for local solid-phase CHNOS budgets throughout the disk via statistical analysis of the behaviour of a large set of individual dust grain models.

We discuss our model and methods in Sect. 2, and benchmark them against a few key earlier works. Subsequently, in Sect. 3, we analyse the coupling behaviour of the various disk processes as a function of vertical and radial position throughout the disk, and we explore the behaviour of individual dust grains under coupled disk processing. We discuss the parameter sensitivities of our model in Sect. 4, and summarize our key conclusions in Sect. 5.

2 Modelling approach

In SHAMPOO (StocHAstic Monomer PrOcessOr), we follow a small tracer dust particle called a monomer, which has size¹ s_m (see Fig. 1). In practice, the monomer will usually be embedded at a certain depth z_m in a larger dust aggregate of effective radius s_a as a consequence of collisions with other monomers and dust aggregates. Therefore, we associate each monomer at any given time with a home aggregate: the dust aggregate that hosts the monomer tracked by our model. Inferences at any given time about the properties of local dust populations can subsequently be made by tracing the evolution of a large number of monomers.

The monomer and its home aggregate interact with the disk environment through a number of processes. Aerodynamic drag and turbulent diffusion may displace the home aggregate including the monomer both radially and vertically throughout the disk (e.g. Weidenschilling 1977; Armitage 2010; Ciesla 2010, 2011). In addition, collisions with other dust aggregates may result in changes in the home aggregate size s_a and monomer depth z_m through coagulation and fragmentation (e.g. Blum & Münch 1993; Dominik & Tielens 1997; Birnstiel et al. 2011). Lastly, gas phase molecules impinging on the home aggregate can be adsorbed as ice on the monomer. Here, z_m plays an important role as the monomer depth determines the probability for gas molecules impinging on the home aggregate to be able to reach the monomer. All these processes depend on the properties of the local disk environment in which the monomer is located. This local disk environment is fully described by the thermochemical disk model ProDiMo² (Woitke et al. 2009; Kamp et al. 2010; Thi et al. 2011, 2013).

We discuss the background disk model in Sect. 2.1, and elaborate on the dynamical model in Sect. 2.2. Subsequently, we introduce our collision model in Sect. 2.3, and our treatment of ice formation in Sect. 2.4.

Fig. 1 Key model idea: a monomer of radius sm embedded at some monomer depth zm inside a home aggregate of effective radius sa.

2.1 Background disk model

The processes which alter the ice mantle of the monomer all depend on the local disk environment as characterized by the spatial physical, thermal and chemical structure of the disk. For this purpose we utilize the thermochemical disk model ProDiMo (Woitke et al. 2009; Kamp et al. 2010; Thi et al. 2011, 2013). This code has been developed to calculate the local physical, thermal, and chemical structure in an azimuthally symmetric disk. We use ProDiMo to determine the local gas and dust density ρ_g, ρ_d, and temperature T_g, T_d, respectively. Time-dependent chemistry in ProDiMo also allows us to infer the local molecule number densities n_x of all species x. We here select molecular species based on their importance for the solid-phase CHNOS mass budgets. Therefore, we will restrict ourselves to volatile H₂O, CO, CO₂, CH₄, NH₃, and H₂S throughout the rest of this work. Lastly, we use the local UV radiation field χ_RT, which is calculated from the 2D radiative transfer model.

The gas density structure of the disk is calculated from a parametrized column density structure (Woitke et al. 2009, 2016). $${\text{Σ}}_{\text{g}}\left(r\right)={\text{Σ}}_0\cdot r^{-ϵ}\cdot \exp \left[{\left(\frac{r}{r_{\text{taper}}}\right)}^{2-\gamma }\right],$$(1)

where ϵ is the density profile power law index, and r_taper denotes the tapering radius beyond which an exponential cutoff occurs. Σ₀ is derived from the normalization condition $$M_{\text{disk}}=2\pi {\displaystyle \underset{0}{\overset{\infty }{\int }}{\text{Σ}}_{\text{g}}}\left(r\right)\text{\hspace{0.17em}\hspace{0.17em}d}r.$$(2)

The vertical gas density structure is given by a Gaussian $$\begin{array}{cc}{\rho }_{\text{g}}\left(r,z\right)={\rho }_0\exp \left(\frac{z^2}{2H_{\text{g}}^2\left(r\right)}\right),& H_{\text{g}}\left(r\right)=H_0{\left(\frac{r}{r_0}\right)}^{\beta }\end{array}$$(3)

Here, H_g denotes the gas scale height, H₀ the reference scale height at a reference distance r₀, and β is the flaring exponent. The normalization ρ₀ is determined via $${\text{Σ}}_{\text{g}}\left(r\right)=2{\displaystyle \underset{0}{\overset{z_{\max }\left(r\right)}{\int }}{\rho }_{\text{g}}\left(r,z\right)\text{d}z},$$(4)

where z_max(r) = 0.5r is the maximum vertical height considered in the model (Woitke et al. 2009).

The initial dust density structure is treated in a similar fashion, where the total unsettled dust density ${\rho }_{\text{d}}^{★}$ is related to the gas density as ${\rho }_{\text{d}}^{★}=\delta {\rho }_{\text{g}}$. Here, δ denotes the dust-to-gas mass ratio. However, we do account for the settling of the largest dust grains via the procedure outlined in Woitke et al. (2016). We assume that the initial unsettled dust size distribution follows the same power law everywhere in the disk: $$f\left(a\right)\{\begin{array}{ll}\propto a^{-a_{\text{pow}}}\hfill & \text{if\hspace{0.17em}}a\in \left[a_{\min },a_{\max }\right]\hfill \\ =0\hfill & \text{elsewhere.}\hfill \end{array}$$(5)

Here, a_min = 5 × 10⁻⁸ m denotes the minimum grain size. Similarly, a_max denotes the maximum grain size, which is usually determined by the behaviour of local collisional growth processes. However, we note that grain growth in the outer regions of older disks may become drift-limited instead (Birnstiel et al. 2012). Assuming dust fragmentation by turbulence-driven relative motion, Birnstiel et al. (2012) expressed the maximum grain size as $$a_{\max }=f_{\text{f}}\frac{2}{3\pi }\frac{{\text{Σ}}_{\text{g}}}{{\rho }_{\text{a}}\alpha }{\left(\frac{{\upsilon }_{\text{frag}}}{c_{\text{s}}}\right)}^2.$$(6)

Here, ƒ_f denotes an offset parameter of order unity, which is assumed $f_{\text{f}}=\frac{1}{2}$ throughout this work. ρ_a denotes the density of dust aggregates, and is assumed equal to the home aggregate density. α represents the turbulence strength (Shakura & Sunyaev 1973). υ_frag denotes the relative velocity above which aggregates undergo fragmentation, while c_s denotes the local isothermal soundspeed.

The proportionality constant in Eq. (5) is determined by requiring the total mass density resulting from integrating over all grain sizes a to be equal to the total unsettled dust mass density ${\rho }_{\text{d}}^{★}$ (Woitke et al. 2016) $${\rho }_{\text{d}}^{\star }=\frac{4\pi }{3}{\rho }_{\text{m}}{\displaystyle \underset{a_{\min }}{\overset{a_{\max }}{\int }}f\left(a\right)a^3\text{\hspace{0.17em}d}a}.$$(7)

Here, ρ_m = 2094 kg m⁻³ denotes the dust grain material density, which is assumed equal to the monomer density. The monomer and dust aggregate density ρ_a are related as ρ_a = ϕρ_m, where ϕ denotes the dust aggregate mass filling factor.

The grain sizes are sampled log-uniformly over N_bins = 100 size bins i between a_min and a_max,D. Here a_max,D denotes the largest value for a_max encountered anywhere throughout the disk, as calculated via Eq. (6). Subsequently, the dust scale height for each size bin is calculated via (Dubrulle et al. 1995) $$H_{\text{d}}=H_{\text{g}}\sqrt{\frac{\alpha }{\text{St\hspace{0.17em}}\sqrt{3}+\alpha }.}$$(8)

Here, St denotes the Stokes number, which is a function of grain size, soundspeed and the gas density in the disk midplane (see also Sect. 2.2). We note that Eq. (8) is similar to the relation between H_d and H_g derived by Youdin & Lithwick (2007). Dust only redistributes vertically, such that the surface density associated with a grain of size a at given radial distance r $${\text{Σ}}_{\text{d}}\left(a,r\right)=2{\displaystyle {\int }_0^{z_{\max }}{\rho }_{\text{d}}\left(a,r,z\right)\text{\hspace{0.17em}d}z},$$(9)

remains constant before and after settling.

The thermal structure of the gas and dust, and the local UV radiation field χ_RT are derived from the local 2D radiation field J_ν, which is calculated with ProDiMo’s radiative transfer module. The gas temperature T_g and dust temperature T_d are treated separately. The dust temperature is found by assuming radiative equilibrium for the dust, while the gas temperature is derived from a detailed heating/cooling balance. For a full description of the radiative transfer and heating and cooling model we refer the reader to Woitke et al. (2009, 2016); Thi et al. (2011); Aresu et al. (2011) and Oberg et al. (2022).

We use the time-dependent chemistry in ProDiMo to calculate the local gas phase molecule number densities n_x of molecular species x. In the models used in this work, we considered the large DIANA chemical network containing 13 elements and 235 species. Within this chemical network, ions and ices of a particular molecular species are treated as different chemical species. Furthermore, we use the adsorption energies listed in the 2012 edition of the UMIST database (McElroy et al. 2013). However, we note that it is not possible to define a single set of adsorption energies which is consistent throughout the entire disk due to the different binding energies associated with different grain surface compositions (Kamp et al. 2017). For a more elaborate discussion of this specific chemical network, we refer the reader to Kamp et al. (2017), whereas the treatment of chemistry in ProDiMo has been discussed in Woitke et al. (2009); Aresu et al. (2011).

We use a two-step approach for the chemistry. The first step involves solving the time-dependent chemistry from initially atomic chemical abundances under molecular cloud conditions. We here choose an integration time of 1.7 × 10⁵ yr, where we follow Helling et al. (2014). This value has previously been derived as the best-fit lifetime for the Taurus Molecular Cloud 1 (McElroy et al. 2013). During the second step, we use the resulting abundances from this first step as initial chemical abundances for our disk model, and subsequently evolve this for another 2 × 10⁵ yr under the local disk conditions to yield representative chemical abundances at the onset of planetesimal formation. Although time-dependent chemistry is used to inform the chemical structure of the background disk models during the simulations, the local chemical abundances in the background models are kept fixed during SHAMPOO simulations.

Throughout this work, we consider four different background disk models, distinguished by the radial behaviour of a_max outlined in Fig. 2. For three of these models, the radial behaviour of a_max is given by Eq. (6) for a fragmentation velocity υ_frag = 1, 5, 10 m s⁻¹. We will refer to these models as vFrag1, vFrag5, and vFrag10. Furthermore, we also consider a background model const where the maximum grain size is fixed as a function of r to a_max = 3 × 10⁻³ m. For all other disk model parameters, we consider the values listed in Table 1 for all four background models. These parameters are primarily setto values of the typical T Tauri model of the DIANA project (Woitke et al. 2016, 2019; Kamp et al. 2017; Dionatos et al. 2019). We also used this DIANA T Tauri model to inform our value for a_max = 3 × 10⁻³ m for the const-model. However, it is thought that the onset of the formation of the first generation of planetesimals starts in the class I disk stage (e.g. Nixon et al. 2018). Therefore, we consider a more massive disk around a younger, more luminous pre-main sequence (PMS) star with respect to the standard T Tauri model parameters (M_disk = 0.1 M_⊙, L_⋆ = 6 L_⊙). This luminosity is consistent with a 0.7 M_⊙ PMS star of 2 × 10⁵ yr and solar metallicity Z_⊕ ≈ 0.02 (Siess et al. 2000).

The above approach gives rise to the disk structures depicted in Figs. A.1–A.4, which show the density (ρ_g, ρ_d) and thermal (T_g, T_d) structures besides the UV radiation field;χ_RT for the four different background models. Furthermore, Figs. A.5–A.8 show n_x as a function of position in the disk for H₂O, CO, CO₂, CH₄, NH₃, and H₂S for the different background disk models, both in the gas and solid phase at the end of the evaluation of the time-dependent chemistry in ProDiMo. Qualitatively, the main differences between the background models are caused by the degree of settling, with background models with larger dust grains (vFrag5 and vFrag10) having their dust settled in a thinner disk compared to the gas density structure, which is identical for all four background disk models. In particular for the vFrag10-model, Fig. A.4 shows that ρ_d is significantly higher in the disk midplane due to settling, at the cost of lower values at higher z/r. As a consequence of dust settling, the midplane region which is shielded from UV radiation is thinner in disk models where a_max is larger. This results in higher temperatures T_g, T_d down to lower z/r, which also limits the ice-forming region towards lower z/r.

Fig. 2 Radial behaviour of amax for the four different background models considered throughout this work.

2.2 Dust dynamics

The key purpose of the dynamical model is to describe the evolution of the radial (r) and vertical (z) position of the monomer during every time step ∆t. Key processes which can displace the monomer and its home aggregate are aerodynamic drag and turbulent diffusion. Our model is based on the work of Ciesla (2010, Ciesla 2011). The dynamical behaviour of the monomer is fully determined by the properties of its home aggregate, and by how well the home aggregate is coupled to the gas. An important measure which characterizes the nature of how the home aggregate interacts with the surrounding gas is the Stokes number (e.g. Armitage 2010; Krijt et al. 2018; Visser et al. 2021): $$\text{St\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}}\{\begin{array}{llll}\sqrt{\frac{\pi }{8}}\frac{{\rho }_{\text{a}}{S}_{\text{a}}}{{\rho }_{\text{g}}{c}_{\text{s}}}\text{Ω}\hfill & \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if}\hfill & S_{\text{a}}<\frac{9{\lambda }_{\text{mfp}}}{4}\hfill & \text{Epstein\hspace{0.17em}drag}\hfill \\ \sqrt{\frac{\pi }{8}}\frac{4}{9{\lambda }_{\text{mfp}}}\frac{{\rho }_{\text{a}}{S}_{\text{a}}^2}{{\rho }_{\text{g}}{c}_{\text{s}}}\text{Ω}\hfill & \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if}\hfill & S_{\text{a}}\ge \frac{9{\lambda }_{\text{mfp}}}{4}\hfill & \text{Stokes\hspace{0.17em}drag.}\hfill \end{array}$$(10, 11)

Here, ρ_a and s_a denote the material density and size of the home aggregate. In addition, ρ_g and c_s denote the local gas density and isothermal soundspeed. Ω denotes the Keplerian orbital frequency $$\text{Ω\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}}\sqrt{\frac{G{M}_{\star }}{{\left(r^2+z^2\right)}^3},}$$(12)

where G is the gravitational constant and M_⋆ the mass of the host star. In addition, λ_mfp denotes the mean free path and is calculated as $${\lambda }_{\text{mfp}}=\frac{{\rho }_{\text{g}}}{\sqrt{2}\mu m_{\text{p}}{\sigma }_{\text{mol}}},$$(13)

where µ = 2.3 denotes the mean molecular weight in atomic mass units, m_p the proton mass, and σ_mol = 2 × 10⁻¹⁹ m² the mean molecular cross section (see e.g. Okuzumi et al. 2012; Krijt et al. 2018).

As the Stokes number is a measure of how decoupled the motion of the home aggregate is from the gas, it can be used to characterize the dynamical behaviour of the home aggregate due to aerodynamic drag. The radial and vertical velocities due to aerodynamic drag υ_r and υ_z are given by (e.g. Armitage 2010) $${\upsilon }_r=-2\eta r\text{Ω}\frac{\text{St}}{1+{\text{St}}^2},$$(14) $${\upsilon }_z=-\text{Ω}z\text{St.}$$(15)

η here denotes the dimensionless gas pressure gradient, given by $$\eta =-\frac{1}{2}{\left(\frac{c_{\text{s}}}{r\text{Ω}}\right)}^2\frac{\partial \ln {\rho }_{\text{g}}}{\partial \ln r}.$$(16)

Altogether, if a monomer is embedded in a home aggregate with significant Stokes number (St ≳ α), the monomer will move towards the disk midplane and radially inward as a consequence of the vertical settling and radial drift of the home aggregate.

In turbulent disks, aerodynamic drift can be countered by turbulent diffusion. Using the single-particle random walk formalism from Ciesla (2010, 2011), we calculate the new radial and vertical positions of the home aggregate and monomer after timestep ∆t as: $$r\left(t+\text{Δ}t\right)=r\left(t\right)+{\upsilon }_r{}_{,\text{eff}}\text{Δ}t+R_1{\left(\frac{2}{\xi }{D}_{\text{d}}\text{Δ}t\right)}^{1/2},$$(17) $$z\left(t+\text{Δ}t\right)=z\left(t\right)+{\upsilon }_z{}_{,\text{eff}}\text{Δ}t+R_2{\left(\frac{2}{\xi }{D}_{\text{d}}\text{Δ}t\right)}^{1/2}.$$(18)

Here, R₁ and R₂ are uniformly drawn random numbers R₁, R₂ ∊ [−1, 1]. ξ denotes the variance of the distribution from which R₁ and R₂ are drawn, which means in our case that ξ = 1/3 (Visser 1997; Ciesla 2010). υ_r,eff and υ_z,eff denote the effective radial and vertical velocities, respectively, and are given by $${\upsilon }_{r,\text{eff}}={\upsilon }_{\text{r}}+\frac{\partial D_{\text{d}}}{\partial r}+\frac{D_{\text{d}}}{{\rho }_{\text{g}}}\frac{\partial {\rho }_{\text{g}}}{\partial r},$$(19) $${\upsilon }_{z,\text{eff}}={\upsilon }_z+\frac{\partial D_{\text{d}}}{\partial z}+\frac{D_{\text{d}}}{{\rho }_{\text{g}}}\frac{\partial {\rho }_{\text{g}}}{\partial z}.$$(20)

D_d denotes the dust diffusivity and is calculated from the gas diffusivity D_g (Youdin & Lithwick 2007) $$D_{\text{d}}=\frac{D_{\text{g}}}{1+{\text{St}}^2}.$$(21)

We here estimate the gas diffusivity D_g with the turbulent viscosity, D_g = ν_turb = αc_sH_g (Shakura & Sunyaev 1973).

As a test for our dynamical model, we consider the vertical motion of 10³ monomers at fixed radius r = 10 AU in Fig. 3. Our testing approach is similar to Ciesla (2010) and Krijt & Ciesla (2016). We perform this test for monomers of three different sizes moving in the const-background disk model: A monomer of radius s_m = 5 × 10⁻⁸ m, s_m = 3 × 10⁻⁴ m, and 3 × 10⁻³ m. At r = 10 AU and z/r = 0, the midplane Stokes numbers associated with these monomer sizes are St ≈ 1.9 × 10⁻⁷, 1.2 × 10⁻³, and 1.2 × 10⁻², respectively. These three monomers correspond to a fully coupled case (St ≪ α), a partially coupled case (St ~ α), and a more decoupled case (St > α). The integration time for each monomer was chosen to be 1.1 × 10⁵ yr, with a constant timestep of 1 yr. This timestep is at least a factor 100 shorter than the turbulent stirring timescale and the settling timescale associated with the largest particles. Since the turbulent mixing timescale at r =10 AU is roughly 10⁴ yr, we omitted the first 10⁴ yr of each monomer trajectory from Fig. 3 (see Sects. 2.5 and 3.1).

The total time that these monomers spend at each height z/r should be proportional to the background density profile of dust particles with the same Stokes number, which is a Gaussian profile with its standard deviation given by the dust scale height (Eq. (8)). Figure 3 shows that this holds for the smallest monomers (5 × 10⁻⁸ m) and largest monomers (3 × 10⁻³ m), whereas Eq. (8) appears to overestimate the width of the density profile of intermediate-sized monomers (3 × 10⁻⁴ m). This deviation was also found by Krijt & Ciesla (2016), and is a consequence of Eq. (8) not taking into account variations in the Stokes number as a function of height z/r (Ciesla 2010), whereas it is in our model.

Fig. 3 Integrated residence time for monomers of sm = 5 × 10−8 m (blue), 3 × 10−4 m (orange), and 3 × 10−3 m (green) as a function of height z at r = 10 AU. The solid lines denote the histogram obtained from the trajectories of 103 monomers combined, whereas the dashed lines are the theoretical vertical density profiles expected from Eq. (8). The vertical dashed lines depict the locations where z = Hd for each monomer size.

2.3 Collisions

With the collision model, we aim to track the evolution of the home aggregate size s_m and monomer depth z_m as the home aggregate is altered by coagulation, fragmentation, and erosion. We here treat collisions in a similar fashion as Krijt & Ciesla (2016), i.e. collisions are allowed to occur randomly, while we track the effect of each collisional interaction on the home aggregate size s_m. In addition, we use a simple probabilistic model to determine the new monomer depth z_m after a collision has occurred. Thus, during every global time step ∆t, we first determine probabilistically whether a collision has occurred, and if so, what the effects are on the size of the home aggregate s_a and on the depth at which the monomer is located z_m.

In order to determine the probability of a collision occurring during time step ∆t, we are required to evaluate the collision rates for a given home aggregate of size s_a with any other particle size bin present in our background dust size distribution. The collision rate for the home aggregate with a particle of size bin i is $$C_i=n_i{\upsilon }_{\text{rel}}{\sigma }_{\text{col}}.$$(22)

Here, n_i denotes the local number density of particles of size bin i, σ_col = π(s_a + s_i)² and υ_rel, the relative velocity between the home aggregate and the collision partner. We relate n_i to the partial disk surface density associated with particles from size bin i, Σ_d,i $$n_i=\frac{{\text{Σ}}_{\text{d,}i}}{H_{\text{d,}i}\sqrt{2\pi }}\exp \left(-\frac{z^2}{2H_{\text{d,}i}^2}\right).$$(23)

Here, H_d,i denotes the dust scale height associated with dust particles from size bin i. We can obtain the total dust surface density Σ from Σ_d,i via $${\text{Σ}}_{\text{d}}={\displaystyle \sum _i^{N_{\text{bins}}}{\text{Σ}}_{\text{d,}i}}.$$(24)

The total relative velocity υ_rel associated with a collision between the home aggregate and a given collision partner constitutes contributions from Brownian motion (υ_BM), differential aerodynamic drag consisting of radial drift (υ_r), vertical settling (υ_z), azimuthal drift (υ_ϕ), and turbulent motion (υ_turb) (Birnstiel et al. 2010; Krijt et al. 2016b). $${\upsilon }_{\text{rel}}=\sqrt{{\upsilon }_{\text{BM}}^2+{\upsilon }_r^2+{\upsilon }_z^2+{\upsilon }_{\varphi }^2+{\upsilon }_{\text{turb}}^2.}$$(25)

Following the approach of Krijt et al. (2016b), we calculate the relative velocity due to Brownian motion as $${\upsilon }_{\text{BM}}=\sqrt{\frac{8k_{\text{B}}{T}_{\text{d}}\left(m_{\text{a}}+m_i\right)}{\pi m_{\text{a}}{m}_i}.}$$(26)

Here $m_{\text{a}}=\frac{4}{3}{\rho }_{\text{a}}\pi s_{\text{a}}^3$ denotes the mass of the home aggregate (with ρ_a denoting the mass density of the home aggregate), and mi denotes the mass of a collision partner from size bin i. All dust grains, whether part of a size bin or a monomer have equal material density.

Differential aeorodynamic drag constitutes a radial (υ_r), vertical (υ_z), and azimuthal (υ_ϕ) component, given by: $${\upsilon }_r=\left|{\upsilon }_{r,\text{a}}-{\upsilon }_{r,i}\right|$$(27) $${\upsilon }_z=\left|{\upsilon }_{z,\text{a}}-{\upsilon }_{z,i}\right|$$(28) $${\upsilon }_{\varphi }=\left|\frac{{\upsilon }_{r,\text{a}}}{2{\text{St}}_{\text{a}}}-\frac{{\upsilon }_{r,i}}{2{\text{St}}_i}\right|.$$(29)

Here, υ_r,n and υ_z,n (n ∊ (a, i)) are calculated from Eqs. (14) and (15), respectively. St_a denotes the Stokes number associated with the home aggregate, while St_i denotes the Stokes number of the collision partner.

For the contribution of turbulence to the total relative velocity, we use the approach derived by Ormel & Cuzzi (2007): $${\upsilon }_{\text{turb}}=\delta \upsilon \cdot Q.$$(30)

Here, $\delta \upsilon =c_{\text{s}}\sqrt{\alpha }$ denotes the mean random velocity of the largest turbulent eddies (Krijt et al. 2016b). The value of Q depends on how the stopping time t_s = St/Ω compares to the orbital period 1/Ω and the turnover timescale of the smallest scale eddies, $t_{\eta }=1/\left(\text{Ω}\sqrt{\mathrm{Re}}\right)$: $$Q=\{\begin{array}{ll}{\text{Re}}^{\frac{1}{2}}\left|{\text{St}}_1-{\text{St}}_2\right|\hfill & \text{if\hspace{0.17em}}{t}_{\text{s1}}\le t_{\eta }\hfill \\ \sqrt{\frac{1}{1+{\text{St}}_1}+\frac{1}{1+{\text{St}}_2}}\hfill & \text{if\hspace{0.17em}}{t}_{\text{s1}}\ge {\text{Ω}}^{-1}\hfill \\ 1.55\sqrt{{\text{St}}_1}\hfill & \text{otherwise.}\hfill \end{array}$$(31)

Here, we associate the Stokes number St₁ with the largest particle, either being the home aggregate or the collision partner, and St₂ with the smallest of the two. Re denotes the turbulent Reynolds number, which we calculate as the ratio between the turbulent viscosity ν_turb and molecular viscosity ν_mol: $$\mathrm{Re}=\frac{v_{\text{turb}}}{v_{\text{mol}}}=\alpha c_{\text{s}}{H}_g{\left(\sqrt{\frac{8}{\pi }\frac{c_{\text{s}}{\lambda }_{\text{mfp}}}{2}}\right)}^{-1}.$$(32)

The above calculations are performed for every size bin i, to find the associated C_i via Eq. (22). Subsequently, the total collision rate C_tot can be obtained via $$C_{\text{tot}}={\displaystyle \sum _i^{N_{\text{bins}}}{C}_i},$$(33)

such that the probability of a single collision event happening during timestep ∆t is $$P_{\text{col}}=1-\exp \left(-C_{\text{tot}}\text{Δ}t\right).$$(34)

Because this expression does not incorporate the possibility of multiple collisions, we are required to set ∆t shorter than the timescale for a single collision, ∆t ≲ 1/C_tot. This can become a problem if our home aggregate is large and located near the midplane. In this situation, the home aggregate will collide frequently with (sub)miron-sized grains, resulting in C_tot becoming very large (Krijt & Ciesla 2016).

To prevent the usage of prohibitively small values of ∆t, we group collisions with collision partners which have a mass m_i which is smaller than a fraction ƒ_c = 10⁻¹ of the home aggregate mass m_a (Ormel & Spaans 2008; Krijt et al. 2015; Krijt & Ciesla 2016). Instead of treating each collision between the home aggregate and small collision partners separately, we only allow the home aggregate to collide with a group of small collision partners at once. The effective collision rate for such a group collision will be (much) smaller than the collision rate between the home aggregate and the individual collision partners which constitute the group. We express the number of particles per group N_col as $$N_{\text{col}}=\frac{1}{f_{\text{c}}}\frac{m_i}{m_{\text{a}}},$$(35)

such that the modified collision rates become $${\tilde{C}}_i=\{\begin{array}{ll}{C}_i\hfill & \text{if\hspace{0.17em}}{m}_i/m_{\text{a}}>f_{\text{c}},\hfill \\ \frac{1}{f_{\text{c}}}\frac{m_i}{m_{\text{a}}}\cdot C_i\hfill & \text{if\hspace{0.17em}}{m}_i/m_{\text{a}}\le f_{\text{c}}.\hfill \end{array}$$(36)

We then calculate the modified total collision rate as $${\tilde{C}}_{\text{tot}}={\displaystyle \sum _i^{N_{\text{bins}}}{\tilde{C}}_i,}$$(37)

such that the probability of a collision occurring during a timestep ∆t is $${\tilde{P}}_{\text{col}}=1-\exp \left(-{\tilde{C}}_{\text{tot}}\text{Δ}t\right).$$(38)

To determine whether a collision event occurs during a timestep, we generate a random number R₃ ∊ [0, 1], such that if R₃ ≤ ${\tilde{P}}_{\text{col}}$, the monomer has undergone a collision event.

In order to determine the effect of a collision on the monomer and home aggregate, we first need to determine the mass of the collision partner, which we determine by generating another random number R₄ ∊ [0, ${\tilde{C}}_{\text{tot}}$], and sum over the various collision size bins up to bin n, which is the first bin which satisfies $${\displaystyle \sum _i^n{\tilde{C}}_i\ge R_4.}$$(39)

We then select a particle from size bin n as the collision partner. Note that this automatically means that size bins which have a higher associated collision rate are more likely to be selected.

The outcome of a collision is fully determined by the mass ratio of the collision partner and the home aggregate, and their sizes. The first step of determining the collision outcome is whether the collision is constructive (net dust growth) or destructive (net dust fragmentation). We assume that the probability of a collision event resulting in fragmentation P_frag can be fully expressed in terms of the relative velocity υ_rel (Birnstiel et al. 2011): $$P_{\text{frag}}=\{\begin{array}{ll}0\hfill & \text{if\hspace{0.17em}}{\upsilon }_{\text{rel}}<{\upsilon }_{\text{frag}}-\delta {\upsilon }_{\text{frag}},\hfill \\ 1\hfill & \text{if\hspace{0.17em}}{\upsilon }_{\text{rel}}\ge {\upsilon }_{\text{frag}},\hfill \\ 1-\frac{{\upsilon }_{\text{frag}}-{\upsilon }_{\text{rel}}}{\delta {\upsilon }_{\text{frag}}}\hfill & \text{otherwise.}\hfill \end{array}$$(40)

We here defined δυ_frag = υ_frag/5, where we follow Birnstiel et al. (2011). After the collision partner and P_frag have been determined, we generate a random number R₅ ∊ [0, 1]. A destructive collision will occur if R₅ ≤ P_frag.

Regardless of collision outcome, we need to re-calculate the home aggregate mass and size, m_a and s_a, after each collision event. In addition, different collision events have an effect on the depth at which our monomer is embedded, z_m. We distinguish between the following collision outcomes:

1. Coagulation. In this case, R₅ > P_frag and the collision event results in the coagulation of the home aggregate and the collision partner(s) of mass m_c. The new home aggregate mass m_a,new then follows from $$m_{\text{a},\text{new}}=m_{\text{a},\text{old}}+m_{\text{c}}{N}_{\text{col}}.$$(41)

Note that if m_c/m_a,old > ƒ_c, the group size is N_col = 1 (Eq. (37)). The new monomer depth z_m,new is calculated as $$Z_{\text{m,new}}=Z_{\text{m},\text{old}}+S_{\text{a},\text{new}}-S_{\text{a,old}},$$(42)

where s_a,old and s_a,new denote the old and new home aggregate size, respectively. Effectively we thus assume that the mass of the collision partner(s) is added to the old home aggregate as a homogeneous surface layer, such that the monomer is buried deeper inside the new home aggregate.

2. Fragmentation. In this case, the home aggregate is catas-trophically disrupted by the collision partner, which must be similar to the home aggregate in terms of mass: $$f_{\text{c}}<\frac{m_{\text{c}}}{m_{\text{a,old}}}<\frac{1}{f_{\text{c}}}.$$(43)

The new home aggregate of the monomer becomes one of the fragments produced during the collision. We assume these fragments follow the distribution (Birnstiel et al. 2010) $$n\left(m\right)\{\begin{array}{ll}\propto m^{-{\xi }_{\text{frag}}}\hfill & \text{if\hspace{0.17em}}{m}_{\text{m}}\le m\le m_{\max },\hfill \\ =0\hfill & \text{otherwise}\hfill \end{array}$$(44)

where ξ_frag = 1.83 (Brauer et al. 2008). In addition, m_max denotes the maximum possible fragment mass. In the case of fragmentation, we set m_max = max(m_a, m_c). We draw the new monomer depth z_m,new from a spherically uniform distribution, such that the probability P(z_m,ew) of the monomer being located at a certain depth after a destructive collision is proportional to the home aggregate mass density ρ_a. This means that the monomer is more likely to be located close to the surface after a fragmentation event since P(z_m) ∝ (s_a − z_m)². In reality, the functional form P(z_m,new) is likely set by the size of the original home aggregate s_a,old size of the collision partner s_c and their relative velocity υ_rel (Dominik & Tielens 1997).

3. Erosion. Erosion occurs when m_c/m_a,old ≤ ƒ_c, as the collision partner is too small to result in full fragmentation of the home aggregate. Instead, the collision partner exhumes and ejects mass of the home aggregate. We here follow Krijt & Ciesla (2016) and assume that the mass ejected is m_ej = 2m_c. In order to determine whether the monomer is ejected, we define the ejection probability P_ej $$P_{\text{ej}}=\frac{m_{\text{ej}}}{m_{a\text{,old}}},$$(45)

and generate another random number R₆ ∊ [0, 1], such that the monomer is ejected if R₆ ≤ P_ej.

• Ejection. In this case the new home aggregate is one of the fragments ejected, whose mass is selected according to Eq. (44). In this case, we set m_max = m_c, and also draw the new monomer depth z_m from a spherically uniform distribution.

• No ejection. If the monomer remains in its old home aggregate, the new lower home aggregate mass m_a,new and smaller monomer depth z_m,new are calculated via $$m_{a,\text{a,new}}=m_{\text{a,old}}-m_{\text{c}}{N}_{\text{col}}$$(46) $$Z_{\text{m,new}}=Z_{\text{m,old}}+S_{\text{a,new}}-S_{\text{a,old}}.$$(47)

  This means that the monomer depth decreases as mass is removed isotropically from the surface.

4. Impact. In this case m_c/m_a,old ≥ ƒ_c, which means that the home aggregate is impacting and eroding the considerably larger collision partner. We here assume that the ejected material only originates from the collision partner, such that the new home aggregate becomes the collision partner, whose mass is $$m_{\text{a},\text{new}}=m_{\text{c}}-m_{\text{a,old}}.$$(48)

Again, we assume that the collision results in major restructuring of the new home aggregate, which is why we again use a spherically uniform distribution to determine a new value for the monomer depth z_m.

We consider a test similar to the one performed by Krijt & Ciesla (2016). Since tracking a monomer effectively entails tracking a unit of mass, the fraction of time which monomers spend inside aggregates of a given size s_a, when summed over many individual monomer histories, must be proportional to the background local dust mass density as a function of grain³ size a = s_a, ρ_d(a, r, z). We compare the total integrated monomer residence time of 6000 monomers to the vFrag5 background model grain size distribution ρ_d(a, r, z) for monomers in aggregates with fragmentation velocities of υ_frag = 1, 5, 10m s⁻¹. All monomer positions are kept fixed at r = 10 AU and z/r = 0, and initially assumed to be bare, such that the initial home aggregate size is always s_a = s_m = 5 × 10⁻⁸ m. Furthermore, the aggregate density is assumed equal to the monomer density (i.e. ρ_a = ρ_m, corresponding to a filling factor ϕ = 1). Each simulation is run for 110 kyr, where we exclude the first 10 kyr of each simulation from the analysis such that we can neglect effects of our initial choice of home aggregate size. The timestep is set to ∆t = 10 yr, which is several orders of magnitude below the total group collision rate $1/{\tilde{C}}_{\text{tot}}$ at 10 AU (see Sects. 2.5 and 3.1).

The integrated monomer residence times shown in Fig. 4 follow the background density distribution for all three fragmentation velocities υ_frag quite well for small aggregates up to ~3 × 10⁻⁴ m. In addition to the integrated monomer time in Fig. 4, we also report the total number of occurrences of the different collision outcomes for all 6000 monomers in Table 2.

For all three fragmentation velocities, the residence time shows a small dip at ~3 × 10⁻⁵ m. This results from the stopping time associated with this aggregate size becoming larger than the turnover timescale of the smallest scale eddies t_η (Eq. (31)), resulting in a sudden increase in the relative velocity due to turbulence (Ormel & Cuzzi 2007), which increases the collision rates associated with aggregates of s_a ~ 3 × 10⁻⁵ m. Therefore, these aggregates are depleted in a similar fashion as the smallest aggregates. We note that both the stopping time and t_η depend on the local disk gas density ρ_g and soundspeed c_s. This means that the location of this dip depends on position in the disk.

Above s_a ~ 3 × 10⁻⁴ m, the integrated residence time starts to deviate significantly from the background grain size distribution for all three values of υ_frag. This is primarily due the fact that although our background grain size distribution incorporates settling, it does not directly account for the effects of collisions. For υ_frag = 1 m s⁻¹ and υ_frag = 10 m s⁻¹, the peak in the residence time occurs respectively for aggregate sizes smaller or larger than the maximum grain size of the background grain distribution. This is expected since the maximum grain size in the background model is calculated for υ_frag = 5 m s⁻¹. For υ_frag = 1 m s⁻¹, aggregates larger than ~3 × 10⁻⁴ m should not be able to form (Fig. 2). This size limit is indeed the location of the peak in monomer residence time in Fig. 4. However, the maximum grain size at 10 AU in the vFгag5-model is a_max ≈ 1 × 10⁻² m, which means that aggregates are still colliding with grains above the maximum aggregate size allowed by the collision model. Such a collision would occur at a relative velocity higher than υ_frag = 1 m s⁻¹, although it is still possible for the monomer to end up in the larger collision partner via an impact event. Table 2 shows that for the monomers in aggregates with υ_frag = 1 m s⁻¹, impact events occur frequently. This can explain why such monomers still spend a considerable amount of time in aggregates with s_a ≳ 3 × 10⁻⁴ m. Such aggregates would have a short lifetime, as in this size regime only destructive collisions are possible. This also explains why the monomers in aggregates with υ_frag = 1 m s⁻¹ undergo significantly more fragmentation, erosion, and ejection events than monomers in aggregates with υ_frag = 5, 10 m s⁻¹ (Table 2).

For monomers in aggregates with υ_frag = 10 m s⁻¹, the aggregate is allowed to grow to sizes beyond the local value of a_max ≈ 1 × 10⁻² m in the background model. In this case, there are no grains of similar size present in the background model at r = 10 AU and z = 0, which means that the aggregate can continue to grow, primarily through the sweep-up of many groups of smaller grains. Similarly, when relative velocities start to exceed υ_frag = 10 m s⁻¹, aggregates have grown to such an extent (up to s_a ~ 4 cm) that collisions with similarly-sized grains almost never occur. Table 2 indeed shows no fragmentation events for monomers in aggregates with υ_frag = 10 m s⁻¹, indicating that only erosion by many groups of small grains is a viable mechanism to limit further growth of the aggregate. However, Fig. 4 shows that this does not prevent monomers from spending a significant amount of time in aggregates larger than a_max.

For the case υ_frag = 5 m s⁻¹, which is consistent with the value for a_max calculated in the background model, the integrated residence time and ρ_d are in good agreement for almost all grain sizes. However, some differences between the integrated monomer residence time and ρ_d(α) still remain. The former shows a deficit with respect to ρ_d(α) in the size regime below a_max = 1 × 10⁻² m, and also shows a small excess of time where monomers reside in aggregates larger than the maximum grain size in the background model a_max. In the collision model, the aggregate size is allowed to peak at sizes for which υ_rel > υ_frag. However, these aggregates are very short-lived, as collisions tend to result in fragmentation or erosion of such a large aggregate. This may also explain the significant number of fragmentation, ejection, and erosion events in Table 2 for υ_frag = 5 m s⁻¹. The formation of aggregates larger than a_max may come at the cost of the number of aggregates between s_a ~ 3 × 10−⁴ m and a_max, as these are preferentially growing to aggregates slightly above a_max with respect to the background model. Altogether these discrepancies appear to result from the fact that the grain size distribution in the background model is approximated with a power law distribution (Eq. (5)), modified by the effects of settling. Furthermore, the grain density drops to zero in the first size bin above a_max(r). However, this may be incorrect near the peak of the grain size distribution, where a dust size distribution limited by turbulence-driven fragmentation generally departs from a power law such as the one assumed in the background model (Birnstiel et al. 2011; Krijt & Ciesla 2016).

Altogether our collisional model is able to reproduce the background dust size distribution which is most consistent with the fragmentation velocity used in our model. As illustrated above, the usage of a background dust size distribution which is not automatically consistent with the collision model can induce errors in the collisional histories of individual monomers. This occurs when there is a difference between the fragmentation velocity in the collision model and the fragmentation velocity which determines the maximum grain size a_max at a given position in the background model (Eq. (6)). The key variable affected is the home aggregate size s_a, specifically the average time a monomer spends in home aggregates of a given size. Furthermore, discrepancies between the value for υ_frag used in the collision model and the value used to calculate a_max in the background model also affect the frequencies of different collision outcomes. Both of these errors propagate into the monomer depth z_m, which in turn can affect the ice evolution of the monomer (see also Sect. 4.2). Therefore, it is key to ensure that the background model dust size distribution ρ_d(a, r, z) is as consistent as possible with the assumptions in the collision model, such as the fragmentation velocity υ_frag.

Fig. 4 Integrated monomer residence time as a function of home aggregate size for different fragmentation velocities (solid coloured lines). The solid black lines on the left and right represent the minimum and maximum grain size occurring throughout the entire vFrag5 background disk model. The dashed lines all denote the same background dust mass distribution as a function of size, ρd(a) at r = 10 AU and z = 0, shifted vertically for comparison with the integrated monomer residence time.

2.4 Ice evolution

The ice evolution model tracks the amount of ices of different species present in the ice mantle around the monomer. We here only track the ice composition of the monomer as it travels through the disk and is constantly changing home aggregate. Therefore, we aim to develop a stochastic approach for how the fluxes of volatile molecules and photons impinging onto the home aggregate affect the ices associated with the monomer. A crucial parameter in regards to both is how deep the monomer is embedded inside the home aggregate, z_m (see Fig. 1). A monomer close to the surface (z_m ~ s_m) can more easily acquire an ice mantle from impinging molecules than a monomer buried beneath many layers of other monomers (z_m ~ s_a ≫ s_m). However, such a buried monomer would also be more protected from photodesorption by the overlying monomers, which would absorb any incident UV radiation. If we assume that the local disk environment allows molecules to approach the home aggregate in straight lines (i.e. if the aggregate is in the Epstein regime), the monomers at z < z_m can be treated as a slab of particles each of size s_m, absorbing individual impinging molecules and photons. From a radiative transfer viewpoint, each geometrical monomer depth z_m can be associated with an optical depth τ defined as $$\tau ={\displaystyle \underset{0}{\overset{z_{\text{m}}}{\int }}\sigma n_{\text{m}}\text{d}z.}$$(49)

Here, σ denotes the collision cross section, and n_m the monomer number density. We approximate the collision cross section with the monomer cross section $\sigma \approx {\sigma }_{\text{m}}=\pi s_{\text{m}}^2$, where we neglect the molecular cross section in case of impinging molecules. However, we note that this may be inappropriate for large molecular complexes which are similar in size to the monomer. Assuming constant monomer density throughout the home aggregate, we can estimate n_m as $$n_{\text{m}}=\frac{\varphi {\rho }_{\text{m}}}{m_{\text{m}}}=\frac{3\varphi }{4\pi s_{\text{m}}^3},$$(50)

where ϕ denotes the home aggregate mass filling factor and $m_{\text{m}}=\frac{4}{3}\pi s_m^3{\rho }_{\text{m}}$ denotes the monomer mass, such that Eq. (49) becomes $$\tau =\frac{3}{4}\frac{z_{\text{m}}}{s_{\text{m}}}\varphi .$$(51)

For an isotropic radiation field, the fraction of radiation which can reach down in the slab to an optical depth τ can be written as (e.g. Whitworth 1975) $$E_2\left(\tau \right)=\tau {\displaystyle \underset{\tau }{\overset{\infty }{\int }}\frac{e^{-t}}{t^2}\text{d}t.}$$(52)

Here, E₂ denotes the exponential integral, and represents the spatially averaged fraction of all impinging molecules and photons which can reach down to depth z_m in the home aggregate. Note that the above equations strictly apply when s_a ≫ s_m, and that the home aggregate has a fractal dimension D_f = 3, such that $m_{\text{a}}\propto s_{\text{a}}^3$.

Although E₂ represents the fraction of molecules and photons which can reach down to a given optical depth τ, the actual flux of molecules and photons received by an individual monomer does not have to scale with E₂(τ). Monomers may be arranged in fractal and clustered configurations inside the home aggregate as a consequence of initial hit-and-stick growth followed by compaction (e.g. Dominik & Tielens 1997; Wurm & Blum 1998; Weidling et al. 2009). This means that some monomers at depth z_m may lie completely in the shadow of other monomers, whereas other monomers at depth z_m may be fully exposed as there are no other monomers blocking incoming molecules and photons that travel inside the home aggregate effective radius s_m. Although full modeling of this spatial flux distribution is beyond the scope of this paper, we account for this behaviour by interpreting Eq. (52) as an exposure probability P_exp: the probability that a monomer, located at a depth z_m, is exposed to the outside environment. Being exposed as a monomer here means that the monomer ice mantle is able to undergo adsorption or photodesorption. Altogether we define the exposure probability as $$P_{\exp }=\{\begin{array}{ll}1\hfill & \text{if\hspace{0.17em}}{z}_{\text{m}}<z_{\text{crit}}\hfill \\ E_2\left(\tau \right);\hfill & \text{if\hspace{0.17em}}{z}_{\text{m}}\ge z_{\text{crit}},\hfill \\ \tau =\frac{3}{4}\frac{z_{\text{crit}}-z_{\text{m}}}{s_{\text{m}}}\varphi .\hfill & \hfill \end{array}$$(53)

We here define the critical monomer depth z_crit = 2s_m = 10⁻⁷ m to ensure that monomers which are closer to the surface than the diameter of a single monomer are always exposed, regardless of filling factor ϕ. Fig. 5 shows the value of P_exp as a function of z_m for different values of ϕ. The filling factor is here varied from ϕ = 1 to ϕ = 10⁻⁵. Such a low value for ϕ is thought to be possible for aggregates grown through efficient sticking and inefficient fragmentation, which results in very fluffy aggregates (Okuzumi et al. 2012; Kataoka et al. 2013). For a filling factor of 10⁻⁵, it appears that molecules and photons have a significant chance to reach monomers which are embedded in the center of even the largest aggregates, amounting to more than 30% for s_a = 10⁻³ m. For larger filling factors, monomers in the center of large aggregates are mostly shielded from gas phase molecules and photons.

To model the evolution of the composition of the ice mantle of the tracked monomer we solve at each time step the following mass balance equation for each chemical species x to find the change in ice mass of that species M_x: $$\frac{\text{d}{M}_x}{\text{d}t}=\{\begin{array}{ll}4\pi s_{\text{m}}^2\left(\frac{1}{2}{ℛ}_{\text{ads},x}-{ℛ}_{\text{tds},x}-\frac{1}{2}{ℛ}_{\text{pds,}x}\right)\hfill & \text{if\hspace{0.17em}}{R}_7\le P_{\exp },\hfill \\ 4\pi s_{\text{m}}^2{ℛ}_{\text{tds,}x}\hfill & \text{otherwise.}\hfill \end{array}$$(54)

Here, R₇ ∊ [0, 1] denotes a random number from a uniform distribution we generate after each collision event to determine whether the monomer is exposed. Furthermore, ℛ_ads,x denotes the specific adsorption rate (in kg m⁻² s⁻¹), ℛ_ads,x the specific thermal desorption rate, and ℛ_pds,x the specific photodesorption rate associated with species x. The factors $\frac{1}{2}$ in Eq. (54) associated with ℛ_ads,x and ℛ_pds,x follow from the assumption that on average, impinging molecules are only effectively adsorbed on the monomer hemisphere aimed towards the aggregate surface, whereas the same applies for photodesorption by impinging UV photons. We treat the ice transport processes as follows:

• Adsorption. Molecules in the gas phase can collide and stick to the monomer, such that the specific adsorption rate (in kg m⁻² s⁻¹) for species x is given by the product of the molecule collision rate and sticking probability: $${ℛ}_{\text{ads,}x}=n_x{m}_x{\upsilon }_{\text{th,}x}S.$$(55)

  Here, n_x denotes the gas phase molecule number density (in m⁻³), mx the molecular mass of species x. υ_th,x denotes the mean thermal speed of species x and is calculated as (Woitke et al. 2009) $${\upsilon }_{\text{th,}x}=\sqrt{\frac{k_{\text{B}}{T}_{\text{g}}}{2\pi m_x}.}$$(56)

  S denotes the sticking probability. We calculate S using the expression of He et al. (2016) $$S={\alpha }_S\left[1-\tanh \left({\beta }_S\left[T_{\text{d}}-{\gamma }_S\frac{E_{\text{ads,}x}}{k_{\text{B}}}\right]\right)\right].$$(57)

  Here, α_S, β_S, and γ_S are fitting parameters, which have been derived for experimental data on the sticking of H₂, D₂, N₂, O₂, CO, CH₄, and CO₂ on amorphous, nonporous water ice (He et al. 2016). We note that α_S = 0.5 for all species to ensure 0 ≤ S ≤ 1, while the values of β_S and γ_S depend on molecule species. However, experimental data on β_S and γ_S for the sticking of different chemical species on different surfaces is limited. Therefore, we follow He et al. (2016) and use β_S = 0.11 K⁻¹ and γ_S = 0.042 for all species as an approximation.

• Thermal desorption. We express the thermal desorption rate k_tds,x for species x as $$k_{\text{tds,}x}=v_x\exp \left(-\frac{E_{\text{ads,}x}}{k_{\text{B}}{T}_{\text{d}}}\right).$$(58)

  Here, ν_x denotes the ice lattice vibration frequency (Tielens & Allamandola 1987; Cuppen et al. 2017): $${\upsilon }_x=\sqrt{\frac{2N_{\text{ads}}{E}_{\text{ads,}x}}{{\pi }^2{m}_x}.}$$(59)

  N_ads = 10¹⁹ m⁻² denotes the density of adsorption sites on the monomer surface. Furthermore, E_ads,x denotes the adsorption energy associated with species x. The specific thermal desorption rate for an ice species x can be expressed as (Cuppen et al. 2017) $${ℛ}_{\text{tds,}x}=k_{\text{tds,}x}{N}_{\text{act}}{f}_x{N}_{\text{ads}}.$$(60)

  Here, N_act is the number of actively desorbing ice molecule monolayers. Typical values for N_act range from 2 to 4 (Cuppen et al. 2017), which motivated using a median value of 3 in this work. We express the ice fraction ƒ_x in terms of mass: $$f_x=\frac{M_x}{m_x{\displaystyle \sum _y{M}_y/m_y}}.$$(61)

  Here, M_x denotes the total mass of species x present in the monomer ice mantle.

• Photodesorption. We treat photodesorption in a similar fashion as Woitke et al. (2009), and write the specific photodesorption rate as $${ℛ}_{\text{pds},x}=Y{f}_x\chi F_{\text{Draine}}{m}_x.$$(62)

  Here Y = 1 × 10⁻³ denotes the UV photon yield, and χ^F_Draine denotes the local UV radiation field strength in photons m⁻² s⁻¹.

In the above treatment, the escape of molecules from the home aggregate is not modelled. However, molecules liberated by thermal desorption from the ice mantles of monomers located at large τ collide with a significant number of other monomers before escaping back into the gas phase, which can result in re-adsorption if the sticking probability S > 0. We explore the timescale associated with molecule escape and the concept of re-adsorption in Appendix B. We find that re-adsorption may be important for particular molecular species in a comparatively small region immediately behind the iceline of the respective molecular species. Thus, our model could underestimate the amount of ice from individual species retained on monomers at large τ in these disk regions.

As a numerical test for our model, we compare the trajectories of pure H₂O monomer particles predicted by our model to the model of Piso et al. (2015). In our comparison, we consider particles which are allowed to undergo radial drift and thermal desorption. We here use the background disk structure and model parameters used by Piso et al. (2015).

The disk gas surface density is calculated as $${\sum }_{\text{P}}\left(r\right)=20000{\left(\frac{r}{\text{AU}}\right)}^{-1}\text{\hspace{0.17em}\hspace{0.17em}kg\hspace{0.17em}}{\text{m}}^{-2},$$(63)

such that the midplane gas density becomes, when assuming a Gaussian vertical density profile $${\rho }_{\text{g,P}}=\frac{{\sum }_{\text{P}}\left(r\right)}{H_{\text{g}}\sqrt{2\pi }}.$$(64)

H_g in this context is calculated as the pressure scale height H_g = c_s/Ω. The soundspeed is obtained from the temperature profile (Piso et al. 2015) $$T_{\text{g,P}}\left(r\right)=120{\left(\frac{r}{\text{AU}}\right)}^{-\frac{3}{7}}K.$$(65)

Note that since collisions, adsorption and photodesorption are not considered by Piso et al. (2015), the evolution of the particle is fully determined by the local gas density and temperature. Since the gas and dust temperature T_g are coupled in the midplane, we assume T_g,P = T_d for this test. For radial drift, the dimensionless pressure gradient η (Eq. (14)) is estimated as $${\eta }_{\text{P}}=\frac{1}{2}{\left(\frac{c_{\text{s}}}{r\text{Ω}}\right)}^2.$$(66)

The sound speed c_s is calculated via $$c_{\text{s}}=\sqrt{\frac{k_{\text{B}}{T}_{\text{g,P}}}{{\mu }_{\text{P}}{m}_{\text{p}}},}$$(67)

where the dimensionless mean molecular weight is set to µ_P = 2.35 (Piso et al. 2015).

We also set the value of the adsorption energy for H₂O, $E_{\text{ads,}}{\text{H}}_2\text{O}$ to 5800 K, the same value as used by Piso et al. (2015). As a default value, we use the UMIST RATE12 value of 4800 K (McElroy et al. 2013). Furthermore, the particles have a density ρ_a = 2000 kg m⁻³, and a filling factor ϕ = 1. We calculate the ice lattice vibration frequency as (Hollenbach et al. 2009; Piso et al. 2015) $${\nu }_{{\text{H}}_2\text{O}}=1.6\times {10}^{11}\sqrt{\frac{E_{\text{ads},{\text{H}}_2\text{O}}{m}_{\text{p}}}{k_{\text{B}}{m}_{{\text{H}}_2\text{O}}}}{\text{s}}^{-1},$$(68)

while we calculate the specific thermal desorption rate ℛ_tds using Eq. (60) with $f_{{\text{H}}_2\text{O}}=1$^.

Figure 6 depicts the evolution of the monomer size as a function of time and H₂O ice mass as a function of radial distance r for pure H₂O grains with sizes of 10⁻³ m up to 10 m. We here compare the evolutionary trajectories predicted by our model with the trajectories depicted in the center row of Fig. 4 of Piso et al. (2015), where the same particles are released in the disk midplane at 10 AU. Note that for this test, the effects of turbulent diffusion on radial migration is neglected, such that aerodynamic drag is the only process affecting the radial position of particles.

Both the model of this work and the model of Piso et al. (2015) predict that regardless of initial particle size, the sublimation is almost instantaneous once a particle crosses the snowline. However, there appears to be a small difference in the time at which sublimation happens. For example, the model of Piso et al. (2015) predicts sublimation of a particle of 10⁻³m after 3 × 10⁶ yr, whereas the sublimation time in our model is 4 × 10⁶ yr. However, a more detailed analysis revealed that the most likely cause of this deviation is inaccuracies in the digitation software used to extract the data points from Fig. 4 of Piso et al. (2015). This is further strengthened by the agreement of the trajectories for the ice masses as a function of radial distance shown in the right panel of Fig. 6. The authors also note the high sensitivity of the sublimation time to changes in model parameters.

Fig. 5 Exposure probability as a function of monomer depth zm for different filling factors ϕ.

Fig. 6 Particle size as a function of time (left panel) and ice mass (scaled with the initial ice mass) as a function of radial distance (right panel) for pure H2O particles undergoing radial drift and thermal sublimation. Different colors denote different particles sizes. The trajectories predicted by the model of this work (solid lines) are compared with the trajectories calculated by Piso et al. (2015; dots). The vertical dashed line in the right panel denotes the position of the H2O iceline at 0.67 AU as calculated by Piso et al. (2015).

2.5 Timescales

To enable us to spatially constrain the regions where disk processes are coupled, we will compare their associated timescales later in this work (Sect. 3.1). In this subsection we outline our approach to estimating these timescales.

For dynamical processes, we compare the radial drift, vertical settling, and turbulent mixing timescales, denoted by τ_rm, τ_zm, and τ_tm, respectively: $$\begin{array}{ccc}{\tau }_{\text{rm}}=\frac{H_{\text{g}}}{{\upsilon }_r},& {\tau }_{\text{zm}}=\frac{H_{\text{g}}}{{\upsilon }_z},& {\tau }_{\text{tm}}=\frac{H_{\text{g}}}{\alpha c_{\text{s}}}\end{array}.$$(69)

Here, υ_r and υ_z denote the drift velocities due to aerodynamic drag, given by Eqs. (14) and (15), respectively, while ν_turb = αc_sH_g. We here choose the gas pressure scale height H_g as the length scale for the migration timescales (τ_rm, τ_zm) since H_g is the typical length scale associated with the turbulent mixing timescale τ_tm. For collision processes, the timescale over which a particle collides with a significant fraction of its own mass worth of collision partners can be estimated as $${\tau }_{\text{col}}=\frac{1}{{\tilde{C}}_{\text{tot}}}.$$(70)

Here the total group collision rate ${\tilde{C}}_{\text{tot}}$ is given by Eq. (37). We note that this definition differs from earlier work, where the timescale for collisional processing has been defined in terms of the time required for a particle to significantly change size as a function of the local dust-to-gas mass ratio (see e.g. Brauer et al. 2008; Birnstiel et al. 2012).

For ice formation, we consider the ice adsorption and desorption timescale for molecular species x as a function of position. We here estimate the average amount of ice present on a monomer from the solid phase number density predicted by the background model, n_x,ice, as (Figs. A.5–A.8) $$M_{x,\text{ice}}=\frac{n_{x,\text{ice}}{m}_x}{{\rho }_{\text{d}}}{m}_{\text{m}}.$$(71)

The adsorption and desorption timescales, τ_ads and τ_des, respectively, are then given by $${\tau }_{\text{ads}}=\frac{M_{x,\text{ice}}}{4\pi s_{\text{m}}^2{ℛ}_{\text{ads,}x}}$$(72) $${\tau }_{\text{des}}=\frac{M_{x,\text{ice}}}{4\pi s_{\text{m}}^2\left({ℛ}_{\text{tds,}x}+{ℛ}_{\text{pds,}x}\right)}.$$(73)

3 Results

As a first step, we use our model to constrain the disk regions where the dynamical, collisional and ice processes are highly coupled (Sect. 3.1). Subsequently, we demonstrate how the coupled behaviour of these processes affects the ice mantle composition of individual monomers (Sect. 3.2). For all monomer models, the parameter values denoted in Table 3 are used, unless noted otherwise.

3.1 Coupling analysis

As a starting point, we investigate the relative importance of the various grain transport processes. We explore the relative importance of turbulent stirring with respect to aerodynamic drift processes as a function of grain size in Fig. 7. Throughout this section, we consider the const-model as the background disk model for our analysis, which means that the maximum grain size a_max = 3 × 10⁻³ m remains constant as a function of r. Furthermore, we always assume a filling factor ϕ = 1, such that the material density is always equal to the monomer density, ρ_a = ρ_m = 2094 kg m⁻³. Figure 7 highlights the physical grain size associated with St = 10⁻³ for the const background disk model discussed in Sect. 2.1. We here inverted Eqs. (10) and (11) to give the grain size as a function of Stokes number. Grains with St = 10⁻³ would undergo significant settling (Eq. (8)) for the turbulence strength parameter value used, α = 10⁻³. It becomes clear that the physical grain size associated with a fixed Stokes number varies greatly throughout the disk, ranging from ~10⁻² m at the transition from the Epstein to the Stokes regime to ≲10⁻⁶ m in the outer disk and above z/r ~ 0.25. This is primarily a consequence of the background model gas density, which varies many orders of magnitude throughout the disk (Fig. A.1). A lower gas density decreases the coupling of dust grains to the gas, and hence increases the Stokes number of a dust grain of given size (Eqs. (10) and (11)). We note that the grain size associated with St = 10⁻³ can also vary as a function of gas temperature due to the dependence of the Stokes number on the soundspeed. However, this effect is much less pronounced in Fig. 7 due to two effects. The dependence of the Stokes number on the gas temperature is weaker than its dependence on the gas density ($\text{St}\propto 1/c_s\propto 1/\sqrt{T_{\text{g}}}$ compared to St ∝ 1/ρ_g). Furthermore, the variation in the gas temperature is small compared to the variation in the gas density, both radially and vertically (Fig. A.1). Despite the limited influence of the gas temperature, the dynamical behaviour of a dust grain of given size and its associated timescales can be expected to vary significantly as a function of position in the disk.

As a next step, we investigate the coupling behaviour of the timescales introduced in Sect. 2.5. As an example, we consider the radial behaviour of the timescales for radial drift, vertical settling, turbulent stirring, collisional growth, ice adsorption, and ice desorption in Fig. 8. The timescales are associated with a dust grain of 10−⁵ m with an H₂O ice mantle. The vertical height is kept fixed at z/r = 0.05.

We find that in the inner disk (r ≲ 1 AU), the collisional growth timescale τ_col becomes shorter than any dynamical timescale by at least an order magnitude, and drops even below ~10⁻¹ yr at r = 0.1 AU. The shorter collision timescale for decreasing r is a result of the increasing dust density closer to the star, which in turn increases collision rates (Eq. (22)).

The turbulent stirring timescale τ_tm, although 1–3 orders of magnitude longer than τ_col, is still far shorter than the aerodynamic drag timescales, τ_rm and τ_zm, which are between ~10⁴ and ~10⁶ yr for r ≲ 1 AU. This is a consequence of the rather small particle size in combination with a high gas density in the inner disk, which together result in a small Stokes number (≪10⁻³ Fig. 7), and hence low drift velocities.

For r ≳ 1 AU, almost all timescales initially increase as a function of r. In case of the dynamical timescales, the increase results from the increase of gas scale height H_g (Eq. (3)). However, beyond 10 AU, the lower gas density results in an increasing Stokes number, and therefore increasing drift velocities result in τ_rm and τ_zm gradually decreasing around r ≳ 100 AU, overtaking turbulent stirring as the fastest dynamical processes around r ~ 150 AU.

For r ≳ 1 AU, the ice formation timescales, τ_ads and τ_des both increase as a function of r. We define the H₂O iceline as the position where τ_ads = τ_des. Although Fig. 8 indicates that τ_ads and τ_des are very similar for r ≲ 4 AU, this condition puts the H₂O iceline at r ≈ 1.7 AU for z/r = 0.05. We note that both beyond the inner and outer boundary of the ice-forming region, τ_des decreases. Interior to the iceline this is due to increased thermal desorption, while in the outermost disk regions this is a consequence of higher photodesorption rates. The radial behaviour of τ_ads is primarily shaped by the local H₂O gas abundance, which explains its rapid decrease interior to the iceline: as H₂O ice sublimates interior to the iceline, the gas phase H₂O abundance increases significantly in the background model (Appendix A.2). The rapid increase of τ_ads for r < 1 AU is a consequence of the decrease in the sticking probability S.

In terms of coupling behaviour, it becomes clear that in the inner disk, the various disk processes are largely decoupled since τ_ads,τ_ads ≪ τ_col ≪ τ_tm ≪ τ_rm, τ_zm, whereas dust dynamics, collisions, and ice processing become more coupled for increasing r up to r ≈ 20 AU. Beyond this distance, both τ_ads and τ_des become much larger than the other timescales, indicating that ice processing is very slow compared to collisional and dynamical processing in the outer disk regions.

It should be emphasized that this result only applies to a dust grain of size s_m = 10⁻⁵ m with a pure H₂O mantle. Since the collisional growth timescale is much smaller than the disk lifetime throughout most of the disk, it is expected that these dust grains will likely be incorporated rapidly into larger dust aggregates. Therefore, it is necessary to explore the coupling and decoupling of the various processes as a function of grain size, ice species, and also for the full extent of height z/r above the disk midplane.

Figure 9 explores the coupling of the various timescales as a function of both r and z/r in case of a dust grain of s_m = 10⁻⁷ m (left panel) and s_m = 10⁻³ m (center panel). For s_m = 10⁻⁷ m, we also explore a case where we only allow a pure CO ice mantle to form instead of a pure H₂O ice mantle (right panel). We highlight the disk regions where the various timescales obey various decoupling conditions. We here categorized the timescales in a dynamical (τ_rm,τ_zm and τ_tm), collisional (τ_col), and ice processing τ_ads,τ_des group. The process with associated timescales τ_i then decouples from the other processes with associated timescales $$f_{\text{dec}}\cdot \max \left({\tau }_i\right)<\min \left({\tau }_j\right),$$(74)

or $$\min \left({\tau }_i\right)>f_{\text{dec}}\cdot \max \left({\tau }_j\right).$$(75)

Here we set f_dec = 10³ to encapsulate that we require a particular process to act much faster (Eq. (74)) or much slower (Eq. (75)) than the other disk processes. For ice processing, we only consider τ_ads interior to the iceline, and τ_des exterior to the iceline in this calculation.

In all three cases, ice desorption decouples from the other processes in any disk region radially interior to the ice forming region. This decoupling originates from the short desorption timescale, since thermal desorption and photodesorption are very efficient outside the ice forming region. We also find in the cases of H₂O, ice formation can decouple from the other disk processes due to the low adsorption rate in the ice forming region (most clearly visible in the center panel of Fig. 9). This is primarily a result of the low gas phase abundance of H₂O inside the ice forming region (Appendix A.2). The center panel of Fig. 9 shows that this effect is more pronounced for a large dust grain since more ice mass must be accumulated onto the grain surface to double the initial ice mantle mass (Eq. (72)).

In all three panels of Fig. 9 we discern a region towards the innermost disk midplane where the three groups of processes fully decouple. In this region, dynamical processes become much slower than collisional growth, which is due to the high collision rates in the inner disk (Fig. 8). The center panel of Fig. 9 also demonstrates that this region of full decoupling in the inner disk is smaller for larger dust grains. The shrinking of this region can be attributed to the fact that more collisions are required to result in a significant change in mass for a massive dust grain. Therefore, the effective collision rates become smaller, and hence τ_col increases.

The most important features of Fig. 9 are the white regions in all three panels around the ice forming regions. For various species, all the three process categories couple to one another. However, the location of this fully coupled region is dependent on the grain size and position of the iceline. Since different grain sizes are connected via coagulation and fragmentation, and the position of the iceline depends on molecular species, there is no ubiquitously decoupled regime in the planet-forming region. In order to understand the ice composition of planetesimals, we thus require a modeling approach such as SHAMPOO which treats dynamical, collisional, and ice processing in a fully coupled fashion.

Fig. 7 Grain size for which the local Stokes number is equal to 10−3, as a function of position in the const background disk model. The red curve encloses where particles of the size for which St = 10−3 are in the Stokes drag regime.

Fig. 8 Behaviour of the timescales associated with various disk processes as a function of radial distance r. Timescales shown are for a monomer of size sm = 10−5 m with a pure H2O ice mantle, located at z/r = 0.05. The dashed line at t = 108 yr indicates the upper limit on timescales relevant for planet formation.

Fig. 9 Process coupling behaviour for different grain sizes and ice mantle compositions. Colored areas indicate the decoupling of the different timescale groups. The filling of the coloured areas with solid lines or dots indicate the decoupling of a process category since it is much faster or slower, respectively. Black lines indicate icelines (τads = τdes).

Fig. 10 Evolution of the vertical positions (upper left panel), radial positions (upper center panel), temperatures (upper right panel), local UV fluxes (lower left), and local molecular number density in the gas phase (lower center) for monomers in home aggregates with different fragmentation velocity υfrag, released from r = 10 AU and z/r = 0. The black dashed line in the lower left panel indicates a UV radiation field strength of 1 FDraine. Linestyles in the lower center panel indicate number densities of different gas phase species.

3.2 Model demonstration and effects of different fragmentation velocity

As a first demonstration of the full SHAMPOO code, we consider the evolution of three monomers, initially not embedded in a home aggregate. The monomers are allowed to be incorporated in aggregates which have a fragmentation velocity of υ_frag = 1, 5, or 10 m s⁻¹. For each of these three scenarios, the appropriate background disk model is used (vFrag1, vFrag5, and vFrag10, respectively), while the monomers are released from z = 0 at r = 10 AU. The initial ice budget of the monomers is informed from the respective background disk model via Eq. (71). This means that the monomers start with a total amount of ice of m_ice = 0.0174m_m, 0.0176m_m, and 0.0113m_m in the vFrag1-, vFrag5-, and vFrag10-model, respectively. The three monomers are allowed to be processed for 10⁵ yr. All home aggregates are assumed to have a filling factor of ϕ = 0.1 for all three scenarios. However, in practice the filling factor can differ for each home aggregate. We therefore explore the effects of different filling factors in Sect. 3.3. The random number generator seed is kept fixed between all three scenarios, such that the fragmentation velocity and the associated background disk model are the only differences between the three simulations.

3.2.1 Dynamical and environmental evolution

Figure 10 shows the evolution of the vertical and radial position of the monomer, the associated monomer temperature and local UV radiation field. The monomer temperature is assumed to be equal to the local dust temperature. Despite the same initial monomer properties and model seed, the three monomers follow clearly distinct dynamical trajectories within the first 10 kyr.

It becomes clear that all monomers are in the dynamical regime dominated by turbulent diffusion (Fig. 7). In the vertical dimension, Fig. 10 shows that in all three cases the monomer makes several vertical excursions to |z/r| ~ 0.2, whereas radially, all three monomers diffuse significantly inwards, with all monomers reaching inside of r = 6 AU in 100 kyr. As the monomers gradually diffuse inwards, all three monomers experience increasing temperatures. However, the monomer in the vFrag1 -model remains cooler than the other monomers between 60 and 80 kyr as it remains around r = 8 AU in this period, while the monomers in the vFrag5 - and vFrag10-model have already diffused interior to r = 6 AU after 60 kyr. We also discern small peaks in the temperature evolutions of the monomers, which can be associated in time with the short-lived vertical excursions.

Figure 10 also reveals that the local UV radiation field strength is usually below 1 FDraine for all three monomers, although exceptions are short peaks in field strength, which can be again attributed to vertical excursions higher in the disk atmosphere. The monomer in the vFragl-model here experiences significantly less peaks in field strength compared to the vFrag5- and vFrag10-model. These latter background models have a thinner dark midplane region (Figs. A.2–A.4), and thus smaller vertical excursions are necessary for a monomer to reach the irradiated surface disk layers. Altogether monomers in the vFrag5- and vFrag10-model are more often in a regime where photodesorption can remove ice. The highest peaks in UV radiation field strength are found between 25 and 30 kyr for the monomer in the vFrag5-model and 40 and 25 kyr for the monomer in the vFragl-model, where the UV radiation field very briefly peaks at 10²F_Draine on both occasions. Although there exist many other peaks in UV radiation field strength, these peaks are usually well below 1 F_Draine and of short duration. Altogether we thus expect photodesorption to affect the ice evolution of monomers only for short periods of time if a monomer is exposed during a peak in the background UV radiation field.

3.2.2 Collisional history

The collisional histories of the three monomers are shown in Fig. 11, where we consider the evolution of the home aggregate size. Figure 11 also shows the type of interactions which result in a change in home aggregate sizes. We also count the number of collision outcomes per monomer history, giving rise to Table 4.

Figure 11 and Table 4 show that fewer collisions occur over 100 kyr for increasing υ_frag. The monomer in the vFrag10-model undergoes less than half as many collisions as the monomer in the vFrag1-model (530 collisions compared to 1165). This is a consequence of the total dust surface mass density Σ_d in the background disk model. In the background disk models with higher υ_frag, larger dust aggregates can form (see also Fig. 2), and thus the same amount of dust mass becomes concentrated in fewer dust particles, which in turn lowers collision rates. This is also confirmed by the typical home aggregate size seen in Fig. 11: s_a can be larger than 1 cm in the vFrag 10-model, while it is smaller than 1 mm in the vFrag1-model. We note that these home aggregate sizes are expected based on the maximum grain sizes allowed between r = 4 AU and r = 10 AU according to Fig. 2. Furthermore, for all three models, the frequency of collisions increases later in the evolution, which is a consequence of the radial diffusion inwards towards higher density disk regions (Fig. 10) for all three monomers.

Although there are thus significant differences between the collision histories of the three monomers, coagulation is the most frequent collision outcome for all three fragmentation velocities, amounting to almost 90% of the collision outcomes in all three cases. Aggregate fragmentation and erosion (comprising erosion, ejection, and impact events) appear to be rather rare compared to coagulation, in particular for the vFrag 5- and vFrag10-model. However, we note that coagulation is the only event allowed for collisions with υ_rel < υ_frag, which is true for most collisions in the monomer history. Once the home aggregate has grown to a size at which collisions with υ_rel > υ_frag become commonplace, destructive collisions usually result in a smaller home aggregate, placing it back in the regime where, υ_rel < υ_frag·

Figure 11 shows that the largest decreases in home aggregate size are associated with monomer ejection, undoing the work of many coagulation events. This is expected since in these events, the monomer we track ends up in a fragment produced by erosion of the previous home aggregate. These fragments are guaranteed to be significantly smaller than the previous home aggregate (Sect. 2.3).

Although coagulation can produce large sudden increases in home aggregate sizes (e.g. the first 5 kyr in the vFrag 10-model), we also discern many episodes in all three models where the home aggregate grows through the sweep-up of many smaller aggregates. We note that the behaviour exhibited in Fig. 11 is similar to the aggregate size history shown as an example in Krijt & Ciesla (2016).

Fig. 11 Collision histories and resulting home aggregate sizes sa as a function of time for monomers in aggregates of different fragmentation velocities υfrag. Note that the vertical position of each collision event coincides with the home aggregate size sa after the collision has occurred.

Fig. 12 Evolution of monomer depth zm for monomers embedded in home aggregates with υfrag = 1 m s−4 (left), 5 m s−1 (center), and 10 m s−1 (right). Green areas denote the times at which the monomer is exposed to gas phase molecules and UV photons, while the monomer is shielded from impinging molecules and UV photons otherwise. In each panel, the dashed line indicates zm = zcrit.

Fig. 13 Evolution of the total ice mass mice of monomers in aggregates with different fragmentation velocity υfrag (left panel), together with snapshots of the monomer ice composition at 10 kyr (center panel), and 100 kyr (right panel).

3.2.3 Ice evolution

The evolution of the amount of ice on monomers depends on the local abundances of gas phase molecules at a given time (lower center panel of Fig. 10). Furthermore, the evolution of the monomer ice mantle crucially depends on the depth at which the monomer is embedded in its home aggregate z_m, which is depicted for the three scenarios in Fig. 12. The monomer depth, which is a product of the collisional history, in combination with dynamical processes gives rise to the ice evolution depicted in Fig. 13. Since aggregates with lower fragmentation velocity remain smaller, it is more likely that the monomer is located close to the surface. Therefore, we find in Fig. 12 that the monomer in the vFrag1-model has the lowest average monomer depth compared to the vFrag5- and vFrag10-models. As a consequence, the monomer in this model is also exposed to the gas phase most often. In the vFrag10-model, the monomer only experiences one short initial period of exposure, which is mainly a result of our initial condition for the monomer to start without a home aggregate. For most of its evolution, the monomer in the vFrag10-model is located below z_m = 10⁻⁴ m. The vFrag5-model represents the intermediate case where the episodes of exposure after the initial phase occur less frequently than for the monomer in the vFragl-model.

Figure 13 shows the evolution of the total amount of ice on the three monomers, along with snapshots of the monomer composition at 10 kyr and 100 kyr. Although the monomer in the vFrag1-model is exposed to gas phase molecules most often, the left panel of Fig. 13 shows that after 100 kyr, the monomer has a similar total amount of ice (2.1 × 10⁻³m_m) as the monomer in the vFrag10-model (1.8 × 10⁻³m_m). Both models lost significant amounts of ice. In the vFrag1-model, this is the consequence of the monomer being exposed around 40 kyr, where the monomer makes a significant vertical excursion (Figs. 10 and 12), allowing all ice to evaporate in the strong UV radiation field via photodesorption. However, the frequent episodes of exposure do allow the monomer to regain a small amount of ice after the photodesorption around 40 kyr.

Although the monomer in the vFrag10-model is protected from photodesorption – as it is almost always unexposed, the left panel of Fig. 13 shows that it still loses a significant amount of ice between 70 kyr and 80 kyr; a period where it also experiences the highest temperatures (~90 K) of all three monomers. Moreover, the snapshot of ice composition at 100 kyr in the right panel of Fig. 13 shows that the monomer in the vFrag10-model has lost all its H₂O ice, while at 10 kyr, all three monomers have similar amounts of water (center panel of Fig. 13). This is striking since the H₂O iceline based on the adsorption and desorption timescales is located around r ≈ 2 AU around the disk midplane for all background models, a radial distance which is not reached by any of the monomers. Therefore, the adsorption rate is larger than the thermal desorption rate. However, since the monomer in the vFrag10-model is always unexposed, all the impinging molecules are adsorbed by monomers closer to the home aggregate surface. As the monomer we track is shielded from adsorption in this case, thermal desorption causes this monomer to lose all H₂O ice on its surface, despite the aggregate being located behind the H₂O ice line. However, this may be an overestimate of the amount of H₂O lost, as the H₂O lost through thermal desorption may in reality be partially replenished by the re-adsorption of H₂O originating from other nearby monomers (see also Appendix B).

Compared to the other two monomers, the evolutionary trajectory for the monomer in the vFrag5-model is less dramatic, with more moderate episodes of adsorption during exposure and desorption at 20 and 90 kyr. We note that this monomer is unexposed at the times of ice loss, suggesting that this loss occurs in a similar although less extreme fashion as for the vFrag10-model.

In terms of composition, the center and right panel of Fig. 13 show that at 10 kyr, the three monomers have ice mantles of similar composition. Here, H₂O is the dominant ice species with smaller amounts of NH₃ and CO₂. It becomes clear that in addition to the H₂O ice loss for the monomers in the vFrag1- and vFrag10-model, all three monomers have lost their CO₂ after 100 kyr. This is a result of the low adsorption energy of CO₂ compared to H₂O and NH₃ ($E_{{\text{ads,CO}}_2}/k_{\text{B}}=2990$ K compared to $E_{{\text{ads,H}}_2\text{O}}/k_{\text{B}}=4800$ K and $E_{{\text{ads,NH}}_3}/k_{\text{B}}=5534$ K; McElroy et al. 2013). As all three monomers gradually diffuse inwards towards warmer disk regions, we expect all CO₂ to be lost through thermal desorption. On the contrary, the amount of NH₃ after 100 kyr remains fairly constant for all three monomers. For the monomer in the vFrag5- and vFrag10-model, this can be explained through the high adsorption energy of NH₃, which means that the role of thermal desorption is limited in comparison to CO₂ and H₂O. We note that lower values for E_ads,NH₃ have been inferred (e.g. Sandford & Allamandola 1993; Brown & Bolina 2007), which would result in NH₃ being less stable against thermal desorption in this disk region. For the vFrag1-model, all NH₃ ice at 100 kyr must have accumulated after the photodesorption event at 40 kyr, indicating that significant NH₃ adsorption has occurred after this catastrophic ice loss.

Fig. 14 Evolution of the total amount of ice (left panel), H2O ice (center), and NH3 (right) averaged over 500 monomers for different home aggregate filling factors. Solid lines indicate the mean of the 500 monomers, whereas dotted lines indicate the standard deviation. The black dots indicate the ice evolution behaviour predicted by ProDiMo for the vFrag5-model. Trajectories of individual monomers are also shown in the background.

3.3 Comparison to background model and effects of filling factor

Throughout this work, the home aggregates were assumed to have a filling factor of ϕ = 0.1 or ϕ = 1, which is consistent with dust aggregates which have undergone some degree of compaction. Experiments have demonstrated that typical filling factors ϕ associated with compact aggregates range from ϕ ~ 0.2 (Langkowski et al. 2008) up to ϕ ≂ 0.4 (Weidling et al. 2009, 2012). Alternatively, modelling work has shown that large, fluffy aggregates with ϕ ≲ 10⁻³ can form if sticking is efficient and fragmentation is inefficient (Okuzumi et al. 2012; Kataoka et al. 2013). Since the value of ϕ affects how often a monomer is exposed and can accumulate ice, it is key to explore the sensitivity of ice formation in our model on ϕ. Furthermore, ice evolution is also considered in the time-dependent chemistry of the ProDiMo background models. Therefore, we can compare the ice evolution predicted for a local monomer population by SHAMPOO to the ice evolution calculated by ProDiMo.

For these purposes, we consider the evolution of the average amount of ice present on 500 monomers over 100 kyr in the vFrag5-model in Fig. 14. Four different scenarios are considered in Fig 14, where monomers are embedded in home aggregates with a constant filling factor of ϕ = 1, 10⁻¹, 10⁻² or 10⁻³. As we also aim to compare the results to the fully local evolution predicted by the time-dependent chemistry in ProDiMo over the same time period, the monomer positions are fixed at r = 10 AU in the disk midplane. This also means that the local disk environment remains the same for all simulations. Tests with dynamical processing enabled were also performed and yielded qualitatively similar results for the total amount of ice as shown in this section. However, the inclusion of dynamics results in a significant increase in the diversity of monomer trajectory and hence the standard deviation.

Initially, all monomers are assigned to a home aggregate of random size s_a which is drawn from the background model dust mass distribution ρ_d(a, r, z). Subsequently the monomer depth z_m is drawn from a spherically uniform probability distribution, such that the initial values of both s_a and z_m are weighted with the distribution of dust mass in the background model. All monomers start with 0.0176 times their bare monomer mass m_m worth of ice, which is calculated according to Eq. (71). The left panel of Fig. 14 shows that for all values of the home aggregate filling factor ϕ, monomers gain ice over the 100 kyr evolution, primarily as a consequence of the adsorption of H₂O and to a lesser extent, NH₃ and CO₂. The amount of ice accumulated, however, depends on the filling factor ϕ. Fig. 14 shows an increase of the total amount of ice of only 0.93% for monomers in aggregates with ϕ = 1, while monomers in very porous aggregates with ϕ = 10⁻³ on average gain ~20% their initial ice mass worth of ices over 100 kyr. The intermediate cases of ϕ = 10⁻² and ϕ = 10⁻¹ result in an average increase in the total ice mass of 6.1 and 2.0%. This decrease in the efficiency of ice adsorption for increasing ϕ is a result from the lower exposure probability P_exp at given z_m for higher filling factors. Figure 5 shows that for ϕ = 10⁻³, 10⁻², 10⁻¹, and 1, gas phase molecules can respectively reach down to depths of approximately 2 × 10⁻⁴ m, 2 × 10⁻⁵ m, 2 × 10⁻⁶ m, and 2 × 10⁻⁷ m. Therefore, monomers spending time in aggregates with higher ϕ spend less time exposed to the gas phase.

The center and right panels of Fig. 14 show that most of the ice accumulated during the 100 kyr evolution is H₂O and NH₃ for all values of ϕ. In all cases, these two species amount to over 90% of the molecules gained as ice, with the remainder being CO₂. This is a result from the fact that H₂O and NH₃ are the most prevalent molecules at 10 AU in the gas phase of our background model, with molecule number densities of $n_{{\text{H}}_2\text{O}}\approx 9.7\times {10}^5{\text{m}}^{-3}$ and $n_{{\text{NH}}_3}\approx 6.6\times {10}^5{\text{m}}^{-3}$ (Appendix A.2). Although CO and CH₄ are much more abundant in the gas phase (n_CO ≈ 2.4 × 10¹⁴ m⁻³ and $n_{{\text{CH}}_4}\approx 2.1\times {10}^{11}{\text{m}}^{-3}$), these molecules do not form any ices since temperatures at 10 AU are too high (see also Fig. 9).

Figure 14 also shows the average amount of ice expected on monomers based on the ice evolution predicted by the time-dependent chemistry of the vFrag5 ProDiMo model. We linearly interpolate between the default snapshot for n_x,ice at t = 2 × 10⁵ yr used in the background model and a new snapshot at t = 3 × 10⁵ yr obtained by evolving the chemistry in the vFrag5 ProDiMo model. The molecule number densities n_x,ice calculated in ProDiMo are here converted to units of ice mass per monomer with Eq. (71), which results in an increase of the total amount of ice of only 0.89% in the background model. This is roughly 10% less than the average total amount of ice expected on monomers in aggregates with ϕ = 1. We also find that in the cases of specific species, the average amount of H₂O and NH₃ on monomers in aggregates of ϕ = 1 calculated with SHAMPOO is ~10% above the increase calculated by ProDiMo. It is expected that monomers in aggregates with ϕ = 1 are in best agreement with ProDiMo, as ProDiMo chemically treats the refractory component of individual dust particles as rigid balls where molecules do not penetrate the aggregate surface upon adsorption. In ProDiMo, ice will only build up in a homogeneous layer on the aggregate surface, such that the rate at which ice of species x is adsorbed on a single dust aggregate is given by $$4\pi s_{\text{a}}^2{ℛ}_{\text{ads,}x}.$$

In SHAMPOO, the effective adsorption rate for the same aggregate can be estimated as $$N_{\exp }\cdot 2\pi s_{\text{m}}^2{ℛ}_{\text{ads,}x},$$

with N_exp being the total number of exposed monomers in the aggregate. In practice, this means that a given aggregate has on average more surface area available for adsorption than a dust particle of the same size in ProDiMo, since monomers below the aggregate surface can have have nonzero exposure probability, allowing for very large N_exp if the filling factor is low. Even for ϕ = 1, Fig. 5 shows that molecules can reach a few monomer radii below the aggregate surface, although the exposure probability drops off very quickly. This nonzero exposure probability down to a few monomer radii overestimates the surface area available for adsorption with respect to ProDiMo, and thus explains ~10% difference.

4 Discussion

4.1 Applications

The SHAMPOO code provides a tool which allows the investigation of the systematic behaviour of dust monomers under the non-local disk processing of a given protoplanetary disk structure. In this section we highlight a number of possible applications of the SHAMPOO code.

The primary goal of SHAMPOO is to reliably predict the CHNOS abundances of dust during the first stages of planet formation, where no planetesimals have formed yet. In case of the Solar System, this would encompass the first 100–500 kyr after the formation CAIs (e.g. Kruijer et al. 2014, 2020). In order to infer local CHNOS abundances in local dust populations, many monomer trajectories would be required. Within the framework of the background disk model, each monomer can be assigned to one spatial grid cells of the background disk model. This means that at any moment in time, every monomer is representing a unit of mass from that spatial grid cell. Predictions for the local CHNOS abundances can then be made by averaging over the individual CHNOS abundances of the monomers associated with that grid cell. A crucial requirement is that the spatial grid resolution is high enough such that the physical, thermal and chemical gradients within a grid cell are small.

Although this work has primarily focussed on the CHNOS-bearing molecules condensed on the monomer as ices, every monomer also consists of a refractory component of mass m_m. This refractory component can also contain a significant fraction of the total CHNOS elemental budget (e.g. Öberg & Bergin 2021). Currently, SHAMPOO does not assign a specific composition to the refractory component, although this could be incorporated in the model for example by assigning a static composition to the refractory monomer component. This refractory composition can be informed from the conditions at which the refractory component of the monomer has condensed. Altogether this would allow for a complete picture of the evolution of local CHNOS abundances.

We note a number of possible applications of SHAMPOO outside its main application. Specifically, SHAMPOO could be used to study the effect of dynamical and collisional processes on the transport of physically or chemically altered minerals throughout the planet-forming disk. Examples here include the possible outward diffusion of processed refractory material such as crystalline sillicates (e.g. Bockelée-Morvan et al. 2002; Watson et al. 2009; Olofsson et al. 2009; Williams et al. 2020), and the transport of hydrated minerals throughout the planet-forming disk (e.g. Brearley 2006; D’Angelo et al. 2019). However, such an application would involve connecting SHAMPOO to models which describe how these processes alter a single monomer given a certain disk environment. Another application lies in predicting the composition of pebble fluxes over time (e.g. Lambrechts & Johansen 2014; Morbidelli et al. 2015). The growth of pebbles involves the compaction of aggregates (i.e. ϕ increases over time) which gradually lowers the sticking probability after a collision. Therefore, such a compact aggregate would ultimately collisionally decouple from the growth-fragmentation cycles shown in Fig. 11 due to the bouncing barrier (Güttler et al. 2010; Zsom et al. 2010; Krijt et al. 2018). The application of SHAMPOO in this context would thus require the inclusion of porosity evolution. Furthermore, an application of SHAMPOO could lie in tracking the effects of non-local disk processing on local isotopic ratios such as D/H, ¹³C/¹²C, ¹⁵N/¹⁴N, and ¹⁸O/¹⁶O (see e.g. Cleeves et al. 2014; Visser et al. 2018; Öberg & Bergin 2021). This would, however, require the incorporation of the isotopologues in the background model chemical network and possibly also the feedback of sublimation into the background gas phase composition.

4.2 Model limitations

In this section we discuss the fundamental limitations of the SHAMPOO code, and highlight how we expect that these limitations are likely to affect predictions made with the SHAMPOO code both in this work and future work.

4.2.1 Feedback of monomer evolution on the disk

The code presented in this work offers a comprehensive description of how disk processes affect the volatile CHNOS abundances of monomer particles. However, the local conditions which affect the monomers are completely set by the background model, and are independent of the evolutionary trajectories of monomers. We highlight a number of effects this assumption may have on the evolutionary trajectories predicted by SHAMPOO.

In Fig. 4 we highlighted the possibility of inconsistencies between the local dust size distribution ρ_d(a, r, z) in the background disk model (Eq. (5)) and the in-situ dust size distribution inferred from the residence time of the monomer in home aggregates of given size. However, these two distributions are not independent since the background dust size distribution determines the collision rates (Eq. (22)). In Fig. 4 we demonstrated that given enough time, the time spent by the monomer in a home aggregate of given size approximately follows the residence time one would expect from the background distribution. However, we also noted significant departures, in particular when the dust size distribution of the background model is not consistent with the steady-state dust size distribution associated with a given fragmentation velocity υ_frag in the collision model (Sect. 2.3). In the background model dust size distribution, the maximum size is specified directly by a_max, whereas in the collision model, the maximum home aggregate size also depends on the fragmentation velocity υ_frag. In reality, however, a_max is a function of the fragmentation velocity. In this work, we accounted for inconsistencies between the maximum aggregate size in the collision model and a_max by using the parametrization from Birnstiel et al. (2012) for a local, steady-state distribution if turbulence is the main source of relative velocity between dust particles. However, a_max does remain a quantity which is not calculated self-consistently and has to be estimated separately from the SHAMPOO-code in the background model.

As a consequence of this discrepancy between the imposed background distribution and the steady-state distribution expected from the collision model, we found in Fig. 11 that monomers in aggregates with too low fragmentation velocity in SHAMPOO (i.e. the background model contains too many large grains) are frequently forced into aggregates which are larger than possible given the value of υ_frag. In these larger aggregates, monomers are more likely to be buried at larger monomer depth z_m, which results in a lower exposure probability. Similarly, a too high fragmentation velocity in SHAMPOO (i.e. the background model does not contain enough large grains) results in monomers spending much time in large aggregates which are not present as grains in the background model. This would thus also result in an underestimation of the exposure probability. For the volatile CHNOS budgets, this means that inconsistencies between the background model and the collision model will likely result in our model underestimating the effects of adsorption and photodesorption on the ice budgets of monomers. The results of Sect. 3.2 show that in particular photodesorption can have extreme effects on the evolution of monomers which are exposed frequently. Inconsistencies between the collision model and the background model would thus underestimate these effects, which highlights the importance of keeping the grain size distribution implied by the collision model and the grain size distribution of the background disk model as consistent as possible.

In our background disk model, we allow for the condensation and sublimation of molecules onto monomers. Since the position of the monomer between the time of condensation and sublimation of a particular molecular changes, the gas phase abundance of the molecule can become depleted in disk regions of systematic condensation onto monomers, and enriched in regions of systematic sublimation. SHAMPOO predicts the effects of a given chemical disk structure on the evolution of individual dust particles. However, it does not incorporate the effects of systematic behaviour of monomers on the disk chemical structure. It is expected that this systematic behaviour does affect the disk chemical structure on ~1 Myr timescales (e.g. Zhang et al. 2020; Krijt et al. 2020), which may in turn affect monomers which are processed in the same region at later times. In future work, we aim to quantify the initial net gas phase depletion-and enrichment fluxes at a certain disk region using the systematic behaviour of many monomers that visit that region. This approach would enable us to constrain the timescales where the condensation or sublimation of volatile molecules have large effects on local gas phase abundances. Altogether it becomes clear that predictions of monomer trajectories made with SHAMPOO are not reliable on timescales comparable to ~1 Myr due to the expected chemical feedback on the gas phase. However, other evolutionary processes such as disk evolution and planet formation are also expected to undermine the assumption of a static background disk environment on these timescales (e.g. Andrews 2020; Raymond & Morbidelli 2022). As a precaution, we have therefore chosen a timescale (100 kyr) significantly shorter than 1 Myr as the upper bound for the simulation time for individual monomers throughout this work.

4.2.2 Optical depth formalism

In Sect. 2.4 we introduced a new stochastic approach for estimating ice adsorption and photodesorption for a monomer embedded at a certain depth zm inside its home aggregate. In this approach, we assumed the monomer to be buried at a depth z_m below a homogeneous slab of monomers associated with an “optical depth” τ. However, the expression for τ derived in Eq. (51) strictly applies to aggregates of constant density (i.e. n_m does not vary as a function of time or position in the aggregate). In reality, however, n_m does vary as a function of position in the home aggregate and evolves over time due to collisional processing.

The most notable effect here is compaction (Weidling et al. 2009; Güttler et al. 2010; Zsom et al. 2010). This is a process which is not explicitly modeled in this work, although the assumed filling factor of ϕ = 0.1 is consistent with aggregates which have undergone some degree of compaction (e.g. Langkowski et al. 2008; Weidling et al. 2009). The fact that n_m depends on position in the aggregate and changes over time due to collisions has important effects on the exposure probability P_exp and the new monomer depth z_m after a collision. The results in Sect. 3.3 suggest that accumulation of new ice is more efficient in porous aggregates, and thus gradual compaction implies that ice formation is a process that becomes more and more confined towards the aggregate surface at later times.

Another simplification made is that the optical depth τ associated with a certain monomer depth z_m is assumed to be the same for impinging gas phase molecules and UV photons. In practice however, UV photons will likely be able to penetrate deeper into the aggregate since the minimum wavelength of UV photons (≳10⁻⁸ m) is either comparable or larger than the monomer size (s_m = 5 × 10⁻⁸). This would imply that the exposure probability for photodesorption at a given depth would be larger, and hence more monomers would be able to lose ice via photodesorption. For r ≂ 10 AU, the results from Sect. 3.2 suggest that the effects of photodesorption are significant for individual monomers in smaller, more fragile aggregates (υ_frag = 1 m s⁻¹) in that disk region, even with no distinction between the optical depth for photons and molecules. For these monomers, our current model may thus underestimate the average amount of ice lost due to photodesorption. Similarly, a lower filling factor also reduces the optical depth, and thus is more likely to leave monomers exposed to impinging UV photons. As a function of position, Figs. A.1–A.4 show that for larger radial distances (r ≳ 100 AU), UV photons can reach deeper in the disk. Therefore, our model may also underestimate the amount of ice lost via photodesorption in small, fragile aggregates at large radial distances. However, the examples in Sect. 3.2 do suggest that the effect of photodesorption is less important for monomers residing in aggregates with higher υ_frag, as the monomer is on average located at considerably larger z_m in these aggregates (Fig. 12). Therefore, we expect that the effect of a distinction between the optical depth for molecules and UV photons will remain limited for monomers in aggregates with higher υ_frag.

5 Conclusions

In this work we introduced and benchmarked the SHAMPOO code, which models the effects of dynamical, collisional and ice processes on the volatile CHNOS abundances of individual dust particles. In addition, we also motivated the need for this model, and used the SHAMPOO code to investigate the effects of different fragmentation velocities and aggregate filling factor on ice evolution at 10 AU. We summarize our conclusions as follows:

• The dynamical behaviour of a dust grain of given size can vary greatly as a function of position throughout the disk. This is a consequence of the fact that the grain size associated with the Stokes number for which aerodynamic drag starts to significantly affect the grain dynamics (St = 10⁻³) varies greatly throughout the disk (Fig. 7). This is primarily a result from the large variation of the gas density throughout the disk, and to a lesser extent also from the gas temperature.

• The disk regions where the timescales associated with dynamical, collisional, and ice processing couple varies, depending on grain size and ice species (Fig. 9). However, different grain sizes are connected through collisional processes, and the ice mantle of a dust particle can consist of multiple ice species. Therefore, we require a fully coupled modeling approach such as SHAMPOO to make physical estimates of the volatile CHNOS during the first stages of planetesimal formation.

• For a fragmentation velocity of υ_frag = 1 m s⁻¹, the monomer spends more time in smaller aggregates, and is therefore on average located closer to the aggregate surface (Figs. 10–13). A lower monomer depth increases the exposure probability P_exp, resulting in more periods of monomer exposure to the gas phase compared to monomers in aggregates of higher υ_frag. On the one hand this means that monomers in fragile aggregates can more frequently accumulate ice through adsorption. On the other hand, the exposure to the gas phase can also result in enhanced ice loss through photodesorption if the period of exposure coincides with a large vertical excursion of the monomer to the disk surface. Altogether, our results suggest that the interiors of fragile aggregates are more likely to have undergone extensive ice processing than aggregates with a higher fragmentation velocity (υ_frag = 5, 10 m s⁻¹).

• In the absence of dynamical transport, the average ice budgets of aggregates with ϕ = 1 are in good agreement with the ice abundance predicted by the time-dependent chemistry of the ProDiMo background disk model. Furthermore, very porous aggregates (ϕ = 10⁻³) have a higher exposure probability P_exp at a given depth. Therefore, gas phase molecules can reach deeper into these aggregates, and these aggregates accumulate more ice (Fig. 14). As a consequence, the ice mantles of monomers inside these aggregates are likely more processed compared to monomers inside more compact aggregates.

Altogether SHAMPOO provides a useful tool for probing the dust processing on a microscopic level. In future work, we aim to use a large set of its individual dust grain simulations to study the influence of non-local disk processing on local CHNOS abundances. Although the model has been developed to study the processing of CHNOS, we also noted applications of SHAMPOO outside its main purpose. Examples included the formation and transport of physically or chemically processed material, such as crystalized sillicates or hydrated minerals, and the formation and evolution of pebbles.

Acknowledgements

The authors would like to thank Rens Waters and Kaustubh Hakim for comments and insightful discussions, and an anonymous reviewer for constructive comments that helped to improve our simulations. This work is part of the second round of the Planetary and Exoplanetary Science Network (PEPSci-2), funded by the Netherlands Organization for Scientific Research (NWO).

Appendix A Structure and composition of the background disk models

The physical disk structure of the const, vFrag1, vFrag5 and vFrag10 background models is shown in Sect. A.1 while Sect. A.2 presents the chemical structures.

Appendix A.1 Physical background disk structures

Overview of the spatial structures of gas temperature T_g, dust temperature T_d, gas mass density ρ_g, dust mass density ρ_d and UV radiation field strength χ_RT provided by the background disk models.

Fig. A.1 Physical disk structure of the const-model.

Fig. A.2 Physical disk structure of the vFrag1-model.

Fig. A.3 Physical disk structure of the vFrag5-model.

Fig. A.4 Physical disk structure of the vFrag10-model.

Appendix A.2 Chemical background disk structure

Overview of the gas (left columns) and ice (right columns) number densities (n_x) of various molecules provided by the background disk models in m⁻³. In each figure from top to bottom: H₂O, NH₃, CO, CO₂, CH₄ and H₂S.

Fig. A.5 Chemical disk structure of the const-model.

Fig. A.6 Chemical disk structure of the vFrag1-model.

Fig. A.7 Chemical disk structure of the vFrag5-model.

Fig. A.8 Chemical disk structure of the vFrag10-model.

Appendix B Diffusion timescales inside the home aggregate

In our treatment of desorption in Sect. 2.4, molecules originating from monomers deep inside aggregates are assumed to escape the aggregate instantaneously. However, this assumption may be inappropriate for significant filling factors, as the mean free path for a molecule traveling inside an aggregate is usually shorter than s_a, which means that molecules liberated by desorption will undergo collisions with other monomers in an aggregate before escaping to the gas phase. This may result in re-adsorption, allowing to possibly retain more ice in the unexposed aggregate interiors than accounted for in the current model. Similarly, this mechanism could also provide a pathway for freshly accumulated ice molecules on exposed monomers to diffuse towards the unexposed monomers deeper in the aggregate. In this section, we develop a simple model to constrain the importance of this process by defining and estimating an upper bound of the timescale associated with the escape of molecules from the aggregate.

We assume that the liberated molecules effectively undergo a random walk as collisions with other monomers occur while escaping the aggregate. If liberated molecules only collide with monomers, the escape timescale τ_esc associated with a molecule depends on two main factors: the number of random walk collisions between the molecule and other monomers n_c before the molecule escapes the aggregate, and the fraction of these collisions which results in re-adsorption, which is given by the sticking factor S (Eq. (57)). Altogether we can thus express τ_esc as $${\tau }_{\text{esc}}=n_{\text{c}}\left({\tau }_{\text{mc}}+S{\tau }_{\text{d}}\right).$$(B.1)

Here, τ_mc denotes the molecule collision timescale, and can be expressed as $${\tau }_{\text{mc}}=\frac{{\lambda }_{\text{mfp}}}{\upsilon }=\frac{4s_m}{3\varphi \upsilon }.$$(B.2)

${\lambda }_{\text{mfp}}=1/\left(\pi s_m^2{n}_m\right)$ in this case denotes the mean free path inside the aggregate. We here use the monomer density n_m as given in Eq. (50). ν denotes the velocity of the molecule while it travels between monomers. Since we aim to obtain an upper bound on the escape timescale, we estimate ν with the molecule ice diffusion speed, which we estimate from the thermal hopping timescale τ_hop (i.e. the average time for a molecule to migrate from one site in the ice lattice to the next) (e.g. Thi 2015; Ruaud et al. 2016; Cuppen et al. 2017), $$\upsilon =\frac{d_{\text{s}}}{{\tau }_{\text{hop}}}=d_{\text{s}}{v}_x\exp \left(-\frac{E_{\text{diff},x}}{k_{\text{B}}{T}_d}\right).$$(B.3)

Here, d_s ≃ 10⁻¹⁰ m (e.g. Caselli & Ceccarelli 2012) denotes the average distance between two molecule binding sites in the ice lattice. v_x is the lattice vibration frequency as given by Eq. (59) and E_dilff,x denotes the diffusion energy, which we estimate as $E_{\text{diff,}x}=\frac{1}{2}{E}_{\text{ads,}x}$ Similarly, the desorption timescale for a single molecule τ_d (in addition to the ice mantle desorption timescale τ_des defined in Eq. (73)) for molecule species χ can be expressed as (e.g. Thi 2015) $${\tau }_{\text{d}}=\frac{1}{v_x}\exp \left(\frac{E_{\text{ads,}x}}{k_{\text{B}}{T}_{\text{d}}}\right).$$(B.4)

Altogether we write Eq. (B.1) as $${\tau }_{\text{esc}}=\frac{n_c}{v_x}\left[\frac{4s_m}{3\varphi d_s}\text{exp\hspace{0.17em}}\left(\frac{E_{\text{ads,}x}}{2k_{\text{B}}{T}_{\text{d}}}\right)+S\exp \left(\frac{E_{\text{ads,}x}}{k_{\text{B}}{T}_{\text{d}}}\right)\right].$$(B.5)

Fig. B.1 Example of a 3D random walk through an aggregate of radius sa = 10−5 m and ϕ = 0.5 (top view). The colorscale indicates the z-coordinate of the molecule. The black dot and cross denote the starting point and the location where the molecule crosses the aggregate surface, respectively.

Fig. B.2 Distribution of escape timescales of 100 model runs in aggregates of ϕ = 0.5 and varying radius. Molecules are always released from the center of the aggregate.

It becomes clear from Eq. (B.5) that τ_esc is very sensitive to T_d. On the one hand, this will result in extremely long τ_esc for aggregates significantly inside the ice-forming region, as thermal desorption becomes a very slow process (see also Sect. 3.1), resulting in extremely long τ_d. Furthermore, S ≃ 1 in these regions, which means that even if a molecule undergoes thermal desorption, it will be re-adsorbed almost instantaneously by another nearby monomer. On the other hand in warm regions where τ_d becomes very short and S → 0, we also expect that re-adsorption will not have a significant effect on monomer ice mantle composition, as τ_esc will be dominated by the molecule collision timescale τ_mc. Altogether Eq. (B.5) suggests that the effects of re-adsorption of species x on monomer ice mantle composition are the largest near the iceline associated with species x. However, another crucial factor in determining the value of τ_esc is the number of collisions before escape n_c, which depends on aggregate size s_a and filling factor ϕ.

In order to constrain n_c and provide an upper bound on τ_esc near the iceline, we consider a simple numerical 3D spherical random walk model. In this random walk model, individual molecules move from one position to the next in steps of size λ_mfp in a spherically random uniform direction. The simulation ends when the molecule distance to the aggregate center is larger than the aggregate radius (i.e. $x_{\text{mol}}^2+y_{\text{mol}}^2+z_{\text{mol}}^2>s_a^2$). This gives rise to trajectories of individual molecules through the aggregate. An example of such a trajectory is shown in Fig. B.1, where a molecule is released from the center of an aggregate with s_a = 10⁻⁵ m and ϕ = 0.5, and undergoes n_c = 7124 collisions before escaping the aggregate. Altogether individual molecules can take significant detours through the aggregate, resulting in a considerable variation in n_c and τ_esc associated with the trajectories of individual molecules. Therefore it is necessary to evaluate Eq. (B.5) for a larger set of molecules.

As a next step therefore we consider the distribution of τ_esc for 100 H₂O molecules escaping from the center of aggregates of s_a = 10⁻⁶, 10⁻⁵, and 10⁻⁴ m at the H₂O iceline (T_d ≈ 140 K, c.f. Fig A.1 and Fig. A.5). As we expect that the effects of readsorption are most important for larger aggregates with a short mean free path, we consider a filling factor ϕ = 0.5. The direct dependence of τ_esc on ϕ in Eq. (B.5) suggests that larger ϕ results in shorter τ_esc. However, a shorter mean free path will also result in larger n_c, which effectively results in the escape timescale being longer for higher ϕ. The resulting escape timescales distributions τ_esc for the different aggregate sizes are shown in Fig. B.2. Molecules escape faster from smaller aggregates due to the lower average number of collisions n_c before escape. Fig. B.2 reveals that the average value of τ_esc increases orders of magnitude per order of magnitude increase in s_m. This suggests that for aggregates larger than 10⁻⁴ m (i.e. millimetre- and centimetre-sized aggregates), the escape timescale may become comparable to the timescales of other disk processes (c.f. Fig. 8). Since the largest aggregates contain most of the dust mass (c.f. Fig. 4), Fig. B.2 implies that re-adsorption after thermal desorption in unexposed aggregate interiors may have significant effects on the average amount of ice retained in monomer ice mantles. Specifically, our model could overestimate the amount of ice lost via thermal desorption on unexposed monomers. However, an increase or decrease of 10 K from T_d was also found to result in a decrease or increase in τ_esc of almost 2 orders of magnitude, respectively. This suggests that the disk region where τ_esc couples to the other disk processes is spatially limited to a small region just behind the H₂O iceline. However, this also means that there always is a region where τ_esc couples to the other disk processes, and thus can affect ice composition on timescales similar as the other disk processes. We expect that this result for H₂O can be generalized to other chemical species, where the disk region where re-adsorption may affect the amount of ice of species x on monomers would be shifted to just behind the iceline of the respective species. Altogether the current model may underestimate amount of ice of species x retained in deeply embedded monomers just behind the iceline of species x.

References

All Tables

All Figures

	Fig. 1 Key model idea: a monomer of radius sm embedded at some monomer depth zm inside a home aggregate of effective radius sa.
In the text

	Fig. 2 Radial behaviour of amax for the four different background models considered throughout this work.
In the text

Fig. 3 Integrated residence time for monomers of sm = 5 × 10−8 m (blue), 3 × 10−4 m (orange), and 3 × 10−3 m (green) as a function of height z at r = 10 AU. The solid lines denote the histogram obtained from the trajectories of 103 monomers combined, whereas the dashed lines are the theoretical vertical density profiles expected from Eq. (8). The vertical dashed lines depict the locations where z = Hd for each monomer size.

In the text

Fig. 4 Integrated monomer residence time as a function of home aggregate size for different fragmentation velocities (solid coloured lines). The solid black lines on the left and right represent the minimum and maximum grain size occurring throughout the entire vFrag5 background disk model. The dashed lines all denote the same background dust mass distribution as a function of size, ρd(a) at r = 10 AU and z = 0, shifted vertically for comparison with the integrated monomer residence time.

In the text

	Fig. 5 Exposure probability as a function of monomer depth zm for different filling factors ϕ.
In the text

Fig. 6 Particle size as a function of time (left panel) and ice mass (scaled with the initial ice mass) as a function of radial distance (right panel) for pure H2O particles undergoing radial drift and thermal sublimation. Different colors denote different particles sizes. The trajectories predicted by the model of this work (solid lines) are compared with the trajectories calculated by Piso et al. (2015; dots). The vertical dashed line in the right panel denotes the position of the H2O iceline at 0.67 AU as calculated by Piso et al. (2015).

In the text

	Fig. 7 Grain size for which the local Stokes number is equal to 10−3, as a function of position in the const background disk model. The red curve encloses where particles of the size for which St = 10−3 are in the Stokes drag regime.
In the text

	Fig. 8 Behaviour of the timescales associated with various disk processes as a function of radial distance r. Timescales shown are for a monomer of size sm = 10−5 m with a pure H2O ice mantle, located at z/r = 0.05. The dashed line at t = 108 yr indicates the upper limit on timescales relevant for planet formation.
In the text

	Fig. 9 Process coupling behaviour for different grain sizes and ice mantle compositions. Colored areas indicate the decoupling of the different timescale groups. The filling of the coloured areas with solid lines or dots indicate the decoupling of a process category since it is much faster or slower, respectively. Black lines indicate icelines (τads = τdes).
In the text

Fig. 10 Evolution of the vertical positions (upper left panel), radial positions (upper center panel), temperatures (upper right panel), local UV fluxes (lower left), and local molecular number density in the gas phase (lower center) for monomers in home aggregates with different fragmentation velocity υfrag, released from r = 10 AU and z/r = 0. The black dashed line in the lower left panel indicates a UV radiation field strength of 1 FDraine. Linestyles in the lower center panel indicate number densities of different gas phase species.

In the text

	Fig. 11 Collision histories and resulting home aggregate sizes sa as a function of time for monomers in aggregates of different fragmentation velocities υfrag. Note that the vertical position of each collision event coincides with the home aggregate size sa after the collision has occurred.
In the text

	Fig. 12 Evolution of monomer depth zm for monomers embedded in home aggregates with υfrag = 1 m s−4 (left), 5 m s−1 (center), and 10 m s−1 (right). Green areas denote the times at which the monomer is exposed to gas phase molecules and UV photons, while the monomer is shielded from impinging molecules and UV photons otherwise. In each panel, the dashed line indicates zm = zcrit.
In the text

	Fig. 13 Evolution of the total ice mass mice of monomers in aggregates with different fragmentation velocity υfrag (left panel), together with snapshots of the monomer ice composition at 10 kyr (center panel), and 100 kyr (right panel).
In the text

Fig. 14 Evolution of the total amount of ice (left panel), H2O ice (center), and NH3 (right) averaged over 500 monomers for different home aggregate filling factors. Solid lines indicate the mean of the 500 monomers, whereas dotted lines indicate the standard deviation. The black dots indicate the ice evolution behaviour predicted by ProDiMo for the vFrag5-model. Trajectories of individual monomers are also shown in the background.

In the text

	Fig. A.1 Physical disk structure of the const-model.
In the text

	Fig. A.2 Physical disk structure of the vFrag1-model.
In the text

	Fig. A.3 Physical disk structure of the vFrag5-model.
In the text

	Fig. A.4 Physical disk structure of the vFrag10-model.
In the text

	Fig. A.5 Chemical disk structure of the const-model.
In the text

	Fig. A.6 Chemical disk structure of the vFrag1-model.
In the text

	Fig. A.7 Chemical disk structure of the vFrag5-model.
In the text

	Fig. A.8 Chemical disk structure of the vFrag10-model.
In the text

	Fig. B.1 Example of a 3D random walk through an aggregate of radius sa = 10−5 m and ϕ = 0.5 (top view). The colorscale indicates the z-coordinate of the molecule. The black dot and cross denote the starting point and the location where the molecule crosses the aggregate surface, respectively.
In the text

	Fig. B.2 Distribution of escape timescales of 100 model runs in aggregates of ϕ = 0.5 and varying radius. Molecules are always released from the center of the aggregate.
In the text