). Hydrous silicate melts and the deep mantle H2O cycle. Earth and Planetary

We report ab initio atomistic simulations of hydrous silicate melts under deep upper mantle to shallow lower mantle conditions and use them to parameterise density and viscosity across the ternary system MgO-SiO 2 -H 2 O (MSH). On the basis of phase relations in the MSH system, primary hydrous partial melts of the mantle have 40-50 mol% H 2 O. Our results show that these melts will be positively buoyant at the upper and lower boundaries of the mantle transition zone except in very iron-rich compositions, where (cid:2) 75% Mg is substituted by Fe. Hydrous partial melts will also be highly inviscid. Our results indicate that if melting occurs when wadsleyite transforms to olivine at 410 km, melts will be buoyant and ponding of melts is unexpected. Box models of mantle circulation incorporating the upward mobility of partial melts above and below the transition zone suggest that the upper mantle becomes eﬃciently hydrated at the expense of the transition zone such that large differences in H 2 O concentration between the upper mantle, transition zone and lower mantle are diﬃcult to maintain on timescales of mantle recycling. The MORB source mantle with ∼ 0.02-0.04 wt% H 2 O may be indicative of the H 2 O content of the transition zone and lower mantle, resulting in a bulk mantle H 2 O content of the order 0.5 to 1 ocean mass, which is consistent with geochemical constraints and estimates of subduction ingassing. © 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
H 2 O exerts a powerful control on the location and composition of partial melts in Earth's deep interior by dramatically reducing the mantle solidus (Hirschmann, 2006).Although relatively insoluble in upper and lower mantle minerals, H 2 O has a high solubility in wadsleyite and ringwoodite in the transition zone (410-670 km) providing the capacity for storing ∼3 ocean masses of H 2 O in this reservoir alone (Bolfan-Casanova et al., 2006;Fei and Katsura, 2020;Hirschmann, 2006;Hirschmann et al., 2005;Ohtani, 2015) (although hydrogen substitutes typically as H + or OH − , hereafter we refer to H 2 O for simplicity).Subducted, H 2 Orich lithologies in oceanic crust undergo dehydration reactions in the shallow upper mantle (< 150 km) leading to the generation of buoyant, hydrous fluids or partial melts responsible for arc volcanism (Poli and Schmidt, 2002).However, in the mantle portion of cooler slabs a significant fraction of the H 2 O dissolved in both hydrous and nominally anhydrous phases may survive to greater depths (Iwamori, 2004;van Keken et al., 2011).Stalling and stagnation of slabs in the transition zone may lead to dehydration of serpentinized mantle in cooler slabs providing a mechanism for deep mantle hydration (Iwamori, 2004;Komabayashi and Omori, 2006), a process implicated in deep focus earthquakes and consistent with the petrology of some sublithospheric diamonds from the deep transition zone or shallow lower mantle (Omori et al., 2004;Pearson et al., 2014;Shirey et al., 2021).
The H 2 O content of the transition zone is not well constrained but has been estimated from geophysical observations.H 2 O contents ranging from <1000 ppm to as much as 0.5 wt% have been suggested on the basis of electrical conductivity constraints (Huang et al., 2005;Karato, 2011;Sun et al., 2015b;Yoshino et al., 2008;Zhang et al., 2021).Global seismic data have also revealed a range from relatively dry to substantial hydration (Houser, 2016;Meier et al., 2009;Wang et al., 2019).For example, modelling of globally averaged seismic velocities and density based on mineral elasticity indicate about 1 wt% H 2 O in the transition zone just below the https://doi.org/10.1016/j.epsl.2022.1174080012-821X/© 2022 The Authors.Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).410 km discontinuity for a pyrolytic mantle composition with 60% olivine/wadsleyite, or a relatively dry (<0.1 wt% H 2 O) transition zone and upper mantle with 50% olivine/wadsleyite (Wang et al., 2019).
If the transition zone is hydrated, then H 2 O can potentially be released at its upper and lower boundaries as hydrous wadsleyite and ringwoodite transform to relatively anhydrous olivine or bridgmanite + ferropericlase during mantle upwelling or downwelling, respectively, leading to the release of H 2 O and the generation of hydrous silicate melts in these regions (Hirschmann, 2006;Panero et al., 2020).This process is the basis of the "transition zone water filter" model, which describes how dense, hydrous silicate partial melts produced by the dehydration of material upwelling through the upper transition zone boundary (410 km depth) may pond and be reabsorbed into the transition zone, efficiently stripping the upwelling mantle of its incompatible trace elements (Bercovici and Karato, 2003;Karato et al., 2020).
Low seismic velocity regions above 410 km and beginning at depths near 350 km have been observed both regionally and globally and have been interpreted as arising from partial melting and potentially indicating neutrally or negatively buoyant melt (Revenaugh and Sipkin, 1994;Song et al., 2004;Tauzin et al., 2010).Conversely, global three-dimensional models coupling shear attenuation and seismic velocity do not show evidence of a global partial melt at depths of 350 km (Debayle et al., 2020).Recent numerical models indicate that a hydrous transition zone containing at least 0.2-0.3wt% H 2 O can explain regions of low seismic velocity observed above and below the Pacific slab by pockets of hydrous melt that are in turn responsible for intraplate and petitspot volcanism in northeast China and offshore Japan (Yang and Faccenda, 2020).Slow seismic velocities have also been seen regionally at about 700-800 km beneath N. America and interpreted as hydrous partial melting (Schmandt et al., 2014).
The gravitational stability of partial melt is dependent on the density contrast between the melt and solid mantle, whereas the ability of melt to migrate depends in part on its viscosity, and both properties are strongly influenced by changes in melt structure at depth.However, the structure and properties of hydrous melts under the high pressure ( P ) and high temperature (T ) conditions above and below the transition zone are poorly constrained, limiting our understanding of their mobility.The density and viscosity of various hydrous silicate melts have been reported at mantle pressures using sink/float (Agee, 2008;Jing and Karato, 2012;Matsukage et al., 2005;Sakamaki et al., 2006), x-ray absorption contrast (Malfait et al., 2014;Sakamaki, 2017;Sakamaki et al., 2009;Seifert et al., 2013), and falling sphere viscometry (Persikov et al., 2017;Poe et al., 2006) but data under transition zone and shallow lower mantle conditions and over a range of appropriately high H 2 O concentrations are sparse.
Ab initio molecular dynamics (MD) simulations, based on density functional theory (DFT), can provide highly accurate models of the atomistic structure and dynamics of liquids and their macroscopic properties (Jahn, 2022).In an early application of DFT to silicate melts, Stixrude and Karki (2005) used ab initio MD with a local-density approximation (LDA) functional to predict the Mie-Grüneisen equation of state, melting curve, and local structure of anhydrous liquid MgSiO 3 .The results predicted a continuous increase in the mean Si-O coordination number from 4 to 6 over the mantle pressure regime resulting in an approximately 5-fold reduction in the liquid-solid density contrast to ∼4% at the coremantle boundary.Wan et al. (2007) reported ab initio MD simulations of MgSiO 3 using a generalised gradient approximation (GGA) functional.While predicting a similar liquid-solid density contrast as Stixrude and Karki (2005), the liquid equation of state is systematically offset by 14 GPa due to the inherent tendency of DFT calculations to overestimate (GGA) or underestimate (LDA) volume.Wan et al. (2007) predicted that negatively buoyant melts from (Mg,Fe)SiO 3 perovskite compositions may reside at deep lower mantle conditions.Subsequent ab initio MD simulations reported by de Koker et al. (2008) predicted a solid-liquid density crossover for anhydrous Mg 2 SiO 4 at 13 GPa, corresponding approximately to the upper boundary of the mantle transition zone.
H 2 O dissolution in silicate melts proceeds by depolymerisation of the melt structure via reaction with bridging oxygen atoms to form hydroxyl (OH − ) units resulting in a profound reduction in viscosity (Kohn, 2000).Various ab initio MD simulations have been reported for hydrous silicate melts under a wide range of P -T conditions (Mookherjee et al., 2008;Karki and Stixrude, 2010a;Karki et al., 2010;Bajgain et al., 2019;Yuan et al., 2020;Karki et al., 2020;Solomatova and Caracas, 2021;Karki et al., 2021).Despite the lower density of hydrous melts relative to anhydrous compositions, Mookherjee et al. (2008) argued that a buoyantly stable melt at the upper boundary of the transition zone would contain approximately 3 wt% H 2 O. Very recently, Karki et al. (2020Karki et al. ( , 2021) ) reported calculations of hydrous MgSiO 3 and Mg 2 SiO 4 melts with up to 25 mol% FeO and predicted that partial melts containing a few weight per cent of H 2 O may be gravitationally stable above and below the mantle transition zone.However, as noted by Yuan et al. (2020), even with iron enrichment, the critical H 2 O content (∼4 wt%) required for neutral buoyancy of hydrous silicate melts relative to pyrolytic solid mantle is much lower than expected for small-degree mantle melts from experimental phase relations (section 1.1), with the consequence that natural hydrous melts cannot experience long-term gravitational stability above, below or within the transition zone.

Melting phase relations and primary hydrous melts at transition zone depths
Shown on Fig. 1 are phase relations in the MSH system at 13 GPa on the basis of experiments and thermodynamic modelling (Myhill et al., 2017;Novella et al., 2017).Partial melts from hydrous peridotite are constrained to the boundary curve where forsterite, enstatite and melt coexist, and for an estimated mantle geotherm at 410 km of ∼1550 ± 50 • C (Katsura et al., 2010), hydrous partial melts will have ∼40 to 50 mol% H 2 O.Such large amounts of H 2 O indicate complete solubility between silicate melt and hydrous fluid at this pressure, consistent with results at 1000-1250 • C for forsterite-rich compositions in the MSH system at 13.5 GPa (Melekhova et al., 2007) (Fig. 1), where melts have 70-80 mol% H 2 O and are beyond the second critical end-point (Fig. 1).The phase relations shown on Fig. 1 are also generally consistent with the experiments for hydrous partial melting of peridotite of Kawamoto and Holloway (1997) and Kawamoto (2004), with melts from those studies generated at 900-1400 • C consistent with the extension of the forsterite+enstatite+melt boundary curve.Also shown in Fig. 1 are the bulk compositions used by previous workers to experimentally assess the densities of a range of hydrous silicate melt compositions to explicitly explore the mobility of hydrous melts around 410 km (Freitas et al., 2017;Jing and Karato, 2012;Matsukage et al., 2005;Sakamaki et al., 2006).In all cases these workers find that model hydrous melt compositions are negatively buoyant relative to the density of PREM at 410 km, and on this basis predict that melts should remain perched at that depth or sink and be reabsorbed back into the transition zone.However, the phase relations of Myhill et al. (2017) at 13 GPa indicate that these chosen melt compositions are either considerably lower in H 2 O than expected in a melt from a hydrous peridotite at 1550 • C (red star in Fig. 1) (Jing and Karato, 2012;Matsukage et al., 2005;Sakamaki et al., 2006) or are far too iron-rich to be a plausible mantle melt (Freitas et al., 2017).Myhill et al. (2017).Red circles show the compositions investigated in this study.The red star shows the composition of partial melt at 1550 • C in equilibrium with model mantle peridotite (green star) based on the mantle geotherm at 410 km of Katsura et al. (2010).Blue diamonds show the compositions used by Jing and Karato (2012) and purple diamonds are from Matsukage et al. (2005) in the CMASF-H 2 O system projected onto MSH (with MgO calculated as total MgO+FeO).Green diamonds are the experimental bulk compositions of Sakamaki et al. (2006) and the yellow diamond shows the extrapolated 0.7% melt fraction model melt of Freitas et al. (2017).Blue, cyan and green squares show the melt compositions used in previous DFT calculations (Karki et al., 2021;Mookherjee et al., 2008;Yuan et al., 2020).The green field shows hydrous silicate melts at 13.5 GPa reported by Melekhova et al. (2007), and the blue field shows melt compositions from hydrous peridotite at 11-24 GPa and 1100-1400 • C (Kawamoto, 2004;Kawamoto and Holloway, 1997) with H 2 O contents estimated from deficiency in microprobe totals.(For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.) In this paper we report the results of ab initio MD simulations of hydrous silicate melts over a comprehensive range of compositions in the MSH ternary system at high temperatures (1800 K, 3000 K) and over a range of pressures (∼1-32 GPa).We parameterise melt density and viscosity as a function of P , T and composition across the entire MSH phase diagram and compare our results with experiments and previous ab initio MD simulations.With reference to the melting phase relations of Myhill et al. (2017), we elucidate the mobility of primary hydrous melts at transition zone depths and discuss implications for the deep mantle H 2 O cycle.

Method
Ab initio MD calculations based on DFT were performed using VASP (Kresse and Furthmüller, 1996a,b).The electronic interactions were described by the projector-augmented wave (PAW) pseudopotentials (Blöchl, 1994;Kresse and Joubert, 1999) using the Perdew-Burke-Ernzerhof (PBE) formulation of the GGA exchange correlation functional (Perdew et al., 1996), and with an energy cut-off of 550 eV sampling the Brillouin zone at the -point.MD trajectories were calculated in the canonical (NVT) ensemble with periodic boundary conditions and a Nosé thermostat (Nosé, 1984).
The final 40 ps of the equilibrated trajectories were used to analyse the melt structure and physical properties.The radial distribution functions and mean square displacements for each system at all P-T conditions were examined to confirm the simulations were representative of the liquid state.The anhydrous configurations were also equilibrated for 15 ps at 1850 K (the enstatite solidus) and 2163 K (the forsterite solidus) for direct comparison with synchrotron x-ray diffraction measurements on these liquids at ambient-P (see the supplementary material).The simulation results are in good overall agreement with the experimental mea- surements (Fig. S2) demonstrating the efficacy of ab initio MD for accurate computation of the atomistic structure of silicate melts.
To account for the under-binding of chemical bonds inherent to the PBE-GGA functional employed in this study, pressures computed at given densities were corrected by comparison between our GGA simulations and additional calculations of anhydrous En and Fo, and hydrous MSH2 melts using the Strongly Constrained and Appropriately Normed (SCAN) meta-GGA functional (Sun et al., 2015a), which provides more accurate volumes (Chen et al., 2017).The SCAN calculations are in good agreement with the equations of state previously reported for liquid MgSiO 3 (Fig. S1a) or Mg 2 SiO 4 (Fig. S1b) at 3000 K from experiment (Mosenfelder et al., 2009;Petitgirard et al., 2015), ab initio MD (de Koker et al., 2008;Stixrude and Karki, 2005), and classical MD (Lacks et al., 2007).An empirical pressure correction, P emp = −2.5 − 0.13P , is required to bring the PBE-GGA results into agreement with the computed SCAN values and anhydrous literature data.The PBE-GGA MSH2 melt data at 3000 K (Fig. S1c) and 1800 K (Fig. S1d) adjusted using the same empirical correction are consistent with the corresponding SCAN values, demonstrating the efficacy of using this correction for the full range of melts simulated in this study.
The Supplementary Material contains details of all the simulations in this study as well as those of En from Stixrude and Karki (2005) and Fo from de Koker et al. (2008) that are included in the parameterisations of density and viscosity described in sections 3.1 and 3.2.The code and the data used to produce the density and viscosity parameterisations can be accessed in the form of an online and interactive Jupyter notebook here: https://mybinder.org/ v2 /gh /olivertlord /hydrous -melts /main ?urlpath =apps /hydrous _melt _ PVTX .ipynb.The source code and data files can be downloaded for offline use here: https://github .com/olivertlord /hydrous -melts.

Results
A full description of the melt structure computed from the simulation trajectories is provided in the Supplementary Material.Consistent with previous studies (Mookherjee et al., 2008;Karki and Stixrude, 2010a;Karki et al., 2010;Bajgain et al., 2019;Yuan et al., 2020;Karki et al., 2020;Solomatova and Caracas, 2021;Karki et al., 2021), H 2 O dissolution proceeds via depolymerisation of the melt structure with little influence on the Si and Mg local coordination environment.Dissolved H 2 O leads to a substantial reduction in both melt density and viscosity, with the hydrous melts highly inviscid over all P -T space.

Melt density
The densities of the simulated melts are shown as symbols in Fig. 2 while the lines describe a global fit to all of the simulations reported in this study as well as those for En from Stixrude and Karki (2005) and Fo from de Koker et al. (2008).During regression, the pressure at each measured V -T -X point is calculated using the Mie-Grüeneisen-Debye (MGD) equation of state as implemented in the Burnman thermodynamic toolkit (Cottaar et al., 2014).The MGD model has six variable parameters: V 0 , the volume at the reference temperature (here 292 K); K 0 , the bulk modulus and K , its derivative with respect to P ; γ 0 , the Grüeneisen parameter at ambient P and an exponent, q, that describes its volumetric dependence, and finally θ 0 , the ambient-P Debye temperature.The value of θ 0 is fixed in our fit at 812 K, the value for pure enstatite at ambient conditions (Yang and Ghose, 1994).This assumption is made because of the lack of information on θ 0 for silicate liquids and has a minimal impact on the quality of the fit.Left to vary, this parameter tends to 0 K, which is physically unreasonable.The five remaining variables are allowed to vary linearly with MgO content (x); there is insufficient data along the MgO-SiO 2 join to require a higher order dependence.The parameters V 0 and q are also allowed to vary linearly with the H 2 O content ( y), so that, for example (1) However, it is necessary to allow the remaining variables K 0 , K and γ 0 to vary quadratically with H 2 O content, so that, for example for a total of 18 fitted parameters.We have tested all 32 possible combinations of linear and quadratic dependence for the five variable MGD parameters as a function of H 2 O content; this combination is the one that produces the most statistically reasonable fit with the lowest reduced chi squared statistic (χ 2 v ).Because of the strong cross-correlation between equation of state variables, other combinations have very similar χ 2 v and the choice has little effect on the computed densities (see Supplementary Material).
Our fit was performed using orthogonal distance regression (ODR) as implemented in the SciPy Python library to take into account uncertainties in both P and T (there is no uncertainty in V as it is fixed in the ab initio MD simulations).The regressed parameters are listed in Table 2; the average residual is 4.5%.Across the wide range of compositions investigated, all hydrous and anhydrous melts are less dense than the solid mantle, as estimated from the Preliminary Reference Earth Model (PREM) (Dziewonski and Anderson, 1981) or any other 1-D seismic model (see Supplementary Material).They are also less dense than the anhydrous solid residue in the MgO-SiO 2 system after extraction of the model melt composition shown in Fig. 1 from a bulk mantle composition modelled after the depleted MORB mantle composition of Workman and Hart (2005).

Melt viscosity
Shear viscosities (η) were computed using the Green-Kubo relation formulated in terms of the integral of the stress autocorrelation function (ACF) where σ denotes the mean of the off-diagonal components (σ xy , σ yz , σ zx , (1/2)[σ xx − σ yy ], (1/2)[σ yy − σ zz ]) of the traceless stress tensor, k B is the Boltzmann constant, t is time, and the angled brackets denote an ensemble average for different time origins t 0 (Cui et al., 1996).All ACFs converge to zero within the equilibrated simulation timescale (40 ps) with convergence times of ∼1 ps at low P and ∼10 ps at the highest P and lowest T (Fig. S4).The statistical error of the viscosity was estimated from the standard deviation of the values obtained for the five independent stress autocorrelation functions obtained from the individual stress tensor components.
The viscosities for the full range of simulated melts are shown as symbols in Fig. 3 while the lines represent a global fit to all of the simulations reported in this study.During regression, the viscosity (η) at each measured V -T -X point is calculated using the modified Vogel-Fulcher-Tammann-Hesse (VFTH) function: where and x denotes the fraction of the oxide components MgO (M), SiO 2 (S) and H 2 O (H) in mol%.Again, ODR was used to take into account uncertainties on P , T and η.The regression coefficients a i,M , a i,S , a i,H are provided in Table 3.The average residual is 19%.
Experimental measurements reveal a reduction in isothermal viscosity of polymerised melts before encountering a viscosity turnover at ∼3-6 GPa, which has been attributed to high compressibility up to the tetrahedral packing limit (Tinker et al., 2004;Wang et al., 2014).An anomalous negative P -dependence in viscosity has also been reported up to 13 GPa in depolymerised melts including MgSiO 3 (Cochain et al., 2017).Previous classical (Lacks et al., 2007) and ab initio MD calculations (Karki and Stixrude, 2010b) predicted an anomalous reduction in viscosity on initial compression of liquid MgSiO 3 at 3000 K, reaching a minimum close to 5 GPa, although the magnitude of the viscosity reduction is of the same order as their reported statistical uncertainty.No such downturn in viscosity is observed for the more depolymerised Mg 2 SiO 4 melt (Lacks et al., 2007;Ghosh and Karki, 2011).Due to the statistical uncertainty in viscosity obtained from the ACF, our results for P -dependence of the viscosity of liquid MgSiO 3 on initial compression are inconclusive.However, we do observe a negative P -dependence in viscosity of liquid MgSiO 3 calculated from the equilibrium diffusion coefficient by using the Stokes-Einstein relation (equation S1) and assuming an O 2− diffusion mechanism (Fig. S5), with viscosity values comparable to previous studies (Lacks et al., 2007;Karki and Stixrude, 2010b).We also observe a low-P viscosity minimum at ∼ 5 GPa for the anhydrous MS melt, which has a similar degree of polymerisation to MgSiO 3 .By comparison, we see no evidence of a viscosity turnover in the more depolymerised Fo melt composition or for any of the hydrous melts, all of which experience a continuous increase in viscosity with increasing P consistent with the free volume theory (Schmelzer et al., 2005).

The effect of iron on the density of hydrous melts
As demonstrated in Fig. 2, all simulated melt compositions in the MgO-SiO 2 -H 2 O ternary are considerably less dense than the mantle as estimated from the PREM model and remain so throughout the upper mantle and transition zone and into the shallow lower mantle.However, densities of melts in the MSH system will be modified considerably when the full complement of elements are considered, especially by the addition of incompatible elements which preferentially enter the liquid phase during partial melting.Iron has a particularly strong influence on melt density due to its high mass compared to magnesium for which it substitutes, and existing data indicate that hydrous melts generated at depth from mantle peridotite will be similar to or enriched in iron relative to the solid mantle (Kawamoto, 2004;Kawamoto and Holloway, 1997;Kushiro and Walter, 1998;Mibe et al., 2006).For example, in the studies of hydrous partial melting of peridotite by Kawamoto and Holloway (1997) and Kawamoto (2004), at pressures of 11 to 24 GPa the Mg# of partial melts range from ∼80 to 90.
On the basis of the experimentally determined partitioning behaviour of iron between olivine and melt at high pressures that considers the effect of melt composition (Mibe et al., 2006), Jing and Karato assume a K D Fe/Mg value of 0.25 to calculate the FeO content used in their studied melt compositions (Fig. 1), yielding melt Mg#s of 61 to 70 that are close to but somewhat more ironrich than the melts in the studies by Kawamoto.In contrast, Freitas et al. (2017) suggest that small degree melts will be much more iron rich.Using their reported analyses of 2 to 25 vol% partial melts from hydrous peridotite, they extrapolate to a melt fraction of 0.7 vol% (Fig. 1) that was chosen to match the shear wave velocity at 410 km to estimate a melt FeO content of ∼33 wt%.
From the compositional information provided, their preferred melt would have ∼1 wt% MgO, yielding an Mg# of ∼5.If such a melt composition were in equilibrium with mantle olivine (Mg#90) this would yield an olivine-melt K D Fe/Mg of < 0.01, which we suggest is entirely unrealistic on the basis of previous experimental data (Kawamoto, 2004;Kawamoto and Holloway, 1997;Kushiro and Walter, 1998;Mibe et al., 2006).We also note that the measured Na 2 O (and Al 2 O 3 ) contents of their partial melts become progressively lower as the degree of melting decreases (Na 2 O = 0.07 wt% at a melt fraction of 2% but is 0.3 wt% in the bulk), opposite to the expectation for an incompatible element and indicative of an erroneous melt composition measurement.
To investigate the effect of FeO substitution on the melts simulated in this study and to compare with experimental data, and in particular to evaluate the potential for a liquid-solid density inversion at the boundaries of the transition zone, densities for all melt compositions were re-calculated based on the difference in atomic mass of Fe and Mg in the system (Mg 1-x Fe x )O-SiO 2 -H 2 O for x = 0.25, 0.5 and 0.75.To confirm the validity of this simple approximation, additional benchmark ab initio MD calculations of Fe-bearing ENH2 composition melts, with all Mg atoms replaced by Fe, were performed at 1800 K using PBE (with and without spin-polarization) and the Hubbard (PBE+U) correction.Regardless of the type of simulation, computed densities are consistent with a simple substitution of the mass of Mg by Fe (Fig. S6).These findings are consistent with previous ab initio MD simulations of Fe-bearing MgSiO 3 liquids in which the valence and spin state of Fe were found to contribute < 2% of changes in melt density at 25% Fe content (Karki et al., 2018).There are some distinct differences in the local coordination environment of Mg and Fe with Fe adopting a more traditional network forming role compared to Mg (Fig. S7).However, these differences in local structure have an inconsequential effect on the bulk density of the melts by comparison to the large effect of atomic mass.Substitution of Mg by Fe has only a minimal influence on the bulk melt diffusive properties (Fig. S8), where the compositional dependence on viscosity is dominated by H 2 O dissolution.2006) and Matsukage et al. (2005).Labels on the X-axis refer to the bulk compositions used in the studies.(b) Comparison of calculated melt densities in the MgO-SiO 2 -H 2 O system from this study (black and red curves) with melt densities from previous DFT-MD simulations on a range of bulk compositions (black and red circles; Karki et al., 2021;Mookherjee et al., 2008;Yuan et al., 2020); bulk compositions are provided on the figure and shown in Fig. 1).

Comparison with previous experimental and theoretical density data on hydrous melts
Fig. 4a shows the density of hydrous silicate melts determined using sink/float experiments as reported in Jing and Karato (2012), Matsukage et al. (2005) and Sakamaki et al. (2006) compared to calculated melt densities using our results for equivalent compositions projected onto the MgO-SiO 2 -H 2 O ternary at the experimental conditions and for the appropriate substitution of FeO for MgO.Five of the ten melt densities are reproduced to better than 2% relative and all are reproduced to better than 6% relative.All of the Matsukage et al. (2005) data are well reproduced (<2%), and the poorest correspondences occur with the melt measured at 14.3 GPa by Jing and Karato (2012) and with the melt by Sakamaki et al. (2006) at 16.8 GPa, 2573 K.It is not surprising that the melt densities of Sakamaki et al. (2006) are the most poorly reproduced because these are in the silica saturated region of the ternary phase diagram where we have no simulations and rely on the extrapolation of our PVTX equation of state (see Fig. 1).On the basis of the phase relations, however, these melts are unlikely to be representative of melts produced in a hydrous mantle at transition zone depths.
Fig. 4b shows the density of hydrous silicate melts calculated previously using DFT-MD for H 2 O-rich compositions along the Mg 2 SiO 4 -H 2 O (Yuan et al., 2020) and MgSiO 3 -H 2 O (Mookherjee et al., 2008;Karki et al., 2021) joins compared to calculated melt densities using our results for equivalent compositions.Our calculated densities are ∼3 to 6% lower than those of Mookherjee et al. (2008) but higher by similar amounts than those of Karki et al. (2021) for hydrous MgSiO 3 melts.These discrepancies are not unexpected given the differences in DFT methodology and corrections applied.Mookherjee et al. (2008) performed LDA simulations, which are known to be less accurate compared to GGA for hydrous systems.The slight increased deviation from our equation of state at higher pressures may be attributed to the P -independent correction to account for the inherent LDA over-binding error.Karki et al. (2021) performed GGA simulations but chose a low plane-wave energy cut-off of 400 eV.Both the Pulay stress and inherent GGA under-binding error are under-corrected, as evidenced by systematic underestimates of the melt density compared to experimental values for anhydrous melts.In contrast, we calculate densities for hydrous Mg 2 SiO 4 melts that are very close to those of Yuan et al. (2020), within 2% in all cases.This is somewhat surprising given that these authors apply no empirical correction to their computed pressures despite their low plane-wave energy cutoff of 450 eV.
Overall, we find good correspondence between our parameterized melt densities and those produced by experiment and theory, with nearly all comparisons within 5%, and most within 3%, relative to our calculated values.

Evaluation of hydrous melt density at the upper and lower transition zone boundaries
The density across a wide range of hydrous melt compositions with Mg#s from 100 to 25 is shown in Fig. 5 at 1800 K and at pressures approximating the upper (13 GPa) and lower (24 GPa) transition zone boundaries.In the iron-free system only the most silica-rich and H 2 O-poor melts are more dense than PREM.At a melt Mg# of 75, only melt compositions with less than about 20 mol% H 2 O are more dense than ambient mantle.Even at a melt Mg# of 50, the estimated 410 km equilibrium melt composition (star on Fig. 5) remains considerably less dense than the ambient mantle.It is not until a melt Mg# of ∼25 that the estimated 410 km melt composition becomes neutrally buoyant.Adopting a solid/melt K D Fe/Mg value of 0.25 for the 410 km model melt composition yields a Mg# of about 70, still far too Mg-rich to be more dense than ambient mantle as defined by PREM at 410 km.While the phase relations for hydrous peridotite are less well known at higher pressures near the base of the transition zone and upper part of the lower mantle, the experiments of Kawamoto (2004) indicate that the H 2 O contents could be as high as in lower pressure melts (Fig. 1) and with similar iron contents.Hydrous melts produced at 24 GPa will also be less dense than surrounding mantle as defined by PREM (Fig. 5b) and there is little change in the relative densities compared to 13 GPa (Solomatova and Caracas, 2021).Only melts that have exceptionally high iron contents (e.g.Mg# ∼25) are likely to be negatively buoyant relative to the surrounding mantle as defined by PREM, either above or below the transition zone.

Implications for the deep mantle H 2 O cycle
We now investigate the implications of positively buoyant hydrous silicate melt above and below the transition zone within the context of deep mantle H 2 O circulation.When transition zone mantle upwells across the 410 km phase boundary, wadsleyite transforms to olivine, and because some amount of H 2 O is soluble in both phases, wadsleyite and olivine will coexist and exchange H 2 O over a range of temperatures (Hirschmann, 2006).For example, Fig. 6a shows a schematic binary two-phase loop in the Mg 2 SiO 4 -H 2 O system.This phase diagram illustrates that for bulk H 2 O contents less than or equal to the solubility limit in olivine, the H 2 O content of olivine when the transition is complete will be the same as in the original wadsleyite; that is, 2 O contents and olivine within the two-phase loop are dictated by H 2 O partitioning but the final olivine H 2 O content is controlled by the bulk H 2 O content.Conversely, when the bulk H 2 O content is greater than the solubility limit in olivine, a H 2 O-bearing fluid phase will form in equilibrium with H 2 O saturated olivine.
In natural mantle peridotite the situation is similar, although clinopyroxene and garnet must also saturate in H 2 O before the storage capacity of the bulk rock is exceeded and a fluid phase forms (Bolfan-Casanova et al., 2006;Férot and Bolfan-Casanova, 2012;Hirschmann, 2006;Hirschmann et al., 2005;Tenner et al., 2012).In this case, a hydrous silicate melt will form with ∼40-50 mol% H 2 O according to the phase relations depicted in Fig. 1 (My hill et al., 2017).Thus, upwelling of transition zone mantle across the phase boundary results in bulk transport of H 2 O from the transition zone to the upper mantle (e.g., here upper mantle refers to mantle above 410 km), either held completely in olivine (and coexisting phases) or in olivine plus a hydrous melt.Since our results show that melts with H 2 O and iron contents expected from melting of peridotite at ∼13 GPa will be positively buoyant, they will migrate upward.The melts will then either react with the upper mantle and hydrate it if it is undersaturated, which could eventually strip the melt of sufficient H 2 O and make it neutrally buoyant, or if the upper mantle is already saturated, the melt will continue to migrate upward, reacting with the mantle, dispersing within it as a stable melt phase, or percolating through it and ponding at the base of the lithosphere.The key point is that the melt will remain within and hydrate the upper mantle.
The low H 2 O storage capacity of lower mantle minerals relative to ringwoodite means that downwelling from the transition zone into the lower mantle can also lead to the release of H 2 O and melting similar to upwelling across 410 km (Karato et al., 2020;Ohtani, 2020;Panero et al., 2020;Schmandt et al., 2014).If the bulk H 2 O content is below the lower mantle storage capacity, ringwoodite dissociation will produce bridgmanite and ferropericlase that are H 2 O-undersaturated (Fig. 6a), but if the bulk H 2 O content is above the storage capacity of bridgmanite and ferropericlase (and Ca-perovskite), H 2 O will be released to produce a hydrous silicate melt.Again, our results suggest that hydrous silicate melts in the lower mantle beneath the transition zone will be positively buoyant such that the melt will migrate immediately back into the transition zone.That is, the net flux of H 2 O into the lower mantle at the phase boundary is through solid phase saturation only.Note that melts are never produced when upper mantle downwells into the transition zone since wadsleyite can always accommodate more H 2 O than the upper mantle.This also applies for lower mantle upwelling into the transition zone.
Positively buoyant, hydrous silicate melts above and below the transition zone have the effect of wholesale transport of H 2 O to the upper mantle (solid + melt) but reduced transport of H 2 O to the lower mantle (solid only).To gain insights into the effects of phase transitions and positively buoyant melts on the deep H 2 O cycle we construct basic box models of mantle H 2 O circulation.We recognize that the H 2 O cycle is considerably more complex in detail than in our conceptual box models, and the simplified models we present here are, therefore, heuristic.We test two whole mantle circulation models: (1) internal convective circulation, and (2) subduction driven circulation, as illustrated in Fig. 6b.

Internal circulation
The objective of this model is to develop an understanding of how H 2 O becomes distributed with time when three reservoirs of different size and H 2 O capacity are mixed through mantle-scale solid-state convection.We assume a mass exchange rate between reservoirs equivalent to the mass of mantle melted to produce the global mid-ocean ridge basalt flux (∼6.8 x 10 14 kg/yr; Bird, 2003) and assuming 10% melting on average.When a mass increment is transferred from one reservoir to another the H 2 O is assumed to For high H 2 O contents above the H 2 O storage capacity (saturation) of olivine and its polymorphs a coexisting H 2 O-rich fluid will form.Wadsleyite and ringwoodite can accommodate considerably more H 2 O than olivine and so may accommodate the excess fluid up to bulk H 2 O contents at their storage capacity (e.g., the blue dashed line).At higher pressures ringwoodite breaks down to form bridgmanite (brg) plus ferropericlase (fp) with a low relative storage capacity.High H 2 O mantle (the blue dashed line) upwelling through the wadsleyite-olivine phase transition or downwelling through the ringwoodite to bridgmanite transition will release a free fluid phase that will cause partial melting.(b) Schematic diagrams illustrating the two simplified whole mantle circulation box models considered here.Internal circulation assumes closed system recycling with a constant mass of material upwelling and downwelling across the phase transition boundaries.Subduction driven circulation assumes all downward flow in the mantle is through subduction of lithospheric plates transporting H 2 O from the hydrosphere into the mantle with return flow resulting in upwelling across phase transitions.Deep slab dehydration occurs in the transition zone and lower mantle and outgassing occurs from the upper mantle to the hydrosphere.Model parameters are provided in the main text.
instantaneously and homogeneously redistribute.Transfer of H 2 O between the transition zone and the upper mantle is modelled in two ways: (1) the total H 2 O in the mass increment is transferred in the form of undersaturated mantle or saturated mantle + excess H 2 O in a positively buoyant melt; (2) for comparison, H 2 O transfer is limited to the amount held in the saturated solid phases in the case of a negatively buoyant melt that cannot rise into the upper mantle, which is the basis for the transition zone filter model (Bercovici and Karato, 2003).Transfer of H 2 O between the upper mantle and transition zone through downwelling is limited to the storage capacity of the solid upper mantle, which we conservatively estimate at 0.1 wt% H 2 O (Férot and Bolfan-Casanova, 2012;Tenner et al., 2012).Transfer of H 2 O to the lower mantle from the transition zone is limited to the storage capacity of the lower mantle, also estimated at 0.1 wt% (Fu et al., 2019), and any excess H 2 O is immediately added back to the transition zone in the positively buoyant melt phase.Transfer of H 2 O from the solid lower mantle to the transition zone upon upwelling is limited to its storage capacity.
Fig. 7a shows the results of the internal circulation model for a case where the initial H 2 O content of the upper mantle and lower mantle are set to zero and the transition zone initially is at its storage capacity of ∼1 wt% H 2 O (∼3 ocean masses).In this model the upper mantle H 2 O concentration rises to its H 2 O storage capacity in ∼200 Ma and to its current H 2 O concentration of ∼0.02-0.04wt% H 2 O (Hirschmann, 2006(Hirschmann, , 2018) ) in ∼50 Ma.The lower mantle H 2 O content lags behind the upper mantle reservoirs due to its larger mass.For comparative purposes, Fig. 7b shows the same model with dense melts at 410 km that are reabsorbed into the transition zone, making H 2 O transport to the upper mantle much less efficient.For example, starting from a dry upper man-tle it would take ∼250 Ma to reach its current H 2 O content and even after 4000 Ma it has still not reached its H 2 O storage capacity.However, in either case, the expectation is that the upper mantle will converge to the transition zone H 2 O content.That is, dense melts at the 410 km discontinuity will eventually lead to the upper mantle reaching its H 2 O capacity because olivine will saturate with H 2 O as material upwells through the phase transition (Fig. 6a).This, however, is contrary to the observation that olivine in the MORB source is undersaturated in H 2 O (Hirschmann, 2006;Hirschmann et al., 2005;Tenner et al., 2012).The zeroth order conclusion from the internal circulation models is that transfer of H 2 O from the transition zone to the upper mantle due to the positively buoyant melt leads to a rapid rise in the H 2 O content of the upper mantle, approaching the H 2 O content of the transition zone.

Subduction-driven mantle circulation
Here we model mixing in the mantle as a combination of subduction of lithospheric plates into the lower mantle accompanied by return flow to the upper mantle.Upward transport of H 2 O across phase transition boundaries separating the three reservoirs is the same as in the model described above.We use the results of the H 2 O subduction model of van Keken et al. (2011) to model downward transport of mass and H 2 O.In their model ∼1 ocean mass (OM) of H 2 O is subducted beyond the volcanic front (> ∼150 km) into the deeper mantle in 4 Ga, a result that is also consistent with a model where the continental freeboard has been roughly constant since the Early Proterozoic (Korenaga et al., 2017).In our model no net H 2 O is added to the upper mantle by subduction; that is, we assume all H 2 O from shallow slab dehydration is returned to the hydrosphere through magmatism and no further slab  (Hirschmann, 2006(Hirschmann, , 2018)).(b) Internal circulation model as in (a) but with negatively buoyant ("dense") melts above the transition zone (e.g., the transition zone H 2 O filter model).(c) Subduction driven circulation with "light" melts above and below the transition zone and with initial H 2 O contents at estimated H 2 O capacities in the three reservoirs (transition zone = 1 wt% H 2 O and the upper mantle and lower mantle = 0.1 wt% H 2 O).(d) Subduction driven circulation with "light" melts above and below the transition zone and with a "dry" initial H 2 O content mimicking a dehydrated post Hadean mantle.dehydration occurs deeper in the upper mantle.For each mass increment, we add 70% of the slab H 2 O content to the transition zone and 30% to the lower mantle mimicking deeper slab dehydration (Iwamori, 2004;Komabayashi and Omori, 2006;Shirey et al., 2021); model outcomes for the upper mantle and transition zone are relatively insensitive to this choice, whereas the lower mantle water content depends critically on this assumption.We also include H 2 O loss from the upper mantle in this model, connecting the hydrosphere to the deep-H 2 O cycle through subduction ingassing and magmatic outgassing.We outgas the upper mantle at a rate calculated assuming an average current oceanic crust mass production of ∼6.8 x 10 13 kg/yr and assuming H 2 O is completely incompatible during melting of hydrated mantle at an average of 10% melt production.
Fig. 7c shows a case in which all three mantle reservoirs be-  (Korenaga et al., 2017).Fig. 7d shows a case where the initial mantle H 2 O content is set to zero, mimicking efficient degassing of an Hadean magma ocean and low solubility of H 2 O in crystallizing phases at very high temperatures as the magma ocean solidifies (Dong et al., 2021).In this model the H 2 O content of all reservoirs increases rapidly once subduction begins.The upper mantle lags behind the transition zone as it becomes hydrated through slab dehydration, reaching its current H 2 O content in ∼1000 Ma, with the two reservoirs converging to a level within uncertainty of the current upper mantle H 2 O content.This model predicts a near steady-state total mantle H 2 O content of about 0.6 OM after ∼ 2000 Ma, and loss from the hydrosphere of the same amount over the age of subduction (4000 Ma).This subduction driven model, beginning with an effectively dry mantle and simplified as it is, is consistent with the current upper mantle H 2 O content and is within current best estimates for total mantle H 2 O content of 0.75 ±0.2 OM based on geochemical mass balance arguments (Hirschmann, 2018), and indicates a net loss of ocean mass through time that can be reconciled with models of continental freeboard (Korenaga et al., 2017).We also note that models for the P-wave, S-wave and density jump across the 410 km discontinuity based on the elastic properties of olivine and wadsleyite produce best fits for low H 2 O contents (e.g.<0.1 wt%) in both reservoirs (Wang et al., 2019).

Summary of mantle circulation models
Our simplified box models of mantle H 2 O circulation show that the upper mantle becomes efficiently hydrated at the expense of the transition zone and that the transition zone cannot maintain a significantly different H 2 O content than the upper mantle over geological timescales of mantle recycling.Accordingly, we suggest that the H 2 O content of the upper mantle (e.g., MORB source) may be similar to the transition zone and likely higher than the bulk of the lower mantle, such that all bulk reservoirs are below their H 2 O storage capacities.Our melt density results are generally inconsistent with the formation of a global melt layer at 410 km.However, we recognize that H 2 O transport into the transition zone and lower mantle through subduction and return to the upper mantle may in reality be spatially heterogeneous and temporally episodic, which could result in localized regions of higher H 2 O content, even H 2 O saturation and melting above 410 km and below 670 km.This may explain local regions where there is seismic evidence of low shear wave velocities above and below the 410 km discontinuity.Our results suggest that if these are regions of hydrous partial melt, they are "snapshots" of melts rising from the lower mantle to the transition zone, and from the transition zone to the upper mantle.

Conclusions
Using ab initio MD simulations, we have parameterised density and viscosity of melts in the system MgO-SiO 2 -H 2 O as a function of pressure, temperature, and composition.All simulated hydrous melts are positively buoyant, inviscid and highly mobile at depths above, throughout and below the transition zone, although this behaviour may be modified with the addition of further components such as Fe.We find that ∼75% substitution of Mg by Fe is required for neutral buoyancy of model hydrous melt compositions at the upper and lower transition zone boundaries.Thus, there is no expectation that hydrous partial melts will pond at 410 km and operate as a filter to incompatible elements.Our box models of mantle H 2 O circulation indicate efficient hydration of the upper mantle at the expense of the transition zone and lower mantle and that all reservoirs tend to converge on geological timescales.We find that subduction-driven whole mantle circulation starting with a dry post-Hadean mantle is generally consistent with the current upper mantle H 2 O content (e.g., the MORB source) and produces a whole mantle H 2 O content consistent with current geochemical constraints.We suggest that the current average H 2 O content of the upper mantle as sampled by MORB may be indicative of the transition zone and perhaps higher than the bulk of the lower mantle.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig.1.Liquidus phase relations in the system MgO-SiO 2 -H 2 O (MSH) at 13 GPa with phase boundaries afterMyhill et al. (2017).Red circles show the compositions investigated in this study.The red star shows the composition of partial melt at 1550 • C in equilibrium with model mantle peridotite (green star) based on the mantle geotherm at 410 km ofKatsura et al. (2010).Blue diamonds show the compositions used byJing and Karato (2012) and purple diamonds are fromMatsukage et al. (2005) in the CMASF-H 2 O system projected onto MSH (with MgO calculated as total MgO+FeO).Green diamonds are the experimental bulk compositions ofSakamaki et al. (2006) and the yellow diamond shows the extrapolated 0.7% melt fraction model melt ofFreitas et al. (2017).Blue, cyan and green squares show the melt compositions used in previous DFT calculations(Karki et al., 2021;Mookherjee et al., 2008;Yuan et al., 2020).The green field shows hydrous silicate melts at 13.5 GPa reported byMelekhova et al. (2007), and the blue field shows melt compositions from hydrous peridotite at 11-24 GPa and 1100-1400 • C(Kawamoto, 2004;Kawamoto and Holloway, 1997) with H 2 O contents estimated from deficiency in microprobe totals.(For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Density of all simulated melts at (a) 3000 K and (b) 1800 K.The solid curves represent a global fit to all of the data (see text for details), where the grey envelopes represent 95% confidence intervals.Mantle density estimated from the Preliminary Reference Earth Model (PREM) of Dziewonski and Anderson (1981) is also shown for comparison (dashed curves).The mantle transition zone is indicated by the cyan shaded area.

Fig. 3 .
Fig. 3. Viscosity of simulated melts along the (a) Fo (Brucite also shown), (b) MSH, and (c) En + H 2 O joins at 1800 K for compositions as indicated in the legends.The solid black curves represent a global fit to the data (see text for details), where the grey envelopes represent 95% confidence intervals of the regression data.The transition zone is indicated by the cyan shaded area.

Fig. 4 .
Fig. 4. (a)Comparison of calculated melt densities (blue circles) from this study with measured melt densities of hydrous silicate melts (red circles) using the sink/float method fromJing and Karato (2012),Sakamaki et al. (2006) andMatsukage et al. (2005).Labels on the X-axis refer to the bulk compositions used in the studies.(b) Comparison of calculated melt densities in the MgO-SiO 2 -H 2 O system from this study (black and red curves) with melt densities from previous DFT-MD simulations on a range of bulk compositions (black and red circles;Karki et al., 2021;Mookherjee et al., 2008;Yuan et al., 2020); bulk compositions are provided on the figure and shown in Fig.1).

Fig. 5 .
Fig. 5. Ternary contour plots showing density of (Mg 1-x Fe x )O-SiO 2 -H 2 O melts at 1800 K at (a) the upper transition zone boundary (13 GPa) and (b) the lower transition zone boundary (24 GPa) for hydrous melts with Mg#s of 100 (x = 0), 75 (x = 0.25), 50 (x = 0.5) and 25 (x = 0.75).The black contour line delineates compositions with density equal to PREM, with a conservative estimate of uncertainty on our parameterized densities shown as 4% relative based on the range observed in previous experimental and theoretical values (e.g., Fig. 4 and Fig. S1).Symbols and boundary curves are as in Fig. 1.

Fig. 6 .
Fig. 6.(a) Schematic phase diagram in the system Mg 2 SiO 4 -H 2 O as a function of pressure at mantle temperatures.For low bulk H 2 O contents (dashed cyan line), olivine, wadsleyite and ringwoodite have the same composition except within the two-phase regions where partitioning determines H 2 O contents during the phase transition.For high H 2 O contents above the H 2 O storage capacity (saturation) of olivine and its polymorphs a coexisting H 2 O-rich fluid will form.Wadsleyite and ringwoodite can accommodate considerably more H 2 O than olivine and so may accommodate the excess fluid up to bulk H 2 O contents at their storage capacity (e.g., the blue dashed line).At higher pressures ringwoodite breaks down to form bridgmanite (brg) plus ferropericlase (fp) with a low relative storage capacity.High H 2 O mantle (the blue dashed line) upwelling through the wadsleyite-olivine phase transition or downwelling through the ringwoodite to bridgmanite transition will release a free fluid phase that will cause partial melting.(b) Schematic diagrams illustrating the two simplified whole mantle circulation box models considered here.Internal circulation assumes closed system recycling with a constant mass of material upwelling and downwelling across the phase transition boundaries.Subduction driven circulation assumes all downward flow in the mantle is through subduction of lithospheric plates transporting H 2 O from the hydrosphere into the mantle with return flow resulting in upwelling across phase transitions.Deep slab dehydration occurs in the transition zone and lower mantle and outgassing occurs from the upper mantle to the hydrosphere.Model parameters are provided in the main text.

Fig. 7 .
Fig. 7. Results of simplified whole mantle H 2 O circulation models.(a) Internal circulation with positively buoyant ("light") melts above and below the transition zone.The initial conditions are a mantle with ∼3 ocean masses (OM) of H 2 O in the transition zone with a dry upper mantle and lower mantle.Estimate of the current upper mantle H 2 O content (e.g., MORB source) is shown as a cyan bar(Hirschmann, 2006(Hirschmann, , 2018)).(b) Internal circulation model as in (a) but with negatively buoyant ("dense") melts above the transition zone (e.g., the transition zone H 2 O filter model).(c) Subduction driven circulation with "light" melts above and below the transition zone and with initial H 2 O James W.E. Drewitt: Formal analysis, Investigation, Methodology, Validation, Visualization, Writing -original draft, Writingreview & editing, Data curation, Software.Michael J. Walter: Conceptualization, Formal analysis, Funding acquisition, Supervision, Validation, Visualization, Writing -review & editing.John P. Brodholt: Methodology, Resources, Validation, Writing -review & editing.Joshua M.R. Muir: Methodology, Validation.Oliver T. Lord: Data curation, Formal analysis, Software, Validation, Visualization, Writing -review & editing.

Table 2
Results of global p-V-T-X fit.

Table 3
Regression coefficients a i, j for Equation (4), where i = 0, 1, 2, 3, 4 and j = M (MgO), OM in this case, as the H 2 O is efficiently distributed to the upper mantle.Efficient transport of H 2 O to the upper mantle also leads to increased outgassing that would result in a large increase in ocean mass over time due to the limited H 2 O subduction rate, which is inconsistent with models for ocean volume evolution with time 2 O content at a high rate to well above its storage capacity (with the excess held as hydrous melt) and equalizes with the transition zone in ∼1000 Ma.Over time all three reservoirs converge and decrease as a result of the efficiency of H 2 O transport to the upper mantle and the rate of magma outgassing.As with internal circulation, this model shows that the current low upper mantle H 2 O content is inconsistent with a high overall man-tle H 2 O content, ∼2.2