Solidification of Earth's mantle led inevitably
to a basal magma ocean

Charles-Édouard Boukaré¹,²*, James Badro¹† and Henri Samuel¹†

¹Université Paris Cité, Institut de Physique du Globe de Paris,
CNRS, 1 rue Jussieu, 75005 Paris, France.

²York University, Department of Physics and Astronomy, 4700 Keele Street, M3J 1P3 Toronto, Canada.

*Corresponding author(s). E-mail(s): boukare@ipgp.fr;

Contributing authors: samuel@ipgp.fr; badro@ipgp.fr;

†These authors contributed equally to this work.

Abstract

Current interpretation of deep-rooted geophysical structures in the man-
tle¹ is that they stem from the top-down solidification of Earth's
primitive basal magma ocean above the core.²,³,⁴,⁵,⁶ However, it remains
debated whether solids first formed at the bottom of the mantle, solid-
ifying upward, or above the melts, solidifying downward. We show
that gravitational segregation of dense, iron-rich melts from lighter,
iron-poor solids drives mantle evolution, regardless of where melting
curves and geotherms intersect. This process results in the accumula-
tion of iron oxide rich melts above the core, forming a basal magma
ocean. We numerically model mantle solidification using a novel multi-
phase fluid dynamics approach that integrates melting phase relations
and geochemical models. This allows estimating the compositional sig-
nature and spatial distribution of primordial geochemical reservoirs,
that may be directly linked to the isotopic anomalies measured in
Archean rocks.⁷,⁸,⁹,¹⁰,¹¹ Interestingly, we find that a significant amount
of solids is produced at the planet's surface, not at depth, inject-
ing geochemical signatures of shallow silicate fractionation in the deep
mantle. This work could serve as a foundation for re-examining the
intricate interplay between mantle dynamics, petrology, and geochem-
istry during the first billion years of the evolution of rocky planets.

001

002

003

004

005

006

007

008

009

010

011

012

013

014

015

016

017

018

019

020

021

022

023

024

025

026

027

028

029

030

031

032

033

034

035

036

037

038

039

040

041

042

043

044

045

046
047

048

049

050

051

052

053

054

055

056

057

058

059

060

061

062

063

064

065

066

067

068

069

070

071

072

073

074

075

076

077

078

079

080

081

082

083

084

085

086

087

088

089

090

091

092

**Keywords:** Mantle Structure and Dynamics, Mantle Differentiation, Magma Oceans, Rocky Planets Formation and Evolution.

Evidence of primordial silicate differentiation and magma ocean hypothesis

Isotopic anomalies in short-lived radiogenic isotopic systems in mantle rocks that record magmatic differentiation processes occurring in the first 100 Myrs shows that the Earth's mantle has preserved chemical heterogeneities$^{7,8,9,10,11}$ dating back to its infancy. These findings are corroborated by the noble gas geochemical record that argues for the preservation of such early formed geochemical reservoirs.$^{12,13,14}$ The solidification of a deep primitive magma ocean alone can explain such an early silicate differentiation event,$^{8,9}$ as it produces the observed fractionation in several trace-element ratios at a global scale. This solidification process can also explain the current seismic structure of the deep mantle, where Large Low Shear Velocity Provinces (LLSVPs) and Ultra-Low velocity zones (ULVZs) can be interpreted as residual products of primordial magma ocean solidification.$^{2,3,4,5}$ As remnants of magma ocean solidification, the two antipodal deep-rooted (above the core-mantle boundary) LLSVPs$^1$ must play a leading role on global mantle and core dynamics,$^{15}$ plate tectonics and hot-spot magmatism$^{16,17,18,19}$ during the entire history of the Earth. Therefore, understanding magma ocean solidification from a dynamical and petrological point of view is essential for our comprehension of the long-term evolution of Earth's mantle and its present-day state.

These geochemical and seismological observations indicate that the last remnants of the terrestrial magma ocean were located deep in the mantle, above the core-mantle boundary, but this remains debated both dynamically and petrologically.$^{3,20,21,22,23,24}$ Classical magma ocean solidification models, akin to those developed for the Moon, stipulate that in a cooling magma ocean, the first solids appear at the bottom of the mantle due to the intersection of liquidus and adiabat at depths, pushing the residual melt upwards.$^{25,21}$ As crystallization proceeds, solidification is expected to occur from the bottom (core mantle boundary) upwards. An alternative scenario, based on the fact that the solids are more buoyant than the melt in the deep mantle, argues that solidification takes place in the middle of the magma ocean, separating it in basal and shallow magma oceans.$^{26,3,20,23}$ The shallow magma ocean solidifies upwards quickly because of efficient cooling at the surface, while the basal magma ocean solidifies slowly, pushing the residual melt down to the core-mantle boundary.

The issue of top-down vs. bottom-up solidification is ultimately controlled by thermodynamic properties that determine: (1) where solidification takes place, i.e., intersection of liquidus with temperature, and (2) where solids and liquids accumulate, i.e., the density contrast between melt and solids. In
addition, the efficiency of solid-liquid phase separation plays a fundamental 093
role because no fractionation can occur if melts cannot gravitationally segre-
gate from solids, regardless of the petrological or geochemical nature of solids 094
and melts. Depending on the efficiency of solid-liquid phase separation,²⁷ the 095
magma ocean can freeze as a homogeneous silicate reservoir (batch crystallisa-
tion scenario) or form strongly fractionated reservoirs of distinct compositions 096
(fractional crystallisation scenario).²⁵ 097
098
099

Magma ocean fluid dynamics and multiphase flow

100

101

102

Previous fluid dynamics studies mostly focused on the issue of solid-liquid phase separation in magma oceans at high melt fraction.²⁵,²⁸,²⁹,³⁰,³¹,³² The question of whether crystals settle or are suspended by the flow have been meticulously investigated. However, in thermal evolution models, crystal fraction is not explicitly transported. In such models crystals can settle, they are deposited at the bottom of the magma ocean at a rate that is a function of the settling velocity corrected by a re-entrainment flux.³³,²⁹ Alternatively, crystal fraction can be determined assuming local thermodynamic equilibrium by comparing the temperature to the solidus and liquidus at a given depth.³¹

103

104

105

106

107

108

109

110

111

Except for a few works,³⁴ less attention has been given to the later stages
of magma ocean differentiation, where the mantle enters a solid-like dynamic
regime due to high crystal fraction. While previous fluid dynamics simulations
have investigated the issue of phase separation efficiency, they did not account
for transport of solids and liquids of distinct composition, crystal production,
and remelting. These assumptions prevent the possibility of scenarios where
(1) a clear distinction can be made between the location of crystal formation
and that of crystal accumulation, and (2) radial magma ocean structure is
controlled by the location where crystals settle, rather than the location where
they form. 121

112

113

114

115

116

117

118

119

120

121

To address these challenges, we use a novel fluid geodynamical approach
with the numerical code Bambari (see Methods and SI Section 1). Bambari
implements a multiphase flow mathematical formalism based on the averaging
method.³⁵,³⁶,²⁷,³⁷ We systematically explore a range of mantle differentiation
scenarios through two independent mechanisms: petrological differentiation
(chemical partitioning) and mechanical differentiation (solid-liquid phase
separation). Our simulations account for a density crossover between melt
and solid at mid-mantle depth (SI Section 1).³⁸,²²,²³ We focus on magma
ocean solidification dynamics when the mantle starts to behave as a solid.

122

123

124

125

126

127

128

129

130

131

We find that an iron-rich Basal Magma Ocean (BMO) always forms above
the core-mantle boundary (Figure 1), and that it corresponds to the last
residue of mantle solidification. This finding is irrespective of the depth where
liquidus and adiabat intersect²⁰,²¹ (i.e., where solidification is expected to
occur). A basal magma ocean forms at the core-mantle boundary (CMB)
through the accumulation of dense, FeO-rich solids and liquids. The FeO-rich

132

133

134

135

136

137

138
4 Solidification of Earth's mantle led inevitably to a basal magma ocean

139 solids, having low melting points, remelt, while the FeO-rich melt, being denser
140 than the solids, migrate downwards due to negative buoyancy. Both processes
141 contribute to the development of the basal magma ocean. We first consider
142 the least favorable scenario for producing a BMO,²¹ where the adiabat inter-
143 sects the liquidus at the base of the mantle. For this scenario, radial models
144 predict that solidification begins at the CMB and progresses upwards.
145

146 Magma Ocean Evolution
147

148 The initial conditions of our simulations represent a compositionally homoge-
149 neous mantle with a uniformly distributed melt fraction of 50%. Large solid
150 production occurs at the surface in the earliest stage of evolution (Fig. 1a-b),
151 because temperature in the shallow mantle drops below the solidus, regardless
152 of where liquidus and adiabat curves intersect at depth. This crystal-rich shal-
153 low layer becomes rapidly gravitationally unstable both thermally - because it
154 is colder than the underlying interior - and compositionally - because it is com-
155 posed of negatively buoyant crystals (see SI 2.1). This crystal-rich unstable
156 layer then sinks as cold downwellings. However, the solids do not accumu-
157 late in the deep mantle because they progressively remelt during their descent
158 (Fig. 1a-c).

159 With further cooling of the magma ocean (Fig. 1d-f), solids formed at
160 shallow depths gradually accumulate in the lower mantle. Our simulations
161 show that these solids, which are still fed from the surface by cold down-
162 wellings, do not remelt anymore (Fig. 1d-f). Notably, these deep-seated solids
163 are not the result of the crystallization of deep mantle melt (where the liq-
164 uidus intersects the adiabat). These solids of shallow origin carry a chemical
165 signature produced by low pressure fractionation. More precisely, these solids
166 are FeO-rich and enriched in incompatible elements. This is at odds with what
167 one would expect solely from the values of partitioning coefficients. However,
168 crystals forming near the cold thermal boundary tend to settle, propelling
169 FeO-rich melt upwards by mass conservation. This FeO-rich melt approaches
170 the cold thermal boundary, cools down by diffusion, and eventually solidifies.
171 As a result, the cold downwellings are relatively enriched in FeO as well as
172 in incompatible trace elements. For trace elements that are very incompati-
173 ble at low pressures (e.g., Sm, Nd, Lu, Hf, W), their relative abundances in
174 the residual liquid remain unchanged until the system reaches very low degree
175 of melting (Fig. SI S3). Consequently, these shallow-origin solids, formed by
176 the rapid quenching of melt with unfractionated trace element ratios, are also
177 expected to exhibit compositions that remain unfractionated relative to the
178 bulk silicate Earth. Later, we demonstrate that the accumulation of these
179 solids partially diminishes the chemical imprint of high-pressure fractionation
180 processes occurring in the deep mantle.

181 The shallow fractionation depth of the Earth's mantle is at odds with
182 all previous assessments based on 1D solidification models, highlighting the
183 importance of lateral variations in temperature and composition that cannot
184
be accounted for in 1D models. The reason why solids do not significantly 185
remelt during their descent is that (1) the average temperature in the mantle 186
(once the rheological transition is reached) is lower than the liquidus at all 187
depths, (2) the slope of the liquidus favors crystallization at larger depths, 188
and (3) the crystal-rich downwellings sink too fast to heat up and remelt by 189
thermal diffusion (see SI Section 2.3). 190

The final stage involves the accumulation of denser solids sinking as downwellings in the deep mantle (Fig. 1e-f). This forms a thermal lid above the core, leading to efficient reheating of the lowermost mantle. This in turn melts mantle fusible and FeO-rich components, that are so negatively buoyant that they accumulate at the base of the mantle, forming a basal magma ocean (Fig. 1g-h). This continuous transport of solids from the surface to the lower mantle, where they remelt and settle above the CMB, gradually enriches the BMO in FeO compared to the Bulk Silicate Earth (BSE). 191
192
193
194
195
196
197
198

Let us now discuss the case of a steeper liquidus curve²⁰ intersecting the 199
adiabat at mid-mantle depths, where a basal magma ocean is the expected 200
final stage of mantle solidification. We find that the dynamical evolution is 201
essentially the same (see SI Fig. S18), which establishes that the crossing of 202
the liquidus and adiabat does not play a major role in the style of solidification 203
of the magma ocean. A liquidus curve that intersects the adiabat in the mid- 204
mantle facilitates remelting in the lower mantle, promoting the formation of 205
a basal magma ocean. Even though the crystal fraction is higher in the mid- 206
mantle, this sluggish shell does not isolate the upper mantle from the lower 207
mantle. Downwellings that form at the surface do not accumulate at mid- 208
mantle depth as they remain thermally negatively buoyant (SI Section 3.1, SI 209
Fig. S18). They pursue their descent and remelt in the deep mantle. 210

We conclude that the formation of the Earth's basal magma ocean is 211
primarily controlled by the intrinsic high density and relatively low melt- 212
ing temperature of FeO-rich silicates, in addition to the density contrast 213
between liquid and coexisting solid silicates (i.e., solids in thermodynamic 214
equilibrium with that melt), and not by the relation between liquidus and 215
adiabat. At the end of magma ocean solidification, the resulting iron-rich 216
thermochemical structures in the lowermost mantle (Fig. 1) are geophysically 217
and mineralogically consistent with the properties of LLSVPs (Fig. 1) and 218
ULVZs.², ³, ³⁹, ⁴, ⁴⁰, ⁵ 219

The role of solid-liquid phase separation 220
221

The solidification sequence described above requires the solid-liquid phase sep- 223
aration to be faster than mantle solidification and mixing by thermochemical 224
convection. Phase separation efficiency is controlled by the dimensionless melt 225
mobility number, δ (SI Section 2.2, 2.2, Fig. S5, and Fig. S11), which is mainly 226
governed by crystal size and melt viscosity. In our simulations, while the vis- 227
cosity of partially or essentially molten regions are unrealistically large, we 228
preserve a realistic balance of the thermochemical convection velocity and 229
230
6 *Solidification of Earth's mantle led inevitably to a basal magma ocean*

231 phase separation velocity - and thus the extent of chemical differentiation - by
232 using a relatively large melt mobility number (see SI Section 2.2). To validate
233 our approach, we explored numerically the Ra-δ parameter space, where Ra
234 is the thermal Rayleigh number and δ is the melt mobility number, that cap-
235 tures the competition between chemical differentiation and convective mixing.
236 For each thermal Rayleigh number we investigated, we were able to identify a
237 critical melt mobility number above which substantial chemical differentiation
238 is observed (see SI Section 2.2).

239 By extrapolating our regime diagram to magma ocean conditions, assum-
240 ing a liquid silicate viscosity of 1 Pa s,⁴¹,⁴² melts will segregate from the
241 solids if the crystal size is larger than 0.01 µm. Even assuming a four orders
242 of magnitude higher viscosity (10⁴ Pa s) the critical crystal size above which
243 melt segregation can occur is 1 µm (SI Section 2.2), consistent with previ-
244 ous estimates.²⁵,⁴³ Since typical grain size estimates in a magma ocean are
245 around 1 mm²⁵ phase separation is expected to be efficient and dominate over
246 remixing under realistic conditions.

247 A fundamental aspect is that the fluid dynamics within or in the vicinity
248 of the top and bottom thermal boundary layer is at the heart of the chemical
249 evolution of the magma ocean. The composition of solids formed at the surface
250 of the planet is controlled by chemical fractionation between melt and low-
251 pressure mineral phases. To transport a fractionated composition into the deep
252 mantle, solids and liquids must segregate from each other in the vicinity of
253 the cold thermal boundary (SI Figures S7 and S8). Otherwise, it is the bulk
254 composition of the solid and liquid mixture that is brought to the deep mantle.
255 Scaling analysis shows that the solid-liquid segregation in the top thermal
256 boundary layer is expected to be faster than the growth rate of Rayleigh-
257 Taylor instabilities (SI section 2.3) in the Earth's magma ocean, allowing for
258 chemical shallow differentiation to occur. Although present, scaling analysis
259 indicates that this process is in fact less pronounced in our fluid dynamics
260 simulations than it would be under real conditions. The symmetric version of
261 the processes described above is responsible for the formation of the BMO.
262 FeO-rich melt migrates downwards within the hot thermal boundary layer
263 before being entrained in upwelling currents (SI Figures S9 and S10).

264
265
Geochemical consequences

266
267 We explore the geochemical consequences of our dynamical model on the
268 production and nature of primordial mantle heterogeneities, inherited from
269 magma ocean solidification. We use experimentally-determined partition coef-
270 ficients of trace elements between melts and the liquidus phase at upper mantle
271 (olivine)⁴⁴ and lower mantle (bridgmanite)⁴⁵ conditions to track a simple (no
garnet and no ferropericlase) evolutionary model of key trace-elements elemen-
272 tal ratios (Sm/Nd, Lu/Hf, Hf/W) during solidification (SI Section 1.2). At the
273 end of mantle solidification, solids in the upper mantle show superchondritic
274 Lu/Hf ratios (in yellow to red, in Fig. 2a), due to olivine crystallisation during
275
276
bottom-up crystallisation of the upper mantle. Conversely, in the lower mantle above the BMO, we observe subchondritic Lu/Hf ratios (in blue to white, in Fig. 2) that indicates bridgmanite fractionation. Furthermore, we observe that solidification and remelting at different depths produce a complex (marble cake-like) geochemical structure in the solid mantle from core to crust (Fig. 2), far from a vision inspired by geochemical two-box modelling of single enriched (residual liquid) and depleted (precipitating solids) reservoirs.⁴⁵,⁴⁶ While the extent and magnitude of trace-element ratio distribution depends dramatically on phase separation efficiency (Fig. S6 and S7), planetary scale solidification systematically generates heterogeneities at all depths (Fig. 2).

The Lu/Hf ratio is a suitable geochemical tracer to quantify the extent of low pressure (i.e., olivine) mantle solidification on a global magma ocean.⁴⁷ This is because Lu is incompatible (i.e., enriched in the melt) both in olivine and bridgmanite, while Hf is incompatible in olivine, but compatible in bridgmanite. In this simple model, solids with superchondritic (high) Lu/Hf ratios stem from olivine crystallisation and denote a shallow origin, while solids with subchondritic (low) Lu/Hf ratios originate from bridgmanite crystallisation and can only have formed at depth. The extent of mixing between these two components on the planetary scale can be seen on a Sm/Nd-Lu/Hf correlation map (Fig. S8). One should note however that the extent and amplitude of fractionation shown in Fig. 2 varies with the efficiency of solid-liquid phase separation (see Fig. S6 and S7 for more extreme cases). The amount of stirring obtained in these 2D Cartesian simulations may also be quantitatively affected by three-dimensional effects and sphericity. Unlike the lunar magma ocean, the Earth's magma ocean does not evolve as a stack of immobile layers each retaining the geochemical characteristics of the depths where they formed. Rather, primordial heterogeneities are scattered at all depths. This marks a shift in our understanding: the way lunar-like magma oceans evolve, based on the differentiation observed in magma chambers on Earth, does not necessarily apply to mantles as deep as the ones of the Earth or Mars. There, despite stirring during solidification caused by thermo-chemical convection, early-formed geochemical reservoirs are preserved, but scattered throughout the mantle.⁴⁸,⁴⁹

We analyzed the statistical distribution of geochemical heterogeneities in the solid and molten (BMO) lower mantle, and more specifically the Lu/Hf population statistics (Fig. 3). If we only account for bridgmanite fractionation (Fig. 3a) and disregard olivine fractionation in the upper mantle, then the distribution of Lu/Hf heterogeneities is centered around 0 for the solid mantle, and the BMO composition is shifted towards higher Lu/Hf ratios (70% higher than the BSE) as one would expect from bridgmanite-melt fractionation. However, when both olivine-melt and bridgmanite-melt fractionations are accounted for (Fig. 3b), the distribution of Lu/Hf heterogeneities for the solid mantle and the BMO are both centered around 0, showing that the opposite effects of low and high pressure solid-melt fractionation cancel each other out.
The role of low pressure chemical fractionation on deep mantle compo-
sition has a fundamental consequence on the quantification of the extent of
magma ocean differentiation from the standpoint of trace element ratios in
mantle rocks. The primitive upper mantle (PUM) has unfractionated refrac-
tory lithophile trace element ratios,⁵⁰ and because bridgmanite crystallisation
strongly fractionates some of these ratios,⁴⁵ this has been used to constrain
the maximum amount of bridgmanite-melt differentiation that took place in
the early Earth. On this basis, and in order to not disturb those ratios in
the PUM (within their uncertainties), no more than 8%⁵¹,⁴⁵ bridgmanite can
be crystallised alone in the magma ocean, which is inconsistent with experi-
mental melt relations and the melting phase diagrams of pyrolite. Our work
relaxes this constraint, because bridgmanite-melt partitioning is not the only
relevant process that describes the chemical differentiation between the PUM
and the lower mantle during magma ocean solidification. Olivine-melt parti-
tioning plays an important role as well, because cold downwellings form at the
surface of the planet, constantly feeding the lower mantle with a PUM-like
signature. These two effects have signatures that partially cancel each other
out, allowing for a much larger extent of lower mantle silicate differentiation
(i.e., bridgmanite-melt) to take place without affecting chondritic ratios in the
PUM (see Fig. 3).

343

Magma ocean outgassing

345

Finally, we conducted simulations to quantify the outgassing of volatile species during magma ocean solidification (see SI Section ??). We used a Lagrangian approach to track and to estimate the amount of mantle materials that would rise to the exsolution depth. We assumed that if mantle parcels remain below the exsolution depth, they remain undegassed, and retain their primordial noble gas signature and volatile content. In agreement with previous work,⁵² our results indicate that a negligible fraction of the mantle is expected to degas during magma ocean solidification beyond the rheological transition. This is consistent with the noble gas geochemical record that suggests the preservation of early formed geochemical reservoirs¹²,¹³,¹⁴ although this topic remains highly debated.⁵³

357

Conclusions

359

Our modelling was performed using the least favourable conditions to pro-
duce a BMO, from both a geochemical (full equilibrium fractionation, see SI
Section 1.2) and a petrological (adiabat intersects the liquidus at the base
of the mantle) standpoint; yet, (1) it systematically produces a basal magma
ocean in the final stage of evolution, which appears to be inevitable on Earth,
and extensible to other large terrestrial planets with Earth-like composition
and (2) the geochemical imprint of this solidification on the solid mantle is
far less marked than that predicted by two-box geochemical models,⁴⁷ due to
extensive vertical mixing during solidification. The composition of the PUM
integrates the signature of a complex mixture of shallow and deep geochemical components. Concomitantly, the signature of olivine-melt trace element fractionation is present in the lower mantle at the end of magma ocean solidification and could be preserved over geological times.

These findings call for a re-interpretation of the available geochemical record and geophysical observations to better reconstruct the Earth's thermal and chemical history from its infancy to the present-day, and more generally to better understand the diversity of terrestrial bodies.

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414
10

Solidification of Earth's mantle led inevitably to a basal magma ocean

<figure><img src="image_1.png" /><figcaption>**Fig. 1** Numerical simulations of the solidification of Earth's mantle from a mushy magma ocean state. Links to the associated videos can be found in the additional information. (a-c) **Initial stage**: After a rapid early stage of solidification, the averaged melt fraction is approximately at the rheological critical melt fraction, i.e., 50 %. (d-f) **Early stage**: Solidification of the mush occurs in concert with thermo-compositional convection at global mantle scale. Even though crystals accumulate in the deep mantle, they are formed at the surface of the planet in cold downwellings plumes. (g-i) **Late stage**: Progressive melt extraction from the cumulates differentiates the mantle. At the end of the upper magma ocean solidification, the mantle is heterogeneous. (j-l) **Final stage**: Fusible iron-rich silicates progressively pile up in the deep mantle, forming a basal magma ocean via downwards melt extraction and remelting of iron-rich solids.</figcaption></figure>

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460
<figure><img src="image_2.png" /><figcaption>**Fig. 2** Geochemical signature of magma ocean solidification for an intermediate ($\delta = 10 \times 10^{-3}$) phase separation efficiency (See Fig. S6 and S7 for more extreme cases). Snapshots at the final stage of the simulations showing the distribution of various geochemical ratios (normalised to their ratio in the BSE) across the mantle (a) Lu/Hf, (b) Hf/W, (c) Sm/Nd. Sm/Nd and Hf/W ratios are always larger than 1 in solids and lower than 1 in residual melts. Lu/Hf on the other hand behaves similarly at low pressure (olivine-melt fractionation). but becomes lower than 1 in solids and larger than 1 in residual melts at high pressures (bridgmanite-melt fractionation). The geochemical signature of both low pressure and high pressure fractionation (Lu/Hf ratio) as well as enriched solids and depleted liquids (Sm/Nd and Hf/W ratios), are preserved in the early solid mantle after magma ocean solidification. These are then expected to be further stirred by ensuing thermochemical convection in the solid mantle over billions of years.</figcaption></figure>

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

## Methods

Recent mathematical and numerical advances<sup>27,54,37,55</sup> now allow extending the multiphase flow formalism to geodynamic regimes relevant to magma ocean

505

506
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

**Fig. 3** Statistical distribution of Lu/Hf heterogeneities (shown as the logarithm of the Lu/Hf normalised to BSE) in the lower mantle (grey) and the basal magma ocean (red) at the end of magma ocean solidification. The histograms are normalized by the maximum number of occurrence. (a) The bridgmanite-melt partition coefficients of trace elements is used at all mantle depths. (b) The olivine-melt partition coefficients are used at low pressures and the bridgmanite-melt partition coefficients are used at high pressures. This second and more realist case highlights the important contribution of low pressure mineral-melt fractionation in the lower mantle, and shows that bridgmanite-melt partitioning coefficients alone are not sufficient to predict the chemical evolution of the lower mantle during magma ocean solidification.

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

solidification dynamics with vigorous convection (expressed by the thermal Rayleigh number, which is set to $10^9$ in this work) and high melt fraction (from 0 to 100%).

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

We model the solid-liquid multiphase physics using the numerical code Bambari.²⁷ Bambari implements a multiphase flow mathematical formalism based on the averaging method⁵⁶, ³⁵, ⁵⁷, ⁵⁸, ⁵⁹, ³⁶, ²⁷, ³⁷, ⁶⁰, ⁶¹ under infinite Prandtl
number and Boussinesq approximations. The original implementation$^{62,63}$ was 553 considerably optimized in recent years, in particular with the use of the numerical stencils implemented for the momentum conservation (Eq. (6)) in the Finite-Volume code StreamV.$^{64,65}$ The latter allow handling reliably large and sharp viscosity contrasts with negligible spurious pressure effects.$^{66}$ 554 555 556 557

Bambari accounts for thermochemical convection in both liquid and solid states, solid-melt phase equilibrium, mineralogical phase change and solid-melt phase separation. Convective motion in a solidifying magma ocean is 558 driven by three types of density differences. These originate from thermal 559 expansivity (due to temperature differences), compositional differences (due 560 to changes in iron content), and phase changes (due to the varying melt fraction).$^{67}$ An additional mechanism that generates motion in a partially molten 561 convective medium is the solid-liquid phase separation driven by shear deformation,$^{68}$ in addition to the density contrasts between the melt and the solid. 562 Phase separation is limited by matrix deformation, i.e., compaction, and viscous drag between the melt and the solid via Darcy law. We implemented 563 a depth-dependent density contrast between liquids and solids,$^{22,23,38}$ and 564 used a thermodynamically self-consistent compositional evolution of melts and 565 solids based on recent experimental melting phase diagrams$^{22,69}$ (SI Section 1, 566 SI Fig. S1 and SI Fig. S2). 567

The model consists of two mechanical phases: a liquid phase and a solid 573 phase. Each mechanical phase is composed of two major compositional phases 574 corresponding to a FeO-rich end-member and a MgO-rich end-member. Trace 575 elements are distributed among the mechanical phases according to their 576 solid-liquid partition coefficients. In addition to the advection and diffusion 577 of temperature, our model tracks the advection of 14 compositional fields: 578 7 liquids (active MgO-rich end-member, FeO-rich end-member, and passive 579 chemical species Hf, Sm, Nd, W, Lu) and the 7 solids counterparts. 580

In the following, we recall the governing equations solved by the code in 581 their dimensionless form. We non-dimensionalize lengths by the thickness of 582 the mantle H, time by the thermal diffusion scale $H^2/\kappa$ ($\kappa$ being the the 583 average coefficient of thermal diffusivity), velocities by $\kappa/H$, temperatures by 584 the total super-adiabatic temperature contrast across the upper and lower 585 thermal boundary layers in the mantle $\Delta T_m$ and pressure and deviatoric stress 586 using viscous pressure scale $\eta_s\kappa/H^2$, where $\eta_s$ is the viscosity of the solid 587 phase. 588

## Mass conservation equations

The mass conservation equations for the four major chemical components 591 write: 592

593 594 595 596 597 598
$$
\begin{align*}
\frac{\partial \phi_1}{\partial t} + \nabla \cdot \mathbf{v}_s \phi_1 &= \Gamma_1 \\
\frac{\partial \phi_2}{\partial t} + \nabla \cdot \mathbf{v}_s \phi_2 &= \Gamma_2 \\
\frac{\partial \phi_3}{\partial t} + \nabla \cdot \mathbf{v}_l \phi_3 &= -\Gamma_1 \\
\frac{\partial \phi_4}{\partial t} + \nabla \cdot \mathbf{v}_l \phi_4 &= -\Gamma_2
\end{align*}
\tag{1}
$$

where $t$ is the time, $\mathbf{v}_l$ is the velocity of the liquid phase, $\mathbf{v}_s$ is the velocity of the solid phase and $\Gamma_i$ are associated with the rate of phase change. The subscripts 1 and 3 refer to the MgO end-members in the solid and liquid phase, respectively. The subscripts 2 and 4 refer to the FeO end-members in the solid and liquid phase, respectively. The densities of the four chemical end-members are:

$$
\begin{align}
\rho_1 &= \rho_0(1 - \alpha T) + \frac{1}{2}\Delta\rho_1 - \frac{1}{2}\Delta\rho_2 \notag \\
\rho_2 &= \rho_0(1 - \alpha T) + \frac{1}{2}\Delta\rho_1 + \frac{1}{2}\Delta\rho_2 \notag \\
\rho_3 &= \rho_0(1 - \alpha T) - \frac{1}{2}\Delta\rho_1 - \frac{1}{2}\Delta\rho_2 \notag \\
\rho_4 &= \rho_0(1 - \alpha T) - \frac{1}{2}\Delta\rho_1 + \frac{1}{2}\Delta\rho_2, \tag{2}
\end{align}
$$

where $\rho_0$ is the reference density, the thermal expansivity $\alpha$ is considered constant in the solid and liquid phase, and $T$ is the temperature. The isochemical density contrast between the liquid and solid, is equal to $\Delta\rho_1 = \rho_1 - \rho_3 = \rho_2 - \rho_4$. Similarly, the chemical density contrast between the dense and light end-members is $\Delta\rho_2 = \rho_2 - \rho_1 = \rho_4 - \rho_3$. We parameterized these density contrasts to fit a self-consistent thermodynamic model$^{22}$ (see SI Fig. S1).

632 Energy conservation equation

We assume local thermal equilibrium between all the phases, i.e., the temperatures at a given location in all the phases are the same ($T_1(x, z, t) = T_2(x, z, t) = T_3(x, z, t) = T_4(x, z, t)$). This leads to a single equation for the conservation of energy:

$$
\partial_t T + \bar{\nabla} \cdot \nabla T = \nabla^2 T + (\Gamma_1 + \Gamma_2)S_t, \quad (3)
$$

where $\bar{\nabla}$ is the velocity of the solid-liquid mixture, $S_t$ is the dimensionless Stefan number that writes,

$$
S_t = \frac{L}{C_p \Delta T_m}, \quad (4)
$$
where $C_p$ is the thermal capacity at constant pressure for the mixture of the four components, and $L$ is the latent heat release upon phase change.

The velocity of the solid-liquid mixture is related to the solid and liquid velocities via the following equation:

$$ \bar{\mathbf{v}} = \phi \mathbf{v}_s + (1 - \phi) \mathbf{v}_l. \qquad (5) $$

## Momentum conservation equations

Since we account for two mechanical phases, two momentum conservation equations are required to describe the two-phase flow mechanical interactions. The following Stokes equation describes the momentum conservation for the averaged solid-liquid mixture:

$$ -\nabla \Pi + \nabla \cdot \bar{\tau} - Ra T \mathbf{g} - Rp \phi \mathbf{g} + \frac{1}{2} Rc (\phi_2 + \phi_4 - \phi_1 - \phi_3) = 0, \quad (6) $$

where $\nabla \Pi$ is the dynamic pressure gradient ($\nabla \Pi = \nabla P - \rho_0 \mathbf{g}$, where $P$ is the total pressure, $\mathbf{g}$ is the gravity vector, and $\rho$ is the reference density), $\bar{\tau}$ is the deviatoric stress tensor for the solid-liquid mixture, the thermal Rayleigh number, $Ra$, and the phase and compositional Rayleigh numbers, $Rp$ and $Rc$ are respectively:

$$ Ra = \frac{\rho_0 \alpha \Delta T_m g H^3}{\kappa \eta_s}, \qquad 667 $$

$$ Rp = \frac{\Delta \rho_1 g H^3}{\kappa \eta_s}, \qquad (7) $$

$$ Rc = \frac{\Delta \rho_2 g H^3}{\kappa \eta_s}, \qquad 671 $$

where $g$ is the constant acceleration of gravity (the magnitude of $\mathbf{g}$). We define the stress tensor of the solid-liquid mixture:

$$ \bar{\tau} = (1 - \phi) \bar{\tau}_s + \phi \bar{\tau}_l = \eta(\phi) \left( \nabla \bar{\mathbf{v}} + [\nabla \bar{\mathbf{v}}]^t \right), \qquad (8) $$

where $\eta(\phi)$ is the viscosity of fluid mixture that varies from $\eta(\phi = 0) = \eta_s$ to $\eta(\phi = 1) = \eta_l$. We use

$$ \eta(\phi) = \eta_s \times 10^{\left(1+\tanh\left(\frac{\phi-0.5}{0.1}\right)\right) \frac{\eta_l}{2\eta_s}}. \qquad (9) $$

For the second momentum conservation equation, we use a retro-action relationship³⁶,²⁷ that describes $\Delta \mathbf{v} = \mathbf{v}_l - \mathbf{v}_s$, the velocity difference between the two phases:

$$ \phi \Delta \mathbf{v} = -\phi^2(1-\phi)\delta^2\left(\zeta \nabla \left[\frac{(1-\phi)}{\phi} \nabla \cdot \phi \Delta \mathbf{v}\right] + X \mathbf{g} + \nabla \cdot \bar{\tau}\right), \qquad (10) $$
691 where

692 $X = Rp + \frac{1}{2} \left( \frac{\phi_3 - \phi_4}{\phi} + \frac{\phi_2 - \phi_1}{1 - \phi} \right) Rc$, (11)

693 $\delta$ is the melt mobility number, and $\zeta$ is the dimensionless compaction viscosity.

694 The melt mobility number and dimensionless compaction viscosity writes:

696 $\delta = \sqrt{\frac{r_0^2 \eta_s}{C_0 H^2 \eta_l}}$, (12)

698 $\zeta = \frac{\eta_c}{\eta_s}$,

700 with $r_0$ is the crystal size, $C_0$ is a constant whose value is discussed in the Supplement, and $\eta_c$ is a constant associated with the compaction viscosity, sometimes referred to as bulk viscosity, that describes the viscous resistance of a solid matrix to compact (i.e., close the porous space)$^{35,36,70,71}$. Here, we follow ref.$^{36}$ where the compaction viscosity, $\eta_{\text{bulk}}$, is defined as $\eta_{\text{bulk}} = \eta_c/\phi$ (the compaction viscosity must be infinite when the melt fraction is 0).

702 The above equations are discretized using conservative finite-difference schemes. The corresponding set of Equations (6) and (10) involves the inversion of two sparse matrices, which are the most time-consuming operations. In Bambari, these matrix inversions are performed using the direct PARDISO library that is parallelized with OpenMP directives.$^{72,73}$ The numerical implementation of the two phases flow physics has been benchmarked against 1D analytical solutions (e.g., Section 3.4 in$^{63}$). The new additions to the numerical implementation were successfully benchmarked against numerical an analytical solutions.$^{66,65}$

718 Phase Change and Chemical Partitioning
719 Methods related to the modeling of phase change and chemical partitioning can be found in the supplementary materials.

722
723
724
725
726
727
728
729
730
731
732
733
734
735
736

722
723
724
725
726
727
728
729
730
731
732
733
734
735
736

722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
1 Dziewonski, A. M., Lekic, V. & Romanowicz, B. A. Mantle anchor structure: an argument for bottom up tectonics. *Earth and Planetary Science Letters* **299** (1-2), 69–79 (2010).

722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
2 Tolstikhin, I., Kramers, J. D. & Hofmann, A. A chemical earth model with whole mantle convection: the importance of a core–mantle boundary layer (d’) and its early formation. *Chemical geology* **226** (3-4), 79–99 (2006).

722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
3 Labrosse, S., Hernlund, J. & Coltice, N. A crystallizing dense magma ocean at the base of the Earth’s mantle. *Nature* **450** (7171), 866–869 (2007). URL http://www.nature.com/nature/journal/v450/n7171/full/nature06355.html. https://doi.org/10.1038/nature06355.
⁴ Lee, C.-T. A. et al. Upside-down differentiation and generation of a ‘primordial’lower mantle. *Nature* **463** (7283), 930–933 (2010). 737

738

⁵ Carlson, R. W. et al. How did early earth become our modern world? *Annual Review of Earth and Planetary Sciences* **42**, 151–178 (2014). 740

741

⁶ Ballmer, M. D., Houser, C., Hernlund, J. W., Wentzcovitch, R. M. & Hirose, K. Persistence of strong silica-enriched domains in the earth’s lower mantle. *Nature Geoscience* **10** (3), 236–240 (2017). 742

745

746

⁷ Willbold, M., Elliott, T. & Moorbath, S. The tungsten isotopic composition of the earth’s mantle before the terminal bombardment. *Nature* **477** (7363), 195–198 (2011). 747

748

749

⁸ Touboul, M., Puchtel, I. S. & Walker, R. J. 182w evidence for long-term preservation of early mantle differentiation products. *Science* **335** (6072), 1065–1069 (2012). 750

752

⁹ Rizo, H., Boyet, M., Blichert-Toft, J. & Rosing, M. T. Early mantle dynamics inferred from 142nd variations in archean rocks from southwest greenland. *Earth and Planetary Science Letters* **377–378**, 324 – 335 (2013). URL http://www.sciencedirect.com/science/article/pii/S0012821X13003877. https://doi.org/http://dx.doi.org/10.1016/j.epsl.2013.07.012. 754

759

¹⁰ Rizo, H. et al. Preservation of earth-forming events in the tungsten isotopic composition of modern flood basalts. *Science* **352** (6287), 809–812 (2016). 761

762

¹¹ Morino, P., Caro, G. & Reisberg, L. Differentiation mechanisms of the early hadean mantle: insights from combined 176hf-142,143 nd signatures of archean rocks from the saglek block. *Geochimica et Cosmochimica Acta* **240**, 43–63 (2018). 764

765

766

¹² Kurz, M. D., Curtice, J., Fornari, D., Geist, D. & Moreira, M. Primitive neon from the center of the galápagos hotspot. *Earth and Planetary Science Letters* **286** (1-2), 23–34 (2009). 769

770

771

¹³ Mukhopadhyay, S. Early differentiation and volatile accretion recorded in deep-mantle neon and xenon. *Nature* **486** (7401), 101–104 (2012). 772

773

¹⁴ Caracausi, A., Avice, G., Burnard, P. G., Füri, E. & Marty, B. Chondritic xenon in the earth’s mantle. *Nature* **533** (7601), 82–85 (2016). 775

776

777

¹⁵ Zhang, N., Zhong, S., Leng, W. & Li, Z.-X. A model for the evolution of the earth’s mantle structure since the early paleozoic. *Journal of Geophysical Research: Solid Earth* **115** (B6) (2010). 778

779

780

781

782
783 16 Thorne, M. S., Garnero, E. J. & Grand, S. P. Geographic correlation between hot spots and deep mantle lateral shear-wave velocity gradients. *Physics of the Earth and Planetary Interiors* **146** (1-2), 47-63 (2004).

786 17 Burke, K., Steinberger, B., Torsvik, T. H. & Smethurst, M. A. Plume generation zones at the margins of large low shear velocity provinces on the core-mantle boundary. *Earth and Planetary Science Letters* **265** (1-2), 49-60 (2008).

791 18 Torsvik, T. H., Burke, K., Steinberger, B., Webb, S. J. & Ashwal, L. D. Diamonds sampled by plumes from the core-mantle boundary. *Nature* **466** (7304), 352-355 (2010).

795 19 French, S. W. & Romanowicz, B. Broad plumes rooted at the base of the earth's mantle beneath major hotspots. *Nature* **525** (7567), 95-99 (2015).

798 20 Fiquet, G. et al. Melting of peridotite to 140 gigapascals. *Science* **329** (5998), 1516-1518 (2010). URL http://www.sciencemag.org/content/329/5998/1516full.pdf.

803 21 Andrault, D. et al. Solidus and liquidus profiles of chondritic mantle: Implication for melting of the earth across its history. *Earth Planet. Sci. Lett* **304** (1-2), 251-259 (2011). URL http://www.sciencedirect.com/science/article/pii/S0012821X1100077X. https://doi.org/10.1016/j.epsl.2011.02.006.

808 22 Boukaré, C.-E., Ricard, Y. & Fiquet, G. Thermodynamics of the MgO-FeO-SiO₂ system up to 140 GPa: Application to the crystallization of Earth's magma ocean. *J Geophys. Res.: Solid Earth* **120** (9), 6085-6101 (2015). URL http://dx.doi.org/10.1002/2015JB011929. https://doi.org/10.1002/2015JB011929, 2015JB011929.

814 23 Caracas, R., Hirose, K., Nomura, R. & Ballmer, M. D. Melt-crystal density crossover in a deep magma ocean. *Earth and Planetary Science Letters* **516**, 202-211 (2019).

818 24 Miyazaki, Y. & Korenaga, J. On the timescale of magma ocean solidification and its chemical consequences: 2. compositional differentiation under crystal accumulation and matrix compaction. *Journal of Geophysical Research: Solid Earth* **124** (4), 3399-3419 (2019).

823 25 Solomatov, V. S. in *9.04 - magma oceans and primordial mantle differentiation* (ed.Schubert, G.) *Treatise on Geophysics* 91-119 (Elsevier, Amsterdam, 2007). URL http://www.sciencedirect.com/science/article/pii/B9780444527486001413.
26 Stixrude, L. & Karki, B. Structure and freezing of MgSiO₃ liquid in Earth's lower mantle. *Science* **310** (5746), 297–299 (2005). URL http://www.sciencemag.org/content/310/5746/297full.pdf.

27 Boukaré, C.-E. & Ricard, Y. Modeling phase separation and phase change for magma ocean solidification dynamics. *Geochemistry, Geophysics, Geosystems* **18** (9), 3385–3404 (2017).

28 Tonks, W. & Melosh, H. The physics of crystal settling and suspension in a turbulent magma ocean. *Origin of the Earth* **1**, 151–174 (1990).

29 Lavorel, G. & Le Bars, M. Sedimentation of particles in a vigorously convecting fluid. *Physical Review E* **80** (4), 046324 (2009).

30 Suckale, J., Elkins-Tanton, L. T. & Sethian, J. A. Crystals stirred up: 2. numerical insights into the formation of the earliest crust on the moon. *Journal of Geophysical Research: Planets* **117** (E8), n/a–n/a (2012). URL http://dx.doi.org/10.1029/2012JE004067. https://doi.org/10.1029/2012JE004067, e08005.

31 Bower, D. J., Sanan, P. & Wolf, A. S. Numerical solution of a non-linear conservation law applicable to the interior dynamics of partially molten planets. *Physics of the earth and planetary interiors* **274**, 49–62 (2018).

32 Maas, C. & Hansen, U. Dynamics of a terrestrial magma ocean under planetary rotation: A study in spherical geometry. *Earth and Planetary Science Letters* **513**, 81–94 (2019).

33 Martin, D. & Nokes, R. Crystal settling in a vigorously converting magma chamber. *Nature* **332** (6164), 534–536 (1988).

34 Hier-Majumder, S. & Hirschmann, M. M. The origin of volatiles in the earth's mantle. *Geochemistry, Geophysics, Geosystems* **18** (8), 3078–3092 (2017).

35 Mckenzie, D. The generation and compaction of partially molten rock. *Journal of Petrology* **25** (3), 713–765 (1984).

36 Bercovici, D., Ricard, Y. & Schubert, G. A two-phase model for compaction and damage: 1. general theory. *Journal of Geophysical Research: Solid Earth* **106** (B5), 8887–8906 (2001). URL http://dx.doi.org/10.1029/2000JB900430. https://doi.org/10.1029/2000JB900430.

37 Keller, T. & Suckale, J. A continuum model of multi-phase reactive transport in igneous systems. *Geophysical Journal International* **219** (1), 185–222 (2018).
875 (2019).

876

877

878

879

880

881

882

883

884

885

886

887

888

889

890

891

892

893

894

895

896

897

898

899

900

901

902

903

904

905

906

907

908

909

910

911

912

913

914

915

916

917

918

919

920

38 Funamori, N. & Sato, T. Density contrast between silicate melts and crystals in the deep mantle: An integrated view based on static-compression data. *Earth Planet. Sci. Lett* **295**, 435–440 (2010). https://doi.org/10.1016/j.epsl.2010.04.021.

39 Wicks, J., Jackson, J. & Sturhahn, W. Very low sound velocities in iron-rich (mg, fe) o: Implications for the core-mantle boundary region. *Geophysical Research Letters* **37** (15) (2010).

40 Bower, D. J., Wicks, J. K., Gurnis, M. & Jackson, J. M. A geodynamic and mineral physics model of a solid-state ultralow-velocity zone. *Earth and Planetary Science Letters* **303** (3-4), 193–202 (2011).

41 Karki, B. B. & Stixrude, L. P. Viscosity of mgsio3 liquid at earth's mantle conditions: Implications for an early magma ocean. *Science* **328** (5979), 740–742 (2010).

42 Dygert, N., Lin, J.-F., Marshall, E. W., Kono, Y. & Gardner, J. E. A low viscosity lunar magma ocean forms a stratified anorthitic flotation crust with mafic poor and rich units. *Geophysical Research Letters* **44** (22), 11,282–11,291 (2017). URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017GL075703. https://doi.org/https://doi.org/10.1002/2017GL075703, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2017GL075703.

43 Elkins-Tanton, L. T. Magma oceans in the inner solar system. *Annu. Rev. Earth Planet. Sci.* **40**, 113–139 (2012). https://doi.org/{10.1146/annurev-earth-042711-105503}.

44 Kennedy, A., Lofgren, G. & Wasserburg, G. An experimental study of trace element partitioning between olivine, orthopyroxene and melt in chondrules: equilibrium values and kinetic effects. *Earth and Planetary Science Letters* **115** (1-4), 177–195 (1993).

45 Corgne, A., Liebske, C., Wood, B. J., Rubie, D. C. & Frost, D. J. Silicate perovskite-melt partitioning of trace elements and geochemical signature of a deep perovskitic reservoir. *Geochimica et Cosmochimica Acta* **69** (2), 485–496 (2005).

46 Rizo, H. et al. The elusive Hadean enriched reservoir revealed by Nd-142 deficits in Isua Archaean rocks. *Nature* **491** (7422), 96–U109 (2012). https://doi.org/{10.1038/nature11565}.

47 Caro, G., Bourdon, B., Wood, B. J. & Corgne, A. Trace-element fractionation in hadean mantle generated by melt segregation from a magma ocean.
Nature 436 (7048), 246–249 (2005). 921

48 Maurice, M. et al. Onset of solid-state mantle convection and mixing during magma ocean solidification. *J. Geophys. Res.: Planets* 122 (3), 577–598 (2017). URL http://dx.doi.org/10.1002/2016JE005250. https://doi.org/10.1002/2016JE005250, 2016JE005250. 922
923
924
925
926

49 Morison, A., Labrosse, S., Deguen, R. & Alboussière, T. Timescale of overturn in a magma ocean cumulate. *Earth and Planetary Science Letters* 516, 25–36 (2019). 927
928
929
930

50 Mcdonough, W. & Sun, S. The composition of the earth. *Chemical Geology* 120 (3-4), 223–253 (1995). https://doi.org/10.1016/0009-2541(94)00140-4. 931
932
933

51 Herzberg, C. & O'hara, M. Origin of mantle peridotite and komatiite by partial melting. *Geophysical Research Letters* 12 (9), 541–544 (1985). 934
935
936

52 Salvador, A. & Samuel, H. Convective outgassing efficiency in planetary magma oceans: insights from computational fluid dynamics. *Icarus* 390, 115265 (2023). 937
938
939
940

53 Parai, R. A dry ancient plume mantle from noble gas isotopes. *Proceedings of the National Academy of Sciences of the United States of America* 119 (29), 1–9 (2022). https://doi.org/10.1073/pnas.2201815119. 941
942
943

54 Oliveira, B., Afonso, J. C., Zlotnik, S. & Diez, P. Numerical modelling of multiphase multicomponent reactive transport in the earth's interior. *Geophysical Journal International* 212 (1), 345–388 (2018). 945
946
947

55 Wong, Y.-Q. & Keller, T. A unified numerical model for two-phase porous, mush and suspension flow dynamics in magmatic systems. *Geophysical Journal International* 233 (2), 769–795 (2023). 948
949
950
951

56 Drew, D. A. Averaged field equations for two-phase media. *Studies in Applied Mathematics* 50 (2), 133–166 (1971). 952
953
954

57 Ribe, N. M. Theory of melt segregation—a review. *Journal of volcanology and geothermal research* 33 (4), 241–253 (1987). 955
956
957

58 Katz, R. F. Magma dynamics with the enthalpy method: Benchmark solutions and magmatic focusing at mid-ocean ridges. *Journal of Petrology* 49 (12), 2099–2121 (2008). 958
959
960
961

59 Rudge, J. F., Bercovici, D. & Spiegelman, M. Disequilibrium melting of a two phase multicomponent mantle. *Geophysical Journal International* 184 (2), 699–718 (2011). 962
963
964
965
966
967 60 Katz, R. F., Jones, D. W. R., Rudge, J. F. & Keller, T. Physics of melt extraction from the mantle: Speed and style. *Annual Review of Earth and Planetary Sciences* **50** (1), 507–540 (2022).

968
969
970
971 Katz, R. F. *The dynamics of partially molten rock* (Princeton University Press, 2022).
972

973
974
975
976
977 62 Šrámek, O., Ricard, Y. & Bercovici, D. Simultaneous melting and compaction in deformable two-phase media. *Geophysical Journal International* **168** (3), 964–982 (2007).

978
979
980
981 63 Šrámek, O. *Modèle d'écoulement biphasé en sciences de la Terre: fusion partielle, compaction et différenciation*. Ph.D. thesis, Université de Lyon - Ecole Normale Supérieure, Lyon (2007).
982
983

984
985
986
987
988 64 Samuel, H. Time-domain parallelization for computational geodynamics. *g-cubed* (2012). Doi:10.1029/2011GC003905.

989
990
991
992
993
994 Samuel, H. A deformable particle-in-cell method for advective transport in geodynamic modelling. *Geophysical Journal International* **214**, 1744–1773 (2018). Doi:10.1093/gji/ggy231.

995
996
997
998
999
1000 66 Samuel, H. & Evonuk, M. Modeling advection in geophysical flows with particle level sets. *g-cubed* **11** (Q08020) (2010). Doi 10.1029/2010GC003081.

1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012 67 Pusok, A. E., Katz, R. F., May, D. A. & Li, Y. Chemical heterogeneity, convection and asymmetry beneath mid-ocean ridges. *Geophysical Journal International* **231** (3), 2055–2078 (2022). URL https://doi.org/10.1093/gji/ggac309. https://doi.org/10.1093/gji/ggac309, https://academic.oup.com/gji/article-pdf/231/3/2055/45807427/ggac309.pdf.

1009
1011
1012
1008
1007
1006
1005
1004
1003
1002
68 Rabinowicz, M. & Vigneresse, J.-L. Melt segregation under compaction and shear channeling: Application to granitic magma segregation in a continental crust. *Journal of Geophysical Research: Solid Earth* **109** (B4) (2004)

1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
70 Rudge, J. F. The viscosities of partially molten materials undergoing diffusion creep. *Journal of Geophysical Research: Solid Earth* **123** (12), 10–534 (2018).
71 Connolly, J. A. & Schmidt, M. W. Viscosity of crystal-mushes and implications for compaction-driven fluid flow. *Journal of Geophysical Research: Solid Earth* **127** (9), e2022JB024743 (2022).

72 Alappat, C. et al. A recursive algebraic coloring technique for hardware-efficient symmetric sparse matrix-vector multiplication. *ACM Trans. Parallel Comput.* **7** (3) (2020). URL https://doi.org/10.1145/3399732. https://doi.org/10.1145/3399732.

73 Bollhöfer, M., Schenk, O., Janalik, R., Hamm, S. & Gullapalli, K. State-of-the-art sparse direct solvers. *Parallel Algorithms in Computational Science and Engineering* 3–33 (2020). URL https://doi.org/10.1007/978-3-030-43736-7_1. https://doi.org/10.1007/978-3-030-43736-7_1.

## Data and Code Availability

The raw files used to produce the figures of this paper are available at https://doi.org/10.18715/IPGP.2024.m42039nd. The code Bambari is available upon reasonable request.

## Acknowledgments

This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 101019965—SEPtiM). Parts of this work were supported by the UnivEarthS Labex program at Université de Paris and IPGP (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02), and Natural Sciences and Engineering Research Council of Canada (RGPIN-2024-06174). We thank Catherine Chauvel for fruitful discussions and Kang Wei Lim for his thorough proofreading of our manuscript. Numerical computations were performed on the S-CAPAD/DANTE platform, IPGP, France.

## Author contributions

The authors declare that there is no financial or non-financial competing interest. Questions regarding reprints and permissions should be addressed to Charles-Édouard Boukaré (boukare@yorku.ca).

**Author 1: C.-E. Boukaré** - Conceived and designed the research - Implemented the numerical code - Ran the numerical simulations - Analysed the data - Wrote the paper. **Author 2: J. Badro** Conceived and designed the research - Designed the implementation of the petrological models in the fluid dynamics simulations - Analysed the data - Wrote the paper. **Author 3: H. Samuel** Code implementation and optimization - Design of the numerical experiments - Benchmarked the numerical implementation - Analysed the data - Wrote the paper.
1059 Additional information

1060

1061 Descriptions of the methods can be found in the Supplementary Information.

1062 Videos corresponding to Figures 1 and 2 can be found via the following links :

1063 • Figure 1 : Temperature and FeO fields,

1064 • Figure 2 : Trace Elements fields.

1065

1066

1067

1068

1069

1070

1071

1072

1073

1074

1075

1076

1077

1078

1079

1080

1081

1082

1083

1084

1085

1086

1087

1088

1089

1090

1091

1092

1093

1094

1095

1096

1097

1098

1099

1100

1101

1102

1103

1104
**Supplementary Information:**

Solidification of Earth's mantle led inevitably
to a basal magma ocean

Charles-Édouard Boukaré<fnref n="1" />,<fnref n="2" />, James Badro<fnref n="1" />† and Henri Samuel<fnref n="1" />†

Contributing authors: boukare@ipgp.fr; samuel@ipgp.fr;
badro@ipgp.fr;

†These authors contributed equally to this work.

We use the label MT to refer to equations introduced in the Methods section of the *Main Text*.

# 1 Phase change

We assume thermodynamic equilibrium between solid and liquid at every time step, everywhere in the model domain. Phase change is therefore not limited by reaction kinetics but only by phase diagrams (which describe the solid-liquid equilibrium) and energy conservation. Phase change is treated differently for major and trace elements. For major elements, phase change is constrained by a parameterized phase diagram (Fig. S2) and an energy conservation equation. For simplicity, we neglect the latent heat of phase change for trace elements. Trace elements are simply distributed among the phases according to their respective partitioning coefficients.

<footnote n="1">¹Université Paris Cité, Institut de Physique du Globe de Paris,
CNRS, 1 rue Jussieu, 75005 Paris, France.</footnote>

<footnote n="2">²York University, Department of Physics and Astronomy, 4700 Keele Street, M3J 1P3 Toronto, Canada.</footnote>
2 Solidification of Earth's mantle led inevitably to a basal magma ocean

<figure><img src="image_43.png" /><figcaption>Fig. S1 Density contrast parameterization. (a) Density of solid and liquid silicate as function of pressure for various composition. (b) Density difference between solid and liquid. The thick grey lines correspond to the exact difference between the density curves shown on the left. The color lines depict the parameterization we used in this study. The iso-chemical density contrast is shown in blue ($\Delta\rho_1$ in MT Eq. (2)) and capture to the effect of highest compressibility of silicate melt compared to their solid couterparts. The chemical density contrast correspond to the density contrast between iron-rich and iron-free silicates ($\Delta\rho_2$ in MT Eq. (2)).</figcaption></figure>

## 1.1 FeO partitioning - Equilibrium thermodynamics

We opted for an idealized binary loop to capture mantle melting. This simple approach allows us to capture the main influence of FeO partitioning on magma ocean evolution (Fig. S2):

* The first-order influence of FeO on both solid and liquid densities, which controls the gravitational stability of FeO-rich layers.

* The first-order role of FeO on the melting temperature, as FeO acts has an anti-freeze effect by shifting the melting curves towards higher values.

* FeO is slightly more compatible in the liquid phase than in the crystals (see binary loop in Fig. S2).

After updating of the temperature field using MT Eq. (3), i.e., prior to phase change, the dimensionless energy conservation in each computational cell (i.e., volume control) of the numerical domain writes:

$$T_f - T_i = (\Gamma_1 + \Gamma_2)S_t, \qquad (1)$$

where $T_f$ is the final temperature, $T_i$ is the temperature before phase change (see MT Eq. (3)). Mass conservation of the melt fraction ($\phi = \phi_3 + \phi_4$) writes,

$$\phi_f - \phi_i = \Gamma, \qquad (2)$$
<figure><img src="image_44.png" /><figcaption>**Fig. S2** Idealized binary loop model used to model solid-liquid melting relations in the Earth's mantle. This parameterization was obtained by fitting diamond anvil cell experiments [1] and self-consistent thermodynamic calculations [2]. This phase diagram has a primary role in the early Earth's mantle evolution. First, it controls the FeO chemical partitioning between solids and liquid, which has a first-order impact on the buoyancy of iron-rich silicates. Second, it controls the effect of FeO content on silicate melting temperature. This plays an important role on the formation of a basal magma ocean by remelting iron-rich solids at the CMB.</figcaption></figure>

where $\phi_f$ is the final melt fraction, $\phi_i$ is the melt fraction before phase change, and $\Gamma = (\Gamma_1 + \Gamma_2)$. For each time step, after solving for the conservation of energy (i.e., after one advection-diffusion time-step for temperature), the solid-liquid mixture is potentially out of equilibrium. We therefore consider three cases to treat this disequilibrium:

1) *The temperature is above the liquidus but the liquid fraction is lower than 1.* If the temperature is hot enough such that it remains above the liquidus after melting of the residual solids, melt fraction is set to 1. The melting temperature is updated using Eq. (1) according to the associated $\Gamma$.

2) *The temperature is below the solidus but the liquid fraction is larger than 0.* If the temperature is cold enough to remain below the solidus after the solidification of the residual melt, melt fraction is set to 0. The temperature is updated using Eq. (1) according to the associated $\Gamma$.

3) All other cases require a special treatment. The challenge is to determine $\Gamma_1$, $\Gamma_2$ and $T_f$ such that $T_f$, the final melt fraction as well as the relative amount of FeO in the melt and the solid correspond to a thermodynamic equilibrium.

We describe below our approach to handle the case 3 mentioned above, which is a variant of the approach proposed in [3].
4 Solidification of Earth's mantle led inevitably to a basal magma ocean

Because the final temperature and the melt fraction must be at equilibrium, we have to solve the following system,

$$ \begin{cases} T_i + \Gamma S_t = T_{eq} \\ \phi_i + \Gamma = \phi_{eq} \end{cases} \qquad (3) $$

where $T_{eq}$ is the equilibrium temperature between the solidus and the liquidus and $\phi_{eq}$ is the equilibrium melt fraction. Both quantities are constrained by the phase diagram. We need additional equations to close this set of equations i.e., we need to solve the system (3) for only one unknown in addition to $\Gamma$).

The branches that describe the liquidus and solidus temperatures in the simplified phase diagram we use (Fig. S2) are parameterized as follows:

$$ T_L = T_{\text{FeO}} + (T_{\text{MgO}} - T_{\text{FeO}}) (1 - x_l^m), \qquad (4) $$

$$ T_S = T_{\text{FeO}} + (T_{\text{MgO}} - T_{\text{FeO}}) (1 - x_s)^n \qquad (5) $$

where $T_L$ is the liquidus temperature, $T_S$ is the solidus temperature, $T_{\text{MgO}}$ and $T_{\text{FeO}}$ are the temperature of the pure end-members, $x_l$ is the composition of the liquid in FeO, $x_s$ is the composition of the solid in FeO, $m$ and $n$ are constant that we used to adjust the shape of the loop. At thermodynamic equilibrium we have:

$$ T_L = T_S = T_{\text{eq}}. \qquad (6) $$

This yields an equation that relates the liquid and solid compositions:

$$ x_l = [1 - (1 - x_s)^n]^{1/m}. \qquad (7) $$

By the lever rule, we have,

$$ \phi_{\text{eq}} = \frac{x_l - x_{\text{bulk}}}{x_l - x_s}, \qquad (8) $$

where $x_{\text{bulk}}$ is the bulk composition of the solid-liquid mixture in FeO. Using Equations (4), (7) and (8), we can rewrite the system of two equations (3) with only two unknown variables, $\Gamma$ and $x_s$,

$$ \begin{cases} T_i + \Gamma S_t = T_S(x_s) \\ \phi_i + \Gamma = \phi_{\text{eq}}(x_l(x_s), x_s) \end{cases} \qquad (9) $$

Due to the non-linearity of Eq. (4) and Eq. (7), we solve this set of equations iteratively by minimizing the following function for $x_s$ and $\Gamma$ using a Newton-Raphson method in two dimensions,

$$ G(x_s, \Gamma) = [T_i + \Gamma S_t - T_S(x_s)]^2 + [\phi_i + \Gamma - \phi_{\text{eq}}(x_l(x_s), x_s)]^2 \qquad (10) $$

$\Gamma_1$ and $\Gamma_2$ can be retrieved a posteriori when the final melt fraction, liquid and solid compositions are known.
While we focused on the limiting case where liquid and solid phases remain at thermodynamic equilibrium, non-equilibrium processes may be important in quantifying the extent of mantle differentiation. Indeed, depending on the efficiency of solid-liquid phase separation, melt migration may occur on a much shorter timescale than the rate of solidification. Several mathematical models have been proposed [4–6] to evaluate the effect of non-equilibrium thermodynamics in the context of geophysical multiphase flows. Further work is required to apply these theoretical approaches to the context of a vigorously convecting magma ocean with larger initial melt fractions. However, it is important to note that non-equilibrium phase change will only increase the allowable extent of silicate-silicate fractionation that is consistent with an non-fractionated PUM; If solids do not have time to re-equilibrate, they will inevitably tend to further retain the geochemical signature of the depth where they originally form.

## 1.2 Trace elements partitioning

At each time step, in every grid cell, we compute the total amount of a given trace element specie i. We distribute the trace element i between the liquid and the solid phase according to its partitioning coefficient, $D_i$ defined as:

$$D_i = \frac{x_i^l}{x_i^s}, \qquad (11)$$

where $x_i^l$ and $x_i^s$ are the liquid and solid composition in trace element i. We use the experimentally-determined partitioning coefficient of [7] and [8]. Note that the trace elements considered here are purely passive, that is they do not affect any solid or fluid properties.

By my mass balance, we can deduce the liquid and solid composition as a function of the melt fraction:

$$x_i^s = \frac{x_0}{1 + \phi(D_i - 1)}, \qquad (12)$$

where $x_0$ is the bulk concentration of the trace element i and $\phi$ is the melt fraction. Figure S3 displays the evolution of liquid and solid composition for the lower mantle and the upper mantle. Note that these plots do not involve any fluid dynamics calculation.

We distribute the trace elements according to Eq. (12) once the equilibrium melt fraction has been reached via the procedure described in Section 1.1. In Fig. S4, we compare the analytical prediction shown in Fig. S3 with compositions obtained using our fluid dynamics simulations. This allows us to benchmark the numerical implementation of trace elements partitioning. The simulations correspond to a case where phase separation is efficient ($\delta = 20 \times 10^{-3}$) such that it maximizes the effect of solid-liquid fractionation. Fig. S4 shows the frequency of occurrence of a given composition. This
6 Solidification of Earth's mantle led inevitably to a basal magma ocean

<figure><img src="image_45.png" /><figcaption>**Fig. S3** Composition of the melt and solids in Sm/Nd (magenta), Lu/Hf (green) and Hf/W (cyan) as function of melt fraction. This plot does not involve any fluid dynamics calculation and severs as a consistency check for the behaviour of our numerical simulations (see Fig. S4). The trace elements ratio are normalized by their ratio in the bulk silicate Earth (BSE). (a and c) In the upper mantle, trace elements are so incompatible in the olivine that the composition departs from the BSE composition only towards the end of solidification. Solids however are very fractionated. (c and d) In the lower mantle, the trace elements considered have a better affinity for the bridgmanite. Liquids are progressively fractionated from the onset of solidification until the end of the crystallization.</figcaption></figure>

snapshot in time captures the end of mantle solidification (corresponding to
Figure 2 in the main text). We observed that the fluid dynamics simulations
(histograms) produce compositions consistent with the partitioning model
(solid and dash lines).

2 Scaling analysis and regime diagrams

2.1 Chemical differentiation near the top and bottom thermal boundary layers

Our fluid dynamics simulation revealed that chemical differentiation in magma oceans differs from what is usually thought. The canonical magma ocean model assumes that solidification and solid-liquid chemical fractionation occur at the depth where liquidus and adiabat intersect. This assumption is based on a 1D description of the temperature in a magma ocean. When considering a 2D temperature field, solids in fact precipitate at the surface of the planet in the top cold thermal boundary layer. This is because the surface temperature remains below the solidus during most of the duration of magma ocean crystallization. The crystal-rich layer is denser than the underlying mantle because both thermal and compositional density contrasts - we detailed below where
<figure><img src="image_46.png" /><figcaption>**Fig. S4** Benchmark of trace elements partitioning in the fluid dynamics simulations. Analytical predictions of both melt and solid composition (green and blue lines) estimated by a simple partitioning model (see Eq. (12) and Fig. S3). (Top) Partitioning coefficients correspond to the one of olivine. (bottom) Partitioning coefficients correspond to the one of bridgmanite. Assuming a melt fraction of 10 %, we can predict the composition of the melt and solid by determining the intersection of the dash (top) or solid (bottom) lines with a horizontal line drawn at y=0.1. This is in good agreement with the frequency distribution of the simulations. However, the fluid dynamics simulations provide a more complex picture of the compositional fields, revealing a distribution of composition rather than a unique composition.</figcaption></figure>

the origin of this compositional differences comes from - act in the same direction. These density contrasts drive the growth of Rayleigh-Taylor instabilities forming crystal-rich downwellings (Figure S7).

As shown in this study, crystals formed at the surface of planet play a major role in the differentiation of the mantle. Indeed, these solids experience shallow depth solid-liquid chemical fractionation governed by low-pressure mineral phases. These crystal then sink as downwellings in the deep mantle without remelting (see Section 2.3). This is the reason why the chemical evolution of
the deep mantle is expected to be significantly affected by low-pressure solid-liquid partitioning. However, to bring a fractionated composition into the deep mantle, solids and liquids must segregate from each other in the vicinity of the surface (cold) thermal boundary. Otherwise, it is the bulk composition of the liquid - that can also be fractionated at other locations - that is brought to the deep mantle.

Here, we perform a scaling analysis to compare the growth timescale of a Rayleigh-Taylor instability with that of solid-liquid segregation. First, we derive a dimensional equation to quickly estimate the timescale ratio using conservative values. Then, we reformulate it with dimensionless numbers to better relate our fluid dynamics simulation to realistic conditions.

The growth timescale for a Rayleigh-Taylor instability driven by a density anomaly of size $R_0$ in a domain of viscosity $\eta$ scales as:

$$ \tau_{RT} \sim \frac{\eta}{\Delta\rho_{RT}gR_0}. \qquad (13) $$

where $\eta$ is a function of viscosity, $\eta = \eta_s \exp(-30\phi)$ [9, 10]. Assuming that the compaction length is smaller than the thickness of the thermal boundary layer - we will come back to this assumption later - the timescale for solid-liquid segregation in this layer that we assume to be mostly solid, is controlled by the melt percolation velocity:

$$ v_{\text{perco}} = \frac{k(\phi) \Delta\rho_{sl}g}{\phi \eta_l}, \qquad (14) $$

where $\Delta\rho_{sl}$ is the density difference between solid and liquid, and the function $k(\phi)$ describes the matrix permeability as function of the melt fraction. We use:

$$ k(\phi) = \frac{1}{C_0} r_0^2 \phi^n, \qquad (15) $$

where $r_0$ is the crystal size, $C_0$ is the tortuosity of the melt network, and $n$ is a power exponent that captures the melt network connectivity. Using the two equations above, the timescale of solid-liquid segregation over a thickness $R_0$ scales as

$$ \tau_{\text{seg}} \sim \frac{C_0 R_0 \eta_l}{\phi^{n-1} \Delta\rho_{sl} g r_0^2}. \qquad (16) $$

Solid can segregate from the top thermal boundary layer if

$$ \frac{\tau_{RT}}{\tau_{\text{seg}}} > 1. \qquad (17) $$
Using Equations (13) and (16) and assuming that the density contrasts $\Delta\rho_{sl}$ and $\Delta\rho_{RT}$ are comparable yields:

$$ \frac{\tau_{RT}}{\tau_{\text{seg}}} \sim \frac{\eta}{\eta_l} \frac{\phi^{n-1}}{C_0} \left(\frac{r_0}{R_0}\right)^2. \qquad (18) $$

For a realistic magma ocean, we use the conservative value of $R_0 = 100$ km (today's Earth's oceanic lithosphere thickness), $r_0 = 1$ cm, $\eta_s = 10^{22}$ Pa s, $\eta_l = 10^2$ Pa s, $C_0 = 300$, $n = 3$ and $\phi = 0.15$. This yields:

$$ \left( \frac{\tau_{RT}}{\tau_{\text{seg}}} \right)_{\text{Earth}} \approx 0.83. \qquad (19) $$

While close to 1, this value remains a conservative estimate, because the thermal lithosphere was likely thinner during the Hadean due to significantly higher heat flux out of the Earth's interior. Additionally, Equation (18) shows that the ratio of these timescales depends on ratio of quantities that are not well-known, such as the grain size in a solidifying magma ocean, $r_0$, the viscosity of the melt $\eta_l$, the thickness of the upper thermal boundary layer, $R_0$, or the effective viscosity of the essentially solid upper thermal boundary, $\eta$.

To gain a more robust estimation of the ratio of the Rayleigh-Taylor and melt segregation timescales, we rely in the following on scaling laws to rewrite Equation (18) as a function of dimensionless numbers that describe our fluid dynamics simulations, e.g., $Ra$ and $\delta$. Finally, this approach allows one to relate the dynamic regime observed in our fluid dynamics simulations to the one of the real Earth.

The ratio of the crystal size to thickness of the thermal boundary layer in our simulations can be estimated using the definition of the melt mobility number (see Method),

$$ \left(\frac{r_0}{R_0}\right)^2 = C_0 \delta^2 \left(\frac{H}{R_0}\right)^2 \frac{\eta_l}{\eta_s}. \qquad (20) $$

Following [11] (their Equation (9.53)), the thickness of the uppermost thermal boundary layer can be related to the mantle thickness in the presence of large viscosity contrasts between the top cold layer and the hot interior:

$$ R_0 \sim H Ra^{-1/3} \left( \frac{\Delta T_m}{\Delta T_R} \right)^{4/3}, \qquad (21) $$

where $Ra$ correspond to the global Rayleigh number computed with the viscosity relevant for the interior. We make a Frank-Kamenetskii approximation [e.g., 12] where $\Delta T_R$ is a temperature scale that captures the variation of viscosity in the top cold boundary layer and in the hot interior (see e.g., Equation
10 Solidification of Earth's mantle led inevitably to a basal magma ocean

(9.36) in [11]),

$$ \frac{\Delta T_m}{\Delta T_R} = \ln \left( \frac{\eta}{\eta_{\text{int}}} \right). \qquad (22) $$

The ratio of the Rayleigh-Taylor instability growth to solid-liquid segregation timescales thus scales as:

$$ \left(\frac{\tau_{RT}}{\tau_{\text{seg}}}\right) \sim \frac{\eta}{\eta_s} \phi^{n-1} \delta^2 Ra^{2/3} \left[\ln\left(\frac{\eta}{\eta_{\text{int}}}\right)\right]^{-8/3}, \qquad (23) $$

where Ra corresponds to the Rayleigh number computed with the viscosity of the partially molten convective interior, $\eta_{\text{int}}$. We rewrite Ra as follows,

$$ Ra = Ra_{\text{ref}} \frac{\eta_s}{\eta_{\text{int}}}, \qquad (24) $$

where $Ra_{\text{ref}}$ is a thermal Rayleigh number computed with the viscosity of the solid, $\eta_s$. This yields:

$$ \left(\frac{\tau_{RT}}{\tau_{\text{seg}}}\right) \sim \frac{\eta}{\eta_s} \phi^{n-1} \delta^2 Ra_{\text{ref}}^{2/3} \left(\frac{\eta_s}{\eta_{\text{int}}}\right)^{2/3} \left[\ln\left(\frac{\eta}{\eta_{\text{int}}}\right)\right]^{-8/3}. \qquad (25) $$

In our fluid dynamic simulations, $\eta \approx \eta_s$ and $\eta_{\text{int}} \approx \eta_l$ when the melt fraction is far from the rheological critical melt fraction (see MT Equation (9)). We use $\delta = 1.12 \times 10^{-2}$ for our reference case (described in the main text), $Ra_{\text{ref}} = 10^6$, $\eta_s = 10^3$, $\phi = 0.15$ (melt fraction in the vicinity of the boundary layer), $n = 2$ (as in our fluid dynamic model, see [13]), and find for our reference case:

$$ \left(\frac{\tau_{RT}}{\tau_{\text{seg}}}\right)_{\text{simulations}} \approx 0.11, \qquad (26) $$

indicating that solid-liquid differentiation in the thermal boundary top thermal boundary layer in our reference simulation is expected to occur, but may not be significant. We performed other simulations with $\delta = 3 \times 10^{-2}$, $Ra = 10^8$, $\eta_s = 10^2$, yielding to,

$$ \left(\frac{\tau_{RT}}{\tau_{\text{seg}}}\right)_{\text{simulations}} \approx 0.88, \qquad (27) $$

which is closer to unity. In our simulations, the maximum viscosity ratio between solids and liquid is $10^4$. However, we achieve a realistic ratio for these two timescales (and thus similar dynamic behavior) by artificially increasing the size of the crystals (see Figure S5a). By doing so, we compensate for the smaller viscosity contrast, ensuring that the dimensionless melt mobility
number is closer to Earth's conditions. As a result, phase segregation velocities become comparable to convective velocities driven by Stokes flow (see Section 2.2 and [13]).

For Earth-like conditions, we have $\eta = \eta_s \exp(-30\phi)$ [9, 10], but the ratio of the solid viscosity to the one of the interior - that controls the thickness of the upper thermal boundary layer - as well as the melt mobility are subject to large uncertainties. To anticipate the expected range of uncertainties, we plot in Figure S5b, the ratio of the Rayleigh-Taylor timescale and solid-liquid segregation as a function of the melt mobility number, and the viscosity contrast between the cold the surface and the interior. For completeness, we assume that the mushy interior is $10^5$ to $10^{15}$ times less viscous than the essentially solid upper thermal boundary layer. Assuming crystal sizes, $r_0$ between 1 mm and 10 cm, and a solid / liquid viscosity contrast between $10^{14}$ and $10^{20}$, we expect the mobility number for a real Earth to range between $10^{-5}$ and $10^2$ (see Figure S5a). For these values, there is a large parameter space where the ratio of the Rayleigh-Taylor timescale is several orders of magnitude larger than the timescale for solid-liquid segregation. This is the reason why we anticipate solid-liquid segregation within the top thermal boundary layer in our simulations. This process is expected to be more significant for the Earth - for sufficiently large melt mobility numbers and viscosity contrasts between the upper thermal boundary layer and the mushy interior - as indicated by this scaling analysis.

Finally, we tested our working hypothesis, that consists in neglecting compaction in the thermal boundary layer, by comparing the thickness of the upper thermal boundary layer (Equation (21)) with the compaction length (Equation (3.3) [14])—which is not directly equivalent to our melt mobility number:

$$L_{\text{comp}} = r_0 \sqrt{\frac{\phi^n \eta_s}{C_0 \eta_l}}. \qquad (28)$$

Figure S6 shows $R_0/L_{\text{comp}}$ as a function of the viscosity ratios. Values larger than one indicate that the thermal boundary layer is thicker than the compaction length, implying that compaction can be neglected. For plausible solid/liquid viscosity ratios, this analysis shows that the viscosity contrast between the thermal boundary layer and the mushy interior cannot be too large (i.e., $< 10^{14}$) for compaction to be negligible in the upper thermal boundary layer.

As depicted in Figure S7, the velocity field for the liquid phase (shown by red arrows) and for the solid phase (shown by blue arrows) near the top boundary layer in our reference simulation ($\delta = 1.12 \times 10^{-2}$, $Ra = 10^9$, $\eta_s / \eta_l = 10^3$). The figure illustrates that the motion of the two phases is primarily driven by phase averaged thermo-chemical density contrast without phase separation. However, we observe that in some regions the liquid velocity is deflected upwards relative to the solid velocity, which is attributed to the upward flux of melt driven by the density contrast between the solid and liquid at shallow depths. This is consistent with our scaling analysis that indicates a ratio of
<figure><img src="image_47.png" /><figcaption>**Fig. S5** (a) Melt mobility number as a function of crystal size and viscosity ratio between the liquid and solid phase. (b) Ratio of the Rayleigh-Taylor to the solid-liquid segregation timescales as a function of melt mobility number and the viscosity ratio between the upper thermal boundary layer and the interior.</figcaption></figure>

the Rayleigh-Taylor growth timescale and solid-liquid segregation smaller but close to 1.

Importantly, the solid-liquid segregation mechanism near the top boundary layer leads to solids in downwelling currents that are enriched in FeO, as shown in Figure S7. This is at odds with what one would expect from the values of partitioning coefficients. These solids should be consistently depleted in FeO compared to liquids. However, any crystal forming near the cold thermal boundary tends to settle, propelling FeO-rich melt upwards by mass conservation. This FeO-rich melt approaches the cold thermal boundary, cools down by diffusion, and eventually solidifies. As a result, the cold downwellings are enriched in FeO, making them thermally and chemically negatively buoyant, as seen in Figure S9.

The symmetric version of the processes described above also occurs in the deep mantle during the formation of upwelling currents at the (hot) bottom thermal boundary (Figure S10). The main difference is that in the deep mantle, melts segregate downwards due to a density crossover between the melt and solid at approximately mid-mantle depth. Near and within the hot thermal boundary layer, solid materials heat up by diffusion and may remelt. Any remaining crystal tend to float, driving FeO-rich melt downwards. Consequently, this denser FeO-rich melt accumulates at the core-mantle boundary (CMB), forming the basal magma ocean (Figure S10).
<figure><img src="image_48.png" /><figcaption>**Fig. S6** Ratio of the upper thermal boundary layer (TBL) thickness and the compaction length. Ratios larger than 1 indicate that compaction can be neglected.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>**Fig. S7** Melt (red) and solid (blue) velocity fields in the vicinity of the top thermal boundary region. In some regions the liquid velocity is deflected upwards relative to the solid velocity, which is attributed to the upward flux of melt driven by the density contrast between the solid and liquid at shallow depths. The top shows closed up views of the fluid parcels highlighted in magenta in the bottom row.</figcaption></figure>
14

Solidification of Earth's mantle led inevitably to a basal magma ocean

<figure><img src="image_50.png" /><figcaption>**Fig. S8** Same as Figure (S7) but for the lower mantle in the vicinity of the basal magma ocean (BMO). In several regions the liquid velocity (red arrows) is deflected downwards relative to the solid velocity (blue arrows), which is attributed to the downwards flux of melt driven by the density contrast between the solid and liquid in the lower mantle. The top panels show closed up views of the fluid parcels highlighted in magenta in the bottom row.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>**Fig. S9** Schematic of fluid dynamics processes involved in the chemical differentiation near the top thermal boundary layer. This is the interpretation of the velocity field of the fluid dynamics simulations shown in Figure S7.</figcaption></figure>
<figure><img src="image_52.png" /><figcaption>**Fig. S10** Schematic of fluid dynamics processes involved in the chemical differentiation near the bottom thermal boundary layer. This is the interpretation of the velocity field of the fluid dynamics simulations shown in Figure S8.</figcaption></figure>

## 2.2 Solid-liquid phase separation vs. entrainment by thermal convection

In the literature, MO thermal evolution models, that track the evolution of the MO melt fraction with time, are generally decoupled from mechanical models of solid-liquid segregation, that track where crystals accumulate. Historically, this choice was motivated by the fact that the first MO thermal models were mostly designed to investigate primordial atmosphere degassing, not mantle differentiation. In these work, the actual structure of the planet interior, i.e., the spatial distribution of crystals, was of secondary importance [15–18]. In fact, if the adiabat is steeper that the liquidus, and crystals remain denser than the melt at all depth, we should always expect an upper magma ocean on the top of solid cumulates. In that case, the location of crystals, e.g., averaged crystal fraction in the magma ocean and the amount of residual liquid in the cumulates, play a relatively minor role on the MO thermal evolution. The basal magma ocean hypothesis [19] changed that paradigm as the presence of a solid shell at mid-mantle depth may thermally decouple the upper magma ocean from the lower magma ocean. In that case, the location where crystals accumulate becomes crucial for the thermal evolution of the magma ocean. [20] proposed the first 1D thermal evolution model capable of accounting for a solid layer at mid-mantle depth. As no crystal spatial distribution is prescribed a priori, substantial numerical challenges arises from the time-integration of the energy equation with a large dynamic range of energy fluxes (see also [21]). However, in this kind of modelling, the temperature (or entropy) is not advected independently from the crystal fraction (but it may be done in theory). Assuming thermodynamic equilibrium, this implies that the spatial
distribution of crystals mimics exactly the relation between temperature and melting curves. In this model, the formation of a basal magma ocean is thus only dictated by temperature, not where crystals go. Here, we aim to propose an alternative scenario where the formation of the basal magma ocean is not dictated by the temperature profile, but by where crystals go. To do that, we must advect the crystal fraction independently from the temperature.

Independently from magma ocean thermal evolution models, the physics solid-liquid phase separation in magma oceans has received considerable attention since the development of the lunar magma ocean hypothesis [22, e.g.,]. Indeed, the physics of solid-liquid phase separation is crucial for the issue of mantle differentiation. Because large melt fractions are required for solid feldspar to segregate from the residual magma to form the primary lunar crust [23], models of magma ocean differentiation mostly focused on the physics of solid-liquid segregation at high melt fraction, i.e., crystal settling and suspension [20, 23–32].

Here we focus on magma ocean differentiation beyond the rheological transition, i.e., once the mantle reaches a solid-state like dynamics (with crystal fractions of about 50 %). Major differences are expected between magma ocean differentiation at low crystal fraction and high crystal fraction. First, the effective thermal Rayleigh number in a mushy magma ocean ($< 10^{15}$) is expected to be substantially lower than that of a fully molten magma ocean ($> 10^{27}$), making the dynamics of mushy magma oceans more closely accessible with current computational capabilities. In this case, solid-liquid segregation is not governed by crystal settling and re-entrainment but by melt percolation, matrix compaction, and shear deformation. Second, during the mushy magma ocean stage, strong feedbacks between the mechanics of solid-liquid segregation and phase change may occur because phase change produces compositional differences that affect solid-liquid density contrast in turn affecting local composition and thermodynamic equilibrium.

To cross-check the relevance of our numerical model that accounts for phase change, chemical fractionation, and solid-liquid viscosity contrast, we analytically predicted the transition between fractional and batch solidification using a simple approach. To do so, we followed previous studies that estimated the transition between the two magma ocean stages regimes mentioned above by comparing the magnitude of velocities associated with thermal convection and phase separation [13, 33]. This simple approach provides an upper bound for the transition between batch and fractional crystallization (see below). Crystals settling at the bottom of the magma ocean (where the vertical component of the convective velocities decreases) can still occur even if crystal velocities are substantially smaller than the magnitude of convective velocities [24, 27, 34]. Therefore, despite its simplicity, comparing the average convective velocity and phase separation appears relevant for validating our numerical approach. We describe below how these two velocities can be extracted from our mathematical description of the magma ocean multi-phase physics.
Mechanical differentiation is governed by competition between solid-liquid phase separation and thermal convection as we consider that the timescale of solidification is not a limiting factor. Indeed, in a sluggish state, the cooling timescale of the magma ocean is larger than the timescale of melt migration. We thus consider two end-member cases:

* **Fractional solidification:** Solid-liquid phase separation generates a compositional layering much faster than thermal convection can stir this layering. This regime promotes fractional crystallization where melt of distinct composition is perfectly extracted from the residual matrix before any large-scale convective motion can occur.

* **Batch solidification:** Thermal convection mechanically stirs any new formed layering faster than phase separation generates this compositional layering. This regime promotes batch crystallization where the mantle is constantly homogenized during solidification.

To estimate the mean velocity of convective motions driven by thermal density contrast, $v_{\text{mean}}^{\text{conv}}$, we use a scaling law of the form [35].

$$v_{\text{mean}}^{\text{conv}} = 0.16Ra^{\frac{3}{2}}. \qquad (29)$$

To estimate the average phase separation velocity, $v_{\text{mean}}^{\text{sep}}$, we use MT Equation (10) assuming that melt extraction is solely driven by density contrasts between melt and solid, and we consider a melt fraction of 50%. We obtain:

$$v_{\text{mean}}^{\text{sep}} = \frac{1}{4}Ra\delta^2. \qquad (30)$$

Importantly, the phase separation velocity is a function of the melt mobility number that accounts for both the melt viscosity and grain size. Fig. S11 maps regions in the $(Ra, \delta)$ space, where batch crystallization or fractional crystallization are expected to occur. The transition between the two regimes corresponds to $v_{\text{mean}}^{\text{conv}} = v_{\text{mean}}^{\text{sep}}$. Despite its simplicity, this prediction is in good agreement with the results from our numerical simulations (see Fig. S11). In our study, we conducted simulations that are in the vicinity of the boundary between the two regimes (Fig. S11). For a maximum thermal Rayleigh number of $10^9$, we observe inefficient differentiation for $\delta = 5 \times 10^{-3}$ (Fig. S12) and substantial differentiation for $\delta = 20 \times 10^{-3}$ (Fig. S13). As expected, we do still observe differentiation when $v_{\text{mean}}^{\text{conv}} > v_{\text{mean}}^{\text{sep}}$, but we do not observe differentiation when $v_{\text{mean}}^{\text{conv}} < v_{\text{mean}}^{\text{sep}}$. This is consistent with our prediction that provides an upper bound for the transition between batch and fractional crystallization as discussed above. This diagram allows to explore magma ocean dynamics in regimes that are accessible with current numerical limitations, but where the critical issue of batch vs. fractional crystallization can still be investigated in a quantitative manner (Fig. S11 and S16).

The regime diagram can also be used to predict the crystal sizes in magma ocean (Fig. S16). If we first consider that melt viscosity is about $10^4$ Pa s, the expected maximum thermal Rayleigh number is $10^{24}$. Assuming that the
<figure><img src="image_53.png" /><figcaption>**Fig. S11** Regime diagram in a Ra (thermal Rayleigh number) - δ (melt mobility) dimensionless space. 22 fluid dynamics simulations (colored diamonds) are used to quantify the amount of chemical differentiation in magma ocean. The maximum concentration of FeO serves as a proxy to measure the degree of chemical differentiation. These numerical simulations that capture several aspects of the magma ocean solidification, are in good agreement with our analytical predictions based purely on mechanical consideration (see text).</figcaption></figure>

Earth's magma ocean solidifies in the intermediate regime, $\delta$ must be around $10^{-5}$, this sets a maximum crystal size to 1 µm. If we consider a much lower melt viscosity value of about 1 Pa s [36], this further decreases the maximum crystal size to 0.01 µm. These values are consistent with previous estimates [27, 28, 37, 38].

Finally, it is important to identify the driving forces for solid-liquid segregation in our simulations. Although body forces, specifically the density contrast between solid and liquid, are generally the main drivers of solid-liquid segregation in the mushy mantle, this is not universally the case across the convective domain (see Figure S14). We observe that deviatoric stresses (the last term on the right side of MT Equation (10)) drive strongly phase separation at the crystallization front, where shear deformation and a large viscosity contrast occur within relatively small regions. Additionally, we note that compaction plays a non-negligible role in regions where there the melt fraction gradient is high. However, this effect is confined to very small regions, consistent with the small values of compaction length used in our simulations.

## 2.3 Crystal remelting during their descent

Solid silicates that form at a planet's surface retain signatures of low-pressure mineral-melt partitioning. In this context, we estimate the extent of remelting expected in a crystal-rich parcel of fluid as it descends into the warmer mantle. These crystal-rich blobs originate as Rayleigh-Taylor instabilities within
<table><caption>**Table S1** Values of the physical parameters used in this study.</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Mean density</td><td>ρ<sub>0</sub></td><td>4500</td><td>kg m<sup>-3</sup></td></tr><tr><td>Thermal expansivity</td><td>α</td><td>10<sup>-5</sup></td><td>K<sup>-1</sup></td></tr><tr><td>Gravity acceleration</td><td>g</td><td>9.81</td><td>ms<sup>-2</sup></td></tr><tr><td>Magma ocean thickness</td><td>H</td><td>3000</td><td>km</td></tr><tr><td>Thermal diffusivity</td><td>κ</td><td>10<sup>-6</sup></td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Solid-liquid isochemical density contrast at the surface</td><td>-</td><td>1400</td><td>kg m<sup>-3</sup></td></tr><tr><td>Solid-liquid isochemical density contrast at the CMB</td><td>-</td><td>100</td><td>kg m<sup>-3</sup></td></tr><tr><td>Maximum chemical density contrast at the surface</td><td>-</td><td>500</td><td>kg m<sup>-3</sup></td></tr><tr><td>Maximum chemical density contrast at the CMB</td><td>-</td><td>1500</td><td>kg m<sup>-3</sup></td></tr><tr><td>Thermal capacity</td><td>C<sub>p</sub></td><td>10<sup>3</sup></td><td>J K<sup>-1</sup> kg<sup>-1</sup></td></tr><tr><td>Latent heat of solidification</td><td>L</td><td>1000</td><td>kJ kg<sup>-1</sup></td></tr><tr><td>Solid viscosity</td><td>η<sub>s</sub></td><td>10<sup>22</sup></td><td>Pa s</td></tr><tr><td>Liquid viscosity</td><td>η<sub>l</sub></td><td>10<sup>0</sup></td><td>Pa s</td></tr><tr><td>Compaction viscosity</td><td>η<sub>c</sub></td><td>10<sup>22</sup></td><td>Pa s</td></tr><tr><td>Crystal size</td><td>r<sub>0</sub></td><td>10<sup>-6</sup> - 10<sup>-2</sup></td><td>m</td></tr><tr><td>Surface temperature</td><td>-</td><td>500</td><td>K</td></tr><tr><td>Core temperature</td><td>-</td><td>5000</td><td>K</td></tr><tr><td>Adiabatic temperature increase over the MO thickness</td><td>-</td><td>1800</td><td>K</td></tr><tr><td>Super-adiabatic temperature jump over the MO thickness</td><td>ΔT<sub>m</sub></td><td>2700</td><td>K</td></tr><tr><td>Mg-rich end-member melting temperature at the CMB</td><td>T<sup>s</sup><sub>MgO</sub></td><td>4700</td><td>K</td></tr><tr><td>Fe-rich end-member melting temperature at the CMB</td><td>T<sup>c</sup><sub>FeO</sub></td><td>3500</td><td>K</td></tr><tr><td>Mg-rich end-member melting temperature at the surface</td><td>T<sup>s</sup><sub>MgO</sub></td><td>2054</td><td>K</td></tr><tr><td>Fe-rich end-member melting temperature at the surface</td><td>T<sup>c</sup><sub>FeO</sub></td><td>1359</td><td>K</td></tr></tbody></table>

<table><caption>**Table S2** Dimensionless quantities used in this study. <sup>1</sup> These values are computed using the parameters presented in Table 2.2. <sup>2</sup> Those values are constrained by current computational limitations. The appropriate solid-liquid viscosity contrast and the melt mobility number cannot be reproduced in our model. A reasonable solid-liquid viscosity contrast is currently out of reach in this kind of fluid dynamics simulations. Even though a realistic melt mobility number cannot be used in our simulations, it is re-scaled in relation to the thermal Rayleigh number so that the competition between entrainment by thermal convection and solid-liquid phase separation remains relevant in our simulations [13]. <sup>3</sup> The depth-dependent iso-chemical Rayleigh number is parameterized as</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Earth's Magma Ocean <sup>1</sup></th><th>This Study <sup>2</sup></th></tr></thead><tbody><tr><td>Solid-state thermal Rayleigh number</td><td>R<sub>a</sub></td><td>3.22 × 10<sup>6</sup></td><td>10<sup>6</sup></td></tr><tr><td>Solid-liquid viscosity contrast</td><td>-</td><td>10<sup>22</sup></td><td>10<sup>3</sup></td></tr><tr><td>Phase Rayleigh number at the surface</td><td>R<sub>p surf</sub></td><td>11.523</td><td>11.523</td></tr><tr><td>Phase Rayleigh number at the CMB</td><td>R<sub>p cmb</sub></td><td>0.823</td><td>0.823</td></tr><tr><td>Iso. Rayleigh exponent<sup>3</sup></td><td>γ</td><td>-</td><td>6</td></tr><tr><td>Compo. Rayleigh number at the surface</td><td>R<sub>c surf</sub></td><td>4.115</td><td>4.115</td></tr><tr><td>Compo. Rayleigh number at the CMB</td><td>R<sub>c cmb</sub></td><td>12.346</td><td>12.346</td></tr><tr><td>Stefan number</td><td>S<sub>t</sub></td><td>0.37</td><td>0.37</td></tr><tr><td>Melt mobility number<sup>2</sup></td><td>δ</td><td>10<sup>-3</sup>-10<sup>1</sup></td><td>5 × 10<sup>-3</sup> - 20 × 10<sup>-3</sup></td></tr><tr><td>Compaction viscosity</td><td>ζ</td><td>1</td><td>1</td></tr><tr><td>CMB temperature</td><td>-</td><td>1.0</td><td>1.0</td></tr><tr><td>Liquidus at the CMB</td><td>-</td><td>0.88</td><td>0.88</td></tr><tr><td>Solidus at the CMB</td><td>-</td><td>0.44</td><td>0.44</td></tr><tr><td>Liquidus at the surface</td><td>-</td><td>0.57</td><td>0.57</td></tr><tr><td>Solidus at the surface</td><td>-</td><td>0.31</td><td>0.31</td></tr><tr><td>Liquidus binary loop exponent</td><td>m</td><td>-</td><td>1.7</td></tr><tr><td>Solidus binary loop exponent</td><td>n</td><td>-</td><td>1.1</td></tr></tbody></table>

$$R_p(z) = \frac{R_p^2 - R_p^{\text{cmb}}}{1 - e^{-\gamma}} (e^{-\gamma z} - 1) + R_p^{\text{surf}}, \text{ where } z \text{ is the dimensionless depth.}$$
<figure><img src="image_54.png" /><figcaption>**Fig. S12** Compositional fields near the end of MO solidification in the case where phase separation is inefficient ($\delta = 5 \times 10^{-3}$) at $Ra_{max} = 10^9$), see Fig. S11</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>**Fig. S13** Compositional fields near the end of MO solidification in the case where phase separation is efficient ($\delta = 20 \times 10^{-3}$) at $Ra_{max} = 10^9$, see Fig. S11).</figcaption></figure>

<table><caption>**Table S3** Trace elements partitioning coefficients. <sup>a</sup> Lower mantle values [7]. <sup>b</sup> Upper mantle values [8]. <sup>c</sup> We use the model of [39]. <sup>d</sup> To our knowledge, solid-liquid partitioning of W in silicates has not been measured experimentally. Therefore, we use the partitioning coefficients of Th that has potentially a similar behavior [40].</caption><thead><tr><th>Element</th><th>D<sub>i</sub> - Lower mantle <sup>a</sup></th><th>D<sub>i</sub> - Upper mantle <sup>b</sup></th><th>Abundances in BSE (in ppb) <sup>c</sup></th></tr></thead><tbody><tr><td>Sm</td><td>0.05</td><td>1.1 &times; 10<sup>-2</sup></td><td>406</td></tr><tr><td>Nd</td><td>0.016</td><td>4.2 &times; 10<sup>-4</sup></td><td>1250</td></tr><tr><td>Lu</td><td>0.79</td><td>3.9 &times; 10<sup>-2</sup></td><td>67.5</td></tr><tr><td>Hf</td><td>1.6</td><td>1.1 &times; 10<sup>-3</sup></td><td>283</td></tr><tr><td>W<sup>d</sup></td><td>0.007</td><td>5.2 &times; 10<sup>-5</sup></td><td>29</td></tr></tbody></table>
<figure><img src="image_56.png" /><figcaption>**Fig. S14** Snapshot of driving forces for solid-liquid phase separation for the simulations shown in Figure S7 and S9.</figcaption></figure>

the top thermal boundary layer of the mantle. Following the mathematical model proposed by [41], we describe melt production in an descending spherical blob, assuming that the blob’s temperature remains between the solidus and liquidus. The conservation of energy within the blob is expressed as (see Equation (5) in [41]):

$$ \rho C_p \frac{\partial T}{\partial t} + \rho L \frac{\partial \phi}{\partial t} - \alpha T \frac{\partial p}{\partial t} = k \nabla^2 T, \quad (31) $$

Here, the first term on the left side accounts for the internal heat change in the blob, the second term represents the energy associated with phase change, and the third term reflects the energy variations due to adiabatic compression (p represents the lithostatic pressure experienced by the blob). These energy sources are balanced only by thermal diffusion on the right side of the equation, where $k = \kappa \rho C_p$ denotes the thermal conductivity.

With the temperature held between the liquidus and solidus, temperature and melt fraction are inherently linked by the phase diagram under the assumption of thermodynamic equilibrium:

$$ \phi = f(T) \text{ for } T_s < T < T_l. \quad (32) $$
<figure><img src="image_57.png" /><figcaption>**Fig. S15** Sm/Nd (a) and Lu/Hf (b) compositional fields near the end of MO solidification shown in Fig. 2 of the main text. We plot in (c) a correlation map of the Sm/Nd and Lu/Hf fields that allows to quickly visualize the origin (e.g., depth of formation) of these compositional anomalies. Correlation fractionnation of Sm/Nd and Lu/Hf ratios (red) is a signature of olivine fractionation, while anticorrelation (blue) is a signature of bridgmanite fractionation.</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>**Fig. S16** Regime diagram of magma ocean dynamics. Large viscosity contrasts between solid ($10^{22}$ Pa s) and liquid (1 Pa s) challenges our understanding of magma ocean evolution. During solidification magma ocean covers a broad range of fluid dynamic regime, ranging from fully molten turbulent convection to solid-state convection. Convective vigor with a thermal Rayleigh number larger than $10^{13}$ is currently not accessible via numerical simulations. However, the question whether magma ocean solidifies in a batch or fractional scenario can still be addressed in convective regimes accessible with current numerical resources by adjusting the melt mobility number.</figcaption></figure>
For simplicity, we assume a linear relationship between melt fraction and temperature (which is a plausible working hypothesis):

$$ \frac{\partial \phi}{\partial t} = \frac{\partial \phi}{\partial T} \frac{\partial T}{\partial t} = a \frac{\partial T}{\partial t}, \qquad (33) $$

where a is a constant. Assuming the spherical blob descends at a constant velocity W, the lithostatic pressure experienced by the blob during its descent, p(t), is given by:

$$ p = \rho g W t \qquad (34) $$

We can now rewrite Equation (31) as:

$$ \rho(Cp + L a) \frac{\partial T}{\partial t} = k \nabla^2 T + \alpha T \rho g W, \qquad (35) $$

We non-dimensionalize Equation (35) using the timescale H/W and R₀, the radius of the spherical blob, as the length scale, yielding:

$$ \frac{\partial \theta}{\partial \tau} = \frac{1}{Pe} \nabla^2 \theta + A, \qquad (36) $$

where θ is the dimensionless temperature, τ is the dimensionless time, and the Péclet number, Pe, is the ratio of diffusion to advection timescales:

$$ Pe = \frac{WR_0^2}{\kappa_e H}, \qquad (37) $$

where κ_e = k/[ρ(Cp + L a)] is an effective diffusivity that accounts for latent heat release. The dimensionless number A corresponds to a dimensionless adiabatic gradient:

$$ A = \frac{\alpha T_0 g W}{C_p + L a}. \qquad (38) $$

Equation (36) shows that temperature (and thus melt fraction) variation is predominantly controlled by the adiabatic gradient when the Péclet number is large (indicative of rapid descent without significant heat exchange with the surrounding mantle). This adiabatic compression is unfavourable to remelting as pressure increases (due to the slopes of liquidus and adiabat) and preserves solids formed at shallow depths in solid form. This process can be seen as the opposite version of decompression melting.

Next, we evaluate the Péclet number for a mushy mantle. We start by assuming that the blob descends at Stokes velocity, primarily driven by thermal density variations (note that in our simulations, compositional factors actually increase the blob's density, therefore this is a conservative estimate):

$$ W = \frac{\rho_0 \alpha \Delta T_m g R_0^2}{2\eta}, \qquad (39) $$
where $\eta$ corresponds to the average fluid viscosity within the mushy mantle, and the temperature difference between the upper thermal boundary layer and the interior is $\Delta T_m/2$. Upon substitution of Equation (39) into Equation (37), we obtain:

$$Pe = \frac{\rho_0 \alpha g \Delta T_m R_0^4}{2 \eta \kappa_e H} \qquad (40)$$

As seen in the equation above, estimating $Pe$ requires the knowledge of the blob size, $R_0$. Assuming the blob matches the size of the upper cold thermal boundary layer, we apply Rayleigh-Nusselt scaling laws (Equation (21)) to link $R_0$ to the mantle thickness, $H$. We assume that the melt fraction varies from 0 to 100 % over a temperature range of 300 K [1, 2], therefore we define $a = \frac{1}{300} K^{-1}$, implying that $\frac{\kappa}{\kappa_e} \sim 1$. Combining Equations (40), and (21) leads to:

$$Pe \approx \frac{1}{2} Ra^{-1/3} \left[ \ln \left( \frac{\eta_s}{\eta} \right) \right]^{16/3} \qquad (41)$$

where $Ra$ is the global thermal Rayleigh number of the partially molten interior. If this requirement is met (i.e, $Pe \gg 1$), solids forming in the top boundary layer cannot remelt via heat diffusion during their descent (though they may remelt at the core-mantle boundary when the vertical advective velocity tapers). These solids will carry low-pressure chemical fractionation signatures into the deep mantle.

Since the thermal Rayleigh number also depends on the effective viscosity of the partially molten interior, we plot in Figure S17 the Péclet number as a function of the viscosity contrast between the surface and the partially molten interior, assuming a reference thermal Rayleigh number of $Ra_{ref} = 10^6$ computed with the viscosity $\eta_s = 10^{22}$ Pa s,

$$Ra = Ra_{ref} \frac{\eta_s}{\eta} \qquad (42)$$

Figure S17 shows that crystal remelting is not expected in a mushy magma ocean when the viscosity contrast between the solid surface and the mush is lower than approximately $10^{21}$. This criterion is also satisfied in our fluid dynamics simulation, where the maximum viscosity contrast between solid and liquid is $10^4$.

It is important to note that this scaling analysis is primarily based on a scaling law originally developed for mantle convection with strong viscosity contrasts driven by thermal variations. In the future, it would be interesting to test the validity of this scaling law for multi-phase systems, which also exhibit a wide range of viscosities due to changes in melt fraction.
<figure><img src="image_59.png" /><figcaption>**Fig. S17** Regime diagram for crystal remelting during their descent based on available scaling laws. A Péclet number value larger than one implies no solids remelting during their descent.</figcaption></figure>

# 3 Additional simulations and videos

## 3.1 Solidification with a super-isentropic melting curve in the deep mantle

We conducted a simulation where the melting curves first intersect the adiabat at mid-mantle depth [42]. The dynamics is very similar to the case where the liquidus first intersects the adiabat at the CMB. It must be noted that crystal-rich downwellings that form at the surface do not accumulate at mid-mantle depth as they remain thermally negatively buoyant. They pursue their descent and remelt in the deep mantle. This further underlines that the crossing of the liquidus and adiabat does not play a major role in the style of solidification of the magma ocean (Fig. S18).

## 3.2 Preservation of a primordial undegassed reservoir

To quantify the amount of degassing in our simulations, we follow the approach used in [43] where we relied on Lagrangian particles to track the minimum depth reached by the magma ocean for the entire duration of our simulation and compare them with exsolution depths for relevant volatiles. Each tracer therefore represents a magma ocean parcel. If a given tracer does not reach the exsolution depth, it remains undegassed, i.e., it retains its primordial noble gas signature and volatile content. In contrast, as soon as a given tracer reaches the exsolution depth it will exsolve all its super saturated volatiles and loose them to the atmosphere. For simplicity, we assume that the magma ocean is undegassed at the beginning of our simulations. We use Ra=10⁶, a solid-liquid
<figure><img src="image_60.png" /><figcaption>**Fig. S18** Case where the melting curves first intersect the adiabat at mid-mantle depth [42].</figcaption></figure>

viscosity contrast of $10^3$ and $\delta = 15 \times 10^{-3}$. Over the duration of the numerical experiments, the tracers travel an average vertical distance of 15 times mantle thickness. The results are shown in Fig. S19 which displays the amount of volatiles that remains in the magma ocean at the end of our simulation for four different exsolution depths. Our results are consistent with those of [43] (performed in the case of a fully molten magma ocean) that show that the exsolution depth must be significantly large to degas the entire mantle on short time interval. For an exsolution depth of 3 km, which is already an extreme upper limit for Earth conditions [43], we observe only a negligible amount of degassing in our simulations (see Fig. S19). These results are expected to depend on the transit times performed by the magma ocean (i.e., the number of times a magma ocean "travels" from one horizontal boundary to the other) [43]. In the present simulations, the number of transit is 15.

Note that in the case considered here the coexistence of the liquid and solid phases along with the possibility of multiple melting and solidification processes limits the applicability of the tracer approach proposed in Samuel and Salvador (2023). For this reason, our results are only a simplification of the reality that reflects the weak ability of mantle convective dynamics to exsolve its volatiles species. A more thorough investigation is therefore required to more finely assess the outgassing capability of mushy magma oceans to exsolve their volatiles but this could require a complete re-implementation of the numerical model to properly track dissolved volatiles in space and time
<figure><img src="image_61.png" /><figcaption>**Fig. S19** Amount of material exsolved after 15 transit times for four different exsolution depths for the case displayed in Fig. 1. The exsolution depth must be significantly large (far beyond plausible values) to allow complete degassing of the Earth's mantle during magma ocean solidification. For an exsolution depth of 3 km, which is already an extreme upper limit [43], we observe a negligible amount of degassing in our simulations.</figcaption></figure>

and across melting and solidification processes, which goes beyond the main purpose of our current study.

Regardless, the preservation of a primordial and undegassed mantle reservoir has recently been questioned [44], which concluded that while the rare gazes isotopic signals collected at hotspot locations may originate from separate reservoirs, these sources may be associated with low concentrations in volatiles.

Overall, these considerations suggest that our results are not incompatible with constraints arising from geochemical analyses of lava samples.

## 3.3 Videos

The videos of numerical simulations described in the main manuscript can be downloaded by using the following links.
We provide three files that correspond to the cases displayed in the Figures 1 (fig1Video.mp4) and 2 (fig2Video.mp4) of the main text, and the additional case that explores the role of a super-isentropic melting curve in the deep mantle (Section S2.1 in the supplementary information, mid_mantleVideo.mp4).

* Figure 1: Temperature and FeO fields.

* Figure 2: Trace Elements field.

* Super isentropic melting (Section S2.1): Temperature and FeO fields.
References

[1] Nabiei, F. et al. Investigating magma ocean solidification on earth through laser-heated diamond anvil cell experiments. *Geophysical Research Letters* e2021GL092446 (2021).

[2] Boukaré, C.-E., Ricard, Y. & Fiquet, G. Thermodynamics of the MgO-FeO-SiO₂ system up to 140 GPa: Application to the crystallization of Earth's magma ocean. *J Geophys. Res.: Solid Earth* **120** (9), 6085-6101 (2015). URL http://dx.doi.org/10.1002/2015JB011929. https://doi.org/10.1002/2015JB011929, 2015JB011929.

[3] Wallner, H. & Schmeling, H. Numerical models of mantle lithosphere weakening, erosion and delamination induced by melt extraction and emplacement. *International Journal of Earth Sciences* **105** (6), 1741-1760 (2016).

[4] Ricard, Y., Bercovici, D. & Schubert, G. A two-phase model for compaction and damage: 2. applications to compaction, deformation, and the role of interfacial surface tension. *Journal of Geophysical Research: Solid Earth* **106** (B5), 8907-8924 (2001).

[5] Oliveira, B., Afonso, J. C., Zlotnik, S. & Diez, P. Numerical modelling of multiphase multicomponent reactive transport in the earth's interior. *Geophysical Journal International* **212** (1), 345-388 (2018).

[6] Keller, T. & Suckale, J. A continuum model of multi-phase reactive transport in igneous systems. *Geophysical Journal International* **219** (1), 185-222 (2019).

[7] Corgne, A., Liebske, C., Wood, B. J., Rubie, D. C. & Frost, D. J. Silicate perovskite-melt partitioning of trace elements and geochemical signature of a deep perovskitic reservoir. *Geochimica et Cosmochimica Acta* **69** (2), 485-496 (2005).

[8] Kennedy, A., Lofgren, G. & Wasserburg, G. An experimental study of trace element partitioning between olivine, orthopyroxene and melt in chondrules: equilibrium values and kinetic effects. *Earth and Planetary Science Letters* **115** (1-4), 177-195 (1993).

[9] Kohlstedt, D. L., Bai, Q., Wang, Z.-C. & Mei, S. Rheology of partially molten rocks. *Physics and chemistry of partially molten rocks* 3-28 (2000).

[10] Costa, A. Viscosity of high crystal content melts: dependence on solid fraction. *Geophys. Res. Lett.* **32** (22) (2005).
30 Solidification of Earth's mantle led inevitably to a basal magma ocean

[11] Jaupart, C. & Mareschal, J.-C. Heat generation and transport in the earth. *Cambridge University Press* **2**, 345–365 (2010).

[12] Solomatov, V. S. & Moresi, L.-N. Scaling of time-dependent stagnant lid convection: Application to small-scale convection on earth and other terrestrial planets. *Journal of Geophysical Research: Solid Earth* **105** (B9), 21795–21817 (2000).

[13] Boukaré, C.-E. & Ricard, Y. Modeling phase separation and phase change for magma ocean solidification dynamics. *Geochemistry, Geophysics, Geosystems* **18** (9), 3385–3404 (2017).

[14] Mckenzie, D. The generation and compaction of partially molten rock. *Journal of Petrology* **25** (3), 713–765 (1984).

[15] Abe, Y. & Matsui, T. Early evolution of the Earth: Accretion, atmosphere formation, and thermal history. *J. Geophys. Res.: Solid Earth* **91** (B13), E291–E302 (1986). URL http://dx.doi.org/10.1029/JB091iB13p0E291. https://doi.org/10.1029/JB091iB13p0E291.

[16] Zahnle, K. J., Kasting, J. F. & Pollack, J. B. Evolution of a steam atmosphere during earth's accretion. *Icarus* **74** (1), 62–97 (1988).

[17] Hamano, K., Abe, Y. & Genda, H. Emergence of two types of terrestrial planet on solidification of magma ocean. *Nature* **497** (7451), 607–610 (2013).

[18] Lebrun, T. et al. Thermal evolution of an early magma ocean in interaction with the atmosphere. *Journal of Geophysical Research: Planets* **118** (6), 1155–1176 (2013).

[19] Labrosse, S., Hernlund, J. & Coltice, N. A crystallizing dense magma ocean at the base of the Earth's mantle. *Nature* **450** (7171), 866–869 (2007). URL http://www.nature.com/nature/journal/v450/n7171/full/nature06355.html. https://doi.org/10.1038/nature06355.

[20] Bower, D. J., Sanan, P. & Wolf, A. S. Numerical solution of a non-linear conservation law applicable to the interior dynamics of partially molten planets. *Physics of the earth and planetary interiors* **274**, 49–62 (2018).

[21] Monteux, J., Andrault, D. & Samuel, H. On the cooling of a deep terrestrial magma ocean. *Earth Planet. Sci. Lett.* **448**, 140 – 149 (2016). URL http://www.sciencedirect.com/science/article/pii/S0012821X16302199. https://doi.org/http://dx.doi.org/10.1016/j.epsl.2016.05.010.
[22] Warren, P. H. The magma ocean concept and lunar evolution. *Annu. Rev. Earth Planet. Sci.* **13** (1), 201–240 (1985).

[23] Suckale, J., Elkins-Tanton, L. T. & Sethian, J. A. Crystals stirred up: 2. numerical insights into the formation of the earliest crust on the moon. *Journal of Geophysical Research: Planets* **117** (E8), n/a–n/a (2012). URL http://dx.doi.org/10.1029/2012JE004067. https://doi.org/10.1029/2012JE004067, e08005.

[24] Martin, D. & Nokes, R. Crystal settling in a vigorously converting magma chamber. *Nature* **332** (6164), 534–536 (1988).

[25] Martin, D. & Nokes, R. A fluid-dynamical study of crystal settling in convecting magmas. *Journal of Petrology* **30** (6), 1471–1500 (1989).

[26] Tonks, W. & Melosh, H. The physics of crystal settling and suspension in a turbulent magma ocean. *Origin of the Earth* **1**, 151–174 (1990).

[27] Solomatov, V. S. & Stevenson, D. J. Suspension in convective layers and style of differentiation of a terrestrial magma ocean. *J. Geophys. Res: Planets* **98** (E3), 5375–5390 (1993). URL http://dx.doi.org/10.1029/92JE02948. https://doi.org/10.1029/92JE02948.

[28] Solomatov, V. S. & Stevenson, D. J. Nonfractional crystallization of a terrestrial magma ocean. *J. Geophys. Res: Planets* **98** (E3), 5391–5406 (1993). URL http://dx.doi.org/10.1029/92JE02579. https://doi.org/10.1029/92JE02579.

[29] Verhoeven, J. & Schmalzl, J. A numerical method for investigating crystal settling in convecting magma chambers. *Geochemistry Geophysics Geosystems* **10** (2009). https://doi.org/{10.1029/2009GC002509}.

[30] Le Bars, M., Wieczorek, M. A., Karatekin, Ö., Cébron, D. & Laneuville, M. An impact-driven dynamo for the early moon. *Nature* **479** (7372), 215–218 (2011).

[31] Maas, C. & Hansen, U. Effects of earth's rotation on the early differentiation of a terrestrial magma ocean. *Journal of Geophysical Research: Solid Earth* **120** (11), 7508–7525 (2015). URL http://dx.doi.org/10.1002/2015JB012053. https://doi.org/10.1002/2015JB012053.

[32] Maas, C. & Hansen, U. Dynamics of a terrestrial magma ocean under planetary rotation: A study in spherical geometry. *Earth and Planetary Science Letters* **513**, 81–94 (2019).

[33] Huppert, H. E. & Sparks, R. S. J. The fluid dynamics of a basaltic magma chamber replenished by influx of hot, dense ultrabasic magma.
32 *Solidification of Earth's mantle led inevitably to a basal magma ocean*
*Contributions to Mineralogy and Petrology 75 (3), 279–289 (1981).*

[34] Lavorel, G. & Le Bars, M. Sedimentation of particles in a vigorously convecting fluid. *Physical Review E* **80** (4), 046324 (2009).

[35] Pandey, A., Verma, M. K. & Mishra, P. K. Scaling of heat flux and energy spectrum for very large prandtl number convection. *Physical Review E* **89** (2), 023006 (2014).

[36] Dygert, N., Lin, J.-F., Marshall, E. W., Kono, Y. & Gardner, J. E. A low viscosity lunar magma ocean forms a stratified anorthitic flotation crust with mafic poor and rich units. *Geophysical Research Letters* **44** (22), 11,282–11,291 (2017). URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017GL075703. https://doi.org/https://doi.org/10.1002/2017GL075703, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2017GL075703.

[37] Solomatov, V. S. & Stevenson, D. J. Kinetics of crystal growth in a terrestrial magma ocean. *J. Geophys. Res: Planets* **98** (E3), 5407–5418 (1993). URL http://dx.doi.org/10.1029/92JE02839. https://doi.org/10.1029/92JE02839.

[38] Elkins-Tanton, L. T. Magma oceans in the inner solar system. *Annu. Rev. Earth Planet. Sci.* **40**, 113–139 (2012). https://doi.org/{10.1146/annurev-earth-042711-105503}.

[39] Mcdonough, W. & Sun, S. The composition of the earth. *Chemical Geology* **120** (3-4), 223–253 (1995). https://doi.org/10.1016/0009-2541(94)00140-4.

[40] Bédard, J. Partitioning coefficients between olivine and silicate melts. *Lithos* **83** (3-4), 394–419 (2005).

[41] Katz, R. F. & Rudge, J. F. The energetics of melting fertile heterogeneities within the depleted mantle. *Geochemistry, Geophysics, Geosystems* **12** (10) (2011).

[42] Fiquet, G. et al. Melting of peridotite to 140 gigapascals. *Science* **329** (5998), 1516–1518 (2010). URL http://www.sciencemag.org/content/329/5998/1516full.pdf. https://doi.org/10.1126/science.1192448, http://www.sciencemag.org/content/329/5998/1516.pdf.

[43] Salvador, A. & Samuel, H. Convective outgassing efficiency in planetary magma oceans: insights from computational fluid dynamics. *Icarus* **390**, 115265 (2023).
[44] Parai, R. A dry ancient plume mantle from noble gas isotopes. *Proceedings of the National Academy of Sciences of the United States of America* **119** (29), 1–9 (2022). https://doi.org/10.1073/pnas.2201815119.