Timescales of chemical equilibrium between the convecting solid mantle and over-/underlying magma oceans

Thank you for accepting our paper! We would like to thank you once again for giving us the opportunity to make changes in the previous version of the manuscript. Once again, we would like to thank the Editor Julien Aubert and the Reviewers for the useful comments that improved the manuscript. We would also like to thank the whole Editorial Board and the Editorial Support Team for handling our paper.

Mush is not considered.

Crystallisation of a magma ocean from the middle 55
If crystallisation of the magma ocean instead starts somewhere at mid-mantle depths and the crystals formed are near-neutrally buoyant (Labrosse et al., 2007;Boukaré et al., 2015), the first layer of solid mantle forms and separates the magma ocean into TMO and BMO (Fig. 1b). Then, two crystallisation fronts move in opposite directions: the TMO front progresses upwards until it reaches the surface of the planet, and the BMO front progresses downwards until it reaches the CMB. In this process, both TMO and BMO, as well as the related cumulates, become progressively enriched in FeO (Fig. 1b). In contrast to the TMO 60 cumulates (see above), BMO cumulates are likely formed over much longer timescales (Labrosse et al., 2007) and are expected to form a stable density profile. By the time the BMO is fully crystallised, a dense stable solid layer may persist at the base of the mantle. This dense layer may explain seismic observations that point to the existence of thermo-chemical piles near the CMB (Masters et al., 2000;Ni and Helmberger, 2001;Garnero and McNamara, 2008;Deschamps et al., 2012;Labrosse et al., 2015;Ballmer et al., 2016).

Motivation
Along these lines, the chemical evolution of the solid mantle depends on the history of early planetary melting and crystallisation. This is a history with either one or two magma oceans, and with convection in the solid mantle driven by unstable thermal and/or chemical stratification, probably while magma ocean(s) at the top and/or bottom are still present. While any such convection would imply re-melting of solid cumulates, the related consequences for mantle evolution are poorly understood. Only 70 a few numerical modelling studies have explicitly incorporated coupled re-melting and crystallisation at the magma ocean mantle boundary or boundaries Morison et al., 2019;Agrusta et al., 2019), and none of these studies have explored the consequences for chemical evolution.
In this paper we use a numerical model to investigate the thermo-chemical evolution of the solid mantle in contact with a TMO and/or a BMO. We consider that convection in the solid mantle starts before the end of magma ocean crystallisation, 75 therefore, dynamic topographies that may form at either or both solid mantle-magma ocean boundaries can melt or crystallise.
We do not explicitly account for the progression of the crystallisation front(s). However, we test several evolution scenarios and different magma oceans thicknesses. We determine the timescales of chemical equilibrium between the magma ocean(s) and the solid mantle, and compare them with those of progression of the crystallisation front (e.g., Lebrun et al., 2013). For simplicity, we hereafter use the term solid-liquid phase changes interchangeably with fractional crystallisation and melting 80 processes at the interface between the solid mantle and TMO and/or BMO.
2 Numerical model

Problem definition
We use the finite-volume/finite-difference method with the convection code StagYY (Tackley, 2008), to model the thermochemical evolution of the solid mantle during magma ocean crystallisation. We test three different evolution scenarios, as the 85 solid mantle may be bounded above by a TMO and/or below by a BMO (Fig. 2). We assume steady crystallisation front (s) and test different magma ocean thicknesses: when only one ocean is present, it can be 100, 500 or 1000 km thick; when both oceans are present, the thickness of each ocean is 100 and/or 500 km.
We assume that the solid mantle is an infinite Prandtl number fluid. We assume mechanical stability between the solid mantle and magma oceans, i.e., ρ TMO < ρ S < ρ BMO , with ρ TMO , ρ S and ρ BMO the densities of the TMO, solid mantle and BMO, 90 respectively. We take gravitational acceleration, g, viscosity, η, thermal diffusivity, κ, heat capacity, C p , thermal expansion coefficient, α, and compositional expansion coefficient, β, as constant. Values of these parameters can be found in Table 1.
We make equations dimensionless to reduce the number of parameters that describe the physical problem. Dimensions of distance, time and temperature can be recovered using, respectively, the thickness of the solid mantle, h S , the thermal diffusive timescale, Ocean (BMO) and right) TMO and BMO. Solid mantle is taken as a spherical shell with density ρS, thickness hS, inner radius R − , and outer radius R + . TMO and BMO are taken with densities ρTMO and ρBMO, and thicknesses hTMO and hBMO, respectively. Superscripts " + " and " − " at the boundaries refer to the boundary between TMO and solid mantle, and solid mantle and BMO, respectively.
We assume incompressibility in the Boussinesq approximation (e.g., Chandrasekhar, 1961). Therefore mass, energy, composition and momentum conservation equations are written as: with u the velocity field, T the lateral average of the temperature field T , t the time, X S FeO the FeO molar content in the solid mantle, X S FeO the lateral average of X S FeO , p the dynamic pressure, Ra the Rayleigh number and B the buoyancy number. The

105
last two are defined respectively as: In this study we consider that magma oceans and solid mantle are made only of (Fe, Mg)O (see section 2.3 for more details).
We set temperature to 1.0 and 0.0, respectively, at the bottom and top solid domain boundaries. Regarding the buoyancy 110 number, although Earth-like models point to a value of B ∼ 3.0 today, the value in the early Earth might be different. In this paper we make a conservative choice and consider a buoyancy number of B = 1.0, to attribute similar weight to compositional and thermal effects on the density.
The solid domain is represented using the spherical annulus geometry (Hernlund and Tackley, 2008), composed of a grid of 512 × 128 cells, in which Eq.
(2) -Eq. (5) are solved. Composition is advected by tracers. We assume that each magma ocean 115 is well-mixed and that its dynamics are fast compared to that of the solid mantle. In our setup, magma oceans are treated as simple 0D compositional reservoirs at solid mantle boundaries. We hereafter use superscripts ' + ' and ' − ' to refer, respectively, to top and bottom solid mantle boundaries. In equations, the sign '±' reads as '+' if a TMO is considered, and '−' if a BMO is considered. The subscript ' MO ' refers to Magma Ocean. Thus, when we introduce a quantity, e.g. ξ, related to a magma ocean, we introduce it as ξ ± MO , with ξ + MO = ξ TMO relating to the TMO, and ξ − MO = ξ BMO relating to the BMO.

Dynamic topography and the phase change boundary condition
Since convection in the solid mantle likely starts before the end of magma ocean crystallisation (Maurice et al., 2017;Ballmer et al., 2017b;Boukaré et al., 2018;Morison et al., 2019;Miyazaki and Korenaga, 2019b), dynamic topographies are supported at either or both solid mantle boundaries. The timescale for producing dynamic topography is noted τ η . This topography can be eroded by solid-liquid phase changes on a timescale related to the transfer of energy and FeO through the magma ocean, 125 from material that is crystallising to material that is melting. We denoted this timescale by τ φ .
The relative values of the two timescales, τ η and τ φ , control the dynamical behaviour of the boundary. If τ η τ φ , dynamic topography can build before being erased by the phase change. In this case, dynamic topography is only limited by the balance between viscous stress in the solid and the buoyancy associated with the topography. In the limit of small topographies, this leads to the classical non-penetrating free-slip boundary condition in which the radial velocity of the solid effectively goes to 0 130 at the boundary (Ricard et al., 2014). On the other hand, if τ η τ φ , the topography is erased faster by phase changes than it can be built by viscous stress in the solid. Consequently, this removes the stress imposed by the topography and the associated limit to the radial velocity. These processes are incorporated into our boundary condition, described by the phase change number, considering that when Φ → ∞, dynamic topography is built (or relaxes) by viscous forces, and when Φ → 0 it is eroded 135 by melting or fractional crystallisation processes (Deguen, 2013;Deguen et al., 2013). The related phase-change boundary condition in dimensionless form, at either or both solid mantle boundaries is: with u r the vertical velocity of the flow in the solid mantle. On one hand, this boundary condition can act like a non-penetrating free-slip boundary condition when Φ → ∞, since vertical velocities of the solid flow tend to 0 at the boundaries. Under this 140 boundary condition, transfer of material across a solid mantle-magma ocean boundary cannot occur. On the other hand, this boundary condition can act as being "open" to phase changes when Φ → 0, since these vertical velocities will be non zero at the solid mantle-magma ocean boundary, and a significant flux of solid and liquid material can cross it to melt and crystallise.
Hence, transfer of material across the phase-change boundary is efficient. In the extreme case of Φ = 0, this boundary condition corresponds to free in-and outflow.

145
The specific value of Φ is difficult to constrain (because τ φ is non-trivial to determine), and also is expected to vary with time (i.e., because τ η depends on the thickness of the solid mantle) (Deguen, 2013;Deguen et al., 2013). However, for a purely thermal case, Morison et al. (2019) and Morison (2019) estimate Φ + ∼ 10 −5 and Φ − ∼ 10 −3 for the Earth. Therefore, significant transfer of material across the solid-mantle magma-ocean boundaries is expected. However, also consider that real multi-phase rocks typically melt over large pressure ranges, unless for truly eutectic bulk compositions. The depleted residue 150 of mantle melting may restrict the efficiency of material transfer across the solid mantle magma-ocean boundaries, depending on the efficiency of melt-solid segregation near the boundaries. In addition to the expected temporal evolution of Φ ± , this potential restriction motivates our exploration of a broad range of Φ ± . In this study we use 7 values of Φ ± that range from 10 −1 to 10 5 . We use Φ = 10 −1 as the lowest value possible for Φ ± because the resolution of the thermal boundary layer is computationally demanding once Φ ± decreases below 10 −1 .
155 Deguen et al. (2013) and Labrosse et al. (2018) found that the critical Rayleigh number, Ra c , for the solid mantle is strongly sensitive to Φ and the setup considered, i.e., having a TMO and/or a BMO, as well as to the thickness of the solid layer. For instance, if the solid mantle is bounded by a TMO of 100 km and Φ → ∞, Ra c is on the order of 10 3 , but for small Φ, is on the order of 10 2 . Ra c can even decrease to arbitrarily small values on the order of ∼ Φ if a TMO and BMO are both considered.
Therefore, we also systematically vary the Rayleigh number, Ra, which controls the convective vigour of the mantle. We choose 160 Ra as multiples of Ra c , according to the super-criticality factor, SC: We use 4 values of SC ranging from 10 2 to 10 5 .

Compositional treatment
In this study we consider a simplified compositional model with only two components, FeO and MgO, which are thought to be 165 the Fe-rich and Mg-rich end-members of mantle silicates. We denote the FeO and MgO molar content in the solid and magma ocean parts, respectively, by X S FeO and X MO FeO , and X S MgO and X MO MgO . We consider mass balance between FeO and MgO in the solid mantle and magma oceans, therefore, X S FeO + X S MgO = 1 and X MO FeO + X MO MgO = 1. Our model simulates melting and crystallisation of material depending on the influx and outflux of material at the solid mantle boundary. Melting of solid material is simulated when dynamic topography develops outside the solid domain, i.e., 170 when there is an outflux of material of the solid domain. It is assumed that no fractionation operates when the solid melts, i.e., all (Fe,Mg)O present in this topography goes into the magma ocean. Therefore, tracers that leave the solid domain pass their information (about mass and composition) to the magma ocean, and are deleted.
We simulate crystallisation of the magma ocean when negative dynamic topography develops in the solid domain, i.e., when there is an influx of mass in the solid domain. When this happens, the influx of material pushes tracers and cells near the 175 boundary are left with no tracers. To ensure mass conservation, new n tracers are introduced in those cells, which simulates solid mantle being created. To determine n, i.e., the amount of tracers that need to be introduced in the solid part, we calculate the influx of mass corresponding to this dynamic topography, and divide it by the tracer ideal mass. We then equally distribute the influx of mass by n tracers. The composition of the solid created is related to that of the liquid by fractional crystallisation, therefore, only a fraction of FeO goes into the solid. This fraction is given by the distribution coefficient, K: We assume K = 0.3 (e.g., Corgne and Wood, 2005;Liebske et al., 2005). The difference between influx and outflux of material through the boundary is of the order of 10 −15 , meaning that conservation of mass in the solid mantle is ensured.
In this paper, we attempt to estimate the characteristic timescale to establish chemical equilibrium between the solid mantle and the magma ocean(s). Assuming a full equilibrium between the solid mantle and magma oceans (superscript " Eq "), the FeO 185 content in the bulk, X Bulk FeO , can be expressed as function of the volumes (V S , V TMO and V BMO ) and the FeO content (X S,Eq FeO , X TMO,Eq FeO and X BMO,Eq FeO ) in the solid mantle and magma oceans, From Eq. (11) and Eq. (12) one can find the FeO content in the solid mantle when it is in chemical equilibrium with the magma ocean(s): where: But, because chemical equilibrium would take too long to reach in a reasonable run time, we look for the timescale to reach chemical half-equilibrium. Starting with a FeO content in the solid mantle equal to X S, Ini FeO , the half-equilibrium is reached when the solid mantle reaches the content X S, Eq/2 FeO , defined as: We denote by t S, Eq/2 the time at which the solid mantle reaches chemical half-equilibrium, t S, Eq/2 = t(X S, Eq/2 FeO ).

200
Previous studies suggest that the Fe content of the present day bulk silicate Earth is 0.113 (Taylor and McLennan, 1985) or 0.107 (McDonough and Sun, 1995). We suppose that some of the Fe could migrate to the core with time (e.g., Nguyen et al., 2018) and therefore, in this study we use X Bulk FeO = 0.120. We start the simulations with a homogeneous FeO content in the solid mantle and magma ocean(s), X S, Ini FeO = X TMO, Ini FeO = X BMO, Ini FeO = 0.120. Although this initial composition is not consistent with the fractional crystallisation assumed in this problem, it serves well our goal of measuring the timescale to reach chemical 205 equilibrium between solid and liquid reservoirs.

Chemical evolution of the mantle bounded on top by a TMO
In this subsection we investigate how the chemical evolution of the solid mantle is affected by the efficiency of mass transfer across the phase-change boundary, as controlled by Φ. As mentioned in the previous section, low values of Φ correspond to 210 efficient material transfer across the phase-change boundaries, and high values of Φ correspond to inefficient material transfer, similar to classical convection. We analyse first the case of a solid mantle bound above by a TMO, as the most straightforward  Because the parameter space explored in this paper is vast, we illustrate here as an example the chemical evolution of a solid 215 mantle bounded by a TMO of 500 km, under three different values of phase change number, Φ + = 10 −1 , 10 2 , 10 3 , at the same super-criticality value of SC = 10 5 . As the critical Rayleigh number, Ra c , decreases as Φ + decreases, these three cases have different values of Rayleigh number, Ra. Hence, for Φ + = 10 −1 , 10 2 , 10 3 , Ra = 100 × 10 5 , 635 × 10 5 , 687 × 10 5 , respectively.
According to Eq. (12) and Eq. (13) the system does not start in chemical equilibrium. We determine the time needed to reach chemical half-equilibrium. Figure 3 shows the chemical evolution in dimensionless time units of these three cases. Our models predict that regardless of the value of Φ + , the FeO content in the solid mantle decreases towards X S,Eq FeO , and the FeO content in the TMO increases towards X TMO,Eq FeO , thereby bringing the solid mantle and the TMO close to chemical equilibrium (but not chemical homogeneity as seen later). However, the lower the value of Φ + , the faster half-equilibrium is reached, since it effectively increases the 225 exchange of material between reservoirs. We calculate the time needed to reach chemical half-equilibrium, and for Φ = 10 −1 , half-equilibrium is reached ∼ 10 times faster than for Φ = 10 2 , and ∼ 200 times faster than for Φ = 10 3 .
In Fig. 4 we present snapshots of FeO content in the solid mantle for these three cases. Our models show that dynamics in the solid mantle is very different between cases. With Φ + = 10 −1 (Fig. 4a), mantle flow is dominated by degree-1 convection, which persists stably for the whole simulation time. With this pattern of convection, there is an upwelling of primordial 230 material (in yellow) that melts on one hemisphere, while material from the TMO crystallises at the boundary and forms a downwelling on the other hemisphere (in blue). This downwelling is FeO depleted, which introduces a strong heterogeneity in the solid mantle. Degree-1 convection involves very little deformation, which explains the existence of a considerable amount of primordial material in the solid mantle, even around the time at which chemical half-equilibrium occurs (snapshot inside red box). As Φ + increases ( Fig. 4b and Fig. 4c), higher degree modes of convection with several convection cells appear. Although . Snapshots of FeO content in the solid mantle, X S FeO , as function of dimensionless time (factor of 10 −6 inside each annulus) for the cases presented in Fig. 3. a) Φ + = 10 −1 , b) Φ + = 10 2 and c) Φ + = 10 3 . In these simulations the top magma ocean thickness is 500 km, and there is no basal magma ocean. Super-criticality is SC = 10 5 in all three cases shown. Snapshots with a red box indicate that the model time is close to the chemical half-equilibrium time, t S, Eq/2 . For Φ + = 10 −1 , 10 2 and 10 3 , t S, Eq/2 = 9.6 × 10 −6 , 114.9 × 10 −6 and 1895.5 × 10 −6 , respectively. Magenta contours correspond to the streamlines of the flow. the composition of the TMO and the average composition of the solid mantle tend to mutual chemical equilibrium in all three cases, chemical homogeneity across the solid mantle is not necessarily reached.
Our models show that for other evolution scenarios, i.e., solid mantle in contact with just a BMO and with a TMO and BMO, the system also evolves to a state close to chemical equilibrium but not chemical homogeneity. In Fig. 5 we present snapshots of FeO content in the solid mantle for different evolution scenarios at about the time of chemical half-equilibrium. When the 240 solid mantle is in contact with just a BMO, material from the magma ocean crystallises at the boundary and forms upwellings (in blue). This material is FeO depleted and, similarly to the TMO case, introduces a strong heterogeneity in the solid mantle around the half-equilibrium time. When the solid mantle is in contact with both oceans, convection occurs with degree-1, i.e., material of the TMO and of the BMO crystallises at the corresponding boundary and forms a downwelling (in blue) and an upwelling (in blue and green), respectively. Note that in this scenario, since the volume of the BMO is smaller than that of the 245 TMO, the BMO composition changes rapidly. Therefore, the composition of the upwelling changes rapidly as well (colours from blue to green). The degree-1 pattern of convection persists stably for the whole simulation time. for a solid mantle in contact with only BMO and with TMO and BMO, did not reach chemical half-equilibrium. Note that the difference between the aspect ratio of each evolution scenario is too small to be noticed in these annuli.
3.2 Timescales of chemical half-equilibrium between the solid mantle and magma ocean(s) Figure 6 shows the timescales of chemical half-equilibrium for the scenarios explored in the previous subsection. For a wide range of SC and Φ ± , these timescales are shown for a solid mantle bounded by a TMO of 500 km thickness (Fig. 6a), by a 250 BMO of 500 km thickness (Fig. 6b), and by a TMO and a BMO of 500 km and 100 km thickness, respectively (Fig. 6c). For all evolution scenarios, models predict that timescales of chemical half-equilibrium decrease for decreasing Φ ± . In other words, chemical half-equilibration is more efficient for efficient material transfer across the solid mantle-magma ocean boundaries. Our results also show that the timescales of chemical half-equilibration are similar (i.e., of the same order of magnitude) for a given SC and Φ ± ranging between 10 −1 and 10 1 , independent of the evolution scenario. This shows that below Φ ± = 10 1 255 a regime with efficient material transfer across the solid mantle-magma ocean boundaries is established. The transition to the regime of inefficient material transfer (i.e., in which mantle flow is limited by viscous building of dynamic topography) occurs somewhere between Φ ± = 10 1 and 10 2 . In this regime, timescales of half-equilibration systematically increase with Φ ± . Our models predict that this transition between regimes occurs over a similar interval of Φ ± for other thicknesses of TMO and/or BMO.

260
To obtain an empirical scaling law, we fit the predicted timescales, t S, Eq/2 pred , for all simulations that reached chemical halfequilibrium. The fitting equation provides t S, Eq/2 pred in dimensionless form as a function of Ra, Φ ± and V S : with V M the volume of the present-day Earth's mantle. Coefficients of this equation can be found in Table 2. In the Appendix A of this paper, we explain the regression method and show a good agreement between the timescales to reach chemical 265 half-equilibrium, obtained with our model predictions and our empirical scaling law (Fig. A1).
Equation (15) presents two branches, each corresponding to a different regime: the left branch corresponds to the regime of efficient material transfer across the solid mantle-magma ocean boundaries, and the right one to the regime of inefficient material transfer. Our models predict that in the regime of efficient material transfer (i.e., for low values of Φ), timescales to reach chemical half-equilibrium are only loosely affected by the volume of the solid mantle (or in other words, by the volume 270 of the magma ocean(s)). The volume of the solid only systematically affects the timescales once the regime shifts to the one of inefficient material transfer. This conclusion is independent of the evolution scenario. One possible explanation for the weak influence of the solid mantle's volume is that at low values of Φ, convection occurs with low degree, so the geometry of the problem is less important.
We find that timescales to reach chemical half-equilibrium are about a factor of 3 larger for a solid mantle in contact with 275 just a TMO (Fig. 6a) than for a solid mantle with just a BMO (Fig. 6b). This can be explained by the fact that the geometry of Table 2. Results of the regressions of the timescales of chemical half-equilibration using the form: t S, Eq/2 the problem is different in both cases. Although the TMO and the BMO have the same thickness, the volume of the TMO is larger than that of the BMO by roughly a factor of 3, which explains the increased duration to reach the half-equilibrium FeO content in the magma ocean.
When it comes to a solid mantle bounded by both oceans (Fig. 6c), models predict that timescales are roughly 2 orders of 280 magnitude smaller than the ones of a solid mantle in contact with just a TMO (Fig. 6a), for a given Rayleigh number. This result is explained by two effects. The critical Rayleigh number is much lower when two magma oceans are present than when only one is present. In principle, when both magma oceans are present, the critical Rayleigh number can even be arbitrarily low as Φ ± decreases towards 0 . Moreover, Agrusta et al. (2019) showed that the heat flow and RMS velocity in the solid mantle vary linearly with Ra when both magma oceans are present, whereas heat flow and RMS velocity in the solid 285 mantle vary as Ra 1/3 and Ra 2/3 , respectively, in the case of only one magma ocean present. This further increases the difference between the two scenarios at a given value of the Rayleigh number. Therefore, one should expect that the timescales to reach chemical half-equilibrium may be arbitrarily low, depending on the efficiency of material transfer across the BMO-mantle and TMO-mantle boundaries.

Chemical half-equilibrium and crystallisation timescales 290
In this subsection we compare the timescales to reach chemical half-equilibrium between a solid mantle and a TMO of a given thickness, with timescales of crystallisation of such a TMO as calculated for the Earth case. The timescales of TMO crystallisation (i.e., before reaching the mush stage) are given by Lebrun et al. (2013), hereinafter denoted by t C L13 . We take the solid mantle bounded on top by a TMO of 100, 500 and 1000 km thickness, and use Eq. (15) to determine the timescales of half-equilibration as a function of phase change number, Φ + . In an attempt to apply our fitting equation to  Figure 7. Timescales to reach chemical half-equilibrium between the solid mantle and a TMO of 100, 500 and 1000 km (respectively orange, blue and purple solid lines, this study, t S, Eq/2 pred ), and timescales of crystallisation of that TMO until it is completely mushy (thicknesses with same colours in dashed lines, from Lebrun et al. (2013), t C L13 ), versus different values of phase change number, Φ + . a) corresponds to a Rayleigh number of the solid part of RaM = 10 8 and b) to RaM = 10 9 . Rose background corresponds to a time higher than the age of the Earth. We note that t C L13 are consistent with the timescales constrained by Nikolaou et al. (2019) and Salvador et al. (2017).
Earth, we assume that the global Rayleigh number of the early-Earth mantle just after solidification of the TMO is between Ra M = 10 8 and Ra M = 10 9 . Ra M is calculated on the basis of the total thickness of the solid mantle, h M . The Rayleigh number, Ra, used in Eq. (15) is then re-scaled to the actual thickness of the solid mantle (i.e., before solidification of the TMO) as follows:

300
This re-scaling neglects the change of various physical parameters (from Eq. (6)), but is sufficient for our discussion.
The comparison between timescales is presented in Fig. 7. The timescale to crystallise the TMO is loosely dependent on its thickness, and this time is around 1 Myr. Our models predict that there are significant chemical exchanges between the TMO and the solid mantle for Φ + <10 and <100, for Ra M = 10 8 and Ra M = 10 9 , respectively. Therefore, for Φ + smaller than these values, the TMO is expected to have reached (at least) chemical half-equilibrium with much of the mantle before reaching 305 the mush stage. Therefore, a very strong enrichment of the final-stage TMO as predicted by fractional crystallisation models (Elkins-Tanton et al., 2003) is not expected to occur for small-to-moderate Φ + .
Increasing the thickness of the TMO, hence, decreasing the thickness of the solid mantle, decreases the Rayleigh number and the ratio V S /V M , which make the dimensionless time increase (see fitting). But also, decreasing the thickness of the solid mantle decreases the scale for time (h 2 S /κ), which partially compensates the increase mentioned previously when recovering the dimensional time. As a result, the thickness of solid in the phase-change regime only loosely affects the dimensional half-equilibration time (Fig. 7). The effect of the thickness of solid is a bit stronger in the high-Φ regime (as mentioned before).

Discussion
Our models address the compositional evolution of the solid mantle bounded by magma oceans above and/or below, and constrain the time needed to chemically (half-)equilibrate these reservoirs. While the concept of a TMO that potentially interacts 315 with the underlying solid mantle is now well accepted, the idea of a long-lived BMO remains controversial. Whether or not a BMO can be stabilised depends on the slope of the adiabat vs. that of the melting curve (Labrosse et al., 2007), and/or on the fate of FeO-rich TMO cumulates that sink to the CMB (Labrosse et al., 2015;Ballmer et al., 2017b). Regardless of these issues, our study can be applied to various scenarios, including those with a solid mantle bounded just by a TMO, and bounded by a TMO and BMO.

320
Classical fractional crystallisation models predict strongly inverse chemical stratification of the initial solid mantle, and consequently a global-scale overturn by the end of a TMO crystallisation (e.g., Elkins-Tanton et al., 2003. The propensity of such a massive density-driven overturn depends on whether or not chemical equilibration between the solid mantle and magma ocean(s) can occur before magma-ocean solidification. Any style of magma-ocean crystallisation leaves a strongly super-adiabatic thermal profile, which should drive convection in the cumulate layers before full solidification (Solomatov and 325 Stevenson, 1993a;Ballmer et al., 2017b;Maurice et al., 2017;Boukaré et al., 2018), and related vertical flow should promote (partial) melting near the TMO-mantle (and BMO-mantle) boundary/ies. In this study we focus on a phase change boundary condition that allows material to flow through the boundary/ies and continuously change the composition of solid and liquid reservoirs. While flow through the boundary/ies remains penalised by the phase-change number, Φ, this approach implies final degrees of melting of at least ∼40%, i.e. beyond the rheological transition, at which crystals become suspended (Abe, 1997;330 Costa et al., 2009). Future work using a more realistic melting model is needed to test whether these high degrees of melting can indeed be reached, or quantify the effects of partial melting on equilibration between the mantle and TMO/BMO.
Our models show that the composition of the solid mantle and magma oceans strongly depends on the phase change number, Φ. In this study we take Φ as being constant through time, but because this number depends on the dynamics and thicknesses of the magma oceans, Φ may change continuously in a more realistic model with moving boundaries. Although we show 335 that chemical equilibration can occur before full crystallisation of the magma oceans, variations of Φ and a moving-boundary scheme should be considered in further studies.
Considering that Φ ± values are low when TMO and BMO (or just TMO) crystallisation starts Morison, 2019), mantle convection would first assume a degree-1 pattern (Fig. 5), possibly with implications for the origin of crustal dichotomy on the Moon (e.g., Ishihara et al., 2009) and Mars (e.g., Roberts and Zhong, 2006;Citron et al., 2018). However, 340 it remains to be shown that such a degree-1 pattern of convection would be able to survive through all stages of magma-ocean crystallisation.
The crystallisation fronts move at different speeds, since the TMO can crystallise in a few Myr years (Lebrun et al., 2013;Salvador et al., 2017;Nikolaou et al., 2019), whereas the BMO may persist for much longer (e.g., Labrosse et al., 2007Labrosse et al., , 2015. Therefore, Φ ± would change accordingly. The efficiency of equilibration during the late-stage magma ocean depends on the timescale of freezing of this final stage, as well as on the efficiency of mass transfer (Φ + ) for a thin and partially mushy TMO.
Once the TMO is fully crystallised, Φ + tends to infinity, while Φ − assumes a finite value as long as the BMO is still present.
Dynamics in the solid mantle would change accordingly: convection in the solid mantle may be either dominated by degree-This implies that the BMO is likely to crystallise much faster than suggested by Labrosse et al. (2007Labrosse et al. ( , 2015 for low Φ − , at least as long as no dense FeO-enriched materials accumulate at the bottom of the solid mantle to prevent efficient mass transfer across the BMO-mantle boundary. The BMO may even be thermally coupled to the relatively fast-cooling TMO for low Φ + and low Φ − . On the other hand, it is conceivable that thermally-coupled TMO and BMO crystallise more slowly than expected for a 355 thermally-isolated TMO (Agrusta et al., 2019). The presence of a BMO makes heat transfer across the mantle and out of the thermally coupled BMO and core more efficient than for cases without a BMO and with a boundary layer at the CMB instead. Such a situation implies that there is a larger heat reservoir available to buffer the temperature of the TMO for a given heat flux through the atmosphere and to space. Note that the mass and heat capacity of the core are similar to that of a 1000 km-thick magma ocean, hence the timescale for full crystallisation could be roughly twice than usually computed (cf. Lebrun 360 et al. (2013)). However, the timescales of TMO, BMO and core cooling would be largely unaffected if the BMO were thermochemically stratified (Laneuville et al., 2018). Whether or not material transfer across the whole mantle, as predicted here for cases with low Φ ± , can efficiently cool the core has important implications for the long-term thermal evolution of terrestrial planets, and the propensity of an (early) dynamo.
Even though timescales of BMO crystallisation are not well constrained, chemical exchange between the two magma oceans 365 (through the solid mantle) is still likely to occur. Note that the same process (i.e. mantle convection) that takes out heat from the BMO and core, is responsible for this chemical exchange. As an example for the Earth, if we take a TMO and a BMO of 100 km thickness each, Φ ± ≤ 10 and a Rayleigh number of 10 8 , we would expect a half-equilibrium between solid mantle, TMO and BMO in less than ∼ 460 ky, i.e. before TMO crystallisation (and even more so, before BMO crystallisation). This chemical exchange, however, does not necessarily imply homogeneity between the TMO and BMO, because the relevant 370 phase diagrams that control fractional crystallisation at the TMO-solid mantle (low pressures) and BMO-solid mantle (high pressures) boundaries are very distinct (e.g., Thomas et al., 2012;Boukaré et al., 2015). For example, while the FeO distribution coefficient, K defined in Eq. (11), is taken as constant in this study, its value is likely to be pressure-dependent (Nomura et al., 2011;Miyazaki and Korenaga, 2019a), potentially causing partitioning of FeO into the BMO. Regardless, any such exchange between the TMO, solid mantle and BMO could be a way to sequester trace elements (including heat-producing elements) into 375 the BMO, particularly if the TMO freezes faster than the BMO.
Once both oceans crystallise and Φ ± → ∞, convection in the solid mantle likely changes to higher degrees of convection (as already seen with Φ = 10 2 , 10 3 in Fig. 4b and Fig. 4c), similar to present-day Earth-mantle dynamics. Our models predict a largely homogeneous solid mantle, with some regions preserving significant primordial heterogeneity for long timescales.
The preservation of heterogeneity is likely to be enhanced once composition-dependent rheology (i.e., a difference in intrinsic 380 strength of mantle materials) is considered (Manga, 1996;Ballmer et al., 2017a;Gülcher et al., 2020). Indeed, primordial cumulates formed in the lower mantle may be strongly enriched in MgSiO 3 bridgmanite (Boukaré et al., 2015), and hence intrinsically strong (Yamazaki and Karato, 2001). In the present-day, LLSVPs are perhaps the most prominent and seismicallyevident large-scale mantle heterogeneities. That they are only rather mildly Fe-enriched (Deschamps et al., 2012) points to rather efficient equilibration between the magma ocean(s) and much of the solid mantle during crystallisation, such as predicted 385 by a subset of our models. For the Earth, the subset of our models with Φ + smaller than ∼100 suggests that chemical (half-)equilibrium between a solid mantle and a (100 to 1000 km-thick)TMO can be accomplished in less than ∼ 1 Myr, i.e., before the TMO is fully solidified or becomes a mush (Fig. 7). Equilibration in such sort timescale makes the solid mantle mostly homogeneous (with some heterogeneities as seen before), which could explain the pyrolitic nature of the mantle (Wang et al., 2015;Zhang et al., 2016;Kurnosov et al., 2017). Therefore, the final-stage TMO and subsequent mush may be efficiently 390 equilibrated with most of the solid mantle. In this case, we expect solid compositions that are by far not as enriched in FeO as predicted by fractional crystallisation models (e.g. Elkins-Tanton, 2012), in which strong enrichment only occurs because the final-stage TMO is fully separated from the solid mantle, with strong disequilibrium between the two reservoirs.
Similarly, we expect moderate enrichment (in FeO and incompatible trace elements) and roughly basaltic-to-komatiitic (i.e., the melting product of a hot ∼pyrolitic mantle) major-element compositions of the primary crust. As our models do not 395 explicitly address the final and mush stages of the TMO, and consider only a strongly simplified compositional model with only two components, (Fe, Mg)O, more detailed studies with a more complex compositional treatment are needed in order to predict the composition of the early crust.
In our models, we consider a simplified initial condition with bulk-planetary TMO and BMO compositions (X Bulk FeO = 0.120 and X TMO, BMO FeO = 0.120). While this condition may be realistic for a formation of the solid mantle due to equilibrium crystalli-400 sation, the TMO and BMO would be significantly more FeO enriched initially if they were formed by fractional crystallisation (Solomatov and Stevenson, 1993a;Xie et al., 2020). In our models, the initial cumulate downwellings formed at the TMOsolid mantle boundary are depleted in FeO and hence buoyant, resisting solid-mantle convection and delaying compositional equilibration, but this effect would be strongly diminished for a more realistic initial condition. Conversely, the initially depleted cumulate upwellings from the BMO-mantle boundary in our models advance convection and equilibration. As these 405 effects that depend on our choice of the initial condition scale with buoyancy number, B, we choose a conservative value of B = 1.0 (see method section). Also note that our models demonstrate that compositional equilibration between the TMO and/or BMO and solid mantle occurs swiftly over a wide range of solid-mantle thicknesses. This result implies that cumulate downwellings/upwellings from the TMO/BMO-solid mantle boundary have a similar composition than the ambient mantle, largely removing the chemical effects on convection that scale with B. Thus, the actual value considered for B is expected to 410 have only a minor effect on the equilibration timescales constrained here.
Smaller planets than Earth are less likely to be chemically equilibrated for a given bulk composition. First, they tend to cool faster, as they contain a smaller total reservoir of heat and volatiles (i.e., stabilising a less massive atmosphere to shield cooling). Moreover, the Ra number is lower for small planets, such that equilibration is expected to take longer according to our results. Thus, the Martian mantle might be less equilibrated (more stratified) than that of Earth (Elkins-Tanton et al., 2003; to be stabilised in their interiors due to high CMB pressures (Stixrude, 2014;Caracas et al., 2019), which has a strong effect on equilibration timescales. Whether or not chemical equilibration during the magma-ocean stage is efficient has important implications for the composition of the primary crust, the propensity of overturn and related stabilisation of a deep dense layer, as well as the long-term evolution of terrestrial planets.

Conclusions
In this work we use a numerical model to investigate the thermo-chemical evolution of the convecting solid mantle bound at top and/or bottom by magma oceans. We parameterise fractional crystallisation and melting processes of dynamic topography at either or both solid mantle boundaries, and determine the timescales to reach chemical half-equilibrium between solid mantle and magma ocean(s). 425 We show that these fractional crystallisation and dynamic melting processes at either or both boundaries play an important role in the chemical evolution of the solid mantle. Efficient transfer of FeO across the mantle-TMO and/or mantle-BMO boundary can prevent strong enrichment of the last-stage magma ocean, and thereby any strong chemical stratification of the early fully-solid mantle. Moreover, this efficient transfer of FeO renders the timescales of chemical (half-)equilibration between the solid mantle and magma ocean(s) shorter than (or on the order of) 1 Myr. Since magma ocean crystallisation occurs in few 430 Myr (Abe, 1997;Lebrun et al., 2013;Salvador et al., 2017;Nikolaou et al., 2019), our study suggests that chemical equilibrium between solid and liquid reservoirs can be reached before the end of magma ocean crystallisation. Therefore, a strong chemical stratification of the solid mantle is unlikely to occur, and the first crust is not expected to be extremely enriched in FeO.
This prediction fundamentally contrasts with that of classical models of fractional crystallisation of the magma ocean (e.g. Elkins-Tanton, 2012).

435
However, more studies are needed to better constrain chemical-equilibration timescales. This could be achieved, for instance, as more realistic compositional models and phase diagrams are accounted for, and/or a moving boundary approach is applied to explicitly model the evolution of either or both crystallisation fronts.

Appendix A: Regression Method
The best fitting coefficients of all regressions are obtained using a simple algorithm. Each free parameter has an initial possible 440 minimum and maximum, chosen here between -1 and 1. All 8 parameters a i in Eq. 15 are scanned between these minimum and maximum boundaries using homogeneous steps. For each point in that 8-D space, we compute the misfit between predicted and observed timescale as log cases t S,Eq/2 pred − t S,Eq/2 case 2 . The set of best fitting parameters are found by selecting the lowest misfit between the analytical formulation and the data.
the initial boundaries in parameter space. Iterations of the search are performed until the solution is converged below fourth digit precision. Figure A1 shows the regression for all cases that reached chemical half-equilibrium.