Settling of inertial particles in turbulent Rayleigh-Benard convection

The settling behaviour of small inertial particles in turbulent convection is a fundamental problem across several disciplines, from geophysics to metallurgy. In a geophysical context, the settling of dense crystals controls the mode of solidification of magma chambers and planetary-scale magma oceans, while rising of light bubbles of volatiles drives volcanic outgassing and the formation of primordial atmospheres. Motivated by these geophysical systems, we perform a systematic numerical study on the settling rate of particles in a rectangular two-dimensional Rayleigh-Benard system with Rayleigh number up to 10^12 and Prandtl number from 10 to 50. Under the idealized condition of spherically-shaped particles with small Reynolds number, two limiting behaviours exist for the settling velocity. On the one hand, Stokes' law applies to particles with small but finite response time, leading to a constant settling rate. On the other hand, particles with a vanishing response time are expected to settle at an exponential rate. Based on our simulations, we present a new physical model that bridges the gap between the above limiting behaviours by describing the sedimentation of inertial particles as a random process with two key components: i) the transport of particles from vigorously convecting regions into sluggish, low-velocity"piles"that naturally develop at the horizontal boundaries of the system, and ii) the probability that particles escape such low-velocity regions without settling at their base. In addition, we identify four distinct settling regimes and analyze the horizontal distribution of sedimented particles. For two of these regimes settling is particularly slow and the distribution is strongly non-uniform, with dense particles being deposited preferentially below major clusters of upwellings.


I. INTRODUCTION
The motion and sedimentation of small particles in a convecting fluid is of great interest to fluid dynamicists, geophysicists, as well as to metallurgists. Dispersion of pollutants, dust [1], and organic material such as pollen in the atmosphere [2], phytoplankton population dynamics in oceans and lakes [3,4], crystallization of magma chambers [5,6], of primordial magma oceans in rocky planets [e.g. 7,8], and of planetary metallic cores [e.g. 9,10] are just a few examples of geophysical processes that to some extent can be described through the settling of inertial particles in a fluid undergoing highly vigorous convection. In industrial settings, the purification of alloys [11] and microfluidic heat transfer technologies [e.g. 12] count into this group.
From a general point of view, the phenomenology of particle-laden turbulent flows has been the subject of extensive studies over the last decades. Researches have focused on the statistical characterisation of particles' dispersion and accumulation upon varying the flow turbulence intensity, the particles' inertia (which is linked to their size and mass density), and their shape [13,14]. For what concerns particle turbulent settling, there are nowadays solid numerical and experimental evidences that particle sinking is enhanced by the effect of turbulence, while the opposite happens for particle rising [for a review, see 15]. However, these results have been obtained in the context of idealized flows: kinematic flows [16,17], unbounded turbulence [18], or channel and pipe flow geometries. Much less explored is the context of thermally-driven flows [for experimental studies, see 5,19].
Despite the relevance of the problem of particles settling in turbulent flows for a variety of natural and industrial systems, this work is largely motivated by the study of the crystallization of primordial magma oceans, with our choice of the parameter space being inspired by this system (Section II A). Planetary-scale volumes of liquid silicates are thought to form during the accretion and differentiation of terrestrial bodies like the Earth, Mercury, Venus, Mars and the Moon [e.g., [20][21][22]. The way magma oceans solidify is of primary importance for the long-term thermo-chemical evolution of the interior of planets [e.g., 23]. Whether or not newly formed crystals settle or remain suspended by turbulent flow determines the initial distribution of the composition of the silicate mantle [e.g. , 7]. This is a difficult problem that depends on the density contrast between crystals and melt, the melt viscosity, the size of the crystals and the convective dynamics of the system. If dense crystals are efficiently maintained in suspension, a magma ocean undergoes equilibrium (or batch) crystallization, which leads to a largely homogeneous composition of the rocky mantle. By contrast, if crystals tend to settle and crystal-melt separation is efficient, fractional crystallization takes place. Residual melts are progressively more and more enriched in so-called incompatible elements such as iron-oxides and heat-producing elements. This process ultimately leads to a compositionally-stratified mantle whose longterm evolution can be dramatically different from that of a homogeneous one [e.g., [24][25][26][27]. Magma chambers are small-scale analogs of magma oceans. Upon cooling and solidification these typically undergo strong fractionation, which is often attributed to crystal settling [e.g., 5,6,28,29].
Similar to the settling of negatively buoyant crystals, floating of light particles is also a fundamental process in the context of the crystallization of magma oceans and magma chambers. In fact, it is the basic mechanism underlying magma degassing, where gas bubbles are released from volatile-saturated magma [e.g., 30]. Greenhouse volatiles such as H 2 O and CO 2 also behave as incompatible species and tend to be strongly enriched in the liquid phase upon magma crystallization. The efficiency with which these are released from a magma ocean controls the formation of primordial atmospheres and the timescale of magma ocean solidification [e.g. [31][32][33].
In the context of a crystallizing magma, Marsh and Maxey [28] modelled the transport of particles by convective motions as a turbulent diffusion process, which was a common approach in studies of mixing in turbulent flows [e.g. 34,35]. Fundamental laboratory experiments aimed at assessing settling rates in a cooling magma were later performed by Martin and Nokes [5], who employed the turbulent diffusion theory to explain their measurements. Assuming the concentration of particles to be spatially uniform, they derived a simple model for particles with a vanishing Stokes velocity according to which the number of suspended particles decays exponentially with time. Although Martin and Nokes [5] anticipated that for particles with a larger Stokes velocity their "diffusion model of turbulent transport will begin to break down and other assumptions will no longer be valid, in particular the assumption of one-dimensionality", surprisingly little effort has been devoted to extend their work. To our knowledge, no experimental or numerical study has been performed that systematically explores the settling mechanism of particles with a non-vanishing Stokes velocity in turbulent, thermally-driven convection.
Differentiation of a cooling magma is a competitive process between generation, sedimentation and reentrainment of crystals. The problem of re-entrainment, in particular, has been addressed by various authors both theoretically [e.g., 36] and experimentally [19,36]. Although the lifting of negatively buoyant particles from the crests of dunes has been recognized as one of the main mechanisms to keep these in suspension [36], in this study we focus entirely on the settling process and completely neglect re-entrainment, which we plan to address in future work. As soon as our particles reach the bottom boundary of the domain (resp. the top boundary for light particles), we eliminate them from the flow, not allowing any accumulation or subsequent lifting of the sedimented material.
We use a modeling approach based on a Eulerian-Lagrangian description of the fluid flow and the particulate phase respectively, and track individual trajectory of each particle. Our approach thus captures how the exact flow structure affects the particle motion. This level of detail is still challenging from the experimental point of view and it brings new results when compared to the one-dimensional turbulent diffusion theory. For example, horizontal variations in the distribution of sedimented particles can be evaluated and linked to the large scale circulation of the fluid (also called the "wind of turbulence", see e.g. Ahlers et al. [37] for a discussion).
Based on the experimental work of Koyaguchi et al. [6], Sparks et al. [38] argued for cyclic sedimentation of crystals in magma chambers caused by the cessation of convection due to the particle concentration exceeding a certain critical value. Similar behaviour was observed by Höink et al. [39] in the context of numerical simulations of metal-silicate separation, and by Verhoeven and Schmalzl [40] in a numerical study that combines a finite volume convection code with a discrete element method [41]. Here we neglect the influence of particles on the convective flow [e.g. 42], i.e. we assume only small particle concentrations. We believe that the dynamics of dilute suspensions is sufficiently rich to warrant a study entirely dedicated to particle settling before considering additional complexities arising from larger solid fractions.
For particles with small radii Verhoeven and Schmalzl [40] obtain statistically stationary suspension in which convective motions keep particles indefinitely entrained. Similar results were obtained in the nonrotating cases of Maas and Hansen [43,44] whose model builds on the one by Verhoeven and Schmalzl [40]. Although it is at odds with the early results of Martin and Nokes [5], Maas and Hansen [44] conclude that "it is generally assumed that vigorous convection would prevent major gravitational segregation in a magma ocean at all latitudes (Andrault et al., 2017)". Here we employ a more elaborate particle model and refute such statement, confirming the experimental results of Martin and Nokes [5] in which particles with a vanishing response time settle at an exponential rate.
The rest of the paper is organized as follows. In Section II we introduce our numerical model and discuss the choice of model parameters. In Section III we present the settling curves of a reference simulation and classify them according to four distinct regimes. In Section III D we then introduce a general model that describes particle settling as a random process. In Section IV A we discuss the horizontal distribution of sedimented particles, showing how it can be strongly non-uniform in some regimes. In Section IV B the focus is on particles lighter than the fluid, including bubbles. These become concentrated in flow vortices, which significantly delays their rising. In Section IV C we analyze how our results depend on the strength of convective vigor and fluid inertia. Finally, in Section V the results are extrapolated to the environment of an extremely vigorous, global magma ocean and that of a large magma chamber.

II. GOVERNING EQUATIONS
Rayleigh-Bénard convection of an incompressible isoviscous fluid is governed by the Boussinesq equations: where U(X, τ ) and T (X, τ ) are respectively the velocity and temperature fields, ν the kinematic viscosity of the fluid, ρ 0 the mean mass density at the reference temperature T 0 , α the volumetric thermal expansion coefficient with respect to the reference temperature, g the gravitational acceleration, κ the thermal diffusivity. The hydrostatic stress, ∇P 0 = ρ 0 g, is already subtracted from Eq. (1), leaving only the dynamic pressure P on its right-hand side (RHS). As shown by the last term in Eq. (1), only temperature-induced variations of density are considered to drive the flow. Equations (1) -(3) are solved in a 2D box with periodic side walls and aspect ratio 2. No-slip conditions are assumed on the top and bottom boundaries, which are are isothermal, with a constant temperature difference ∆T driving thermal convection.
We non-dimensionalize the governing equations by scaling the length with the height of the box H, x := X/H, the velocity with the characteristic velocity u * := √ αg∆T H, u := U/u * , the density with the reference density ρ 0 , and we introduce non-dimensional temperature θ := (T − T 0 )/∆T . For the time and pressure it then follows: t := τ u * /H, and p := P/(ρ 0 u * 2 ).
In terms of non-dimensional quantities, the governing equations read: where Ra and P r are the Rayleigh and Prandtl number that control the flow characteristics: The fluid carries inertial particles, whose trajectory is governed by friction from the surrounding fluid in combination with particle buoyancy. Under idealized conditions of spherically-shaped particles with small Reynolds number, the Lagrangian equation of motion for a massive particle reads [e.g. 45]: Settling of particles in thermal convection Submitted to Phys. Rev. Fluids where V is the particle velocity and the first term on the RHS denotes the material derivative of the fluid velocity. The modified density ratio β = 3ρ f /(ρ f + 2ρ p ) relates the density of the fluid ρ f with the particle density ρ p . Eq. (8) is a truncated version of the original equation derived independently by Maxey and Riley [46] and Gatignol [47]. Due to the small Reynolds number and size of the particles we neglect here both the unsteady drag term, known as history term, and the Faxén corrections. We systematically vary the ratio ρ f /ρ p and assume that it is constant for each particle throughout the simulation. This assumption does not necessarily neglect the density variations of material: in view of the Boussinesq approximation ρ f = ρ 0 [1 − α(T − T 0 )], employed in Eq. (1), the assumption simply means that each particle is always at the same temperature as the surrounding fluid and has the same thermal expansivity [for effects resulting from keeping the particles at a different temperature than the fluid, see 48]. The particle response time τ D = r 2 c /(3νβ) depends quadratically on the particle radius r c , which we also vary systematically.
After non-dimensionalizing it with the same scales introduced above, Eq. (8) takes the form: leaving three non-dimensional parameters to control the particle dynamics, namely: The Stokes number St is the particle response time τ D divided by the characteristic time H/u * . It describes the viscous friction acting on each particle due to the difference between particle and fluid velocity. The parameter Λ (hereafter buoyancy ratio) expresses the relative importance of particle buoyancy with respect to the thermally-induced buoyancy of the fluid (the unit vectorẑ points vertically upward). The first term on the RHS of Eq. (9) is the so called added-mass as estimated by [49]. We do not consider any feedback mechanism with respect to the flow: the fluid velocity u is obtained from Eqs. (4)-(6) and does not depend on the particle velocity v, i.e. we adopt a one-way coupling [for the distinction between one-way and two-way coupling, see the review of 50]. In a turbulent flow, the adopted particle model can be considered appropriate as long as the particle size, r c (resp. r c /H in dimensionless units) is up to the same order of magnitude as the spatial dissipative scale of turbulence, η. In Rayleigh-Bénard flow the global value of such scale, in the current dimensionless units, goes as η = P r 1/2 (Ra(N u − 1)) −1/4 [e.g. 51], meaning that it decreases at increasing the thermal forcing Ra (and so the Nusselt number N u) but increases at increasing the Prandtl number P r (see also Discussion).
For particles suspended in a fluid at rest, i.e. with u ≡ 0, Eq. (9) can be solved analytically, yielding: In the limit t→∞, the so-called terminal or Stokes' velocity v t is reached: where the terminal velocity is defined positive for sinking particles and negative for rising particles. We inject 10 6 particles of 301 different types into a fully developed, two-dimensional, statistically-steady thermal convection, with each particle type represented by three values: St, Λ, and β. Since we are primarily interested in the dynamics of the particles, we refer to the thermal flow of the carrier as the "background" flow. For particles with v t > 0 (⇔ Λ < 0 ⇔ β < 1), i.e. those denser than the fluid (labeled as heavy), we measure the time it takes until they reach the bottom boundary. For particles with v t < 0 (⇔ Λ > 0 ⇔ β > 1, labeled as light) we do the same with respect to the top boundary. For brevity, both these cases are referred to as "settling". Initially, all particles are distributed uniformly across the domain and their velocity is set equal to the local velocity of the fluid. 300 different types of particles are obtained by evenly sampling ρ f /ρ p and r 2 c ; one particle type is reserved for fluid tracers. The above described model system is numerically simulated by means of the Eulerian-Lagrangian code ch4-project [52]. The code adopts a lattice Boltzmann (LB) algorithm for the computation of the fluid and temperature dynamics, while it uses a second order time-stepping and grid-to-particle bi-linear interpolation for the computation of particles' trajectories. This code has been already extensively employed in studies involving turbulent thermal convection and inertial particle dynamics [53].

A. Model parameters
We aim to map the settling behaviour of particles over the entire St, Λ, and β space, while focusing on highly vigorous convection (Ra = [10 8 , 10 10 , 10 12 ]), with moderate to small importance of inertia (P r = [10,50]). As such, our results are applicable to a range of natural systems (see Section I). Throughout the paper we strictly use non-dimensional control parameters, but it is instructive to demonstrate how these are linked to physical parameters of a particular system, namely the thermal convection of a large reservoir of crystallizing magma. In this section we inspect how the parameter space is mapped and discuss intrinsic limitations of our numerical approach.
In Table I we list the physical parameters that roughly describe the thermal convection of a primordial, mantle-deep magma ocean for the Earth. A relativley large uncertainty is in the value of the viscosity of high-pressure and -temperature magma. First-principles simulations suggest that the kinematic viscosity of MgSiO 3 , one of the major mantle silicates, over the temperature and pressure range relevant for a global magma ocean (∼ 2000 − 4000 K and 0 − 130 GPa) is on average of the order of 10 −5 − 10 −6 m 2 /s [54]. Since the Prandtl number is defined as ν/κ, the lower and upper bounds of ν define the range of interest of P r and we indicate ν directly as P r × κ in Table I. The temperature contrast ∆T driving convection is also difficult to determine precisely. The reference value of only 1 K reported in the table roughly corresponds to the contrast predicted by parameterized models of the thermal evolution of the Earth's magma ocean in the presence of an atmosphere [e.g., 32,33]. Such a low value is also representative for planetary cores, where a large volume of low-viscosity metallic liquid undergoes thermal convection [e.g., 55]. ρp/ρ f [0, 2] a from Solomatov [8]; b from Ni et al. [56]; c see Karki and Stixrude [54] for typical viscosities of silicate liquids at high pressure and temperature; d see e.g. Lebrun et al. [32] and Nikolaou et al. [33] for typical temperature contrasts during the evolution of magma oceans.
Sampling the ranges of r ref c and ρ p /ρ f from Table I results in a sampling of the non-dimensional parameter space β, St, Λ (set A in Fig. 1). The y-axis in Fig. 1 represents the absolute value |Λ| rather than Λ in order to fit both light (Λ > 0) and heavy (Λ < 0) particles into a compact plot. The modified density ratio β is marked by color only. Later we will show that, apart from the effect described in Section IV B, the first term on the RHS of Eq. (9) has secondary importance on the settling behaviour, which sidelines the relevance of β.
The terminal velocity v t and the response time St can be used to a priori estimate the number of time steps that are required to evaluate the settling time of a given particle type. In the time 1/|v t |, each particle would cross the model domain vertically if sinking (or rising) at the speed v t , making 1/|v t | a proxy for the minimum required duration of a simulation (and thus CPU time). Due to constraints arising from the numerical integration of Eq. (9), which we explain below, we use max(10∆t/St, 1)/|v t | to estimate the minimum required CPU time (green line in Fig. 1). Here, ∆t denotes the maximum time step allowed by the Courant-Friedrichs-Lewy condition, i.e. by the advective and diffusive time scales of the background flow.
For illustration purposes, let us assume u ≡ 0 and discretize Eq. (9) via explicit, first-order Euler scheme: It follows that ∆t/St must be smaller than 2 in order to avoid numerically unstable solutions. Demanding numerical accuracy limits the admissible values of ∆t/St even further -only a small fraction of v t must be added at each time step to ensure that convergence to the stationary solution v t is smooth. In a turbulent flow (u = 0), Eq. (9) yields accurate trajectories only when ∆t/St < 0.1 (the exact value depends on the employed numerical scheme, with 0.1 resulting from our experience with the second-order Adams-Bashforth formula that we use to advect the particles). This constraint increases the CPU time of each simulation (resp. decreases the allowed time step) by an additional factor, 10∆t/St, where St is the smallest Stokes number in the respective set of particles. Red lines in Fig. 1 mark isolines of v t . In Section III, we show that to first-order the settling behaviour of particles can be described by their terminal velocity only: for a given background flow, particles with the same v t settle in a similar manner. Since the green and red lines in Fig. 1 have different slopes, it is convenient to modify the parameters from Table I  For the purpose of this study, particle sets B and C can be simply understood as the selected coverage of model parameter space (later we argue that based on these sets we map the entire St, Λ, β space reasonably well). We note, however, that both sets also have a certain geophysical interpretation. When compared to the original set A, set B has enlarged r c and a narrowed ρ p /ρ f range. As such, its parameters roughly correspond to large clusters of crystals with density similar to that of the surrounding magma. Set C has a reduced thermal expansivity and a larger temperature contrast that lies in the range of temperature contrasts characteristic of magma oceans that cool in absence of an atmosphere [32,33].
Based on the parameter values listed in Table I, the Rayleigh number Ra ref of a mantle-deep magma ocean would be of the order of 10 27 . This is far from being reachable with any numerical method because the thickness of the thermal boundary layer in a convecting system scales approximately as Ra −1/3 , demanding higher resolution for higher Ra. Here we model a series of Rayleigh numbers up to Ra = 10 12 , using up to 4096×2048 grid points. The possibility of extrapolating our results to higher Rayleigh numbers is analyzed in Section V.
We note that our code is based on a dimensional formulation. Therefore, in order to reduce the Rayleigh number we need to modify some of the parameters in Table I For each simulation set, we first wait for thermal convection to develop into a statistically steady state and then we inject all the particles at once, distributing them uniformly in space and assigning them the velocity of the carrier fluid, i.e. v(t=0, x, z) = u(t=0, x, z). In Fig. 2a, we show the average root mean square velocity of the background flow for all the tested values of Ra and P r. Fig. 2b shows the corresponding Reynolds number, Re := U rms H/ν. For P r = 10 we run simulations with Ra equal to 10 8 and 10 10 , while for P r = 50 we test three values, Ra = 10 8 , 10 10 , and 10 12 (lowering Prandtl number increases the resolution demandssee the Reynolds number for two simulations with the same Ra but different P r). Our simulation sets are labeled as B or C, depending on the range of particle parameters (see Fig. 1), and by upper and lower indices we label the exponent of Ra and the value of P r (e.g., C 10 50 stands for simulation set C with Ra = 10 10 and P r = 50).

III. RESULTS: SETTLING CURVES
In Section II A we anticipated that the terminal velocity v t is capable of sorting the settling behaviour of particles in a flow of given Ra and P r. In order to compare the settling behaviour also across flows with a different convective vigor, one more parameter is needed. Similarly to the experimental study of Martin and Nokes [5], who divide v t by the average vertical velocity of the flow and use the resulting ratio to organize their results, we use u rms to account for the properties of the background flow. In this section, we show that there are four distinct regimes describing the sinking or rising of particles, and that the ratio |v t |/u rms determines to which regime a given particle type belongs.
In Fig. 3 we plot the temporal evolution of the settling process for the simulation set C 10 50 , which is taken as a reference case. When heavy particles (v t >0) reach the bottom, resp. light particles (v t <0) reach the top, we eliminate them from the flow and mark them as settled. Each line in Fig. 3 represents a different particle type, although we do not show the respective values of St, Λ, and β. Instead, we mark each line by the value of |v t |/u rms -this single parameter uniquely orders the obtained settling curves. We identify four distinct groups: i) "stone-like" regime, |v t |/u rms > 2.0 (Fig. 3a); ii) bi-linear regime, 0.3 < |v t |/u rms < 2.0 ( Fig. 3b); iii) transitional regime, 0.02 < |v t |/u rms < 0.3 (Fig. 3c); and iv) "dust-like" regime |v t |/u rms < 0.02 (Fig. 3d).
The time on the x-axis of Fig. 3 is multiplied by v t for each settling curve individually. The x-axis thus represents distance rather than time. The "terminal distance", t v t , corresponds to the distance a particle with a given v t would travel in a fluid at rest by the time t (i.e. t v t = 1 represents sinking with the Stokes velocity through the entire container). This means that, even though the particles in Fig. 3d seemingly take only approximately 5 times longer than those in Fig. 3a to completely settle, the actual time differs by more than two orders of magnitude because the corresponding value of v t differs by more than a factor hundred in both subplots. The same applies for different settling curves within each subplot: two settling curves that overlap but have different colors correspond to different settling rates with respect to time t.
FIG. 3. Settling curves obtained from the simulation set C 10 50 . All particle types are separated into four groups based on the |vt|/urms ratio, which is shown in color for each subplot. The x-axis represents the terminal distance t vt. Black and red solid lines represent the analytic solutions (14) and (15) respectively. Fig. 3 is the theoretical prediction

The black line in
where N s and N 0 are the number of settled particles and the initial number of particles respectively. The velocityṽ is given by Eq. (11), and zero initial conditions are considered,ṽ(x, t=0) ≡ 0. Eq. (14) thus expresses the percentage of particles that would settle at the time t if initially they were distributed uniformly in a still fluid. The shape of the black curve is nearly identical to simply min(t v t , 1) because it takes a negligible time for the particles to accelerate from 0 to v t (see Section III A for a further discussion). Below we analyze the settling regimes individually and explain underlying mechanisms. In Section IV we discuss how the regimes' properties and boundaries depend on the characteristics of the background flow.
A. "Stone-like" regime (|vt|/urms 2) The simplest regime corresponds to the case with a high |v t |/u rms ratio. For a particle with |v t | > 2u rms , the average convective velocities are more than twice smaller than the speed at which the particle would be sinking if there was no convection. This implies that particles with this property are little affected by the flow -they sink almost as if the fluid was at rest because thermal convection is slow relative to the particle's vertical drift.
Therefore, when |v t |/u rms 2, Eq. (11) provides a good prediction of the settling behaviour ( Fig. 3a). As analyzed later, this result is very robust with respect to the values of Ra and P r because the background flow is nearly irrelevant in this regime.
For even higher |v t |/u rms ratios, the fit to Eq. (14) becomes perfect, and the acceleration from 0 to v t begins to play a role in the shape of obtained solutions. In Fig. 4 we demonstrate this effect on the simulation set xC 10 50 , constructed using particles with 10 times larger radii than in the reference set C 10 50 . For clarity of the figure, we only show a few particle types from the set, with the |v t |/u rms ratio ranging from 2 to 50. As |v t | and St increase in value, the analytic solution (14) loses its symmetry with respect to light and heavy particles because its second term gains in relative importance. The second term in Eq. (14) is not symmetrical with respect to v t : particles with the same terminal velocity but different modified density ratio β have different values of St (recall the definitions of St, Λ, and that v t := −StΛ). As a result, light particles (β>1) have a shorter response time St when compared to heavy ones (β<1), and accelerate to v t faster (cf. the last term in Eq. 11).
We label this regime stone-like. Although particles with |v t |/u rms 2 still interact with the structure of the flow (see Section IV A below), their vertical speed is close to the free-fall speed |v t |. Moving to lower |v t |/u rms ratios, the settling curves become approximately piecewise linear, with two distinct settling rates. The two different rates correspond to different initial positions of the particles.
In Fig. 5 we depict particles with 0.3 < |v t |/u rms < 2.0 that are still suspended at the time t = 0.4/(0.3 u rms ) ≈ 13. This snapshot corresponds to the time at which the settling curves of particles with |v t |/u rms = 0.3 change their slope (cf. the darkest settling curve in Fig. 3b). The particles in Fig. 5 form a cloud, centered above a cluster of upwellings, and the larger the value of v t , the smaller is the respective cloud. The surroundings of major downwellings are free of particles. Note that we depict heavy particles only; light particles are located above the central cluster of downwellings.
After t ≈ 13, all particles with 0.3 < |v t |/u rms < 2.0 settle at a reduced rate, while up to this time their settling is well captured by Stokes' law ( Fig. 3b). It follows that if a heavy particle (with 0.3 < |v t |/u rms < 2.0) is initially injected close to a major downwelling, or below the top boundary layer where horizontal velocities are large, it settles quickly. Perhaps surprisingly, the settling rate of such particles is represented well by Stokes' law and does not exceed it (compare the first linear segment of all settling curves in Fig. 3b). One could expect the downwellings to sediment the carried particles downstream, speeding up their settling beyond the rate predicted by Stokes' law. This, however, does not happen. In the first stage of the bi-linear regime, particles touch both horizontal boundaries with little to no lateral preference (i.e. the x-coordinates FIG. 4. Selected settling curves from the simulation set xC 10 50 that consists of particles with a high |vt|/urms ratio (2 to 50 for the depicted selection). Solid lines show the theoretical prediction (14), while colored points are obtained from the numerical simulation.
of the settling events have a uniform distribution), and the percentage of settled particles grows linearly in time, with a slope that matches the Stokes velocity. We pay further attention to the horizontal distribution of settled particles in Section IV A.
In the second stage of the bi-linear regime, the settling rates are significantly reduced. This is because bursts of upwelling flow act against the particles' tendency to settle. The terminal velocities studied in this section are still sufficiently large (|v t | > 0.3 u rms ) for the particles to efficiently penetrate through the fluid flow, but at the same time the existence of plumes alters particle trajectories significantly, in particular by lifting particles that get caught in strong conduits. The higher the |v t |/u rms ratio, the less particles survive after the first-stage settling, and the closer they are to the central axis of the major upwelling structure (compare black and white dots in Fig. 5). As a result, the settling rate is smaller for higher |v t |/u rms ratios (Fig. 3b). Note that throughout the text we refer to the relative settling rate, i.e. to the slopes of the settling curves in Fig. 3, where the x-axis represents distance rather than time. With respect to time t, the settling is generally faster for higher |v t |/u rms ratios.
For |v t |/u rms ≈ 2, the settling curves are close to being flat in the second stage of the bi-linear regime. This is in agreement with expectations: the velocities of plume heads are typically close to 2u rms and the carried particles thus could, in principle, be indefinitely suspended in a fixed point in space by the action of a stationary plume (cf. Eq. (9) with u ≡ −v t ). For Ra = 10 10 the flow is highly non-stationary and such situation never occurs, but the idealized scenario helps explaining the very slow settling rates. The settling curves are nearly identical when |v t |/u rms 0.02. The very existence of a limit is a nontrivial result. While one can a priori expect the applicability of Stokes' law when |v t |/u rms → ∞, when |v t |/u rms → 0 the particles should behave as fluid tracers. This can be seen directly by inspecting Eq. (9): for a fixed Λ, if the terminal velocity tends to zero, the Stokes number St also tends to zero. In the limit St → 0, the second term on the RHS dominates Eq. (9), yielding v = u (i.e. fluid tracers). For this reason, τ D is called the response time -particles quickly adapt to the carrier fluid velocity when τ D is small. Fluid tracers, however, never touch the bottom nor the top boundary of the model domain and it is thus a priori unclear how particles with a small |v t |/u rms ratio should settle. A simple theoretical model for small particles (i.e. with a small Stokes number) was developed by Martin and Nokes [5]. They proposed that at the base of the model domain, where convective velocities vanish, all particles are free to settle from the fluid with a speed equal to their terminal velocity v t . Therefore, the rate at which the number of particles in the flow N decreases with time is given by where A is the area of the base of the domain, c(z) the horizontally averaged concentration of particles at height z above the bottom boundary, and N 0 the initial number of particles. The exponential solution in Eq. (15) is obtained by assuming c(0) to be the current average concentration N/(AH).
Indeed, Fig. 3d shows that 1 − exp(−v t t) fits the settling curves well for |v t |/u rms → 0, confirming the theoretical and experimental conclusions of Martin and Nokes [5]. The derivation of Eq. (15), however, is based on a counter-intuitive assumption: particles must be uniformly distributed throughout the entire model domain by convection, and yet there must be a boundary layer with little to no mixing, thick enough for the particles to separate from the fluid and accelerate to v t . Moreover, the concentration of particles in the boundary layer is assumed to be the same as in the bulk of convection, without assessing the mutual transport between the two regions.
In the next section, we focus on the statistics of particle transport between convection cells and boundary layers of the flow. We interpret particle settling as a random process, allowing us to provide a quantitative description of the settling regimes' boundaries, and to explain in detail why the exponential law (15) fails for particles with |v t |/u rms 0.02. Most importantly, our approach serves as a unifying theory, capable of containing all four settling regimes in a new equation that estimates the time required for a complete sedimentation of the particles.

D. Particle settling as a random process
The velocity structure of the simulation C 10 50 is shown in Fig. 6. On top of the temperature field, shades of green show regions with |u| < 0.5 u rms . Below, these regions are referred to as the "low-velocity piles" (the definition may seem rather arbitrary here -later we investigate piles defined generally as the regions where |u|/u rms < p f and vary the pile factor p f ). The centres of the large scale convection rolls, which also show small velocities, are not considered as the low-velocity piles (thick green line in Fig. 6, explained later). Grey points in Fig. 6 show suspended particles, this time we select only those with |v t |/u rms < 0.3. For dust-like terminal velocities (|v t |/u rms < 0.02), the particles appear uniformly distributed, while for |v t |/u rms → 0.3 the spatial distribution resembles the one from Fig. 5. In this section we explore the rate at which particles enter the low-velocity piles, and analyze how likely it is that a particle enters one but does not settle at its base and returns to fast convection cells instead. In Fig. 7 we take a particular instant in time and compare Eq. (15) with the number of particles that have actually settled (red line vs green dots). In agreement with Fig. 3, there is a good match with the exponential law for |v t |/u rms 0.02, but for 0.02 |v t |/u rms 1.5 the settling is slower than Eq. (15) predicts.
Similarly as in Fig. 3, the black line represents |v t |t, i.e. the prediction based on simple Stokes' settling. Note, however, that while Fig. 3 maps the temporal evolution of the settling, Fig. 7 captures only a snapshot in time and does not contain any information about the settling rate. For instance, particles with |v t |/u rms 1.5 have completely settled at time t = 50 -both the exponential and Stokes' laws match the observation and Fig. 7 cannot be used to distinguish between the two, although Fig. 3a shows that particles with a high |v t |/u rms ratio follow Stokes' law.
Blue stars in Fig. 7 mark the percentage of particles that have entered the low-velocity piles at least once by the time t = 50. For heavy particles we only consider the piles located in the bottom fourth of the model domain (thick green line in Fig. 6), for light particles only the upper fourth of the domain is considered. The center of the graph (v t = 0) represents fluid tracers (treated as heavy for the purpose of this analysis).
It is important to notice two things: First, the probability of ever entering the low-velocity piles decreases with the |v t |/u rms ratio, no matter how small the ratio is. This means that, even in the transitional and dust-like regimes, the trajectories of inertial particles differ from those of fluid tracers: particles governed by the Maxey-Riley equation are more likely to cross sluggish regions of the flow. Second, the probability of escaping the low-velocity piles and returning to the flow increases with the |v t |/u rms ratio. This can be seen from the difference between the percentage of particles that have entered the piles and the percentage of particles that have settled: while for higher values of |v t |/u rms the two nearly coincide, for lower values they differ, with the difference growing as |v t |/u rms approaches zero (compare blue and green symbols in Fig. 7).
In Fig. 8 we plot the probability of escaping the low-velocity piles. It is computed as follows: when a particle reaches a region with |u| < 0.5 u rms , we mark it as captured. If a captured particle, instead of settling at the wall of the container, is transported back to a region with |u| > u rms , we mark it as escaped. Each particle can be captured and escape multiple times and we store the record for each individual particle. The escape probability P e is simply the sum of all escapes divided by the sum of all captures, taken over all particles of a given type. As expected, P e goes to one for fluid tracers, because tracers never settle and eventually always escape any low-velocity regions. On the other hand, P e tends to zero for particles with a large terminal velocity, because such particles always settle upon encountering a slow region. Now we have the means to describe particle settling as a random process. The exponential law (15) and its underlying differential equation dN = −|v t |N dt are, in a way, calling for such approach: in the terminology of stochastic processes, the equation says that particle settling is a Poisson counting process with intensity |v t |.
A random process is defined through an event of a given probability. Here, the event is the settling to the base of a low-velocity pile, with the probability equal to 1−P e . For each particle type, the number of particles that settle over a time dt can be written as where dN c is the number of particles that enter the low-velocity piles in the time interval t, t+dt . In other words, dN c describes the supply of particles into our random process, and, as outlined by Fig. 7, it is a function of |v t |. Under the assumption (indicated by the question mark above the equality sign) that particle settling is a Poisson counting process with intensity |v t |, we may further write Eq. (17) can be validated numerically. However, due to intrinsic fluctuations of thermal convection, it is convenient to integrate Eq. (17) in time and compare the integral quantities instead, The RHS of Eq. (18) is plotted in black in Fig. 8a. The ratio of the two integrals is labeled as t c because it is the average time between particle captures, that is, between repeated entries into the low-velocity piles. By comparing 1−P e and |v t | t c , we observe that for particles with sufficiently small terminal velocities (|v t |/u rms 0.02), the governing equation of the Poisson process is well satisfied (compare blue and black symbols in the inset of Fig. 8a). This confirms the idea underlying the derivation in Eq. (15). Note that for a small value of p f our description reduces to the one by Martin and Nokes [5]: piles with a vanishing p f factor become only thin layers at the top and bottom boundaries, from which there is no escape (i.e. P e =0).
When P e = 0, Eq. (17) reduces to N (t) |v t | dt ? = dN c , where the left hand side expresses the particle flux that can be expected through any horizontal plane if uniformly distributed particles drift vertically with the Stokes velocity. In the limit of p f → 0, the right hand side is the supply of particles into the boundary layer. According to the assumptions of Martin and Nokes [5], these two particle fluxes are equal (Eq. (15)).
Our results show that this assumption is not valid generally: the value of |v t | t c reaches up to 4 for light particles in the set C 10 50 (black circles in Fig. 8a). Indeed, unlike the probability 1−P e , the value of |v t |t c is not limited by 1: its value merely evaluates how frequent the transport of particles between fast and slow regions of the flow is. The difference between the ratio |v t |t c /(1−P e ) and 1 is then the observed deviation from a Poisson counting process (Fig. 8b).
Already for terminal velocities exceeding ∼ 0.02 u rms , the RHS of Eq. (18) is larger than the probability 1−P e . This explains why the exponential law (15) fails: it is too difficult for particles to enter the lowvelocity piles. In other words, the transport of particles into the sluggish regions, where separation from the fluid flow takes place, is much slower than what one would predict when simply assuming that the particles drift vertically with the speed |v t |.
Comparing |v t | t c and 1−P e provides a unified description of all four regimes. Dust-like: Both quantities are equal and N = N 0 exp(−|v t | t) describes particle settling well. Transitional: The supply of particles into the piles is too small (|v t | t c > 1−P e ), which implies that settling is slower than that predicted by Eq. (15), and the disagreement increases as |v t |/u rms increases. Bi-linear: The supply of particles is still too small, but increases quickly as |v t |/u rms further rises, effectively reducing the difference between |v t | t c and 1−P e . This is because the particles are increasingly efficient in penetrating the fluid flow. Stone-like: Eventually, |v t | t c becomes smaller than 1−P e , and the settling curves become faster than exp(−|v t | t), soon reaching the Stokes's law instead. These newly defined regime boundaries do not exactly overlap those used in Fig. 3, but both definitions are in a rough agreement.
While P e and t c depend on the value of p f that is used in the definition of the low-velocity piles (i.e. |u|/u rms < p f ), their ratio does not (Fig. 8b). The quantity |v t | t c /(1 − P e ) can thus be used to construct a general model, extending the exponential law (15). First, we fit the ratio with a skew normal distribution f : where x = |v t |/u rms and erf is the error function. The amplitude A is prescribed as in order to get the observed match between |v t |t c and (1−P e ) in the limit of zero terminal velocity (i.e. f = 1 for |v t |/u rms → 0). The remaining parameters, µ, σ, and λ, are obtained by fitting the data (see Fig. 8b). The analytic prescription (19) allows us to continue in describing particle settling as a random process. Setting dN c = N dt/t c in Eq. (16), we have Eq. (21) has exponential solution: Equation (22) establishes an extension to the solution (15). It is valid also for |v t |/u rms > 0.02, and bridges the gap between existing analytic solutions for small and large terminal velocities (i.e. the |v t | → 0 and |v t | → ∞ limits). Note that the domain depth H must be added to the denominator when |v t |t are to be replaced by their dimensional counterparts |v t u * |τ . Apart from Eq. (15), Martin and Nokes [5] also develop a more sophisticated theory, in which c(0) is not simply the average concentration N/AH. By assuming a depth-dependent concentration c = c(z) whose temporal changes are governed by the diffusion equation [e.g. 34], Martin and Nokes [5] find solutions for c(0)/c for several flow and particle parameters, withc being the average concentration N/AH (see their Table 2). Note that inserting c(0) =c into Eq. (15) is exactly analogous to our Eq. (22), with c(0)/c being analogous to our 1/f . While Martin and Nokes [5] predict c(0)/c > 1 for particles with a non-vanishing Stokes velocity, we obtain the exact opposite, f > 1. However, the experimental measurements of Martin and Nokes [5] are in agreement with our results, as they systematically measure the settling rates to be slower than Eq. (15) predicts, especially for particles with v t /u rms ≈ 0.5. Martin and Nokes [5] acknowledge the discrepancy between their theoretical prediction and measurements, and speculate that it may be related to the breakdown of the assumption of one-dimensionality. In particular, they anticipate that the large scale circulation in the fluid could be responsible for the failure of the one-dimensional turbulent diffusion theory, which is exactly what we observe.
The misfit between Eq. (22) and the observed settling curves never exceeds 30%, with the largest error occurring for particles with 0.3 < v t /u rms < 1.0. This is not surprising -already from Fig. 3b it is clear that the settling curves are not exponential when v t /u rms 0.3. Nevertheless, Eq. (21) is still useful for estimating the characteristic time of complete sedimentation for all particle types (see Fig. 10 at the end of this section and the accompanying discussion).
The imperfect fit of the observed settling curves is caused by the fact that t c is not a function of |v t | only, but it is also a function of time, t c = t c (|v t |, t). For instance, in the stone-like regime, the number of particles N follows N (t) ≈ N 0 × (1 − v t t), and dN c ≈ N 0 v t dt. Eq. (18) then yields t c ≈ max(1/v t −t/2, 0.5/v t ), with 0.5/v t being the end value that is reached after all particles have settled (see the black line in Fig. 9, resp. the dotted line in Fig. 8b). Eq. (22), on the other hand, assumes that the average time between particle captures, t c , does not vary in time (see the definition of t c in Eq. (18) versus its use in Eq. (21)).
In Fig. 9 we show the temporal evolution of t c for selected v t /u rms ratios. Regardless of the value of v t /u rms , t c equilibrates before t ≈ 100 (note that the x-axis is logarithmic -for most of the depicted timewindow t c is steady). By t ≈ 100, however, most of the particles in the stone-like and bi-linear regimes have already settled (dashed lines). Refining Eq. (22) would thus require accounting for the time dependence of t c .
With the black line in Fig. 9 we plot max(1/0.3−t/2, 0.5/0.3), which is the theoretical value of t c obtained for Stokes's settling of particles with v t = 0.3 (see above). The initial, short-lived drop of the actual evolution of t c (t) (dark blue line) is caused by accelerating to v t from zero (the average vertical velocity of the fluid, u z dV , and thus also the initial average velocity of the particles, is zero). In the transitional and dust-like regimes, on the other hand, t c (t) has only a mild temporal evolution and the transient phase is also less relevant, because the particles take longer to settle. Most of the particles with small terminal velocities settle when t c (t) is completely steady, allowing Eq. (22) to provide a good fit to the observed settling curves (< 10 % error for v t /u rms < 0.1). In Fig. 10 the above results are summarized. In the literature, the terminal velocity v t is used as a measure of the settling rate of particles in a convective flow. This implies that the time required for a complete sedimentation of all particles is 1/(v t ), or H/(v t u * ) in terms of dimensional quantities. Fig. 10 shows the factor F by which 1/(v t ) has to be multiplied in order to obtain the correct settling time. Since Eq. (22) is an exponential law, the predicted time of complete sedimentation is, in principle, infinite. To circumvent such inconvenience, we compute the time until 95% of particles have settled and divide the resulting value by 0.95/v t (i.e. normalize by the respective terminal time). The plotted factor, F , is therefore equal to: The settling time / terminal time ratio, F , is computed with the use of Eq. (22) when v t /u rms < 2 and from Eq. (14) when v t /u rms > 2. We restrict the applicability of our analytic model, because for v t /u rms > 2 the Stokes' formula (11) is more accurate and physically appropriate. The parameters µ, σ, and λ used in Fig. 10 come from the set C 10 50 (see Table II for the respective values, Fig. 10 is plotted for heavy particles only).
Solid and dashed red lines in Fig. 10 are isolines of v t that correspond to the regime boundaries. Since u rms = 0.1, the respective values are v t = 0.002, 0.04, and 0.2. For a more vigorous flow these values must be adjusted accordingly (see Section V).
On top of the analytic prediction (22) we plot the settling times observed in the simulation sets B 10 50 , C 10 50 , and xC 10 50 (pentagons, squares, and circles respectively). Only the particle types for which the simulations have reached at least 95% settling are plotted. The difference between Eq. (22) and the numerical simulations is typically less than 3%. One exception is the vicinity of |v t |/u rms ≈ 0.6, where the settling rates are very small during the second-stage of the bi-linear regime, which significantly prolongs the settling time with respect to expectations (up to 40% discrepancy in the value of F ). Nevertheless, N s (t 95 )/N 0 − 0.95 stays below 7% for all particle types, ensuring reasonable accuracy of Eq. (22) for general applications. Here, t 95 is the settling time computed by setting (N 0 − N )/N 0 = 0.95 in Eq. (22).
The region labeled "slow belt" is characterized by increased settling times, with the average settling rates being up to 13 times slower than Stokes's law predicts in the present conditions Ra = 10 10 and P r = 50. With the exception of the above mentioned discrepancy between the theoretical prediction and measurements of Martin and Nokes [5], its existence is not reported in previous literature, which calls for experimental and 3D numerical confirmation of our findings.

IV. RESULTS: IMPORTANCE OF THE BACKGROUND FLOW
In the dust-like and stone-like regimes, the settling curves are robust with respect to properties of the background flow. In the transitional and bi-linear regimes, on the other hand, the above results suggest that settling curves depend on large-scale circulation of the fluid (Fig. 5). Thermal convection of an isoviscous fluid naturally leads to the formation of low-velocity piles similar to those depicted in Fig. 6 (see also Discussion), but their coherence and erosion depend on the values of Ra and P r. In this section we analyze the interplay between particle settling and the velocity structure of thermal convection.

A. Horizontal distribution of settled particles
The escape probability P e is close to zero already for |v t |/u rms 0.3 (Fig. 8), and the settling problem is thus mostly reduced to measuring the transport of particles from the bulk of convection into the low-velocity piles. Here we show that the near-boundary regions depicted in green in Fig. 6, i.e. the piles with p f = 0.5, act as dominant sinks for most particle types.
In Section III B we demonstrated that in the second stage of the bi-linear regime, heavy particles hover inside a cluster of upwellings (Fig. 5). Eventually, these particles settle in the underlying low-velocity pile, whose edges are a continuous source of the plumes that keep lifting the particles. Note that settling below upwellings is somewhat counter-intuitive -one could naturally expect heavy particles to concentrate below major downwellings.   11 shows the horizontal distribution of settled particles. The distribution combines information from both stages of settling, but we note that the non-uniformity is produced in the second stage only. In the first stage of the bi-linear regime, the relatively fast sedimentation along with the horizontal drag associated with large-scale convection rolls ensure that particles are distributed evenly across the bottom and top boundaries (see also Video S1 in Supplementary material).
For brevity, and because the horizontal distribution of light particles is analogous to the distribution of heavy particles, we will discuss only the heavy particles in this section. Top panel of Fig 11 shows the vertically-and time-averaged vertical component of the velocity field. Peaks of the function correspond to the edges of large-scale convection cells, whose centres are indicated by the grey symbols in panels b) and c). 12. As Fig. 11, only here we show the horizontal distribution of settled particles in the simulation set xC 10 50 . Fig. 12 shows that also stone-like particles settle preferentially in the large low-velocity pile located at the edges of the model domain (as long as the v t /u rms ratio is 10 -see the dark blue lines). Due to their large sinking velocity, stone-like particles efficiently cut through the bursts of upwellings. Yet, the large low-velocity pile still acts as a sink. This is because, due to the large-scale rolls, the fluid velocities in the lower-half of the model domain generally point towards the side edges, which sways the sinking particles into that direction.
For the particle trajectories to be close to "ballistic", i.e. to follow the analytic solution (11) exactly, the |v t |/u rms ratios of much larger than 2 are necessary. Even when Eq. (11) effectively governs particle dynamics, the horizontal distribution of settled particles is not uniform. The fluid flow then enters Eq. (11) through the initial velocities,ṽ(x, t=0) = u(x, t=0). Because of large-scale circulation, ballistic particles in the upper part of model domain are generally injected in the direction of major downwellings and have sufficient time to move laterally. Particles injected into the lower half of the domain, on the other hand, do not have the time to move to the sides. As a result, the horizontal distribution reverses for |v t |/u rms 30, with most of the heavy particles settling below the major downwelling (light green and yellow curves in Fig. 12).
For particles with 0.02 < |v t |/u rms < 0.3, we plot the horizontal distribution of settling in Fig. 13. In the transitional regime, particles still "see" the low velocity piles, but the settling events become horizontally uniform as |v t | → 0.02. The particles experience on average over twenty flow overturns between the repeated crossings of the low-velocity piles (see Fig. 9). Typically, they circulate in large convection cells, waiting to enter the slow regions through small-scale irregularities of the flow that are produced by births of new plumes (see also the concentration of grey and black dots inside the green regions in Fig. 6). As discussed in Section IV C, for Ra = 10 10 the low-velocity piles are particularly coherent and do not move horizontally (see also the top panels of the figures in this section), which makes it difficult to enter them. Under such conditions, the particle supply into the piles can be particularly small and the horizontal distribution of settling events can be highly non-uniform. FIG. 13. As Fig. 12, only here we analyze the transitional regime, 0.02 < |vt|/urms < 0.03.
In the dust-like regime, the horizontal distribution of settled particles is uniform. These, almost tracer-like particles remain suspended in the fluid for very long times, repeatedly entering and leaving the low-velocity piles (Fig. 7). The escape probability P e approaches 1 as |v t |/u rms → 0, which means that for particles with small terminal velocities the piles become transparent. Convective motions inside the piles, though relatively slow, are still fast enough to drag along the particles with a vanishing response time St (i.e. the dust-like particles). The only structure where such particles can separate from the fluid becomes the thin, laterally uniform part of the no-slip boundary layer.

B. Motion of light particles towards vortices (β-effect)
Describing the settling behaviour with the help of |v t |/u rms only is a crucial reduction of the fivedimensional model parameter space. Given the complexity of the problem at hand, such description can only be used as a first-order approximation.
One difficulty appears already in Fig. 8a, since there is clearly an asymmetry between light and heavy particles that have the same amplitude of the terminal velocity. For Ra = 10 10 , the differences are small and may partially result from a random asymmetry of the up-and downwelling regions. For Ra = 10 12 , the asymmetry becomes a prominent feature (Fig. 14). Its underlying mechanism is explained below.
Due to the term βDu/Dt in Eq. (9), heavy particles (β<1) have a tendency to move away from strong flow vortices [57]. Indeed, in Fig. 15, where we show a snapshot from the simulation C 12 50 , there is a reduced concentration of heavy particles at the edges of the particle cloud, i.e. close to the centres of the two largest convection rolls (in Fig. 5 this effect is not observed because the Reynolds number is too small).
Light particles, on the other hand, move toward flow vortices thanks to the βDu/Dt term [16]. In Fig. 16, we show only the particles with v t /u rms ≈ −0.3 (the light/heavy asymmetry seems largest for this value, see Fig. 14). Around strong and long-lived vortices, there is an increased concentration of light particles. These particles are trapped until the respective vortices vanish, which explains the enhanced flattening of the settling curves. Due to the way particle sets B and C are constructed, the range of β is approximately 0.99, 1.01 for one and 0.6, 3.0 for the other. Since the outward (resp. inward) motion of heavy (resp. light) particles depends on how much β departs from unity, the effects encountered in Figs 15 and 16 for the set C 12 50 do not occur in B 12 50 . In terms of the settling curves, B 12 50 shows a symmetry between light and heavy particles, while C 12 50 exhibits a slightly enhanced settling of heavy particles and a significantly delayed settling of light particles (Fig. 14). Each simulation set contains several particle types with roughly the same |v t |/u rms ratio. For C 12 50 and v t /u rms = −0.3 ± 0.02, the modified density ratio β ranges from 1.4 to 2.5 (see Fig. 1). The respective particle types are shown in Fig. 16. For these few particle types, the concentration in the vicinity of stable vortices seems similar, and also their settling curves are comparable. Based on this particular example, we can crudely conclude that for β ∈ 1.4, 2.5 the β-effect is similarly strong and for β < 1.01 no effect is observed. A more detailed analysis should be the subject of future work.
One way to quantify the β-effect as a function of v t /u rms is through the v t t c /(1−P e ) ratio. Already in Fig. 8 there was a difference between light and heavy particles, and the gap further increases as the Reynolds number increases. In Fig. 17 we plot v t t c /(1−P e ) for the sets B 10 10 , C 10 10 , and C 12 50 . Consistently with the . Light particles (circles) experience slower settling than the heavy ones, as long as the β-effect comes into play (sets labeled as C). The separation between the ratios for light and heavy particles increases with the Rayleigh number (cf. the orange and green symbols).
analysis above, light particles settle very slowly in C 12 50 , with v t t c /(1−P e ) going up to 14 for particle types that are trapped inside vortices for a particularly long time, while there is little difference between the settling rates of light and heavy particles in the set B 10 10 . Note that for heavy particles the value of v t t c /(1−P e ) may drop below 1 due to the β-effect (green triangles in Fig. 17). This result can be interesting in view of the debate between Wang and Maxey [18] and Mei [58]: the former observed faster than Stokes' settling of heavy particles due to preferential sweeping in downward moving fluid, while the latter did not observe the effect [see also 19,59]. Here we find an increased settling rate that is higher than predicted by the exponential law of Martin and Nokes [5], but it never exceeds the Stokes' velocity on average. Note, however, that our P r is always relatively high, limiting turbulence effects.

C. Effects of convective vigor and fluid inertia
In this section we analyze how the above results are affected by Ra and P r of the background flow. In Table IIa we show the parameters µ, σ, and λ of the skew normal distribution f for all the simulation sets labeled as B, along with the maximum value of f , denoted as m, and the value of |v t |/u rms at which the maximum is reached. The maximum m can be used as an estimate for how much slower the particle settling can be when compared to the exponential decay N = N 0 exp(−|v t |t).  The most important outcome of the comparison in Table IIa is that the critical ratios |v t |/u rms that mark the regime boundaries are largely independent of Ra and P r: the maximum of the function f , i.e. the boundary between the transitional and bi-linear regimes, always lies at |v t |/u rms ≈ 0.4 ± 0.1. With a similar accuracy, the stone-like regime is always obtained for |v t |/u rms ≈ 2.0 (recall that Stokes' settling satisfies v t t c /(1−P e ) = 0.5, see the dotted line in Fig. 8b).
The settling behaviour in the limits |v t |/u rms → ∞ and |v t |/u rms → 0 can be derived analytically regardless of the values of Ra and P r: in the first case, the flow is irrelevant. In the second case, the only requirements are those discussed above in relation to Eq. (15). The critical values of |v t |/u rms for which the settling curves start to substantially deviate from these limits, could, however, strongly depend on the Rayleigh and Prandtl numbers. The fact that there seems to be no such dependence makes the possibility to extrapolate our results to arbitrary thermal flows promising.
The only notable difference between the various simulation sets is the width of the bi-linear and transitional bands, i.e. the spread of the settling curves in the transitional and bi-linear regimes. The width of the band is directly linked with the maximum value m (compare e.g. the sets B 12 50 and C 12 50 in Figs 14 and 17), and m depends non-trivially on Ra and P r. Generally, there is a trend between m and the Reynolds number, with m being the largest for Re ∈ [10 3 , 10 4 ].
Increasing the Reynolds number Re enhances the short-wavelength content of the velocity field, and alters the stability of large-wavelength structures. In particular, the low-velocity regions become less stable (see Fig. 15). High convective vigor is capable of tearing the sluggish, boundary-based piles into chunks that are advected into the rest of the fluid and mixed. This results in spatial variations of large-scale circulation and thus in faster settling because the sinking particles spread across a broader area (see Videos S2 and S4 in Supplementary material). Surprisingly, particles with 0.02 < |v t |/u rms < 2.0 show faster settling also when Ra, resp. Re decreases. As expected, in the simulation C 8 50 the plumes are thicker and live longer than those in the simulations C 10 50 and C 12 50 . The thicker and well separated plumes, however, allow particles to sink in between them, which results in an increased settling rate in the second stage of the bi-linear regime when compared to the reference set C 10 50 (see also Video S5 in Supplementary material). For sets labeled as C and Re ≥ 10 2 , the β-effect splits the v t t c /(1−P e ) ratio into two clearly distinct functions. In Table IIb we provide the corresponding sets of parameters µ, σ, and λ. The difference between m for the heavy and light particles increases with the Reynolds number.

V. APPLICATION TO CRYSTALLIZING MAGMA
The Rayleigh number of a magma chamber is of the order of 10 9 -10 17 [e.g. 60], making our results directly applicable to relatively small volumes of magma. For a global magma ocean, however, Ra ∼ 10 27 (Table I).
First of all, in order to apply our results to this system, it is necessary to estimate u rms for such an extreme flow regime.
In terms of non-dimensional control parameters, the v t /u rms ratio can be expressed as: where the Reynolds number is defined through the volume-averaged root mean square velocity, Re := U rms H/ν. Upon employing a Re = Re(Ra, P r) scaling, Eq. (24) can be used to compute the |v t |/u rms ratio for various Rayleigh and Prandtl numbers. Relationships between Reynolds number and Rayleigh and Prandtl numbers have been obtained based on various experimental, numerical, and theoretical work. Here we adopt the Grossmann-Lohse theory [61]. The theory defines four regimes for isoviscous thermal convection, depending on whether kinetic and thermal energy dissipation takes place dominantly in the boundary layer region, or in the convective bulk.
For the ranges of Ra and P r investigated here, the energy dissipation is dominated by the convective bulk, and the thermal boundary layer is nested inside the kinetic one. In an idealized case such situation yields Re ∝ Ra 4/9 P r −2/3 (see Table 2, regime IV u in Ahlers et al. [37]). Based on our 5 data points, we observe Re = (0.07±0.01) (Ra 4/9 P r −2/3 ) (1.31±0.02) , i.e. our exponents differ by ca. 30% from the idealized case. However, the Grossmann-Lohse theory is derived for a 3D box with thermally insulating side-walls, while we perform 2D simulations with periodic sides, which may be a source of the discrepancy (see also Discussion). When extrapolating to Ra = 10 27 according to our Re(Ra, P r) relationship, Eq. (24) gives u rms ≈ [5.2, 2.9] for P r = [10,50].
In Fig. 18a we show the settling time of crystals in a global magma ocean (Table I). The results are obtained by multiplying the terminal time t t := 1/|v t | (resp. H/(v t u * ) in dimensional units) with the factor F = F (|v t |/u rms ). The black triangles in Fig. 18a correspond to heavy particles in the particle set A, i.e. to the black triangles from Fig. 10. The mean velocity of the flow is assumed to be u rms = 4.0, a value representative of the relevant range [2.9, 5.2]. Since u rms is now larger than in Fig. 10, the slow belt moves to the right with respect to the positions of the particle types, which now fall into the dust-like and transitional ranges (see the red lines and black symbols in Figs 10 and 18a).
The factor F is computed using the a i coefficients derived for the particle set C 12 50 , as this simulation has the highest Ra and includes the β-effect. F reaches a maximum of 19 for light particles, while for the heavy ones it only slightly exceeds 3 (see the light/heavy asymmetry in Fig. 18a).
Note that there are three independent parameters in the simplified Maxey-Riley equation (9), but in the dimensional version (8) there are only two, τ D and β, because the gravitational acceleration g is usually fixed. Considering St Λ, and β as independent is thus a generalization for arbitrary gravity. In Fig. 18 gravity is fixed again, g = 9.8 m/s 2 , and the two independent variables are chosen as ρ p /ρ f and r c , because these quantities are typically measured.
The settling times that we obtained tend to be small compared to the typical lifetimes of magma oceans [e.g. 32,33]. This suggests that a magma ocean will likely solidify via fractional crystallization. An exception is for particle types with extremely small density contrasts. For example, for r c = 1 mm, the value of |ρ p /ρ f −1| must be smaller than 4 × 10 −6 in order to obtain settling time longer than 1 Myr, while the typical values are ca.   Table I, for which we estimate urms to be [2.9, 5.2]. The black triangles correspond to the black triangles in Fig. 10, circles are the respective light particle types. We show log 10 of the value of ttF , where tt is the terminal time H/(vtu * ) and F accounts for the characteristic settling behaviour of the various particle types (see Eq. (23) and Table II, the parameters µ, σ, and λ correspond to the set C 12 50 , i.e. to the simulation with the highest Ra). (b) Same as above, only this time urms = 0.5 instead of 4.0 and the size of the system is H = 2.9 km instead of H = 2890 km. The particle set now spreads across the slow belt. We show only the heavy particles from the set (the position of particle types and the red lines are symmetrical with respect to 1 on the x-axis).
While the position of the slow belt in the (St, Λ)-diagram depends on the value of u rms only, the positions of crystals of given sizes and radii in that diagram depend on various other parameters (see Eq. (10)). Similarly as for a global magma ocean, we can provide first-order estimates also for magma chambers. The key difference between the two is in the value of H. When H = 3 km instead of H ≈ 3 × 10 3 km, crystals with the same radius range r ref c = [0.5, 10] mm move by 3 orders of magnitude to the right in Fig. 10. At the same time, the lower value of H reduces the Rayleigh number to ≈ 10 18 , yielding u rms ≈ 0.5, which shifts the slow belt by a factor of 5 to the right in Fig. 10. As a result, the crystals span over all the settling regimes, with 0.5 mm crystals being in the dust-like regime for ρ p /ρ f up to 1.5, and 10 mm crystals being in the stone-like regime for ρ p /ρ f larger than 1.1 (Fig. 18b). Again, the expected density contrasts result in settling times significantly smaller than the life-span of the system [typically more than 100 kyr, see e.g. the review 62], indicating fractional crystallization.
The position of the studied particles with respect to the slow belt has important consequences, as it is directly related to the horizontal distribution upon sedimentation. Thus, for a mantle-deep magma ocean, the crystal radius must be 10 mm in order to experience a horizontally non-uniform accumulation of sediments (see the black triangles that fall into the transitional regime in Fig. 18a). For a 3-km-deep magma chamber, on the other hand, the overlap between [0.5, 10] mm crystals and the slow belt indicates that the majority of suspended particles will eventually settle in the low-velocity piles (under the assumption that a large scale circulation is present and the low-velocity piles form, see also Discussion).
Note that the rhythmic sedimentation suggested by Sparks et al. [38] intrinsically relies upon the assumption that precipitated crystals settle much faster in a non-convecting fluid than when convection is present. The existence of the slow-belt presented here is thus in favour of the scenario proposed by Sparks et al. [38], for which the authors found petrological evidence in fully solidified chambers.

VI. DISCUSSION
The crystallization of a primordial molten mantle is a complex system in which the generation, settling, and re-entrainment of crystals are competing processes. Here we only focused on one of these components: the settling of crystals.
Typical time scales for the solidification of a whole-mantle terrestrial magma ocean range from ∼ 10 3 years in the absence of atmosphere, up to ∼ 10 6 years in the presence thereof [e.g. 32,33]. When compared to these time scales, our results indicate a fast settling (Fig. 18a), and thus support the idea of a fully fractional crystallization. In the series of papers by Solomatov [for a review, see 8], it is instead argued for equilibrium crystallization of the majority of the primordial mantle. This is largely because re-entrainment of sedimented particles from the bottom of the fluid is claimed to be the dominant process. Solomatov et al. [63] derive a formula for the equilibrium crystal fraction, Φ eq = 18 αρνQ/(gc p ∆ρ 2 r 2 c ), where Q is the surface heat flux and denotes the fraction of available convective energy that goes into re-entrainment, estimated to be ≈ 0.1 − 1 % [63]. For the simulations presented here, e.g. for Ra = 10 10 , P r = 50 and crystals with ∆ρ/ρ = 0.1 and r c = 1 cm, the resulting Φ eq is only around 3 % (the formula for Φ eq is designed for a single type of crystals only; for a range of crystal properties it must be decided how much of the available energy goes into the lifting of the various types). For the significantly larger heat fluxes that accompany the early stages of a global magma ocean solidification, it quickly reaches 100%, indicating full suspension [64].
More recently, the scaling law of Solomatov et al. [63] was confirmed by the experimental study of Lavorel and Le Bars [19], who systematically varied the density ratio ρ p /ρ f and the temperature contrast that drives thermal convection. An important finding of their study is that in a highly turbulent flow the molecular viscosity that appears in the formula for the Stokes' velocity must be replaced by an apparent viscosity in order to account for turbulent eddies smaller than the particles. Such approach was successfully used to describe the dynamics of finite-size particles in turbulent flows [65,66], but it is difficult to generalize in non-homogeneous flows such as for thermal convection between parallel plates. In the context of particle settling, Lavorel and Le Bars [19] found this approach viable, and the apparent viscosity that they measured at Ra = 3 × 10 9 was only ca. 2.7 times larger than the molecular viscosity. We note, however, that replacing molecular viscosity with apparent (or turbulent) viscosity in their formula for the decay of the number of suspended particles (their Eq. (9)) is similar to dividing the terminal velocity |v t | by f as in our Eq. (22), resp. it is mathematically identical if f can be treated as constant for the investigated particles. Moreover, the amplitude 2.7 is within the range that we obtain for f . In other words, the decrease of the Stokes' velocity that Lavorel and Le Bars [19] computed to reconcile their measurements may have been caused by the effects of large-scale circulation analyzed in this paper as well as by the effect of sub-particle sized turbulence.
Similarly to the previous works, Lavorel and Le Bars [19] used a cuboid tank filled with salty water and spherical, polymethyl methylacrylate particles. We note that the conditions of the experiments [5,63,67] on which the energetic analysis of Solomatov and Stevenson [36] is based differ significantly from the environment of non-spherical silicate crystals that accumulate at the bottom of a cooling magma. In nature, the sedimented crystals may be subject to chemical and petrological altering, possibly binding the crystals together, i.e. making them prone to re-entrainment. Note that fractional crystallization is often reported in exposed plutons [e.g. 38]. While we do not argue against re-entrainment as such, we merely point out that its workings should be thoroughly investigated in future work in the context of magma environments.
In fact, the non-uniform horizontal distribution of settling events that we observe in the transitional and bi-linear regimes is slightly in favour of re-entrainment. As discussed in Solomatov et al. [63], embedded particles may be lifted by the tangential stresses caused by rising plumes, and the respective stresses increase with distance from the domain boundary. We observe a large concentration of sedimented particles in the low-velocity piles. Inside these regions, the sediments would thus build tall dunes. Since most of young plumes are born at the edges of the low-velocity piles, the crests of these dunes should be exposed to large tangential stresses.
In terms of the predicted settling time, some of our results differ dramatically from those of Verhoeven and Schmalzl [40]. Roughly speaking, for particles satisfying (ρ p − ρ f ) < ρ f α∆T (see their Fig. 12 and Eq. 5), they obtain a temperature-dominated convection mode (T-regime), in which the flow is thermally driven and all particles are held indefinitely in suspension (similar results are obtained in the non-rotating cases of Maas and Hansen [43,44]). This is because in their formulation the momentum equation is solved for the volumetric average of the fluid and particle velocity, i.e. their particles by definition follow the fluid as described in Eq. 6 in Verhoeven and Schmalzl [40] (but they can also invoke fluid motion, see the next paragraph). Therefore, in case of vigorous thermal mixing of the fluid, their particles never settle. In our study, particles can have different velocities from the surrounding fluid, which allows fluid-particle separation and thus sedimentation regardless of convective vigor.
For stronger density contrasts, roughly for (ρ p − ρ f ) > ρ f α∆T , Verhoeven and Schmalzl [40] obtain a particle-driven convection in which a layer of sediment is segregated from the rest of the fluid (C-regime). We note that the condition (ρ p − ρ f ) > ρ f α∆T simply means that the critical concentration that is required for the formation of a settling front, as described earlier by Koyaguchi et al. [6] and Sparks et al. [38], is less than 100%. In other words, the C-regime is established whenever the formation of a settling front (and thus cessation of convection due to particle motion) can take place for some critical particle concentration C * < 100 %. As analyzed by Solomatov and Stevenson [36], for the crystals of interest in magma oceans and chambers, C * is typically less than 100% (and this is also the case for most of the particle types investigated in the present paper).
Within the C-regime, Verhoeven and Schmalzl [40] develop a model that is based on the theory of Martin and Nokes [5]. In particular, they complement the theory by accounting for the volume occupied by particles that have already sedimented. This is not to be confused with our model, in which the function f is a measure of the rate of particle transport into the low-velocity piles, normalized by the probability of not escaping from these regions. As such, our function f depends non-trivially on the structure of the background flow, while the factor f in Eqs. 26-32 of Verhoeven and Schmalzl [40] represents packing of sedimented material. Verhoeven and Schmalzl [40] then verify their model on a set of simulations in which the employed particles fall into the dust-like regime according to our classification (see their Eq. 8, Fig. 10, and the parameters listed below Fig. 15). This should be understood as yet another confirmation of the applicability of Eq. (15) for particles with a small value of the |v t |/u rms ratio.
At the beginning of our simulations, we inject the particles uniformly throughout the entire model domain. In systems where the particulate phase is a product of a chemical reaction or phase change, this is typically not the case. In particular, in a cooling magma, the solid crystals nucleate in the relatively cold downwellings and degassing takes place due to decompression in hot rising plumes. The importance of the initial positions of newly formed particles depends on their v t /u rms ratio. On one hand, in the dust-like regime the particles are likely to get thoroughly mixed and the settling rates would not be affected. In the stone-like regime, on the other hand, the particles' trajectories and settling times strongly depend on the particles' initial positions and velocities. Given the particle positions associated with the second stage of the bi-linear regime (Fig. 5), faster settling than reported here is to be expected on average if the majority of heavy crystals form preferentially in downwellings.
For Ra = 10 12 and P r = 50, i.e. for the highest investigated Reynolds number, the Kolmogorov length scale η is ≈ 2×10 −3 , while the upper bound of the non-dimensional particle radius r c /H is 1.5×10 −3 (i.e. the smallest turbulent eddies are only slightly larger than the investigated particle radii in the respective simulation set). For the simulation set xC 10 50 the assumed crystal sizes even exceed the Kolmogorov length, and the particle Reynolds numbers are in turbulent rather than laminar regime (i.e. breaching the range of validity of Eq. (9)). We note, however, that the set xC 10 50 is performed only for illustrative purposes in the context of explaining some aspects of the stone-like regime.
Our study shows that the existence of a stable large-scale flow structure has a clear signature on the settling of particles. First, it delays the settling on average. Second, it is responsible for the non-uniform horizontal distribution of settling events. Since the large-scale circulation is observed at the highest Rayleigh numbers reached so far in numerical simulations [68] and experiments [37], we speculate that this can be an important feature for the extreme regime of a cooling magma ocean. Superstructures were analyzed in detail in 3D systems with large aspect ratios by Pandey et al. [69]. To test the influence of a larger aspect ratio, we performed the reference simulation C 10 50 also with aspect ratio 4, and the resulting function f as well as the amplitude of the horizontal variations of settling events were nearly identical.
In Section V we performed an extrapolation to Ra = 10 27 , assuming Re = 0.07 Ra 0.582 P r −0.873 . Grossmann and Lohse [61] distinguish four idealized regimes of convection depending on whether kinetic and thermal dissipation rates are dominated by the convective bulk or the boundary layers. For our Rayleigh and Prandtl numbers, the idealized case ("pure power-law") yields Re ∝ Ra 4/9 P r −2/3 , but already for P r ≈ 1 we would be on the boundary with the Re ∝ Ra 1/2 P r −1/2 regime, in which the kinetic energy dissipates in the velocity boundary layer of the flow. While these are idealized cases, real convection is a mixture of these regimes, and for a detailed treatment one must employ the full theory of Grossmann and Lohse [61]. Using scaling laws to extrapolate to very high Rayleigh numbers, such as done in Section V, is, however, still subject to an open debate [e.g. 37].
It is interesting to note that the amplitude of u rms does not depend on Ra and P r in the limit of Re ∝ Ra 1/2 P r −1/2 , i.e. in the idealized scenario for low Prandtl and high Rayleigh numbers. In such a case, the settling behaviour could be predicted simply by computing the terminal velocities of particles of interest since the v t /u rms ratio would not depend on the exact values of Ra and P r (see Eq. (24)). However, for low values of P r new settling regimes may exist, which we plan to investigate in the future.
Convection in magma oceans occurs in the presence of rotation, which is not included in our model. Using estimated values for Earth's magma ocean, the convective Rossby number lies in a range of 0.03 -100 [43], and the large-scale circulation thus may be disrupted due to rotation [70,71]. Both the settling rates and the horizontal distribution of sedimented material would be affected in such scenario. Recently, the effects of rotation on the distribution of crystals and the rate of settling in the early stages of Earth's primordial magma ocean were analyzed in spherical geometry by Maas and Hansen [44]. In a rotation dominated scenario (Rossby 0.3) they find a much pronounced settling as convection and thus vertical entrainment of particles is suppressed.

VII. SUMMARY
We evaluate the settling rate of inertial particles that are injected into statistically steady state thermal convection. Previously, the number of suspended particles in such system was assumed to follow either the relation N = N 0 exp(−v t t) or N = N 0 (1 − v t ), with v t being the Stokes' velocity. We observe a new regime with particularly slow sedimentation, in which large-scale circulation prevents particles from reaching the boundary layers of the fluid. By introducing a new framework that treats the settling mechanism as a random process, we develop a model that unifies the observed settling rates into the general equation N = N 0 exp(−v t t/f ), where f is a function of the ratio of Stokes' and mean characteristic velocity of the flow, v t /u rms . We investigate f over a broad range of Reynolds numbers and show that the function is relatively robust. It reaches its maximum for v t /u rms ≈ 0.4, the maximum value ranging approximately from 1.5 to 3 for particles with mild density contrasts with respect to density of the fluid (Table IIa).
We also analyze the horizontal distribution of settled particles. Within the regime of slow settling, heavy particles accumulate preferentially below major clusters of upwellings. These are located at edges of largescale convection rolls.
For Reynolds numbers larger than ∼ 1000 and particles with a stronger density contrast, additional complexity arises because of the preferential concentration phenomenon, i.e. light particles have a tendency to move towards flow vortices, while heavy particles move away from them. As a result, light particles get captured inside long-lived vortices, which significantly prolongs their sedimentation at the top boundary.
The maximum value of f is close to 10 (resp. the normalized settling time F ≈ 30) for our simulation with the highest Rayleigh number, Ra = 10 12 .
When extrapolated to the extreme conditions of solidifying magma chambers and oceans, our results predict fractional crystallization. For a better understanding of such complex systems, it is possible (and necessary) to extend our method to account for 3D geometry, rotation, re-entrainment of sedimented particles, self-consistent nucleation of solid crystals, and the coupling between particle concentration and momentum conservation of the fluid.

VIII. ACKNOWLEDGEMENT
We thank Tina Rückriemen for useful discussions. VP and NT acknowledge support from the Helmholtz Association (project VH-NG-1017) and the DLR-DAAD research fellowship program.