© ESO 2020

1 Introduction

More than 4000 extrasolar planets have been confirmed to date, but only a handful of them can be studied in detail (see e.g. Nikolov et al. 2016; Huitson et al. 2017; Birkby et al. 2017; Arcangeli et al. 2018). Indirect observations, for instance transmission spectroscopy, have demonstrated the presence of clouds (Sing et al. 2016; Nikolov et al. 2016; Pino et al. 2018; Gibson et al. 2017; Kirk et al. 2018; Tregloan-Reed et al. 2018). Far easier targets for atmosphere studies are brown dwarfs, which are very similar to planets with respect to their physical parameters and atmospheric processes. The coolest brown dwarfs (Y dwarfs) reach effective temperatures as low as 250 K (Leggett et al. 2017; Luhman 2014). The observation of brown dwarfs allows us to identify the vertical cloud structures (Apai et al. 2013; Buenzli et al. 2015; Yang et al. 2016; Helling & Casewell 2014). To date, between 1500 and 2000 brown dwarfs are known (depending on whether late-M dwarfs and/or early-L dwarfs are included; Gagné et al. 2015; Best et al. 2018) and are relatively well studied compared to the ~4000 extrasolar planets, for which the era of spectral analysis has only just begun.

Cloud formation has a profound impact on the remaining gas-phase abundances and radiative transfer effects, but cloud particles will also affect the ionisation state of the atmosphere, which is well known for Solar System objects (Helling et al. 2016a,b). Efforts are therefore ongoing to construct physical models describing the formation of clouds in exoplanet and brown dwarf atmospheres. Detailed models like this are necessary tools to provide the context for observations and to uncover processes that are not directly accessible by observations. Part of this effort is the consistent coupling of cloud formation with 1D atmosphere models with radiative transfer and convection (Tsuji et al. 1996; Ackerman & Marley 2001; Tsuji 2002; Witte et al. 2009; Allard et al. 2012; Juncher et al. 2017; also Helling et al. 2008a), but also in 3D in order to study the time-dependent climate of extrasolar planets (Lee et al. 2016; Lines et al. 2018) and to understand observational implications beyond 1D (Lee et al. 2017; Lines et al. 2018).

As our understanding of cloud formation progresses (e.g. Lee et al. 2015b; Krasnokutski et al. 2017; Hörst et al. 2019), including its implication for habitability (Narita et al. 2015), we start to refine our approaches. One long-standing discussion is how to model the element replenishment in 1D cloud-forming atmospheres, because without replenishment, a quasi-static atmosphere must be cloud free (Appendix A in Woitke & Helling 2004). Parmentier et al. (2013) used passive tracers to study the atmospheric mixing in 3D (shallow water approximation) simulations for irradiated dynamic but convectively stable atmospheres of (gas giant) planets. They observed that cloud particles are distributed throughout the whole atmosphere.

Parmentier et al. (2013) pointed out that in their relaxed hydrodynamical models, the upward transport by mixing statistically balances the downward transport due to particle settling. The mixing does not require convection in particular, but any large-scale vertical motions in the stably stratified atmosphere can provide such mixing. These flow patterns need to be revealed by global circulation models, for example driven by the day-night heating contrast. This assessment confirms our conclusion that the upward transport of condensable elements through the atmosphere by mixing is indeed the key for understanding cloud formation. However, challenges arise from the choice of the inner boundary condition (Carone et al. 2019), chemical gradients (Tremblin et al. 2019), and the need to include cloud particle feedback in order to test mixing parametrisations. A particular interesting case will be the ultra-hot Jupiters, where day and night sides can be expected to have very distinct (vertical) mixing patterns and scales. In this paper, we consider self-luminous gas giant planets, for which the irradiation from their host stars is negligible, such as young gas giant planets and brown dwarfs. Brown dwarf atmospheres are by now understood to be rather similar to gas giant planets, in particular atmospheres from low-gravity brown dwarfs and young gas giants (Charnay et al. 2018).

Moses et al. (2000) pointed out that large-scale mixing helps to homogenise a gravitationally stratified atmospheres consisting of different types of molecules. However, this only prevails up to a certain altitude above which gas-kinetic diffusion starts to dominate mixing (Zahnle et al. 2016). Different approaches have been chosen to represent this vertical mixing in 1D atmosphere models (Ackerman & Marley 2001; Woitke & Helling 2004; Helling et al. 2008a; Allard 2014; Juncher et al. 2017) and in 3D models (Lee et al. 2015a; Lines et al. 2018) by measuring vertical velocity fluctuations and deriving mixing parametrisations from 2D or 3D radiation-hydrodynamics simulations (Ludwig et al. 2002; Freytag et al. 2010; Parmentier et al. 2013; Zhang & Showman 2018).

Cloud formation modelling also becomes an increasingly important part of exoplanet and brown dwarf retrieval approaches for which, however, computational speed is an essential limitation. As part of the ARCiS¹ retrieval platform, Ormel & Min (2019) presented a new forward-model that consistently solves diffusive mixing and cloud particle growth for exoplanet atmospheres.

In this paper, we present a new theoretical approach that consistently combines cloud formation modelling with diffusive transport for element replenishment. After presenting the main formula body of our model in Sects. 2.1–2.3, we summarise our ansatz for handling the diffusion coefficient in Sect. 2.4, before we present our main results in Sects. 4 and 5. We conclude in Sect. 6. An overview of quantifying diffusion coefficients in the literature is provided in Appendix D.

2 Cloud formation with diffusive transport of gas and cloud particles

Cloud formation involves at least seed particle formation (nucleation), surface growth and evaporation, element depletion, gravitational settling, and element replenishment. During their descent through the atmosphere, cloud particles may change phase, or more general, chemical composition, and may collide with each other, which leads to further growth. These cloud formation processes have been described previously (Woitke & Helling 2003; Helling & Woitke 2006; Helling & Fomins 2013), and different cloud formation models have been compared by Helling et al. (2008a) with an update by Charnay et al. (2018). We therefore only provide a short summary here. A recent review can be found in Helling (2019).

Clouds are made of particles (aerosols, droplets, and solid particles). The formation of these particles requires condensation seeds, which are produced, for the case of the Earth atmosphere, by volcano eruptions, ocean sprays, and wild fires. In absence of these crucial processes, which all require the existence of a solid planet surface, cloud formation needs to start with the formation of seed particles through chemical reactions in the gas phase, involving the formation of molecular clusters. The formation of seed particles requires a highly supersaturated gas. Once such seed particles become available, many materials are already thermally stable and can simultaneously condense onto these surfaces. Nucleation and growth reduce the local element abundances and have a strong feedback on the local composition of the atmospheric gas. As macroscopic cloud particles form, they display a spectrum of sizes as well as a mixture of condensed materials. The local particle size distribution and the material mixture change as the cloud particles move through the atmosphere (hence, encounter different thermodynamic conditions), for example by gravitational settling (rain). Particle–particle collision will continue to shape the size distribution function. Cloud particles may break up into smaller units (shattering) or stick together to form even larger units (coagulation). Cloud particles may also be transported upward and downward by macroscopic mixing processes. Particle–particle processes are not part of our present model, which focuses on the formation of cloud particles and their feedback on the local chemistry through element depletion or enrichment. We note that the surface growth does shut off the nucleation process due to efficient element depletion (Lee et al. 2015b) such that a simultaneous treatment of nucleation and growth is required in order to calculate the number of cloud particles that form in the first place.

2.1 Cloud formation as a reaction-diffusion system

As introduced in Woitke & Helling (2003), we consider the evolution of the size distribution function f(V) [cm⁻⁶] of cloud particles in the particle volume interval V... V + dV as affected by advection, settling, surface reactions, and (new) by diffusion according to the following master equation: $$\frac{\partial (f(V)\text{d}V)}{\partial t}+\nabla (v(V)f(V)\text{d}V)={\sum }_k{R}_k\text{\hspace{0.17em}d}V - \nabla {\varphi }_{\text{d}}\text{\hspace{0.17em}d}V.$$(1)

R_k are the various gain and loss rates due to surface chemical reactions, which lead to growth and evaporation of the particles (see Eqs. (59)–(62) for large Knudsen numbers, and Eqs. (68)–(71) for small Knudsen numbers in Woitke & Helling 2003). A list of symbols with physical meaning and units is provided in Table 1. The volume of the particles V is chosen as size variable to formulate the material deposit by surface reactions in the most straightforward way. The last term in Eq. (1) accounts for the additional gains and losses due to diffusive mixing. The cloud particle velocity v(V) is assumed to be given by the hydrodynamical gas velocity v_gas plus a vertical equilibrium drift velocity ${\stackrel{°}{v}}_{\text{dr}}(V),$ $$v(V)=v_{\text{gas}}+{\stackrel{°}{v}}_{\text{dr}}(V).$$(2)

When we apply Fick’s first law (see e.g. Bringuier 2013), the diffusive flux ϕ_d of the cloud particles in volume interval V... V + dV (ϕ_d d V has units [cm⁻² s⁻¹]) is given by the concentration gradient of these particles, $${\varphi }_{\text{d}}\text{\hspace{0.17em}d}V=-\rho D_{\text{d}}\nabla \left(\frac{f(V)\text{d}V}{\rho }\right),$$(3)

where D_d [cm² s⁻¹] is the diffusion coefficient for these cloud particles and ρ [g∕cm³] is the gas density. We introduce moments of the cloud particle size distribution as $$\rho L_j={\displaystyle {\int }_{V_{\text{}\mathcal{l}}}^{\infty }\text{}}\text{}f(V)V^{j/3}\text{\hspace{0.17em}d}V.$$(4)

Multiplying Eq. (1) with V^j∕3 and integrating over volume, we obtain the following system of moment equations for large Knudsen numbers (see details in Woitke & Helling 2003): $$\begin{array}{lll}\frac{\partial (\rho L_j)}{\partial t}\hfill & +\hfill & \nabla (v_{\text{gas}}\rho L_j)=\underset{\text{nucleation}}{\underbrace{V_{\text{}\mathcal{l}}^{j/3}{J}_{\star }}} + \underset{\text{growth}}{\underbrace{\frac{j}{3}\chi \rho L_{j-1}}}\hfill \\ \hfill & -\hfill & \underset{\text{drift}}{\underbrace{\nabla {\displaystyle {\int }_{V_l}^{\infty }{v}_{\text{dr}}}(V)f(V)V^{j/3}\text{\hspace{0.17em}d}V}} - \underset{\text{diffusion}}{\underbrace{\nabla {\displaystyle {\int }_{V_l}^{\infty }{\varphi }_{\text{d}}}{V}^{j/3}\text{\hspace{0.17em}d}V}},\hfill \end{array}$$(5)

where J_⋆ [s⁻¹ cm⁻³] is the nucleation rate and χ [cm s⁻¹] the net growth velocity. For large Knudsen numbers and subsonic velocities (Epstein regime), the equilibrium drift velocity, also called final fall speed, is given by Schaaf (1963), $${\stackrel{°}{v}}_{\text{dr}}=-\frac{\sqrt{\pi }g{\rho }_{\text{d}}a}{2\rho c_{\text{T}}}\widehat{r},$$(6)

where a is the particle radius, ρ_d is the cloud particle material density, $\widehat{r}$ is the unit vector pointing away from the centre of gravity, and g is the gravitational acceleration. $c_{\text{T}}=\sqrt{2kT/\overline{\mu }}$ is an abbreviation, T is the temperature, k is the Boltzmann constant, and $\overline{\mu }$ is the mean molecular weight of the gas particles.

Using Eq. (6) with $a={(\frac{3V}{4\pi })}^{1/3}$ and assuming that the particle diffusion coefficient D_d is independent of size, we can carry out the integrations in Eq. (5). The final result is $$\begin{array}{lll}\frac{\partial (\rho L_j)}{\partial t}+\nabla (v_{\text{gas}}\rho L_j)\hfill & =\hfill & V_{\text{}\mathcal{l}}^{j/3}{J}_{\star }+\frac{j}{3}\chi \rho L_{j-1}\hfill \\ \hfill & \hfill & +\nabla \left(\xi \frac{{\rho }_{\text{d}}}{c_T}{L}_{j+1}\widehat{r}\right)+\nabla \left(D_{\text{d}}\rho \nabla L_j\right),\hfill \end{array}$$(7)

with abbreviation $$\xi =\frac{\sqrt{\pi }}{2}{\left(\frac{3}{4\pi }\right)}^{1/3}g.$$(8)

A size-dependent diffusion coefficient, D_d(V), would lead to an open set of moment equations. This has been discussed by Helling & Fomins (2013).

2.2 Generalisation to mixed materials

We assume that all cloud particles are perfect spheres with a well-mixed material composition that is independent of size, but depends on time and location in the atmosphere (Helling & Woitke 2006; Helling et al. 2008c). Using the index s = 1 ... S to distinguish between the different solid materials, for example Al₂ O₃[s], TiO₂ [s], Mg₂ SiO₄[s] and Fe[s], we write $$V={\displaystyle \sum _s{V}^s}\text{\hspace{0.17em},}{L}_3={\displaystyle \sum _s{L}_3^s}\text{\hspace{0.17em},}{J}_{\star }={\displaystyle \sum _s{J}_{\star }^s}\text{\hspace{0.17em},}\chi ={\displaystyle \sum _s}{\chi }^s\text{\hspace{0.17em},}{b}^s=\frac{L_3^s}{L_3},$$(9)

where b^s is the volume fraction of material s in the cloud particles². The mean cloud particle material density is given by ρ_d = ∑_sb^sρ_s, where ρ_s is the mass density of a pure material s. Most materials will not nucleate themselves ($J_{\star }^s=0$), but will use alien nuclei to grow on. Using this approximation, we can split the third-moment equation into a set of third-moment equations for single materials as follows: $$\begin{array}{lll}\frac{\partial (\rho L_3^s)}{\partial t}+\nabla (v_{\text{gas}}\rho L_3^s)\hfill & =\hfill & V_{\text{}\mathcal{l}}{J}_{\star }^s+{\chi }^s\rho L_2\hfill \\ \hfill & \hfill & +\nabla \left(\xi \frac{{\rho }_{\text{d}}}{c_T}{b}^s{L}_4\widehat{r}\right)+\nabla \left(D_{\text{d}}\rho \nabla L_3^s\right).\hfill \end{array}$$(10)

When we sum Eq. (10) for all solids s, we again obtain Eq. (7) for j = 3. The different materials grow at different speeds that depend on the amount of available atoms and molecules in the gas phase and on the supersaturation ratio. Islands of some materials may grow, whereas others are thermally unstable and shrink. This behaviour is obtained by summing the contributions of all surface reactions r = 1... R (for examples, see Table 1 in Helling et al. 2008c) $$\chi ={\displaystyle \sum _s}{\chi }^s={\displaystyle \sum _s{\displaystyle \sum _r{c}_r^s}}{V}_0^s,$$(11)

where $V_0^s[{\text{cm}}^{\text{3}}]$ is the volume of one unit of material s in the solid state and $c_r^s[{\text{cm}}^{\text{-2}}{\text{s}}^{\text{-1}}]$ is the effective surface reaction rate, $$c_r^s=\sqrt[3]{36\pi }\frac{{\nu }_r^s{n}_r^{\text{key}}{\text{v}}_r^{\text{rel}}{\alpha }_r}{{\nu }_r^{\text{key}}}\left(1-\frac{1}{S_{\text{}r}}\right)\times \{\begin{array}{cc}1& \text{if}{S}_{\text{}r}\ge 1\\ b^s& \text{if}{S}_{\text{}r}<1\end{array}\text{}\text{},$$(12)

where $n_r^{\text{key}}$ is the particle density [cm⁻³] of the key species of surface reaction r, α_r is the sticking probability, ${\nu }_r^{\text{key}}$ is its stoichiometric factor in that reaction, ${\text{v}}_r^{\text{rel}}=\sqrt{kT/(2\pi m^{\text{key}})}$ is its thermal relative velocity, and m^key is its mass. These growth rates are derived from a simple hit-and-stick model, where we commonly assume α_r = 1. The impact of the limited number of known α_r ≠ 1 has been studied by Helling & Woitke (2006). S_r is the reaction supersaturation ratio as introduced in Appendix B of Helling & Woitke (2006). For example, in the reaction $$\text{2}\text{\hspace{0.17em}CaH} \text{+} \text{2}\text{\hspace{0.17em}Ti} \text{+} \text{6}{\text{H}}_{\text{2}}\text{O} ⟷ \text{2}{\text{CaTiO}}_{\text{3}}[\text{s}] \text{+} \text{7}{\text{H}}_{\text{2}}$$(13)

the key species is either CaH or Ti, depending on which species is less abundant. ${\nu }_r^{\text{key}}=2$ in either case, s = CaTiO₃[s] is the solid species, and ${\nu }_r^s=2$ units of CaTiO₃[s] are produced by one reaction.

We note that Eq. (12) differs slightly from our previous definition (Eq. 4 in Helling et al. 2008c). The new growth/evaporation rates now always change sign at S_r = 1 as they should, independent of the value of b^s. When supersaturated (S_r > 1), we assume that the total surface of the particles acts as a funnel to collect the impinging molecules from the gas phase, followed by fast hopping to find a place on a matching island of kind s (see Fig. 1 in Helling & Woitke 2006). For under-saturation (S_r < 1), however, we assume that the molecules triggering the evaporation processes must hit one of the islands of the matching kind, the probability of which is b^s.

2.3 Element conservation with diffusive replenishment

An integral part of our cloud formation model is the element conservation. We must identify a replenishment mechanism that is able to compensate for the losses of elements due to cloud particle formation and settling in the upper atmosphere. In this paper, we include diffusion of gas particles along their concentration gradients. As cloud particles form, they consume certain elements in the upper atmosphere, creating a negative element abundance gradient. Thus, gas particles containing these elements will ascent diffusively in the atmosphere. Analogously to the formulation of the master equation for the dust particles, we formulate the element conservation as diffusion-reaction system, $$\begin{array}{lll}\frac{\partial (n_{\langle \text{H}\rangle }{\varepsilon }_k)}{\partial t}+\nabla (v_{\text{gas}}{n}_{\langle \text{H}\rangle }{\varepsilon }_k)\hfill & =\hfill & -{\displaystyle \sum _s{\nu }_k^s}{N}_{\mathcal{l}}{J}_{\star }^s - \rho L_2{\displaystyle \sum _s{\displaystyle \sum _r{\nu }_k^s}}{c}_r^s\hfill \\ \hfill & \hfill & + \nabla \left(D_{\text{gas}}{n}_{\langle \text{H}\rangle }\nabla {\varepsilon }_k\right),\hfill \end{array}$$(14)

where ɛ_k is the abundance of element k with respect to hydrogen. The chemical reactions leading to nucleation and growth appear as negative source terms here because they consume elements. We chose n_⟨H⟩ as density variable in Eq. (14), the total hydrogen nuclei particle density, which is proportional to ρ in hydrogen-dominated atmospheres. n_⟨H⟩ɛ_k is the total number density [cm⁻³] of nuclei of element k in any chemical form. ${\nu }_k^s$ is the stoichiometric factor of element k in solid s, for example ${\nu }_{\text{Ti}}^{\text{TiO2}[\text{s}]}=1$, and D_gas [cm² s⁻¹ ] is the gas diffusion coefficient. For simplicity, we assumed that all molecules are transported by the same diffusion coefficient, which is valid within a factor 2 or 3 for gas-kinetic diffusion (sometimes called the binary diffusion coefficient, see Eq. (16) and Appendix D), and is entirely justified when eddy diffusion dominates. The involved diffusive gas element flux ${\varphi }_k^{\text{diff}}$ [cm⁻²s⁻¹] is given by $${\varphi }_k^{\text{diff}}=-D_{\text{gas}}{n}_{\langle \text{H}\rangle }\nabla {\varepsilon }_k.$$(15)

2.4 Gas diffusion coefficient

The diffusion coefficients provide the kinetic information for calculating the transport rates from concentration gradients (e.g. Lamb & Verlinde 2011). In general, gas and cloud particles diffuse with different efficiencies because of their different inertia and collisional cross-sections with the surrounding gas.

The determination of the gas diffusion coefficient D_gas is of crucial importance for our model. We included gas-kinetic diffusion and large-scale turbulent (eddy) diffusion as mixing processes. The gas-kinetic diffusion coefficient is given by $$D_{\text{micro}}=\frac{1}{3}{\text{v}}_{\text{th}}\mathcal{l},$$(16)

where ℓ = 1∕(σ n) is the mean free path, n is the total gas particle density, and σ ≈ 2.1 × 10⁻¹⁵ cm² is a typical cross-section for collisions between the molecules under consideration with H₂. The thermal velocity is defined as ${\text{v}}_{\text{th}}=\sqrt{8kT/(\pi m_{\text{red}})}$, where m_red is the reduced mass for collisions between the molecule and H₂ (Woitke & Helling 2003). This gas-kinetic diffusion ∝ 1∕n is negligible in the lower high-density layers of brown dwarf and planetary atmospheres, where instead mixing by large-scale turbulent orconvective motions is the dominant mixing process (Ackerman & Marley 2001; Ludwig et al. 2002; Woitke & Helling 2004; Allard et al. 2012; Parmentier et al. 2013; Lee et al. 2015a). The large-scale (turbulent, convective, or eddy) gas diffusion coefficient is given by $$D_{\text{mix}}\approx \langle {\text{v}}_z\rangle L\text{with}L=\alpha H_{\text{p}},$$(17)

where ⟨v_z⟩ is the root-mean-square average of the fluctuating part of the vertical velocity in the atmosphere,considering averages over sufficiently large volumes and/or long integration times, L is the mixing length, and H_p is the local pressure scale height. α is a dimensionless parameter of the order of one. We used α = 1 in this work,but α can be fine-tuned to describe the actual mixing scales as revealed by detailed hydrodynamic modelling. Inside the convective part of the atmosphere ⟨v_z⟩≈v_conv was assumed, where v_conv is the convective velocity, which is an integral part of stellar atmosphere models, derived from mixing length theory in 1D models. Above the convective atmosphere, where the Schwarzschild criterion for convection is false, ⟨v_z ⟩ decreases rapidly with increasing z, but never quite reaches zero due to convective overshoot (see e.g. Brandenburg 2016). We applied a power law in log p to approximate this behaviour, $$\log \langle {\text{v}}_z\rangle =\log {\text{v}}_{\text{conv}}-{\beta }^{\prime }\cdot \max \{0,\log p_{\text{conv}}-\log p(z)\},$$(18)

with a free parameter β′≈ 0.0... 2.2 (Ludwig et al. 2002). The total gas diffusion coefficient is then $$D_{\text{gas}}=D_{\text{mix}}+D_{\text{micro}}.$$(19)

At high altitudes, the gas density n is low and hence D_micro is large, whereas D_mix is small when β′ > 0. Therefore, at some pressure level in the atmosphere, the gas-kinetic diffusion will start to dominate. Figure 1 shows a typical structure as assumed in our models. The minimum of D_gas around 10⁻³ mbar correspondsto the cross-over point (called the homopause), upward of which D_micro dominates and the atmospheres is not well mixed. Moses et al. (2000) drew similar conclusions concerning Saturn’s atmosphere. The maximum of D_gas around p_conv = 0.2 bar results from the start of the convective layer, within which both ⟨v_z⟩ and D_gas are approximately constant. Appendix D summarises some of the formulas currently applied in the literature for the gas diffusion coefficient.

Fig. 1 Gas diffusion coefficient Dgas (Eq. (19)) in the new DIFFUDRIFT model for a brown dwarf atmosphere model with Teff = 1800 K, log g =3 [cm s−2] and solar abundances. The grey line shows the inverse mixing timescale τmix as assumed in the previous DRIFT model. τmix is calculated according to Eq. (9) in Woitke & Helling (2004). Both quantities are computed for β = β′ = 1, and both y-axes show exactly eight orders of magnitude.

2.5 Cloud particle diffusion coefficient

The diffusion of solid particles due to turbulent gas fluctuations was studied, in consideration of protoplanetary discs, by Dubrulle et al. (1995), Schräpler & Henning (2004), Youdin & Lithwick (2007), and Riols & Lesur (2018). All works applied mean field theory (also called Reynolds decomposition ansatz), where the densities and velocities of both the particles and the gas are decomposed into a mean component (that depends only on z) and a small fluctuating part.

The response of the solid particles to the turbulent gas variations is then determined by comparing two timescales. The stopping or frictional coupling timescale is given by $${\tau }_{\text{stop}}=\frac{a{\rho }_{\text{d}}}{v_{\text{th}}\rho }\text{with}{v}_{\text{th}}=\sqrt{\frac{8kT}{\pi \overline{\mu }}},$$(20)

where a is the particle radius. Equation (20) follows from a general relaxation ansatz ${\tau }_{\text{stop}}=m{(\partial F_{\text{fric}}/\partial v_{\text{dr}}{|}_{v_{\text{dr}}^{{}^{°}}})}^{-1}$, see Eq. (21) in Woitke & Helling (2003), for the special case of large Knudsen numbers in a subsonic flow (the so-called Epstein regime), which we assume is valid here.

The second timescale is the eddy turnover or turbulence correlation timescale τ_eddy(l) in consideration of a spectrum of different turbulent modes associated with different wave-numbers k or different spatial eddy sizes l. In a Kolmogorov type of power spectrum P(k) ∝ k^−5∕3, any given cloud particle of size a tends to co-move with all sufficiently large and slow turbulent eddies, whereas its inertia prevents it from following the short-term, small turbulence modes.

In order to arrive at an effective particle diffusion coefficient, the advective effect of all individual turbulent eddies has to be averaged and thereby transformed into a collective diffusive effect. This procedure is carried out with different procedures and approximations. The result of Schräpler & Henning (2004, see their Eq. (27)) reads $$D_{\text{d}}=\frac{D_{\text{mix}}}{1+\text{St}}\text{with}\text{St}=\frac{{\tau }_{\text{stop}}}{{\tau }_{\text{eddy}}},$$(21)

where D_d is the size-dependent cloud particle diffusion coefficient and St is the Stokes number of the particle in consideration of the largest eddy size L. The eddy turnover timescale of the largest turbulence mode is given by $${\tau }_{\text{eddy}}=\frac{L}{\langle {\text{v}}_z\rangle }.$$(22)

Both the size of the largest eddy L and the average of the fluctuating part of the vertical velocity ⟨v_z⟩ are assumedto be identical to the mixing length and velocity appearing in Eq. (17). Combining the above equations, we find $$\text{St}=\frac{a{\rho }_{\text{d}}{D}_{\text{mix}}}{v_{\text{th}}\rho L^2}.$$(23)

The impact of the size dependence of D_d on the cloudparticle moments was explored by Helling & Fomins (2013), who showed that this leads to an open set of moment equations, which seems impractical for an actual solution. In the frame of this work, we therefore only explore the two limiting cases of very large and very small Stokes numbers throughout the atmosphere. For small particles with St ≪ 1, we have D_d → D_mix and for huge particles St →∞, we have D_d → 0. $$\begin{array}{lll}\text{case 1:}\hfill & D_{\text{d}}=D_{\text{mix}}\hfill & \text{if all cloud particles are small,}\hfill \\ \text{case 2:}\hfill & D_{\text{d}}=0\hfill & \text{if all cloud particles are large}\text{.}\hfill \end{array}$$(24)

Our results show that both approximations lead to rather similar cloud structures in the models explored so far, that is, the inclusion of turbulent cloud particle motions does not seem to be a critical ingredient to our present model. However, in preliminary models for hot Jupiters, where D_mix(z) is flatter or even increases with height, this might be different.

3 Static plane-parallel atmosphere

Before we proceed with the numerical solution of the full time-dependent model of cloud formation in diffusive media, we first discuss the 1D static case in order to better understand the expected results from these equations. When we consider the plane-parallel (∇→d∕dz), static (v_gas = 0), and stationary case (∂∕∂t = 0), our Eqs. (7), (10) and (14) simplify to $$\begin{array}{lll}0\hfill & =\hfill & V_{\text{}\mathcal{l}}^{j/3}{J}_{\star }+\frac{j}{3}\chi \rho L_{j-1}+\xi \frac{\text{d}}{\text{d}z}\left(\frac{{\rho }_{\text{d}}}{c_T}{L}_{j+1}\right)+\frac{\text{d}}{\text{d}z}\left(D_{\text{d}}\rho \frac{\text{d}{L}_j}{\text{d}z}\right)\hfill \\ 0\hfill & =\hfill & V_{\text{}\mathcal{l}}{J}_{\star }^s+\rho L_2{\displaystyle \sum _r{c}_r^s}{V}_0^s+\xi \frac{\text{d}}{\text{d}z}\left(\frac{{\rho }_{\text{d}}}{c_T}{b}^s{L}_4\right)+\frac{\text{d}}{\text{d}z}\left(D_{\text{d}}\rho \frac{\text{d}{L}_3^s}{\text{d}z}\right)\hfill \\ 0\hfill & =\hfill & -{\displaystyle \sum _s{\nu }_k^s}{N}_{\mathcal{l}}{J}_{\star }^s-\rho L_2{\displaystyle \sum _s{\displaystyle \sum _r{\nu }_k^s}}{c}_r^s+\frac{\text{d}}{\text{d}z}\left(D_{\text{gas}}{n}_{\langle \text{H}\rangle }\frac{\text{d}{\varepsilon }_k}{\text{d}z}\right).\hfill \end{array}$$

3.1 Total element fluxes

In the hydrostatic stationary case, the total vertical flux of elements (due to vertical settling of cloud particles and diffusive transport) must be zero everywhere in the atmosphere and for each element. This conclusion can be derived formally by summing Eq. (27) and ∑_s (Eq. (26)) $\cdot {\nu }_k^s/V_0^s$, using $V_{\mathcal{l}}=N_{\mathcal{l}}{V}_0^s$. The chemical source terms (nucleation and growth terms) cancel out exactly, and in case D_d = 0, we find $$\begin{array}{lll}\hfill & \hfill & \frac{\text{d}}{\text{d}z}\left(D_{\text{gas}}{n}_{\langle \text{H}\rangle }\frac{\text{d}{\varepsilon }_k}{\text{d}z}\right)+\xi \frac{\text{d}}{\text{d}z}\left(\frac{{\rho }_{\text{d}}}{c_T}{L}_4{\displaystyle \sum _s\frac{{\nu }_k^s{b}^s}{V_0^s}}\right)=0\hfill \\ \hfill & \hfill & \Rightarrow \underset{-{\varphi }_k^{\text{diff}}}{\underbrace{D_{\text{gas}}{n}_{\langle \text{H}\rangle }\frac{\text{d}{\varepsilon }_k}{\text{d}z}}} + \underset{{\varphi }_k^{\text{settle}}}{\underbrace{\xi \frac{{\rho }_{\text{d}}}{c_T}{L}_4{\displaystyle \sum _s\frac{{\nu }_k^s{b}^s}{V_0^s}}}}={\text{const}}_k.\hfill \end{array}$$(28)

${\varphi }_k^{\text{diff}}$ [cm⁻² s¹ ] is the upward element flux by diffusion in the gas phase and ${\varphi }_k^{\text{settle}}$ is the downward flux of elements contained in the settling cloud particles at this point. Equation (28) would still allow for solutions with constant (i.e. time-independent and height-independent) total element fluxes throughout the atmosphere, but this would require matching feeding and removing rates at the bottom and the base of the atmosphere, which does not seem to be physically plausible.

3.2 Equlibrium between gas diffusion and particle settling

From Eq. (28), assuming const_k = 0, we find $${\varphi }_k^{\text{diff}}={\varphi }_k^{\text{settle}}\Rightarrow \frac{\text{d}{\varepsilon }_k}{\text{d}z}=-\frac{\xi {\rho }_{\text{d}}{L}_4}{c_T{D}_{\text{gas}}{n}_{\langle \text{H}\rangle }}{\displaystyle \sum _s\frac{{\nu }_k^s{b}^s}{V_0^s}} \le 0.$$(29)

The element abundance gradients in a cloudy, static (v_gas = 0), and stationary (∂∕∂t = 0) atmosphere must be negative. The downward transport of elements through the precipitation of cloud particles must be balanced by an upward-directed diffusive flux of elements in the gas phase, which requires a negative concentration gradient in the gas phase. This conclusion is correct whenever cloud particles are present (L₄ > 0) and gravity is active (ξ > 0), otherwise, the gas element abundances are constant. The abundance gradients of different elements are proportional to the element composition of the settling cloud particles at this point. The abundances of all elements k involved in cloud formation must decrease monotonically toward the top of the atmosphere.

4 Numerical solution of the time-dependent cloud formation problem

Equations (25)-(27) form a system of (3 + S + K) coupled second-order differential equations, which can be transformed into a system of 2 × (3 + S + K) first-order ordinary differential equations (ODEs). Unfortunately, we have not been able to solve this ODE system directly. The boundary conditions are partly given at the lower and partly at the upper boundary of the model, see Sect. 4.4. The integration direction must be downward in order to model the nucleation of new cloud particles. Hence, we tried a shooting method where ɛ_k (z_max) was varied at the top of the atmosphere until ɛ_k(z_min) was met, that is, the given values in the deep atmosphere. We found it impossible to proceed this way. A tiny change of ɛ_k (z_max) in the 12th digit was still observed to change ɛ_k(z_min) by a factor of two. The reason for this extreme sensitivity seems to be the nucleation rate with its threshold character as a function of supersaturation, and hence as a function of ɛ_k.

Looking for alternatives, we found that a simulation of the time-dependent equations on a given vertical grid can be performed by means of the operator splitting method as explained in Sect. 4.2. We evolved the atmospheric cloud structure $\{L_j(z,t),L_3^s(z,t),{\varepsilon }_k(z,t)\}$ for a sufficiently long time, until it approached the time-independent case $\{L_j^{°}(z),L_3^{s,°}(z),{\varepsilon }_k^{°}(z)\}$, which is the stationary structure we are interested in. Assuming a plane-parallel (∇ →d∕dz) and static (v_gas = 0) atmosphere,Eqs. (7), (10), and (14), including the time-dependent terms, read $$\begin{array}{lll}\frac{\text{d}}{\text{d}t}(\rho L_j)\hfill & =\hfill & V_{\text{}\mathcal{l}}^{j/3}{J}_{\star }+\frac{j}{3}\chi \rho L_{j-1}\hfill \\ \hfill & \hfill & +\xi \frac{\text{d}}{\text{d}z}\left(\frac{{\rho }_{\text{d}}}{c_T}{L}_{j+1}\right)+\frac{\text{d}}{\text{d}z}\left(D_{\text{d}}\rho \frac{\text{d}{L}_j}{\text{d}z}\right)(j=1,2,3)\hfill \\ \frac{\text{d}}{\text{d}t}(\rho L_3^s)\hfill & =\hfill & V_{\text{}\mathcal{l}}{J}_{\star }^s+{\chi }^s\rho L_2\hfill \\ \hfill & \hfill & +\xi \frac{\text{d}}{\text{d}z}\left(\frac{{\rho }_{\text{d}}}{c_T}{b}^s{L}_4\right)+\frac{\text{d}}{\text{d}z}\left(D_{\text{d}}\rho \frac{\text{d}{L}_3^s}{\text{d}z}\right)(s=1\mathrm{...}S)\hfill \\ \frac{\text{d}}{\text{d}t}(n_{\langle \text{H}\rangle }{\varepsilon }_k)\hfill & =\hfill & -{\displaystyle \sum _s{\nu }_k^s}{N}_{\mathcal{l}}{J}_{\star }^s-\rho L_2{\displaystyle \sum _s{\displaystyle \sum _r{\nu }_k^s}}{c}_r^s\hfill \\ \hfill & \hfill & +\frac{\text{d}}{\text{d}z}\left(D_{\text{gas}}{n}_{\langle \text{H}\rangle }\frac{\text{d}{\varepsilon }_k}{\text{d}z}\right)(k=1\mathrm{...}K).\hfill \end{array}$$

4.1 Closure condition

The moment Eqs. (30) and (31) are not closed because L₄ appears twice on the right side, a consequence of larger particles settling faster (Eq. (6)). Therefore, a numerical solution requires a closure condition as $$L_4=F(L_0,L_1,L_2,L_3).$$(33)

We used theclosure condition explained in Appendix A.1 of Helling et al. (2008c). The idea is to approximate the particle size distribution f by a double δ-function that has four parameters. These parameters are determined by matching the given moments L₀, L₁, L₂, and L₃, and result in the fourth-moment L₄ according tothe definition of the dust moments (Eq. (4)).

4.2 Operator splitting method

Figure 2 visualises our numerical approach using the operator splitting method (Klein 1995). We list the steps below.

• 1

  We update L_j and $L_3^s$ only according to the settling source terms (the terms on the right side of Eqs. (30) and (31) containing L_j+1 and L₄), applying half a time step Δt∕2, see details in Appendix B.

• 2

  We call the diffusion solver for half a time step Δt∕2 to update ɛ_k, and optionally, L_j and $L_3^s$, if the cloud particles are to be diffused as well, see Appendix A.

• 3

  We integrate the chemical source terms (nucleation, growth, and evaporation) for a full timestep Δt. These equations are stiff at high densities and require an implicit integration scheme. We use the implicit ODE-solver LIMEX 4.2A1 (Deuflhard & Nowak 1987). The computation of the chemical source terms on the right-hand side proceeds as follows: (i) for a given temperature T, density n_{⟨H ⟩}, and element abundances ɛ_k, we call the equilibrium chemistry code GGCHEM (Woitke et al. 2018) to calculate all molecular concentrations; (ii) these results are used to calculate the reaction supersaturation ratios S_r ; (iii) the nucleation rates $J_{\star }^s$ and net surface reaction rates $c_r^s$ are determined.

• 4

  We complete the timestep by calling again the diffusion solver for Δt∕2 and the settling solver for Δt∕2 in this order.

• 5

  Checkpoint and output files are written for visualisation. The method is computationally limited by the time consumption for the implicit integration of the chemical source terms, which requires numerous calls of the equilibrium chemistry. Therefore we do not split CF (Fig. 2), but set it with a full time step in the centre of the operator splitting calling sequence. The cloud formation part of the code is parallelised and can be executed for all atmospheric layers independent of each other.

Fig. 2 Operator splitting calling sequence. Sett = gravitational settling, Diff = diffusion, CF = cloud formation (nucleation, growth and evaporation), and OP = output. 1/2 means half a time step, and 1 is a full time step.

4.3 Timestep control

In order to produce accurate second-order solutions, the timestep must be limited to ensure that each operator remains in the linear regime. For example, the sole application of a cloud-chemistry timestep must not change the amount of dust or the element abundances substantially in any computational cell. In order to achieve code stability and accuracy, we limit the timestep as follows:

• 1

  The cloud particles must not jump over layers by settling $$\Delta t<0.5\frac{\Delta z}{v_{\text{dr},j}},$$(34)

  where Δz is the vertical grid resolution and v_dr,j is the mean drift velocity affecting moment ρL_j, as given by Eq. (B.3). This is the usual Courant-Friedrichs-Lewy (CFL) criterion to stabilise explicit advection scheme, with an additional safety factor 1/2.

• 2

  Nucleation and cloud particle growth and evaporation, as integrated over Δt, must not change any of the gas element abundances by more than a given maximum relative change (the default accuracy is 15%).

• 3

  The time step must not exceed the maximum explicit diffusion timestep (Eq. (A.27)).

If one of these criteria becomes false during the simulation, the time step is discarded and Δt reduced. If the criteria are met easily, however, Δt is increased for the subsequent time step.

4.4 Boundary conditions

As our upper boundary condition, we assumed that no cloud particles settle into the model volume from above v_dr,j(z_max) = 0. In the diffusion solver, we used a zero-flux (closed box) upper boundary condition, that is, the gradients of ɛ_k were assumed to be zero at z = z_max. The same applies to the cloud particle moments $\frac{\text{d}}{\text{d}z}{L}_j(z_{\text{max}})=0$ if they are to be diffused as well.

The lower boundary was placed well below the main silicate cloud layer to ensure that the temperature was too high to allow any cloud particles to exist near the lower boundary L_j (z_min) = 0. We also demanded that the element abundances at the lower boundary equal the given values as present deep in the atmosphere ${\varepsilon }_k(z_{\text{min}})={\varepsilon }_k^0$, where the ${\varepsilon }_k^0$ are considered as free parameters, for example solar abundances (Asplund et al. 2009).

4.5 Initial conditions

We started all simulations from a cloud-free atmosphere L_j (z, t = 0) = 0. We experimented with two cases for the element abundances: (1) an “empty” atmosphere ɛ_k (z, t = 0) = 0 or (2) a “full” atmosphere ${\varepsilon }_k(z,t=0)={\varepsilon }_k^0$, where the index k is applied to all elements that might be transformed completely into solids (k = Si, Mg, Fe, Al, Ti, etc.), but not H, He, C, N, O, etc. For the latter elements we set ${\varepsilon }_k(z,t=0)={\varepsilon }_k^0$ in both cases. We found an identical final structure in both cases (see Appendix C), which is very reassuring. The models calculated from initial condition (2), however, need much more computational time to complete. In this case, the nucleation rate is initially huge and a very large number of tiny cloud particles are created shortly after initialisation, which take a long time tosettle down in the atmosphere.

In order to reach the final relaxed, time-independent state, the model must be evolved until (i) the atmosphere is completely replenished several times by fresh elements ascending diffusively from the lower boundary to the very top and (ii) new grains formed high in the atmosphere have sufficient time to settle down to the cloud base and evaporate. In comparison, the chemical processes are typically quite fast. We need to evolve one model for about 10⁶ time steps, which, depending on global parameters such as T_eff and log g, translates into real evolutionary times between a few months to a few dozen years. On a parallel cluster, one such model can be completed in a few days’ time when 16 processors (about 500 CPU hours) are used, where the chemical equilibrium solver GGCHEM is called a few 10⁹ times.

5 Results

5.1 Comparison to our previous cloud formation model

In the previous Helling & Woitke cloud formation models (Woitke & Helling 2004; Helling & Woitke 2006; Helling et al. 2008c), henceforth called the DRIFT models, the replenishment of elements was treated in a different way, using a prescribed mass-exchange timescale τ_mix(z) to replenish the atmosphere with fresh elements from the deep as $n_{\langle \text{H}\rangle }({\varepsilon }_k^0-{\varepsilon }_k)/{\tau }_{\text{mix}}$. The mass-exchange timescale was approximated by a power law logτ_mix(z) = const − βlogp(z) with power-law index β = 2.2 to describe convectional overshoot, see Eq. (9) in Woitke & Helling (2004). This simple approach led to an ODE-system that can be solved within about 2 CPU minutes.

Figure 3 compares the results of a previous DRIFT model with the new diffusive model, henceforth called the DIFFUDRIFT model. Both approaches model seed formation, kinetic surface growth and evaporation of cloud particles, and gravitational settling in the same way, but they differ in the treatment of the mixing that enters the cloud formation and the element conservation equations. The underlying temperature and pressure structures for all models discussed in this paper are taken from the DRIFT-PHOENIX atmosphere grid (Dehn et al. 2007; Helling et al. 2008b; Witte et al. 2009, 2011). In Fig. 3 we selected a model with effective temperature T_eff= 1800 K, surface gravity log g = 3, and metallicity Z = 1 (i.e. solar abundances are assumed deep in the atmosphere). The DRIFT-PHOENIX models solve the complete 1D model atmosphere problem including convection, radiative transfer, and hydrostatic structure, coupled to our previous DRIFT model, where the cloud opacities are calculated by Mie and effective medium theory. The resulting atmosphericstructure were frozen for this study, that is, the feedback of the new cloud formation model on the (p, T)-structure was not included.

The chemical setup for this comparison has 16 elements (H, He, Li, C, N, O, Na, Mg, Si, Fe, Al, Ti, S, Cl, K, Ca), one nucleation species (TiO₂), 12 solid species (TiO₂[s], Al₂ O₃[s], CaTiO₃[s], Mg₂ SiO₄[s], MgSiO₃[s], SiO[s], SiO₂[s], Fe[s], FeO[s], MgO[s], FeS[s], Fe₂ O₃[s]), and 60 surface reactions. The molecular setup in the new models is not quite identical because the DRIFT model uses a previous version of the chemical equilibrium solver GGCHEM, which has been replaced by the latest version (Woitke et al. 2018)in the DIFFUDRIFT model. We used 189 molecules in DRIFT and 308 molecules in DIFFUDRIFT to find the molecular concentrations in chemical equilibrium. We do not see any substantial differences in molecular concentrations caused by this data update, unless the local temperature falls below about 400 K. The thermochemical data for theselected solids are not entirely identical either, but these differences are likewise not substantial because the local temperatures remain above 700 K in this test. We assumed the mixing powerlaw index to be β = β′ = 1 for both τ_mix(z) in DRIFT and D_gas(z) in DIFFUDRIFT, see Eq. (18) and Fig. 1, although the meaning of β and β′ is slightly different. We note that using β′ > β would likely produce results that are more similar to each other than those presented in this paper. The lower volume boundary for the size integration of the cloud particle moments was set to V_ℓ = 10 × V_TiO2, where V_TiO2 = 3.14 × 10⁻²³ cm³ is the assumed volume of one unit of solid TiO₂[s].

The resulting cloud structures, as predicted by our previous DRIFT and the new DIFFUDRIFT models, are compared in Fig. 3. The diffusive transport of condensable elements up into the high atmosphere with DIFFUDRIFT is much less efficient than the assumed replenishment in the DRIFT model. As these elements are slowly mixed upwards by diffusion, they can collide with existing cloud particles to condense on, and hence far fewer of these elements reach the high atmosphere where the nucleation takes place. This is the main difference between the DRIFT and the DIFFUDRIFT models. In the previous DRIFT models, the mixing process was assumed to take place instantly.

Fig. 3 Comparison of cloud formation models for Teff = 1800 K, log g =3, metallicity Z = 1, and β = β′ = 1. The previous DRIFT model is shown by the thick grey lines. Two DIFFUDRIFT models are overplotted assuming pure gas diffusion (dashed lines) and gas + particle diffusion (black solid lines). dust/gas = ρdL3 is the dust-to-gas mass ratio, nd = ρL0 the number density of cloud particles, ${\langle a^3\rangle }^{1/3}={(3L_3/(4\pi L_0))}^{1/3}$ the mass-mean particle radius, and $\langle {\text{v}}_{\text{dr}}\rangle ={\langle a^3\rangle }^{1/3}\sqrt{\pi }g{\rho }_{\text{d}}/(2\rho c_T)$ the corresponding drift velocity according to Eq. (6).

Cloud structure

Consequently, the new DIFFUDRIFT model is featured by up to five orders of magnitude lower nucleation rates (Fig. 3) and less cloud particles high in the atmosphere. At intermediate pressures (10⁻⁶... 10⁻³ bar), we find that the fewer cloud particles in the DIFFUDRIFT model grow quickly and reach particle sizes of about 10 μm at 1 mbar, whereas in the DRIFT model, the cloud particles remain smaller, about 0.3 μm, because there are so many of them. The growth of the cloud particles is limited by the amount of condensable elements available per particle, and therefore, this effect is expected.

The dust-to-gas mass ratio, ρ_d∕ρ_gas, increases more steeply in the DIFFUDRIFT model, but reaches about the same maximum of order 10⁻³ at 1 mbar as in the DRIFT model. Thus, overall, the cloud formation process is about equally effective, but the clouds are spatially more confined in the DIFFUDRIFT model, reaching up just a few scale heights above the cloud base.

Table 2 lists vertically integrated cloud column densities for the three models. We find values of a few milligrams of condensates per cm², where the DIFFUDRIFT model without cloud particle diffusion is found to be the most dusty one. When an order of magnitude estimate of cloud particle opacities (see Appendix E) is used, values between several 100 cm² g⁻¹ to several 1000 cm² g⁻¹ at λ= 1μm are expected, depending on material and particle size distribution, that is, a column density of 1 milligram of condensate per cm² roughly corresponds to an optical depth of one at λ = 1 μm. This impliesthat all three models discussed here have optically thick cloud layers.

The computation of more realistic cloud particle opacities will need to take into account the height-dependent material composition, size distribution, and possibly the shape of the cloud particles, as done, for example, by Dehn et al. (2007), Witte et al. (2009, 2011) and Helling et al. (2019). In comparison to the previous DRIFT models, the particles in the DIFFUDRIFT models are larger, which is likely to cause the optical depths to be somewhat smaller, although the cloud mass column densities are similar. In addition, molecular opacities need to be added to calculate the spectral appearance of the objects and feedback onto the (p, T)-structure, which goes beyond the scope of the present paper.

The resulting particle properties in the main cloud layer below 1 mbar depend not only on the treatment of mixing, but also on whether we switch on the dust diffusion in the DIFFUDRIFT model. In this region, the cloud particles stepwise purify chemically as they decent in the atmosphere (see Fig. 5). Thermodynamically less stable materials such as Fe₂ O₃[s], FeO[s], and FeS[s] sublimate off the cloud particles sooner. Subsequently, the abundant magnesium-silicates MgSiO₃[s] and Mg₂SiO4[s] sublimate as well, which causes the cloud particles to shrink significantly around 1 mbar. Close to the cloud base, at about 1800 K in this model, only the most refractory materials remain, in particular metal oxides such as Al₂ O₃[s], CaTiO₃[s] and TiO₂[s], before even these materials eventually sublimate and the cloud particles evaporate completely.

As one material sublimates, the liberated elements may re-condense into different condensates, which are thermodynamically more stable, leading to rapid changes in particle size and material composition. There is also a dynamical effect. When the particles shrink, their fall speeds decrease, which leads to spatial accumulation, hence the number density of cloud particles n_d increases. While these effects and the general behaviour of the cloud particles are similar in all three models, the steps of sublimation are more pronounced in the DIFFUDRIFT model without dust diffusion. Dust diffusion tends to smooth out the variations of particle size and density.

Element abundances

Figure 4 compares the resulting gas element abundances. We see a strong depletion of condensable elements in the main cloud layer in all three models, by up to five orders of magnitude, concerning elements Ti, Al, Mg, Si, Mg, and Fe. However, the details are different. The previous DRIFT model is featured by minimums of ɛ_k that are similar in depth as compared tothe overall decrease of ɛ_k in the DIFFUDRIFT models. High up in the atmosphere, where cloud particles are virtually absent, there is no surface to condense on, and so the instantaneous mixing assumption in the DRIFT models causes a re-increase of ɛ_k toward the top of the atmosphere, unless the element can form nuclei. In the extremely low density gas at these heights, the nuclei simply fall through the atmosphere without much interaction, whereas elements that cannot nucleate, accumulate.

In contrast, in the new DIFFUDRIFT models, the abundances of all elements involved in cloud formation decrease with height in a monotonic way. This behaviour is expected in the final, time-independent, relaxed state of the atmosphere, as discussed in Sect. 3.1. In the stationary case, the downward transport of condensable elements via the falling cloud particles must be compensated for by an upward diffusive transport of these elements in the surroundinggas, which implies negative element abundance gradients throughout the cloudy atmosphere, see Eq. (29). The total drop of element abundances is deepest for Ti, but more shallow for Si and Fe than for Mg.

Freytag et al. (2010) performed 2D radiation hydrodynamical simulations of substellar atmospheres that included a time-dependent description for the formation of a single type of cloud particles for a fixed concentration of seed particles. The authors discussed substellar objects with T_eff = 900−2800 K, log(g) = 5 and solar element abundances. Their Fig. 9 (bottom left panel) shows the fraction of condensable gas in the atmosphere as a function of pressure, very similar to our Fig. 4 (lower plot). These results of Freytag et al. (2010) support our new DIFFUDRIFT results, where abundances of condensable elements in the gas phase decrease fast in the cloud layers and remain about constant above the clouds. We note that Freytag et al. (2010) prescribed the number of seed particles and considered only one generic condensate in their simulations.

Fig. 4 Effect of our assumptions on the mixing processes in the atmosphere on the resulting gas element abundances, in models with the same parameters as in Fig. 3. Upper plot: ${\varepsilon }_k^{\text{gas}}$ for instantaneous mixing as assumed in the previous DRIFT model. Lower plot: results according to the new DIFFUDRIFT models, where the full lines show the results for gas and dust diffusion, and the dashed lines show the results for gas diffusion only.

Fig. 5 Material volume composition of the cloud particles $b^s=L_3^{\text{s}}/L_3=V_s/V_{\text{tot}}$ for the samethree models as discussed in Figs. 3 and 4.

Cloud particle composition

Figure 5 shows the corresponding solid material composition (by volume) of the cloud particles. In all three models we assumed that the cloud particles have a well-mixed material composition that is independent of size, but that composition changes as the particles fall through the atmosphere, hence material composition depends on height. All three models show the same basic vertical structure.

• 1

  A layer containing only the most stable metal-oxides at the cloud base, in this model Al₂ O₃[s], TiO₂[s], and CaTiO₃[s]. The position of the cloud base, which depends on T_eff and log g, is located ataround 1800 K in this model.

• 2

  A thin layer of cloud particles around 1500 K that mainly consist of metallic Fe[s].

• 3

  Main silicate cloud layer composed of Mg₂ SiO₄[s], MgSiO₃[s], MgO[s], SiO[s], and SiO₂[s], mixed with metallic iron, upward of about 1400 K in this model.

• 4

  Lessstable solid materials such as FeS[s], FeO[s], and Fe₂ O₃[s] are incorporated into the silicate cloud particles at temperatures lower than about 1100 K, 1000 K, and 850 K, respectively, in this model.

• 5

  Pure nuclei at the top, here TiO₂[s], which fall through the atmosphere so quickly that they practically do not grow.

Further inspection shows, however, that the material composition of the main silicate cloud layer (3) differs substantially between the DRIFT and the DIFFUDRIFT models. In the new diffusive DIFFUDRIFT models, the first magnesium-silicate to form is fosterite Mg₂ SiO₄[s], which has a stoichiometric ratio Mg : Si = 2 : 1. The formation process of Mg₂ SiO₄[s] stops when the reservoir of Mg is exhausted, still leaving about half of the available Si in the gas phase. Because the mixing is diffusive, very little Mg can be mixed upwards through these Mg₂ SiO₄[s] clouds. Thus, the remaining amount of Si above the Mg₂ SiO₄[s] clouds preferentially forms other silicate materials, in particular SiO₂[s] and SiO[s], but only very little MgSiO₃[s]. This is different in our previous DRIFT model, where the depleted elements are assumed to be instantly replenished at similar rates, in which case both Mg₂ SiO₄[s] and MgSiO₃[s] are found to be about equally abundant condensates in the main silicate cloud layer.

Another difference concerns FeS[s] (troilite). FeS[s] is found to form in large quantities in the previous DRIFT models, causing ɛ_S to drop significantly, see the upper part of Fig. 4. However, this depends on our assumptions about how the elements are replenished. In the new diffusive models, upward mixing of gaseous Fe is rather inefficient because the Fe atoms have many opportunities to condense in form of Fe[s] on existing cloud particles along their way upwards in the atmosphere. When the temperature is low enough to allow FeS[s] to form, so little Fe is left in the gas phase that the S abundance is more or less unaffected by FeS[s] formation, and therefore sulphur remains available for other condensates to form.

5.2 Cloud structures as a function of T_eff

In this section, we study the results of a sequence of the new DIFFUDRIFT cloud formation models with decreasing effective temperature T_eff. We used a slightly different chemical setup here that allows us to discuss secondary cloud layers. We considered four nucleation species ${({\text{TiO}}_{\text{2}})}_N$, ${(\text{SiO})}_N$, ${(\text{KCl})}_N$, and ${(\text{C})}_N$. The nucleation rates of TiO₂, KCl, and C were calculated by modified classical nucleation theory (Helling et al. 2017; Gail et al. 1984), with a surface tension value for KCl from Lee et al. (2018). The nucleation rate of SiO is calculated according to Gail et al. (2013). We included 16 elements in this setup (H, He, Li, C, N, O, Na, Mg, Si, Fe, Al, Ti, S, Cl, K, and Ca), 14 condensed species (TiO₂[s], Al₂ O₃[s], MgSiO₃[s], Mg₂ SiO₄[s], SiO[s], SiO₂[s], Fe[s], FeS[s], FeO[s], MgO[s], KCl[s], NaCl[s], Na₂S[s], and C[s]), 308 molecules and 50 surface growth reactions. Molecular equilibrium constants and Gibbs free energies of the condensates were all taken from Woitke et al. (2018). Dust diffusion was included in all models.

Figure 6 shows the total column densities of selected cloud materials Σ_s [g cm⁻²] in a series of DIFFUDRIFT models with constant logg = 3 and mixing index β′ = 1, but decreasing T_eff. The column densities of the condensed species were computed as $${\Sigma }_s={\displaystyle \int {\rho }_s}\rho L_3^s\text{\hspace{0.17em}d}z,$$(35)

where ρ_s [g cm⁻³] is the material density of the pure condensate of kind s and $\rho L_3^s[{\text{cm}}^{\text{3}}{\text{cm}}^{\text{-3}}]$ is the volume of condensed kind s per volume of atmosphere. For example, for T_eff = 1800 K we find about 10 mg condensates per square centimetre, mostly made of Mg₂SiO₄[s], Fe[s] and Al₂O₃[s], followed by SiO[s] and MgSiO₃[s].

On the left side of this plot, the first model that shows condensation appears at T_eff= 2800 K. Here, the temperatures are too high for condensates other than just the most stable metal oxides in form of Al₂O₃[s] and TiO₂[s]. In the next few models down to T_eff = 2000 K, the main silicate layer forms, mixed with iron. In this range of effective temperatures, Al₂O₃[s] still has the highest column density because the metal-oxide layer is situated deeper in the atmosphere where the densities are higher. Only for T_eff < 2000 K does the silicate-iron layer start to dominate by mass. At the very end of the sequence, for T_eff < 1500 K we find the first models that host a third cloud layer made of di-sodium sulfide Na₂S[s].

Figure 7 shows a few more details from this T_eff series of new cloud formation models. The upper left plot shows the atmospheric density and temperature structures we assumed, taken from the DRIFT-PHOENIX atmosphere grid (Dehn et al. 2007; Helling et al. 2008b; Witte et al. 2009, 2011). The kinks in deep layers (T ~ 2500−3000 K) indicate the beginning of the convective layer (Schwarzschild criterion for convective instability). The upper right plot shows the assumed diffusion coefficient in the atmosphere, which decreases with T_eff, because the convective layer sinks into deeper layers, hence the spatial distance to the source, which causes the mixing motions in the atmosphere, increases.

The left middle plot in Fig. 7 shows the dust-to-gas mass ratio, which has its maximum in the main silicate-iron layer, and a shoulder on the right caused by the deeper metal-oxide clouds, which are made of the rarer elements with the highest condensation temperatures: aluminium, calcium, and titanium. As T_eff decreases, both layers move inward to deeper layers and become successively more narrow, until finally, for T_eff = 1400 K, a new cloud layer occurs that mainly consists of di-sodium sulfide Na₂S[s]. The right middle plot shows the effect on the silicon abundance in the gas phase. All curves decrease monotonically towards the surface, with higher Si depletions for lower T_eff where the silicate cloud particle formation is more complete.

The nucleation rates of ${({\text{TiO}}_{\text{2}})}_{\text{N}}$ and ${(\text{SiO})}_{\text{N}}$ particles are depicted in the lower left plot. A complicated double-peaked pattern emerges that has a minimum around the main peak of the dust-to-gas ratio (at the peak position of the main silicate-iron layer). ${({\text{TiO}}_{\text{2}})}_N$ is usually the most significant nucleation species, but cooler models show additional contributions by ${(\text{SiO})}_N$. The resulting mean particle sizes are plotted on the lower right, with a tendency to produce larger particles deep in the atmosphere for lower T_eff. An in between minimum in particle size occurs where the main silicate material evaporates. Only the coolest model has a second minimum where Na₂S[s] evaporates. Interestingly, the hottest and coolest models in Fig. 7 show about equally large cloud particles at high altitudes, whereas all other models show smaller particles.

Fig. 6 Column densities [g cm−2] of different condensates in the atmosphere for a sequence of models with decreasing Teff but constant logg = 3 and β′ = 1. A value of 10−3 g cm−2 roughly corresponds to an optical depth of one at a wavelength of λ = 1 μm (see Appendix E).

6 Summary and discussion

This paper has introduced a new cloud formation model that is applicable to the atmospheres of brown dwarfs and gas giant (exo-)planets. We combined our previous cloud particle moment method (Woitke & Helling 2004; Helling & Woitke 2006; Helling et al. 2008c) with a diffusive mixing approach, according to which fresh condensable elements are diffusively transported upwards to replenish the upper atmosphere via a combination of turbulent (eddy) mixing and gas-kinetic diffusion in the final relaxed time-independent state of the atmosphere. Our formulation of the problemarrives at a system of about 30 second-order partial differential equations of reaction-diffusion type, where the formation and growth of the cloud particles follows from a kinetic treatment in phase-non-equilibrium.

Model setup

The new cloud formation model was applied to a given 1D (p, T) atmosphericstructure in this paper. The model was calculated time-dependently, using an operator splitting technique. All models are found to relax toward a time-independent stationary solution, where the condensable elements are constantly mixed up diffusively, cloud particles nucleate from the gas phase high in the atmosphere, grow by the simultaneous condensation of different solid materials on their surface, and then descent through the atmosphere due to gravitational settling, before the particles stepwise purify and eventually sublimate completely at the cloud base.

Timescales

The real-time simulation time required to reach this stationary solution varies between a few months to several dozen years, depending on log(g) and T_eff. The relaxation is quicker when models are started from an atmosphere that is devoid of any condensable elements at t = 0. These relatively long simulation times make the models computationally expensive (about 500 CPU-hours per model) because the intrinsic nucleation and growth reactions are very fast, which means that the models need to be advanced on short computational time steps of the order of seconds to guarantee numerical stability. The long physical timescales involved in the simulations are (i) the overall settling time for small particles inserted high in the atmosphere, and (ii) the overall mixing time for gas parcels to diffusively reach the highest point in the model from the cloud base. This implies that 3D simulations of cloud formation (GCM models, e.g. Freytag et al. 2010; Lee et al. 2016; Lines et al. 2018; Powell et al. 2018; Charnay et al. 2018) must be advanced for similar real-time simulation times before a relaxed physical structure can be expected. However, how long these physical timescales actually are will depend on the exact formulation of mixing and setting in the GCM models.

Fig. 7 Sequence of cloud forming models with decreasing Teff at constant logg = 3 and mixing power-law index β′ = 1. Top row: gastemperature T and diffusioncoefficient Dgas as a functionof pressure (both assumed). Middle row: resulting dust-to-gas mass ratio and element abundance of silicon in the gas phase ɛSi. Lower row: resulting nucleation rates J⋆ and mean particle sizes ${\langle a^3\rangle }^{1/3}$.

Cloud density and particle sizes

In comparison to our previous DRIFT models, the DIFFUDRIFT models show fewer but larger cloud particles, which are more concentrated towards the cloud base. However, the physical properties of the cloud particles in the main silicate-iron layer towards the bottom of the clouds (dust-to-gas mass ratio, particle sizes, optical depth, chemical composition, etc.) are found to be similar to the results of the previous models. The dust-to-gas ratio in the main silicate-iron layer reaches a peak value of about 0.002 to 0.003, quite independent of T_eff, for moderately hot models (T_eff > 2500 K). This is close to the maximum value of 0.0045 that is expected from complete condensation of a gas with solar abundances (Woitke et al. 2018). The physical reason for the stronger concentration of the cloud particles around the cloud base in the DIFFUDRIFT models is that the diffusive element replenishment is less effective for the upper atmosphere because the molecules carrying the elements diffusively upwards have a high probability of colliding with existing cloud particles on their way up the atmosphere. This effect was not accounted for in our previous models.

Element abundances

The concentration of condensable elements in the gas phase shows a steep decline in the DIFFUDRIFT models above the cloud base, followed by a monotonous decrease towards a plateau that continues on that level toward the upper boundary of the model. This behaviour is expected in the time-independent relaxed state because the downward flux of condensable elements via the falling cloud particles must be compensated for by an upward diffusive flux of elements in the gas phase. This in turn requires a negative concentration gradient. We find that the abundances of the condensable elements high above the cloud layers strongly depend on effective temperature, in agreement with the results of 2D radiation hydro-models by Freytag et al. (2010). For example, the silicon abundance is reduced by about 2.5 orders of magnitude for T_eff = 2000 K, but by 6 orders of magnitude for T_eff = 1300 K in our models.

Fig. 8 Concentration of condensed species in a model with Teff = 1300 K, log g =3, and β′ = 1, showing a secondary cloud layer almost entirely made of di-sodium sulfide Na2 S[s]. $n_{\text{cond}}^s=\rho L_3^s/V_{\text{mat}}^s[{\text{cm}}^{\text{-3}}]$ is the number density of solid units of condensed species s, $V_{\text{mat}}^s=m_s/{\rho }_s$ is the volume occupied by one unit of solid s in the pure material, and ms is the mass of one such unit, for example 100.4 amu for s = MgSiO3. $n_{\text{cond}}^s/n_{\langle \text{H}\rangle }$ is directly comparable to element abundances.

Cloud material composition

The chemical composition of the cloud particles is characterised by (i) a deep layer with only the most stable metal oxides at the cloud base (Al₂ O₃[s], TiO₂[s], and CaTiO₃[s] in the current setup), (ii) the main silicate-iron layer (Mg₂SiO₄[s] and Fe[s]), which is enriched by other silicates such as SiO[s] and SiO₂[s] and MgSiO₃[s] at increasing heights. Some less stable condensates are also found in smaller quantities, in particular FeS[s] and FeO[s]. The condensation in both these cloud layers leads to a removal of certain elements from the gas phase, and the stoichiometry of the condensates decides upon which elements remain for further condensation higher in the atmosphere. In particular, the formation of Mg₂ SiO₄[s], with stoichiometry Mg : Si = 2 : 1, causes the Mg abundance to drop quickly, whereas roughly half of the originally available Si remains in the gas phase. This then favours the formation of SiO[s] and SiO₂[s] above the Mg₂ SiO₄[s] layer, rather than the formation of MgSiO₃[s], which is a rather unimportant condensate in our new DIFFUDRIFT models. The low abundance of MgSiO₃[s] in the main silicate-iron layer is a result that differs from the results obtained with our previous DRIFT models, and from phase-equilibrium models starting from complete solar abundances (Woitke et al. 2018).

Na₂S[s] clouds

Our coolest DIFFUDRIFT models show the occurrence of a secondary cloud layer almost entirely made of di-sodium sulfide Na₂ S[s], see Fig. 8. The presence of Na₂S-clouds in brown dwarf atmospheres has been proposed by Morley et al. (2012) to fit the optical appearance of two red T dwarfs. The formation of Na₂ S-clouds requires the presence of sulphur and sodium in the gas phase at low temperatures. In phase-equilibrium models starting from solar abundances,this combination is prevented by the formation of FeS[s] (troilite), which consumes the sulphur. However, in our new diffusive kinetic cloud formation models, iron is depleted by the formation of metallic iron Fe[s] at high temperatures, so that FeS[s] cannot form in large quantities. Consequently, sulphur remains available to eventually form Na₂ S[s] at lower temperatures. The condensation of Na₂S[s] then reduces the possibility of forming NaCl[s] at even lower temperatures, and so on. Therefore, the new diffusive DIFFUDRIFT models reveal new details about the condensation sequence in cloudy atmospheres, and we need more experiments with our selection of condensates during model initialisation to arrive at more distinct conclusions.

7 Conclusions

The physical description of the replenishment mechanism for condensable elements in planetary atmospheres seems crucial for realistic cloud formation models. This paper has used a quasi-diffusive approach in 1D to simulate the turbulent eddy-mixing processes in cloudy atmospheres, using the new DIFFUDRIFT models. This approach can be considered as the limiting case of small-scale mixing. At the other extreme, large-scale hydrodynamic motions (convection, Hadley cells, etc.) might be able to dredge up these elements in a possibly more direct way, which was the idea in our previous DRIFT models. In reality, there is not only vertical, but also horizontal mixing, which is likely to be quite fast, for example in super-rotating horizontal jets as known from Jupiter (Schneider & Liu 2009). However, horizontal mixing alone cannot couterbalance vertical settling on the long run, it rather stimulates vertical mixing via e.g. Kelvin-Helmholtz instabilities. More 3D numerical experiments are required to quantify the efficiency of mixing to inform our cloud formation models.

Acknowledgements

Ch.H. thanks Will Best and Jonathan Gagné for the discussion on the number of known brown dwarfs. We thank Robin Baeyens and Ludmila Carone for insightful discussions on mixing regimes in gas giant planets. The computer simulations were carried out on the UK MHD Consortium parallel computer at the University of St Andrews, funded jointly by STFC and SRIF.

Appendix A Diffusionsolver

We used a self-developed 1D explicit/implicit diffusion solver in this paper that has second-order accuracy in both the formulation of the differential equations and the boundary conditions³. In case of a plane-parallel static atmosphere, the diffusion equation for element k is given by $$\frac{\text{d}(n_{\langle \text{H}\rangle }{\varepsilon }_k)}{\text{d}t}=\frac{\text{d}}{\text{d}z}\left(n_{\langle \text{H}\rangle }{D}_{\text{gas}}\frac{\text{d}{\varepsilon }_k}{\text{d}z}\right),$$(A.1)

where the diffusive element flux is $$j_k^{\text{diff}}=-n_{\langle \text{H}\rangle }{D}_{\text{gas}}\frac{\text{d}{\varepsilon }_k}{\text{d}z}.$$(A.2)

A.1 Vertical grid and discretisation of derivatives

We introduce an ascending vertical grid z_i (i = 1,..., I). The first and second derivatives of any quantity f(z_i) = f_i at grid point z_i are approximated as $$\begin{array}{lll}\frac{\text{d}{f}_i}{\text{d}z}\hfill & =\hfill & d_i^{l,1}{f}_{i-1}+d_i^{m,1}{f}_i+d_i^{r,1}{f}_{i+1}\hfill \\ \frac{d^2{f}_i}{\text{d}{z}^2}\hfill & =\hfill & d_i^{l,2}{f}_{i-1}+d_i^{m,2}{f}_i+d_i^{r,2}{f}_{i+1},\hfill \end{array}$$

that is, as linear combinations of the function values on the neighbouring grid points, where $d_i^{l,1}$, for example, is the coefficient for the first derivative on the point left of the grid point i, $d_i^{m,1}$ the same on the mid-point, and $d_i^{r,1}$ the same on the point right of grid point i. Similar, for the second derivative, the coefficients are $d_i^{l,2}$, $d_i^{m,2}$ and $d_i^{r,2}$. When we use a second-order polynomial approximation for function f(z), the coefficients are given by $$\begin{array}{lll}{d}_i^{l,1}\hfill & =\hfill & -\frac{h_i^r}{(h_i^r+h_i^l)h_i^l}\hfill \\ d_i^{m,1}\hfill & =\hfill & +\frac{h_i^r-h_i^l}{h_i^l{h}_i^r}\hfill \\ d_i^{r,1}\hfill & =\hfill & +\frac{h_i^l}{(h^r(i)+h_i^l)h_i^r}\hfill \\ d_i^{l,2}\hfill & =\hfill & +\frac{2}{(h_i^r+h_i^l)h_i^l}\hfill \\ d_i^{m,2}\hfill & =\hfill & -\frac{2}{h_i^r{h}_i^l}\hfill \\ d_i^{r,2}\hfill & =\hfill & +\frac{2}{(h_i^r+h_i^l)h_i^r},\hfill \end{array}$$

where $h_i^l=z_i-z_{i-1}$ and $h_i^r=z_{i+1}-z_i$ are the left- and right-hand side grid point distances. For the special case of an equidistant grid, we have $h=h_i^l=h_i^r$ and hence $$\begin{array}{lll}\frac{\text{d}{f}_i}{\text{d}z}\hfill & =\hfill & \frac{f_{i+1}-f_{i-1}}{2h}\hfill \\ \frac{d^2{f}_i}{\text{d}{z}^2}\hfill & =\hfill & \frac{f_{i+1}-2f_i+f_{i-1}}{h^2}.\hfill \end{array}$$

The above equations are valid for grid points i = 2,..., I − 1. For the first derivative at the boundaries, we write $$\begin{array}{lll}\frac{\text{d}{f}_1}{\text{d}z}\hfill & =\hfill & d_1^{l,1}{f}_1+d_1^{m,1}{f}_2+d_1^{r,1}{f}_3\hfill \\ \frac{\text{d}{f}_1}{\text{d}z}\hfill & =\hfill & d_I^{l,1}{f}_{I-2}+d_I^{m,1}{f}_{I-1}+d_I^{r,1}{f}_I,\hfill \end{array}$$

which is also second-order accuracy by using the information on the three leftmost or three rightmost grid points, respectively. The coefficients are given by $$\begin{array}{lll}{d}_1^{l,1}\hfill & =\hfill & -\frac{h_2+h_3}{h_2{h}_3}\hfill \\ d_1^{m,1}\hfill & =\hfill & \frac{h_3}{h_2(h_3-h_2)}\hfill \\ d_1^{r,1}\hfill & =\hfill & -\frac{h_2}{h_3(h_3-h_2)}\hfill \\ d_I^{r,1}\hfill & =\hfill & \frac{h_{I-1}+h_{I-2}}{h_{I-1}{h}_{I-2}}\hfill \\ d_I^{m,1}\hfill & =\hfill & -\frac{h_{I-2}}{h_{I-1}(h_{I-2}-h_{I-1})}\hfill \\ d_I^{l,1}\hfill & =\hfill & \frac{h_{I-1}}{h_{I-2}(h_{I-2}-h_{I-1})},\hfill \end{array}$$

where h₂ = z₂ − z₁, h₃ = z₃ − z₁, h_I−1 = z_I − z_I−1, and h_I−2 = z_I − z_I−2.

A.2 Spatial derivatives

The diffusion term at grid point z_i (i = 2... I − 1) is numerically resolved, with abbreviation D_gas(z_i) = D_i, as $$\begin{array}{l}{\frac{\text{d}}{\text{d}z}\left(n_{\langle \text{H}\rangle }{D}_{\text{gas}}\frac{\text{d}{\varepsilon }_k}{\text{d}z}\right)|}_{z_i}={\frac{\text{d}(n_{\langle \text{H}\rangle }{D}_{\text{gas}})}{\text{d}z}|}_{z_i}\cdot {\frac{\text{d}{\varepsilon }_k}{\text{d}z}|}_{z_i}+{n_{\langle \text{H}\rangle }{D}_{\text{gas}}\frac{{\text{d}}^2{\varepsilon }_k}{\text{d}{z}^2}|}_{z_i}\\ =\left(d_i^{l,1}\text{nHim}1D_{i-1}+d_i^{m,1}{n}_{\langle \text{H}\rangle ,\text{i}}{D}_i+d_i^{r,1}\text{nHip}1D_{i+1}\right)\\ \cdot \left(d_i^{l,1}{\varepsilon }_{k,i-1}+d_i^{m,1}{\varepsilon }_{k,i}+d_i^{r,1}{\varepsilon }_{k,i+1}\right)\\ +n_{\langle \text{H}\rangle ,\text{i}}{D}_i\left(d_i^{l,2}{\varepsilon }_{k,i-1}+d_i^{m,2}{\varepsilon }_{k,i}+d_i^{r,2}{\varepsilon }_{k,i+1}\right),\end{array}$$(A.21)

and the diffusive fluxes across the lower and upper boundaries are $$\begin{array}{lll}{\varphi }_{k,1}\hfill & =\hfill & -D_1{n}_{\langle \text{H}\rangle ,\text{1}}\left(d_1^{l,1}{\varepsilon }_{k,1}+d_1^{m,1}{\varepsilon }_{k,2}+d_1^{r,1}{\varepsilon }_{k,3}\right)\hfill \\ {\varphi }_{k,I}\hfill & =\hfill & -D_I{n}_{\langle \text{H}\rangle ,I}\left(d_I^{l,1}{\varepsilon }_{k,I-2}+d_I^{m,1}{\varepsilon }_{k,I-1}+d_I^{r,1}{\varepsilon }_{k,I}\right).\hfill \end{array}$$

A.3 Boundary conditions

As boundary conditions, we implemented three options, for example considering the lower boundary:

• 1.

  fixed concentration: ɛ_k,1 is a given constant;

• 2.

  fixed flux: ϕ_k,1 is a given constant;

• 3.

  fixed outflow rate: the flux through a boundary is assumed to be proportional to the concentration of species k at the boundary, for instance,

$${\varphi }_{k,1}={\beta }_k{n}_{\langle \text{H}\rangle ,\text{1}}{\varepsilon }_{k,1}{\text{v}}_k[{\text{cm}}^{\text{-2}}{\text{s}}^{\text{-1}}],$$(A.24)

where β_k is a given probability (fixed value)and v_k is the speed at which element k moves through the boundary (also fixed value).

Fig. A.1 Test problem with fixed concentrations on the left and right boundaries (ɛ = 0). The thin black lines overplot the analytic solution ɛ(z, t) = exp(−ωt)sin(kz) with k = π and ω = Dk2.

Fig. A.2 Diffusive evolution of an initial δ peak with analytic solution overplotted. The analytic solution is $\varepsilon (z,t)=A(t)\exp (-{(z-0.5)}^2/w^2(t))$ with $A(t)=\sqrt{t_0/t}$ and $w(t)=2\sqrt{Dt}$.

A.4 Explicit integration

A straightforward way to integrate Eq. (A.1) for a time step Δt is the following explicit scheme: $$f_i^{(n)}=f_i^{(n-1)}+\Delta t\frac{\text{d}{f}_i^{(n-1)}}{\text{d}t},$$(A.25)

where $f_i^{(n)}$ is some quantity on grid point i at time tⁿ and $f_i^{(n-1)}$ is the quantity on grid point i at time tⁿ⁻¹ with tⁿ = tⁿ⁻¹ + Δt. In consideration of Eq. (A.1), this leads to $$\begin{array}{lll}{\varepsilon }_{k,i}^{(n)}\hfill & =\hfill & {\varepsilon }_{k,i}^{(n-1)}+\Delta t{\displaystyle \sum _{j=1}^I{A}_{ij}}{\varepsilon }_{k,j}^{(n-1)},\hfill \end{array}$$(A.26)

where A is a tri-diagonal matrix, the elements A_ij of which are given by Eq. (A.21). Equation (A.26) applies to the grid points i = 2,..., I − 1, but not to the boundaries. On the boundary points, the following equations are applied depending on the choice of boundary conditions, here for example the lower boundary,

• 1.

  fixed concentration: ${\varepsilon }_{k,1}^{(n)}={\varepsilon }_k^0$

• 2.

  fixed flux: ${\varepsilon }_{k,1}^{(n)}=\frac{1}{d_1^{l,1}}\left(-\frac{{\varphi }_{k,1}}{D_1{n}_{\langle \text{H}\rangle ,\text{1}}}-d_1^{m,1}{\varepsilon }_{k,2}^{(n)}-d_1^{r,1}{\varepsilon }_{k,3}^{(n)}\right)$

• 3.

  fixed outflow rate:

$$\begin{array}{l}{\beta }_k{n}_{⟨\text{H}⟩,1}{\varepsilon }_{k,1}{\text{V}}_k=-D_1{n}_{⟨\text{H}⟩,1}\left(d_1^{l,1}{\varepsilon }_{k,1}+d_1^{m,1}{\varepsilon }_{k,2}+d_1^{l,1}{\varepsilon }_{k,3}\right)\\ \Rightarrow {\varepsilon }_{k,1}^{\left(n\right)}=\frac{-d_1^{m,1}{\varepsilon }_{k,2}^{\left(n\right)}-d_1^{r,1}{\varepsilon }_{k,3}^{\left(n\right)}}{d_1^{l,1}+\frac{{\beta }_k{\text{V}}_k}{D_1}}\end{array}$$

These assignments are applied at time tⁿ, that is, after an explicit diffusion time step has been completed on grid points i = 2... I − 1. To guarantee numerical stability, the explicit time step must be limited by α ≤ 0.5 according to $$\Delta t=\alpha \cdot {\min }_{i=2,\mathrm{...},I}\frac{{(z_i-z_{i-1})}^2}{\frac{1}{2}(D_i+D_{i-1})}.$$(A.27)

A.5 Implicit integration

To avoid the time-step limitations in the explicit solver and to guarantee numerical stability for much larger time steps, an implicit integration scheme can optionally be applied, $$\begin{array}{lll}{f}_i^{(n)}\hfill & =\hfill & f_i^{(n-1)}+\Delta t\frac{\text{d}{f}_i^{(n)}}{\text{d}t},\hfill \end{array}$$(A.28)

which is a system of linear equations for the unknowns $f_i^{(n)}$. In consideration of Eq. (A.1), we have $$\begin{array}{lll}{\varepsilon }_{k,i}^{(n)}\hfill & =\hfill & {\varepsilon }_{k,i}^{(n-1)}+\Delta t{\displaystyle \sum _{j=1}^I{A}_{ij}}{\varepsilon }_{k,j}^{(n)}.\hfill \end{array}$$(A.29)

We re-write this equation more generally, including the boundary conditions, by means of another unit-free matrix as $$\begin{array}{lll}B{\varepsilon }_k^{(n)}\hfill & =\hfill & R_k,\hfill \end{array}$$(A.30)

where we have $$B_{ij}={\left(\text{mathds}1-\Delta tA\right)}_{ij}\text{and}{R}_{k,i}={\varepsilon }_{k,i}^{(n-1)}\text{mbox}for\$i=2\mathrm{...}I-1\$,$$(A.31)

where some of the B_ij and R_k,i depend on the selected boundary conditions, for example at the lower boundary,

• 1.

  fixed concentration: B₁₁ = 1 and $R_{k,1}={\varepsilon }_k^0$

• 2.

  fixed flux: B₁₁ = 1, $B_{12}=d_1^{m,1}/d_1^{l,1}$, $B_{13}=d_1^{m,1}/d_1^{l,1}$, and $R_{k,1}=-{\varphi }_{k,1}/(n_{\langle \text{H}\rangle ,\text{1}}{D}_1{d}_1^{l,1})$

• 3.

  fixed outflow rate: $B_{11}=1+{\beta }_k{\text{v}}_k/(D_1{d}_1^{l,1})$, $B_{12}=d_1^{m,1}/d_1^{l,1}$, $B_{13}=d_1^{m,1}/d_1^{l,1}$, and R_k,1 = 0.

We can now perform an implicit time step according to Eq. (A.30) as $${\varepsilon }_k^{(n)}=B^{-1}{R}_k,$$(A.32)

where B⁻¹ is the inverse of the matrix B. As long as the spatial grid points z_i, the densities n_{⟨H ⟩,i} and diffusion constants D_i, the constants involved in the boundary conditions (e.g. ϕ_k,1 or β_k), and the time step Δt do not change, we need to perform the matrix inversion only once. Successive time steps are then performed by simply incrementing n, re-computing the vector R_k, and applying again Eq. (A.32). B⁻¹ is also usually the same for all elements k to be diffused.

Fig. A.3 Testing the model convergence for different initial conditions: gas element abundances ɛk as a functionof pressure p in two cloud formation models with different initial conditions, “empty” (left) and “full” (right), see text for further explanations. Panels from top to bottom: respective results after t = 0, t = 10 days, and t = 8000 days.

This favourable property of B makes the computation of implicit time steps very fast. We note, however, that B⁻¹ is in general a full I × I matrix where all entries are positive ${(B^{-1})}_{ij}>0$. This leads to a very stable numerical behaviour for arbitrary time steps. In contrast, the matrix A has positive entries along the main diagonal, but negative entries along both semi-diagonals, which leads to numerical instabilities when the time step is too large.

Appendix B Settling solver

For the 1D vertical settling in the Epstein regime we solve $$\frac{\text{d}(\rho L_j)}{\text{d}t}=\frac{\text{d}}{\text{d}z}{\displaystyle {\int }_{V_{\text{}\mathcal{l}}}{v}_{\text{dr}}^{{}^{°}}}(V)f(V)V^{j/3}\text{\hspace{0.17em}d}V=\frac{\text{d}}{\text{d}z}\left(\xi \frac{{\rho }_{\text{d}}}{c_T}{L}_{j+1}\right),$$(B.1)

according to Eqs. (6) and (8). The settling flux for cloud particle moment ρL_j is hence $${\varphi }_j=-\xi \frac{{\rho }_{\text{d}}}{c_T}{L}_{j+1}=-\rho L_j{\text{v}}_{\text{dr},j},$$(B.2)