Atomic Scale Modelling of Icy Surfaces: A Best Practice for Validating Interatomic Potentials and Ice Substrates in Extreme Environments

Introduction: Several current and upcoming missions will be focusing on understanding the evolution of icy surfaces and their interaction with the exosphere. For example, Europa Clipper and JUICE aim to study the icy Moons of Jupiter [1], and ongoing work on the Moon and Mercury has identified the icy permanently shadowed regions as being of scientific and operational importance. However, while significant computational and laboratory research has been conducted on volatile interactions with silicate surfaces, there is comparatively less work on icy surfaces. Understanding the interaction of icy surfaces with the space environment on these bodies is necessary for interpreting much of the upcoming observation data. For example, binding energies and diffusion characteristics of volatiles on icy surfaces can be incorporated into exospheric models to better understand the exosphere formation. While experiments are costly and time-consuming, molecular dynamics (MD) offers a theoretical alternative by simulating the behavior of atoms in extreme environments. It has proven valuable in understanding surface diffusion and surface binding energies of key volatiles on silicate surfaces [2,3]. These approaches commonly use reactive interatomic potentials (ReaxFF) that are capable of modelling dynamic bond breaking and reformation found during chemical reactions. However, while ReaxFF potentials have been well studied and validated for silicates they remain relatively untested for the conditions and compositions relevant to icy bodies. A validation of available interatomic potentials (IPs) for ice is needed before we can apply these MD-methods to exosphere modelling.

Here, we aim to validate ReaxFF for icy surfaces by comparing the mechanical properties of ice to a standard, and well validated, water-based potential called TIP4P/Ice (that cannot model chemical reactions) and available experimental results to help build a computational and methodological framework for future surface studies. We perform MD simulations of crystalline and amorphous water-ice and focus on validating against diffusion, Youngs Modulus (E_Y), isothermal compressibility (k_T), and density (ρ).

Methodology: We have studied crystalline and amorphous ice at 25 K, 100 K, and 264 K, the first two temperatures due to the relevance in shadowed craters on the Moon and Mercury and the latter for comparison against experiment. Crystalline proton disordered, non-polar, 1h ice was first created using GenIce. Amorphous ice was then made by melting the crystalline structure at 360 K and then equilibrating at 264 K while for lower temperatures for the amorphous ice was equilibrated at 245 K before quenching to 25 K and 100 K with rates ranging from 80 ns to 80 ps. Following the work by Baran et al. [4], we calculated the diffusion of Oxygen in amorphous ice as a function of temperature for 40 ns in a constant volume and temperature ensemble. For the E_Y the crystalline and amorphous substrates we produced stress-strain curves using a minimum strain rate [5]. Finally, the isothermal compressibility was computed from volume variations.

Results: First, we demonstrate a notable methodological advancement by validating the ability to convert equilibrated ice structures from TIP4P/Ice to ReaxFF formats. This allows researchers to leverage the faster TIP4P/Ice for equilibration of the surface and then switch to the more computationally intensive ReaxFF potential for the chemically reactive simulations that will be found during volatile interactions.

Due to limited experimental data, the E_Y of the two IPs was only compared to experimental data at 264 K [6]. For both IPs there was strong agreement when compared to experiment, a 5.2% and 5.6% difference for TIP4P/Ice and ReaxFF potentials respectively. As temperature increases from 25 K to 264 K, we found that the difference in E_Y between the two potentials decreased from 23% to 0.3%, suggesting that at temperatures below melting ReaxFF is well optimized. When comparing amorphous to crystalline ice, we found that the E_Y is lower by ~50% for ReaxFF for both 25 K and 100 K cases whereas for TIP4P/Ice the E_Y decreases by 28% and 42% for 25 and 100 K respectively. As expected at 264 K the sample is melted and has a Young Modulus of zero. The E_Y values for amorphous ice were found to increase with increasing the quenching time. The third ReaxFF potential showed similar trends but was more diffusive translating to a higher melting point and lower E_Y.

The calculated diffusion values for both tested IPs compared well to previous simulations that used TIP4P/Ice [4]. Isothermal compressibility for both potentials is consistent with each other but is underestimated compared to previous studies [7]. Density falls within expected ranges for both crystalline, 0.88 < ρ < 0.94, and amorphous ice, 1.05 < ρ < 1.16.

Conclusion: We tested ReaxFF potentials for key properties of icy surfaces and suggest a new validation methodological approach for future simulations. Hence, we provided a clear framework that can be reliably applied by other researchers to assess new and emerging potentials relevant to space applications. This advancement not only ensures accurate simulation of mechanical behavior of ice but also opens pathways for further exploration into icy surface chemistry. This study lays the groundwork for accelerating our understand of surface exosphere connections on icy bodies.

[1] Magnanini et al. (2024) Astronomy & Astrophysics 687 A132. [2] Morrissey, et al. (2022) Icarus, vol. 379, article no. 114979. [3] Morrissey, et al. (2022) The Astrophysical Jounral Letters 925. [4] Baran et al (2023) J. Chem. Phys. 158 (6): 064503. [5] Santos-Flórez, et al. (2018) The Journal of Chemical Physics 149.16. [6] Schulson. (1999) Jom 51 21-27. [7] Neumeier. (2018) J. Phys. Chem. 47 (3).