Stability and Solid Solutions of Hydrous Alumino-Silicates in the Earth’s Mantle

The degree to which the Earth’s mantle stores and cycles water in excess of the storage capacity of nominally anhydrous minerals is dependent upon the stability of hydrous phases under mantle-relevant pressures, temperatures, and compositions. Two hydrous phases, phase D and phase H, are stable to the pressures and temperatures of the Earth’s lower mantle, suggesting that the Earth’s lower mantle may participate in the cycling of water. We build on our prior work of density functional theory calculations on phase H with the stability, structure, and bonding of hydrous phases D, and we predict the aluminum partitioning with H in the Al2O3-SiO2-MgO-H2O system. We address the solid solutions through a statistical sampling of site occupancy and calculation of the partition function from the grand canonical ensemble. We show that each phase has a wide solid solution series between MgSi2O6H2-Al2SiO6H2 and MgSiO4H2-2δAlOOH + SiO2, in which phase H is more aluminum rich than phase D at a given bulk composition. We predict that the addition of Al to both phases D and H stabilizes each phase to higher temperatures through additional configurational entropy. While we have shown that phase H does not exhibit symmetric hydrogen bonding at high pressure, we report here that phase D undergoes a gradual increase in the number of symmetric H-bonds beginning at ∼30 GPa, and it is only ∼50% complete at 60 GPa.


Introduction
The storage and cycling of water in the Earth's deep interior can be greatly enhanced if hydrous silicates are stable under the high pressures and temperatures of the Earth's mantle. In particular, phases D and H are likely candidates for the storage and transport of water in the transition zone and deep mantle (e.g., [1,2]). Each phase is a hydrous magnesium silicate, MgSi 2 O 6 H 2 and MgSiO 4 H 2 , respectively. Due to the similarity in ion size, a common coupled substitution of Mg 2+ + Si 4+ = 2Al 3+ predicts the stability of both aluminous solid solutions and equivalent aluminous phases, like aluminous phase D, Al 2 SiO 6 H 2 , and aluminous phase H, otherwise referred to as δ-AlOOH.
Synthesis experiments of phase D in the absence of aluminum demonstrates that it is stable at pressures greater than 16 GPa, with increasing thermal stability with increasing pressure, such that under pressures of the base of the Earth's transition zone, phase D is stable to nearly 1700 K, suggesting it could even be stable in a mantle assemblage at moderately cool temperatures [3,4]. Phase D then breaks down to phase H with excess silica above 40-48 GPa [2,5], while phase D + brucite transforms to phase H at pressures as low as 33 GPa. Similarly, Walter et al. [6] demonstrate that increasing silica in a MgO-SiO 2 -H 2 O system increases the stability of phase H relative to phase D in equilibrium with bridgmanite [6].
Adding aluminum to the system increases the thermal stability of both phases D [5] and H [7]. Solid solutions in phase H are well documented, in which [8] demonstrated a complete solid solution series between MgSiO 4 H 2 and δ-AlOOH, with more limited mixing with SiO 2 . The entropy of mixing on this solid solution series is sufficient to explain an increase in stability by as much as 800 K. Phase D, in contrast, has a potentially incomplete solid solution series between the MgSi 2 O 6 H 2 and Al 2 SiO 6 H 2 end-member compositions, in which case crystal refinements of solid solutions show Al substituting for both Mg and Si but with partial occupancy of sites, a complexity not likely in the phase H system [9].
Multicomponent systems increase the pressure stability range of phase H [7] with a wide mutual stability region of bridgmanite, phase D, and phase H assemblage between 33 and 55 GPa [6], with some work suggesting that phase H can be stable as low as 25 GPa [10,11]. The partitioning of aluminum between phases D and H when in equilibrium, however, is unclear. Nishi et al. [5] demonstrate that Al preferentially partitions into phase H over phase D at 50 GPa and 1273 K by a factor of 6. In contrast, Bindi et al. [12] show under similar conditions that the aluminum instead favors phase D over phase H by nearly the same amount.
In this work, we focus on the stability, structure, and hydrogen bonding of hydrous phase D in the MgO-Al 2 O 3 -SiO 2 -H 2 O system, following the approach developed by [8]. Together with the results in [8], we examine the phase boundary between phases D and H, with implications on the partitioning of Al between the two phases.

Mineral Structures
For the MgSi 2 O 6 H 2 phase D composition, single crystal refinement [1] places both the Mg and Si in octahedral sites of the P31m space group; the SiO 6 units form an edge-sharing sheet of octahedra in which one in three octahedra are vacant. The MgO 6 octahedra are linked by their corners to the silica layer above and below each vacant octahedral site in the silica layer. Hydrogen atoms bond to the underbonded corners of the silica octahedra in the c-direction, extending into the gaps between MgO 6 units ( Figure 1a). As with most hydrous silicates, hydrogen atoms bond in one of two potential wells between neighboring oxygen atoms. This leads to a short OH bond and a longer H..O bond. With increasing pressure, the OH bond length increases at the expense of the H..O hydrogen bond as the two potential wells begin to merge, eventually leading to a symmetric hydrogen bond at about 40 GPa [13]. Such symmeterization is associated with a shift in the material properties.
The aluminum end-member of phase D, Al 2 SiO 6 H 2 can be thought of as a systematic substitution of two Al octahedra for pairs of Mg and Si octahedra, predicting an ordered arrangement of alternating Al and Si octahedra and an interlayer made up of AlO 6 units ( Figure 1c). However, and with analogy to phase Egg [14] and δ-AlOOH-phase H solid solutions [8], the cations are likely disordered on octahedral sites. This is supported by a structural analysis of iron-bearing aluminous phase D [9]. Pure Al 2 SiO 6 H 2 , however, may have additional disordering amongst all nominally vacant octahedral sites, which may explain the further thermal stability of this phase at intermediate pressures [15].
In contrast to phase D, phase H is a superstructure of the CaCl 2 -structure, in which all octahedra are corner sharing, with an ordering of Mg and Si on the octahedra sites [8]. The fully aluminous structure, δ-AlOOH, is similar to phase H in which each octahedral site is filled with Al. The solid solutions lead to a complete disordering between the octahedral sites with a near-zero excess enthalpy of solution for pressures between 20 and 60 GPa. A consequence of the cation-site disorder is an ensemble average bonding environment in which OH..O bonds form an asymmetric, single-well potential at midmantle pressures. This is in contrast with the end-member compositions, in which the ordered arrangements of the charge distribution leads to an effective merging of two OH..O energy wells [2,8,13].

Methods
We calculate the structure and energetics of Mg-and Al-phase D end-member compositions and Mg 1−x Al 2x Si 2−x O 6 H 2 solid solutions. To do this, we calculate the energetics of at least 10 configurations in 10 GPa steps of Mg 1−x Al 2x Si 2−x O 6 H 2 from 0 to 60 GPa for 9 compositions ranging from x = 0 to x = 1, for a total of 630 structural calculations. In each configuration, Al, Si, and Mg are randomly assigned to octahedral sites, and H atoms are placed in the stable arrangement as found in the ordered, MgSi 2 O 6 H 2 composition structure. Additionally, end-member compositions also include the fully ordered structure as an additional configuration at each pressure.

Calculations
We use the Vienna Ab-initio Simulation Package (VASP) package [16] to perform density functional theory calculations with the projector augmented wave method (PAW) [17] and generalized-gradient approximation (GGA) in the Perdew-Burke-Ernzerhof [18] formulation for the exchange-correlation part, with a kinetic energy cutoff of 800 eV. The electronic density is sampled in the reciprocal space on grids of regular and high-symmetry 2 × 2 × 2 k-points to achieve better than 0.004 eV convergence in enthalpy. Each phase D composition is modeled in a 2 × 2 × 2 supercell (88 atoms), generated in the same approach as in [8].

Thermodynamics
For each pressure and composition, we calculate the partition function, Z, as: where E i is the calculated energy of the configuration i and β is (k B T) −1 . The probability of each state, P i , is therefore a function of both the energy difference relative to other states and the temperature, according to The ensemble average energy is then calculated directly from the weighted probabilities of each configuration and represents the expectation value of the energy of that composition.
We treat the configurational entropy, S con f ig , as a counting of microstates in which each octahedral site is independent; we neglect the effects of compositional variations in the vibrational component of the entropy, and work only with the "configurational" Gibbs free energy of the system, defined as where x is a compositional index representing the Al-component of phase D.
where the end-member compositions, A and B, are MgSi 2 O 6 H 2 and Al 2 SiO 6 H 2 , respectively

Phase D-Phase H Phase Boundary
In the MgO-Al 2 O 3 -SiO 2 -H 2 O system, we solve for the equilibrium coexistence of solid solution of the phases D and H by identifying the compositions for each in which chemical potentials of the two components in the coexisting phases are equal, thereby minimizing the Gibbs free energy as a function of pressure, temperature [19,20], in the system, where phD, phH, and St refer to the phase D, phase H, and thermodynamically stable SiO 2 structures, respectively. St is in the stishovite structure below 45 GPa and the CaCl 2 structure at higher pressures. The energetics for the phase H system along the MgSiO 4 H 2 -δ-AlOOH binary are drawn from [8].
We neglect the further distribution of Al and/or H into the silica phase as the mixing enthalpies are significantly greater than in phase H (e.g., [8,21]).

Phase D Structure
We confirm the structure of phase D as determined through single-crystal x-ray diffraction techniques, with an ensemble average trigonal structure with space group P31m. Within the MgSi 2 O 6 H 2 -Al 2 SiO 6 H 2 phase D binary, we find structural consistency across all relaxed structures, in which randomly generated configurations relax towards octahedral coordinations for each Mg, Si, and Al cation with H atoms forming hydrogen bonds across the layer occupied by the MgO 6 octahedra.
Expectation values for the lattice constants and the unit-cell volume are calculated at 1500 K according to Equation (3) and fully reported in Table 1, where 1500 K is chosen to be most relevant to the synthesis conditions of many of the experiments. At zero pressure, we find unit-cell volume decreasing with increasing aluminum content dominated by the decrease in c-axis (Table 1 and Figure 2) consistent with several synthesis experiments, noting significant scatter in the synthesis results. The decrease in unit-cell volume is accommodated almost entirely along the c-axis.

Phase D Disordering and Energetics
A single exchange of neighboring Mg and Si sites within the eight molecular units of MgSi 2 O 6 H 2 leads to an excess formation enthalpy about 10 kJ/mol greater than the fully ordered structure at all pressures between 0 and 60 GPa ( Figure 1). As the configurational entropy for the single defect calculation stabilizes the system by 4.7 kJ/mol at 1500 K, we proceed with the assumption that the Mg-end member composition of phase D is minimally disordered up to its melting temperature. We further assume that Mg does not enter the layer composed of edge-sharing sheets for subsequent solid solution calculations.
In contrast, the single exchange of neighboring Al and Si sites within the eight molecular units of Al 2 SiO 6 H 2 in this calculation leads to a change of energy indistinguishable from zero up to 40 GPa, increasing to 2 kJ/mol at 60 GPa. The disordered system, therefore, becomes more stable than the ordered system for all pressures at mantle-relevant temperatures, consistent with previous experiments [15]. As a confirmation, for multiple randomly generated configurations of the aluminous phase D, we find that the ensemble average enthalpy at a given pressure is indistinguishable from the one of the nominally ordered structure.
For each configuration at each composition and pressure, we calculate the probability of each state from Equation (2) and find that the probabilities of configurations across all compositions can be described by a normal distribution, in which 67% of configurations have a probability between 2% and 14% at 1500 K. The configurations with lower probabilities are associated with a greater number of AlO 6 -MgO 6 -AlO 6 links, whereas higher probability configurations minimize such clusters. We therefore consider the configurational entropy of the system to be broadly ideal for nearly random mixing of Al on Mg-sites with associated random mixing of Al on Si-sites.
We model the binary solution in which the excess enthalpy has the form of a regular solution (Figure 3), where x Al is the mole fraction aluminum in phase D, and W is the Margules interaction parameter between the end-member compositions. Between 20 and 40 GPa, there is little variation in the interaction parameter of about 23 kJ/mol (Equation (7)). In contrast, the mixing enthalpy is indistinguishable from zero in the same pressure range along the δ-AlOOH-phase H binary.

Hydrogen Bonding
In the fully ordered, MgSi 2 O 6 H 2 structure, we find consistent structural results as [13], with a lengthening of the OH distance and associated shortening of the OH..O distance along a nearly straight OH..O bond, becoming symmetric at about 40 GPa when O-O distances decrease below 2.34 Å. This is consistent with bond symmetrization for phase H and δ-AlOOH (Figure 4). The Al-phase D, however, does not display consistent bond symmetrization: at 60 GPa only about 25% of the OH..O bonds are symmetric.
In the solid solutions, bond distances as a function of pressure exhibit a pattern intermediate to the end-member compositions. In general, with increasing pressure, shortening O-O distances tend to lengthen the OH distance and shorten the OH..O distance; eventually the hydrogen sits in a single-well potential. Symmetrization occurs when the O-O distance reaches about 2.4 Å with H locate equidistant to the oxygen atoms. In contrast to the ordered MgSi 2 O 6 H 2 structure, the solid solutions display a symmetrization that is protracted across more than 20 GPa in pressure (Figure 4). At 30 GPa, 12% of bonds are within 0.06 Å of equidistance between neighboring oxygen atoms; at 40 GPa, the number of symmetric bonds increases to 40%, with 48% of the bonds symmetric at 60 GPa ( Figure 4).
These results are in contrast to those of phase H solid solutions, which displays a failure to undergo bond symmetrization. That is due to an asymmetric potential well that develops and deepens under pressure, thus leading to a near-constant OH bond distance with increasing pressure (Figure 4), about 0.2 Å shorter than the hydrogen bond.

Phase D-Phase H Phase Boundary
Under static conditions where the configurational entropy is neglected, the idealized phase D, MgSi 2 O 6 H 2 , breaks down to phase H, MgSiO 4 H 2 with excess silica at 41 GPa. This is consistent with prior computational work [2] and consistent with synthesis experiments [5] (Figure 5). In contrast, when in the presence of excess brucite, the breakdown occurs at 23.5 GPa, in agreement with other previous computational studies [2] but about 10 GPa lower than identified in synthesis experiments [5]. In the latter, phase D is stable at 1100 • C and 32 GPa but breaks down to phase H at lower temperatures, suggesting a positive Clapeyron slope. This is consistent with a decrease in entropy as hydrogen is incorporated in the denser phase H upon breakdown of brucite with significant proton disorder [24]. In contrast to the Mg-rich composition, the ordered, aluminous phase D, Al 2 SiO 6 H 2 , is not stable at any pressure relative to 2δ-AlOOH + silica under static conditions ( Figure 5).
The aluminous phase D is likely disordered at mantle relevant temperatures, suggesting a thermal stabilization of phase D relative to δ-AlOOH. Assuming a random population of Al and Si in each of the octahedral sites and associated randomization of the OH bonding site, we find that the OH bonds relax to one of four configurations. This possible disordering of the H and the OH bonds could in principle be captured by in situ vibrational observations, like Raman, or by finite-temperature lattice dynamics calculations [25]. While these are both envisageable, they fall outside the purpose of this article. Therefore, for the time being, we approximate the configurational entropy in the aluminous phase D as S con f ig = R * (ln(4) + ln(2)) in which the first term represents the number of H bonding configurations, while the second represents the number of configurations for the Al in the Si octahedral layer. If instead there is a random distribution of all Si and Al across all octahedral sites as suggested in [15], the configurational entropy increases to S con f ig = R * (ln(4) + 1/3ln(2) + ln (3)). This latter model successfully reproduces the observed, high-temperature stability of aluminous phase D ( Figure 6) at 26 GPa, whereas the lower entropy model requires temperatures 200-300 K higher to stabilize phase D relative to δAlOOH and stishovite.  (4) (8) and (ii) an ideal configurational entropy for Mg 1−x Al 2x Si 1−x O 4 H 2 phase H as in [8], In the absence of other phases (e.g., bridgmanite, superhydrous B, phase Egg, and free water), we predict that the Mg-component of phase D is stable to 41 GPa, breaking down to phase H + SiO 2 . The addition of aluminum to the system stabilizes an alumina-rich phase H in equilibrium with an alumina-poor phase D, in which the phase H/phase D Al partitioning is 17.5 at 1300 K and 20 GPa, decreasing to about 5 at 41 GPa. At 1800 K and 20 GPa, the phase H/phase D Al partitioning is 2.2, increasing slightly with increasing pressure to 3.5 at 41 GPa Figure 7.
A second phase loop develops in the phase D + phase H mutual stability field, such that at 2100 K phase D is stable for all compositions up to 30.3 GPa, with the aluminous phase D stable to 40 GPa in equilibrium with, the relatively Al-poor phase H + SiO 2 , with a phH/phD Al partitioning of 0.95 at 2100 K between 30-39 GPa. However, the high-temperature, high-pressure phase stability predicted here of aluminous phase D cannot explain the experimental observation of mixtures of phases D and H in an alumina-rich assemblage at 45-50 GPa and 1273 K [5,12]. At these temperatures the aluminous phase D is not stable relative to δ-AlOOH+SiO 2 ; at these pressures the stability of phase D occurs only above ∼2100 K, suggesting either nonequilibrium process or an unexplored structural component may be active.

Discussion
We find that the structure of phase D can be explained by an ordered arrangement of Si and Mg in octahedral layers in the P31m space group, undergoing bond symmetrization at pressures of about 40 GPa. The aluminous phase D, in contrast, exhibits disordering between all octahedral sites at moderate temperatures, consistent with previous XRD refinements [1,15]. Through comparison with the temperature stability field of synthesized aluminous phase D, further disorder of aluminum and silicon must exist across all octahedral sites as suggested by [15].
The broad, mutual stability field of phases D and H in a silica-rich system emphasizes the complexity of these hydrous phases should they persist in the mantle. While broadly consistent with many synthesis experiments in the MSH and MASH system, significant discrepancies persist. Many synthesis experiments exploring phase D report greater hydrogen concentrations than expected from stoichiometry [15] with less than ideal Mg/Si or Al/Si ratio. They imply an additional, unexplored defect mechanism that should include either a hydrogarnet substitution of Si 4+ = 4 H + or a vacancy on the Mg-site as Mg 2+ = 2 H + , two of the most common defects stable under transition zone pressures [26,27]. Such additional mechanisms may also help explain the mismatch with multiple synthesis experiments in the a-axis as a function of aluminum content (Figure 2) as well as explain the scatter in structural results.
A few test calculations reveal that the hydrogarnet substitution in the MgSi 2 O 6 H 2 end-member phase D in the presence of brucite is more stable than the hydrogarnet-free phase D system. The difference is 10 kJ/mol decrease in Gibbs free energy at 30 GPa ( Figure 5), similar to results for synthesis experiments with excess brucite [6]. This also helps to explain chemical analyses with excess H with a Si deficit [23]. Should this stabilization energy be broadly constant as a function of Al content, this will increase the pressure at which Mg-phase D breaks down but still does not predict the stability of Al-phase D at 50 GPa, while increasing the Mg/Si ratio, counter to what is observed.
If such additional defect mechanisms are active in phase D they may enhance the hydrogen carrying capacity through the Earth's transition zone, while limiting the thermal stability of this phase.
Funding: This research was funded by NSF-EAR grant number 1724693 to Wendy R. Panero. Razvan Caracas was funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 681818 IMPACT to Razvan Caracas) and by the Research Council of Norway through its Centres of Excellence funding scheme, project number 223272. Razvan Caracas acknowledges access to the GENCI supercomputers (Occigen, Ada, Jean-Zay, and Curie) through the stl2816 series of eDARI computing grants and the TGCC supercomputers (Irene-AMD) through the PRACE grant RA4947.