The fate of particles in a volumetrically heated convective fluid at high Prandtl number

The dynamics of suspensions plays a crucial role on the evolution of geophysical systems such as lava lakes, magma chambers and magma oceans. During their cooling and solidification, these magmatic bodies involve convective viscous fluids and dispersed solid crystals that can form either a cumulate or a floating lid by sedimentation. We study such systems based on internal heating convection experiments in high Prandtl fluids bearing plastic beads. We aim to determine the conditions required to produce a floating lid or a sedimented deposit. We show that although the sign of particles buoyancy is the key parameter, it is not sufficient to predict the particles fate. To complement the model we introduce the Shields formalism and couple it with scaling laws describing convection. We propose a generalised Shields number that enables a self-consistent description of the fate of particles in the system, especially the possibility to segregate from the convective bulk. We provide a quantification of the partition of the mass of particles in the different potential reservoirs (bulk suspension, floating lid, settled cumulate) through reconciling the suspension stability framework with the Shields formalism. We illustrate the geophysical implications of the model by revisiting the problem of the stability of flotation crusts on solidifying rocky bodies.


Introduction
According to the classical scenarii of planetary formation, terrestrial bodies were likely partially or totally molten, forming a magma ocean (Taylor & Norman 1992;Tonks & Melosh 1993;Abe 1995Abe , 1997. This initial stage in planetary history is due to two major phenomena. In the first few million years of the solar system, during the accretion of planetesimals, the decay of short lived radioactive elements such as 26 Al and 60 Fe is an important heating source (Urey 1955;Neumann et al. 2012;Weidenschilling 2019;Kaminski et al. 2020). Later, collisions between planetary embryos and giant impacts converted gravitational energy into heat and produced massive melting events (Tonks & Melosh 1992;Safronov & Ruskol 1994;Canup & Asphaug 2001). During the cooling of such a system, the temperature evolves from the liquidus to the solidus, a melting interval which is several hundred of degrees for silicate systems. If the cooling rate controls the 2 C. Sturtz,É. Kaminski, A. Limare and S. Tait. pace of solidification of the system, solidification and the fate of crystals can also introduce a feedback on the thermal history of the system. For example, the anorthosite crust of the Moon formed by flotation of light plagioclase crystals (Wood 1972;Warren 1985;Shearer et al. 2006). This thick crust has an insulating effect on the convective system underneath (Lenardic & Moresi 2003;Grigné et al. 2007) hence slowed down its thermal evolution (Elkins-Tanton et al. 2011;Maurice et al. 2020).
If the fluid were quiescent, crystals would settle down or float according to the sign of their buoyancy. Convection may prevent this behavior by maintaining particles in suspension. The crucial issue lies in the determination of the stability of such a suspension. Sparks et al. (1984) pointed out that the Stokes velocity of particles in magma chambers, which stands for the typical settling velocity, is small compared to the mid-depth vertical root-mean square fluid velocity. In that case, suspensions in magmatic reservoirs should always be stable, as particles behave like passive tracers. However, Martin & Nokes (1988) illustrated experimentally that negatively buoyant particles initially in suspension eventually settle down and form deposits. Lavorel & Le Bars (2009) furthermore highlighted that deposition occurs at Stokes velocity even though convection is turbulent. These observations can be interpreted considering the interaction between particles and the dynamical boundary layers that develop at the borders of the reservoir, where velocities vanish because of rigid boundary conditions (Sparks et al. 1984). As a matter of fact, dealing with suspension sustainability requires a criterion that involves both convection and sedimentation.  proposed a description based on the energy balance of fluid-particle interaction. These authors assumed that the fluid can transfer an amount of its convective energy to particles. If the suspension gravitational energy that drives settling exceeds this quantity, the suspension is not stable and deposits form. Even though = 0.1 − 1 % had been evaluated from experiments Lavorel & Le Bars 2009), it has been underlined that this ad-hoc parameter is not well constrained ). If it were, this model could make possible a totally self-consistent mass budget, without involving any ad-hoc parametrization.
Instead of dealing with suspension stability, other authors prefer to consider the stability of the beds of particles that may form. Erosion and entrainment of particles is usually described by the formalism proposed by Shields (1936). Basically, particles are locked on the bed surface because of frictional forces that have the same order of magnitude as the buoyancy force according to Coulomb's law. Entrainment is possible if the ratio of the stress acting on particles over the particles buoyancy exceeds a critical value that lies between 0.1 − 0.2 (Charru et al. 2004). This model was used to determine sediment transport upon bedloads (Lajeunesse et al. 2010) and the equilibrium height of the settled bed (Leighton & Acrivos 1986). However, these studies refer to experiments that involve controlled flow, with isoviscous and isothermal conditions, that are not consistent with the geophysical applications considered here. Convection in magmatic reservoirs involves destabilization of thermal boundary layers (TBL) that complexifies both temperature and velocity fields.  adapted the previous reasoning by considering that a single particle in the TBL can not be lifted, but is moved horizontally at the surface of the bed by the tangential stress. Then, particles form dunes that make the stress quasi-vertical, enabling entrainment. However, these authors do not study the equilibrium thickness of the underlying bed, or its influence on the thermal state.
In the present study, we aim at quantifying the behavior of such a lid embedded within the TBL of a convective system. We will consider the case of an internally heated system cooled from above, displaying only one TBL which is at the upper surface of the Erosion and deposition in an internally convective fluid.
3 convective layer, and we will focus on the stability and the equilibrium thickness of a floating lid. The first part of the study tackles the dimensionless equations that outline the problem, and underlines the key dimensionless numbers that describe the system. In the second part, we develop an experimental approach to study this issue. The setup used is composed of a tank containing the convective fluid internally heated by microwave absorption and plastic beads that represent crystals. We examine the erosion of the floating lid, and we introduce a model that predicts the equilibrium thickness of the crust according to the thermal state. We identify and test experimentally the stability condition to form a cumulate. We emphasize one dimensionless parameter, the Shields number, as the key parameter to deal both with the lid stability and with the suspension sustainability. Finally, we discuss a geological case illustrating the relevance of our model: the condition of stability of a floating lid on terrestrial bodies.

Internally heated convective systems
Geophysical problems considered here involve convective systems driven by internal heating and secular cooling, phenomena that are mathematically equivalent. Thermal convection in the Boussinesq approximation is governed by the following equations representing the conservation of momentum, energy and mass (see, e.g.: Jaupart & Mareschal 2010, pp.114):

2)
∇ · u f = 0, (2.3) where u f is the fluid velocity field, ρ 0,f is the reference fluid density, η f is its dynamic viscosity, α f is its thermal expansion, c p,f is its specific heat, λ f is its thermal conductivity, g is the acceleration of gravity, H is the rate of internal heat generation, P is the pressure, T is the temperature field, and θ = T − T is the thermal anomaly with T the horizontal average temperature.
Internally heated convective systems are characterized by the internal temperature scale ∆T H defined as: where h is the vertical thickness of the system. This temperature scale introduces a new definition of the Rayleigh number, called the Rayleigh-Roberts number (Roberts 1967): where κ f = λ f /ρ 0,f c p,f is the fluid's thermal diffusivity. This dimensionless number enables the characterization of the vigor and patterns of convection (e.g.: Vilella et al. 2018).
Using h as the length scale, ∆T H as the temperature scale, W = ρ 0,f α f g∆T H h 2 /η f as the velocity scale that represents the Stokes' velocity of a laminar thermal, h/W as the time scale, and η f W/h as the pressure scale, one can provide dimensionless form of the Boussinesq equations as follows (see, e.g.: Jaupart & Mareschal 2010, pp.114):

6)
Ra H ∂T ∂t + u f · ∇T = ∇ 2 T + 1, (2.7) ∇ · u f = 0, (2.8) where P r is the Prandtl number, which compares viscous and thermal diffusion: where ν f = η f /ρ f is the kinematic viscosity. The Prandtl and Rayleigh-Roberts numbers characterize the regime of convection occurring in the system. They scale inertia in equation (2.6). The high-P r low-Ra H limit corresponds to laminar flows, whereas low-P r high-Ra H yields turbulent inertial flows. In geophysical systems, these dimensionless numbers span a large range of values. For instance, solid-state convection in the current Earth's mantle verifies P r ≈ 10 23 and Ra ≈ 10 7 , whereas magmatic reservoirs are rather evolving with P r ≈ 10 3 − 10 8 and Ra ≈ 10 11 − 10 16 . In the case of magma oceans, huge variation of P r and Ra H over the thermal history are expected. Massol et al. (2016) estimated for the terrestrial magma ocean a Rayleigh number that goes from 10 31 at the very beginning if the planet is totally molten, to 10 14 when crystals are in suspension. This drop can partly be explained by the fact that the presence of crystals increases the apparent viscosity of the mixture by several orders of magnitude when the rheological transition is reached (Lejeune & Richet 1995;Guazzelli & Pouliquen 2018). As a consequence, the P r number increases drastically during the solidification, from an initial value around 10 1 − 10 2 to the current value of 10 23 for solid mantle convection. In comparison, experiments carried out in the present study lie in the following ranges: P r ≈ 10 3 and Ra H ∈ [3.10 6 , 10 8 ]. According to the theory by Grossmann & Lohse (2000a,b), partially crystallized magma oceans and our experiments occur in the same regime of convection.
Internally heated convective systems are characterized by a single upper TBL generating cold instabilities. The temperature drop across the TBL ∆T T BL and the velocity of downwellings W i have been studied experimentally and numerically, and can be expressed based on local scaling analyses (Limare et al. 2015;Vilella et al. 2018): (2.10) where the pre-factors C T and C W depend only on the mechanical boundary condition at the top of the system (see : table 1). Introducing particles into the convective fluid add significant complexity in the problem. Indeed, depending on the density, shape and size of the particles, but also depending on the fluid flow, multiple phenomena are likely to occur, which are described in an abundant literature (Andreotti et al. 2011;Boyer et al. 2011;Houssais et al. 2015). We summarize below the theoretical framework and fundamental parameters that will later be used to model these interacted phenomena.

Particles in suspension -two-phase flow
The presence of crystals dispersed in a magma reservoir makes it a two-phase system. To describe the dynamics of two-phase flow, a set of complementary equations has to Erosion and deposition in an internally convective fluid.

5
Conditions C T C W
be added to the conservation equations (2.1)−(2.3) to describe the two phases and the interactions between them (see, e.g.: Andreotti et al. 2011, pp.306). Considering the fluid phase (f ) and the solid phase (p), and assuming that the volume fraction φ of the solid phase is uniform, i.e. no chemical or mass exchanges between the phases, equations of motion for the fluid and the particles are respectively: where ∆ρ = ρ f − ρ p is the density difference between the fluid and the particles, and f is the fluid-particle interaction force. All other forces are neglected, such as Van der Waals interactions or frictional forces between particles. There are many approaches describing the interaction force f depending on the flow regime and particle properties (Maxey & Riley 1983). Here, we consider that at high P r number, the interaction between the fluid and particles is dominated by the viscous drag, which is written as: where r is the particle radius and β(φ) is a dimensionless function that refers to the contribution of the other particles to the drag. It increases as φ increases (Andreotti et al. 2011, pp.306).
Using the same scales as before, the dimensionless set of equations becomes: Two other dimensionless numbers appear. The C parameter compares the gravitational potential energy of particles to the inertial drag: (2.17) In our experiments, this parameter is of order unity. The Stokes number characterizes the interaction between the fluid and particles: For reservoirs much larger than the size of particles, which is a limit relevant for magma oceans or magma reservoirs, and/or for laminar flow, the Stokes number is likely to be much smaller than unity. In this case, particles are statistically passive tracers, following fluid motions (e.g.: Crowe et al. 2011, pp.25). Consequently, (2.16) becomes: We emphasize that this limit describes the average behavior of particles but does not imply that particles never settle. Over a short timescale compared to the convective timescale, the particles are indeed passive tracers and the equality (2.19) holds true. Nevertheless, particles are still buoyant, so there is a small component of the particles velocity that participates to settling. In this way, equation (2.19) rewrites u p = u f + u s , with ||u s || ||u f ||. Thus, sedimentation actually occurs with a low probability (Patočka et al. 2020), and deposits form at long timescales compared to the convective one.
A direct consequence of particles sedimentation is the formation of settled cumulates or floating lids. This implies in turn that erosion and particles re-entrainment from these layers must be taken into account in the modeling of convective systems.

Settling and re-entrainment
Erosion and re-entrainment from settled cumulates and/or floating lids bring particles back in suspension . However, the framework that describes this phenomenon is different from the one used to treat suspensions as it depends on local mechanical equilibrium of particles (Charru et al. 2004;Lajeunesse et al. 2010). Particles at the surface of the bed are submitted to two forces: (i) the frictional force between the particles and the underlying bed that captures particles at the surface of the bed and is proportional to the particles buoyancy according to Coulomb's law and (ii) the shear stress induced by the flow. A dimensionless number, called the Shields number, compares these two forces (Shields 1936): where τ is the shear stress at the surface of the bed. This ratio enables the definition of a critical value ζ c that describes the threshold behavior of particles on the bed. If ζ < ζ c , particles are locked on the bed by frictional forces, whereas if ζ > ζ c , the shear stress is strong enough to erode particles. For spherical plastic particles homogeneously sheared by a viscous, laminar flow, Charru et al. (2004) estimated ζ c = 0.12. The force driving particles (re-)entrainment does not take into account any vertical pressure effects. This follows the comment made by  that a single particle can not be lifted by a vertical pressure gradient by comparing the pressure force exerted on the particle to its buoyancy. The authors proposed that the mechanism for entrainment of particles is strongly linked to shear stress acting at the interface, that manage to displace particles and form dunes that enables entrainment.
The framework presented above can be used to describe the dynamics of a suspension and the coupled stability of cumulates and/or floating layers. To our knowledge, this coupling has not been studied yet in internally heated convective systems relevant for magma reservoirs. To fill this gap, we present below lab-scale experiments.
Erosion and deposition in an internally convective fluid.

Convection with internal heating
Achieving experimental convection driven by homogeneous internal heating at high Rayleigh numbers was challenging until Limare et al. (2013) developed a unique experimental set-up based on microwave absorption. A 30 × 30 cm 2 wide and 5 cm high tank is introduced in a specially designed microwave oven (Surducan et al. 2014). The top of the tank is a thermostated aluminium plate whose temperature is fixed and monitored. The other walls and the base of the tank are made of poly(-methyl methacrylate) (PMMA), transparent to visible light and microwave radiation, and ensuring adiabatic thermal boundary conditions.
A laser sheet scans half of the tank (15 cm), and we acquire images at a spacing of 1 cm. Two CCD cameras register images in different spectral ranges allowing noninvasive measurement of the temperature field via a two-dye Laser Induced Fluorescence (LIF) method. The velocity field is measured by Particle Image Velocimetry (PIV). The temperature and velocity spatial resolution are 0.2 and 0.8 mm respectively. Further details on the experimental setup and calibration can be found in Fourel et al. (2017). The same set-up and methods are used in the following study, but the fluid is adapted to study the sedimentation of beads. Typical 2D velocity and temperatures field are shown in figure 1 for experiments without beads. Panels (a.1) and (b.1) show the 2D horizontal and vertical velocity fields and their correspondent root mean square (RMSvertical profiles (panels ((a.2) and (b.2)). The velocity is zero on the boundaries since they are all rigid. Negative vertical velocities are associated with thermal instabilities generated at the top boundary of the tank. Positive vertical velocities corresponds only to return flow. Figure 1 (c.2) displays the temperature vertical profile showing that the thermal structure of the convective layer can be split into an upper boundary layer and a convective interior. An important feature is that the fluid interior has slightly negative Properties are all measured in the lab, except those marked with (*) which are taken from Mark (2007). See Supplementary Materials for further information on the way properties measurements have been carried on.

Fluid and particles
The working fluid used in experiments is a mixture of 44 wt% glycerol and 56 wt% ethylene glycol. Particles are PMMA spherical beads. Two sets of beads are used, corresponding to two different radii (r 1 = 290 µm, r 2 = 145 µm). The main properties are summarised in table 2. Particles have a different thermal expansion coefficient α p than the fluid allowing the investigation of the full range of particle behaviours. For both phases, the density is linked to the temperature according to the thermal equation of  (table 3). In some cases, due to progressive heating of fluid, beads become heavier than the fluid and settle to form a cumulate. state: where T 0 is the reference temperature, ρ 0,i is the reference density at the reference temperature T 0 , and i refers to the fluid or particles. In this case, the density difference between the fluid and particles is: If the fluid is cold enough, particles are lighter than the fluid and float, whereas at higher temperature, beads become heavier and can sink. Thus, an inversion of buoyancy exists at an "inversion temperature" T inv (figure 2). Furthermore, when the thermal state in the experimental tank spans a large range of temperature, from the cold surface temperature T s < T inv to the bulk temperature T bulk > T inv , the system can display simultaneously both a floating lid and cumulate formation. In our case, T inv = 37.4 o C. Experimental conditions are summarized in table 3. Fluid Prandtl number is high (P r ≈ 1000) and experiments reached high Rayleigh-Roberts numbers (Ra H ∈ [3.10 6 , 10 8 ]). Particles Stokes number is about 10 −5 − 10 −4 , which makes them passive tracers.
As particles are made of PMMA, they are transparent to microwave radiation, so internal heating only occurs in the fluid. Comparing the diffusive time scale in one particle τ d ∼ r 2 /κ p and the convective time scale τ c ∼ h/W leads to τ c /τ d ≈ 10 1 − 10 2 , which shows that the thermal equilibration of particles is rapid. Thus, we assess that the local temperature difference between the particles and the fluid is negligible.
Unfortunately, in this configuration, we could not achieve refractive index matching between the beads and particles (see Supplementary Materials). Consequently, plastic beads are light-scattering objects and the images recorded by the cameras are blurred. This has little influence on the velocity measurement, as particles behave like passive tracers. In order to check if it affects the temperature measurements, we carried out a sedimentation experiment using an homogeneous suspension in a controlled isothermal environment. We monitored the fluorescent signal whilst particles settled, and confirmed that the mean temperature measured was consistent with the imposed temperature. Hence, the presence of beads does not affect the measurement of the average properties of convection, on which further reasoning is based.

Experiments
Experiments were conducted as follows. The tank was filled with a mixture of fluid and particles. One has to avoid introducing air into the system, in order to limit surface tension effects (see Supplementary Materials for more details). The system is thermostated at the surface temperature T s generally bellow T inv , so that particles form initially a floating lid (see figure 3). Then, the microwave power is turned on and convection starts within a time lapse of few tens of seconds. As convection proceeds, the floating lid is eroded, and particles are put in suspension. In some experiments, the inversion temperature was reached, and particles could further form a cumulate (Table  3). We waited until the thermal steady state was reached (which corresponds to a period of time spanning from 2 to 6 hours).
In the following section, we will study in detail these two aspects: the erosion of the floating lid and the formation of the cumulate.

Thermal steady state
The presence of a floating lid is likely to influence the thermal state of the system as it is situated in the TBL (figure 4). In the following section, we will quantify the thickness of the steady lid and the thermal state of the system. Figure 5 reveals that the dimensionless drop of temperature between the bulk and the surface is systematically higher than the one predicted by the scaling laws (2.10). We displayed the scalings for both mechanical conditions (rigid and free-slip) since this boundary condition is not well defined beneath a granular lid. The lid reduces the efficiency of heat transfer that occurs at the top of the reservoir and causes an insulating effect related to the lid thickness.
To quantify this effect, we use the model developed to study convection under a conductive lid (Guillou & Jaupart 1995;Lenardic & Moresi 2003;Grigné et al. 2007). The reasoning is based on two hypotheses. First, we consider the lid as a homogeneous conductive layer with an average conductivity λ, and an averaged thickness δ, such that the heat flow through the lid is: 10 7 ) Table 3. Experimental characteristics: the beads radius r, the initial floating bed thickness δ 0 , the imposed surface temperature T s , the mean bulk temperature at steady state T bulk , the Rayleigh-Roberts number calculated at steady state Ra H . The two last columns inform about the erosion of the floating lid (whether partial or total), and the formation of a cumulate at the bottom of the tank. Note that two families of beads with different radius are used. The last 8 rows correspond to experiments done without beads.
where T lid is the basal temperature of the lid and Q s = Hh is the surface heat flux at steady state. The lid's mean conductivity is estimated based on the individual values of each component weighted by their respective volume fraction: λ = φ RLP λ p +(1−φ RLP )λ f , with φ RLP the beads packing in the lid, assumed to be the random loose one: φ RLP = 56%. We further consider that the scaling laws for thermal convection hold true beneath the conducting lid. It has been verified experimentally for homogeneous internally heated systems that the drop of temperature between the surface and the bulk T bulk − T s scales like the drop of temperature through the TBL ∆T T BL given in (2.10) with a pre-factor C T = 3.38 ± 0.16 at steady state (Limare et al. 2019). Thus, thanks to the continuity of   Table 1).
temperature between the lid and the fluid, one can get: which combined with (4.1) yields: showing that the larger the lid thickness, the hotter the bulk. To validate this model that uses the thermal state of the system to estimate the lid thickness, an independent way to measure δ is required. In that aim we explore the Erosion and deposition in an internally convective fluid.  Consequently, by symmetry, the vertical velocity field has a maximum value at middepth -figure 6 (a.2). Also by symmetry, the horizontal profiles have a minimum value at mid-depth between two maxima corresponding to 1/4 and 3/4 of the total thickness of the tank as can be seen in figure 6 (a.1). This shape is due to the confined environment that creates a horizontal recirculation flow. As emphasized in figure 6 (b.1) and (b.2), the presence of the floating lid shifts vertically the point where the horizontal velocity is minimal by an amount ∆δ, which is set by the "mechanical" thickness of the lid δ m : Comparison between the thickness deduced from the velocity-shift method δ m and the inverted thermal thickness δ th calculated thanks to (4.3) is shown in figure 7. In this plot, we only use experiments where the lid is partially eroded and where no cumulate appears at steady state, in order to limit unwanted shifts due to the presence of a deposit of particles at the base of the tank. The agreement between the measurements and the model is fair, despite some scatter due to the simplifications of the model used. We assessed that the floating lid can be approximated by an homogeneous conductive lid whilst it is composed of packed particles containing interstitial fluid. Moreover, its thickness is not strictly uniform as dunes form during erosion/deposition processes (see Supplementary Materials). In the following we will use the thermal thickness δ th to characterize the lid thickness.
Besides, the symmetry observed in the horizontally averaged vertical profile U x,RM S is clearly related to the return flow due to the rigid lateral boundaries and is not a general feature of convection driven by internal heating. For instance, in the case of a spherical shell, this symmetry has a priori no reason to exist. This experimental feature was used as an additional measurement of the floating lid thickness. This assumption does not affect the validity of the reasoning as our model does not require a symmetry of the lateral flow. . Comparison between the lid thickness inverted from the thermal state δ th and the one measured by the velocity-shift method δ m . Error-bars on δ m correspond to one pixel of the PIV-grid, and those on δ th correspond to one particle radius.

Local mechanical equilibrium
As illustrated in figure 4 (a), local equilibrium of the beads is set by the balance between erosion forces and the bead buoyancy. To quantify such an equilibrium, we rely on the threshold theory of mechanical erosion (Glover & Florey 1951;Métivier et al. 2016) and we define the Shields number ζ lid at the base of the floating lid as: where τ = η fγ is the characteristic convective shear acts on the bottom of the lid, andγ the corresponding strain rate. We consider only the experiments with partial erosion of the floating lid at steady state, which means that the lid thickness δ = 0. In these experiments, the Shields number reaches the critical threshold (ζ lid = ζ c in (4.5) and δ = δ th in (4.1)). At steady state, we assume that the temperature filed in the floating lid varies linearly with depth, so that the temperature at the base of the lid is T lid = T s + Q s δ th /λ. Introducing this temperature in (3.2), some algebra yields the following set of equations: where ζ s is the Shields number calculated at the surface temperature T s . With (4.6), ζ s appears as the control parameter that determines the critical thickness δ c , as T s and Q s are known. The problem thus boils down to the determination of the characteristic convective shear stress τ . For instance, in experiments done by Charru et al. (2004), the shear is experimentally controlled and homogeneously applied to the bed, which facilitates its characterization. Here, the bed lies in the unstable cold TBL. Thus the flow field is Erosion and deposition in an internally convective fluid.  Figure 8. Scaling laws derived from experiments without beads. Horizontal (a) and vertical (b) velocity, horizontal (c) and vertical (d) strain rate. All these physical properties are evaluated thanks to their root mean square value calculated in the entire volume of the tank and are represented by the dots. Blue lines represents best fit law with fixed exponent 3/8 and the orange dashed line are those with the exponent left to vary. All parameters of these laws are summarised in Table 4. complex and contains spatial and temporal fluctuations. Hence, characterizing the shear stress in a homogeneously heated convective system is required.

Scaling laws for velocities and shear stresses
Equation (4.5) emphasizes the importance of the horizontal shear stress, hence of the velocity field, on the erosion process. The strain rateγ scales as follows: where U L is the characteristic horizontal velocity, and δ v is the characteristic length over which velocity vanishes. By definition, the latter corresponds to the dynamical boundary layer (DBL) thickness. First, the Reynolds numbers Re reached in our experiments are low which implies laminar flows and thus, δ v ∼ h (see appendix A). As a consequence, the strain rateγ becomes:γ (4.10) We consider the volume of fluid in the TBL that is drained by one downwelling (Figure  4 (b)). On one side, fluid from the TBL is drained at the characteristic velocity W i by the downwelling whose cross section area is A i . On the other hand, fluid is converging at the characteristic horizontal velocity U L through the lateral surface of the cylinder of thickness δ T BL and area ∆S. This reasoning is based on the fact that the fluid drained Using the same scalings as in Vilella et al. (2018): H , valid for 10 6 Ra H 10 9 , we obtain the horizontal velocity scale: with C u a pre-factor which depends on the boundary conditions. To verify these scaling laws, experiments without beads have been carried out using the same fluid and the same methods as those described previously. We recorded the horizontal and vertical velocities using the PIV method, and we calculated the horizontal and vertical strain rate. As we are interested in the average shear rate, we determined their RMS values calculated over the entire volume of the tank. Results are displayed in figure 8 (a), (b) for the horizontal and vertical velocities. The scaling laws pre-factors determined experimentally are summarized in table 4. In our experiments, Re ≈ 10 −1 −1, so the convection is laminar. The scaling law for the strain rate is therefore C γ H . Results are displayed in figure 8 (c) and (d) for the horizontal and vertical strain rates respectively. In both cases, the predicted scaling law is in good agreement with experimental data.
We can thus estimate the Shields' number as follows: (4.13)

Critical Shields number and stability of deposits
Thanks to the previous scaling analyses, equation (4.6) provides a way to measure experimentally the critical Shields number. By considering that the lid thickness at steady state corresponds to the critical thickness δ th = δ c , one can determine the critical Shields number ζ c for each experiment: (4.14) with ζ s given by (4.13). The critical number is calculated for each experiment, and results are displayed in figure 9. We get: ζ c = 0.29 ± 0.17. This value is of the same order of magnitude of the one estimated by Charru et al. (2004). Error bars represent a variation of one bead radius of the floating lid thickness. The discrepancy is due to the sensitivity of the prediction to the value of the lid thickness, and the higher the heat flux, the steeper the thermal gradient in the lid and the greater the error bars. This is the reason why one experiment in particular (IHB13) does not appear in the plot because the uncertainties are too large.
With the critical number ζ c determined, we can compare the stability of the different deposits predicted by the Shields approach with the experimental observations (figure 10 (a) and (b)). For the floating lid, we calculate the surface Shields number ζ s . If ζ s > ζ c , the floating lid is unstable and erosion should be total. But if ζ s < ζ c , a floating lid of some thickness is stable. Results are displayed in figure 10 (a) and the transition between total erosion and partial erosion is well described by ζ c . Similarly, we calculate the bulk Shields number ζ bulk , which is also the Shields number at the base of the reservoir: and we compare it to ζ c . If ζ bulk > ζ c , convection prevents deposition at the base of the tank. Otherwise, a basal deposit is stable. This transition is also verified experimentally in figure 10 (b).

Suspension stability
The scalings derived above allow predicting whether a deposit can form at the base of the reservoir and a stable lid at the top of it, but we need another framework to describe the full dynamics of the suspension containing the particles eroded from the deposits. In a fluid at rest, particles with negative buoyancy will all eventually settle down. Observations show that in a convective fluid, even negatively buoyant particles can remain in suspension at steady state (Lavorel & Le Bars 2009).  proposed that the dynamics of the suspension can be described based on the equilibrium between buoyancy and shear forces, described by the balance: or in compact form, where B p is the integral on the left-hand side of (5.1) related to beads buoyancy, V f is the integral on the right-hand side of (5.1) referring to the bulk viscous dissipation, (V ) is the total volume of the reservoir, and is the percentage of viscous energy being used to maintain particles in suspension (also called efficiency parameter ). This description can be further used to determine the maximal concentration of particles φ max that can be maintained in suspension by convection. Assuming to be constant,  get the following law for φ max : where τ and ∆ρ stand for the volume averaged values of τ and ∆ρ respectively, and C s = 9/2 . Basically, if the concentration of particles in the convective bulk φ is below this limit, particles stay in suspension. Otherwise, the convective fluid only sustains the quantity of particles corresponding to the maximal concentration φ max , and the rest settles, forming a cumulate. Interestingly, this law can be expressed with a Shields parameter ζ similar to the one used previously. The major difference lies in the fact that ζ is a global parameter, comparing volume average properties, whereas previously ζ has been estimated locally. However, by neglecting the thin TBL at the top of the reservoir, one can approximate that ζ ≈ ζ bulk = ζ(T bulk ), and thus φ max is described by previous scaling laws.
In order to verify the validity of this criterion in our experiments, we considered the formation of the cumulate at the base of the tank. As the direct measurement of particle concentration φ is not possible due to the refractive index mismatch, we calculate a proxy ϕ of the quantity of particles that is eroded from the floating lid. To do so, taking δ 0 as the initial thickness of the bed, and δ th as the thickness at steady state, the quantity of particles that is eroded is proportional to δ 0 − δ th . The coefficient of proportionality is linked to the packing of particles inside the floating lid. The packing is assumed to be constant, as all experiments are prepared similarly. Thus, the concentration of particles in the bulk is calculated as if all the eroded material stay in suspension: where a is a constant which depends potentially on the packing of beads inside the lid (a = 0.56 for a random loose packing). As the measurement of the total quantity of deposit particles is subject to large uncertainties, we detect φ max as the limit between partial sedimentation and absolute suspension defined by the limit φ max ∼ ζ 2 . The blue line in figure 11 represents the empirical boundary between the two regimes. Now, we have a complete framework that describes the steady state of any suspension in a convective environment driven by internal heating. The Shields's formalism enables us to quantify the thickness of the lid located in a boundary layer, and Solomatov's approach enables us to describe the maximal quantity of particles that can stay in suspension.

Lids formation in a convective system bearing particles
Our results can be used to describe the fate of particles in a convective fluid by splitting the system into three reservoirs: (i) the floating lid situated in the TBL, (ii) the bulk fluid containing suspended particles, and (iii) the deposit at the base. Our model quantifies the volume of particles that can be stored in each reservoir (figure 12). Buoyant crystals first fill the floating lid reservoir. We can define the maximal capacity of the floating lid V c determined by (4.6). Particles exceeding this critical volume of the lid remain in suspension or form a cumulate. In the same way, for negatively buoyant crystals, if the bulk Shields number is below the critical value, a deposit may form. The concentration of crystals that stay in suspension is limited to φ max , the rest settles and forms the basal deposit.

Validity of the experimental results
The approach adopted here takes into account the main mechanisms that deal with suspension and deposits stability. The main goal was to highlight how to use the framework treating of fluid/crystal interaction in the case of convective systems bearing 20 C. Sturtz,É. Kaminski, A. Limare and S. Tait. Figure 12. Synopsis that fully determines the system's regime at steady state. The surface Shields number ζ s is given by (4.13) and depends on the surface temperature T s , the Rayleigh-Roberts number Ra H and the initial volume of particles V 0 which are all input parameters. V c is the volume of the critical floating lid that can form at the surface of the convective fluid, corresponding to a lid of thickness δ c given by (4.6). Stability of the floating lid: if ζ s > ζ c , the floating lid is unstable and all particles are put in suspension. If ζ s < ζ c , the erosion is partial. In this case, if V 0 < V c , the initial crust is stable, and very few particles are put in suspension. If V 0 > V c , only the volume of particles V c stay at the surface, the rest of the volume V sus = V 0 − V c is put in suspension. Cumulate formation : we consider the stability of the suspension in the case of negatively buoyant suspended particles. If the basal Shields number ζ bulk is greater than ζ c , the suspension is stable and no cumulate forms. If ζ bulk > ζ c , a basal cumulate can settle down. The volume of this cumulate is given by the maximal concentration of particle that can stay in suspension φ max given by (5.3). The convective fluid can only sustain a maximal volume V max = φ max V in suspension, where V is the volume of the convective layer. V max has to be compared to the volume of particle that have been eroded V sus . The surplus V d = V sus − V max forms a basal cumulate.
crystals. It is based on hypotheses that allows to express the problem in terms of scaling laws that can easily be re-scaled to conveniently describe a geophysical system as it will be illustrated in next section.
Even though experimental results are well described by our models at first order, certain discrepancies exist and should not to be ignored. They result from different factors. First, our experiments enable the estimation of the floating lid thickness by an indirect measurement based on the average temperature state of the system. This method might increase the uncertainty on the determination of δ and, thus, all parameters that are linked to it, such as the entrainment threshold. One way to improve this measurements could be the use of fluid and beads that have a near-perfect optical index matching. In this way, the bottom of the lid would be directly observable, and the thickness measurable more precisely. This prerequisite enabled the study of the bed surface motions and even of the motions of beads inside a granular bed (Mouilleron et al. 2009;Houssais et al. 2015). In this way, the bottom of the lid will be directly observable, and the thickness measurable more precisely.
Second, our theory is based on averaged physical properties, which leads to the general predicament of describing a 3D problem by a 1D theory. This crude approximation leads to results that are fairly consistent with data at first order, but a local description of the convection and the erosion mechanism can improve this approach. Downwellings develop from the TBL, which can have an influence on the local behavior of beads at the bottom of the floating lid. Equivalently, we describe the floating lid as a quiescent solid medium of homogeneous thickness, but this assumptions is a strong simplification. We point out the existence of a topography (see supplementary materials) that can have an influence on the local flow in the TBL, which, in return, impacts the erosion mechanism that occurs at the interface (Charru & Hinch 2006). Local recirculation flow that induces dune formation are intimately linked to downwellings in our case, and understanding the interaction between dune growth and downwelling dynamics can upgrade the present study. Furthermore, the suspension is also considered in our model as evenly distributed in the bulk. This assumption represent only a particular case and numerical simulations underline that crystals can concentrate depending on the flow field and their buoyancy (Patočka et al. 2020). Further investigations are required to observe experimentally the behavior of particles in suspension as a function, for instance, of the Shields number, with a set-up that allow the study of a large range of values for ζ.
Dealing with the erosion mechanism, we only treated the threshold in terms of one critical value of the Shields' number. Our estimate is consistent with data but the deviation from this value can be due to an oversimplifications. First, as timescales to reach the steady state are long, the steady state might not be reached for the erosion of the floating lid when we stopped the experiment. This might induce an over-estimate of the equilibrium thickness of the bed, and consequently an under-estimate of the threshold value. Besides, the Shields number may depend on parameters such as the grain size and the flow regime. Buffington (1999) showed for instance a variation for turbulent flow between 0.03 − 0.1, whereas Ouriemi et al. (2007) underlined that the threshold is stable in laminar flow but it depends on the way the bed is packed (Agudo & Wierschem 2012). The latter authors measured changes in the critical Shields number of a factor 2 up to a factor 5 depending on the exposure of the beads to the flow and the packing and orientation of the substrate. Thus, a dependance of ζ c on other physical parameters related to the flow (such as the Ra H ), or related to the substrate (such as the packing) has to be verified with an experimental device that enables a study over several orders of magnitude.
All these effects could explain the outliers that appears in our data. Nonetheless, our results are proposed here as a first step in the path of using the physical framework used in granular media in order to describe the behavior of particles in convective systems.

Characterization of fluid/particle coupling in our experiments
We now illustrate a way to draw a parallel between the two-phase flow framework and the above energy balance approach developed in section 5. The purpose is to show qualitatively how the efficiency parameter could be refined and how it relates to the different regimes that are reached in our experiments. This parameter is supposed to depend on many factors such as fluid and particles densities, the particle size, the concentration and the heat flux involved ) but this dependency has not yet been clarified. In a convective system without particles, there is a bulk equilibrium between the total viscous dissipation V f and the total buoyancy of the fluid B f . In the presence of particles, part of the viscous dissipation corresponds to the work done to maintain the particles in suspension: with: Combining (2.12) with (7.1) leads to: This relation links quantitatively the efficiency parameter to the coupling force between the fluid and particles. In our case, for laminar convection, the coupling force (2.14) can be used and we obtain: As discussed earlier, as St 1, particles are statistically passive tracers. Nevertheless, we observed experimentally the formation of cumulates, so there is a weak component of particle motion that enables settling. In this way, we assume the following decomposition for the particle velocity: where u s stands for the settling velocity, which is considered to be the Stokes velocity, and verifies ||u f || ||u s || in our case. This hypothesis is also corroborated by the experiments done by Lavorel & Le Bars (2009), who showed that particles in turbulent convection settle at a speed that scales with the Stokes velocity. Equation (7.4) thus suggest the Erosion and deposition in an internally convective fluid.
To verify this expression in our experiments, we rewrite (5.3) by substituting with (7.7) and deduce a new scaling law for the maximal concentration of particles in the suspension: As a consequence, the proxy ϕ of the particle concentration should also follow this law: ϕ = C ϕ h/rζ with C ϕ a constant that is constrained experimentally C ϕ = 8.10 −4 ( Figure  13). To verify qualitatively the consistency of this expression with , we compare (5.3) and (7.8) to get an expression for the efficient coefficient: Assuming that particles are packed randomly in the floating lid (a = 0.56), we get the average value of = 4 ± 3% in our experiments. This estimate is slightly higher than the values reported in the literature but still consistent with them Lavorel & Le Bars 2009). However, the main purpose of this physical argument is to show how to reconcile Solomatov's ad-hoc expression with the complete physical framework that governs the two-phase flow and to verify it qualitatively. Experiments where the bulk concentration of beads is precisely measured are required to corroborate quantitatively this reasoning.

Crystal segregation in magma ocean
Based on the model developed above, we propose an insight into the segregation process of particles in the thermal history of a magma ocean, and especially the flotation of plagioclase. This type of crystals is known to be lighter than the residual liquid from which they form (Elkins-Tanton et al. 2011). The flotation of plagioclase is the main scenario invoked to explain the formation of the light anorthosite crust of the Moon (Wood 1970;Solomon & Longhi 1977). Convection driven by secular cooling is the relevant regime for convective magma ocean. Secular cooling is mathematically equivalent to internal heating treated in our model (see, e.g.: Krishnamurti 1968). Deschamps et al. (2012) demonstrated that scaling laws that describe thermal convection in volumetrically heated cartesian boxes is also valid in spherical geometry, provided we take into account a geometrical factor describing the shell curvature. In this case, the surface heat flux is linked to the secular cooling and the internal heating rate as follows: No stable crust  This study Planetesimals 200 km LMO Figure 14. Diagram of flotation versus suspension for two planetary bodies. Each solid line corresponds to the critical radius r c given by (7.11). Particles with radius greater than the critical one are allowed to float, otherwise they stay in suspension. Dashed line show the critical radius of crystals predicted by Solomatov's law (7.12).
where f = (R − h)/R is the ratio between the depth of the magma ocean h and the planetary radius R. Further details can be found in Appendix B.
Q s ranges from 10 6 W m −2 for molten planetary bodies that release heat by radiation to space (Elkins-Tanton et al. 2011;Massol et al. 2016), down to 10 −3 − 10 −2 W m −2 when the flotation crust is present (Maurice et al. 2020). All physical quantities used in the model are summarized in table 5. We deduce the critical particle radius enabling crystal flotation for two types of magma ocean: a 200 km deep shallow lunar magma ocean (R = 1737 km, f = 0.88, g = 1.6 m s −2 ) and a fully molten planetesimal (R = 300 km, f = 0.42, g = 0.23 m s −2 ). This critical radius r c is calculated from the critical Shields number (4.13): Results are displayed in figure 14. As crystal size during crystallization is estimated in the range r = 0.1−10 mm (Solomatov 2000), we deduce that crystals float during magma ocean cooling. Magma ocean episodes enable efficient crystal segregation. This segregation is less efficient for smaller bodies such as planetesimals (with R 300 km), suggesting that crystallization history might be more complex in these cases.We compared our results to those displayed in Elkins-Tanton (2012) based on the law of : (7.12) with η = 0.1 Pa s. (7.12) was considered by these authors to be applicable in both laminar and turbulent experiments but established for system heated from below. Here we propose a model for convective systems driven by secular cooling. Our model enables smaller particles to participate to crustal formation.

Crustal thickness and flotation efficiency
Our model not only provides criteria for the formation of stable crystal reservoirs, floating crust or cumulates, but can also be used to quantify the partitionning of the crystals between the suspension and these deposits. We illustrate how to quantify this partition of crystals between the deposits and the convective bulk in the particular case of the lunar magma ocean (LMO).
In the literature, it is widely assumed that all plagioclase crystals settle instantaneously once they nucleate in the LMO (Elkins-Tanton et al. 2011;Maurice et al. 2020), which yields a system composed of a flotation crust growing above a purely liquid magma ocean. The main results of these two studies are indicated in the two highlighted bands in figure  15. This hypothesis was relaxed using the flotation efficiency which is the ratio between the amount of crystals that float and the total volume of crystals that nucleate in the system. This parameter can be estimated from field data. For instance, Namur et al. (2011) determined the value for plagioclase flotation by studying the Sept Iles layered intrusion, and they found an efficiency of 50%. They assumed that this flotation efficiency would be the same in the case of the LMO, and they use it to constrain the amount of plagioclase that composed the anorthosite crust and its thickness. Charlier et al. (2018) used the crustal thickness measured by spacecraft gravity data (Wieczorek et al. 2012) and a flotation efficiency in the range between 40 and 100% in order to constrain the thickness of the LMO.
This ad-hoc parameter can be discussed in the light of our model. According to the present study, if the initial volume of crystals is smaller than the critical volume that can be sustained at the surface, all crystals form a flotation crust. If the crust has reached its maximum thickness, the crystals "in excess", i.e. that can not be incorporated in the crust, remain in suspension. This is a fundamentally different way of thinking which links crystal flotation to physical parameters.
We consider that the LMO is a shallow shell corresponding to 20 vol% of the total volume of the Moon. This high level of solidification is necessary to trigger plagioclase crystallization (Elkins-Tanton et al. 2011). The evolution of the maximal thickness sustainable at the surface of the LMO is displayed in figure 15. The crystal size has almost no influence on the steady state crustal thickness except on the low limit. This is due to the sub-critical value of the Shields number, which is essentially constrained by the low value of the melt viscosity. In order to form a crust with a thickness of 40 km corresponding to current estimates of the lunar crust, our model predicts that the surface heat flux should be below 3 W m −2 . In the cases of Elkins-Tanton et al. (2011) and Maurice et al. (2020), 100% flotation efficiency would imply a thicker crust. Our model adds some physical constraints on the LMO evolution in general, by coupling quantitatively the segregation of crystals to the thermal history.
We emphasize that our reasoning is dealing with stability of the crust and its steady state. The predicted crust thickness is reached after a time-dependent deposition process which is not discussed in this paper. It requires a complete description of the transient regime of the flotation crust.

Limitation of the model
We aim to illustrate how the criterion developed in this paper can be applied to natural systems and how they can provide new constraints on conditions under which a deposit can form or not. Of course, our illustration is a crude approximation of the complexity of magma ocean. We use simplified hypotheses, but the framework we propose here is robust enough to enable refinements that can be added at will in order to include complexity.
Our first strong assumption is the bulk crystallization that was abundantly considered in the literature as a common end-member (Solomatov 2000;Elkins-Tanton et al. 2011). This is justified by the nature of convection driven by secular cooling, which implies that crystals nucleate within the cold downwellings that end up in an isentropic convective bulk. For small-size bodies, this is equivalent to an isothermal interior.
We also assume that crystals are non-cohesive, as in experiments. This hypothesis is compatible with batch crystallization where crystals nucleates in the convective bulk. This regime of crystallization is consistent with magma oceans, but may not be for smaller magmatic reservoirs such as magma chambers or lava lakes. In these cases, crystals may mainly form at the boundary of the system, directly where they are supposed to settle. This adds a complexity to the system as cohesion between crystals becomes important. One way to take cohesion into account in the present model may be by increasing the erosion threshold ζ c . The stronger the interaction between crystals are, the higher the value of ζ c .
In addition, we did not consider in detail the effect of crystallization series on the evolution of composition and its consequences on all the physical parameters such as the number of crystals, crystal size and liquidus and solidus temperatures. We expect that in limit of small crystal fractions, our results would still apply, but the transient behavior of crystallizing systems needs to be studied specifically. Moreover, our model does not take into account the influence of pressure on all the physical parameters. This hypothesis is consistent with small planetary bodies such as planetesimals, or with the latest stages of shallow magma oceans, but it is less justified for larger bodies such as the Earth. In the latter case, pressure would affect the whole thermodynamics, as it influences the geotherm, the solidus and liquidus of crystals, the viscosity of the liquid magma, the formation of crystals, etc. In these cases, our stability criterion based on the Shields number is still valid, but the threshold value and the scaling laws that are involved should be adapted accordingly.
Furthermore, we described the system by adopting stability criterion as did Solomatov. But magma oceans are evolving systems where temperature varies a lot during the thermal history. Hence, describing their complex behavior over time thanks to a stationary criterion is a strong assumption that has to be relaxed. Complete thermal model that takes into account the kinetics of sedimentation and erosion is required to deal with suspension stability in magma oceans. Erosion/deposition mechanism commonly used in geomorphology (Charru et al. 2004;Lajeunesse et al. 2010) could be further revisited in terms of crystal erosion and sedimentation in magmatic reservoirs.
We considered here the high Prandtl limit where inertia is negligible compared to viscous forces. Nonetheless, magma oceans have experienced turbulent episodes (Solomatov 2000) that imply high Rayleigh numbers that are not yet attainable experimentally or numerically. This seems to be all the more true for large bodies. In this case, scaling laws used for the thermal state have to be adapted. For instance, Kraichnan (1962) expressed the different scaling law for turbulent flow, depending on the Prandtl number. It would lead to other scaling laws for the shear stress and thus other expressions for the Shields number. But fundamentally, the bases of our approach would remain true even in the turbulent case.

Conclusion
Particles sheared in a convective fluid can form deposits, floating crusts or cumulates, depending on the sign of their buoyancy. Using laboratory experiments, we proposed a model that estimates the partitioning of particles in such systems at high Prandtl number and for Rayleigh-Roberts number up to 10 9 . Our model is based on the estimation of one dimensionless parameter that encapsulates the main physical ingredients required to describe the system: the generalized Shields number ζ that compares the buoyancy of particles to the shear stress generated by convection. The value of this parameter quantifies both the thickness of the flotation crust, and the maximal concentration of crystals that can be sustained in suspension. The volume of particles that forms a cumulate can be deduced from these two pieces of information by a mass budget. This unifying framework can be adapted to understand transient episodes of the thermal history of natural systems driven by secular cooling and/or internal heating. In order to do so we need to complete the model by the detailed description of the dynamics of erosion and deposition, and in particular of the characteristic timescales that are involved. These timescales have to be compared to the convective system cooling timescale. Further studies are needed to completely understand the feedback existing between these intricate timescales. This is a necessary step in order to explain the wealth of the observed structures in planetary bodies and in their vestiges that arrive to us as meteorites.    which leads to δ v ∼ h Re −1/2 where Re = ρ f hU L /η f is the Reynolds number. This relation justifies quantitatively that the dynamical boundary layer occupies the entire reservoir in the laminar limit, which is the regime reached in our experiments ( figure 16). The scaling law that characterizes velocity in internally heated convective systems can be rewritten in terms of Reynolds number Re using (4.11): which has been verified experimentally as illustrated in figure 16. The pre-factors for horizontal and vertical Reynolds number are given in table 6.

Appendix B. Secular cooling and internal heating
In geophysical systems evolving over long periods of time, the thermal state is transient, so one can argue that our reasoning is not applicable as scaling laws that are used here hold only at steady state. Here, we will show that it is possible to adapt them to describe the transient state.
First, we treat the floating lid as a homogeneous layer where Fourier's law can be applied. In the transient state, the heat flux departs from (4.1), which is valid only in steady-state conditions. Hence, we can define the diffusive time scale in the floating lid τ diff = δ 2 /κ, with κ = λ/ρ c p the thermal diffusivity of the floating lid where λ, ρ and c p are respectively its average thermal conductivity, density and specific heat. This time scale is compared to the time scale related to the evolution of the thermal state of a convective fluid heated from inside. The latter is given by: Limare et al. 2019). Thus, the criterion that ensures that the lid's thermal state evolves quasi-statically is based on the ratio: If R τ 1, the lid reaches its steady state quickly and (4.1) holds true. In our experiments, R τ ≈ 0.1 − 0.5, so the approximation is justified. In geophysical systems, the approximation holds true if Ra H is moderate, and/or if the floating lid is thin compared to the depth of the convective mantle. For instance, for lids that represent δ/R = 1% equivalent to the ratio between the lithosphere thickness and the Earth's radius, the quasi-static approximation is correct for Ra H up to 10 12 , which is relevant for geophysical applications. Second, the model is based on scaling laws that describe the thermal state of convective systems heated internally. In the transient state, the secular cooling term can be considered as an additional source of internal heat. We define the modified rate of internal heat generation as: In this way, the Rayleigh-Roberts number is also modified: Thus, the drop of temperature across the TBL becomes: Limare et al. (2021) showed that this model describes well the transient state in experimental convective systems. The authors measured C * T = 3.58 ± 0.15 similar to the value retrieved from experiments in homogeneous, steady state conditions. This model was further used to describe the thermal evolution of parent bodies of iron meteorites (Kaminski et al. 2020). It shows that convective systems where convection is mainly driven by secular cooling can also be treated in the framework presented in this paper.
Consequently, we are able to estimate the evolution of the Shields number from the thermal history of such systems, enabling the study of the formation of floating crust and/or deposits.

Fluid
The fluid used in this study is a mixture of 44wt% glycerol and 56wt% ethylene glycol. Density and thermal expansion have been determined thanks to a DMA 5000 Anton Paar densimeter. Thermal properties, such as thermal diffusivity κ f , thermal conductivity λ f and specific heat c p,f , have been determined thanks to a photothermal method described by Dadarlat and Neamtu [2009]. The viscosity η f and its dependance on temperature have been monitored thanks to an RS600 Anton Paar rheometer. The refractive index n f and its temperature dependence have been obtained thanks to the Abbemat 350 Anton Paar refractometer. These properties are shown in table 1 and their tempertaure dependence is plotted in figure 2 (a), (b) and (c).

Beads
As beads are small solid spheres, the measure of previous properties is much more complicated. To measure the density, we introduced the beads in a tank containing a stratified fluid. This density stratification is achieved by the double bucket method described by Fortuin [1960] and illustrated in figure 3. The buckets contain two mixtures of different proportions of water and glycerol having different densities, called (A) and (B) for the high and low value respectively. As the heavy mixture (A) is pumped into the working tank (S), the lighter one (B) is introduced in (A) thanks to a basal linking tube, reducing the density of (A) and, hence, the density of the fluid deposed in the tank (S).
Step by step, as the pumping is going on, the density of fluid introduced in (S) becomes less important. The stratification so created is sampled at different heights, and controlled apart, thanks to the DMA 5000 Anton Paar densimeter. Beads are introduced in the stratified liquid, and they stop at a level corresponding to their neutral buoyancy heigh (NBH). The determination of this height is a way to measure their density.
In order to measure the thermal expansion α p , this type of set up has been carried on for different thermostated environments. Achieving a stable and isothermal stratified environment is extremely complicated, because even a difference of a few degrees between the wall of the tank and the bulk might trigger convection, and thus, destroy the stratification. In order to prevent this effect, the fluids (A) and (B) have been prepared in a thermostated room, whose temperature can be tuned within a certain range (between 18 and 27 o C). In this way, we get three measurements of density at three different temperatures (Figure 2 (a)), and assuming a constant thermal expansion, we get the following value for beads : α p = 3.17.10 −4 K −1 .
We work with two families of beads of different radius. Their size has been determined by image analysis of a sample of hundreds of beads of each family. The size distribution is plotted in Figure 2 (d), and it leads to the estimation of the averaged radius : r 1 = 290 µm and r 2 = 145 µm.
All the other thermal properties are taken in Mark [2007] : λ p = 0.21 W/m/K, c p,p = 1765 J/kg/K and thus κ p = 1.10 −7 m 2 /s. Concerning the refractive index, the handbook exposes a wide range of values, from 1.445 for bead to 1.56 (Mark [2007], p.830). As we noticed a strong index mismatch in experiments, the refractive index of the beads is likely to be closer to the upper bound of this range.

How to introduce beads without bubbles ?
Achieving an appropriate experimental setup that involves a fluid bearing beads is actually challenging. Especially, we have to avoid introducing air bubbles in the mush, as it might induce unwanted surface tension effects, and, thus, increase cohesion between particles.
To do so, we systematically applied the following protocol. First, the working tank is filled with the fluid. Then, we prepared a mixture of fluid and beads in the home-made stir-machine -see Figure 5 (a) and (b). The stir-machine includes a vessel and a rotor that is composed of an horizontal grid aiming to maintain the beads in the lower part of the stir-machine, and let   the bubbles escape. This protocol is done at room temperature, therefore beads float and form a lid against the grid. In order to prevent settling and to maintain the beads in suspension, three blades are also attached to the rotating grid. Degassing the mixture consists in stirring the fluid with beads, and letting it quiescent in order to let the bubbles escape. This procedure is repeated several times, until no more bubbles escaped from the separation grid. Then, the stir-machine is connected to the working tank. Isostatic pressure enables the beads to go from the stir-machine to the tank. The flow has to be controlled in order not to let particles be entrained and escape the working tank. Once the amount of particles introduced is sufficient, the tank is closed, shaken to homogenize the mixture, and placed into the oven. We waited until particles settle before beginning any experiment. The total beads mass introduced in the tank is measured at the end of the set of experiments, when the tank is drained, the mixture is sieved, and the beads are washed, dried and then weighed. Figure 5 shows two pictures of the deposits studied : (a) the floating crust that forms at the surface, and (b) the cumulate that appears when some conditions enable its formation. Picture (a) has been taken at the end of an experiment, when all particles have settled at the surface. Note the egg-box-like shape of the floating lid, very similar to shape that can be observed in geomorphology and in experiments by . Concerning the basal cumulate (d), note that deposits form patterns like dunes. We must emphasize that these dunes move and deform as the convection goes on.  Figure 4 -Stir-machine that enables the introduction of beads without air bubbles. (a) The grid in the rotor has a mesh smaller than the beads diameter, that retains them in the lower part of the cylinder while letting the bubbles escape. (b) Set up to introduce beads in the tank in the vertical position. The grid and the rotor are immersed in the fluid. Blades have also been fixed on the rotor, in order to keep beads in suspension, to avoid their settling against the grid and to allow bubbles detachment and escape through the grid. 5 Data for the scaling laws dealing with strain rates and velocities Refractive index mismatch between the fluid and beads affected the velocity field as soon as the laser sheet scanned larger distances from the front wall and the concentration of particles became significant. We performed 8 experiments without any beads in order to determine accurately the scaling laws for the velocity and the strain rate. We measured the steady-state velocity field in order to calculate either the characteristic velocities (vertical and horizontal) and the characteristic strain rate (vertical and horizontal). All these values are evaluated by their root mean square (RMS) values calculated over the entire volume of the tank. The results are shown in Table 2. These data were used to determine the scaling laws shown in section 4.2.2.

Fluid
The fluid used in this study is a mixture of 44wt% glycerol and 56wt% ethylene glycol. Density and thermal expansion have been determined thanks to a DMA 5000 Anton Paar densimeter. Thermal properties, such as thermal diffusivity κ f , thermal conductivity λ f and specific heat c p,f , have been determined thanks to a photothermal method described by ?. The viscosity η f and its dependance on temperature have been monitored thanks to an RS600 Anton Paar rheometer. The refractive index n f and its temperature dependence have been obtained thanks to the Abbemat 350 Anton Paar refractometer. These properties are shown in table 1 and their tempertaure dependence is plotted in figure 2 (a), (b) and (c).

Beads
As beads are small solid spheres, the measure of previous properties is much more complicated. To measure the density, we introduced the beads in a tank containing a stratified fluid. This density stratification is achieved by the double bucket method described by ? and illustrated in figure 3. The buckets contain two mixtures of different proportions of water and glycerol having different densities, called (A) and (B) for the high and low value respectively. As the heavy mixture (A) is pumped into the working tank (S), the lighter one (B) is introduced in (A) thanks to a basal linking tube, reducing the density of (A) and, hence, the density of the fluid deposed in the tank (S).
Step by step, as the pumping is going on, the density of fluid introduced in (S) becomes less important. The stratification so created is sampled at different heights, and controlled apart, thanks to the DMA 5000 Anton Paar densimeter. Beads are introduced in the stratified liquid, and they stop at a level corresponding to their neutral buoyancy heigh (NBH). The determination of this height is a way to measure their density.
In order to measure the thermal expansion α p , this type of set up has been carried on for different thermostated environments. Achieving a stable and isothermal stratified environment is extremely complicated, because even a difference of a few degrees between the wall of the tank and the bulk might trigger convection, and thus, destroy the stratification. In order to prevent this effect, the fluids (A) and (B) have been prepared in a thermostated room, whose temperature can be tuned within a certain range (between 18 and 27 o C). In this way, we get three measurements of density at three different temperatures (Figure 2 (a)), and assuming a constant thermal expansion, we get the following value for beads : α p = 3.17.10 −4 K −1 .
We work with two families of beads of different radius. Their size has been determined by image analysis of a sample of hundreds of beads of each family. The size distribution is plotted in Figure 2 (d), and it leads to the estimation of the averaged radius : r 1 = 290 µm and r 2 = 145 µm.
All the other thermal properties are taken in ? : λ p = 0.21 W/m/K, c p,p = 1765 J/kg/K and thus κ p = 1.10 −7 m 2 /s. Concerning the refractive index, the handbook exposes a wide range of values, from 1.445 for bead to 1.56 (?, p.830). As we noticed a strong index mismatch in experiments, the refractive index of the beads is likely to be closer to the upper bound of this range.

How to introduce beads without bubbles ?
Achieving an appropriate experimental setup that involves a fluid bearing beads is actually challenging. Especially, we have to avoid introducing air bubbles in the mush, as it might induce unwanted surface tension effects, and, thus, increase cohesion between particles.
To do so, we systematically applied the following protocol. First, the working tank is filled with the fluid. Then, we prepared a mixture of fluid and beads in the home-made stir-machine -see Figure 5 (a) and (b). The stir-machine includes a vessel and a rotor that is composed of an horizontal grid aiming to maintain the beads in the lower part of the stir-machine, and let the bubbles escape. This protocol is done at room temperature, therefore beads float and form   a lid against the grid. In order to prevent settling and to maintain the beads in suspension, three blades are also attached to the rotating grid. Degassing the mixture consists in stirring the fluid with beads, and letting it quiescent in order to let the bubbles escape. This procedure is repeated several times, until no more bubbles escaped from the separation grid. Then, the stir-machine is connected to the working tank. Isostatic pressure enables the beads to go from the stir-machine to the tank. The flow has to be controlled in order not to let particles be entrained and escape the working tank. Once the amount of particles introduced is sufficient, the tank is closed, shaken to homogenize the mixture, and placed into the oven. We waited until particles settle before beginning any experiment. The total beads mass introduced in the tank is measured at the end of the set of experiments, when the tank is drained, the mixture is sieved, and the beads are washed, dried and then weighed. Figure 5 shows two pictures of the deposits studied : (a) the floating crust that forms at the surface, and (b) the cumulate that appears when some conditions enable its formation. Picture (a) has been taken at the end of an experiment, when all particles have settled at the surface. Note the egg-box-like shape of the floating lid, very similar to shape that can be observed in geomorphology and in experiments by ?. Concerning the basal cumulate (d), note that deposits form patterns like dunes. We must emphasize that these dunes move and deform as the convection goes on.  Figure 4 -Stir-machine that enables the introduction of beads without air bubbles. (a) The grid in the rotor has a mesh smaller than the beads diameter, that retains them in the lower part of the cylinder while letting the bubbles escape. (b) Set up to introduce beads in the tank in the vertical position. The grid and the rotor are immersed in the fluid. Blades have also been fixed on the rotor, in order to keep beads in suspension, to avoid their settling against the grid and to allow bubbles detachment and escape through the grid. 5 Data for the scaling laws dealing with strain rates and velocities Refractive index mismatch between the fluid and beads affected the velocity field as soon as the laser sheet scanned larger distances from the front wall and the concentration of particles became significant. We performed 8 experiments without any beads in order to determine accurately the scaling laws for the velocity and the strain rate. We measured the steady-state velocity field in order to calculate either the characteristic velocities (vertical and horizontal) and the characteristic strain rate (vertical and horizontal). All these values are evaluated by their root mean square (RMS) values calculated over the entire volume of the tank. The results are shown in Table 2. These data were used to determine the scaling laws shown in section 4.