Hydrogen Diffusion in the Lower Mantle Revealed by Machine Learning Potentials

Hydrogen may be incorporated into nominally anhydrous minerals including bridgmanite and post‐perovskite as defects, making the Earth's deep mantle a potentially significant water reservoir. The diffusion of hydrogen and its contribution to the electrical conductivity in the lower mantle are rarely explored and remain largely unconstrained. Here we calculate hydrogen diffusivity in hydrous bridgmanite and post‐perovskite, using molecular dynamics simulations driven by machine learning potentials of ab initio quality. Our findings reveal that hydrogen diffusivity significantly increases with increasing temperature and decreasing pressure, and is considerably sensitive to hydrogen incorporation mechanism. Among the four defect mechanisms examined, (Mg + 2H)Si and (Al + H)Si show similar patterns and yield the highest hydrogen diffusivity. Hydrogen diffusion is generally faster in post‐perovskite than in bridgmanite, and these two phases exhibit distinct diffusion anisotropies. Overall, hydrogen diffusion is slow on geological time scales and may result in heterogeneous water distribution in the lower mantle. Additionally, the proton conductivity of bridgmanite for (Mg + 2H)Si and (Al + H)Si defects aligns with the same order of magnitude of lower mantle conductivity, suggesting that the water distribution in the lower mantle may be inferred by examining the heterogeneity of electrical conductivity.

throughout the lower mantle.Therefore, we investigate the hydrogen diffusion in Al-bearing bridgmanite considering (Al + H) Si defects using the on-the-fly machine learning force field (MLFF) trained from ab initio data.
We consider three supercells containing 80 (2 × 2 × 1 for Pbnm bridgmanite), 360 (3 × 3 × 2 for Pbnm bridgmanite), and 2880 (6 × 6 × 4 for Pbnm bridgmanite) atoms.The water content of the system is controlled by substituting different numbers of Si atoms with Al atoms and hydrogen atoms, generating (Al + H) Si defects that are evenly distributed in the supercell.All ab initio calculations and MLFF generations are performed using VASP (Kresse and Furthmuller, 1996).For ab initio calculations, we adopt the PBEsol approximation (Perdew et al., 2008) with the projector augmented wave (PAW) method (Kresse and Joubert, 1999).The core radii of O (2s 2 2p 4 ), Si (3s 2 3p 2 ), Mg (2p 6 3s 2 ), Al (3s 2 3p 1 ), and H (1s 1 ) are 0.820 Å, 1.312 Å, 1.058 Å, 1.402 Å, and 0.370 Å, respectively.We employ a 500-eV energy cutoff and sample the Brillouin zone at the Γ point.All molecular dynamics (MD) simulations are performed in the NVT ensemble at 2000 K controlled by the Nosé-Hoover thermostat (Hoover, 1985) with a time step of 0.5 fs.During MD simulations, MLFFs are trained using total energies, forces, and stress tensors from automatically selected DFT structures.Similar to the DeePMD method in the main text, the potential energy of the system is decomposed as the sum of local energies derived from the local environments of atoms.A Bayesian error estimate is applied to decide between DFT or MLFF forces for the MD simulation.If the predicted errors exceed the threshold, a new reference structure will be chosen and recalculated, based on which the coefficients of the MLFF will be adjusted.A detailed description of the MLFF approach can be found in Jinnouchi et al. (2019a) and Jinnouchi et al. (2019b).
We initiate our on-the-fly training with a small system (Mg 16 Si 14 Al 2 O 48 H 2 ) and extend the training set by continuing the on-the-fly training using larger systems (3 × 3 × 2 for Pbnm bridgmanite) with various water contents (0.12-1.68 wt%) at different pressures (24,25,and 26 GPa).The Al content in the DFT training set ranges from 0.37-5.04wt% and is proportional to the water content.The final training set consists of 1800 DFT configurations and the corresponding MLFF is refitted for prediction-only MD runs.
We compare the energies, atomic forces, and stresses from the MLFF to those from DFT calculations for 200 configurations that are not included in the training set, and the root-mean-square errors of the energies, atomic forces, and stresses are 0.55 meV atom −1 , 0.09 eV Å−1 , and 0.09 GPa, respectively.Since the training set only contains one mineral phase with one hydrogen substitution mechanism, and the MLFF is limited to ∼2000 K and 24-26 GPa, this error is much smaller than our MLP introduced in the main text.= 0.12 wt%).The simulation time is 100-200 ps for each run and 10 independent runs with different initial velocity distributions are performed for each system to evaluate the uncertainties.Two 100-ps trajectories of hydrogen in the largest systems with the (Al + H) Si defect and the (Mg + 2H) Si defect are shown in Fig. S7a and Fig. S7b, respectively.The hydrogen diffusion coefficients in all those systems are derived from the mean square displacements and plotted in Fig. S8.For Mg 17 Si 15 O 48 H 2 with the (Mg + 2H) Si defect, we used both the MLP and MLFF methods to perform simulations under completely identical conditions, and obtained very close diffusion coefficients (with only a 3.4% difference, Fig. S8a), verifying the consistency of the results from these two methods.Consistent with the MD trajectory of hydrogen associated with (Al + H) Si reported by Muir and Brodholt (2018), hydrogen atoms are not strongly bound to a specific oxygen atom in a fixed lattice site, but are able to diffuse far from the Al 3+ ion.In addition, at 2000 K and 25 GPa, the hydrogen diffusion coefficients for (Al + H) Si and (Mg + 2H) Si both fall within the order of magnitude of 10 8 m 2 s −1 , exceeding (2H) Mg and (4H) Si by at least an order of magnitude (see Fig. 5).We conclude that (Al + H) Si and (Mg + 2H) Si exhibit remarkably similar hydrogen diffusion patterns with closely matching hydrogen diffusion coefficients at the uppermost lower mantle conditions, considering all three water and aluminum contents.
Owing to this intrinsic similarity between (Al + H) Si and (Mg + 2H) Si , we propose that the hydrogen diffusion in these two defects can be treated analogously, and our discussions about (Mg + 2H) Si might well apply to (Al + H) Si too.

Text S2 Defect distribution in the supercell
To accurately quantify the distribution of defects in the supercell, we used the average inter-defect distance of nearest neighbors (ξ ) as a collective variable defined as where N is the total number of defects in the system, and d min is the distance between a defect and its nearest neighboring defect.Taking (2H) Mg as an example, given positions of Mg atoms that will be substituted by hydrogen atoms, i.e., the positions of defects in the initial configuration for MD simulations, we compute the distance between all defects and their respective nearest neighboring defects, i.e., d min , which are then averaged to derive ξ .This parameter can quantify how far away the different water-containing defects are.
Specifically, small ξ means that defects tend to gather into clusters, and large ξ means that the defects are dispersed with the largest ξ corresponding to the uniform distribution of defects.
In order to determine the stability of different defect distributions and find the most energy-favored one, we conduct a simple Monte Carlo simulation.We generate 100 different hydrous bridgmanite (Mg 60 Si 64 O 192 H 8 ) configurations with (2H) Mg defects, among which 98 configurations are randomly generated using Atomsk (Hirel, 2015), that is, 4 of the 64 Mg atoms are randomly selected to be substituted by hydrogen atoms.The rest two configurations include one with the largest ξ where defects are evenly distributed in the supercell, and the other with the smallest ξ where 4 nearest Mg atoms within a unit cell were replaced with H atoms.
We calculated the enthalpy of each configuration using both MLP and DFT methods, as shown in Fig. S9.
For the MLP results, although the different defect distributions lead to some variation in the enthalpy of the system, the standard deviation is only 6.06 meV atom −1 enthalpies are almost identical considering the uncertainty of our MLP, which is smaller than the RMSE of our MLP (RMSE = 6.39 meV atom −1 ).This is more pronounced in the DFT results, where the standard deviation of the enthalpy is as low as 0.46 meV atom −1 , which is far lower than the theoretical uncertainty of calculating enthalpy using the DFT method (tens of meV atom −1 , Wang et al., 2021).Therefore, we argue that the enthalpy variation in Fig. S9 is negligible and all random configurations with different defect distributions may have almost the same enthalpy.Therefore, it may be true that the probability of various defect distributions occurring in nature is approximately equal, and it is difficult to find a specific defect configuration that is preferred in terms of energy.Thus, we can directly use any type of defect distribution configuration for diffusion simulation, instead of searching for the most stable one.
We investigate the impact of defect distribution on hydrogen diffusion by simulating and calculating the hydrogen diffusion coefficients using multiple different configurations with randomly distributed defects.
We then compared the results with those for which defects are uniformly distributed, in order to quantify how the defect distribution ultimately affects the hydrogen diffusion coefficient.We choose two representative bridgmanite systems at 25 GPa and 2000 K: Mg 16368 Si 16384 O 49152 H 32 (0.0175 wt% water, Fig. S10a) and Mg 8128 Si 8192 O 24576 H 128 (0.140 wt% water, Fig. S10b).The blue square represents the average value of hydrogen diffusion coefficients calculated from 10 independent simulations with different initial velocity distributions, considering the same supercell containing uniformly distributed defects.The orange square represents the average result obtained from 10 different supercells with randomly distributed defects.The results suggest that while using a random defect distribution may slightly increase the uncertainty of the calculated diffusivity, the impact on the mean is negligible compared to the calculation error, validating the result obtained from the uniform defect distribution in the manuscript.Table S2: Hydrogen diffusion coefficient and proton conductivity of bridgmanite (Brg) and post-perovskite (pPv) along an adiabatic geotherm Katsura et al. (2010).

Figure S1 :Figure S2 :Figure S3 :
Figure S1: Comparisons of mean square displacements (MSD) and hydrogen diffusion coefficients (D H ) derived with the density functional theory (DFT) and the machine learning potential (MLP) for hydrous bridgmanite (Mg 15 Si 16 O 48 H 2 ) at 3000 K and 25 GPa.The shading and the error bar represent the standard deviation of 10 independent simulations.(a) MSD as a function of the simulation time using double logarithmic axes; (b) MSD as a function of the simulation time using linear axes; (c) D H calculated using DFT and MLP.

Figure S5 :
Figure S5: The atomic positions along the three crystallographic axes as a function of time in hydrous bridgmanite for (a) the (2H) Mg defect (Mg 960 Si 1024 O 3072 H 128 ) and (b) the (Mg + 2H) Si defect (Mg 1088 Si 960 O 3072 H 128 ) at 2000 K and 25 GPa.The trajectories of these atoms are shown in Fig. 3.

Figure
Figure S6: (a) Color mapping of the trajectory of 32 hydrogen atoms over time in hydrous bridgmanite for the (2H) Mg defect (Mg 16368 Si 16384 O 49152 H 32 , C water = 0.0175 wt%) at 2000 K and 25 GPa.As the simulation time increases from 0 ps to 1000 ps, the color of the trajectory changes along the visible spectrum.All atoms are omitted for clarity.(b) The atomic position along the three crystallographic axes as a function of time of one hydrogen atom (pointed by gray arrows).

Figure S7 :
Figure S7: Color mapping of the trajectories of two hydrogen atoms over time in hydrous bridgmanite for (a) the (Al + H) Si defect (Mg 576 Si 568 Al 8 O 1728 H 8 , C water = 0.12 wt%) and (b) the (Mg + 2H) Si defect (Mg 580 Si 572 O 1728 H 8 , C water = 0.12 wt%) at 2000 K and 25 GPa.As the simulation time increases from 0 ps to 100 ps, the color of the trajectory changes along the visible spectrum.Orange, blue, red, and white spheres represent magnesium, silicon, oxygen, and hydrogen, respectively.

Figure S9 :
Figure S9: The static enthalpy of 100 different defect distributions of hydrous bridgmanite (Mg 60 Si 64 O 192 H 8 ) calculated using DFT and MLP at 25 GPa.The error bars represent the 2RMSE of our MLP.