Structure and Density of H2O‐Rich Mg2SiO4 Melts at High Pressure From Ab Initio Simulations

Water has a strong effect on silicate melt properties, yet its dissolution mechanism in depolymerized melts, typical for mantle composition, remains poorly understood. Here we report results of first‐principles molecular dynamics simulations for hydrous Mg2SiO4 melts with 6, 16, and 27 wt% H2O at pressure and temperature conditions relevant to the upper mantle and mantle transition zone. The results show that hydrogen bonds not only to the network‐forming cation Si but also to Mg which—nevertheless—remains the most important network modifier. There is no evidence to support the hypothesis based on experimental data that water may cause an increase in melt polymerization for ultramafic magmas; the ratio of non‐bridging oxygen per Si increases with the addition of the oxygen from H2O. The partial molar volume of water is independent on concentration in our simulations which allows us to examine the density of hydrous melt systematically. The critical water content—at which melts are neutrally buoyant compared to the surrounding mantle—is ~4 wt% H2O for a pyrolite composition, much lower than the high water content (>10 wt%) observed in petrological experiments and estimated thermodynamically for low‐degree partial melts formed in the vicinity of the mantle transition zone.


Introduction
Water is thought to play an essential role in the Earth's dynamics, as it fundamentally alters the physical properties of materials in the planet's interior. Geophysical observations have revealed the existence of low-velocity and high conductivity layers above the 410-km (Alex Song et al., 2004;Revenaugh & Sipkin, 1994;Tauzin et al., 2010;Toffelmier & Tyburczy, 2007) and below the 660-km seismic discontinuities (Liu et al., , 2018Schmandt et al., 2014), presumably indicating the occurrence of hydrous silicate melts. A water-rich mantle transition zone (MTZ) is supported by the observation of ice-VII (H 2 O) (Tschauner et al., 2018) and hydrogen within the ringwoodite crystal lattice (Pearson et al., 2014), both found as inclusions in diamonds formed in this depth range. In addition, hydrous minerals such as phase egg (AlSiO 3 OH) and δ-AlOOH, stable at pressure (P) and temperature (T) conditions of the MTZ, have also been discovered (Wirth et al., 2007). These lines of evidence suggest that the MTZ may contain a substantial amount of water regionally. Alternatively, based on constraints from mineral viscosity, Fei et al. (2017) suggested that the MTZ contains 1-2 wt% water globally.
As water preferentially partitions into melts compared to solid peridotite residues (Aubaud et al., 2004;Novella et al., 2014), partial melts of hydrous material from the MTZ are expected to be water-rich (up to 16.5 wt% H 2 O), supported by observations from melting experiments (Freitas et al., 2017;Litasov & Ohtani, 2002;Zhang et al., 2017) and insight from thermodynamic modeling (Hirschmann et al., 2009;Novella et al., 2017). Whether these water-rich melts stay at depth-and therefore provide a viable explanation for the geophysical observations-is an open question, as they appear to be far too water-rich to be neutrally buoyant. However, iron also becomes enriched in partial peridotite melts (Kawamoto & Holloway, 1997;Zhang & Herzberg, 1994), which increases melt density (ρ) substantially. Experiments determining ρ suggest that melts with up to 6-7 wt% H 2 O may be gravitationally stable at the 410-km seismic discontinuity (Matsukage et al., 2005;Sakamaki et al., 2006). On the other hand, magmas contributing to the formation of large igneous provinces (LIPs)-though they contain less water (~2 wt% H 2 O) (Kuritani et al., 2019)-are still buoyant enough to propagate to the surface from near-MTZ storage Zhao et al., 2009).
Hydrogen incorporation in silicate melts has been generally understood in terms of the reaction: The water molecule breaks up the network structure by converting one bridging oxygen, Si-O-Si, to two non-bridging oxygen (Stolper, 1982). Computations on the SiO 2 -H 2 O systems over a wide range of H 2 O contents (4-23 wt% H 2 O) (Karki & Stixrude, 2010;Pöhlmann et al., 2004;Spiekermann et al., 2016) have provided valuable insight into melt structure and dynamics and confirmed that water dissolution mechanism in polymerized melts mainly involves dissociation of H 2 O following Reaction 1.
The addition of water to depolymerized melts, on the other hand, has been proposed to potentially result in an increase of polymerization (Mysen & Cody, 2005;Xue & Kanzaki, 2004). In particular, structure factors determined by in situ X-ray diffraction on water-rich silicate melts with 15-20 wt% H 2 O support this hypothesis due to the interaction between water species and Mg cations indicated by the reaction (Yamada et al., 2007(Yamada et al., , 2011: In carbonated silicate melts, Moussallam et al. (2016) found that CO 2 − 3 -ions can alter the role of Ca and Mg cations from their usual depolymerizing role in breaking up the covalent silicate network, causing a significant polymerization of the melt.
Solid evidence of the polymerization effect of H 2 O dissolution on the network structure, however, has been lacking. To quantify the degree of melt polymerization, one needs to know the abundance of bridging oxygen. This is usually approached from experiments by deconvoluting 29 Si nuclear magnetic resonance (NMR) or Raman spectra as outlined, for example, by Moussallam et al. (2016), which may result in ambiguities in band assignment and non-unique solution, however. A detailed analysis of the local atomic structure is required, which is currently beyond the reach of experimental studies based on total X-ray and neutron structure factors or NMR and Raman spectroscopy. The structure of water-rich silicate melts and hydrogen speciation further remain elusive in part because these melts are difficult to quench to glass, or it is difficult to account for H 2 O lost during T-quench (Agee, 2008;Gavrilenko et al., 2019;Hirschmann et al., 2009;Mysen & Cody, 2005;Mysen & Wheeler, 2000).
Considerable experimental efforts have been made to investigate the role of water in controlling silicate melt properties, for example, density (Matsukage et al., 2005;Mookherjee et al., 2008;Ochs & Lange, 1999;Sakamaki et al., 2006), viscosity (Dingwell et al., 1996;Fei et al., 2017;Poe et al., 2006;Richet et al., 2000), diffusivity, and electrical conductivity (Gaillard, 2004;Ni et al., 2011;Pommier et al., 2008). Experiments have been supplemented by molecular dynamics (MD) simulations based on density functional theory (DFT) that investigate the high-P dissolution and incorporation behavior of H 2 O (Bajgain et al., 2015;Karki & Stixrude, 2010;Mookherjee et al., 2008;Pöhlmann et al., 2004;Spiekermann et al., 2016) and CO 2 Moussallam et al., 2016) in silicate melts. They provide complementary information to experiments, as the DFT-MD approach gives insight into the structure of melts at an atomic scale and allows access to conditions that are not readily explored by experiments. While DFT-MD simulations on the silica-H 2 O system (Spiekermann et al., 2016) did explore water-rich compositions (23 wt% H 2 O), the vast majority of previous theoretical investigations have focused on probing the behavior of silicate-H 2 O systems containing hydrogen with water contents typically not exceeding 10 wt % (e.g., Bajgain et al., 2015;Karki & Stixrude, 2010;Mookherjee et al., 2008). Here we go beyond this H 2 O content, following the suggestions of high water concentrations in silicate melts formed near the MTZ (Freitas et al., 2017;Hirschmann et al., 2009;Novella et al., 2017), and perform DFT-MD simulations on hydrous forsterite (Mg 2 SiO 4 ) melts with 16.6 wt% and 26.7 wt% H 2 O and compare the results to those of a more typical composition considered in simulations with 5.5 wt% H 2 O. With these simulations, we explore whether H 2 O has a polymerizing effect on depolymerized silicate melts as proposed by previous experiments (Mysen & Cody, 2005;Xue & Kanzaki, 2004;Yamada et al., 2007Yamada et al., , 2011, and we constrain the influence of H 2 O on melt densities-and therefore buoyancy-in the MTZ and its vicinity.

Computational Details
We conduct first-principles electronic structure calculations based on Kohn-Sham (KS) DFT to obtain total energy, Hellmann-Feynman stresses (diagonal terms only) and forces with a plane wave approach and periodic boundary conditions using the Vienna ab initio simulation package (Kresse & Furthmüller, 1996;Kresse & Hafner, 1993) and a planewave energy cutoff of 450 eV which is sufficient to obtain converged results for both energy and pressure. We use the projector augmented wave implementation and atomic files (Kresse & Joubert, 1999) and the generalized gradient approximation to exchange and correlation (Perdew et al., 1996). Electronic KS-DFT states are calculated at the Brillouin zone center only.
Molecular dynamics simulations are performed in the canonical N-V-T ensemble, where the number of atoms (N) and the volume (V) of the cell is held fixed, and T is controlled by a Nosé-Hoover thermostat (Hoover, 1985). After convergence of the electronic cycle, atoms are advanced with a time step (Δt) of 0.25 fs due to the presence of hydrogen in the system. Run durations are varied from 10 to 20 ps. We double the weight of hydrogen to that of deuterium, following the approach of Pan and Galli (2016) and Schwegler et al. (2008), as this permits the integration of the equations of nuclear motion with larger time steps. We find that assigning deuterium mass to hydrogen affects H dynamics (section 3.5) but has no influence on other transport or structural properties in the liquid (supporting information Figure S1). However, we calculate the density of the melts using the atomic mass of hydrogen. Initial melting and equilibration are performed at 8,000 K and followed by quenching to the desired T of 1,800 and 2,200 K at constant volume. During structure analysis, the configurations until equilibration (initial 2.5 ps) were discarded to ensure no correlation with the structure at 8,000 K. A volume range of V = 1.2V X to 0.7V X (with V X = 1,033.5 Å 3 ) for the simulation cells is generated to cover P relevant for the upper mantle and the MTZ between 0 and 15 GPa (Table S1). Empirical corrections (e.g., Mookherjee et al., 2008;Stixrude & Karki, 2005) to the resulting P computed by Hellmann-Feynman stresses are not applied.
We set up H 2 O-bearing forsterite (Mg 2 SiO 4 ) liquids in cubic supercells with compositions of 22Mg 2 SiO 4 + 10H 2 O (184 atoms), 18Mg 2 SiO 4 + 28H 2 O (210 atoms) and 13Mg 2 SiO 4 + 37H 2 O (202 atoms), corresponding to water contents of 5.5, 16.6, and 26.7 wt% ("fo06H 2 O," "fo17H 2 O," and "fo27H 2 O"), respectively (Table 1). To compare the compressibility of water in silicate melt to that of pure water, additional simulations on a supercell containing 70 H 2 O molecules are performed. As Mg 2 SiO 4 forsterite-with a melting point of 2,163 K at ambient P and an initial melting slope of 314 K/GPa (Ohtani & Kumazawa, 1981)-is not stable under many of the conditions considered in this study, we cannot compute a dry reference by DFT-MD directly. To compare to dry forsterite liquid, we use the thermodynamic description of de Koker et al. (2008) for Mg 2 SiO 4 liquid and extrapolate it to the metastable low-T region.
The high H 2 O-content in our simulations raises the question on whether demixing occurs. An inspection of the equilibrated liquid structure and radial distribution functions (RDF) g(r) shows a homogeneous phase even for fo27H 2 O. This observation is supported by two experimental inferences: (i) Novella et al. (2017) suggested that Mg 2 SiO 4 + 20.4 wt% H 2 O is a single liquid phase at high P and T based on quench textures of their samples; (ii) phase relations of the peridotite-H 2 O system (22-63 wt% H 2 O) studied by in situ X-ray radiography (Mibe et al., 2007) indicate that at P > 4 GPa the miscibility gap between silicate melt and fluid is closed; that is, hydrous silicate melt and aqueous fluid in the deep upper mantle become indistinguishable from each other. However, our observation of a homogeneous phase is not sufficient to claim that fo27H 2 O is thermodynamically stable: Phase separation rarely occurs spontaneously in regular DFT-MD simulations, a problem that partly arises from the fact that an interface between two phases is not stable in small systems (Hong & Van De Walle, 2013). For a quantitative evaluation of phase stability, it is necessary to compute the Gibbs energy of solvation of liquid H 2 O in Mg 2 SiO 4 melts by using thermodynamic integration (e.g., Yuan & Steinle-Neumann, 2020), for example. However, this approach is prohibitively expensive for the problem at hand and beyond the scope of this study.
Our simulated T-range is close to the melting point of the H 2 O-bearing system (e.g., Novella et al., 2017), permitting direct comparison to experiments, and considerably lower than previous DFT-MD simulations on hydrous silicate melts (Bajgain et al., 2015;Karki & Stixrude, 2010;Mookherjee et al., 2008;Stixrude & Karki, 2005). Simulations at such low T may not readily reach equilibrium, but we find in our results that they are relevant for liquid properties. We check that the cells considered are indeed in the molten state by monitoring the mean-square displacement (MSD) of all species and g(r) between them. The MSD-t curve displays three regimes ( Figure S2): (i) at short times t < 0.1 ps, ballistic and vibrational motion of atoms around their local equilibrium dominate, (ii) followed by a range where the atoms are temporarily confined in a cage made of their nearest neighbors, (iii) in the t range of >1 ps, the atoms leave their cages and display a typical diffusional displacement, following the linear Einstein relationship. For fo06H 2 O at 1,800 K, we have looked beyond these indicators and performed additional simulations and analysis to support this notion: 1. By performing a series of simulations from 1,600 to 3,600 K at comparable P ( Figure S3a), we find that the calculated self-diffusion coefficients of all species closely follow an Arrhenian relation. This behavior suggests that the liquid-glass transition does not occur over this T-range, as diffusion is expected to deviate from an Arrhenian relationship near the glass transition temperature (Hui & Zhang, 2007;Karki et al., 2013). 2. Extending the simulation duration to 100 ps, we find no significant changes in terms of liquid structure and transport properties ( Figure S3b).

Equation of State
We describe the P-V-T results obtained for fo06H 2 O, fo17H 2 O, fo27H 2 O, and H 2 O by combining an isothermal (T 0 = 1,800 K) Birch-Murnaghan equation of state (P iso ) of second order (third order for H 2 O) with the Mie-Grüneisen model for the thermal contribution (B TH ) to P (Bajgain et al., 2015): The thermal P coefficient B TH (V) is described by where a, b, c, and d are constants and u = V/V 0 = ρ 0 /ρ, with V 0 (ρ 0 ) the zero pressure/equilibrium volume (density) at 1,800 K. Fit parameters are summarized in Table 1; results and fits are displayed Fixed for the fitting of second-order Birch-Murnaghan equation of state.
The equation of state fit for the various hydrous Mg 2 SiO 4 melts (Table 1) allows us to calculate the partial molar volume of H 2 O in the silicate melts ( Figure 1d). Differences between fo06H 2 O, fo17H 2 O, and fo27H 2 O at 1,800 K yield molar volume of water in silicate melt (V H2O ) that are within 3% of one another for the various possible combinations; that is, V H2O does not depend on H 2 O concentration in the melt, and we display an average value in Figure 1d. These V H2O differ from those calculated using the volume difference between hydrous Mg 2 SiO 4 liquids (this study) and extrapolated dry Mg 2 SiO 4 melt (de Koker et al., 2008) by~5% at ambient P and <1% at P > 10 GPa. Higher T = 2,200 K results in insignificantly larger values for V H2O . Values for V H2O decrease from~15 cm 3 /mol at ambient P to~10 cm 3 /mol at 13.4 GPa. At low P, V H2O is significantly smaller than that computed for molecular H 2 O (V H2O ), but pure H 2 O is initially more compressible than H 2 O dissolved in the silicate melt. Therefore, for P > 4 GPa, V H2O > V H2O , but both show a similar P-dependence (Figure 1d). Experimental measurements of V H2O show a large variation, and our calculated values for V H2O compare favorably with data at 3-4 GPa by Agee (2008) and at 12-14 GPa by Jing and Karato (2012), but are higher by~20% than the values by Matsukage et al. (2005) and Sakamaki et al. (2006) in hydrous mafic/ultramafic magmas. A plausible explanation for this discrepancy is that these experiments suffered from water loss at high P and T, an aspect that was examined by Agee (2008) who found that after 30 s at T = 2,023-2,133 K the experimental charge retained only 2 wt% H 2 O of the 5 wt% in the starting material. Matsukage et al. (2005) and Sakamaki et al. (2006) did not report the H 2 O content of their recovered samples as they found their samples to consist of dendritic quenched crystals, interstitial glasses, and pores, which makes any compositional analysis challenging.
Our predicted ambient-P density ρ 0 = 2.52 g/cm 3 for fo06H 2 O (5.6 wt% H 2 O) is comparable to that of a multicomponent realistic melt analog to MORB (4.9 wt% H 2 O, iron-free) simulated previously (Bajgain et al., 2015), despite the difference in composition. H 2 O leads to a ρ-decrease in the melt by Δρ = 22-27 mg/cm 3 per 1 wt % H 2 O at 1,800 K in the P-range of 0-15 GPa, in good agreement with experimental data (Matsukage et al., 2005), despite the difference in absolute values for V H2O mentioned. The calculated Δρ due to H 2 O dissolution is slightly larger than that induced by CO 2 with Δρ~15 mg/cm 3 per 1 wt% CO 2 .

Hydrogen Speciation
The speciation of hydrogen in silicate melts can be characterized in terms of coordination environments of hydrogen by oxygen (HO n ); the coordination environments of oxygen by hydrogen (OH n ) provides complementary information. The O-H RDF is governed by a narrow maximum at a distance of~1.0 Å, followed by a first minimum. The position of this minimum (r min ) is frequently used for local atomic structure analysis, in this case to obtain HO n and OH n . With increasing P, r min , however, becomes poorly defined. The abundance of HO n in the hydrous silicate systems (and also OH n in the pure water system) strongly depends on the cutoff value which determines the distance between O and H atoms for the first coordination sphere of atoms that are directly bonded to the reference atom ( Figure S4). A small proportion of the hydrogen atoms are not bonded to oxygen when r min values are used for the analysis of H-speciation which is physically and chemically unreasonable. Starting from the position of the first sharp peak, we slowly increase the cutoff distance until all H atoms are bound to the O atoms; that is, the proportion of HO 0 reduces to zero, determining a critical value (r c ). The r c values show little or no dependence on P ( Figure S5). Average r c values determined for the pure H 2 O system at 1,800 K (r c = 1.41 Å) and 2,200 K (r c = 1.48 Å) are therefore used to examine hydrogen-oxygen coordination (HO n and OH n ).
An inspection of the partial RDFs reveals that H atoms exclusively form bonds with the only anion in the melt, that is, the O atom. We observe a peak in the RDF for H-H at~0.7 Å in a few simulations ( Figure S6 and Table S1) for the fo17H 2 O composition, which suggests the presence of metastable H 2 molecules. The simulation results show that hydrogen in the silicate melts typically occurs as OH 1 hydroxyl groups, OH 2 water molecule, and bridging HO 2 (-O-H-O-). The abundances of HO n speciation in the hydrous silicates with 5.5, 16.6, and 26.7 wt% H 2 O are very similar for a given P and T (Figure 2). The speciation of hydrogen is sensitive to variations in P and shows a weak dependence on T. Between 0 and 15 GPa, the proportion of hydroxyl groups (HO 1 ) in the three hydrous silicate systems gradually decreases from 0.95 to 0.85 at 1,800 K (Figures 2a, 2e, and 2i) and from 0.90 to 0.80 at 2,200 K (Figures 2b, 2f, and 2j). Concurrently, the proportion of bridging -O-H-O-(HO 2 ) gradually increases with P, reaching 0.15 at 1,800 K and 0.20 at 2,200 K. Our results are consistent with observations from prior DFT-MD simulations on MgSiO 3 -H 2 O systems that show a changing role of hydrogen from network modifier to network former (Karki & Stixrude, 2010;Mookherjee et al., 2008).
In the water-poor fo06H 2 O liquid, the abundance of molecular water OH 2 is negligible over the whole P-range considered (Figure 2), and more than 80% of O atoms do not bond to H at all (OH 0 ). With increasing bulk water content and increasing O/cation (especially O/Si) ratio, both OH 1 and OH 2 complexes increase in abundance. The total proportion of H-bonded oxygen (OH 1 and OH 2 ) becomes comparable to that of unbound oxygen (OH 0 ) as the water content increases to 17 wt%. At ambient P, for fo27H 2 O, OH 1 is most abundant (>0.4) and unbound oxygen (~0.3) still outnumber OH 2 molecules (~0.2). With compression, the proportion of OH 1 groups stays constant, and OH 2 groups become slightly more abundant at the expenses of the OH 0 group-their abundances are comparable at the highest P explored in this study.

Oxygen Coordination Environment and Polymerization
Oxygen is the most abundant element in silicate melts, where it occurs as "bridging" (BO), "non-bridging" (NBO), and as "free oxygen." BOs link two SiO n polyhedra (n = 4 and 5 for the P-range considered in this study), providing strong bonds between the smallest structural units of the silicate network (left-hand side of Reaction 1 and right-hand side of Reaction 2). The degree of polymerization of a silicate composition is often parameterized by (i) the number of NBO per Si cation, NBO/Si, that can vary from 0 (fully polymerized) to 4 (fully depolymerized) (Mysen & Richet, 2005) and (ii) the mole fraction of NBO (x NBO , the number of NBO per oxygen) (e.g., Lee et al., 2020); NBO/Si and x NBO are critical in determining the thermodynamic and dynamic properties of melts (Mysen & Richet, 2005;Ni et al., 2015). In our analysis, we count an O atom with two or more Si atoms as nearest neighbors as BO, an O atom with one Si atom as the nearest neighbor as 10.1029/2020JB020365

Journal of Geophysical Research: Solid Earth
NBO (i.e., OSi 1 ), and an O atom with no attachment to Si atoms as free oxygen (OSi 0 ) (which includes groups such as molecular water and Mg-O-H) ( Figure S7).
In highly polymerized melts such as SiO 2 , water dissolves by reacting with bridging oxygen to form depolymerizing Si-O-H complexes (Reaction 1), whereas for depolymerized melts, there are competing water solution mechanisms that may include both network-forming and network-modifying cations. NMR and Raman measurements on depolymerized Ca,Mg-silicate melts (Xue & Kanzaki, 2004) and peralkaline melts in the Na 2 O-SiO 2 -H 2 O system (Mysen & Cody, 2005) have been interpreted in terms of Ca-O-H and Na-O-H complexes playing a critical role in polymerizing the melt. Similarly, in the Mg 2 SiO 4 -H 2 O systems, the Mg cations may be scavenged from their network-modifying role in anhydrous melts to form Mg-O-H complexes (Reaction 2). If Reaction 2 is P-dependent, the equilibrium proceeds toward one side upon compression. When the reaction proceeds to the right, the relative proportion of the Si-O-Si complexes (i.e., the degree of melt polymerization) increases which potentially provides an explanation for the observed anomalous P-evolution in the position of the first sharp diffraction peak (FSDP) in the structure factors determined by in situ X-ray diffraction (Yamada et al., 2007(Yamada et al., , 2011, a subject we further discuss in section 3.4.  Hydrogen-oxygen speciation and abundance of species for different water contents in our simulations at 1,800 K (a, e, i, c, g, k) and 2,200 K (b, f, j, d, h, l). HO 1 and HO 2 indicate that a hydrogen atom is bonded to one and two oxygen atoms, respectively; OH 0 , OH 1 , and OH 2 indicate that an oxygen atom is bonded to zero, one, and two hydrogen atoms, respectively.  (Figures 3a-3d), and the Si-O-Mg complex remains important in fo27H 2 O despite a high water concentration (Figures 3e and 3f). In line with the experimental observation that water species can be linked either to Si or metal cations in depolymerized melts (Mysen & Cody, 2005;Xue & Kanzaki, 2004), our results show that a large proportion of Mg cations are bonded to OH. While in the experimental studies the formation of free hydroxyls is interpreted as an increase in polymerization degree (e.g., Xue & Kanzaki, 2004)

X-Ray Structure Factor of Hydrous Mg 2 SiO 4 Melts
Our DFT-MD simulation results can be compared directly to experiments by calculating the scattering structure factor. We compute partial structure factors S αβ (q) between species α and β from the Fourier transform of the corresponding partial pair distribution function by means of where c α(β) = N α(β) /N is the concentration of α(β) species with α,β ∈ {Si, Mg, O, H} (Gutiérrez & Johansson, 2002). The total X-ray scattering structure factor can be obtained from the partial structure factors by weighting them with X-ray form factors f α (s) (Winkler et al., 2004): For fo17H 2 O, we compare the results of the dominant partial S(q) and S X-ray (q) at 6.4 GPa and 1,800 K (Figure 4) with the experimental data for dry glassy Mg 2 SiO 4 (Kohara, 2004) and hydrous liquid Mg 2 SiO 4 (Yamada et al., 2011). The structure factor and the FSDP, in particular, provide information on the intermediate-range order with a spacing larger than the interatomic distance (Elliott, 1991). The position of the FSDP at zero P is predicted at q~2.18 Å −1 for fo06H 2 O, comparable to that in vitreous Mg 2 SiO 4 (Kohara, 2004) and shifts to q~2.22 Å −1 for fo17H 2 O and to q~2.32 Å −1 for fo27H 2 O, representing a general trend of a larger q for the FSDP with increasing H 2 O content at constant P. This is consistent with experimental observations regarding the effect of hydrogen (Zotov et al., 1996(Zotov et al., , 1992 and the alkali elements K and Na (Du & Corrales, 2006) on the FSDP position in melts and glasses. We find a second peak close to the FSDP at~3.0 Å −1 , consistent with the one observed in experiments (Yamada et al., 2011). This second peak is correlated with S MgO .
The FSDP of glassy forsterite (Mg 2 SiO 4 ) determined by X-ray diffraction is located at q~2.2-2.3 Å −1 (Kohara, 2004;Kohara et al., 2011), larger than for MgSiO 3 enstatite glass with q~2.1-2.2 Å −1 (Funamori et al., 2004;Kohara et al., 2011). In hydrous Mg 2 SiO 4 melts, H interrupts the covalent network connectivity in a similar way as Mg and should therefore lead to a decrease of medium-range order, resulting in a shift of FSDP to larger q values. By contrast, the measured FSDP position (q~2.0 Å −1 ) in hydrous Mg 2 SiO 4 melt (Yamada et al., 2011) is at lower q than that in anhydrous Mg 2 SiO 4 glass (Kohara, 2004;Kohara et al., 2011), and similarly, the FSDP position in hydrous MgSiO 3 melt (Yamada et al., 2011) is at lower q than that in dry MgSiO 3 melt (Funamori et al., 2004;Kohara et al., 2011) at P~6 GPa. These results seem to be at odds with the previous experimental conclusion for hydrous glass (Zotov et al., 1996(Zotov et al., , 1992. For a quantitative comparison of the FSDP positions from different experiments, one needs to take T-effect into account, however. Our DFT-MD results suggest that the FSDP position in hydrous Mg 2 SiO 4 is positively correlated to H 2 O content and P ( Figure S10) and there is no evidence to support an increasing degree of melt polymerization with P to 15 GPa and H 2 O content to 27 wt%. We emphasize that the variation of water content and Mg/Si ratio of the melts may lead to different conclusions. Melts studied by Xue and Kanzaki (2004) contain less MgO compared to the Mg 2 SiO 4 melt considered in this study; although the Mg 2 SiO 4 melts studied previously (Yamada et al., 2007(Yamada et al., , 2011 are water-rich, the diamond capsule with platinum caps used in their experiments may not fully retain water. The anomaly in the FSDP position with respect to P observed by Yamada et al. (2011), if true, more likely occurs at water-poor conditions. This highlights the need to conduct further simulations on water-poor and less depolymerized melts such as MgSiO 3 .

Diffusivities
The diffusion constants of a species α in the liquid can be calculated from the MSD using the Einstein rela- . As deuterium mass was assigned to hydrogen atoms, the actual calculations of the dynamic properties are performed for deuterium instead of hydrogen. Hydrogen diffusivity (D H ) is larger than deuterium diffusivity (D D ) with a D H /D D ratio of~1.2 for the one example (fo27H 2 O at 0 GPa and 1,800 K) we have explored explicitly ( Figure S1). To determine the uncertainty on diffusivity, we performed multiple simulations (n = 50 runs) starting from the same configuration but different initialized velocities both at 1,800 and 2,200 K, following the suggestions of Maginn et al. (2019). For each trajectory, we performed a linear fitting of MSD-t curve in a time range between 1 and 10 ps and obtained a distribution of diffusivity (expressed as histograms in Figure S11). The standard deviation of this distribution for diffusivity is in the range of 30-40% at 1,800 K and 20-30% at 2,200 K. Simulations at high T have faster dynamics and better sampling whereas those at lower T sample fewer linear diffusion events and therefore are associated with larger statistical uncertainty.
As expected, the computed diffusivity sequence is D H/D > D Mg ≥ D O > D Si for all conditions considered ( Figure 5). Hydrogen diffuses fastest in the melts; Si is the slowest species as it is confined within SiO n polyhedra; network modifiers (i.e., Mg) are weakly bonded to oxygen and therefore relatively mobile. Pressure leads to decreased mobility of species and consequently to smaller D α -values for Si, Mg, and O, while H/D does not experience a significant reduction in D α over the entire P-range considered. Compared to dry Mg 2 SiO 4 liquid in previous molecular dynamics simulations (Adjaoud et al., 2011;Ghosh & Karki, 2011), the diffusion coefficients of all species (Si, Mg, and O) are larger over the entire P-range (Figures 5b-5e). Self-diffusion coefficients increase approximately linearly with water content (C H2O in wt%) both at ambient and high P (Figure 5a).
Mass transport properties of liquids, such as diffusivity, are well known to be strongly related to melt structure and can therefore be used as an indicator for changes in it (e.g., Posner et al., 2017). Assuming that there is a structural change in response to P as proposed by Yamada et al. (2007Yamada et al. ( , 2011 in hydrous Mg 2 SiO 4 melts

10.1029/2020JB020365
Journal of Geophysical Research: Solid Earth (Reaction 2), an anomaly in diffusivity would be expected. However, our predicted self-diffusion coefficients display a typical P-induced decrease for depolymerized melts Mookherjee et al., 2008;Ni et al., 2015), supporting the interpretation that H 2 O dissolution does not cause an increase in the degree of polymerization of the network structure for mafic-ultramafic melts (here Mg 2 SiO 4 ).

Implications for the Gravitational Stability of Hydrous Melts Near the MTZ
The MTZ is potentially a significant water reservoir in the deep Earth due to the high solubility of water in its dominant minerals, wadsleyite and ringwoodite (Inoue et al., 1995;Kohlstedt et al., 1996;Kudoh et al., 1996). When material moves away from the MTZ, dehydration melting may occur as the amount of water held in MTZ minerals probably exceeds the saturation limit of the phase assemblages in the upper and lower mantle (Bercovici & Karato, 2003;Schmandt et al., 2014). Moreover, melts are unlikely to be generated in the dryperidotite system at comparable conditions, as the solidus T increases significantly with P (Takahashi, 1986;Zhang & Herzberg, 1994). The origin of low-velocity zones near the MTZ has therefore been widely attributed to hydrous partial melting.
The gravitational stability of melt in the deep mantle is primarily controlled by the ρ-contrast between silicate melt and surrounding peridotite rock, which depend primarily on T and composition. Temperatures at which partial melting occurs at the 410-km and 660-km seismic discontinuities are estimated to be in the range of~1,800-1,900 K (Chust et al., 2017;Ito & Katsura, 1989), which is lower than T used in most previous determinations of melt ρ by experiments (e.g., Matsukage et al., 2005;Sakamaki et al., 2006) and especially simulations (e.g., Bajgain et al., 2015;Mookherjee et al., 2008).
The generation of melt in ultramafic composition with a ratio of Mg/Si = 2 used in our simulations may appear unlikely in the mantle. However, melts produced by incongruent melting of peridotite move to ultramafic compositions with increasing P, and the ratio (Mg + Fe)/Si is approaching a value of 2 or larger for hydrous compositions (Inoue, 1994;Kawamoto, 2004;Kawamoto & Holloway, 1997) at MTZ conditions, making our results relevant for considerations of partial melting there. Using the compressibility behavior of dissolved water in silicate melt ( Figure 1 and Table 1) and the compressibility of the dry silicate liquid, ρ of hydrous silicate melts can be estimated based on our results. We examine ρ of melts with idealized compositions of Mg 2 SiO 4 and (Mg 1−x Fe x ) 2 SiO 4 , and a multicomponent realistic peridotitic melt at 1,800 K, and compare them with ρ-profiles of seismic models in and around the MTZ, especially just above the 410-km seismic discontinuity.
We determine the density of Mg 2 SiO 4 -based melts by combining the compressibility parameters of dry Mg 2 SiO 4 melt (de Koker et al., 2008) with those of the water component in silicate melt from our DFT-MD results (Tables 1 and S3). Calculated ρ of both dry and hydrous Mg 2 SiO 4 melts are significantly lower than ρ-values of seismic models (Figure 6a). For instance, ρ of dry Mg 2 SiO 4 melt is smaller bỹ 10-12% compared to the Preliminary Reference Earth Model (PREM) (Dziewonski & Anderson, 1981) and the average value of the 99 pyrolite models from Monte Carlo inversion for mantle structure (99PREFum) (Cammarano et al., 2005) throughout the MTZ; above the 410-km seismic discontinuity, they are lower by~6% and~10% compared to these profiles, respectively.
A ρ-crossover between melts and crystals is only possible by taking the effect of the iron component into account. Under anhydrous conditions (Mibe et al., 2006), a distribution coefficient K  (Kawamoto & Holloway, 1997). Using K FeO=MgO D = 0.3 and assuming the melt is in equilibrium with a pyrolite mantle which has an Mg# of 89 (i.e., the molar ratio of Mg/(Mg + Fe) × 100), the hydrous melt has an Mg# of 71. Using a Vegard-type mixing model in which the molar V of a silicate melt is given by a linear combination of partial molar V of MgO and FeO oxide components using the ratio from Lange and Carmichael (1987) and appropriate molar weights, we estimate ρ 0 for (Mg 0.71 Fe 0.29 ) 2 SiO 4 melt at 1,800 K (Table S3). While ideal mixing of volumes has proven to be a reliable method of determining ρ 0 of silicate melts (Lange & Carmichael, 1987), there is evidence that the partial molar volume FeO (V FeO ) varies with composition (e.g., Guo (Jing & Karato, 2008), we use the linear mixing approach. We further assume that values of K T and K′ of (Mg 0.71 Fe 0.29 ) 2 SiO 4 melt are the same as those used for Mg 2 SiO 4 , insensitive to change of Mg# (Jing & Karato, 2008). For 5 wt% of H 2 O, this melt follows the seismic ρ-profile throughout the MTZ to first order (Figure 6b) but is positively buoyant above both the 410-km and 660km seismic discontinuities. The critical water content-at which the melt is neutrally buoyant relative to mantle peridotite-at 13.4 GPa and 23.8 GPa is~9 wt%.
Finally, we evaluate ρ of hydrous melt with a representative mantle composition by using ρ H2O from our DFT-MD simulations (Tables 1 and S3) in conjunction with ρ of (dry) peridotitic melt. For that, we use an experimental composition in the CaO-MgO-Al 2 O 3 -SiO 2 pyrolite system with 2 wt% H 2 O (Litasov & Ohtani, 2002)-a water content consistent with mineral viscosity considerations by Fei et al. (2017)-that results in a partial melt composition of 49.2 wt% SiO 2 , 3.4 wt% Al 2 O 3 , 34.5 wt% MgO, and 12.9 wt% CaO. Other oxide components such as K 2 O and Na 2 O are not considered because of their low concentrations and therefore negligible effects on melt ρ (Freitas et al., 2017). To include FeO in the melt chemistry, an Mg# of 71 was used following the arguments given above, and we estimate ρ 0 based on partial molar volume ratios by Lange and Carmichael (1987). We ignore non-ideal mixing effects between SiO 2 , Al 2 O 3 , FeO, MgO, and CaO expressed by an excess volume term (e.g., Xin et al., 2019), corresponding to the possible interactions between them, as it results in only a marginal shift of ρ 0 . For compressibility of the dry peridotitic melt, we use data by Ohtani and Maeda (2001) and Sakamaki et al. (2009) (Table S3). This melt model results in ρ-profiles at MTZ conditions shown in Figure 6c: The dry peridotite closely follows the seismological ρ-profile throughout the MTZ, and the critical water content is~4 wt% H 2 O just above the 410-km seismic discontinuity. This water content is lower than estimated previously, for example, 6-7 wt% H 2 O in Matsukage et al. (2005), who have adopted V H2O < 8 cm 3 /mol, at the low end of their experimental constraints and smaller than the values found in our study ( Figure 1d).
Consequently, the key consideration for the fate of hydrous partial melts formed at the edges of the MTZ is the H 2 O content. Both melting experiments (Litasov & Ohtani, 2002;Zhang et al., 2017) and thermodynamic modeling (Hirschmann et al., 2009;Novella et al., 2017) suggest that partial melts above the 410-km seismic discontinuity form with a minimum of~15 wt% dissolved H 2 O, which is much larger than the critical water content at the 410-km (~4 wt%) and 660-km seismic discontinuities (~0 wt%) estimated here. This small critical H 2 O content suggests that water is a significant buoyancy source for triggering deep plumes ascending from lowermost upper mantle storage as imaged by seismic tomography Xia et al., 2016;Zhao et al., 2009).

Conclusion
We have simulated hydrous Mg 2 SiO 4 melts with 5.5-26.7 wt% H 2 O using ab initio molecular dynamics methods to study their structure, density, and dynamic properties over P-T conditions of the upper mantle and the mantle transition zone. Despite high water contents in the melts, the Si-O-Mg complex constitutes a large proportion of the X-O-X complexes (X represents Mg, Si, and H cations), and it is the major cause for depolymerizing melt structure. Water hardly changes the depolymerization role of the Mg cations, even though a strong interaction between water species and Mg can be observed from both our simulations and previous spectroscopy studies (e.g., Xue & Kanzaki, 2004). However, our results based on combined evidence from static and dynamic properties of the melt do not support the hypothesis inferred from NMR data and X-ray structure factors that ultramafic magmas may be polymerized by the presence of dissolved water. This result is quite different from the polymerizing effect of CO 2 on the melt network shown in previous simulations (Moussallam et al., 2016).
We have examined the buoyancy stability of water-rich silicate melts generated through dehydration melting near the mantle transition zone. Taking iron partitioning into account, the critical water content at which hydrous melts have the same density as the mantle (density crossover) is less than 10 wt% despite iron enrichment in the melt. Therefore, water-rich melts cannot experience long-term gravitational stability in the mantle transition zone and its vicinity. Low velocity zones imaged near the 410-km and 660-km seismic discontinuities are likely transient features resulting from continuous downward/upward flow across the transition zone boundaries.

Data Availability Statement
Computations were performed at the Leibniz Supercomputing Center of the Bavarian Academy of Sciences and the Humanities. The Vienna Ab Initio Simulation Package (VASP) is proprietary software for which licenses are issued online (https://www.vasp.at/). In compliance with the journal's FAIR data availability requirements, the simulation files are described in computational methods and are archived in Figshare (https://doi.org/10.6084/m9.figshare.11803113).