J Monteux
email: julien.monteux@uca.fr

B Qaddah

D Andrault

Conditions for segregation of a crystal-rich layer within a convective magma ocean

Keywords: Early Earth, thermal evolution, magma ocean, numerical modeling, two-phase flow, convection

. Hence, during the early crystallization of a magma ocean, suspension should be the dominant process.

We then investigate the implications of our models on the mechanism of magma ocean solidification after a major melting event on an Earth-like planet.

We define the relevant set of parameters; including the P-V-T equations of state of coexisting melt and bridgmanite in the mushy MO. We observe that bridgmanite grains are unlikely to segregate in a mantle of pyrolite composition.

Introduction

During its early history, the Earth has experienced several episodes where a significant volume fraction of its mantle was molten. These melting events were enhanced by radiogenic heating [START_REF] Yoshino | Core formation in planetesimals triggered by permeable flow[END_REF], viscous dissipation during core formation [START_REF] Monteux | A model of metal-silicate separation on growing planets[END_REF][START_REF] Samuel | Heterogeneous lowermost mantle: Compositional constraints and seismological observables[END_REF] and kinetic energy converted into heat during large impacts [START_REF] Tonks | Magma ocean formation due to giant impacts[END_REF]. The giant impact leading to the Earth/Moon system could even be responsible for the complete melting of the Earth mantle [START_REF] Nakajima | Melting and mixing states of the Earth's mantle after the Moon-forming impact[END_REF] and the subsequent formation of a nearly 3000 km thick magma ocean.

Extremely vigorous convection is likely to occur within such a molten reservoir [START_REF] Solomatov | Magma Oceans and Primordial Mantle Differentiation[END_REF][START_REF] Monteux | On the cooling of a deep terrestrial magma ocean[END_REF]. The variations in temperature are strong enough to generate a significant natural convection flow. Moreover the potentially large thickness of a molten-silicate reservoir coupled with low viscosity material [START_REF] Karki | Viscosity of MgSiO 3 Liquid at Earthś Mantle Conditions: Implications for an Early Magma Ocean[END_REF][START_REF] Xie | Direct Viscosity Measurement of Peridotite Melt to Lower-Mantle Conditions: A Further Support for a Fractional Magma-Ocean Solidification at the Top of the Lower Mantle[END_REF] leads to Rayleigh numbers (comparing advective and diffusive timescales) ranging between 10 20 to 10 30 [START_REF] Patočka | Settling of inertial particles in turbulent Rayleigh-Bénard convection[END_REF] while the current Rayleigh number of the Earth mantle is ranging between 10 6 and 10 8 [START_REF] Ricard | Physics of Mantle Convection[END_REF]. The cooling and solidification dynamics are difficult to monitor within such an environment. However, recent numerical models agree on the characteristic cooling timescales that is constrained to several thousands of years [START_REF] Lebrun | Thermal evolution of an early magma ocean in interaction with the atmosphere[END_REF][START_REF] Monteux | On the cooling of a deep terrestrial magma ocean[END_REF].

The study of particle sedimentation in vigorously convecting fluids and, by extension, the study of crystal behaviour in convecting magmas is a debated subject in fluid dynamics [START_REF] Lavorel | Sedimentation of particles in a vigorously convecting fluid[END_REF] and in Earth Sciences [START_REF] Jaupart | Dynamics of differentiation in magma reservoirs[END_REF]. A large number of study has been dedicated to the crystal settling in magma chambers using laboratory experiments [START_REF] Martin | Crystal settling in a vigorously converting magma chamber[END_REF][START_REF] Sturtz | The fate of particles in a volumetrically heated convective fluid at high Prandtl number[END_REF], numerical modelling [START_REF] Höink | Formation of compositional structures by sedimentation in vigorous convection[END_REF][START_REF] Verhoeven | A numerical method for investigating crystal settling in convecting magma chambers[END_REF] or both [START_REF] Weinstein | Evolution of crystalsettling in magma-chamber convection[END_REF]. All these studies point out the importance of particle deposition and the subsequent differentiation mechanisms that may occur within large plutonic reservoirs. In magma oceans, the physical context is quite different with potentially more vigorous convection related to much thicker reservoirs. Recent numerical models from [START_REF] Patočka | Settling of inertial particles in turbulent Rayleigh-Bénard convection[END_REF] constrained the settling of inertial particles in turbulent Rayleigh-Bénard convection without considering the re-entrainement process.

In particular, [START_REF] Patočka | Settling of inertial particles in turbulent Rayleigh-Bénard convection[END_REF] identified different settling regimes and derived characteristic residence timescales. For their range of particule sizes and densities, their results show that the settling should be rapid, much shorter than the typical time scales for the solidification of a whole-mantle terrestrial magma ocean supporting the idea of a fully fractional crystallisation. However, using analog models, [START_REF] Solomatov | Entrainment from a bed of particles by thermal convection[END_REF] introduced the conditions for particle re-entrainment after sedimentation at the bottom of a convecting reservoir.

For sufficiently high Rayleigh numbers, their models illustrate that the particles might be re-entrained by the viscous stress produced by thermal plumes emphasising that entrainement of crystals cumulated within dunes was likely at the bottom magma oceans. Such a re-entrainment process could hence favour equilibrium crystallisation. More recent laboratory experiments from [START_REF] Lavorel | Sedimentation of particles in a vigorously convecting fluid[END_REF] confirmed the founding from [START_REF] Solomatov | Entrainment from a bed of particles by thermal convection[END_REF] and characterized the influence of the density ratio between fluid and crystals and the temperature difference driving thermal convection.

Several parameters likely influence the settling dynamics: (1) the density difference between the settling crystal and the surrounding convecting magma ocean, (2) the crystal size that is a strong function of the magma ocean cooling rate and (3) the viscosity of the magma ocean which value may cover several orders of magnitude depending on the pressure/temperature conditions within the early magma ocean. Recently, [START_REF] Caracas | Melt-crystal density crossover in a deep magma ocean[END_REF] proposed that the accumulation of crystals might occur near the depth of neutral buoyancy between crystals and the coexisting melt. They localised this neutral buoyancy depth for different partition coefficients iron partitioning and solid fractions. They show that the crystals may cumulate and form a mushy layer inside the magma ocean separating a shallow MO and a basal MO. Here, we monitor the stability of a thin mushy layer cumulated within a magma ocean. We consider grain sizes, the density differences between the crystals and the surrounding magma ocean and viscosities of the magma ocean derived from petrological experiments. In particular, we aim at determining if the crystal layer will be sedimented at the bottom of the reservoir or will be efficiently mixed in the convective reservoir.

The hypothesis of a fully molten mantle is difficult to reconcile with geochemical data and the need for hidden reservoirs [START_REF] Boyet | 142 Nd Evidence for Early (¿4.53 Ga) Global Differentiation of the Silicate Earth[END_REF]. However, the cooling and solidification processes occurring within large partially molten silicate reservoirs may lead to important chemical differentiation events and to chemical segregation. The depth at which solidification initiates is governed by the crossing between the early geotherm and the melting curves. Crystal settling could then occur towards neutral buoyancy depths where crystals are paradoxically more easily entrained by the flow [START_REF] Caracas | Melt-crystal density crossover in a deep magma ocean[END_REF]. Both chemical composition and temperature contribute to the buoyancy of a crystal layer within a convecting reservoir. The chemical differentiation occurring after a large melting event is strongly related to the ability of the crystal fraction to sink or float within the convective reservoir. This ability is itself governed by the propension of the crystallised minerals to integrate in their structure heavy elements such as Fe [START_REF] Nomura | Spin crossover and iron-rich silicate melt in the Earth's deep mantle[END_REF][START_REF] Andrault | Solid-liquid iron partitioning in Earth's deep mantle[END_REF]. Depending on the crystal buoyancy relative to the molten magma, different scenarios may emerge with a bottom/up solidification [START_REF] Monteux | On the cooling of a deep terrestrial magma ocean[END_REF] or a mid-mantle solidification making a basal magma ocean possible [START_REF] Labrosse | A crystallizing dense magma ocean at the base of the Earth's mantle[END_REF]. In our study, we will discuss the implications of our model on the mechanism of magma ocean solidification after a major Moon forming impact on an Earth-like planet.

In the following study, we have developed a numerical model to monitor the cooling, segregation and chemical evolution of the early Earth mantle. Our study is organised as follows. In section 2, we present our physical and numerical models. In section 3, we present the main numerical results from our systematic parametric study emphasising the influence of the melt viscosity, the crystal size and the density difference between crystals and melt. Section 4 details a criterion for the stability of a crystallising layer. Finally section 5 describes the evolution of mushy material in the context of a progressive magma ocean crystallisation.

We discuss our results and conclude in section 6.

Physical and numerical models

Governing equations

We monitor the stability of a horizontal mushy layer composed of solid crystals within an initially convecting fully molten magma reservoir. For this purpose, we consider the Euler-Euler model which is an accurate dispersed multiphase flow model describing both the crystal and the liquid phases. We distinguish two non-miscible phases in our models: the continuous phase (with c subscript) that represents the liquid magma ocean and the dispersed phase The dynamical and thermal evolution of the mushy layer and ambient liquid are governed by the Navier-Stokes and heat transfer equations. We initially solve the equations for the continuous phase heat in 2D rectangular geometry using the Rayleigh-Bénard system. Once it reaches a steady state, we then solve the dynamical evolution of the mushy layer within the convective reservoir. The Euler-Euler flow model defines one set of Navier-Stokes equations for each phase:

• The mass conservation:

∇.(φ c u c + φ d u d ) = 0, (1)
∇.(φ d u d ) = 0 (2)
with u d and u c represent of the dispersed phase and the continuous phases fluid velocity vector respectively (m.s -1 ).

• The momentum conservation: • The heat conservation:

ρ c ( ∂u c ∂t + u c .∇u c ) = ∇.[-P I + η c (∇u c + (∇u c ) T )] + ρ c g + F m,c φ c (3) ρ d ( ∂u d ∂t + u d .∇u d ) = ∇.[-P I + η d (∇u d + (∇u d ) T )] + ρ d g + F m,d φ d (
ρ c C p [ ∂T ∂t + u.∇T ] + ∇.[-k∇T ] = 0 (5)
with C p the heat capacity at constant pressure (J.kg -1 .K -1 ), T the fluid temperature (K), u the fluid velocity vector (m.s -1 ), k the thermal conductivity (W.m -1 .K -1 ). No heat source is considered in our models (e.g. radiogenic, tidal or latent heat). For simplicity, we consider that both the continuous and dispersed phases have the same heat capacities and thermal conductivities.

For each fluid parcel, the viscosity ranges between the viscosity of the fully molten phase η c that is imposed for each simulation (see Tab. 1 for values) and the viscosity of the solid dispersed phase. In between the mixture viscosity η depends on the dispersed solid phase volume fraction φ d . For the mixture, we consider a Krieger-type dynamic viscosity:

η = η c (1 -M IN ( φ d φ d,max , 0.999)) -2.5φ d,max (6)
where φ d,max is the maximum packing concentration. φ d,max depends on the geometry of the crystals within the system and on the shear rate. For monodispersed systems with spheres and low shear rates, φ d,max ranges between 0.5 and 0.7 [START_REF] Dörr | A discrete model for the apparent viscosity of polydisperse suspensions including maximum packing fraction[END_REF]. We consider in our models that φ d,max is fixed to 0.6 meaning that φ c can locally range between 40 and 100% while φ d can locally range between 0 and 60%.

To compute the drag force F m , we consider the Schiller-Naumann drag model [START_REF] Schiller | A drag coefficient correlation[END_REF][START_REF] Ibrahim | Coupled Flow Modelling in Geotechnical and Ground Engineering: An Overview[END_REF] that is relevant for dispersed rigid spheres and for the range of Reynolds number (comparing inertial and viscous forces) values in our study (lower than 10 4 ). F m is defined as:

F m,c = -F m,d = β(u d -u c ) (7)
with β the momentum transfer coefficient (in kg m -3 s -1 ) between fluid and solid particles [START_REF] Ibrahim | Coupled Flow Modelling in Geotechnical and Ground Engineering: An Overview[END_REF])

β = 3 4 ρ c φ d C D (u d -u c ) D d ( 8
)
D d is the particle diameter (m) and C D is the drag coefficient for a single particle. For Re < 1000:

C D = 24(1 + 0.15Re 0.687 ) Re (9)
For Re > 1000, C D = 0.44 [START_REF] Wen | Mechanics of Fluidization[END_REF].

The density of the magma ocean is supposed to vary linearly with temperature in the buoyancy term of the momentum equation following the Boussinesq approximation. The density calculation accounts for the contrast related to the thermal effects and the phase influence. It is calculated as:

ρ c = ρ c,0 (1 -α(T -T 0 )) ( 10
)
where T 0 is the reference temperature (in K) and α is the thermal expansion coefficient (K -1 ). The initial reference temperature of the reservoir is highly influenced by the full thermal history of the early earth as well as its depth within the Earth mantle. Here we fix T 0 = 4000 K. The particle density is calculated as:

ρ d = ρ c + ∆ρ (11)
with ∆ρ the intrinsic density difference between the molten silicates and the solid crystals.

Non dimensional parameters

Before considering the stability of a crystal layer within the molten reservoir, we first solve only the equations relative to thermal convection within the reservoir. In our geophysical investigation, the main parameter is the magma ocean viscosity which ranges widely in value. Here we use viscosities for the continuous phase (η c,0 ) ranging between 0.1 and 10 Pa.s [START_REF] Karki | Viscosity of MgSiO 3 Liquid at Earthś Mantle Conditions: Implications for an Early Magma Ocean[END_REF][START_REF] Xie | Direct Viscosity Measurement of Peridotite Melt to Lower-Mantle Conditions: A Further Support for a Fractional Magma-Ocean Solidification at the Top of the Lower Mantle[END_REF]. The convection flow induced within the magma ocean is characterised by the thermal Rayleigh number (Ra) describing the vigor of convection as a consequence of basal heating:

Ra = αρ c g∆T H 3 η c κ (12)
with κ the heat diffusivity of the convective fluid (see Tab. 1 for value).

This leads to Ra values ranging between 1.8 × 10 6 and 1.8 × 10 8 . From our purely thermally convective models, we monitor the surface heat flux F which is proportional to the average value of the temperature gradient at the surface of the reservoir. We obtain F as a function of time for the range of Ra values considered in our study. We stop our thermally convective simulations once the surface heat flux has converged to a constant value and use this value to calculate the corresponding Nusselt number N u:

N u = HF k∆T (13)
Fig. 1 shows the N u-Ra relation obtained from our models and compares our numerical results with the scalings obtained from [START_REF] Turcotte | Geodynamics, Applications of continuum physics to geological problems[END_REF] and [START_REF] Wolstencroft | Nusselt-Rayleigh number scaling for spherical shell Earth mantle simulation up to a Rayleigh num-ber of 10 9[END_REF]. Our results show an agreement with these 2 scalings and we obtain N u = 0.14Ra 0.34 . For each Ra value used in our models we obtain a temperature field that will be used as the initial temperature field for the Euler-Euler model in order to study the stability of a thin mushy layer with different characteristics (D d and ∆ρ).

Numerical model

We used the commercial software COMSOL Multiphysics (version 5.4) previously validated for two phase flow applications [START_REF] Qaddah | Dynamics and stability of an iron drop falling in a magma ocean[END_REF][START_REF] Qaddah | Thermal evolution of a metal drop falling in a less dense, more viscous fluid[END_REF]).

We consider a 2D domain with a thickness H and a width of 4H (See Fig. 2). As a first attempt to characterise the stability of a thin newly formed mushy layer within a convective magma ocean, we arbitrarily consider that the thickness of this mushy layer is h = H/10. The mush thickness is set constant in all our models. Within the mushy layer, the crystals are considered as spheres with a fixed diameter D d . In our models, this diameter ranges between 10 -4 and 10 -2 m. The initial volume fraction of crystal φ d,0 is also set to a constant value (φ d,0 = 0.1). For the dynamics, we use free-slip conditions at all the system boundaries for both the dispersed and continuous phases. For the temperature, we consider no flux conditions at lateral boundaries. We fix the temperature to 4100 K at the bottom boundary and to 4000 K at the top boundary (i.e.

∆T = 100K).

The COMSOL Multiphysics software employs the finite-element method on a grid mesh to solve the set of equations (Eqs. 1 to 5). In order to study the sensitivity analysis of mesh size on the main results of this paper, we run several models varying the size of the mesh (between 1.6 mm and 50 mm). In these resolution tests we consider the overturn of a mushy layer with h = H/10, φ d,0 = 0.5, η c = 10 (Ra = 1.8 × 10 6 ) and ∆ρ = 30 kg.m -3 . Fig. 3 shows the average value φ d,average along a horizontal profile located at the bottom of the reservoir at t = 600s when the mushy layer is fully sedimented. φ d,average

increases when decreasing the mesh size from 0.05 m to 0.01m. Then, below a mesh size value of 0.01 m, φ d,average converges to a constant value of 0.45. In order to save computational time, we consider a mesh size of 0.08 which leads to an error of 7.7% relative to the converged value obtained for the smallest mesh size. In addition, we decrease the size of the mesh by a factor 2 in the lowermost part of the reservoir where crystal-rich patches are likely to settle or within the whole reservoir for the most vigorously convective contexts. Therefore, the calculation domain is composed of 100 000 to 180 000 cells. We list the details of our models in Tab. 2. Each run presented in this study represents (for the Euler-Euler flow model only) 12 to 72 hours computation time on a bi-processor, 8-cores, 3.2-GHz workstation.

Reference case

In this subsection, we detail the case corresponding to simulation no. 47 in Tab. 2. For this reference case, we consider the case with Ra = 4.4 × 10 6 (corresponding to η c = 4 Pa.s), ∆ρ = 50 kg.m -3 (∆ρ/ρ c,0 = 1.7% ) and D d = 10 -3 m. Fig. 4 illustrates the thermal and dynamical time evolutions of the convective magma ocean and the mushy crystal layer. Because of negative buoyancy, the crystal layer rapidly sinks. Fig. 4 shows that an overturn rapidly occurs concentrating the crystals in the lower part of the reservoir during the first hundreds of seconds after the layer formation. The crystal fraction is first diluted in the lower reservoir then starts to organise in patches where the concentration is much higher than the initial concentration. This early stage limits the heat transfer from the bottom surface to the reservoir leading to a transient cooling of the magma.

To monitor the evolution of the average concentration C (particles/m 3 ) in crystals in the upper mid reservoir for the reference case, we define:

e = H H/2 4H 0 Cdxdz H 0 4H 0 Cdxdz (14)
Höink et al. (2005) refers to this parameter as the relative entrainment. Fig. 5 illustrates the time evolution of e for our reference case and shows that during 2000 s, the upper reservoir is depleted in crystals. After nearly 10 000 s, large hot plumes formed at the bottom of the reservoir (see Fig. 4) transport the crystals back in the upper reservoir. Hence the concentration in crystals in the upper reservoir increases. This second phases leads to an efficient mixing in the whole reservoir and the crystal concentration gets back and slightly oscillates around its initial value. Fig. 5 shows that the initial concentration in the upper reservoir is slightly larger than its theoretical value of 0.01. This is due to the mesh properties that encompasses a larger volume of crystals leading to an initial value of e = 0.0105.

Parametric study

In the following section, we monitor the influence of the density difference between the crystals and the magma, the crystal size, and the magma viscosity on the settling dynamics within a convecting magma ocean.

Influence of the density difference ∆ρ

During the solidification of the magma ocean, the density of the crystal phase potentially increases (or decreases) relative to the molten phase depending on its ability to integrate dense chemical elements in its structure. Within the fully molten mantle, the first crystal to form is a crystal of bridgmanite. This first crystal is Fe-poor, and becomes chemically neutrally buoyant close to the core mantle boundary [START_REF] Caracas | Melt-crystal density crossover in a deep magma ocean[END_REF]. Then, the relative Fe content increases in the melt potentially leading to the formation of a basal magma ocean overlaid by a thick layer of bridgmanite. The front of neutral buoyancy moves towards the shallow mantle as the solidification proceeds [START_REF] Solomatov | Fluid Dynamics of a Terrestrial Magma Ocean[END_REF][START_REF] Boukaré | Thermodynamics of the MgO-FeO-SiO 2 system up to 140 GPa: Application to the crystallization of Earth's magma ocean[END_REF][START_REF] Caracas | Melt-crystal density crossover in a deep magma ocean[END_REF]. This scenario potentially leads to a fractional crystallisation of the early mantle [START_REF] Patočka | Settling of inertial particles in turbulent Rayleigh-Bénard convection[END_REF]. However, when accounting for re-entrainment enhanced by vigorous convection, crystal dispersion might be an efficient process that should favor batch (or equilibrium) crystallisation of the early mantle [START_REF] Solomatov | Entrainment from a bed of particles by thermal convection[END_REF]. The melt-crystal density crossover in both these scenarios is discussed by [START_REF] Caracas | Melt-crystal density crossover in a deep magma ocean[END_REF]. Here, we monitor the influence of the chemical density difference on the settling dynamics considering relatively small chemical density differences (10-100 kg.m -3 ) which is relevant in the early stages of mantle solidification. The chemical density difference is fixed in each simulation and is superimposed to the thermal density difference which can evolve during the cooling. We note here that in our models, the density difference related to temperature differences (αρ∆T = 3 kg.m -3 ) is smaller that the density difference between the magma and the crystals.

In Fig. 6, we considered the reference case illustrated in Fig. 4 in which we increased the density difference from ∆ρ = 50 kg.m -3 to ∆ρ = 100 kg.m -3 .

As already illustrated in Fig. 4, Fig. 6 shows that because of negative buoyancy, the crystal layer rapidly sinks. After a rapid settling at the bottom of the reservoir, particles entrainment occurs from the dunes (10 000 s after the initiation of the settling). However, because of a larger density difference between the dispersed phase and the continuous phase, efficient mixing is prevented.

Particle entrainment nearly vanishes and the layer is fully sedimented at the bottom of the reservoir (28 000 s after the initiation of the settling). Again, the sedimented layer acts as an insulating layer by limiting the heat transfer from the bottom of the reservoir. As a consequence, the mean temperature is significantly lower due to bottom insulation when full sedimentation is achieved.

t overturn = 4π 2 η lay,0 (φ d,0 ∆ρ -αρ lay,0 ∆T )gδ ( 15
)
where η lay,0 is the initial viscosity of the crystal rich layer (see Eq.6) and ρ lay,0 is the initial density of the crystal rich layer (= φ d,0

ρ d + (1 -φ d,0 )ρ c ).
The overturn timescales decreases with increasing density differences as shown in Fig. 5.

Influence of the crystal size

Upon solidifying, crystals form following the nucleation and growth mechanisms [START_REF] Solomatov | Magma Oceans and Primordial Mantle Differentiation[END_REF]. The crystal size reached at equilibrium depends on the nature of the mineral and on the cooling rate within the convecting reservoir.

The grain size may range between 10 -4 m and 10 -2 m in typical magma ocean conditions [START_REF] Solomatov | Magma Oceans and Primordial Mantle Differentiation[END_REF]. Here, we monitor the influence of the crystal size in our models. To that aim, we compare the reference case (with

D d = 10 -3 m)
to cases with crystals larger (D d = 10 -2 m) and smaller (D d = 10 -4 m). Fig. 7 shows the time evolution of the crystal fraction in the upper part of the reservoir for the three particle diameters D d considered. Fig. 7 illustrates that, even for the wide range of crystal sizes considered (2 orders of magnitude), the overturn timescale is poorly sensitive to the particle diameter D d (factor 3 in difference).

This result is in agreement in agreement with Eq. 15.

However, Fig. 7 shows that the final behaviour of the reservoir in terms of sedimentation/suspension is dependent on the crystal size. Indeed, for the parameters chosen for the reference case ( η c = 4 Pa.s and ∆ρ = 50 kg.m -3 ), Fig. 7 shows that when the particle size is equal to 10 -2 m, the suspension regime is not reached after the overturn. For smaller particle diameters (D d = 10 -4 and D d = 10 -3 ), the overturn process is rapidly followed by a suspension regime where the particles are efficiently mixed within the convecting reservoir.

Influence of the viscosity

The viscosity of the molten phase (η c ) is a key parameter governing the cooling dynamics of the magma ocean and in particular the solidification timescales [START_REF] Monteux | On the cooling of a deep terrestrial magma ocean[END_REF]. According to [START_REF] Karki | Viscosity of MgSiO 3 Liquid at Earthś Mantle Conditions: Implications for an Early Magma Ocean[END_REF] Fig. 8 shows that the continuous phase viscosity governs the ability of the solid phase to mix or sediment at the bottom of the reservoir. The interesting feature arises from the non linearity of the process. Indeed Fig. 8 shows that for the lowest viscosity, the crystal layer rapidly mix with the convecting reservoir.

Then for 0.4 ≤ η c ≤ 1, the rapid overturn process is not followed by a mixing process and the final concentration within the upper part of the reservoir does not reach the typical value synonymous of efficient mixing. Finally, for η c > 1 Pa.s, the overturn process is followed by a vigorous mixing process similarly to the case with η c = 0.1 Pa.s.

Fig. 8 shows that the molten phase viscosity also governs the early overturn timescale (i.e. the time needed for the particle concentration in the upper part of the reservoir to reach the minimal value). This overturn time increases from nearly 40s (for η c = 0.1 Pa.s) up to 2000s (for η c = 10 Pa.s). We represent these overturn timescales in Fig. 9 and compare them with the scaling law from Eq. 15. Fig. 9 shows that the overturn timescales from our models is 5 to 40 times larger than the overturn timescale from Eq. 15. We attribute this discrepancy to two points: (1) Eq. 15 is related to an overturn timesecale while settles with a timescale that is larger than if the layer could maintain its initial structure.

A criterion for the stability of a crystallising layer

We have illustrated the separated influence of the density difference between crystals and melt, grain size and melt viscosity on the behaviour of a crystallising layer within a magma ocean. We now aim at anticipating the ability of the flow generated from the convection within magma ocean to initiate the motion and suspension of sedimented crystals and to lead to efficient mixing within the magma ocean. We introduce the Shields number Sh that is the ratio between the tangential stress from convective flow over the buoyancy stress relative to density difference between the crystals and the melt [START_REF] Solomatov | Entrainment from a bed of particles by thermal convection[END_REF]: ∆ρ crit is derived from Eq. 16 by calculating:

Sh = η c αgF C p 1/2 1 ∆ρgD d ( 16
∆ρ crit = η c αgF C p 1/2 1 Sh c gD d ( 18
)
and considering Sh c = 0.1 which is a conservative value regarding a magma ocean viscosity of 0.1 Pa.s or less [START_REF] Xie | Direct Viscosity Measurement of Peridotite Melt to Lower-Mantle Conditions: A Further Support for a Fractional Magma-Ocean Solidification at the Top of the Lower Mantle[END_REF]. At a given depth (z), a ∆ρ M O larger than the critical value ∆ρ crit would yield an efficient Bg sedimentation to the bottom of the MO. By extension, we also consider that Bg would efficiently segregate into a solid layer above a basal magma ocean if ∆ρ M O /∆ρ crit is smaller than -1. Calculating ∆ρ M O requires a number of assumptions:

1. We use pyrolite as a model of reference for the mantle composition. Pyrolite is obtained by adding some crust component to the present-day depleted (peridotitic) upper mantle [START_REF] Ringwood | Pyrolite and the Chondritic Earth Model[END_REF]. To implement flexibility on the possible mantle compositions, we add or subtract a component of bridgmanite (Bg) to pyrolite; it yields less or more mafic mantle material, respectively. Later, we will argue that the composition of the primitive upper mantle (PUM) should be close to pyrolite, as this is widely accepted. 3. The Bg composition is calculated by equilibrating it with the melt in the mushy MO, using partitioning coefficients K(Al) Bg/M elt = 0.5 and K(F e) Bg/M elt = 0.1 or 0.25. The value of K(Al) Bg/M elt does not impact our calculations significantly. K(F e) Bg/M elt remains controversial as it is reported between ∼ 0.07 to ∼ 0.5 [START_REF] Andrault | Solid-liquid iron partitioning in Earth's deep mantle[END_REF][START_REF] Nomura | Spin crossover and iron-rich silicate melt in the Earth's deep mantle[END_REF]. The former value could be underestimated as it was measured in a sample presenting extreme chemical segregation in a temperature gradient, which artificially promotes the Fe-depletion of Bg at the highest temperature. In the latter study, the experimental conditions corresponded to solidus-like temperatures, at which Fe is more compatible in Bg. Relevant temperatures for MO solidification evolve from liquidus-like to intermediate between liquidus and solidus. Therefore, we perform our calculations for K(F e) Bg/M elt values of 0.1 or 0.25.

The composition of the bulk silicate

4. As a reference model for the mushy MO, we use Bg grain diameter (D Bg ) of 0.1 to 10 mm, viscosity (η) of 0.1 Pa.s [START_REF] Xie | Formation of bridgmaniteenriched layer at the top lower-mantle during magma ocean solidification[END_REF], intrinsic conductivity (k) of 5 W/m/K [START_REF] Hofmeister | Mantle Values of Thermal Conductivity and the Geotherm from Phonon Lifetimes[END_REF] and heat capacity (Cp) of 1800 J/Kg/K [START_REF] Asimow | Melts under extreme conditions from shock experiments[END_REF]. We consider that the melt has the same composition at all depths because of efficient mixing in the turbulent MO.

5. We consider that the term ∆T used in F calculation from Eq. 18 is the temperature difference between the Earth's surface and the bottom of the MO, after correcting for the adiabatic increase of temperature with mantle depth. We first fix the surface temperature (T surf ) to 500 K [START_REF] Lebrun | Thermal evolution of an early magma ocean in interaction with the atmosphere[END_REF][START_REF] Monteux | A mushy Earth's mantle for more than 500 Myr after the magma ocean solidification[END_REF], however, we will investigate later T surf of 1000 and 1500 K. We consider that the bottom of the MO (z M O ) is the depth where the mantle solidification is achieved for a mixture of 60% of Bg and 40% of melt. During the mantle solidification, z M O evolves from the CMB depth towards the surface. We assume that the temperature at a depth z M O follows the solidus temperature profile of pyrolite (T sol (z M O )) (Pierru et al., subm.). We note that using the solidus profile yields a conservative ∆T , compared to the true temperature that should be between solidus and liquidus. Then, we calculate ∆T (z M O ) of a MO extending from a depth z M O to the Earth surface using the following equation:

∆T (z M O ) = T sol (z M O ) -T surf -∆T ad (z M O ) (19)
where ∆T ad (z M O ) is the temperature difference between T sol (z M O ) and the corresponding potential surface temperature extrapotated from T sol (z M O )

to the surface along the MO adiabatic temperature profile. The adiabatic gradient in the MO is calculated using the equations of state available for different silicate melt compositions [START_REF] Asimow | Melts under extreme conditions from shock experiments[END_REF]. ∆T (z M O ) is found to range between ∼ 1300 K and ∼ 2000 K for MO extending from 660 km (24 GPa) to 2800 km (135 GPa), respectively (see Fig. 11).

6. We calculate the density of Bg (ρ Bg ) along the mantle solidus, using the equations of state (EoS) of MgSiO 3 , FeSiO 3 and Al 2 O 3 end-members [START_REF] Matas | On the bulk composition of the lower mantle: predictions and limitations from generalized inversion of radial seismic profiles[END_REF]. The decomposition in end-members assumes Bg of (Mg [START_REF] Asimow | Melts under extreme conditions from shock experiments[END_REF]. These calculations yield the Bgto-melt density contrast (∆ρ M O = ρ Bg -ρ M elt , Fig. 12) along the mantle solidus: ∆ρ decreases with increasing pressure whatever is the composition of the mushy mantle, because the melt is more compressible than Bg.

7. We do not attempt to model the latest stages of MO crystallization, when the bottom of the MO is shallower 660 km depth, because the phase relations rely on the melting diagrams at upper mantle P-T conditions.

Conditions for segregation of Bg in the mushy magma ocean

Effect of melt composition

We now calculate the ratio of density differences χ ρ = ∆ρ M O /∆ρ crit at all lower mantle depths using the parameters defined above. First, we fix the grain size of Bg to 1 mm, K(F e) Bg/M elt to 0.1 or 0.25, and vary the MO composition from pyrolite -40% of Bg to pyrolite +40% of Bg (i.e. 70% of pyrolite and 30% of Bg). χ ρ is found to be positive or negative (Fig. 13), depending on the solidmelt density difference (Fig. 12). A ∆ρ M O value larger than ∆ρ crit (i.e. χ ρ > 1)

indicates an efficient Bg segregation to the bottom of the MO. Our calculations show that Bg segregation happens in shallow MOs for all chemical compositions investigated. For pyrolite, this is in agreement with a previous work dedicated to the determination of the MO viscosity [START_REF] Xie | Formation of bridgmaniteenriched layer at the top lower-mantle during magma ocean solidification[END_REF]. The higher is the Bg content in the MO, the earlier (the deeper) Bg segregation occurs during the MO solidification; this is because ∆ρ M O is larger when the bulk Fe-content is lower.

As expected, lower K(F e) Bg/M elt favors a denser melt, and thus lower ∆ρ M O .

Still, our calculations show a moderate impact on the results of K(F e) Bg/M elt variations from 0.1 to 0.25. Finally, efficient Bg segregation could occur above a dense MO for a MO composition highly depleted in Bg; this is suggested by χ ρ < -1 for pyrolite -40%Bg at pressures larger than ∼ 120 GPa, for a Bg grain size of 1 mm (Fig. 13).

Effect of grain size

We now investigate the role of the Bg grain size, for five different mantle compositions (Fig. 14). We first focus on a melt of pyrolite composition. χ ρ is found positive at most of the mantle depths, because ∆ρ M O is also positive (Fig. 12), except in the lowermost mantle for K(F e) Bg/M elt = 0.1. Still, χ ρ can only be larger than 1 for D Bg larger than 0.5 mm; Bg-grains of e.g. 2 mm segregate when pressure (depth) is larger than 70 GPa (1650 km). Grain size in a crystallizing MO remains a complicated issue [START_REF] Solomatov | Magma Oceans and Primordial Mantle Differentiation[END_REF]: (i) as long as the MO is fully liquid close to the planetary surface, grains are often remolten within the large scale convection movements and, therefore, grain size is controlled by nucleation processes; (ii) When the crystallization zone reaches the surface, it enables some grains to grow indefinitely. According to the shape of mantle solidus and liquidus, this change of mechanism could happen when the bottom of the MO reaches a depth of 1000 km ( 40 GPa) (see Xie et al.

(2020)). In this article, we aim at modelling the chemical segregation in the first stages of the MO solidification, when grain size is controlled by nucleation processes. D Bg in these conditions depends largely on the cooling rate of the MO, with estimated values of 0.1 to 1 mm [START_REF] Solomatov | Magma Oceans and Primordial Mantle Differentiation[END_REF]. For 1 mm grain size, χ ρ is larger than unity at pressures (depths) larger than 45GPa (1100 km) (Fig. 13). As a result, it appears unlikely that Bg grains could segregate in a crystallizing MO of pyrolite composition for a depth larger than ∼ 1100 km. Alternatively, a BSE composition depleted in Bg yields a denser melt in the mushy MO, as discussed in several works [START_REF] Boukaré | Thermodynamics of the MgO-FeO-SiO 2 system up to 140 GPa: Application to the crystallization of Earth's magma ocean[END_REF]Ballmer et al., 2017a;[START_REF] Caracas | Melt-crystal density crossover in a deep magma ocean[END_REF]. Our calculations of χ ρ for melt compositions of pyrolite -20% and -40% of Bg show that Bg is indeed less dense than the melts in a large fraction of the lower mantle (Fig. 12). Bg remains always denser at low pressures, due to its very high bulk modulus. For Bg grain size of 1 mm, values of χ ρ are found lower than -1 at pressures (depths) larger than 120 GPa (2600 km) for a melt composition of pyrolite -40% of Bg (Fig. 14 e). The pressure (depth) at which Bg could segregate upwards decreases with increasing the Bg grain size. From a geochemical point of view, an early segregation close to the CMB of a Bg-depleted melt requires a BSE composition depleted in Bg, if the shallow MO were to be pyrolitic.

Effect of temperature

Another important parameter in ∆ρ crit is ∆T (Eq. 19), which is intrinsically related to the rate of MO cooling. During the mushy stage of a MO solidification, the surface temperature is ≈ 500 K, which yields to final MO solidification in 100-200 thousands years after the Moon forming impact [START_REF] Lebrun | Thermal evolution of an early magma ocean in interaction with the atmosphere[END_REF][START_REF] Monteux | A mushy Earth's mantle for more than 500 Myr after the magma ocean solidification[END_REF]. Presence of a blanketing atmosphere atop the MO would induce a larger surface temperature and a lower rate of cooling [START_REF] Sleep | Terrestrial aftermath of the Moon-forming impact[END_REF][START_REF] Salvador | The relative influence of H 2 O and CO 2 on the primitive surface conditions and evolution of rocky planets[END_REF]. We note that the formation of a blanketing atmosphere is unlikely to affect the first stages of the MO solidification for different reasons: (i) H 2 O (a most efficient blanketing compound) is efficiently vaporized to the atmosphere only after a high level of MO solidification [START_REF] Elkins-Tanton | Linked magma ocean solidification and atmospheric growth for Earth and Mars[END_REF][START_REF] Salvador | The relative influence of H 2 O and CO 2 on the primitive surface conditions and evolution of rocky planets[END_REF]. (ii) a surface temperature larger than 1500K yields an atmosphere dominated by rock-vapor, which is inefficient to stop the surface cooling by thermal radiation [START_REF] Sleep | Terrestrial aftermath of the Moon-forming impact[END_REF]. For these reasons, we also perform calculations of ∆ρ M O for surface temperatures of 500K, 1000K and 1500 K (Fig. 15).

Calculations show that larger surface temperatures yield more extreme positive and negative values of χ ρ and therefore easier Bg segregation. This is because larger T surf corresponds to smaller temperature difference ∆T (Eq. 19) and thus a less turbulent MO. Still, the impact remains modest; for a pyrolite mantle, variation of T surf from 500 to 1500 K yields Bg segregation at pressures (depths) lower than 45 GPa (1100 km) and 67 GPa (1600 km), respectively. is very likely that it would be enriched in SiO 2 , compared to pyrolite [START_REF] Allègre | The chemical composition of the Earth[END_REF][START_REF] Javoy | The integral enstatite chondrite model of the Earth[END_REF][START_REF] Mcdonough | The composition of the Earth[END_REF]. Such enrichment raised long-term discussions about how much the core could contain some Si and how much the deep mantle could be enriched in SiO 2 , compared to the present day upper mantle. A high MgSiO 3 /MgO ratio in the lower mantle is compatible with results from inversions of seismic profiles based on elastic parameters of minerals [START_REF] Samuel | Heterogeneous lowermost mantle: Compositional constraints and seismological observables[END_REF][START_REF] Matas | On the bulk composition of the lower mantle: predictions and limitations from generalized inversion of radial seismic profiles[END_REF][START_REF] Murakami | A perovskitic lower mantle inferred from high-pressure, high-temperature sound velocity data[END_REF].

Our model shows that increasing the Bg content in BSE favors the segregation of Bg in a deep MO (Fig. 14 a andb). Bg grains of 1mm diameter segregate efficiently (i.e. χ ρ > 1) at pressures (depths) lower than 45GPa (1100 km), 55 GPa (1350 km) and 65 GPa (1550 km) for pyrolite, pyrolite with addition of +20%Bg and +40%Bg, respectively. Still, Bg grains would segregate at the early stages of the solidification if the grain size would be 5 mm or more.

Such big grains are not expected in a very thick MO [START_REF] Solomatov | Magma Oceans and Primordial Mantle Differentiation[END_REF]. Bg segregation becomes possible when the MO is shallower and the greater the Bg enrichment, the earlier (the deeper) the Bg-segregation. For a grain size of 1mm, a Bg layer could form at MO depths of 1100 to 1550 km, depending on the BSE composition. Such thin MO could also favor grain growth, because the MO is partially crystallized up to the Earth's surface [START_REF] Solomatov | Magma Oceans and Primordial Mantle Differentiation[END_REF][START_REF] Xie | Formation of bridgmaniteenriched layer at the top lower-mantle during magma ocean solidification[END_REF]. Such Bg-layer could yield silica-enriched domains (Ballmer et al., 2017b). Unfortunately, our models are not adapted to address the thickness of layer, because of too many adjustable parameters. Segregation of some Bg in the lower mantle would drive the MO toward a pyrolite mantle, for a BSE is enriched in Bg. Decreasing the Bg-content makes the segregation of Bg grains less efficient (Fig. 13).

Solidification of a pyrolite-type bulk silicate Earth

A pyrolite-like BSE would solidify with no Bg segregation from the melt, as long as the Bg grain size always remains lower than 0.5 mm (Fig. 14 c). Bg segregation could still happen in the uppermost lower mantle if the Bg grain size gets larger. We note, however, that the grain growth of olivine and other (iii) We calculate the composition of the coexisting melt and Bg grains, fixing K(F e) Bg/M elt = 0.25. Anyway, the value of K(F e) Bg/M elt has no significant impact on this calculation. (iv) Starting from the CMB, we cool the mantle and calculate the composition of a mixture of 60% of Bg grains and 40% of melt which solidifies at the bottom of the MO. The chemical content of this bottom layer is then subtracted to the composition of the overlying mushy MO, before a next step of the solidification is performed within a thinner MO.

This calculation shows significant changes of the SiO 2 , FeO and CaO contents with depth in the solidified mantle (Fig. 16). They occur due to Bg presenting higher SiO 2 and lower FeO and CaO contents, compared to the melt.

The shallow MO presents a chemical composition close to pyrolite, when the BSE composition is adjusted to 80% of pyrolite and 20% of Bg. The SiO 2 content in the deep mantle is found to decrease from ∼46 to 39%, which translate into a decrease from ∼90 to 70% of the Bg-content, with decreasing depth from 2900 to 660 km, respectively. The Fe-content also increases with decreasing depth, because of its incompatible character; this effect is found even larger when using K(F e) Bg/M elt = 0.1, instead of K(F e) Bg/M elt = 0.25. Still, the overlying mushy MO does not become denser than the solidified layers, because Bg is denser than the melt at all mantle depth for Bg-enriched BSE composition (Fig. 12).

Such a chemical layering would only survive as such in a mantle that would be stagnant after its solidification. Instead, one could expect important vertical mixing in mushy mantle presenting a relatively low viscosity and with a major heat flux from the core to the surface [START_REF] Monteux | On the cooling of a deep terrestrial magma ocean[END_REF]. Another effect that can moderate the chemical stratification is a simultaneous mantle solidification in a large range of mantle depths. Upon cooling, the temperature profile becomes almost parallel to the solidus [START_REF] Miller | The equation of state of a molten komatiite: 2. Application to komatiite petrogenesis and the Hadean Mantle[END_REF][START_REF] Stixrude | Thermodynamics of silicate liquids in the deep Earth[END_REF][START_REF] Monteux | On the cooling of a deep terrestrial magma ocean[END_REF], which could yield to a bulk solidification.

Still, these effects may not suppress entirely the chemical gradients, in particular between a shallow MO and the lower mantle. The progressive enrichment of Bg with mantle depth is in good agreement with the mineralogical models [START_REF] Samuel | Heterogeneous lowermost mantle: Compositional constraints and seismological observables[END_REF][START_REF] Matas | On the bulk composition of the lower mantle: predictions and limitations from generalized inversion of radial seismic profiles[END_REF][START_REF] Murakami | A perovskitic lower mantle inferred from high-pressure, high-temperature sound velocity data[END_REF], as well as with the difference between shallow the mantle composition today (pyrolite) and a "chondritic-type" BSE (see above).

We note that the scenario of Bg-segregation in the lowermost mantle does not help explaining the LLSVPs atop the CMB. Several origins were proposed for the LLSVP, in particular some primitive mantle that did not melted upon the Moon forming impact. Primitive mantle compositions imply a relatively low proportion of ferropericlase, which is compatible with the comparison between the LLSVP seismic signature and the elastic parameters of the minerals [START_REF] Vilella | Constraints on the composition and temperature of LLSVPs from seismic properties of lower mantle minerals[END_REF]. Still, LLSVPs represent a quite small fraction of the mantle volume, typically 1.5 to 2.5 % [START_REF] Burke | Plume Generation Zones at the margins of Large Low Shear Velocity Provinces on the core mantle boundary[END_REF]. Therefore, our models would not be much different if some primitive mantle would have survived to melting at the Moon Forming impact.

Conditions for the formation of a primary basal magma ocean

The segregation of large Bg grains to the bottom of a pyrolite mantle would yield a mafic Fe-enriched overlying MO. Following this scenario, it has been proposed that the melt could become denser than the underlying mantle and induce early mantle overturns [START_REF] Boukaré | Thermodynamics of the MgO-FeO-SiO 2 system up to 140 GPa: Application to the crystallization of Earth's magma ocean[END_REF]Ballmer et al., 2017a;[START_REF] Caracas | Melt-crystal density crossover in a deep magma ocean[END_REF]. Such a reversal could result in a secondary MO with a composition typical of a shallow MO. Our results do not strongly support this scenario, because the lower the Bg-content in the mushy MO, the less efficient the Bg segregation (Fig. 14). Therefore, bulk solidification should rapidly be favored, after some Bg-segregation may have occurred in a thick MO. Mantle overturns may still occur due to the crystallization processes in the upper mantle; this should be investigated in further works based on analyses of the melting diagrams in the upper mantle.

Formation of a basal MO during the solidification process (i.e. a primary basal-MO) would require the efficient segregation of Bg in an overlying solid layer [START_REF] Labrosse | A crystallizing dense magma ocean at the base of the Earth's mantle[END_REF]. This conditions translates into χ ρ < -1 in our model. Such condition is reached for mafic BSE composition, which present a lower Bg content, compared to pyrolite (Fig. 14 c andd). Still, our calculations show that Bg segregation upwards is only expected (i) in the lowermost part of a MO (ii) for a BSE composition depleted of more than -20% of Bg, compared to pyrolite and (iii) for Bg grain size larger than 0.5 mm. We note that low K(F e) Bg/M elt value at these conditions help making the melt denser than Bg (Fig. 13, [START_REF] Nomura | Spin crossover and iron-rich silicate melt in the Earth's deep mantle[END_REF]). For a MO composed of pyrolite -40%Bg, Bg grains of 1 mm diameter would aggregate efficiently to a Bg-rich layer that would be located pressure (depth) of 120 GPa (2600 km). For a MO composed of pyrolite -20%Bg, only Bg grains of 5 mm diameter could segregate efficiently upwards. The simultaneous formation of a basal MO and a solid layer of Bg may maintain mafic the composition of the overlying MO. No further Bg segregation should be expected above ∼2500 km (i.e. 100-110 GPa) in a mafic MO, for a reasonable Bg grain size. In case the MO would be further depleted in Bg due to the mechanism of "at equilibrium fractionation" described above (Fig. 16), it would make it even more mafic. At this point, it is difficult to imagine which geological process could lead to the disappearance of such a major volume of early mafic MO. The path to pyrolite-type upper mantle seems difficult to achieve. It makes the formation of a primary basal MO simultaneous to MO solidification relatively unlikely.

Conclusion

The solidification within a deep convecting magma ocean (MO) is a key process that significantly shaped the Earth's interior (Ballmer et al., 2017a) and likely the early mantle evolution for several planetary contexts : Mars (Elkins- et al., 2005), Mercury [START_REF] Mouser | Experimental Investigation of Mercury's Magma Ocean Viscosity: Implications for the Formation of Mercury's Cumulate Mantle, Its Subsequent Dynamic Evolution, and Crustal Petrogenesis[END_REF] or the Moon [START_REF] Elkins-Tanton | The lunar magma ocean: Reconciling the solidification process with lunar petrology and geochronology[END_REF][START_REF] Dygert | A Low Viscosity Lunar Magma Ocean Forms a Stratified Anorthitic Flotation Crust With Mafic Poor and Rich Units[END_REF]. This process is responsible for the separation between compatible elements that prefer the solid phase (e.g., Mg, Si) and incompatible elements which prefer the liquid phase (e.g., Al, Na, Fe). The solidi-fication depth is governed by the geotherm and the melting curves. With density profiles of the melt and bridgmanite crystals crossing in the lower mantle, crystal settling could then occur towards neutral buoyancy depths where crystals are paradoxically more easily entrained by the flow. Both chemical composition and temperature contribute to the buoyancy of a crystal layer within a convecting reservoir. Hence different scenarios may emerge: Bottom/up solidification or mid-mantle solidification making a basal magma ocean possible.

Tanton
In our study we have shown that the ability of the crystal fraction to sediment or to disperse within the domain is the result of a combination of the influence of the crystal size, the density difference between the crystals and the magma and the magma viscosity. We have shown that this ability is a fonction of the ratio of the convection over the buoyancy stress (expressed as the Sh dimensionless number). From a relevant range of parameters we observe that the critical Sh c separating suspension and sedimentation regimes exhibits values lower than proposed by [START_REF] Solomatov | Entrainment from a bed of particles by thermal convection[END_REF]. A a consequence, our results strengthen the idea that suspension mechanism should be the dominant process during the magma ocean solidification compared to sedimentation and, as a consequence, equilibrium crystallisation should be favoured compared to fractional crystallisation.

In our models, we consider a particule size range (10 -4 -10 -2 m) and a density difference range between the magma and the particles (10 -100 kg.m -3 ), which are probably closer to the relevant MO conditions (see e.g. ∆ρ = -200 to 200 kg/m 3 in Fig. 12) than the range of parameters proposed by Patočka et al.

(2020) (0.5 × 10 -3 -10 -2 m and 10 -3000 kg.m -3 ). Moreover, their models neglect the re-entrainement process by eliminating the particles that reach the bottom of the reservoir. Hence, it is not surprising that their models favour a rapid crystal settling and a a fractional crystallisation. [START_REF] Patočka | Settling of inertial particles in turbulent Rayleigh-Bénard convection[END_REF] identified three settling regimes: "stone-like", bi-linear and "dust-like" regimes depending on the ratio between the crystal terminal velocity and the average vertical velocity of the flow. It appears that the models we run mostly fall out in the "dust-like regime" where the terminal velocity of a particle is smaller than the average velocity within the convective reservoir. In the "dust-like" regime the particles from the models developed by [START_REF] Patočka | Settling of inertial particles in turbulent Rayleigh-Bénard convection[END_REF] appear uniformly distributed in the reservoir in agreement with our models.

We have also considered in our models different end-members for the chemical composition of the early mantle adding or subtracting up to 40% of Bg to a pyrolite mantle reference. Exploring the early mantle composition within such a range leads to a wide range of density difference between the molten and the solid phases. Hence the chemical composition also influences the ability of the crystals to cumulate at the bottom or in the middle of the magma ocean. Our models emphasize that (1) Bg grains are unlikely to segregate in a mantle of pyrolite composition, except in a shallow MO when the grain size may increase to 1 mm and more, (2) Bg grains are likely to segregate at the bottom of a MO enriched in SiO 2 , compared to pyrolite and (3) the formation of a basal MO simultaneous to mantle solidification could be possible for mafic BSE compositions (difficult to reconcile with the Earth's geochemical models). We also note that a homogeneous bottom-up solidification of the MO is expected to produce a significant SiO 2 enrichment with increasing mantle depth in a primitive mantle (compatible with seismic observations) and a higher SiO 2 content in BSE compared to pyrolite (compatible with Earth's geochemical models).

This study is a supplementary step toward a more realistic approach that should account for several other major processes. In our models we have not consider the time evolution of the crystal size, of the crystal fraction and of the density of the crystals during their growth. If crystal settling is likely to occur close to mechanical boundaries (CMB or solidified wall) where the velocity is lower than in the middle of the magma ocean ( within "calm" conditions ), future models should explore parameters compatible with the center of the magma ocean (i.e. larger Rayleigh numbers). This would imply much higher Reynolds number and would require to consider the inertial forces and the turbulence equations in our model. The likeliness of a thick magma ocean and the presence of a large Moon much closer to the Earth than today should also lead us to consider te influence of the Coriolis forces [START_REF] Solomatov | Fluid Dynamics of a Terrestrial Magma Ocean[END_REF]. 9. We note we have tested the case with ηc = 0.1 Pa.s with two different meshes (a 6 and 8 mm mesh) and the difference between the two curves is not significative.

Figure 9: Overturn timescale as a function of the molten phase viscosity. Black symbols illustrate the results from our numerical models. The dashes lines show the overturn timescale obtained from Eq. 15 (with 1 to 100 toverturn). We note we have tested the case with ηc = 0.4 Pa.s with two different meshes (a 4 and 8 mm mesh) and the difference between the two timescales is not significative.

Figure 10: Shields number value (Sh) as a function of the molten continuous phase viscosity (ηc) (or Ra for the top x-axis). Square symbols correspond to the models run in this study (see Tab. 2 for values). Empty black symbols represent cases where the suspension regime is reached. Filled green symbols represent cases where the sedimentation regime is reached. The two regimes are separated by dashed lines (green for ηc ≤ 1, blue for ηc ≥ 1) illustrating the critical Shields number (Shc). Red symbols and dotted lines illustrate the values for the cases represented in Fig. 5, Fig. 7 and Fig. 8 respectively. We calculate the Sh values considering the combination of Eq. 12, Eq. 13 and the expression N u = 0.14Ra 0.34 obtained from Fig. 1 to compute the heat flow F within Eq. 16.

Figure 11: Temperature profiles corresponding to (grey thick line) the solidus of the pyrolite mantle (Pierru et al., subm.) and (colored lines) adiabatic gradients of MO for potential surface temperature from 1800 to 2500 K.

Figure 12: Solid-melt density difference (∆ρ M O = ρ Bg -ρ melt ) calculated using the P-V-T equation of state of Bg [START_REF] Matas | On the bulk composition of the lower mantle: predictions and limitations from generalized inversion of radial seismic profiles[END_REF] and MO [START_REF] Asimow | Melts under extreme conditions from shock experiments[END_REF], after decomposition in end-members. Calculations are performed for several MO compositions by adding or subtracting a fraction of Bg to a pyrolite-based mantle reference, and for Fe partition coefficient between Bg and MO (K(F e) Bg/M elt ) of 0.1 or 0.25. Positive value indicates lighter melt. (Eq. 19). Thick, intermediate and thin lines correspond to 500, 1000 and 1500 K, respectively. The calculations are also performed for various mantle compositions and K(F e) Bg/M elt values.

Figure 16: Chemical profiles obtained in the case of no Bg segregation from the melt. The chemical gradient comes from the solidification of ∼60% of Bg and 40% of the melt at the bottom of the MO. These profiles assume no mantle mixing once the mantle layers are eventually solidified. "BSE" and "s-MO" are bulk silicate Earth and a shallow MO extending up to 660 km depth, respectively. Here, the BSE composition is adjusted to 80% pyrolite and 20% of Bg. Table 2: Parameters for all performed simulations used in this study. We calculate the Sh value from Eq. 16 considering the combination of Eq. 12, Eq. 13 and the expression N u = 0.14Ra 0.34 obtained from Fig. 1

(

  with d subscript) that represents the crystal phase. Both phases behave as Newtonian and incompressible fluids. The local volume fraction of the fluid continuous phase φ c and the dispersed solid phase φ d are the local average volume of melt and crystals respectively (φ c = 1 -φ d ).

  4) with ρ c the melt density (kg.m -3 ), ρ d the solid phase density, η c the melt viscosity (Pa.s), η d the solid phase viscosity (Pa.s), t the time (s), P the fluid pressure (Pa), g the gravitational acceleration (m.s -2 ) and I the identity matrix. F m,c and F m,d are the drag forces (N/m 3 )

Fig. 5

 5 Fig. 5 compares the time evolution of the concentration of particles in the upper part of the reservoir for the two density differences illustrated in Fig. 4 and 6. During all the simulation, the concentration for the case with ∆ρ = 50 kg.m -3 remains larger than the concentration for the case with ∆ρ = 100 kg.m -3 . This feature is in agreement with Lavorel and Le Bars (2009) who monitored the relative number of particles in suspension as a function of the fluid/particle density difference during the first 1000 s after being uniformly distributed within the reservoir. Fig. 5 also shows that increasing the density difference decreases

  , η c for liquid anhydrous MgSiO 3 may decrease down to 0.048 Pa.s at mid-mantle conditions (P = 70 GPa and T = 4000 K along the magma ocean isentrope). For molten peridotite, Xie et al. (2021) obtained a viscosity ranging between 0.017 and 0.038 Pa.s. To maintain a laminar flow, (i.e. with low Reynolds number), we consider 0.1 ≤ η c ≤ 10 Pa.s.

Fig. 9

 9 Fig. 9 represents the time from which the concentration in the upper part of the reservoir re-increases after the overturn which is closer to a re-entrainement timescale and (2) the mixing dynamics tends to dilute the concentration in particles within the overturning layer during its sinking. Initially the particle concentration in the layer is φ d,0 . As illustrated in Figs. 4 and 6 (second line), the vigorous convection redistribute the particles within the reservoir even during the early overturn. As a consequence the density of the mushy layer relative to the ambient fluid decreases from its initial value. The mushy layer

  ) Large Sh values favor a suspension regime while small Sh values favor a sedimentation regime. The two regimes are separated by a critical Shields number value Sh c . From their analog experiments, Solomatov et al. (1993) have estimated a critical value Sh c = 0.1 -0.2. We have run a total of 76 numerical models (with values listed in Tab. 2) considering different sets of D d , ∆ρ and η c with value ranges in Tab. 1. Each model corresponds to a Sh value calculated considering the combination of Eq. 12, Eq. 13 and the expression N u = 0.14Ra 0.34 obtained from Fig. 1 to compute the heat flow F (See Tab. 2 for values). Monitoring the concentration of grains in the upper reservoir (e(t)) in each models enables to decipher if the models end up in the suspension or sedimentation regime. In Fig. 10, we illustrate the Sh values as a function of the molten phase viscosity η c . A case where the grain concentration in the upper part of the reservoir (e) never reaches back its initial value is considered in the sedimentation regime. Oppositely, a case where the grain concentration in the upper part of the reservoir gets back to its initial value is considered in the suspension regime. As shown in Fig. 8, some cases may end up in an intermediate regime where the final concentration is 0.005 < e < 0.01. We consider that these intermediate cases fall out in the sedimentation regime. Sedimentation and suspension regimes are separated by the critical Shields number value (Sh c ). According to Fig. 10, the Sh c value is a non-linear function of the continuous molten phase viscosity. For η c > 1, Sh c seems to be constant and independent of the η c value. However, when η c < 1, the Sh c value decreases with decreasing values of the viscosity following Sh c = 0.1η c 0.45 . As a consequence Sh c is not influenced by the increase of Rayleigh number up to Ra = 1.8 × 10 7 . For Ra values larger than 1.8 × 10 7 the increase of the kinetic energy induced by the vigorous convection in the system results in a decrease of Sh c following: Sh c = 184.Ra -0.45 (17) The non-linear dependence of Sh c on the viscosity explains why the crystal settling behaviour may vary back and forth between sedimentation and suspension regimes when the viscosity increases. Indeed Fig. 10 shows with the red dotted line the path followed for the models detailed in Fig. 8. Moreover, Fig. 10 shows that Sh c is smaller than the range of values proposed by Solomatov et al. (1993) (0.1 -0.2). This feature is particularly visible for the lowest range of viscosities used in our models. For a magma ocean viscosity value of 0.1 Pa.s, the Sh c value may decrease down to 0.04 which means that the crystal settling is difficult to envision in such a convecting environment especially regarding the typical viscosity values for a magma ocean (Karki and Stixrude, 2010; Xie et al., 2021). Decreasing the viscosity down to the lowest values proposed by Xie et al. (2021) may enhance this feature. 5. Fate of mushy material in the context of a progressive magma ocean crystallisation 5.1. Parameters controlling the solidification of a magma ocean on a Earth-like planet We now investigate the implication of our model on the mechanism of magma ocean (MO) solidification on Earth. More precisely, we investigate the possibility of Bridgmanite (Bg) segregation depending on its buoyancy (positive or negative). To that aim, we compare the density difference between the MO and the solid crystal ∆ρ M O calculated for different compositions and the critical density contrast ∆ρ crit above which crystal segregation is likely to occur.

  Earth (BSE) is still not well established, because a number of deep mantle reservoirs could have a composition different than pyrolite: The large low seismic velocity provinces (LLSVP) seem to present a different mineralogy (e.g. Vilella et al. (2021)),the deep lower mantle could be more bridgmanitic than pyrolite[START_REF] Murakami | A perovskitic lower mantle inferred from high-pressure, high-temperature sound velocity data[END_REF] and there could be relics above the core-mantle boundary (CMB) of a basal magma ocean[START_REF] Labrosse | A crystallizing dense magma ocean at the base of the Earth's mantle[END_REF] or of early mantle overturns (e.g.[START_REF] Boukaré | Thermodynamics of the MgO-FeO-SiO 2 system up to 140 GPa: Application to the crystallization of Earth's magma ocean[END_REF]). In our models, we give a large freedom of size of the possible Bg-rich or mafic reservoirs that could have segregated in the deep mantle during the MO solidification. The flexibility of BSE composition is achieved by adding or subtracting up to 40% of Bg to the pyrolite mantle reference. This is a much broader range of chemical variation than the uncertainty remaining today on the major elements concentration in BSE.

A

  BSE composition enriched in Bg yields a lighter melt on the first stages of the MO solidification, and thus a larger ∆ρ M O . It makes it easier the segregation of Bg at the bottom of the mushy MO. For a melt composition enriched with 20% or 40% of Bg, Bg grains of 1 mm could efficiently segregate at pressures (depths) of 55 GPa (1350 km) or 65 GPa (1550 km), respectively (Fig.14 a and b). From a geochemical point of view, an early segregation of Bg to the CMB requires a BSE composition enriched in Bg if the late shallow MO should be of pyrolite composition; we will discuss this matter later in the article.

5. 3 .

 3 Possible scenarios of MO solidification 5.3.1. Conditions for Bg segregation at the bottom of the MO Geochemical scenarios of the Earth's formation involve building blocks of chondritic composition. Uncertainties remain about the composition of the primitive "chondritic-type" mantle after the core-mantle segregation. Still, it

  minerals of the upper mantle may have a moderate impact on the grain size of Bg in the lower mantle. Modelling the mechanism of MO solidification in the upper mantle is beyond the scope of this article. Interestingly, absence of Bg segregation does not imply a chemically homogeneous mantle once it is entirely solidified. The origin of this effect is that the solidified mantle contains ∼60% of Bg and ∼40% of melt, while the overlying mushy MO contains a different fraction of the two phases. Here we assume the rate of mushy-mantle compaction (with extraction of the melt) significantly slower than the rate of MO solidification, following a previous work (Solomatov , 2015). To model the chemical stratification resulting from homogeneous solidification, we proceed as follows: (i) we choose a BSE composition close to pyrolite. Still, we keep its Bg content adjustable, to enable a final shallow MO (thickness of ∼660 km) composition as close as possible from pyrolite. (ii) At each time step of the progressive MO solidification, we estimate the fraction of Bg grains in the MO above the bottom of the MO (φ M O ) based on our previous model of progressive MO solidification (Fig. 3 of Monteux et al. (2016)).

Figure 1 :

 1 Figure1: Nusselt number as a function of the Rayleigh number during the phase the first phases where only the thermal convection is considered (no settling of particles). Results from our models are represented with black circles. For comparison, we plot the scaling laws proposed by[START_REF] Wolstencroft | Nusselt-Rayleigh number scaling for spherical shell Earth mantle simulation up to a Rayleigh num-ber of 10 9[END_REF] and[START_REF] Turcotte | Geodynamics, Applications of continuum physics to geological problems[END_REF] with blue and green dashed lines respectively. The red dashed line represents the result from a power-law fitting of our numerical results. Corresponding scaling laws are detailed in the legend of the figure.

Figure 2 :

 2 Figure 2: Schematic representation of the computational domain of our study. Blue material represents the continuous phase (molten silicates) and green material represent the initially mixed continuous and dispersed phases. A close-up view is illustrated in the right panel.

Figure 3 :

 3 Figure3: Dispersed fraction averaged along an horizontal profile at the bottom of the computational domain as a function of the mesh size. In these models, h = H/10, φ d,0 = 0.5, ηc = 10 Pa.s and ∆ρ = 30 kg.m -3 . The red dashed line represents the asymptotic value at which φ d,average is converging for high-resolutions (i.e. small mesh size) models.

Figure 4 :

 4 Figure 4: Time evolution (from top to bottom) of the dispersed volume concentration (left) and the temperature (right) for the reference case (with Ra = 4.4 × 10 6 (corresponding to ηc = 4 Pa.s), ∆ρ = 50 kg.m -3 and D d = 10 -3 m).

Figure 5 :

 5 Figure 5: Time evolution of the particles concentration in the upper part of the reservoir with Ra = 4.4 × 10 6 (corresponding to ηc = 4 Pa.s) and D d = 10 -3 m. The red line illustrates the case where ∆ρ = 50 kg.m -3 and the green line the case where ∆ρ = 100 kg.m -3 . Vertical dashed lines illustrate the time when the minimum value of e is obtained, which is referred in the text as the overturn time.

Figure 6 :

 6 Figure 6: Time evolution (from top to bottom) of the dispersed volume concentration (left) and the temperature (right) for the case with Ra = 4.4 × 10 6 (corresponding to ηc = 4 Pa.s), ∆ρ = 100 kg.m -3 and D d = 10 -3 m.

Figure 7 :

 7 Figure 7: Time evolution of the particles concentration in the upper part of the reservoir with Ra = 4.4 × 10 6 (corresponding to ηc = 4 Pa.s) and ∆ρ = 50 kg.m -3 . The red line illustrates the case with D d = 10 -3 m, the green line the case with D d = 10 -4 m and the blue line the case with D d = 10 -2 m. Vertical dashed lines illustrate the time when the minimum value of e is obtained, which is referred in the text as the overturn time.

Figure 8 :

 8 Figure 8: Time evolution of the particles concentration in the upper part of the reservoir with ∆ρ = 50 kg.m -3 and D d = 10 -3 m. The different color lines illustrate with cases with different viscosities of the continuous phase ranging from 0.1 to 10 Pa.s. Vertical dashed lines illustrate the time when the minimum value of e is obtained, which is referred in the text as the overturn time. The corresponding values of the overturn times are plotted in Fig.9. We note we have tested the case with ηc = 0.1 Pa.s with two different meshes (a 6 and 8 mm mesh) and the difference between the two curves is not significative.

Figure 13 :

 13 Figure 13: Density contrast ratio ∆ρ M O /∆ρ crit in the MO. Values larger than 1, or smaller than -1, indicate efficient Bg segregation below, or above, a mushy MO, respectively.

Figure 14 :

 14 Figure14: Density contrast ratio ∆ρ M O /∆ρ crit in the MO calculated as a function of the Bg grain size, from 0.1 to 5 mm, for 5 different melt compositions, from pyrolite -40%Bg to pyrolite +40%Bg. Values larger than 1 or lower than -1 indicate efficient Bg segregation below, or above, the mushy MO, respectively. χρ = 0 is found for ∆ρ M O = 0 (Fig.12). It indicates infinite value of Sh and perfect crystal suspension (Fig.10). Continuous and dashed lines indicate K(F e) Bg/M elt of 0.1 or 0.25, respectively.

Figure 15 :

 15 Figure15: Density contrast ratio ∆ρ M O /∆ρ crit calculated for various surface temperatures (Eq. 19). Thick, intermediate and thin lines correspond to 500, 1000 and 1500 K, respectively. The calculations are also performed for various mantle compositions and K(F e) Bg/M elt values.

  Thermal diffusivity κ = k/(ρ c,0 C p ) 1.7 × 10 -6 m 2 .s

  Figure 1.

  Figure 2.

  Figure 3.

  Figure 4.

  Figure 5.

  Figure 6.

  Figure 7.

  Figure 8.

  Figure 9.

  Figure 12.

  1-x-y Fe x Al y )(Si 1-y Al y )O 3 stoichiometry with a dominant fraction of Fe 2+ over Fe 3+ at lower mantle conditions, in agreement with pre-SiO 4 forsterite, Fe 2 SiO 4 fayalite, CaMgSi 2 O 6 diopside, and CaAl 2 Si 2 O 8 anorthite

	same P-T conditions, we calculate the MO density (ρ M elt ) and thermal
	expansion (α M elt ) using the EoS of the end-member molten MgSiO 3 en-
	statite, Mg 2
	vious works (e.g. Andrault et al. (2018)). This assumption has little
	impact on the Bg density; the dominating parameter is the global Fe con-
	tent in Bg, itself controled by K(F e) Bg/M elt (discussed above). At the

Table 1 :

 1 Symbol definitions and values of the physical parameters used in this study.

  to compute the heat flow F .

		η c (P a.s)	Ra	D d (m) ∆ρ (kg/m 3 )	Sh
	#1	0.1	1.8 × 10 8	10 -4	5	4.16
	#2	0.1	1.8 × 10 8	10 -4	10	2.08
	#3	0.1	1.8 × 10 8	10 -4	20	1.04
	#4	0.1	1.8 × 10 8	10 -4	30	0.69
	#5	0.1	1.8 × 10 8	10 -4	40	0.52
	#6	0.1	1.8 × 10 8	10 -4	50	0.42
	#7	0.1	1.8 × 10 8	10 -3	5	0.41
	#8 #9 #10 #11 #12 #13 #14 #15 #16 #17 #18 #19 #20 #21 #22 #23 #24 #25 #26 #27 #28 #46 #47 #48 #49 #50 #51 #52 #53 #54 #55 #56 #57 #58 #59 #60 #61 #62 #63 #64 #65	0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.4 0.4 0.4 0.4 0.4 0.4 0.7 0.7 4. 4. 4. 4. 4. 4. 10 10 10 10 10 10 10 10 10 10 10 10 10 10	1.8 × 10 8 1.8 × 10 8 1.8 × 10 8 1.8 × 10 8 1.8 × 10 8 1.8 × 10 8 1.8 × 10 8 1.8 × 10 8 1.8 × 10 8 1.8 × 10 8 1.8 × 10 8 9. × 10 7 9. × 10 7 4.4 × 10 7 4.4 × 10 7 4.4 × 10 7 4.4 × 10 7 4.4 × 10 7 4.4 × 10 7 3. × 10 7 3. × 10 7 4.4 × 10 6 4.4 × 10 6 4.4 × 10 6 4.4 × 10 6 4.4 × 10 6 4.4 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6 1.8 × 10 6	10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -2 10 -2 10 -2 10 -2 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -4 10 -3 10 -3 10 -3 10 -3 10 -2 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -4 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 10 -3 0.015 10 0.2 20 0.1 30 0.07 40 0.052 50 0.042 80 0.026 100 0.02 5 0.041 10 0.02 20 0.01 30 0.007 50 0.052 60 0.044 5 0.65 10 0.33 20 0.16 30 0.11 40 0.082 50 0.065 40 0.1 50 1.4 50 0.14 70 0.1 80 0.0875 100 0.07 50 0.014 5 19 10 9.51 20 4.76 30 3.17 40 2.37 50 1.9 100 0.95 5 1.9 10 0.95 20 0.47 30 0.31 40 0.24 50 0.19 10 -3 100 0.095
	#66	10	1.8 × 10 6	10 -2	5	0.19
	#67	10	1.8 × 10 6	10 -2	10	0.095
	#68	10	1.8 × 10 6	10 -2	15	0.063
	#69	10	1.8 × 10 6	10 -2	20	0.047
	#70	10	1.8 × 10 6	10 -2	30	0.032
	#71	10	1.8 × 10 6	10 -2	50	0.019
	#72	10	1.8 × 10 6	10 -2	100	0.0095

Acknowledgements

The authors thank Pierre Bonnand and Maylis Landeau for fruitful discussions.This research received funding from French PNP program (INSU-CNRS) and the French Government Laboratory of Excellence initiative No. ANR-10-LABX-0006. This paper is Laboratory of Excellence ClerVolc Contribution No. XXX.

Data availability

The data that support the findings of this study were obtained using the commercial software COMSOL Multiphysics (version 5.4). COMSOL Multiphysics is a simulation platform that provides fully coupled multiphysics and singlephysics modeling capabilities. COMSOL Multiphysics (version 5.4) has been previously validated for two phase flow applications [START_REF] Qaddah | Dynamics and stability of an iron drop falling in a magma ocean[END_REF][START_REF] Qaddah | Thermal evolution of a metal drop falling in a less dense, more viscous fluid[END_REF]).

To compute our simulations we used the Heat Transfer (https://www.comsol.com/heattransfer-module) and Computational Flow Dynamics (https://www.comsol.com/cfdmodule) modules in addition to the main Multiphysics platform (https://www.comsol.com/comsol-

Code availability

The software used for this study is the commercial software COMSOL Multiphysics (version 5.4) previously validated for two phase flow applications [START_REF] Qaddah | Dynamics and stability of an iron drop falling in a magma ocean[END_REF][START_REF] Qaddah | Thermal evolution of a metal drop falling in a less dense, more viscous fluid[END_REF]). A trial version of COMSOL Multiphysics may be requested (see https://www.comsol.com).


Pierre-Yves Ancelin

Julien Berthe

Alain Devos

When Yesterday'

When yesterday's war generates life today: example of the ecological and territorial resilience of

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

(pierre-yves.ancelin@univ-reims.fr)