ENVIRONMENT AND ENERGY

The solid Earth is mainly composed of rocks consisting of both minerals, i.e. crystalline phases of natural occurrence, and amorphous or glassy material. In geologically active regions, e.g., continental margins, subduction zones or volcanic environments, high temperatures and pressures and/or volatile chemical components lead to the formation of molten rock, mineral alteration processes, or the release of magmatic-hydrothermal fluids. Such processes have shaped our present-day Earth and triggered chemical differentiation that is reflected in various rock types, ranging from common basalts or granites to technologically important ore deposits. The solid, liquid and gaseous phases found in the Earth’s interior are generally chemically complex, and their properties are often not directly accessible at the relevant thermodynamic conditions. Therefore, interpretation of field observations and understanding of the underlying mechanisms of geological processes require insights from laboratory experiments and numerical simulations.

In the context of this long-term supercomputing project, our research focused on the following two goals:

• development of reliable structure models for non-crystalline phases such as oxide and silicate melts and glasses or aqueous fluids in a wide range of pressures and temperatures. This is an important computational task since the three-dimensional atomic structure of non-crystalline phase is not directly accessible experimentally. The validation of the simulation models was primarily done by comparison to experimental fingerprints of the respective samples. This included data from X-ray and neutron diffraction, vibrational spectroscopy (IR and Raman), as well as X-ray spectroscopy (XANES, EXAFS, X-ray Raman spectroscopy).
   
• prediction of thermodynamic properties of crystalline and non-crystalline phases. Knowledge of these properties is required to understand phase behavior in the Earth’s interior and the (re-)distribution of major and minor chemical elements and their isotopes during geochemical reactions. In addition, the simulations help to parameterize advanced thermodynamic models and equations of state, which are key components for larger scale geophysical, geochemical and geodynamic models.

To achieve these goals, proven computational chemistry methods were used, namely molecular dynamics simulations in combination with electronic structure or classical interaction potential methods. This numerical modeling approach makes it possible to simulate the motion of atoms and molecules at defined volume (or pressure) and temperature in real time, i.e. on the time scales of femto- to nanoseconds. Depending on the system, the particle interactions are either described by rather accurate but computationally expensive first-principles methods based on density-functional theory, or by more efficient force fields, which usually have a limited scope. In both cases, the system sizes and simulation times are limited by the available computational resources. The allocated computational time allowed us to perform extensive simulations with dedicated MPI-OpenMP-parallelized codes (e.g. CPMD, CP2K, ABINIT) on the JURECA and JUWELS CPU supercomputers at JSC.

One of our research interests is the study of speciation and thermodynamic behavior of metal ions in supercritical aqueous fluids. Such fluids are geologically important in various environments such as magmatic-hydrothermal systems, hot geothermal reservoirs, or subduction zones. We considered pure H₂O and chloride and fluoride solutions as simplified solvents for the metal cations (e.g. Sahle et al. 2013; Jahn et al. 2015; Sahle et al. 2017; Stefanski et al. 2020). Significant structural changes in solvents are observed when pressure and temperature vary from near-surface conditions to conditions in the lower crustal or upper mantle (e.g. Sahle et al. 2018; Stefanski et al. 2018; Elbers et al. 2021). Solvent properties determine speciation and solubility of metal ions. For example, the low dielectric constant of water destabilizes highly charged species at high temperatures and pressures, leading to the formation of (near) neutral contact ion pairs.

Fig. 1 shows, as an example, structural information of Y speciation in mixed fluoride and chloride supercritical solutions from ab initio molecular dynamics simulations using the CPMD code. In order to obtain the stability constants of different complexes, the dissociation constants were derived using the potential of mean force approach of thermodynamic integration. This involves determining the change in free energy during stepwise removal of individual ligands (e.g. Cl or F). In practice, this is done by performing constraint molecular dynamics simulations in which the distance between the cation and the dissociating anion is held constant and the average force between these two ions is determined.

Repeating this procedure for different distances between the equilibrium distance of the associated complex and a distance at which the complex is considered dissociated, the potential of mean force is obtained. The free energy of dissociation is given by the integral of this potential. Stability constants describe the reverse process of stepwise association, but they are easily derived from the (standardized) dissociation energies. Finally, species distribution diagrams are constructed based on the computed species stability constants (Fig. 2). In the example presented here, the fluoride species become significant and increasingly dominant at fluoride concentrations above a few millimole per kg H₂O. Further details of this study are published in Stefanski and Jahn (2020).