We follow the modeling strategy of Elkins-Tanton et al. (2003) and Charlier et al. (2018), who developed a numerical code to simulate the crystallization of the Martian and lunar magma oceans, respectively. Our code (Saracino et al., 2026) is rooted in experimentally determined phase equilibria and traces the solidification steps of the MMO through a stepwise approach in which we calculate the mineralogy, composition, and density of cumulates while tracking the compositional evolution of the residual silicate liquids. The initial compositions of the silicate liquid are the primitive compositions Mer8 and Mer15. Major elements and heat-producing elements (U, Th, K) were considered. FeO was excluded from the input compositions because MESSENGER data have demonstrated the paucity of FeO on the surface (Nittler et al., 2020). Fractional crystallization is described as (Charlier et al., 2018; Namur et al., 2012; Zhang et al., 2024):

c_{i,Liq}^{0} = (1 - z)c_{i,Liq}^{1} + z\left[\left(\sum_{j = 1 \rightarrow n}^{}{X^{j}c_{i}^{j}}\right)X_{Sol}^{Mush} + c_{i,Liq}^{0}X_{Liq}^{Mush}\right] \tag{1}

where z is the crystallization increment (fixed at 1 % of the residual liquid; i.e. z becomes increasingly small as crystallization proceeds), c_{i,Liq}^{1}is the concentration of element i, bonded to oxygen, in the liquid at each step of fractionation, c_{i,Liq}^{0} is that during the previous step of fractionation, c_{i}^{j} is the concentration of element i in solid phase j, X^{j} is the proportion of phase j in the cumulate assemblage, and X_{Sol}^{Mush} and X_{Liq}^{Mush} are the bulk proportions of the solid phases and trapped liquid in the solidifying crystal mush.

To determine the density of the silicate liquid produced during MMO crystallization, we first calculated the density at zero pressure as:

\rho_{Liq} = \sum_{}^{}\frac{Χ_{i}M_{i}}{Χ_{i}V_{i}} \tag{2}

where Χ_{i} is the molar fraction of component i in the melt, M_{i} is its molar mass, and V_{i} is its partial molar volume. Density is then adjusted to account for the effect of pressure and temperature with a third order Birch-Murnaghan equation of state:

P = \frac{3}{2}k_{T}\left[ \left( \frac{V_{0}}{V_{f}} \right)^{\frac{7}{3}} - \left( \frac{V_{0}}{V_{f}} \right)^{\frac{5}{3}} \right] \left\{ 1 - \left( \frac{3}{4} \right) \left( 4 - k_{T}' \right)\left[ \left( \frac{V_{0}}{V_{f}} \right)^{\frac{2}{3}} - 1 \right] \right\} \tag{3}

where P is pressure, k_{T} is Young’s modulus, V_{0} is the initial molar volume, V_{f} is the molar volume at pressure P, and k_{T}' is the pressure derivative of k_{T}. k_{T} and k_{T}' were determined via multiple linear regression of literature data as a function of major element melt compositions (Guillot and Sator, 2007; Ohtani and Maeda, 2001). The temperature dependence of the liquid density is calculated using the compositional expression of thermal expansion from Lange and Carmichael (1987). The parameters used to calculate the melt density are reported in the Supplementary Material (Section S.1). Other equations of state used for the calculation of the liquid density are shown for comparison in Figure S1 (in the Supplementary Material). The effect of S on the density of the silicate liquid is poorly constrained in the literature, and a specific discussion of this aspect is presented in Section 6.4 MMO density and potential mineral flotation.

Magma oceans are considered to crystallize from the bottom up and to strongly convect, such that no chemical gradient is assumed (Brown and Elkins-Tanton, 2009; Walker et al., 1975). Starting from the density of the cumulus phases, the thickness of the different layers can be calculated as crystallization proceeds at each increment z.

The initial depth of the MMO was calculated using gravity models for Mercury. A smaller polar moment of inertia (MOI; 0.333 ± 0.005; Genova et al., 2019) yields a deeper core-mantle boundary (CMB), and therefore a deeper interface between the magma ocean and the metallic core (\sim 485 ± 20 km depth in Goossens et al., 2022; \sim 463 km in Steinbrügge et al., 2021). Conversely, models with a higher MOI (0.346 ± 0.014; Margot et al., 2012) return a shallower CMB interface at \sim 430 km depth (Bertone et al., 2021). In this work, we simulate the crystallization of the MMO from an initial depth of 480 km, corresponding to the maximum possible CMB depth. Accordingly, we account for a higher volume of the silicate magma ocean, which is considered to have been \sim 17 % greater than the present-day solid silicate shell.

Literature data on near-liquidus experiments representative of the BSMe provide precious insights into the phase equilibria of early forming minerals (Anzures et al., 2020; Saracino et al., 2025; Xu et al., 2024). Xu et al. (2024) and Saracino et al. (2025) showed that enstatite is the liquidus phase in sulfide-saturated Mer8 (Mg/Si = 0.88) whereas forsterite is the liquidus phase in sulfide-saturated Mer15 (Mg/Si = 1.16) between 1 and 7 GPa. Interestingly, Saracino et al. (2025) showed that the residual melt compositions of sulfide-saturated Mer8 and Mer15 converge, eventually reaching a single cotectic curve after sufficient differentiation. The initial fractionation of the MMO can therefore be simulated by fractionating only enstatite for Mer8 and only forsterite for Mer15 until the residual MMO liquids reach the forsterite-enstatite cotectic (Fig. 2). After reaching this cotectic curve, all Mercurian melts, including low- and high-Mg/Si melts, evolve towards the same eutectic where crystallization ends (Namur and Charlier, 2017). Accordingly, during the first stage of crystallization, the Mg/Si ratios of the Mer8 and Mer15 residual liquids converge until they cross at \sim 2.6 GPa (Fig. 3). To further evolve the MMO, we consider an average of the two residual liquid compositions (Fig. 3).

The second stage of the forward model is the co-crystallization of forsterite and enstatite. The instantaneous relative proportions of forsterite and enstatite in the cumulate assemblage can be obtained from the tangent to the cotectic curve. We thus plotted experimental melts saturated with both forsterite and enstatite relevant to Mercury’s mantle compositions (Saracino et al., 2025; Xu et al., 2024) in the diopside-SiO₂-forsterite ternary system (Fig. S2 in Section S.2). The plotted melts were produced from the same starting materials (Mer8 and Mer15) over a wide pressure range (1.5 GPa–7.0 GPa). Calculated mineral proportions from these experimental melts returned a range of relative proportions, from around Fo₄₂En₅₈ to Fo₁₀En₉₀, with higher pressures producing a higher Fo/En ratio due to the expansion of the enstatite stability field relative to that of forsterite with increasing pressure (Bowen, 1914; Chen and Presnall, 1975; Falloon and Green, 1988; Kushiro, 1969; Weng and Presnall, 2001: S2). As an approximation, we therefore considered a rough average of the relative proportions: Fo₂₅En₇₅.

Our experimental study aims to simulate the evolution of the silicate liquid in the MMO during cooling and with decreasing pressure as the cumulate pile thickens while the melt fraction (F) decreases. We used a stepwise approach to simulate the fractional crystallization of the BSMe. Starting composition Mer_1 corresponds to the cotectic composition at which forsterite and enstatite coexist, following the crystallization of forsterite from Mer15 and enstatite from Mer8 using the forward crystallization model described above (Fig. 3). Experiments were performed with the Mer_1 composition until clinopyroxene joined the mineral assemblage with forsterite and enstatite. We then used the new residual melt composition in equilibrium with clinopyroxene, Mer_2, in further experiments until the appearance of quartz, which also corresponded to the destabilization of forsterite due to the high activity of SiO₂ (a_SiO2). Lastly, a final composition produced from the Mer_2 residual melts, Mer_3, was used to reach the most evolved stages of magma ocean solidification. These three starting compositions are reported in Table 2. Each new series of experiments pertaining to new starting materials were performed at progressively lower pressure-temperature conditions to simulate the incipient cooling and thinning of the residual MMO. The P-T-F profile and the complete crystallization sequence were obtained by combining the experimental results and the forward crystallization model. In total, our experiments span temperatures of 1525 ℃–1125 ℃ and pressures of 1.5 GPa–0.5 GPa (Fig. 4).

The starting materials (Table 2) were produced from high-purity oxide powders: SiO₂, TiO₂, Al₂O₃, Cr₂O₃, MnO, MgO, CaSiO₃, Na₂SiO₃, K₂Si₄O₉, and AlPO₄. CaSiO₃, Na₂SiO₃, and K₂Si₄O₉ were produced by decarbonating mixtures of carbonates and SiO₂. Sulfur was added as FeS, FeS + S, or only as elemental S powder. Powders were mixed in an agate mortar with methanol and then stored in an oven at 120 ℃. Different oxygen fugacity (fO₂) conditions were achieved by using different metallic Si/SiO₂ molar ratios in the starting materials: 0, 0.10, 0.20 (referred to as Mer(0), Mer(10), Mer(20)). Although fO₂ may not affect the stability of silicate phases in the Fe-free system, it affects S solubility in the melt and therefore the timing and nature of sulfide saturation (Namur et al., 2016b). With this approach, we were able to track the range of sulfides that can be produced using different Fe/S ratios. Neither S nor Si were added to Mer_3(0) because the solubility of S was observed to be low (< 0.5 wt % S) in Mer_2 at temperatures below 1250 ℃, in accordance with the decreasing solubility of S with decreasing T and increasing SiO₂ content (Namur et al., 2016b).

Experiments were conducted with a Voggenreiter Mavo LPC 250-300/50 end-loaded piston-cylinder apparatus at the University of Liège (Belgium). Half inch (1/2") assemblies were used in the pressure range 1.0 GPa–1.5 GPa, and 3/4" assemblies for runs at 0.5 GPa. For the 1/2" assembly, a graphite capsule with a MgO spacer was placed in a graphite furnace, and a BaCO₃ cell was used as the pressure medium. In the 3/4" assembly, a talc + pyrex cylinder was used as the pressure medium. Details on the cell assemblies employed in this study are reported in the Supplementary Material (Fig. S3). Temperature was monitored with a W₇₅R₂₅/W₉₇R₃ D-type thermocouple for both assemblies. A correction was applied for the temperature gradient between the tip of the thermocouple and the center of the capsule, which we estimate to be \sim 25 ℃. The assemblies were first pressurized at room temperature to 0.7 GPa in the 1/2" setup and 0.4 GPa in the 3/4" setup. Then, temperature was increased to 865 ℃ at 100 ℃/min while keeping the pressure constant. Temperature was then held for 6 min to pressurize to the target pressure. Finally, temperature was increased again at 50 ℃/min until the target temperature was attained. Friction correction coefficients of 9.3 % for the 1/2" assembly and 20 % for the 3/4" assembly were applied, respectively (Condamine et al., 2022). Experiments were quenched by switching off the power. Experimental samples were cut in half using a diamond wafer saw, mounted in epoxy, and ultimately polished with an alcohol-based polycrystalline 1 µm diamond suspension preserving all sulfides.

We performed time-series experiments to test the attainment of equilibrium and the potential loss of volatiles through the graphite capsule (Na, S). Four experiments at 1300 ℃ and 1 GPa were run at different durations (30 min, 1 h, 3 h, 6 h) using Mer_2(20) as the starting material. In all of these samples, enstatite, clinopyroxene, and quartz coexisted with the silicate melt, and these three crystalline phases produced small euhedral crystals regardless of the duration investigated. As commonly observed in experiments with Si/SiO₂ = 0.20, a slight excess of Si metal (\sim 1 wt %) was observed in all four experiments. The compositions of the silicate glasses (Fig. S4 in the Supplementary Material; Saracino et al., 2026) show that sulfur was readily incorporated into the melt after 30 min and slightly decreased after one hour (from 3.4 to 2.7 wt %; Fig. S4d). No substantial loss of Na was observed at durations up to 6 h (Fig. S4c). We also observed that the concentration of Si in the silicate glass slightly increased at durations longer than 1 h (Fig. S4a), explained by the likely oxidation of residual Si metal in longer duration runs. The progressive oxidation of the system is corroborated by the slight depletion of Si in the metallic phase (Fig. S5) over durations longer than 1 h and the lower solubility of S in the silicate melt in longer experiments. Based on these results, we conclude that we approached thermodynamic equilibrium within \sim 1 h.

Imaging and phase identification were performed with the TESCAN MIRA 4^th Generation Scanning Electron Microscope (SEM) at the Department of Geology, KU Leuven (Belgium). Quantitative analyses of our experimental products were conducted with the JEOL JXA-8530F electron probe micro-analyser (EPMA) at the Department of Materials Engineering, KU Leuven. Silicate crystals and sulfide and metallic phases were measured with a 15 kV accelerating voltage and 10 nA beam current. Silicate melts (glass) were measured with a 15 kV accelerating voltage and a 15 nA beam current. Crystals were analyzed with a focused beam and glasses with a beam defocused to 5 µm–20 µm depending on the size of the glass pool. Sulfides were measured with a focused beam and metals with either a focused or a defocused beam (5 µm–10 µm) depending on the size of the metallic phase and the presence of heterogeneities and quench textures. Peak counting times were 10 s–30 s and background counting times on each side of the peak were 5 s–15 s for each element. We employed both natural and synthetic standards. For silicate crystals we used albite for Na, orthoclase for K and Al, diopside for Ca, olivine for Si and Mg, rutile for Ti, MnO for Mn, Cr₂O₃ for Cr, and fayalite for Fe. For silicate glass we used albite for Na, obsidian for K and Al, olivine for Mg, diopside for Ca, rutile for Ti, rhodonite for Mn, Cr₂O₃ for Cr, hematite for Fe, and BaSO₄ for S. For sulfide and metallic phases we used albite for Na, Si metal for Si, diopside for Ca, olivine for Mg, obsidian for Al, rutile for Ti, Cr metal for Cr, Fe metal for Fe, Mn metal for Mn, and pyrite for S. Replicate measurements of international standards indicate analytical errors of less than 5 % for major elements and 10 %–15 % for minor elements. For sulfur measurements in the silicate glass, reliability was attained through accurate calibration and reproducibility on BaSO₄, which suggest errors below 10 %–15 %. Raw data were corrected for matrix effects with the CATZAF software.

Under reduced conditions, the oxygen fugacity of the experimental charges can be calculated considering the equilibrium between Fe metal and FeO-bearing silicate melts as (Corgne et al., 2008):

{Fe}_{Metal} + \ \frac{1}{2}O_{2}\ = \ {FeO}_{Melt} \tag{4}

The activity of Fe in the metal a_{Fe} is given by:

a_{Fe} = Χ_{Fe}\gamma_{Fe} \tag{5}

where Χ_{Fe} is the molar content of Fe in the Fe metal and \gamma_{Fe} is the activity coefficient of Fe. Analogously, the activity of FeO in the silicate melt a_{FeO} is given by:

a_{FeO} = Χ_{FeO}\gamma_{FeO} \tag{6}

where Χ_{FeO} is the molar content of FeO in the silicate melt and \gamma_{FeO} is the activity coefficient of FeO. However, typical Mercurian melt compositions have extremely low FeO contents (< 1 wt %) and are close to or below the detection limit of analytical instruments. Measurements uncertainties for low FeO concentrations may therefore limit the accuracy of fO₂ estimates (Anzures et al., 2020; Cartier et al., 2014; Namur et al., 2016a; Pirotte et al., 2023; Steenstra et al., 2020a). Nonetheless, under FeO-poor conditions, fO₂ can be determined using the equilibrium between Si-rich metal and SiO₂ in the silicate melt as (Cartier et al., 2014):

{SiO}_{2,\ Melt} = {Si}_{Metal} + \ O_{2} \tag{7}

where {SiO}_{2,\ Melt} is the SiO₂ concentration in the melt and {Si}_{Metal} is the Si concentration in the metal. The activity of Si in the metal is given by:

a_{Si} = Χ_{Si}\gamma_{Si} \tag{8}

with Χ_{Si} being the molar Si content in the metal and \gamma_{Si} the activity coefficient of Si. Here we use the Si activity coefficients calculated following Ma (2001). Similarly, the activity of SiO₂ in the silicate melt is calculated as:

a_{{SiO}_{2}} = Χ_{{SiO}_{2}}\gamma_{{SiO}_{2}} \tag{9}

where Χ_{{SiO}_{2}} is the molar SiO₂ content in the melt and \gamma_{{SiO}_{2}} is the activity coefficient of SiO₂. For experiments that do not saturate silica phases (A399, A397, A396, A404), a_{{SiO}_{2}} was calculated based on the reaction Mg₂SiO₄ + SiO₂ = Mg₂Si₂O₆ (Saracino et al., 2025). In experiments that saturate silica phases, a_{{SiO}_{2}} can be approximated as \sim 1 (Cartier et al., 2014). Only for runs in which no Si-bearing Fe metal was available (A461, A462, A469) did we calculate oxygen fugacity following the method of Namur et al. (2016b), where the oxygen fugacity can be calculated based on the sulfur content at sulfide saturation (SCSS) as:

\ln\lbrack S\rbrack_{SCSS} = a + \frac{b}{T} + \frac{cP}{T} + d \log(fO_{2}) + \sum e_{i} \frac{X_{i}}{X_{SiO_{2}}} \tag{10}

where \left[S\right]_{SCSS} is the SCSS (in wt%), T is temperature (in K), P is pressure (in bar), fO₂ is the oxygen fugacity, X_i are the molar fractions of oxides in the silicate melt (renormalized to 100 % on a S-free basis), and a, b, c, d, eTiO₂, eMgO, and eNa₂O are derived coefficients.

The fO₂ in experiments containing Si-bearing Fe metal spans the range from IW-5.8 to IW-8.4. The coexistence of excess Si metal and quartz suggests that some experiments are buffered at Si–SiO₂ equilibrium. For experiments that did not contain Si-bearing metal, we calculated the oxygen fugacity to be between IW-3.7 and IW-3.9.

Representative backscattered electron images of the experimental products are shown in Figure 5. All experimental runs, P-T-fO₂ conditions, and phase assemblages are reported in Table 3. Mer_1 contained enstatite as the liquidus mineral at T = 1500 ℃, so slightly offsets from the predicted enstatite-forsterite cotectic. Forsterite and clinopyroxene then appeared at T = 1450 ℃. Enstatite and forsterite occurred as relatively large (5 µm–50 µm) euhedral crystals (Fig. 5a), and clinopyroxene as small (< 10 µm), bright crystals often associated with enstatite. For Mer_1, we also ran experiments using a mix of different sulfur sources. Mer_1 + FeS,S still contained enstatite as the liquidus mineral at T = 1500 ℃, and the melt later saturated with enstatite + forsterite + clinopyroxene at T = 1450 ℃. Lastly, anhedral quartz occurred at T = 1400 ℃ (Fig. 5b). In contrast, Mer_1 + S contained silicate melt saturated with enstatite + forsterite + clinopyroxene + quartz at T = 1500 ℃. In Mer_2(20), enstatite + forsterite + clinopyroxene + quartz co-crystallized at T = 1350 ℃. Forsterite later disappeared at T = 1325 ℃. Forsterite was small, rare, and usually surrounded by small (2 µm–20 µm), euhedral crystals of enstatite (Fig. 5d). Clinopyroxene was small and often displayed elongated, needle–like textures (Fig. 5e). Quartz appeared as relatively large crystals (5 µm–50 µm), always with rounded edges. In Mer_2(10), enstatite + forsterite + clinopyroxene + quartz were the cotectic assemblage at T = 1300 ℃. Forsterite was destabilized at T = 1250 ℃. Lastly, Mer_3 runs contained an assemblage of quartz, plagioclase, and clinopyroxene starting at T = 1250 ℃ (Table 3). Quartz crystals were large (5 µm–70 µm), whereas clinopyroxene crystals were small (< 10 µm; Fig. 5f). Plagioclase crystals were small (< 10 µm), tabular, and always found as clusters (Fig. 5f). All SiO₂ phases featured in our experiments fall in the quartz stability field (Presnall, 1995). Our experimental samples also contained metal and sulfide phases. Mer_1 contained Si-bearing Fe metal (Fig. 5a) surrounded first by a thick liquid FeS layer, then by a small (Mg,Fe,Ca)S veneer. Among experiments with different sulfur sources, Mer_1 + FeS,S mostly contained sulfide globules of unmixed FeS and (Mg,Fe,Ca)S. At 1400 ℃, we also observed pure Si metal. Conversely, Mer_1 + S only contained excess Si metal globules. We also identified some small euhedral Cr sulfide crystals at 1450 ℃ (Fig. 5c). We are not aware of any other experimental work reporting this phase, although Cr-bearing sulfides have been reported in chondritic meteorites (Ma et al., 2011). Mer_2(20) mostly exhibited liquid FeS and (Mg,Fe,Ca)S droplets. Si-bearing Fe metal (Fig. 5d) and Si metal globules appeared at T = 1400 ℃. Mer_2(10) contained mostly Si metal and liquid FeS. Fe metal was very rare in our experiments, with only one small Si-free globule observed in only one sample (A462).

Quantitative analyses performed on the silicate glasses are shown in Table 4. Experimental liquid lines of descent (LLDs) are shown in Figures 6 and S6. The Si content constantly increased (26 wt %–37 wt %, Fig. 6a) from 1550 to 1100 ℃, whereas the Mg content steadily dropped (15 wt %–0 wt %, Fig. S6a). Na was also progressively enriched in the glass (1 wt %–6 wt %, Fig. 6d) with decreasing temperature. S content decreased from 1550 to 1250 ℃ (no S was added to experiments run at < 1250 ℃, Fig. S6f). The difference in temperature between the appearance of quartz when different S sources are added in Mer_1 (Table 3) might be a cause of the differing Si content in the melt (Table 4), which is affected by sulfur speciation. The glass in the more reduced experiment Mer_2(20) contained 2.0 wt %–4.5 wt % S, whereas that in the more oxidized experiment Mer_2(10) contained 0 wt %–1 wt % S. Both Al and Ca were enriched in the glass from 1550 to 1250 ℃, but were then depleted at lower temperatures (Fig. 6b–c) due to the saturation of plagioclase and clinopyroxene. More Ca was present in Mer_2 when 10 wt % Si was added than when 20 wt % Si added (Fig. 6c). This was likely the result of more Ca partitioning into the sulfide phases at the more reduced conditions in Mer_2(20). In all experiments, olivine and orthopyroxene matched their Mg endmember compositions (forsterite and enstatite, respectively). We only observed variations in CaO and Al₂O₃ contents in the enstatite composition over the range of starting materials investigated: both CaO (0.2 wt %–2.7 wt %) and Al₂O₃ contents (1.2 wt %–6.6 wt %) increased in enstatite with decreasing temperature (Saracino et al., 2026). In Mer_1 and Mer_2, most clinopyroxene crystals were too small for compositional analysis, with the exception of a few crystals at 1250 ℃ (A469). In Mer_3, clinopyroxene crystals were larger and contained 3 wt %–10 wt % Al₂O₃, 11 wt %–19 wt % MgO, 20 wt %–22 wt % CaO, and 0.7 wt %–1.8 wt % Na₂O. The average clinopyroxene composition was Wo₄₈En₅₂ (Fig. S7). Plagioclase was only observed in Mer_3(0) and became more sodic with decreasing temperature (An₈₄ to An₅₆, where An = molar Ca/[Ca + Na]). The compositions of metals and sulfides are reported in Saracino et al. (2026). Fe metallic alloys in Mer_1 + FeS contained 83 wt %–87 wt % Fe, 6 wt %–10 wt % Si, and 0.2 wt %–0.8 wt % S. Metal alloy in Mer_2(20) + FeS contained 82 wt %–87 wt % Fe, 13 wt %–20 wt % Si, and 0.1 wt %–0.3 wt % S. In Mer_2(10) + FeS we observed a single small Fe metal globule containing \sim 93 wt % Fe, 0 wt % Si, and \sim 1.7 wt % S in a single run at 1300 ℃. FeS in Mer_1 contained 51 wt %–57 wt % Fe, 36 wt %–42 wt % S, 2 wt %–6 wt % Cr, and 0.7 wt %–4.7 wt % Ti. Troilite in Mer_2 had a higher Fe content (60 wt %–64 wt %) and lower S (34 wt %–38 wt %), Cr (0.2 wt %–0.9 wt %), and Ti contents (< 1.7 wt %).

The experimentally determined phase relations (Table 3) were used to mimic the fractional crystallization of the MMO. We calculated the cotectic proportions of mineral phases in order to model melt differentiation consistent with the experimental LLDs (Table 4). The cotectic proportions of phases used in our modeling are listed in Section S.8 (Table S2 of the Supplementary Material). We used temperature as the main criterion to define the saturation or the destabilization of solid phases in the growing cumulate pile. We developed a S-free liquid thermometer based on experimental studies of Mercury-like mantle compositions (Charlier et al., 2013; Namur et al., 2016a; Namur and Charlier, 2017; Saracino et al., 2025, Figs S8, S9 in Section S.7) and applied a correction for the depression of the liquidus temperature in the presence of sulfur (Saracino et al., 2025, Section S.7). The phase assemblages considered in our magma ocean modeling as a function of temperature are reported in Section S.9 (Fig. S10, Table S3). The densities of crystals in the cumulate pile were determined at the solidus temperature derived by fitting the solidi of the evolving MMO liquid as determined using MAGEMin (Riel et al., 2022). The effects of temperature and pressure on mineral densities were calculated with a third-order Birch-Murnaghan equation of state. The parameters used in the calculation of mineral densities are reported in the Supplementary Material (Table S1 in Section S.1).

The evolution of the modeled MMO liquids for both Mer8 and Mer15 and a comparison with the experimental melts are reported in Figure 7. As a first approximation, we assumed that no trapped liquid is present in the growing cumulate assemblage so that X_{Liq}^{Mush} = 0 in Equation 1. The influence of Mg- and Ca-rich sulfides on the evolution of the MMO liquid is also shown (discussed further in Section 6.3 Sulfides in the mantle of Mercury). The crystallization sequence of Mer8 and Mer15 as MMO bulk compositions is shown in Figure 8 as a function of depth. In Mer8, the MMO solidification starts with the formation of a \sim 90 km thick basal orthopyroxenitic layer (Fig. 8a), whereas Mer15 forms a \sim 300 km thick dunitic basal layer (Fig. 8b). The residual melt of Mer8 becomes co-saturated in forsterite and enstatite at F = 0.81, whereas Mer15 reaches the forsterite-enstatite cotectic at F = 0.38. As also shown in Figure 2 for Mer15, a thick dunite layer is produced before the MMO residual liquid saturates enstatite too. Here, the compositions of the residual melts from sulfide-saturated Mer8 and Mer15 become similar (Fig. 7). Over this range of F (1.00–0.81 for Mer8 and 1.00–0.38 for Mer15), Si, Al, and Ca contents increase in the melt and the Mg content decreases due to the formation of enstatite and/or forsterite (Fig. 7b). Clinopyroxene is the next phase to appear at F = 0.40 in Mer8 and F = 0.35 in Mer15. Forsterite-enstatite-clinopyroxene assemblages are then followed by the appearance of quartz at F = 0.24 in Mer8 and at F = 0.28 in Mer15. Forsterite is stable until F = 0.15 in Mer8 and until F = 0.21 in Mer15, after which it disappears. Plagioclase finally crystallizes at F = 0.14 in Mer8 and F = 0.19 in Mer15, when enstatite becomes unstable. The formation of plagioclase is accompanied by a depletion of both Ca and Al in the residual liquids of both Mer8 and Mer15 (Fig. 7c,d). Crystallization models were stopped at F = 0.12, as we do not have experimental constraints at lower F.

Here, we compare our models with the crystallization sequences modeled by Brown and Elkins-Tanton (2009) and Mouser and Dygert (2023). We rescaled their crystallization sequences to our updated MMO depths to directly compare the models (Fig. 9; the unscaled figure is Figure S11 in the Supplementary Material).

The model of Mouser and Dygert (2023) shows a basal dunite layer, as we obtained for our Mer15 composition (Fig. 9b,c). This is caused by the similar starting Mg/Si ratio of Mer15 and the modified silicate composition of the CH chondrite (Fig. 1). The forsterite cumulate in Mer15 is nonetheless more than double the thickness of their basal layer. We also observe that the most evolved products differ substantially (Fig. 9a–c). The lithologies that we obtained based on our experiments include quartz at depths shallower than \sim 140 km–120 km, whereas their model includes quartz-free assemblages at those depths. Our modeled sequence shows that a magma ocean composition with a high starting Mg/Si ratio (> 0.8) can produce a large variety of residual liquid compositions (45 wt %–75 wt % SiO₂ and 40 wt %–0 wt % MgO) through fractional crystallization (Fig. 7).

Although the models of Brown and Elkins-Tanton (2009) were based on an older estimate of the size of Mercury’s core, a few considerations can be drawn from our comparison. First, none of their models produce a basal orthopyroxenite layer as produced from our Mer8 composition (Fig. 9d–g), although the Mg/Si ratios of their starting compositions (0.9–1.1 with one exception) are comparable to ours (0.8–1.2). The only major difference is therefore the absence of S in their models, which allows olivine to appear in all their basal layers because sulfur tends to favor orthopyroxene over olivine (Namur et al., 2016a; Saracino et al., 2025). Second, clinopyroxene saturation occurs at around F = 0.50–0.40 in their models, similar to clinopyroxene occurring at F = 0.40–0.35 in our models. Sulfur speciation has been shown to decrease the activity of CaO in silicate melts, delaying the formation of clinopyroxene (Anzures et al., 2020). Compared to S-free models from the literature, our S-bearing models indicate that this effect might not be as strong as previously suggested.

The different scenarios in the models of Brown and Elkins-Tanton (2009) also suggest that garnet and dense opaque phases (chromite, ilmenite) may have been stable in the MMO. Experiments at 7 GPa on the same compositions used herein (Xu et al., 2024) feature pyrope garnet at low melting fractions of F = 0.20–0.30. Garnet is thus a near-solidus phase and is not expected to fractionate from an evolving MMO near its liquidus. The production of high density, opaque minerals (ilmenite, chromite) at lower liquid fractions (Brown and Elkins-Tanton, 2009) is also not supported by our experiments. The paucity of Cr and Ti at the surface of Mercury (Cartier et al., 2020; Nittler et al., 2011, 2018, 2023; Peplowski et al., 2011), coupled with the more chalcophile behavior displayed by both Cr and Ti in reduced conditions (Cartier et al., 2020; Steenstra et al., 2020b), preclude those elements as being important to mineral formation in the MMO within the P-T-fO₂ range investigated herein.

The lithological evolution of Mercury’s primordial mantle impacts mantle melting conditions and the production of partial melts, eventually influencing the compositions of surface lavas. Indeed, widespread volcanism on Mercury generated its diverse geochemical regions (Head et al., 2011; Weider et al., 2015). These regions have been shown to be the product of different degrees of partial melting of a heterogeneous mantle (Charlier et al., 2013; Namur et al., 2016b; Wang et al., 2022). Charlier et al. (2013) originally proposed different melting sources for two chemically distinct surface regions approximately in the southern hemisphere (termed the ‘intercrater plains-highly cratered terrains’ or IcP-HCT) based on the first XRS data from the MESSENGER spacecraft. They suggested that the two geochemically distinct units G1 (low Al₂O₃ and high CaO, MgO) and G2 (high Al₂O₃ and low CaO, MgO) were the product of low-pressure (< 10 kbar) partial melting of lherzolitic and harzburgitic mantle sources, respectively. As new geochemical data became available (Weider et al., 2015), Namur et al. (2016a) proposed that the geochemical terrains had equilibrated at low pressure (< 15 kbar) following adiabatic decompression of lherzolitic mantles differing by their clinopyroxene fraction. In particular, they found that MgO-rich lavas were generated by a high degree of partial melting (F = 0.35 ± 0.03) starting deeper in the mantle (\sim 360 km depth, corresponding to \sim 4.5 GPa) and at higher temperatures (\sim 1720 ℃) compared to Northern Volcanic Plain lavas (F = 0.27 ± 0.04) that formed at shallower depths (\sim 160 km) and lower temperatures (\sim 1435 ℃). Similarly, Wang et al. (2022) showed that the geochemical region characterized by high Mg/Si (the ‘high Mg region’ or HMR) may have originated from multiple mantle sources characterized by jagged, irregular boundaries, whose melting was likely favored by local thermal anomalies early in the planet’s history (> 4.1 Ga). Melting experiments on a synthetic EH4 composition at 0.5 GPa–5.0 GPa (relevant to the petrogenesis of Mercury’s surface lavas) showed that melting of a single enstatite-rich chondritic composition can produce a variety of silicate melt compositions, from high-Mg melts at high pressure (2 GPa–5 GPa) to silica-rich silicate melts at low pressure (0.5 GPa–1 GPa; Boujibar et al., 2025). They also argued that Mercury’s mantle is mainly pyroxenitic in composition due to the large stability field of orthopyroxene in sulfur-bearing systems.

Studying the primordial mantle of Mercury as produced from a fractionating magma ocean is paramount to understanding the melting sources responsible for the formation of the surface lavas. Importantly, previous melting studies overlooked that certain portions of the primordial mantle would have been more refractory to melting. Our results show that a refractory mantle is produced in the lower cumulate pile, which could be solely forsterite, solely enstatite, or include both phases. In the upper-MMO sequence, a fertile mantle forms when clinopyroxene becomes part of the cumulus assemblage. The formation of clinopyroxene-bearing mantle lithologies starts at around 160 km–190 km depth (Fig. 10). Experimental studies have shown that a refractory mantle made of an olivine + orthopyroxene assemblage would melt congruently between 1680 ℃ at 1.2 GPa and 1780 ℃ at 2.5 GPa (Chen and Presnall, 1975). Conversely, a mantle containing olivine + orthopyroxene + clinopyroxene would start melting between 1310 ℃–1320 ℃ at 1.2 GPa and 1520 ℃ at 2.5 GPa (Walter and Presnall, 1994).

Because of the contrasting fusibility of the refractory and fertile mantles, it is expected that the fertile mantle was the dominant contributor to crustal production via partial melting. Here, we evaluate if the fertile mantle produced by MMO crystallization can actually produce the relatively thick volcanic surface of Mercury by partial melting alone. We consider the thickness of the fertile mantle to be 140 km. We then calculate the volume of the silicate shell corresponding to the fertile reservoir as:

V_{Fertile\ Mantle} = \frac{4}{3}\pi\left( R_{Mercury}^{3} - R_{Core + Refractory\ Mantle}^{3} \right) \tag{11}

where V_{Fertile\ Mantle} is the volume of the fertile mantle, R_{Mercury} is the radius of Mercury (2440 km), and R_{Core + Refractory\ Mantle} is the radius of the core plus the refractory reservoir of Mercury (2300 km). From this, we can calculate the volume of melt produced at a fixed melt fraction F, and the relative thickness it would produce as lavas after extraction and cooling at the surface. We first consider a minimum degree of partial melting (F = 0.30) among estimated values in the literature (Beuthe et al., 2020; Boujibar et al., 2025; Namur et al., 2016a; Wang et al., 2022). Then we can calculate the volume of partial melt produced, and therefore the thickness of the layer that such a melt would produce at the surface (or within the crust; i.e. intrusive magmatism). We find the thickness of the volcanic crust to be 40 km, in good agreement with the estimated crustal thickness of 35 ± 17 km of Padovan et al. (2015). Higher degrees of partial melting (up to 0.70 for the HMR region; Boujibar et al., 2025) have also been suggested. Higher values of F (0.50–0.60) in our models would indeed yield a crustal thickness in the range 60 km–80 km. The range of crustal thicknesses produced by low degrees of partial melting are well within the various estimates of the extent of the Mercurian crust (Beuthe et al., 2020; Padovan et al., 2015). Our results therefore suggest that Mercury’s volcanic crust may have originated from the partial melting of the fertile mantle portion alone. It is important to note, however, that the mantle of Mercury, as formed post-MMO crystallization, might have undergone one or multiple episodes of density-driven mantle overturn (Brown and Elkins-Tanton, 2009; Mouser and Dygert, 2023), which may have refertilized deeper mantle sources capable of producing partial melts by adiabatic decompression melting (Namur et al., 2016a).

SCSS can be calculated at each incremental step of our crystallization model to track the formation of Mg- and Ca-rich sulfides. Studies investigating the MgS-CaS system show that a (Mg,Ca)S solid solution is stable over the MMO’s temperature range (Pitsch et al., 2025). Under reduced conditions, SCSS sensibly decreases with temperature and is affected by the composition of the silicate melt (Namur et al., 2016b). Two end-member scenarios are usually considered when investigating the fate of sulfur in early differentiation processes on Mercury (e.g. Boukaré et al., 2019; Pirotte et al., 2023). In the first, the magma ocean was initially saturated with sulfides, allowing for sulfide phases to form at the CMB, perhaps together with an FeS layer atop the core. In the second, the MMO reached sulfide saturation later (at lower temperature) during crystallization, preventing sulfides from forming during the early stages.

To explore the potential sulfur content of the cumulate pile, we calculated SCSS using the model of Namur et al. (2016b), assuming an average mantle oxygen fugacity of log(fO₂) = IW-5.4 (Namur et al., 2016b) for both Mer8 and Mer15 (Fig. 8). At CMB conditions, Mer15 (SCSS \sim 8.8 wt % at 2040 ℃) has a higher SCSS than Mer8 (SCSS \sim 7.8 wt % at 1960 ℃) due to the higher temperatures and MgO concentrations of Mer15 melts. We use a range of BSMe S concentrations (1 wt %–8 wt % S) and track the onset (fraction of residual liquid) of sulfide crystallization. We treat sulfur as a fully incompatible element in the silicate melt (D_{S}^{crystal/melt} = 0). The concentration of sulfur in the melt can therefore be determined as:

C_{S,Liq} = \ C_{S,0}F^{- 1} \tag{12}

where C_{S,Liq} is the concentration of sulfur in the residual liquid and C_{S,0} is the starting concentration of sulfur in BSMe. When the path of increasing sulfur content with decreasing depth (or the fraction of residual liquid) reaches the modeled SCSS, the magma ocean becomes saturated with sulfide. We observe that for 1 wt % S in the starting BSMe, sulfide saturation is attained at F = 0.33 (\sim 160 km depth) for Mer8 and F = 0.45 (\sim 210 km depth) for Mer15. At 5 wt % S in the BSMe, sulfide saturation occurs closer to the CMB (\sim 350 km depth for Mer8 and \sim 360 km depth for Mer15). At 8 wt % S, the Mer8 MMO would be saturated with S at the onset of crystallization, whereas the Mer15 MMO would only reach sulfide saturation at \sim 450 km–460 km depth (Fig. 8b). At > 8.8 wt % S, the Mer15 MMO would be saturated with S at the CMB.

In Figure 8a,b, we show the potential sulfide content in the cumulate pile. As an approximation, (Mg_0.9Ca_0.1)S is considered as the stable sulfide phase that would form in this fO₂ range (ΔIW < -4.5; Anzures et al., 2020, 2025). Different starting BSMe S contents would affect the amount of (Mg,Ca)S in Mercury’s primordial mantle. For example, in Mer8, 1 wt % S in the starting silicate magma ocean would produce \sim 5.5 wt % (Mg,Ca)S at 160 km depth, whereas 5 wt % S in the magma ocean would crystallize \sim 12.0 wt % sulfides at 350 km depth. Lark et al. (2022) calculated the abundance of Mg- and Ca-rich sulfides in Mercury’s mantle based on the inferred mantle S content of Namur et al. (2016b). They found that 13 wt %–20 wt % (Mg,Ca)S would crystallize from the MMO.

We also evaluate the role that 10 wt % (Mg,Ca)S could play on the evolution of the residual MMO liquid in Figure 7. Sulfide-free residual liquids (blue lines in Fig. 7), after an initial slight increase in Si (up to 25 wt %–28 wt %), become depleted in silica with progressive crystallization in both the Mer8 and Mer15 cases (with Si down to 20 wt %–15 wt %). In contrast, due to decreased proportions of crystallizing silicate minerals, sulfide-bearing liquids (red lines) show enrichments of up to \sim 32 wt %–33 wt % Si under the investigated conditions (Fig. 7a). This is supported by the composition of the experimentally determined S-bearing silicate melts (Fig. 7a). Importantly, the fate of Al and Ca changes dramatically depending on the presence or absence of sulfides (Fig. 7c,d). In accordance with our experiments, Al and Ca concentrations in sulfide-bearing liquids will first increase (up to 7 wt %–10 wt % Al and 2 wt %–4 wt % Ca) because they are hardly consumed by early forming phases (enstatite, forsterite). Then, once plagioclase crystallizes, their concentrations decrease (down to 5 wt %–9 wt % Al and 0 wt %–3 wt % Ca). Conversely, in sulfide-free liquids, we observe an enhanced enrichment of Al and Ca before plagioclase formation (> 10 wt % Al and > 5 wt % Ca). In addition, we also observe that Mg behaves similarly in the sulfide-free MMO as in its sulfide-bearing counterparts (Fig. 7b).

In this study, we modeled the crystallization of MMO considering Mg-dominated sulfides (plus Ca), which is considered to be the main sulfide phase stable in the range of oxygen fugacities investigated in this study (Anzures et al., 2020, 2025). It is, however, important to add that Mercury’s magma ocean might have segregated other sulfides containing Na, K, Ti, and Cr (Anzures et al., 2020, 2025; Mouser and Dygert, 2023). Such exotic sulfides have been also found in aubrites, enstatite-rich achondrites that might represent analogs of Mercury’s surface (e.g. Cartier and Wood, 2019; Wilbur et al., 2022).

Determining the density of silicate melts containing sulfur requires knowledge of the nature of bonding between S and typical silicate melt cations, for which the uncertainties remain large. Silicate melts are considered to dissolve ionic S (e.g. O’Neill, 2021), although several studies have demonstrated the bonding of S with cations like Mg, Ca, Si, and Na in silicate glasses using a variety of analytical approaches (Anzures et al., 2020, 2025; Namur et al., 2016b; Pommier et al., 2023). Some of these analyses, however, were performed on quenched products, which may show differences in speciation and elemental distributions compared to silicate melts (Mysen, 1983). Under reduced conditions, anionic oxygen is partly replaced by sulfur (Fincham and Richardson, 1954); therefore, in terms of molar mass, the density of a silicate liquid containing dissolved S is slightly higher than that of a S-free silicate liquid because sulfur has roughly twice the mass of oxygen (32 vs. 16 g mol ^-1). However, the ionic radius of S^2- (1.84 Å) is larger than that of O^2- (1.41 Å), implying a greater molar volume and thus a lower specific density of sulfur-bearing melt components. Indeed, S has a lower specific density than other typical melt components. At room temperature, the partial specific density (M_i / V_i, with M_i being the molar mass and V_i the molar volume of melt component i) of S (\sim 2060 kg m ^-3) is within the range of the specific densities of other melt components like SiO₂ or K₂O (2000 kg m ^-3–2200 kg m ^-3), and higher than those of common volatiles (H₂O, CO₂ at 900 kg m ^-3–1400 kg m ^-3; Lesher and Spera, 2015). In contrast, if we assume MgS and CaS as melt species in sulfur-bearing silicate melts, the bonding of S with Mg and Ca would decrease the activities, and therefore the molar fractions, of MgO and CaO species in the silicate melt, which have higher specific densities (3300 kg m ^-3–3600 kg m ^-3) compared to S and other sulfides. Furthermore, the occurrence of MgS and CaS melt species will increase, which have specific densities of \sim 2500 kg m ^-3–2700 kg m ^-3 in their solid form (Lark et al., 2022). The densities of these species are higher than dominant melt components like SiO₂ (2200 kg m ^-3) but lower than the specific densities of MgO and CaO.

Therefore, the densities of silicate liquids containing dissolved S should, in principle, slightly decrease. We calculated the effect of sulfur on the density of the MMO liquid (details on the procedure are reported in Section S.12). The resulting densities of the MMO calculated for Mer8 and Mer15 are shown in Figure 11. To our knowledge, this study is the first attempt to explore the impact of S on the density of reduced, Mg-rich silicate melts, although we stress that further studies specifically dedicated to understanding the volumetric properties of S and other sulfides at high temperatures and pressures are needed. Nonetheless, we importantly find that sulfur plays a major role in decreasing the density of silicate liquids. For example, around 10 wt % S dissolved in a silicate liquid would decrease the density by \sim 210 kg m ^-3. This effect may have serious consequences on the ability of newly formed solid phases to float in a vast magma ocean.

To investigate the possibility that minerals formed in the MMO could become buoyant, we compare the densities of the S-bearing residual liquids for both the Mer8 and Mer15 BSMe compositions with the density of solid phases as a function of depth (Fig. 11). Mineral phase densities are calculated with third-order Birch-Murnaghan equations of state in terms of density (Birch, 1947; Vander Kaaden and McCubbin, 2015) using parameters listed in Section S.1 (Table S1). We confirm that graphite has a lower density than any residual MMO melt (brown dash-dotted line in Fig. 11), and may have thus continuously contributed to the formation of the primordial graphite crust. We also observe that quartz, appearing at F = 0.24–0.28, is denser than the residual melts (black dashed line in Fig. 11). Plagioclase is the last mineral appearing in our models at F = 0.14–0.19 (purple dashed line in Fig. 11) and is also shown to be denser than the residual MMO liquid. Although we did not model the very last MMO crystallization products, we expect plagioclase to become more albitic in the last stages of MMO crystallization. We therefore plotted albite in Figure 11 (dark green dashed line), which is also denser than the residual MMO liquid. We also plot potential sulfide minerals (niningerite, MgS; blue dash-dotted line; oldhamite, CaS; orange dash-dotted line) that could have segregated from the MMO (Boukaré et al., 2019; Mouser and Dygert, 2023). For Mer8, since the beginning of MMO crystallization, the density of the MMO liquid is lower than the densities of the sulfides. If the MMO was initially sulfide saturated, we would then expect the sulfides to sink into the cumulate pile with the early formed silicates like enstatite and/or forsterite (last panel in Fig. 8a). Even in the case that sulfide saturation was attained later during MMO crystallization, the newly formed sulfides would still have been denser than the residual liquid (Fig. 8a). As for Mer15, we observe that at CMB conditions, sulfides are slightly positively buoyant with respect to the MMO liquid (Δ\rho^{sulfides,\ MMO} = −3 to −13 kg m ^-3; last panel of Fig. 8b). Gravitational phase segregation depends on several variables such as convective regimes, crystal-melt density contrast, crystal size and crystal fraction, and melt viscosity (Elkins-Tanton, 2012; Mouser et al., 2021; Solomatov, 2015). Sulfide flotation may have been hindered by turbulent convection in the magma ocean (Schmidt and Kraettli, 2022). In the case of Mer15, the density contrast between the sulfides and the magma ocean is too small to effectively segregate sulfides; in this case, sulfides would be trapped amongst the co-crystallizing forsterite minerals.

Considering a MMO of Mer15 composition saturated with sulfides at the onset of crystallization, the sulfides, whether formed before or after the sulfide density crossover, would be stored in the lower and intermediate sections of the mantle. In the alternative case where the MMO Mer15 composition is not saturated with sulfides from the onset of crystallization, sulfides would likely start forming at intermediate depths in the mantle.

As mentioned in Section 6.3 Sulfides in the mantle of Mercury, sulfides other than niningerite and oldhamite may have segregated from the MMO. Sulfides containing Na, K, Cr, Ti may have formed during MMO solidification (Anzures et al., 2025). Their appearance is further corroborated by the presence of natural sulfide phases like troilite, caswellsilverite, heideite, daubréelite, alabandite, and djerfisherite in aubrites (e.g. Keil, 2010; Wilbur et al., 2022). These differentiated meteorites are generally dominated by FeO-free enstatite, minor Na-rich plagioclase, FeO-free diopside, and forsterite (Keil, 2010). Aubrites formed under reduced conditions, which makes them good analogs of Mercury’s surface (Cartier and Wood, 2019). Mineral densities of aubritic sulfides (at 1 atmosphere) are listed in Table 1 of Mouser and Dygert (2023), for which parameters needed for density calculations are lacking. Nonetheless, one important point should be noted. Except for oldhamite and niningerite modeled here (Fig. 11), these sulfides would be much denser (\rho > 3200 kg m ^-3) than the residual MMO liquid over the entire range of pressures investigated here, plotting outside Figure 11. They would likely sink at the bottom of the residual MMO.

The primordial structure of the Mercurian mantle affected the distribution of HPEs like U, Th, and K. The partitioning of these elements between the silicate, metallic, and sulfide phases may have also been influenced by the reduced conditions during Mercury’s differentiation (Boujibar et al., 2019; Pirotte et al., 2023). Elements tend to become more chalcophile and siderophile at increasingly reduced conditions (Boujibar et al., 2019; McCubbin et al., 2012; Pirotte et al., 2023; Steenstra et al., 2020b; Wohlers and Wood, 2017). A FeS layer atop the metallic core may also affect the HPE distribution and heat production (Boujibar et al., 2019; Boukaré et al., 2019), although the existence of such a layer remains debated (Cartier et al., 2020; Pirotte et al., 2023; Smith et al., 2012). Oxygen fugacity and the initial S content of the MMO have been shown to affect the storage of U in the Mercurian mantle (Boukaré et al., 2019). At low initial S contents (\sim 1 wt %) and relatively oxidizing conditions (IW-4 to IW-2), sulfides remain in sufficient abundances that HPEs could be preferentially hosted in a sort of KREEPy (K, REE, and P-rich) layer in the upper mantle or, if present, in the FeS layer at the CMB. Under more reducing conditions (IW-6), sulfides are segregated in the shallow mantle of Mercury (< 96 km depth), making this region richer in HPEs. At higher initial S contents (\sim 5 wt % S) and oxidizing conditions (> IW-4), SCSS is low, and an FeS layer would form at the CMB, storing important amounts of HPEs (Boukaré et al., 2019). At high S contents and reduced conditions, HPEs would mostly be concentrated in the intermediate and upper mantle (< 220 km depth; see Fig. 2 in Boukaré et al., 2019).

Importantly, large uncertainties surround the HPE distribution between Mg- and Ca-rich sulfides and the silicate melt (e.g. Boukaré et al., 2019). Only Pirotte et al. (2023) recently determined the partition coefficients of HPEs between MgS and the silicate melt under reduced conditions (IW-8.5), showing the incompatible behavior of U and Th. Differences in the lattice structure between Mg- and Ca-bearing sulfides might affect the incorporation of HPEs, and particularly U and Th (Blundy and Wood, 2003). Indeed, Ca-rich sulfides have been shown to be significant U and Th carriers in enstatite achondrites (Murrell and Burnett, 1982). Rare earth elements (REEs) also partition into oldhamites found in aubrites and enstatite chondrites (Dickinson and McCoy, 1997; Ingrao et al., 2019). This hints at a potentially higher partitioning of HPEs in Ca-rich sulfides. Here we quantify HPE partitioning in the MMO for both Mer8 and Mer15 compositions considering two end-member scenarios: (i) incompatible HPEs, using D_{U}^{MgCaS/Melt} = 0.5, D_{Th}^{MgCaS/Melt} = 0.1 from Pirotte et al. (2023), and considering for D_{K}^{MgCaS/Melt} an average (0.05) from the range found in Pirotte et al. (2023); (ii) compatible U and Th, assuming D_{U}^{MgCaS/Melt} = D_{Th}^{MgCaS/Melt} = 2, following the case of moderately high partition coefficients in Boukaré et al. (2019, here K was not considered as we believe it would still behave very incompatibly). The partitioning of element i between a solid phase and the silicate melt is expressed by the partition coefficient D_{i}, as: D_{i}^{Crystal/Melt} = C_{i}^{Crystal}/C_{i}^{Melt} \tag{13}

where C_{i}^{Crystal} is the concentration of element i in a given crystal phase, and C_{i}^{Melt} is the concentration of element i in the melt. The concentrations of trace elements, like HPEs, in the MMO are calculated following Rayleigh’s equation for fractional crystallization (Eq. 1). We calculate the distribution of HPEs based on the estimated (Mg,Ca)S content as shown in Figure 8a,b. Partition coefficients for the silicate minerals are from White (2020). The starting U, Th, and K concentrations in the BSMe are from Pirotte et al. (2023). The concentrations of U, Th, and K in the MMO as a function of depth are reported in panel (a) of Figure 12. Furthermore, after calculating the concentration of HPEs in the evolving MMO liquid, we could also determine the HPE abundances in the cumulate assemblages. The results are reported in Figure 12b. For clarity, we only show Mer8 (Fig. 12). Plots for Mer15 are shown in Section S.13 (Fig. S13).

In scenario (i: incompatible HPEs) U, Th, and K are all incompatible, and their concentrations increase in the residual silicate liquid as the cumulate grows in thickness (red lines in Fig. 12a). We also show sulfide-free magma ocean models for comparison (grey lines in Fig. 12). The residual sulfide-bearing liquids are slightly more depleted than compared to the equivalent corresponding sulfide-free liquids, given the more compatible behavior displayed by U and Th (Fig. S14). Silicate liquids trapped in the growing cumulate might affect the trace element budget of the residual magma ocean liquid (Elkins-Tanton et al., 2011; Snyder et al., 1992). We therefore modeled the HPE budget considering 5% trapped liquid in the cumulate assemblages (dark green lines in Fig. 12). HPEs are still enriched in the residual MMO liquid as crystallization proceeds. We only observe slightly lower HPE abundances compared to sulfide-bearing liquids with no trapped melt. In scenario (ii: compatible HPEs), U and Th are still enriched in the residual MMO liquid (orange and light green lines in Fig. 12a), although to a lesser extent as compared to scenario (i), given the higher bulk partition coefficients of the sulfide-bearing cumulates (Fig. S14a). The most noticeable difference between the two scenarios is found in the HPE abundances in the cumulates (Figs 12b and S13b for Mer15). We show an increased content of U (4 ppb–8 ppb in sulfide-bearing Mer8 with no trapped liquid; 8 ppb–10 ppb in sulfide-bearing Mer8 with trapped liquid; 2 ppb–9 ppb in sulfide-bearing Mer15 with no trapped liquid; 6 ppb–10 ppb in sulfide-bearing Mer15 with trapped liquid) and, particularly, Th (32 ppb–36 ppb in sulfide-bearing Mer8 with no trapped liquid; 37 ppb–42 ppb in sulfide-bearing Mer8 with trapped liquid; 11 ppb–36 ppb in sulfide-bearing Mer15 with no trapped liquid; 24 ppb–40 ppb in sulfide-bearing Mer15 with trapped liquid) as compared to scenario (i) in the lower and intermediate mantle of Mercury (< 150 km–200 km depth) caused by both the higher compatible behavior of HPEs in the sulfides coupled with the higher content of sulfides calculated in the lower portions of the mantle (Fig. 8).

In general, considering both scenarios, we observe that the residual liquid of the MMO would likely contribute to forming a shallow layer that is rich in incompatible elements (Fig. 12a, Fig. S13a). A KREEP layer in the upper mantle of Mercury is therefore hypothesized, similar to the Moon (Warren, 1988; Warren and Wasson, 1979; Zhang et al., 2024). On Mercury, such a layer would be found between the late cumulates and the primary graphite crust. As regards the cumulates, we notice that the presence of sulfides combined with that of trapped melt have an effect on HPE distribution in the cumulates (Fig. 12b, Fig. S13b). Although we observe that sulfides alone will only cause a small U, Th enrichment in early cumulates (roughly a 2 ppb–5 ppb difference), the presence of trapped liquid will cause a steeper U, Th enrichment, especially in the late cumulate sequence (< 200 km depth). K is almost absent in early formed cumulates (only trapped liquid has a minor effect), and it will only be stored when plagioclase saturates.