Processing math: 100%
Mazzucchelli, M. L., Moulas, E., Kaus, B. J. P., & Speck, T. (2024). Fluid-mineral Equilibrium Under Nonhydrostatic Stress: Insight From Molecular Dynamics. American Journal of Science, 324, 2. https://doi.org/10.2475/001c.92881
Download all (8)
  • Figure 1. (A) Schematic illustration of the decomposition of a general, nonhydrostatic stress state into a mean-stress (hydrostatic part) and a deviation (deviatoric stress). In an elastically isotropic material, the hydrostatic and the deviatoric part of the stress control the changes in volume and shape, respectively. The differential stress (σmaxσmin) describes the difference between the largest and smallest normal stresses and it is proportional to the deviation from hydrostatic conditions. (B) Under hydrostatic stress, all the normal stresses acting on a point are equal to each other. If the stress state is the same at every point in a system, the stress is said to be homogenous. If the magnitude of the hydrostatic stress varies across the system, the hydrostatic stress is nonhomogeneous. (C) Under a nonhydrostatic stress state, the normal stress components acting on a point are not equal one another. If the nonhydrostatic stress is the same at every point in the system, the system is under a homogeneous, nonhydrostatic stress. In a system under a non-homogenous, nonhydrostatic stress, the normal components can vary in magnitude, sign and orientation across the system.
  • Figure 2. Mechanical stress state of a porous system composed of a solid porous matrix and a fluid (geometry after fig. 1 of Dahlen, 1992). All quantities are presented in dimensionless form (indicated by a star, see appendix A for details). Despite the system is assumed to be under its own weight and without any far-field stress (lithostatic), the solid-fluid coexistence led to the development of local nonhydrostatic stress in the solid matrix. (A) Pressure of the interconnected fluid (Pf- colored region). (B) Negative mean stress (pressure) of the solid grains (Ps - colored region). (C) Differential stress of the solid matrix (σsd- colored region). The non-zero differential stress in the solid grains indicates that their state of stress is not hydrostatic.
  • Figure 3. (A) Snapshot of the simulated coexisting LJ system equilibrated at hydrostatic conditions. Atoms belonging to the fcc solid and fluid layers are shown in orange and blue, respectively (see appendix C for details on the distinguishing procedure). Rendering performed with the OVITO software (Stukowski, 2010). The crystallographic directions [1 0 0] and [0 1 0] of the solid are aligned with the reference axes x1 and x2, respectively. (B) Stress configuration of the solid in the MD simulations. The strain applied to the solid is axially-symmetric (ε11=ε22) resulting (because of the crystallographic symmetry and orientation of the solid) in an axially-symmetric stress (σs11=σs22). Because of mechanical equilibrium σs33=Pf. The shear stress components of the solid are zero.
  • Figure 4. Schematic representation of solid-fluid equilibrium in a PT diagram. (A) The P,T points at which the solid and the fluid coexist at equilibrium under a homogeneous and hydrostatic stress state define the usual hydrostatic phase boundary in the PT space. Examples of snapshots of the atomic configuration of the LJ fluid (p1) and fcc solid (p3) phases, and of the solid-fluid coexistence at equilibrium (p2) as calculated from MD are shown. Atoms belonging to the solid and fluid phases (see appendix C) are shown in orange and blue, respectively. The rendering was performed using the OVITO software (Stukowski, 2010). (B) When the solid is stressed in a nonhydrostatic manner, the equilibrium between the solid and the fluid is achieved at a fluid pressure Pf higher than the pressure on the hydrostatic boundary at the same temperature (e.g., p2 in A and p4 in B). For a fixed stress state of the solid, the variation of the equilibrium Pf with temperature would define a nonhydrostatic phase boundary that is shifted to higher pressures (and lower temperatures) compared to the hydrostatic boundary. Therefore, a solid can be melted under nonhydrostatic stress conditions because of the shift of the phase equilibria (e.g., p2 is a coexistence state in A but it melts in B).
  • Figure 5. Potential energy-temperature hysteresis, generated by heating and cooling at constant rate in an NPT ensemble for 30106 timesteps at the hydrostatic PH=0.5. The T+=0.737 and T=0.435 are found as the T of melting and crystallization, respectively. Because of the cooling rate and the effect of boundary conditions, defects are formed during the crystallization along the cooling path. Therefore, the solid formed from supercooled liquid has higher enthalpy than the initial solid, explaining the small difference in energy between the heating and cooling paths at T<T. The equilibrium temperature of melting at hydrostatic pressure TH=0.606 is calculated as described in Luo et al. (2004, eq. (14)). All quantities are reported as LJ reduced units.
  • Figure 6. Results of MD simulations as a function of the nonhydrostatic stress applied to the solid, with components qs11=qs22 (fig. 3B, table 4) at a constant temperature T0=0.60935 (LJ reduced units). (A) Pressure of the LJ fluid at nonhydrostatic equilibrium. (B) Deviation of the mean stress (¯σs) of the nonhydrostatic LJ solid layer from its initial mean stress (PH) at hydrostatic equilibrium. (C) Deviation of ¯σs from the negative of the equilibrium pressure (Pf) of the fluid. The solid lines in (A) and (C) are calculated from equation (14) using the properties at hydrostatic conditions reported in tables 2 and 3. The ¯σs along the solid line in (B) is calculated from the strain applied to the solid and the Pf, using linear elasticity (elastic properties in table 3). All the stress and pressure quantities are reported as LJ reduced units.
  • Figure 7. Equilibrium quantities of the LJ system under nonhydrostatic stress, as obtained from the theoretical relationship of equation (14) by varying the components q11 and q22 independently. (A) Variation of the equilibrium pressure of the fluid Pf with stress (grey surface). (B) Deviation of the mean stress (¯σs) of the nonhydrostatic solid layer from its initial mean stress (PH) at hydrostatic equilibrium, shown as contours. (C) Deviation of ¯σs from the negative of the pressure of the fluid (Pf) at nonhydrostatic equilibrium, shown as contours. Three deformation paths are shown as lines. Red solid line: path with q11=q22, as explored in our MD simulations (see fig. 6). Red dashed line: path with q11=q22; along this path the ¯σs = Pf (C), but deviates from the initial PH (B). Black dashed line: path along which the ¯σs = PH (B), but deviates from Pf (C). All the shown stress and pressure quantities are reported as LJ reduced units. The properties used in the calculations are reported in tables 2 and 3.
  • Figure 8. Melt-solid equilibrium as a function of the nonhydrostatic stress applied to the solid with components qs11=qs22 at constant T. The initial hydrostatic PH, TH conditions of Mg2SiO4 (2.70 GPa,2304 K), SiO2 (1.49 GPa,2262 K), NaCl (0.95 GPa,1264 K) and the LJ material (0.5ε/σ3,0.609ε/kB) are determined from their respective melt lines at the nondimensional pressure P=P/K=0.0253, as explained in the main text. All of the quantities are normalized on the initial hydrostatic PH of the system. (A) Variation of the melt equilibrium pressure calculated from equation (B2) assuming isotropic elastic properties for the solid. The dashed line (anisotropic LJ) has been calculated with equation (14) assuming anisotropic elastic properties for the solid. (B) Deviation of the mean stress (¯σs) of the nonhydrostatic solid layer from its initial mean stress (PH) at hydrostatic equilibrium. (C) Deviation of ¯σs from the negative of the equilibrium pressure (Pf) of the fluid. The properties of each phase at hydrostatic conditions are reported in tables 2, 3 and 5.

Abstract

The interpretation of phase equilibria and reactions in geological materials is based on standard thermodynamics that assumes hydrostatic and homogeneous stress conditions. However, rocks and minerals in the lithosphere can support stress gradients and nonhydrostatic stresses. Currently, there is still not an accepted macroscopic thermodynamic theory to include the effect of nonhydrostatic stress on mineral reactions, and the use of several thermodynamic potentials in stressed geological system remains under debate. In experiments under nonhydrostatic stress, it is often difficult to resolve the direct effect of differential stress on phase equilibria because pressure gradients may be developed. Such gradients can affect the metamorphic equilibria at the local scale. Here, we investigate the direct effect of a homogeneous, nonhydrostatic stress field on the solid-fluid equilibrium using molecular dynamics simulations at non-zero pressure and elevated temperature conditions. Our results show that, for simple single-component systems at constant temperature, the equilibrium fluid pressure of a stressed system is always larger than the value of fluid pressure at hydrostatic stress conditions. The displacement of the equilibrium value of the fluid pressure is about an order of magnitude smaller compared to the level of differential stress in the solid crystal. Thus, phase equilibria can be accurately predicted by taking the fluid pressure as a proxy of the equilibration pressure. On the contrary, the mean stress of the solid can deviate substantially from the pressure of the fluid in stressed systems at thermodynamic equilibrium. This has implications on the use of thermodynamic pressure in geodynamic models since the fluid pressure is a more accurate proxy for predicting the location of metamorphic reactions, while the equilibrium density of the solid has to be determined from its mean stress.

1. INTRODUCTION

Our understanding of geological processes requires that we can accurately model how minerals, rocks and fluids behave chemically and mechanically in a wide range of pressure, temperature and stress conditions. Commonly, such conditions are characterized by increased pressure and temperature (P-T) compared to the Earth’s surface, thus making the direct observation problematic for processes like mineral deformation, fluid-mineral interaction and phase transformations in general. For example, the understanding of geodynamic processes such as plate subduction and mountain building largely relies on indirect observations such as the investigation of rocks that have recrystallized in subduction zones and are now found at the Earth’s surface. The preserved (quenched) mineral assemblages in metamorphic rocks are frequently interpreted in the framework of equilibrium thermodynamics, based on the experimentally produced mineral assemblages at laboratory conditions (e.g., Bucher & Frey, 2002; Spear, 1993). This approach has given us important insights on the range of conditions that are experienced by igneous and metamorphic rocks in natural environments. However, a major assumption of experimental petrology is that the stress state is hydrostatic and spatially homogeneous (e.g., fig. 1A, B), i.e., all the normal stresses in coexisting minerals and fluids are equal. It is important to note that by “hydrostatic stress state” we refer to the isotropy of the stress tensor (fig. 1A) and we do not imply that the stress can be calculated using the lithostatic stress gradient.

Figure 1
Figure 1.(A) Schematic illustration of the decomposition of a general, nonhydrostatic stress state into a mean-stress (hydrostatic part) and a deviation (deviatoric stress). In an elastically isotropic material, the hydrostatic and the deviatoric part of the stress control the changes in volume and shape, respectively. The differential stress (σmaxσmin) describes the difference between the largest and smallest normal stresses and it is proportional to the deviation from hydrostatic conditions. (B) Under hydrostatic stress, all the normal stresses acting on a point are equal to each other. If the stress state is the same at every point in a system, the stress is said to be homogenous. If the magnitude of the hydrostatic stress varies across the system, the hydrostatic stress is nonhomogeneous. (C) Under a nonhydrostatic stress state, the normal stress components acting on a point are not equal one another. If the nonhydrostatic stress is the same at every point in the system, the system is under a homogeneous, nonhydrostatic stress. In a system under a non-homogenous, nonhydrostatic stress, the normal components can vary in magnitude, sign and orientation across the system.

An interesting realization is that all Gibbs-energy minimization algorithms and mineral equilibrium software (e.g., Connolly, 2005, 2017; de Capitani & Petrakakis, 2010; Riel et al., 2022) rely on this hydrostatic approximation. This means that in these algorithms the pressure of each phase is the same (i.e., it is homogeneous) and the stress tensor is assumed to be isotropic independently of the nature of the material and the loading configuration (fig. 1B). In mechanically weak mineral assemblages, the hydrostatic approximation seems to provide an adequate description of the state of stress in the Earth, because all of the normal stress components will be almost equal to each other (e.g., Moulas et al., 2019). However, because of its rheological properties, the Earth’s lithosphere can sustain significant levels of nonhydrostatic stress with macroscopic differential stresses that may reach up to 0.5–1 GPa (Brace & Kohlstedt, 1980; Prieto et al., 2012). In this work, we use the term nonhydrostatic stress to refer to the general state of stress where the stress tensor is not isotropic. For isotropic materials it is common to refer to the isotropic part of the stress tensor as hydrostatic stress (i.e., the pressure, when adjusted for sign conventions), whereas the anisotropic part is referred to as the deviatoric stress (fig. 1A). In principle, for elastically isotropic materials, the pressure is solely related to volume changes while the deviatoric stress induces only shape changes. However, for a general anisotropic material (i.e., for any crystal with crystallographic symmetry lower than cubic) even a purely hydrostatic state of stress will not result in a purely isotropic strain state. Therefore, in solid minerals, volume changes are generally influenced also by the deviatoric part of the stress. In this case, the decomposition of the stress tensor into a hydrostatic and deviatoric part does not hold a precise physical meaning (Augusti et al., 1969). To avoid confusion, in this work we use the term “nonhydrostatic stress” to describe a general state of stress in which the stress loading is not equal from every side of the crystal. By differential stress we refer to the difference between the maximum and the minimum normal components (eigenvalues) of the stress tensor, which scales proportionally to the degree of anisotropy of the stress tensor (fig. 1A).

In presence of fast tectonic deformations, the differential stress in the Earth’s lithosphere can reach several hundreds of MPa in magnitude (e.g., Andersen et al., 2008). Nonhydrostatic stresses can also be observed when relatively strong viscous creep flow laws are employed in the modeling of subduction settings (e.g., Reuber et al., 2016). At the mineral scale, nonhydrostatic stress states with differential stress up to 1 GPa are observed when strong mineral hosts enclose minerals with different elastic properties (Campomenosi et al., 2018; Gonzalez et al., 2021; Kohn et al., 2023; Moulas et al., 2020, 2023). In general, any change in the pressure and temperature experienced by a solid mineral aggregate will generate stresses at the grain-scale level simply because different minerals respond differently to the applied stress/strain conditions. In recent years, an increased number of studies have highlighted the importance of volume changes during mineral reactions and phase transitions (Gillet et al., 1984; Van der Molen & Van Roermund, 1986; Zheng et al., 2018).

These effects become particularly significant in a confined environment such as the one expected in the deep Earth. In such cases, mechanical models predict the development of GPa-level stress and pressure differences between the host mineral and the inclusion (Plümper et al., 2022; Zhukov & Korsakov, 2015). Moreover, in multiphase systems with coexisting minerals and fluids, pressure and stress gradients are intrinsically generated due to changes in the contact area of the solid grains (e.g., Dahlen, 1992; Gallagher et al., 1974; Pollard & Fletcher, 2005, p. 194). However, despite the quantification of the interplay between stress, thermodynamic equilibria and reactions being essential for understanding many key geological processes, the magnitude of the effect of nonhydrostatic stresses on fluid-mineral and mineral-mineral equilibria is still actively debated.

Gibbs (1876) first realized that the classic expression of free energy and phase equilibria was only applicable to hydrostatic condition and discussed the influence of a nonhydrostatic stress on the equilibrium between a solid and a fluid of the same composition. However, the magnitude of the effect of nonhydrostatic stresses on fluid-mineral and mineral-mineral equilibria is still actively debated among researchers who, sometimes, have incompatible views about the use of various thermodynamic formulations (e.g., Hess & Ague, 2021; Hobbs & Ord, 2015; Tajčmanová et al., 2015; Wheeler, 2014). Several experiments have tried to shed light on such controversies by imposing nonhydrostatic stress on pure substances that experience phase transitions (Hirth & Tullis, 1994; Hobbs, 1968; Richter et al., 2016; Vaughan et al., 1984; Zhou et al., 2005). However, recent studies that combine data from experiments and models of rock deformation suggest that the stress field in such experiments is nonhomogeneous, and therefore the experimental results are not diagnostic (Cionoiu et al., 2019, 2022; Ji & Wang, 2011; Moulas et al., 2013). The development of pressure gradients in these experiments is directly proportional to the applied far-field differential stress (Cionoiu et al., 2022). Thus, in mechanical experiments, it is difficult to resolve whether the observed phase transitions are controlled by the nonhydrostatic stress state (which in this work we call the direct effect of nonhydrostatic stress) or if they are caused by local pressure concentrations due the presence of a nonhomogeneous stress within the sample (i.e., the indirect effect of nonhydrostatic stress, fig. 1C). Apart from the problem of nonhomogeneous stress distribution, typical rock-deformation experiments do not allow us to resolve the stress state within the sample during the occurrence of phase transitions and mineral reactions. The products of such experiments are typically studied after the sample has been recovered and the only available information about its evolving stress state is related to the far-field loading conditions. Therefore, the actual phase distribution at the end of the experiment may not correspond to an instantaneous equilibrium phase distribution but may be influenced by the entire stress/strain history of the sample (e.g., Cionoiu et al., 2022). For this reason, without knowledge of the stress distribution within the sample, it is extremely difficult to assess the direct effect of stress on mineral equilibria.

The presence of nonhydrostatic stress conditions and stress gradients is not just limited to solid aggregates, but they can also be developed in porous system where solid grains are in contact and react with fluid phases, such as aqueous fluids and melts (Gallagher et al., 1974). Fluid flow and associated hydration and dehydration reactions play a fundamental role in many phenomena in the lithosphere, such as the dynamics of shear zones (e.g., Austrheim, 1987), the generation of earthquakes (e.g., Brantut et al., 2017; Ferrand et al., 2017), or rheological weakening of rocks (e.g., Jolivet et al., 2005; Kaatz et al., 2021; Moulas et al., 2022). In magmatic systems, the production of melt due to phase transitions can have a dramatic influence on the rheological properties of porous rocks and on the overall magma dynamics (e.g., Hirth & Kohlstedt, 1995; Schmeling et al., 2012; Yarushina & Podladchikov, 2015). The accurate prediction of fluid-solid reactions is also central in industrial processes with examples that range from the prediction of mineral reactions in geological carbon sequestration (e.g., Kelemen & Hirth, 2012), to the assessment of the volume changes during geothermal energy extraction (Erfani et al., 2019) and to the engineering of underground long-term nuclear waste repositories (Spycher et al., 2003).

While in hydrostatic thermodynamics the Gibbs potential can be adopted to find equilibrium assemblages, there is no experimentally verified, thermodynamic theory for nonhydrostatic reactions in deforming geomaterials. To assess which thermodynamic description is the most appropriate in describing mineral phase equilibria under nonhydrostatic stresses, we performed atomistic computer simulations with molecular-dynamics (MD). With this approach, we can evaluate the evolution of the trajectories of the individual atoms that constitute a two-phase solid-fluid system, under prescribed stress and temperature conditions. With MD simulations, the direct influence of several controlled nonhydrostatic stress states on the mineral-fluid equilibrium can be assessed independent of any assumption regarding the thermodynamic potential of the system. Moreover, in MD simulations the stress in the solid phase can be kept homogeneous, allowing us to isolate the direct effect of the nonhydrostatic stress on the perturbation of the equilibrium, avoiding the development of local stress concentrations which, as deduced by experiments, may indirectly affect the equilibrium. For our MD simulations we adopted the Lennard-Jones (LJ) method (Lennard-Jones, 1924a, 1924b), which provides a numerical representation of an idealized network of atoms and the interatomic potential that controls the dynamics of the individual atoms. This allowed us to explore the equilibrium states of a coexisting solid and fluid in an idealized one-component system under controlled stress conditions.

We begin our paper by introducing the problem of stress gradients in fluid-mineral systems and then proceed with a short summary of the previous thermodynamic analysis related to the single-component solid-fluid equilibration. We then introduce our MD models in relation to the problem of solid-fluid equilibration under stress. Our approach follows similar studies from the literature but considers the effects of nonhydrostatic stress under confined conditions (i.e., initial pressure larger than zero). Our simulations show that, at constant temperature, a positive change of the fluid equilibrium pressure is always to be expected, whether the applied nonhydrostatic stress is compressive or tensile. For small applied differential stresses, the change in fluid pressure (Pf) required to maintain the solid-fluid equilibrium at constant temperature is small and approximately one order of magnitude less than the differential stress generated in the solid. On the other hand, our results show that, for the particular loading configuration, the mean stress of the solid (¯σs) changes considerably in response to the nonhydrostatic stress applied to it. Therefore, when the solid is under differential stress, its ¯σs deviates from the Pf with which it is in equilibrium, leading to a decoupling between the “pressure” (i.e., ¯σs) of the solid phase and the pressure of the solid-fluid thermodynamic equilibrium. With this approach, we can therefore assess the magnitude of the deviation of phase equilibria from the hydrostatic approximation in simple systems under nonhydrostatic stress. Finally, we discuss the magnitude of the stress inhomogeneities that arise in solid-fluid multiphase systems and which stress measure should be taken as a proxy of the thermodynamic pressure in order to accurately describe reactions.

2. STRESS IN FLUID-SOLID MULTIPHASE SYSTEMS

Many processes in the Earth are intrinsically related to the evolution of multiphase systems in which solid mineral phases are coupled with fluid phases, such as aqueous fluids and melts. Theoretical considerations and experimental studies have verified that in porous and granular rocks stresses are not distributed homogeneously (Gallagher et al., 1974). In particular, a granular rock having a fluid-filled porosity is expected to develop stress gradients because of the inability of the fluid to sustain large deviatoric stresses. Even without far-field tectonic loads, stress gradients are developed due to the presence of gravity and the occurrence of stress concentration along grain contacts. We refer to such conditions as quasi-static since true lithostatic conditions require, inter alia, that the rock is stratified with respect to its density (e.g., Moulas et al., 2019). Porous rocks with lateral density gradients do not satisfy this condition but, in the absence of far-field stresses can be approximated as static, at least macroscopically. In order to calculate and show the mean-stress (pressure) gradients developed in a granular rock under the influence of gravity we consider a schematic two-dimensional section of a model rock as it was provided in the work of Dahlen (1992, his fig. 1). In addition, we adopt a viscous mechanical model with a gravity field (the details of the calculation are reported in appendix A). Figure 2A shows that when the porosity is interconnected, the fluid will always experience a hydrostatic stress gradient. On the contrary, the solid grains will experience different levels of stress due to the changes in the contact area between each grain (e.g., Pollard & Fletcher, 2005, p. 194). As a consequence, in rocks that are porous or have a granular structure, local stress gradients are always developed, even if their macroscopic pressure gradient is “lithostatic” (fig. 2B, C). The presence of stress and pressure gradients may have important implications on our interpretations of processes in the Earth.

Figure 2
Figure 2.Mechanical stress state of a porous system composed of a solid porous matrix and a fluid (geometry after fig. 1 of Dahlen, 1992). All quantities are presented in dimensionless form (indicated by a star, see appendix A for details). Despite the system is assumed to be under its own weight and without any far-field stress (lithostatic), the solid-fluid coexistence led to the development of local nonhydrostatic stress in the solid matrix. (A) Pressure of the interconnected fluid (Pf- colored region). (B) Negative mean stress (pressure) of the solid grains (Ps - colored region). (C) Differential stress of the solid matrix (σsd- colored region). The non-zero differential stress in the solid grains indicates that their state of stress is not hydrostatic.

Current numerical models that describe the mechanical and chemical evolution of such systems assume that the fluid pressure is equal to the rock pressure (e.g., Yang & Faccenda, 2020) and that the chemical reactions are controlled by the pressure of the fluid phase (e.g., Schmalholz et al., 2020). From a thermodynamic perspective, these assumptions imply that the solid matrix and the coexisting fluid are assumed to be under hydrostatic-stress conditions. Therefore, within this approximation, the resulting phase equilibria and reactions are described by the minimization of Gibbs free energy at constrained PT conditions. However, as we have shown, the presence of a nonhomogeneous-stress distribution at the grain scale breaks the fundamental assumptions behind the use of Gibbs potentials, i.e., a hydrostatic stress tensor and homogeneous stress field. The presence of nonhydrostatic stress can lead to changes of the physical conditions at which phase equilibria and reactions take place. These changes can have a large effect on the chemical and mechanical evolution of the system. For example, the production of fluid/melt during mineral reactions and phase transitions can have a dramatic influence on the rheological properties of porous rocks. This is because the effective bulk and shear viscosity of partially molten rocks are inversely dependent on porosity and melt fraction (e.g., Schmeling et al., 2012; Yarushina & Podladchikov, 2015). Such effects become even more relevant if we consider cases in which the fluid/melt flow is driven by pressure gradients in tectonically stressed, porous rocks (e.g., Connolly & Podladchikov, 2004; Spiegelman & McKenzie, 1987).

Therefore, the quantification of the effect of nonhydrostatic stress on mineral phase equilibria is crucial to describe the mechanical properties and the mineralogy of porous, multiphase, reactive rocks. However, such perturbations cannot be evaluated directly in natural systems since, as shown in figure 2, they do not provide simple example of homogeneous loading in the solid phase. This implies that, similarly to what observed in experiments, also in natural samples the displacements of thermodynamic equilibria due to the nonhydrostatic stress state (direct effect) cannot be separated from the effect due to stress concentrations and pressure gradients (indirect effect). Such perturbations cannot be evaluated either by existing numerical models which are purely mechanical, and therefore cannot account for the two-way feedback between the deformation and the nonhydrostatic thermodynamic response. Models that are able to assess the microscopic scale behavior of multiphase systems under homogeneous stress in the solid phase, such as those that can be performed with MD, are therefore required to assess the direct effect of a homogeneous, nonhydrostatic stress on the thermodynamic phase equilibria in deformed systems.

3. THERMODYNAMICS OF SOLID-FLUID EQUILIBRIUM AT NONHYDROSTATIC CONDITIONS

The thermodynamics of systems which have a nonhydrostatic stress state has been the topic of several recent studies (Hess & Ague, 2021; Hobbs & Ord, 2016; Tajčmanová et al., 2015; Wheeler, 2014). In those studies, different interpretations are proposed regarding the role of pressure (also called “thermodynamic pressure”). In this work, we will highlight the difference between the fluid and solid pressures in a thermodynamically equilibrated system. As it will become clear later, the fluid and solid pressures can both be considered “thermodynamic” as they can be derived from thermodynamic principles. The main difference is that the equilibrium fluid pressure is related to the solid-fluid coexistence conditions while the solid “pressure” (¯σs) is related to the volumetric response of the solid.

In systems under nonhydrostatic conditions that comprise solids in contact with fluids, chemical reactions can involve dissolution of solid phases in the fluid, chemical reactions in the fluid and subsequent precipitation of the newly reacted material on the stressed solid. The equilibrium between a solid under nonhydrostatic stress and fluids was first investigated by Gibbs (1876, pp. 34–520). Assuming that temperature is homogeneous in the system and that the deformations are small, Gibbs derived the equilibrium between a solid and a fluid in a single component system (his eq 386 and eq 3 in Sekerka & Cahn, 2004):

ˆUsˆΩs0TˆSsˆΩs0+PfˆΩsˆΩs0=ˆμfˆρs

where T (K) is the absolute temperature of the system at equilibrium. Superscript s refers to variables of the solid, where ˆΩs0 (in m3) is the reference volume of the solid in the undeformed configuration, ˆUs its energy (J), ˆSs its entropy (J/K), and ˆρs=m^Ω0s (kg/m3) the mass density of the considered component in the solid at the reference state (notation is summarized in table 1). Superscript f refers to variables of the fluid, with Pf (Pa) being the pressure of the fluid, and ˆμf (J/kg) being the chemical potential of the substance of the solid in the fluid. By considering the atom density of the substance of the solid in the undeformed state ρsa=NˆΩs0 (m-3) and expressing the per-atom chemical potential as μf, we can divide equation (1) by ρsa to obtain:

usTss+Pfωs=μf

where us=ˆUsˆΩs0 ρsa , ss=ˆSsˆΩs0 ρsa , ωs=ˆΩsˆΩs0 ρsa , μf=ˆμfˆρsρsa are now per-atom quantities.

Table 1.List of Symbols and Abbreviations.
Symbol/Abbreviation Meaning Units
P Pressure Pa or
ˆa/ˆb3 in LJ systems
T Temperature K or ˆa/kB in LJ systems
η Effective viscosity Pa s
ρ Mass density kg/m3
K Bulk modulus Pa or
ˆa/ˆb3 in LJ systems
σij Cauchy stress tensor Pa or
ˆa/ˆb3 in LJ systems
τij Deviatoric stress tensor Pa
Vi Velocity m/s
PH Pressure of the solid and fluid at hydrostatic equilibrium Pa
TH Temperature at hydrostatic equilibrium K
Pf Fluid pressure Pa
ˆUs Solid energy J
ˆSs Solid entropy J/K
ˆΩs0 Solid reference volume m3
ˆρs Mass density of the substance of the solid in the undeformed state kg/m3
ˆμf Chemical potential of the substance of the solid in the fluid J/kg
ρsa Atom density of the substance of the solid in the undeformed state atoms/m3
us Solid energy, per atom J/atom
ss Solid entropy, per atom J/K/atom
ωs Solid volume, per atom m3/atom
μf Chemical potential of the substance of the solid in the fluid, per atom J/atom
ψs Helmholtz free energy of the solid, per atom J/atom
gsH Gibbs free energy of the solid at hydrostatic conditions, per atom J/atom
μfH Chemical potential of the substance of the solid in the fluid at hydrostatic conditions, per atom J/atom
εsij Small-strain tensor in the solid None
sfH Fluid entropy at hydrostatic conditions, per atom J/K/atom
ssH Solid entropy at hydrostatic conditions, per atom J/K/atom
ωfH Fluid volume at hydrostatic conditions, per atom m3/atom, or
ˆb3/atom in LJ systems
ωsH Solid volume at hydrostatic conditions, per atom m3/atom, or
ˆb3/atom in LJ systems
Csijkl Elastic stiffness tensor of the solid Pa, or
ˆa/ˆb3 in LJ systems
Dsijkl Elastic compliance tensor of the solid Pa-1 or
(ˆa/ˆb3)1 in LJ systems
qsij = σsij+δijPf, deviation of the stress of the solid from the Pf Pa or
ˆa/ˆb3 in LJ systems
¯σs Solid mean stress Pa or
ˆa/ˆb3 in LJ systems
MD molecular dynamics --
LJ Lennard-Jones interaction --
a Parameter of the LJ interaction (minimum potential energy) ˆa
b Parameter of the LJ interaction (distance at which the potential energy is zero) ˆb
Rc Cutoff distance of the LJ interaction ˆb
m Particle mass of the LJ substance ˆm
Δt Time step of the MD simulations (ˆmˆb2ˆa)12
˜l = l/b, reduced units of length in the LJ system None
˜t =t(am b2)12, reduced units of time in the LJ system None
˜P =Pb3a, reduced units of pressure in the LJ system None
˜σ =σb3a, reduced units of stress in the LJ system None
˜T =kBTa, reduced units of temperature in the LJ system. None
NPT MD ensemble with constant: N, number of particles; P, hydrostatic pressure; T, temperature --
NPx3AT MD ensemble with constant: N, number of particles; Px3, total stress along the x3 axis; A, area in the plane normal to x3; T, temperature --
NPH MD ensemble with constant: N, number of particles; P, hydrostatic pressure; H, enthalpy --
NVT MD ensemble with constant: N, number of particles; V, volume; T, temperature --
P’ =P/K, nondimensional pressure used to rescale result of different solid-melt systems None

Starting from the work of Gibbs (1876), Frolov and Mishin (2010) have explored in more detail the case of a single component solid in equilibrium with its liquid (i.e., a single-component fluid with the same composition of the solid). Their analysis has extended the theory to include also the anisotropic elastic properties of the solid phase, which is particularly relevant for minerals that are elastically anisotropic (e.g., Nye, 1985). Here we follow Sekerka and Cahn (2004) and Frolov and Mishin (2010) and develop a relationship that describes the effect of nonhydrostatic stresses on the solid-fluid equilibrium in single-component systems. We start by recognizing the expression of the Helmholtz free energy of the solid, ψs, in equation (2):

ψs= usTss

Therefore, using the previous equation, equation (2) can be written as:

ψs+Pf ωs=μf

When both the solid and the fluid are under a homogeneous and hydrostatic stress, the pressure of the fluid Pfis equal to the pressure of the system PH. In that case, the Gibbs free energy (per atom) of the solid is:

gsH= ψsH+PHωsH

where the subscript H refers to the properties in the system where both the solid and the fluid are under the same hydrostatic conditions. By rearranging equation (4) and substituting the result in equation (5), we arrive to the usual equilibrium condition for a single component solid in contact with its pure melt under hydrostatic stress:

gsH=μfH

When the solid is taken from the hydrostatic state to nonhydrostatic conditions its energy changes, and therefore, in order to preserve the equilibrium, the chemical potential of the component of the solid in the fluid must change accordingly. The variation of the Helmholtz energy of the solid with stress/strain can be evaluated by expressing its differential with respect to the applied strain and the temperature:

dψs=ωsHi,j=1,2,3σsijdεsijssdT

where ωsH is the volume per atom of the solid in the hydrostatic reference state used to define the strain, and σsij and εsij are the Cauchy stress tensor and the small strain (i.e., Lagrangian infinitesimal) tensor of the solid, respectively. Following the usual convention in continuum mechanics, σsij and εsij are positive in tension. Equation (7) shows that, for the single component solid, when the system is kept at constant temperature, the variation of the Helmholtz free energy is only due to the applied strain. Therefore, for an isothermal process, the Helmholtz free energy of the solid is equal to the elastic energy, and any change in the strain state of the solid will increase its free energy. This increase in the specific (per atom) elastic energy occurs whether the strain in the solid is tensile or compressive.

Assuming that the solid is linear elastic and that the applied strain is small, Hooke’s law can be used to express the linear relationship between the Cauchy stress and the applied strain: σij= Cijklεkl, with Cijkl the 4th-rank elastic stiffness tensor of the solid (with components in pressure units, Pa). Under this assumption, the specific Helmholtz free energy of the solid can be found by integrating equation (7) for small variations in strain (assuming constant temperature):

ψs=ωsH2i,j=1,2,3σsijεsij=ωsH2i,j,k,l=1,2,3Csijklεsijεskl= ωsH2i,j,k,l=1,2,3Dsijklσsijσskl

where Dsijkl= Cs1ijkl is the elastic compliance tensor of the solid. At constant temperature, a solid can increase its energy by increasing its elastic stress, and therefore the Helmholtz energy (eq 8) of the linear-elastic solid is an implicit function of its elastic stress state (not necessarily hydrostatic) and of the absolute temperature, independent from the mechanism that has generated the stress. In other words, given a specific temperature and stress state, the specific Helmholtz energy of the solid is the same whether the solid is in contact with a fluid or with another solid.

Since we are interested in describing the equilibrium between a solid and a fluid, we can exploit the requirements of mechanical equilibrium, which imposes that the normal stress of the solid orthogonal to the solid-fluid interface is equal in magnitude and opposite in sign to the pressure of the fluid. For example, we can consider a configuration in which a rectangular block of the solid is in contact with a fluid along two faces and is oriented so that the principal stress component σ33, along the Cartesian x3 direction, is orthogonal to the solid–fluid interface (e.g., fig. 3). In this case it follows that σ33=Pf. If we now assume that the shear stresses are zero, and the stress state is homogeneous in the solid, we can represent the specific Helmholtz energy of the solid as a function of this stress and the absolute temperature:

ψs=ψs (T, σ11, σ22, σ33 )=ψs (T, σ11, σ22, Pf )

Figure 3
Figure 3.(A) Snapshot of the simulated coexisting LJ system equilibrated at hydrostatic conditions. Atoms belonging to the fcc solid and fluid layers are shown in orange and blue, respectively (see appendix C for details on the distinguishing procedure). Rendering performed with the OVITO software (Stukowski, 2010). The crystallographic directions [1 0 0] and [0 1 0] of the solid are aligned with the reference axes x1 and x2, respectively. (B) Stress configuration of the solid in the MD simulations. The strain applied to the solid is axially-symmetric (ε11=ε22) resulting (because of the crystallographic symmetry and orientation of the solid) in an axially-symmetric stress (σs11=σs22). Because of mechanical equilibrium σs33=Pf. The shear stress components of the solid are zero.

We can now consider the equilibrium between a nonhydrostatically stressed solid and a fluid. Gibbs (1876) showed that the condition outlined above (i.e., eq 4) for the equilibrium between a single-component, homogeneous, hydrostatic solid and a fluid is formally valid even when the solid is loaded in a nonhydrostatic manner. However, the quantities that appear in the equation assume different values in the hydrostatic and the nonhydrostatic configuration, since they directly depend on the stress configuration. Assuming that the temperature T is uniform in the system and that the geometrical configuration is as in figure 3, we can rewrite the equilibrium condition to highlight which quantities are functions of the stress state of the solid (eq 9 of Sekerka & Cahn, 2004):

μf (T, Pf)= us(T, σ11, σ22, Pf) T ss(T, σ11, σ22, Pf)+ Pf ωs(T, σ11, σ22, Pf)= ψs(T, σ11, σ22, Pf)+ Pf ωs(T, σ11, σ22, Pf)

where T is the temperature of equilibrium between the fluid and the nonhydrostatic solid. Equation (10) highlights that the volume of the solid, and therefore its density, is controlled by its full stress state and not only by the pressure of the fluid to which it is in contact. It also shows that, in agreement to equation (8), the equilibrium between the solid and the fluid is determined by the complete stress state of the solid and not by the Pf alone. For example, in a system at constant temperature, a nonhydrostatic stress applied to the solid will always increase its specific Helmholtz energy, whatever the strain/stress state (see eq 7). Therefore, one can already expect qualitatively from equation (10) that at constant temperature the chemical potential of the component of the solid in the fluid (μf) must vary in order to preserve the equilibrium between the stressed solid and the fluid. In a single-component system at constant T, the variation of μf can only be generated by a variation of the pressure Pf. The latter implies that the fluid must change its pressure in order to preserve the equilibrium with a stressed solid.

Frolov and Mishin (2010) have derived a full description that relates the changes in temperature of the system, pressure of the fluid, and nonhydrostatic stress of an elastically anisotropic solid for a solid-fluid system at equilibrium. Such a relationship can be obtained by starting from the expressions of the differentials of the internal energy of the solid and of the chemical potential of the fluid for a variation of stress state and entropy of the solid and of temperature and pressure of the fluid, respectively:

dus=THdss+ωsHi,j=1,2,3 σsij dεsijdμf= sfHdT+ωfHdPf

where the subscript H refers to quantities evaluated in the hydrostatic state as before. The equation that describes the change in the equilibrium Pf and temperature of the system for small departures from the hydrostatic state, is obtained by taking the full differential of equation (10) and substituting the relationships of equation (11). By linearizing the differentials of the resulting equation (first-order Taylor expansion), and assuming a linear elastic relationship between strain and stresses, Frolov and Mishin (2010) integrated the resulting equation from the initial hydrostatic state to an (nonhydrostatic) equilibrium state:

ΔsH(TTH)+12(ΔsT)H(TTH)2ΔωH(PPH)12(ΔωP)H(PPH)2(ΔωT)H(TTH)(PPH)+ωsH2i,j,k,l=1,2Dsijklqsijqskl=0

ΔsH= sfHssH and ΔωH= ωfHωsH are, respectively, the differences between the atomic entropy and volume of the fluid and the solid in the reference hydrostatic state. The tensor qsij= σsij+δijPf expresses the deviation of the stress components of the solid from the Pf. We note that the tensor qsij does not represent the deviatoric stress tensor of the solid phase, since it relates the components of the stress tensor of the solid to the Pf, and not to the ¯σs. As it will be shown later, under nonhydrostatic stress the ¯σs and the Pf can assume different values even at equilibrium, leading to a deviation of the components qsij from those of the deviatoric stress tensor of the solid phase.

For our geometric configuration (fig. 3) and since the fluid cannot support shear stresses, the shear stress components in the x1x2 plane are parallel to the solid-fluid interface implying that the components σs13, σs23, and their symmetric equivalents σs31, σs32, are zero. Moreover, as discussed above, mechanical equilibrium prescribes that one of the normal components of the stress of the solid must be numerically equal to the negative of the fluid pressure (in our case σs33=Pf, see fig. 3B). Such mechanical constraints limit the summation in equation (12) to i,j,k,l = 1,2, since all the components of the tensor qsmn with indices m = 3 and/or n = 3 become zero. We note that, despite the simplification, this remains a completely general 3D problem because we do not impose any arbitrary restriction on the value of the stress components.

Equation (12) is obtained by assuming that the solid is anisotropic and linear elastic, the deformations are small (i.e., negligible volume changes) and the elastic properties of the solid are constant with stress and strain. Starting from equation (12), one can evaluate the effect of the stress state of the solid on the T and the Pf at equilibrium for any thermodynamic path. For example, at isobaric conditions (i.e., the pressure of the fluid is kept constant) and neglecting second order terms, equation (12) becomes:

TTH=  ωsH2 ΔsH i,j,k,l=1,2Dsijklqsijqskl

Equation (13) implies that the temperature of equilibrium changes from the initial hydrostatic state (TH) to the new equilibrium state (T) at which the solid is a nonhydrostatically stressed. The change in equilibrium temperature is quadratic in the qsij components, and the resulting equilibrium surface is a paraboloid in the coordinates T, q11, q12, and q22. In absence of shear stress in the solid (q12=0), the coexistence surface becomes a two-dimensional paraboloid.

In the same way, a relationship can be found between the change of equilibrium Pf and the stress of the solid at isothermal conditions. Again, neglecting second-order terms and at constant temperature, equation (12) becomes:

PfPfH= ωsH2 ΔωHi,j,k,l=1,2Dsijklqsijqskl

As for equation (13), equation (14) also implies that the pressure of the fluid changes when the solid is nonhydrostatically stressed in order to maintain the thermodynamic equilibrium between the solid and the fluid. The pressure change is quadratic in nonhydrostatic stress, and in absence of shear stress in the solid, the coexistence surface becomes a two-dimensional paraboloid.

Sekerka and Cahn (2004) derived relations equivalent to equations (13) and (14) valid for elastically isotropic solids (appendix B). While minerals are elastically anisotropic (Nye, 1985), the anisotropic elastic tensors of minerals at high P and T conditions are often not known, implying that equations (13) and (14) generally cannot be directly evaluated in mineral systems. On the other hand, the elastic properties, the molar volumes and the entropies of minerals can be calculated at high hydrostatic P and T from thermodynamic databases (e.g., Holland & Powell, 2011), under the approximation that the phases involved are elastically isotropic. In this case, as it will be discussed later in section 6.1, the isotropic equations (B1) and (B2) in appendix B can be adopted to estimate the variations of the equilibrium T and P in stressed geological systems. In appendix B we also assess the deviations expected from the use of the isotropic approximation for linear-elastic anisotropic minerals.

4. MOLECULAR DYNAMICS SIMULATIONS

For further insights on the influence of nonhydrostatic stress on the solid-fluid equilibrium, we have investigated the equilibrium conditions (PfT) of a stressed solid in contact with a fluid of the same composition. We used the software package Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS; Thompson et al., 2022) to perform classical molecular dynamics simulations where the coexistence of the two phases (solid and fluid) is monitored as a response to changes in the thermodynamic variables T, σ11= σ22 and Pf.

In our simulations we adopted the Lennard-Jones (LJ) pair potential, which is a simplified model that describes the interactions between simple atoms (Lennard-Jones, 1924a, 1924b). The LJ potential can be viewed as describing the interaction of a fictive chemical element (sometimes referred to as ‘Lennard-Jonesium’), whose thermodynamic properties and phase transitions have been extensively investigated by analytical and numerical computations (e.g., Schultz & Kofke, 2018; Stephan et al., 2019). The LJ model is a special case of the Mie potential (Grüneisen, 1926; Mie, 1903; Schwerdtfeger et al., 2021), from which the Mie-Grüneisen equation of state was developed and that is used to describe the PT dependence of volume of solid materials and of minerals at geological conditions (e.g. Anderson, 1995; Angel et al., 2022; Heuzé, 2012; Speziale et al., 2001). Despite its simplicity, this model represents physically realistic interactions for fluid and condensed systems, making it a common choice for the development of theories and methods concerning the properties of matter (e.g., Jorgensen et al., 1996; Luo et al., 2004; Shukla & Cowley, 1985; Yu et al., 2022).

In this work, we adopt the LJ potential in order to model a single-element system with a phase change between a face centered cubic (fcc) crystalline solid and a disordered liquid melt. The potential energy ϕ due to the LJ interaction is defined as:

ϕ(r)=4 a [(br)12(br)6],

where r is the distance between two interacting particles and a is the minimum potential energy of the interaction (i.e., the depth of the potential well) and b is the distance at which the potential energy is zero. Our simulations use a “force-shifted” interaction which combines the standard Lennard-Jones function (eq 15) and subtracts a linear term based on the cutoff distance Rc:

u(r)=ϕ(r) ϕ(Rc)(rRc)dϕdr |r=Rc r <Rcu(r)= 0r>Rc

where u(r) is the value of the potential energy at distance r, and ϕ(r) is the standard LJ potential at distance r. This modified relationship provides a computationally efficient approximation that allows the adoption of relatively small cutoff distances in the simulations while avoiding the effects of the discontinuity in the potential and forces at the cutoff distance that would arise by using the standard equation (15) (Toxvaerd & Dyre, 2011). We adopt a constant value for Rc = 2.5 b. In this work we express all the quantities obtained by MD simulations in reduced units rescaled on the a (energy) and b (length) parameters of the LJ potential. The a, b parameters and the particle mass (m) are set to unity, and the reduced units of time (˜t), pressure (˜P), stress (˜σ) and temperature (˜T) of the LJ system are defined as ˜t=t(aˆmb2)12 , ˜P=Pb3a , ˜σ=σb3a and ˜T=kBTa (with kB the Boltzmann constant), respectively (table 1).

The LJ substance can be stable as different physical phases depending on the P, T and density conditions. In this work, we focus on the equilibration between the fluid and the fcc solid LJ phases (e.g., Schultz & Kofke, 2018). The fcc phase is elastically anisotropic because of its crystallographic cubic symmetry (Quesnel et al., 1993). Therefore, the response to the applied strain is orientation-dependent. In all of our simulations, we assume that the crystallographic orientation of the solid phases is defined by the crystallographic directions [1 0 0] and [0 1 0] parallel to the Cartesian x1 and x2 of the simulation block, respectively (fig. 3A). Moreover, we used an integration time step Δt=0.003 in reduced time units, which ensures numerical stability and accuracy in the energy conservation in our simulations.

4.1. Solid-fluid equilibration at hydrostatic conditions

For a solid phase coexisting with a fluid of the same component, the hydrostatic equilibrium is achieved when the system, which is kept under a uniform hydrostatic pressure PH, reaches the corresponding temperature TH on the solid-fluid thermodynamic phase boundary (fig. 4A). However, the direct assessment of the thermodynamic equilibrium by MD simulations is often challenging, since the melting/crystallization process is a first order phase transition and the energy barriers may induce appreciable kinetic effects, i.e., deviations from thermodynamic equilibrium.

Figure 4
Figure 4.Schematic representation of solid-fluid equilibrium in a PT diagram. (A) The P,T points at which the solid and the fluid coexist at equilibrium under a homogeneous and hydrostatic stress state define the usual hydrostatic phase boundary in the PT space. Examples of snapshots of the atomic configuration of the LJ fluid (p1) and fcc solid (p3) phases, and of the solid-fluid coexistence at equilibrium (p2) as calculated from MD are shown. Atoms belonging to the solid and fluid phases (see appendix C) are shown in orange and blue, respectively. The rendering was performed using the OVITO software (Stukowski, 2010). (B) When the solid is stressed in a nonhydrostatic manner, the equilibrium between the solid and the fluid is achieved at a fluid pressure Pf higher than the pressure on the hydrostatic boundary at the same temperature (e.g., p2 in A and p4 in B). For a fixed stress state of the solid, the variation of the equilibrium Pf with temperature would define a nonhydrostatic phase boundary that is shifted to higher pressures (and lower temperatures) compared to the hydrostatic boundary. Therefore, a solid can be melted under nonhydrostatic stress conditions because of the shift of the phase equilibria (e.g., p2 is a coexistence state in A but it melts in B).

To obtain an initial estimate of the melting temperature at a given pressure PH, we performed nonequilibrium melting and crystallization MD simulations following the hysteresis method (Luo et al., 2004). To this aim we first create an initial configuration corresponding to a homogeneous solid LJ material by placing the atoms on the nodes of a cubic fcc crystal lattice. The unit cell is replicated 6 times along x1 and x2 and 36 times along x3, and contains a total of 5184 atoms. The boundary conditions of the MD simulations are periodic along x1, x2 and x3. The simulation was performed in a NPT ensemble (with constant properties at every timestep: N, number of particles; P, hydrostatic pressure; T, temperature). The homogenous crystal was heated incrementally at constant hydrostatic PH until complete melting, and then cooled to the initial temperature. The heating and the cooling have an equal duration and the T is varied following a linear ramp. The evolution of the energy (E) of the system as a function of T is discontinuous because of the latent heat of melting/crystallization (fig. 5). Since the kinetic barrier induces superheating and supercooling effects, the T+ of complete melting and the T of complete crystallization are not equal, and they can be determined from the hysteresis defined in the ET space (fig. 5). Following Luo et al. (2004), we determined from the hysteresis the equilibrium melting TH at hydrostatic conditions (TH = 0.606 in LJ reduced units, see fig. 5). However, such nonequilibrium approach does not provide us with a configuration of the system in which the solid and the fluid phases are coexisting in a stable equilibrium.

Figure 5
Figure 5.Potential energy-temperature hysteresis, generated by heating and cooling at constant rate in an NPT ensemble for 30106 timesteps at the hydrostatic PH=0.5. The T+=0.737 and T=0.435 are found as the T of melting and crystallization, respectively. Because of the cooling rate and the effect of boundary conditions, defects are formed during the crystallization along the cooling path. Therefore, the solid formed from supercooled liquid has higher enthalpy than the initial solid, explaining the small difference in energy between the heating and cooling paths at T<T. The equilibrium temperature of melting at hydrostatic pressure TH=0.606 is calculated as described in Luo et al. (2004, eq. (14)). All quantities are reported as LJ reduced units.

In order to obtain an equilibrated and stable system with coexisting solid-fluid, we implemented the two-phase method in which the solid and fluid phases are kept in direct contact (Dozhdikov et al., 2012; Luo et al., 2004; Yoo et al., 2009). We used the same initial configuration as in the previous simulations with a solid supercell with 5184 atoms. The extension of the simulation block along the x3 axis prevents phase boundary effects during the simulation (Noya et al., 2008). This solid system is equilibrated in an NPT ensemble at the hydrostatic pressure PH at a temperature Ts below the equilibrium melting temperature TH that was estimated previously with the hysteresis method. The atoms in the simulation block are then grouped into three layers, each containing an equal number of atoms. The interfaces separating the layers are normal to the x3 axis of the simulation block (fig. 3). The atoms of the central layer are kept “frozen”, at the temperature Ts, and therefore they remain in the solid state. The other two layers are heated up to a temperature Tm above the T+ determined from the hysteresis method. The heating is performed in a NPx3AT ensemble (with constant: N, number of particles; Px3, total stress along the x3 axis; A, area in the plane normal to x3; T, temperature) which ensures that the total stress along the x3 axis is kept close to the prescribed value PH during the entire heating (Yoo et al., 2004). The Steinhardt’s local bond-order parameters (Auer & Frenkel, 2005; Steinhardt et al., 1983) are determined from the coordinates of the atoms in order to assess for each individual atom whether it is in the solid (more ordered) or liquid state (less ordered). The details of the calculation of the order parameter are reported in appendix C. With this procedure, we can assess if the two lateral layers (fig. 3) are completely melted.

Then, the temperature of the molten layers is rapidly dropped to the level Ts of the “frozen” solid layer in a NPx3AT ensemble. At the final stage, the entire system composed of the coexisting solid and fluid layers is allowed to equilibrate. The equilibration is performed in the NPH (i.e., constant particle number, hydrostatic pressure and enthalpy) ensemble by maintaining the prescribed PH and allowing all the sides of the simulation box to be adjusted independently in order to relax the anisotropic stress that might arise because of the initial small differences of temperature and stress between the solid and the fluid layers (Morris & Song, 2002; Yoo et al., 2009). During this final equilibration in the NPH ensemble, the temperature is not constrained and evolves according to the energy and the configuration of the system. If the solid-fluid system is at T higher than the melting T, the Gibbs free energy of the crystalline phase is higher than that of the liquid phase and the crystalline phase will melt due to absorption of heat, resulting in a decrease in kinetic T of the system (Yoo et al., 2009). The opposite situation takes place if the T of the system is lower than the melting T, leading to the crystallization of the melt and to the increase in the kinetic T of the system.

The solid-fluid equilibrium is achieved when T fluctuates around a constant mean value (i.e., the temperature distribution is Gaussian) in time. The number of solid and liquid atoms (as determined with the analysis in appendix C) must also remain steady in time, which indicates that the rate of melting and crystallization are equal. On the contrary, if the initial energy was too high or too low, the coexistence of the solid and fluid phases would not be stable, and it would lead to a complete crystallization or melting of the system. In this case, the procedure is repeated by varying the initial Ts or Tm to find a final state where both the solid and fluid are present and are in mutual equilibrium. When the coexisting system has reached equilibrium, the average pressure and temperature of melting (PH, TH) are calculated as the time-averages over a production run of 5 106 timesteps.

The TH = 0.60935 (15) (LJ reduced units) of equilibrium melting determined at PH=0.5 (LJ reduced units) with the two-phase methods agrees with the estimate obtained from the hysteresis method (TH = 0.606, LJ reduced units) to better than 0.5%, as expected for an equilibrated system (Luo et al., 2004). Moreover, as an independent benchmark of the reliability of this approach we used our simulation strategy to reproduce the equilibrium melting T reported in the literature for the standard LJ formulation (eq 15). Mastny and de Pablo (2007) calculated accurately the melting T (0.7793 ± 0.0004, LJ reduced units) at a pressure of 1 (reduced units) assuming an infinite-size and full-potential. By adopting the same standard LJ formulation, and a large cutoff distance of 10b to approximate the full-potential, we obtained an equilibrium melting T of 0.7787 ± 0.0003 at a pressure of 1 (reduced units), which closely approximates the results of Mastny and de Pablo (2007) despite the finite-size of our system.

4.2. Simulations under nonhydrostatic stress at constant temperature

The variation of the pressure of the fluid at equilibrium with a stressed solid was investigated in a series of simulations in which various strain states are applied to the solid. To this aim, the configuration obtained after the previous hydrostatic equilibration at TH, PH (tables 2, 3) was used as the initial configuration. During the equilibration in the NPH ensemble discussed above, the lateral dimensions of the simulation box are adjusted independently, and they fluctuate randomly around their equilibrium value under hydrostatic pressure. Therefore, the sides lx1 and lx2 of the LJ solid (along x1 and x2 , respectively) in the final configuration may have a small difference (less than 0.03%), which results in a small virtual lateral strain of the solid. In order to have a reference configuration with a completely laterally unstrained solid, we fix lx1 = lx2 and perform a further isothermal equilibration in the canonical (NVT, constant number of particles, volume and temperature) ensemble by employing a Nosé–Hoover thermostat (Shinoda et al., 2004) at the equilibrium TH = 0.60935. The final configuration (shown in fig. 3A) still preserves the solid-fluid equilibrium at hydrostatic conditions, and the stress of the system is calculated as time-averages during a production run of 50 106 timesteps (table 2).

Table 2.Properties of the reference LJ solid-fluid system equilibrated at hydrostatic conditions in a MD simulation.
TH PH Stress of the solid layer Average ωsH Average ωfH
 σS11 σs22 σs33 ¯σs
0.60935(1) 0.4976(13) -0.5000(17) -0.5000(18) -0.4971(27) -0.4990(36) 1.067(3) 1.204(5)

The properties are calculated as averages over the simulation time. The number in parentheses for each entry indicates the 68% confidence intervals for the rightmost digits of the reported value. Stress, pressure, temperature and volume values are reported as LJ reduced units.

Table 3.Elastic properties of the solid LJ at hydrostatic pressure at the PH,TH conditions of solid-fluid equilibrium (see table 2).
C1111 C1122 C3232 D1111 D1122 D3232 K E ν
33.4(8) 19.4(5) 21.9(3) 0.0521 -0.0191 0.0456 24.0 30.5 0.2884

The components of the stiffness tensor (Cijkl, in LJ reduced stress units) are calculated with MD for a cubic fcc LJ solid with orientation [1 0 0] // x1 and [0 1 0] // x2, following the procedure outlined in Appendix E. The elastic compliance components are obtained as Dijkl= C1ijkl. The K (bulk modulus, in LJ reduced stress units), E (Young’s modulus, in LJ reduced stress units) and ν (Poisson’s ratio) are calculated from the compliance tensor. The number in parentheses for each entry indicates the 68% confidence intervals for the rightmost digits of the reported value.

To generate the nonhydrostatic stress states in the solid, the previous hydrostatic configuration was subject to lateral tensile or compressive strains parallel to the x1 and x2 coordinate axes that are normal to the “free” surfaces of the solid, by scaling the respective dimensions of the box. The x3 dimension is left unchanged since the stress in this direction will be determined by Pf. During this loading, the volume of the solid may change in response to the applied strain, and it is not constrained in the simulation (see section 5 for a discussion on the evolution of the mean stress of the solid under stress). Upon loading, the applied compressive or tensile strain determines an initial increase or decrease of the Pf, respectively. At the same time, a nonhydrostatic stress is developed in the solid phase due to the applied strain. However, such virtual configuration is not in equilibrium since the stress component of the solid normal to the solid-fluid interface (σs33 in fig. 3) may not be equal to the Pf, as would be expected from force balance. An isothermal equilibration was therefore performed in the NVT keeping the new x1, x2 and x3 dimensions of the simulation box and the TH constant, for 50 106 timesteps.

At this point we would like to point out that with MD simulations the equilibration is not found by imposing a specific thermodynamic relationship, such as that of equation (4). Instead, the equilibration is achieved by letting the trajectories of the atoms evolve in response to their potential and kinetic energy, which are in turn affected by the prescribed strain and temperature conditions. This results in the spontaneous crystallization or melting of a limited amount of material due to difference in energy between the solid and fluid phases in the strained system. The atomic volume of the LJ fcc solid is smaller than that of the liquid (table 2). During the equilibration in the NVT ensemble, the crystallization of part of the fluid or the melting of part of the solid leads to a decrease or increase of the Pf, respectively. The stress state of the solid also varies in response to its own crystallization/melting and to the variation of the Pf. If the equilibrium is restored at the prescribed strain conditions, the stress component σs33 becomes equal to the Pf and the energy configuration of the atoms in the two phases counteracts further crystallization or melting.

After the equilibration is achieved, a 50 106 timesteps production run is performed to compute the resulting equilibrium Pf and to produce the snapshots with the per-atom positions and stresses required to calculate the order parameters (appendix C) and the stress of the solid layer (calculate as described in appendix D) in the post-processing. For each simulation, the pressure of the fluid and the stress of the solid in the equilibrium stage were computed as ensemble averages and are reported in table 4.

Table 4.Lateral strains (ε11=ε22) applied to the solid and the corresponding stress state of the solid and solid-fluid equilibrium pressure (PfMD) obtained from MD simulations.
lateral strain Stress of solid layer from MD Stress of solid layer calculated from strain (assuming linear elasticity) PfMD Pftheory
 σS11 σS22 σS11=σS22
-0.00800 -0.7726(13) -0.7734(13) -0.7509 0.5122(17) 0.5150
-0.00600 -0.6971(13) -0.6966(13) -0.6865 0.5057(10) 0.5069
-0.00500 -0.6587(13) -0.6589(13) -0.6544 0.5036(12) 0.5037
-0.00400 -0.6265(14) -0.6261(14) -0.6231 0.5008(24) 0.5016
-0.00300 -0.5929(13) -0.5928(13) -0.5922 0.4998(24) 0.4998
-0.00200 -0.5575(14) -0.5577(14) -0.5611 0.4985(20) 0.4984
-0.00100 -0.5238(13) -0.5234(13) -0.5299 0.4980(11) 0.4977
0.00000 -0.5000(17) -0.5000(18) -0.5000 0.4976(13) 0.4976
0.00100 -0.4684(10) -0.4684(10) -0.4693 0.4977(10) 0.4978
0.00200 -0.4429(15) -0.4424(15) -0.4402 0.4983(12) 0.4983
0.00300 -0.4132(14) -0.4136(14) -0.4101 0.4993(11) 0.4994
0.00400 -0.3921(15) -0.3917(14) -0.3809 0.5009(15) 0.5006
0.00500 -0.3687(15) -0.3684(15) -0.3521 0.5040(9) 0.5023
0.00600 -0.3455(15) -0.3463(14) -0.3234 0.5066(16) 0.5042
0.00800 -0.3041(15) -0.3040(15) -0.2663 0.5127(11) 0.5087

The stress state of the solid expected from the applied strain assuming linear elasticity is reported. The theoretical solid-fluid equilibrium pressure Pftheory was calculated using the eq. (14) starting from the stress state computed from MD. The number in parentheses for each entry indicates the 68% confidence intervals for the rightmost digits of the reported value. Stress/pressure values are reported in LJ reduced units.

5. SOLID-FLUID EQUILIBRIUM AT NONHYDROSTATIC CONDITIONS ASSESSED BY MOLECULAR DYNAMICS SIMULATIONS

For a solid-fluid system at hydrostatic conditions, all the P,T points at which the solid and the fluid phases coexist in thermodynamic equilibrium define the usual thermodynamic phase boundary in PT space (e.g., fig. 4A). As discussed above, equation (14) describes the expected variation of the equilibrium Pf when a single-component solid-fluid system is taken from a hydrostatic equilibrium to a nonhydrostatic state at constant temperature. According to equation (14), the pressure of a fluid in equilibrium with a nonhydrostatic solid is expected to be always higher than the pressure of the same fluid equilibrated with a hydrostatic solid due to the positive contribution of deviatoric stress to the specific Helmholtz energy of a solid. As a consequence, the equilibrium phase boundary, defined by all the Pf, T points representing the pressure of the fluid in equilibrium with the stressed solid at varying temperature, would always be shifted to higher pressures (and lower temperatures) compared to the hydrostatic (classic) phase boundary (fig. 4B).

In order to confirm quantitatively the magnitude of the change in equilibrium fluid pressure and the mean stress of the solid under various nonhydrostatic stress conditions, we performed several atomistic MD simulations by adopting the procedure outlined above. Compared to the previous study by Frolov and Mishin (2010), we investigated the conditions of nonhydrostatic equilibria under large (non-zero) initial pressure. In our simulations we assume that the strain along the x1 and x2 coordinate axes are always equal one another (ε11=ε22) with no shear components. We define the resulting strain configuration as axially-symmetric (fig. 3B). As will be discussed later, such a configuration produces a volume strain in the solid phase. Because the internal pressure of the fluid phase is independent of the shape of the simulation cell, the stress state of the fluid is hydrostatic in the entire range of deformation. On the other hand, the applied axially-symmetric strain together with the chosen orientation of the solid (fig. 3) induces a nonhydrostatic stress in the solid with equal lateral stress components (σs11=σs22), which are aligned along the x1 and x2 axes of the simulation cell. Because of the previous relation, we have the following equality between the qs11 and qs22 tensor components: qs11= σs11+Pf=qs22=σs22+Pf . Finally, the mechanical equilibrium along the x3 direction imposes that σs33=Pf and, since the principal axes of strain and stress coincide with the coordinate axes, the shear components are zero (qij=0 when ij).

The results of the MD simulations (fig. 6A) confirm that the Pf at equilibrium increases for any stress state of the solid away from the initial hydrostatic reference pressure PH. Figure 6B shows that also the ¯σs varies from the initial hydrostatic pressure when a nonhydrostatic stress is applied to it. This is expected, because of the strain field used in our simulations (ε11=ε22), combined with the crystallographic symmetry and orientation of the solid, which results in a volumetric strain of the solid. At the maximum deformation explored in our simulations, the variation of ¯σs with respect to the initial hydrostatic PH ( 37%, fig. 6B) is one order of magnitude larger than the variation of the Pf ( 2.8%, fig. 6A). Moreover, when the solid is placed under nonhydrostatic conditions, its “pressure” (i.e., ¯σs) is decoupled from the pressure of the fluid Pf with which it is in thermodynamic equilibrium and can deviate largely from it (fig. 6C).

Figure 6
Figure 6.Results of MD simulations as a function of the nonhydrostatic stress applied to the solid, with components qs11=qs22 (fig. 3B, table 4) at a constant temperature T0=0.60935 (LJ reduced units). (A) Pressure of the LJ fluid at nonhydrostatic equilibrium. (B) Deviation of the mean stress (¯σs) of the nonhydrostatic LJ solid layer from its initial mean stress (PH) at hydrostatic equilibrium. (C) Deviation of ¯σs from the negative of the equilibrium pressure (Pf) of the fluid. The solid lines in (A) and (C) are calculated from equation (14) using the properties at hydrostatic conditions reported in tables 2 and 3. The ¯σs along the solid line in (B) is calculated from the strain applied to the solid and the Pf, using linear elasticity (elastic properties in table 3). All the stress and pressure quantities are reported as LJ reduced units.

The solid line shown in figure 6A represents the Pf calculated from equation (14) using the elastic constants determined on our solid at the hydrostatic reference state (tables 2 and 3). A very good agreement is found between the theoretical solution (eq. 14) and the results of our MD simulations under small deformations (fig. 6A). The deviations observed under large deformations are due to the assumption of linear elasticity in the derivation equation (14), which assumes that the deformations are small and neglects the variation of the elastic properties of the solid with stress. Therefore, the analytical solution (eq. 14) does not account for the stiffening of the solid under compressive deformations, nor for its softening under tensile stress. As a consequence, it underestimates the value of stress (and therefore of qs11) in compression, while it overestimates it in tension. This is also reflected in Fig. 6B, which shows that by using constant elastic properties (solid line), the variation of ¯σs is underestimated under compression and overestimated under tension.

To further investigate this effect, for each simulation we used linear elasticity to recalculate the “linearized” stress of the solid, by knowing the components ε11 and ε22 of the applied strain and imposing that, because of mechanical equilibrium, σs33=Pf(fig. 3). Figure 6A shows an excellent agreement of the recalculated points with the predictions of the analytical solution (blue dots in fig. 6A), confirming that the small deviations with respect to our MD results are due to the assumptions of small strains behind equation (14). However, the errors of the analytical solution on the determination of Pf are less than 1% even at the largest deformation, despite the very large differential stress of the solid (50% of the initial hydrostatic pressure) at these conditions.

We note that not all of the possible strain configurations lead to the same evolution of Pf, and to the same deviation of ¯σs from PH and Pf, as those shown in figure 6. Figure 7 shows these equilibrium quantities as obtained from the theoretical relationship of equation (14), by varying q11 and q22 independently. In our simulations we explored the path with q11=q22 (solid red line in fig. 7A, B, C). This path produces the largest deviation of the ¯σs from the reference pressure PH (fig. 7B) and from Pf (fig. 7C). Since ¯σs=q11+q22 3Pf3, the ¯σs remains equal to Pf only when q11= q22 (red dashed line in fig. 7C). Because ¯σs=Pf along this particular path, the tensor qsij= σsij+δijPf is equal to the deviatoric stress of the solid. With our choice of crystallographic orientation of the solid, and since the LJ fcc solid has a cubic crystallographic symmetry, ¯σs stays equal to Pfwhen the applied strain has components ε11= ε22. However, even this strain configuration induces a volume strain of the solid phase with respect to the initial hydrostatic reference state at PH, and therefore a variation of its ¯σs (red dashed line, fig. 7B). This is due to the fact that, despite ε11= ε22, the normal strain component along the third direction (ε33) is controlled by the equilibrium Pf, which varies in response to the strain applied to the solid (fig. 7A). A more complex path in q11 and q22 can be found that maintains the ¯σs equal to the initial PH (black dash-dotted line, fig. 7B). However, along this path the ¯σs does not equal the equilibrium Pf in the stressed system (black dash-dotted line, fig. 7C).

Figure 7
Figure 7.Equilibrium quantities of the LJ system under nonhydrostatic stress, as obtained from the theoretical relationship of equation (14) by varying the components q11 and q22 independently. (A) Variation of the equilibrium pressure of the fluid Pf with stress (grey surface). (B) Deviation of the mean stress (¯σs) of the nonhydrostatic solid layer from its initial mean stress (PH) at hydrostatic equilibrium, shown as contours. (C) Deviation of ¯σs from the negative of the pressure of the fluid (Pf) at nonhydrostatic equilibrium, shown as contours. Three deformation paths are shown as lines. Red solid line: path with q11=q22, as explored in our MD simulations (see fig. 6). Red dashed line: path with q11=q22; along this path the ¯σs = Pf (C), but deviates from the initial PH (B). Black dashed line: path along which the ¯σs = PH (B), but deviates from Pf (C). All the shown stress and pressure quantities are reported as LJ reduced units. The properties used in the calculations are reported in tables 2 and 3.

6. DISCUSSION

6.1. Equilibrium pressure of mineral phases in equilibrium with their melt under nonhydrostatic conditions

The results of our MD simulations reported in fig. 6A show that the theoretical relationship based on linear elasticity (eq 14) describes correctly the change of the fluid-solid equilibrium pressure up to compressions and tensions with a qs11=qs22 of at least ± 20% of the initial pressure at hydrostatic conditions. Therefore, we applied this theory to calculate the equilibrium pressure between several solid mineral phases (Mg2SiO4, forsterite; SiO2, quartz; NaCl, halite) and their pure melt. Initially the mineral-melt equilibrium is evaluated at hydrostatic conditions. Since each of these systems reaches the solid-melt equilibrium at different P and T conditions, in order to compare consistently their behavior, we introduce a nondimensional pressure defined as the pressure PH at which a solid phase reaches the hydrostatic equilibrium with its pure melt at a given TH, normalized over the bulk modulus (K) of the solid phase at that PH,TH. For the LJ system at the hydrostatic conditions PH,TH of our simulations, the nondimensional hydrostatic equilibrium pressure is P=PHK=0.0253 (table 3). We used the Perple_X software (Connolly, 2005, 2009) with data from the Holland and Powell (2011) thermodynamic database to map the variation in the P-T space of the nondimensional solid pressure P=PK for each mineral phase. Then, in order to compare the solid-melt equilibrium pressure of all the systems, we found for each of them the P and T condition along their melting line where the condition P= 0.0253 is satisfied. The P and T of hydrostatic solid-melt equilibrium obtained in this way for each phase are reported in table 5, and they are used as reference states to evaluate the perturbation due to nonhydrostatic stress.

Table 5.Properties of the solid minerals at solid-melt hydrostatic equilibrium calculated with the Perple_X software (Connolly, 2005, 2009) with data from the Holland and Powell (2011) thermodynamic database.
Phase P (GPa) T (K) Poisson's ratio Bulk modulus (GPa) Young modulus (GPa) Solid molar volume (cm3/mol) Melt molar volume (cm3/mol)
NaCl 0.95 1264 0.163 37.5 75.9 30.08 33.91
SiO2 1.49 2262 0.217 58.9 100.0 23.22 25.16
Mg2SiO4 2.70 2304 0.270 106.7 146.9 46.02 48.60

For all the phases the properties are calculated at the non-dimensional pressure P=P/K = 0.00253 as explained in the main text. The corresponding dimensional pressures and temperatures on the respective melt lines are reported.

Since the anisotropic properties of these mineral phases are not available at the required P and T conditions, we used their isotropic properties calculated from the Perple_X software (Young’s modulus and Poisson’s ratio, table 5) together with equation (B2) to investigate the variation of the equilibrium pressure of the melt expected from the application of a nonhydrostatic stress to the solid mineral. We applied an axially-symmetric strain state as in the MD simulations (i.e., ε11=ε22, and no shear components) assuming the same orientation of the Cartesian axes as in figure 3. Since the minerals are assumed to be elastically isotropic, the exposition of the solid layer to these lateral strains along the x1and x2 directions and to the Pf along the x3 direction, results in a nonhydrostatic stress in the solid with qs11=qs22 (i.e., the same path as the solid red line in fig. 7).

Figure 8 shows the change in the equilibrium Pf and in the ¯σs for a nonhydrostatic compression and tension up to ± 20% of the initial hydrostatic pressure. In order to compare the results obtained for minerals at solid-melt equilibrium at different P and T, for each system the applied nonhydrostatic stress and the variations in Pf(fig. 8A) and ¯σs (fig. 8B, C) are all normalized to the initial PH of the system.

Figure 8
Figure 8.Melt-solid equilibrium as a function of the nonhydrostatic stress applied to the solid with components qs11=qs22 at constant T. The initial hydrostatic PH, TH conditions of Mg2SiO4 (2.70 GPa,2304 K), SiO2 (1.49 GPa,2262 K), NaCl (0.95 GPa,1264 K) and the LJ material (0.5ε/σ3,0.609ε/kB) are determined from their respective melt lines at the nondimensional pressure P=P/K=0.0253, as explained in the main text. All of the quantities are normalized on the initial hydrostatic PH of the system. (A) Variation of the melt equilibrium pressure calculated from equation (B2) assuming isotropic elastic properties for the solid. The dashed line (anisotropic LJ) has been calculated with equation (14) assuming anisotropic elastic properties for the solid. (B) Deviation of the mean stress (¯σs) of the nonhydrostatic solid layer from its initial mean stress (PH) at hydrostatic equilibrium. (C) Deviation of ¯σs from the negative of the equilibrium pressure (Pf) of the fluid. The properties of each phase at hydrostatic conditions are reported in tables 2, 3 and 5.

Figure 8A shows that for all the investigated minerals, the equilibrium pressure of the melt increases for any compressional or tensile stress state of the solid. In general, the equilibrium melt pressure increases more in systems in which the solid phase is stiffer. Indeed, solids that are elastically stiffer store more elastic energy for the same amount of applied strain (eq 8). Therefore, in order to preserve the solid-melt equilibrium in systems with a stiffer solid phase, the pressure of the melt must also increase for the chemical potential of the melt (μf) to keep up with the increase in Helmholtz energy of the solid (see eq 10). Moreover, figure 8A shows a difference in the change of the equilibrium Pf for the anisotropic and isotropic LJ, calculated respectively with equations (14) and (B2). The discrepancy between the isotropic and anisotropic curves can be estimated from equation (B3), which shows that for the LJ material the change in the Pf is underestimated by 29% when using isotropic elastic properties. This is due to the LJ solid being highly anisotropic, with a universal anisotropic index (Ranganathan & Ostoja-Starzewski, 2008) AU=1.73 (for reference, an isotropic solid would have AU= 0). The anisotropic calculation cannot accurately be performed for the other phases, since their anisotropic elastic constants at the chosen melting conditions are not available. However, we can roughly estimate the discrepancy between the anisotropic and the isotropic solution, by evaluating equation (B3) at ambient conditions (1 bar, 298K) and assuming that this discrepancy stays constant at increased P and T conditions.

For example, for quartz (anisotropic index AU=1.30 at room conditions), the isotropic approximation underestimates the equilibrium pressure by 27% when using the anisotropic elastic properties of quartz from ambient conditions (cf. Lakshtanov et al., 2007). However, even allowing for the contribution of elastic anisotropy, the change in the normalized melt-solid equilibrium pressure (ΔPf/PH) would still be at least one order of magnitude less than the applied normalized nonhydrostatic stress (qs11/PH). On the other hand, the normalized variation of the mean stress of the solid (¯σs+PH)/PH; fig. 8B) is of the same order of magnitude as the applied normalized nonhydrostatic stress (qs11/PH). Under this deformation, the ¯σs also deviates substantially from the Pf of the fluid to which it is in equilibrium with (fig. 8C).

6.2. Implications for solid-fluid multiphase systems in porous rocks

As shown in figure 2 and discussed above, pressure gradients are expected in multiphase systems when solid minerals are in direct contact with fluids, even under static conditions and without external tectonic stresses. While fluid will always experience a hydrostatic stress gradient (fig. 2A), solid grains will experience nonhomogeneous and nonhydrostatic stress (fig. 2B, C) due to the changes in the contact area between the grains and the requirements imposed by mechanical equilibrium. Therefore, nonhydrostatically stressed minerals and hydrostatic fluids can coexist, a condition that breaks the fundamental assumption behind the use of Gibbs potentials to describe reactions in such solid-fluid systems. With MD simulations we can quantify the magnitude of the errors expected if a hydrostatic Gibbs potential is used to describe phase equilibria between a fluid and a solid for nonhydrostatic conditions. Our results are in agreement with previous theory (Frolov & Mishin, 2010; Gibbs, 1876; Sekerka & Cahn, 2004), and confirm that for single component systems, the equilibrium pressure of the fluid changes quadratically with the nonhydrostatic stress applied to the solid matrix (fig. 6).

However, for the nonhydrostatic stresses expected in the lithosphere (with differential stresses of several hundreds of MPa; Andersen et al., 2008; Brace & Kohlstedt, 1980; Kiss et al., 2023; Prieto et al., 2012), the variation of the equilibrium pressure of the fluid, and therefore the displacement of the thermodynamic pressure of phase equilibria, is almost negligible. For example, we can consider the worst case among those reported in figure 8A of a forsterite in hydrostatic equilibrium with its pure melt at a pressure of 2.7 GPa and temperature of 2304 K (table 5). For an applied axially-symmetric nonhydrostatic stress with qs11=qs22= ± 0.5 GPa, the equilibrium pressure of the melt increases by 0.02 GPa, i.e., an increase of less than 1% with respect to the initial pressure in the purely hydrostatic system. This implies that, for cases similar to ours, the pressure of a melt that is in thermodynamic equilibrium with a nonhydrostatically stressed solid is close to melt pressure at the hydrostatic limit. Therefore, the thermodynamic pressure that controls the solid-melt equilibria is equal to the pressure of the fluid, which remains almost constant even when a differential stress is applied to the solid matrix. This has relevant implications for models that describe multiphase systems under deformation, since, as a first approximation, they can use the pressure of the melt in a standard hydrostatic thermodynamic model (e.g., Perple_X software; Connolly, 2005, 2009) in order to estimate the phase equilibria of the stressed system.

However, it is important to note that when the system is under nonhydrostatic stress, even at thermodynamic equilibrium, the pressure of the melt Pf may not be equal to the “pressure” of the solid (i.e. -¯σs), and their deviation is often non negligible (fig. 6C, 7C, 8C). Therefore, in a stressed system, the use of the solid “pressure” instead of the fluid pressure to calculate the equilibrium assemblages may lead to large errors. As an example, we use once again the Mg2SiO4 system with an initial hydrostatic solid-melt equilibrium at 2.70 GPa and 2304 K. When an axially-symmetric nonhydrostatic stress qs11=qs22=0.5 GPa is applied, the Pf at thermodynamic equilibrium becomes 2.72 GPa, i.e., close to its value at the initial hydrostatic equilibrium. However, the solid “pressure” reaches 3.04 GPa (i.e., +13% variation) in the same deformation, and diverges significantly from the Pf.

The deviation of the ¯σs from the Pf with increasing applied stress in porous systems (fig. 2A and B), has important implications for understanding the interactions between fluids and rocks under deformation. In current hydro-mechanical-chemical (HMC models (e.g., Schmalholz et al., 2020), the density and the energy of the solid matrix are approximated using the values of the Pf. This is justified when the difference between the ¯σs and the Pf is small. However, since nonhydrostatic stress can be developed in the solid matrix in such systems, the use of Pf as a proxy of ¯σs can lead to inaccurate estimates of the density of the solid, even when the solid and the fluid are in thermodynamic equilibrium. The magnitude of such errors must be assessed on a case-by-case basis, because they depend on the specific deformation applied to the system (fig. 7C).

7. CONCLUSIONS

The interpretation of phase equilibria and reactions in geological materials has long been approximated by standard thermodynamics which assumes that the stress in the systems is hydrostatic and homogeneous (i.e., the same for all the phases involved). However, stress gradients and nonhydrostatic stresses are common in rocks and can occur even in porous systems with coexisting solid minerals and fluids (e.g., fig. 2). Currently there is still no accepted theory that extends thermodynamics to include the general effects of nonhydrostatic stress on reactions, and the use of several thermodynamic potentials in stressed geological system is still debated.

By means of MD simulations, we have investigated the effect of a homogeneous nonhydrostatic stress on the equilibrium between a solid and its pure melt. The advantages of this approach are: (i) the stress of the solid phase can be kept homogeneous. Therefore, the direct effect of a homogeneous, nonydrostatic stress on thermodynamic equilibria can be assessed, and the indirect perturbations due to local stress concentrations, that may be present in experiments, are avoided. (ii) It can be used to explore the equilibration of fluid-solid systems under known and controlled strain/stress states. (iii) The energy of the system and the local stress are monitored simultaneously while the system evolves toward equilibrium. (iv) It intrinsically accounts for the elastic anisotropy of the solid phase. (v) The equilibrium condition can be defined without any a priori assumptions of a specific thermodynamic theory. Such an approach is therefore extremely useful to validate the applicability of existing thermodynamic theories.

Our MD simulations show that for small applied strains, the equilibrium pressure of the fluid increases quadratically with the nonhydrostatic stress generated in the solid. These findings agree with the existing thermodynamic theory for single component systems (Frolov & Mishin, 2010; Gibbs, 1876; Sekerka & Cahn, 2004) and extend previous MD simulations performed at zero initial pressure (Frolov & Mishin, 2010).

We then use the analytical relationship derived by Sekerka and Cahn (2004) to explore the equilibrium between common mineral phases and their pure melts at the differential stresses expected in the lithosphere. We observe that, by applying an axially-symmetric load, the Pfat equilibrium increases by less than 1% even when the solid matrix is subject to a differential stress that is 20% above (compression) or below (tension) the initial hydrostatic pressure. However, the variation of the ¯σs is one order of magnitude larger than the variation of Pf under these stress conditions (fig. 8). In general, a decoupling of the “pressure” of the solid (¯σs) from the pressure of the fluid Pf at the thermodynamic equilibrium is expected whenever the solid is put under a nonhydrostatic stress (fig. 7C). This has important implications for hydro-mechanical-chemical (HMC) models that describe the interactions between fluids and minerals under deformation. Our results imply that, at least for simple single-component porous systems, phase equilibria can be accurately predicted by taking the Pf as a proxy of the equilibration pressure and assuming that it is not perturbed significantly by the typical levels of differential stress expected in the solid matrix in such systems (fig. 2). Despite the fact that the ¯σs can deviate substantially from the Pf in stressed systems, ¯σs still remains a valuable thermodynamic quantity that is required for the accurate calculation of density and energy of the solid phase.

Our preliminary MD simulations at high pressure and temperature conditions show that this approach can provide quantitative estimates on the direct effect of nonhydrostatic stresses on phase equilibria and reactions in multiphase systems at geologically relevant conditions. While experimental studies of multicomponent systems suggest a significant role of the Pf (e.g., Llana-Fúnez et al., 2012), a complete thermodynamic description of such systems under stress is still lacking. In this respect, our framework may be extended to real mineral systems by adopting appropriate force fields for geological materials. While our analysis does not take into account the effect of stress heterogeneities in the sample at larger scales, which can indirectly affect phase equilibria, the knowledge of the direct effect of deviatoric stress on thermodynamic equilibria gained from MD simulations can be incorporated into chemical-mechanical models that describe macroscopic geological systems under deformation.


ACKNOWLEDGEMENTS

All the authors acknowledge the Mainz Institute of Multiscale modelling (M3ODEL) for providing funding at the early parts of this research. M.L.M. was supported by an Alexander von Humboldt research fellowship. Boris Kaus received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (ERC CoG grant agreement No 771143 – MAGMA). We thank M. Alvaro, R.J. Angel, M. Prencipe and two anonymous reviewers for constructive comments that allowed us to improve this manuscript. The editorial handling and insightful comments by M. Brandon and A. Smye are greatly appreciated.

AUTHOR CONTRIBUTIONS

All authors contributed to the development of the ideas and writing of the manuscript. M.L.M and T.S. developed the strategy for the molecular dynamics simulations. M.L.M performed the molecular dynamics simulations and the post-processing. E.M. provided the calculation shown in figure 2. M.L.M, E.M. and B.K. derived the geological implications.

DATA AVAILIBILITY

All the sources of the data used are explicitly mentioned in the text.

Editor: Mark T. Brandon, Associate Editor: Andrew Smye

Accepted: January 24, 2024 EDT

References

Andersen, T. B., Mair, K., Austrheim, H., Podladchikov, Y. Y., & Vrijmoed, J. C. (2008). Stress release in exhumed intermediate and deep earthquakes determined from ultramafic pseudotachylyte. Geology, 36(12), 995–998. https://doi.org/10.1130/g25230a.1
Google Scholar
Anderson, O. L. (1995). Equations of state of solids for geophysics and ceramic science. Oxford University Press. https://doi.org/10.1093/oso/9780195056068.001.0001
Google Scholar
Angel, R. J., Gilio, M., Mazzucchelli, M., & Alvaro, M. (2022). Garnet EoS: A critical review and synthesis. Contributions to Mineralogy and Petrology, 177(5), 54. https://doi.org/10.1007/s00410-022-01918-5
Google Scholar
Auer, S., & Frenkel, D. (2005). Numerical Simulation of Crystal Nucleation in Colloids. In C. Holm & K. Kremer (Eds.), Advanced Computer Simulation: Approaches for Soft Matter Sciences I (pp. 149–208). Springer. https://doi.org/10.1007/b99429
Google Scholar
Augusti, G., Martin, J. B., & Prager, W. (1969). On the decomposition of stress and strain tensors into spherical and deviatoric parts. Proceedings of the National Academy of Sciences, 63(2), 239–241. https://doi.org/10.1073/pnas.63.2.239
Google ScholarPubMed CentralPubMed
Austrheim, H. (1987). Eclogitization of lower crustal granulites by fluid migration through shear zones. Earth and Planetary Science Letters, 81(2–3), 221–232. https://doi.org/10.1016/0012-821x(87)90158-0
Google Scholar
Bokeloh, J., Rozas, R. E., Horbach, J., & Wilde, G. (2011). Nucleation Barriers for the Liquid-To-Crystal Transition in Ni: Experiment and Simulation. Physical Review Letters, 107(14), 145701. https://doi.org/10.1103/physrevlett.107.145701
Google Scholar
Brace, W. F., & Kohlstedt, D. L. (1980). Limits on lithospheric stress imposed by laboratory experiments. Journal of Geophysical Research: Solid Earth, 85(B11), 6248–6252. https://doi.org/10.1029/jb085ib11p06248
Google Scholar
Brantut, N., Stefanou, I., & Sulem, J. (2017). Dehydration-induced instabilities at intermediate depths in subduction zones. Journal of Geophysical Research: Solid Earth, 122(8), 6087–6107. https://doi.org/10.1002/2017jb014357
Google Scholar
Bucher, K., & Frey, M. (2002). Petrogenesis of Metamorphic Rocks. Springer Science & Business Media. https://doi.org/10.1007/978-3-662-04914-3
Google Scholar
Campomenosi, N., Mazzucchelli, M. L., Mihailova, B., Scambelluri, M., Angel, R. J., Nestola, F., Reali, A., & Alvaro, M. (2018). How geometry and anisotropy affect residual strain in host-inclusion systems: Coupling experimental and numerical approaches. American Mineralogist, 103(12), 2032–2035. https://doi.org/10.2138/am-2018-6700ccby
Google Scholar
Chorin, A. J. (1997). A Numerical Method for Solving Incompressible Viscous Flow Problems. Journal of Computational Physics, 135(2), 118–125. https://doi.org/10.1006/jcph.1997.5716
Google Scholar
Cionoiu, S., Moulas, E., Stünitz, H., & Tajčmanová, L. (2022). Locally Resolved Stress-State in Samples During Experimental Deformation: Insights Into the Effect of Stress on Mineral Reactions. Journal of Geophysical Research: Solid Earth, 127(8), e2022JB024814. https://doi.org/10.1029/2022jb024814
Google Scholar
Cionoiu, S., Moulas, E., & Tajčmanová, L. (2019). Impact of interseismic deformation on phase transformations and rock properties in subduction zones. Scientific Reports, 9(1). https://doi.org/10.1038/s41598-019-56130-6
Google ScholarPubMed CentralPubMed
Connolly, J. A. D. (2005). Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation. Earth and Planetary Science Letters, 236(1–2), 524–541. https://doi.org/10.1016/j.epsl.2005.04.033
Google Scholar
Connolly, J. A. D. (2009). The geodynamic equation of state: What and how. Geochemistry, Geophysics, Geosystems, 10(10). https://doi.org/10.1029/2009gc002540
Google Scholar
Connolly, J. A. D. (2017). A primer in gibbs energy minimization for geophysicists. Petrology, 25(5), 526–534. https://doi.org/10.1134/s0869591117050034
Google Scholar
Connolly, J. A. D., & Podladchikov, Y. Y. (2004). Fluid flow in compressive tectonic settings: Implications for midcrustal seismic reflectors and downward fluid migration. Journal of Geophysical Research: Solid Earth, 109(B4). https://doi.org/10.1029/2003jb002822
Google Scholar
Cormier, J., Rickman, J. M., & Delph, T. J. (2001). Stress calculation in atomistic simulations of perfect and imperfect solids. Journal of Applied Physics, 89(1), 99–104. https://doi.org/10.1063/1.1328406
Google Scholar
Dahlen, F. A. (1992). Metamorphism of nonhydrostatically stressed rocks. American Journal of Science, 292(3), 184–198. https://doi.org/10.2475/ajs.292.3.184
Google Scholar
de Capitani, C., & Petrakakis, K. (2010). The computation of equilibrium assemblage diagrams with Theriak/Domino software. American Mineralogist, 95(7), 1006–1016. https://doi.org/10.2138/am.2010.3354
Google Scholar
Dozhdikov, V. S., Basharin, A. Y., & Levashov, P. R. (2012). Two-phase simulation of the crystalline silicon melting line at pressures from –1 to 3 GPa. The Journal of Chemical Physics, 137(5), 054502. https://doi.org/10.1063/1.4739085
Google Scholar
Duretz, T., Räss, L., Podladchikov, Y., & Schmalholz, S. (2019). Resolving thermomechanical coupling in two and three dimensions: Spontaneous strain localization owing to shear heating. Geophysical Journal International, 216(1), 365–379. https://doi.org/10.1093/gji/ggy434
Google Scholar
Erfani, H., Joekar-Niasar, V., & Farajzadeh, R. (2019). Impact of Microheterogeneity on Upscaling Reactive Transport in Geothermal Energy. ACS Earth and Space Chemistry, 3(9), 2045–2057. https://doi.org/10.1021/acsearthspacechem.9b00056
Google Scholar
Fenley, A. T., Muddana, H. S., & Gilson, M. K. (2014). Calculation and Visualization of Atomistic Mechanical Stresses in Nanomaterials and Biomolecules. PLOS ONE, 9(12), e113119. https://doi.org/10.1371/journal.pone.0113119
Google ScholarPubMed CentralPubMed
Ferrand, T. P., Hilairet, N., Incel, S., Deldicque, D., Labrousse, L., Gasc, J., Renner, J., Wang, Y., Green, H. W., II, & Schubnel, A. (2017). Dehydration-driven stress transfer triggers intermediate-depth earthquakes. Nature Communications, 8(1), 15247. https://doi.org/10.1038/ncomms15247
Google ScholarPubMed CentralPubMed
Frolov, T., & Mishin, Y. (2010). Effect of nonhydrostatic stresses on solid-fluid equilibrium. I. Bulk thermodynamics. Physical Review B - Condensed Matter and Materials Physics, 82(17), 1–14. https://doi.org/10.1103/physrevb.82.174113
Google Scholar
Gallagher, J. J., Jr, Friedman, M., Handin, J., & Sowers, G. M. (1974). Experimental studies relating to microfracture in sandstone. Tectonophysics, 21(3), 203–247. https://doi.org/10.1016/0040-1951(74)90053-5
Google Scholar
Gibbs, J. W. (1876). On the Equilibrium of Heterogeneous Substances—Part II. Transactions of the Connecticut Academy of Arts and Sciences, 3, 108–248. https://doi.org/10.5479/sil.421748.39088007099781
Google Scholar
Gillet, P., Ingrin, J., & Chopin, C. (1984). Coesite in subducted continental crust: P-T history deduced from an elastic model. Earth and Planetary Science Letters, 70(2), 426–436. https://doi.org/10.1016/0012-821x(84)90026-8
Google Scholar
Gonzalez, J. P., Mazzucchelli, M. L., Angel, R. J., & Alvaro, M. (2021). Elastic Geobarometry for Anisotropic Inclusions in Anisotropic Host Minerals: Quartz-in-Zircon. Journal of Geophysical Research: Solid Earth, 126(6), e2021JB022080. https://doi.org/10.1029/2021jb022080
Google Scholar
Goodier, J. N. (1934). XLVI.An analogy between the slow motions of a viscous fluid in two dimensions, and systems of plane stress. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 17(113), 554–576. https://doi.org/10.1080/14786443409462415
Google Scholar
Grüneisen, E. (1926). Zustand des festen Körpers. In C. Drucker (Ed.), Thermische Eigenschaften der Stoffe (pp. 1–59). Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-642-99531-6_1
Google Scholar
Hess, B. L., & Ague, J. J. (2021). Quantifying the Effects of Non-Hydrostatic Stress on Single-Component Polymorphs. Journal of Geophysical Research: Solid Earth, 126(5), 2020021594. https://doi.org/10.1029/2020jb021594
Google Scholar
Heuzé, O. (2012). General form of the Mie–Grüneisen equation of state. Comptes Rendus Mécanique, 340(10), 679–687. https://doi.org/10.1016/j.crme.2012.10.044
Google Scholar
Hirth, G., & Kohlstedt, D. L. (1995). Experimental constraints on the dynamics of the partially molten upper mantle: Deformation in the diffusion creep regime. Journal of Geophysical Research: Solid Earth, 100(B2), 1981–2001. https://doi.org/10.1029/94jb02128
Google Scholar
Hirth, G., & Tullis, J. (1994). The brittle-plastic transition in experimentally deformed quartz aggregates. Journal of Geophysical Research: Solid Earth, 99(B6), 11731–11747. https://doi.org/10.1029/93jb02873
Google Scholar
Hobbs, B. E. (1968). Recrystallization of single crystals of quartz. Tectonophysics, 6(5), 353–401. https://doi.org/10.1016/0040-1951(68)90056-5
Google Scholar
Hobbs, B. E., & Ord, A. (2015). Dramatic effects of stress on metamorphic reactions. Geology, 43(11), e372. https://doi.org/10.1130/g37070c.1
Google Scholar
Hobbs, B. E., & Ord, A. (2016). Does non-hydrostatic stress influence the equilibrium of metamorphic reactions? Earth-Science Reviews, 163, 190–233. https://doi.org/10.1016/j.earscirev.2016.08.013
Google Scholar
Holland, T. J. B., & Powell, R. (2011). An improved and extended internally consistent thermodynamic dataset for phases of petrological interest, involving a new equation of state for solids: THERMODYNAMIC DATASET FOR PHASES OF PETROLOGICAL INTEREST. Journal of Metamorphic Geology, 29(3), 333–383. https://doi.org/10.1111/j.1525-1314.2010.00923.x
Google Scholar
Ji, S., & Wang, Q. (2011). Interfacial friction-induced pressure and implications for the formation and preservation of intergranular coesite in metamorphic rocks. Journal of Structural Geology, 33(2), 107–113. https://doi.org/10.1016/j.jsg.2010.11.013
Google Scholar
Jolivet, L., Raimbourg, H., Labrousse, L., Avigad, D., Leroy, Y., Austrheim, H., & Andersen, T. B. (2005). Softening trigerred by eclogitization, the first step toward exhumation during continental subduction. Earth and Planetary Science Letters, 237(3–4), 532–547. https://doi.org/10.1016/j.epsl.2005.06.047
Google Scholar
Jorgensen, W. L., Maxwell, D. S., & Tirado-Rives, J. (1996). Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. Journal of the American Chemical Society, 118(45), 11225–11236. https://doi.org/10.1021/ja9621760
Google Scholar
Kaatz, L., Zertani, S., Moulas, E., John, T., Labrousse, L., Schmalholz, S. M., & Andersen, T. B. (2021). Widening of Hydrous Shear Zones During Incipient Eclogitization of Metastable Dry and Rigid Lower Crust—Holsnøy, Western Norway. Tectonics, 40(3), e2020TC006572. https://doi.org/10.1029/2020tc006572
Google Scholar
Kelemen, P. B., & Hirth, G. (2012). Reaction-driven cracking during retrograde metamorphism: Olivine hydration and carbonation. Earth and Planetary Science Letters, 345–348, 81–89. https://doi.org/10.1016/j.epsl.2012.06.018
Google Scholar
Kiss, D., Moulas, E., Kaus, B. J. P., & Spang, A. (2023). Decompression and Fracturing Caused by Magmatically Induced Thermal Stresses. Journal of Geophysical Research: Solid Earth, 128(3), e2022JB025341. https://doi.org/10.1029/2022jb025341
Google Scholar
Kohn, M. J., Mazzucchelli, M. L., & Alvaro, M. (2023). Elastic Thermobarometry. Annual Review of Earth and Planetary Sciences, 51(1), 331–366. https://doi.org/10.1146/annurev-earth-031621-112720
Google Scholar
Lakshtanov, D. L., Sinogeikin, S. V., & Bass, J. D. (2007). High-temperature phase transitions and elasticity of silica polymorphs. Physics and Chemistry of Minerals, 34(1), 11–22. https://doi.org/10.1007/s00269-006-0113-y
Google Scholar
Lechner, W., & Dellago, C. (2008). Accurate determination of crystal structures based on averaged local bond order parameters. The Journal of Chemical Physics, 129(11), 114707. https://doi.org/10.1063/1.2977970
Google Scholar
Lennard-Jones, J. E. (1924a). On the determination of molecular fields. —II. From the equation of state of a gas. Proceedings of the Royal Society of London. Series A, 106(738), 463–477. https://doi.org/10.1098/rspa.1924.0082
Google Scholar
Lennard-Jones, J. E. (1924b). On the determination of molecular fields.—I. From the variation of the viscosity of a gas with temperature. Proceedings of the Royal Society of London. Series A, 106(738), 441–462. https://doi.org/10.1098/rspa.1924.0081
Google Scholar
Llana-Fúnez, S., Wheeler, J., & Faulkner, D. R. (2012). Metamorphic reaction rate controlled by fluid pressure not confining pressure: Implications of dehydration experiments with gypsum. Contributions to Mineralogy and Petrology, 164(1), 69–79. https://doi.org/10.1007/s00410-012-0726-8
Google Scholar
Luo, S.-N., Strachan, A., & Swift, D. C. (2004). Nonequilibrium melting and crystallization of a model Lennard-Jones system. The Journal of Chemical Physics, 120(24), 11640–11649. https://doi.org/10.1063/1.1755655
Google Scholar
Mastny, E. A., & de Pablo, J. J. (2007). Melting line of the Lennard-Jones system, infinite size, and full potential. The Journal of Chemical Physics, 127(10), 104504. https://doi.org/10.1063/1.2753149
Google Scholar
Menon, S., Leines, G. D., & Rogal, J. (2019). pyscal: A python module for structural analysis of atomic environments. Journal of Open Source Software, 4(43), 1824. https://doi.org/10.21105/joss.01824
Google Scholar
Mie, G. (1903). Zur kinetischen Theorie der einatomigen Körper. Annalen Der Physik, 316(8), 657–697. https://doi.org/10.1002/andp.19033160802
Google Scholar
Morris, J. R., & Song, X. (2002). The melting lines of model systems calculated from coexistence simulations. The Journal of Chemical Physics, 116(21), 9352–9358. https://doi.org/10.1063/1.1474581
Google Scholar
Moulas, E., Burg, J.-P., & Podladchikov, Y. (2014). Stress field associated with elliptical inclusions in a deforming matrix: Mathematical model and implications for tectonic overpressure in the lithosphere. Tectonophysics, 631, 37–49. https://doi.org/10.1016/j.tecto.2014.05.004
Google Scholar
Moulas, E., Kaus, B., & Jamtveit, B. (2022). Dynamic pressure variations in the lower crust caused by localized fluid-induced weakening. Communications Earth & Environment, 3(1). https://doi.org/10.1038/s43247-022-00478-7
Google Scholar
Moulas, E., Kostopoulos, D., Podladchikov, Y., Chatzitheodoridis, E., Schenker, F. L., Zingerman, K. M., Pomonis, P., & Tajčmanová, L. (2020). Calculating pressure with elastic geobarometry: A comparison of different elastic solutions with application to a calc-silicate gneiss from the Rhodope Metamorphic Province. Lithos, 378–379, 105803. https://doi.org/10.1016/j.lithos.2020.105803
Google Scholar
Moulas, E., Podladchikov, Y. Y., Aranovich, L. Ya., & Kostopoulos, D. (2013). The problem of depth in geology: When pressure does not translate into depth. Petrology, 21(6), 527–538. https://doi.org/10.1134/s0869591113060052
Google Scholar
Moulas, E., Podladchikov, Y., Zingerman, K., Vershinin, A., & Levin, V. (2023). Large-strain Elastic and Elasto-Plastic Formulations for Host-Inclusion Systems and Their Applications in Thermobarometry and Geodynamics. American Journal of Science, 323. https://doi.org/10.2475/001c.68195
Google Scholar
Moulas, E., Schmalholz, S. M., Podladchikov, Y., Tajčmanová, L., Kostopoulos, D., & Baumgartner, L. (2019). Relation between mean stress, thermodynamic, and lithostatic pressure. Journal of Metamorphic Geology, 37(1), 1–14. https://doi.org/10.1111/jmg.12446
Google Scholar
Noya, E. G., Vega, C., & de Miguel, E. (2008). Determination of the melting point of hard spheres from direct coexistence simulation methods. The Journal of Chemical Physics, 128(15), 154507. https://doi.org/10.1063/1.2901172
Google Scholar
Nye, J. F. (1985). Physical properties of crystals: Their representation by tensors and matrices. Oxford University Press.
Google Scholar
Plümper, O., Wallis, D., Teuling, F., Moulas, E., Schmalholz, S. M., Amiri, H., & Müller, T. (2022). High-magnitude stresses induced by mineral-hydration reactions. Geology, 50(12), 1351–1355. https://doi.org/10.1130/g50493.1
Google Scholar
Pollard, D., & Fletcher, R. C. (2005). Fundamentals of Structural Geology. Cambridge University Press.
Google Scholar
Prieto, G. A., Beroza, G. C., Barrett, S. A., López, G. A., & Florez, M. (2012). Earthquake nests as natural laboratories for the study of intermediate-depth earthquake mechanics. Tectonophysics, 570–571, 42–56. https://doi.org/10.1016/j.tecto.2012.07.019
Google Scholar
Quesnel, D. J., Rimai, D. S., & DeMejo, L. P. (1993). Elastic compliances and stiffnesses of the fcc Lennard-Jones solid. Physical Review B, 48(10), 6795–6807. https://doi.org/10.1103/physrevb.48.6795
Google Scholar
Ranganathan, S. I., & Ostoja-Starzewski, M. (2008). Universal elastic anisotropy index. Physical Review Letters, 101(5), 055504. https://doi.org/10.1103/physrevlett.101.055504
Google Scholar
Räss, L., Utkin, I., Duretz, T., Omlin, S., & Podladchikov, Y. Y. (2022). Assessing the robustness and scalability of the accelerated pseudo-transient method. Geoscientific Model Development, 15(14), 5757–5786. https://doi.org/10.5194/gmd-15-5757-2022
Google Scholar
Reuber, G., Kaus, B. J. P., Schmalholz, S. M., & White, R. W. (2016). Nonlithostatic pressure during subduction and collision and the formation of (ultra)high-pressure rocks. Geology, 44(5), 343–346. https://doi.org/10.1130/g37595.1
Google Scholar
Richter, B., Stünitz, H., & Heilbronner, R. (2016). Stresses and pressures at the quartz-to-coesite phase transformation in shear deformation experiments. Journal of Geophysical Research: Solid Earth, 121(11), 8015–8033. https://doi.org/10.1002/2016jb013084
Google Scholar
Riel, N., Kaus, B. J. P., Green, E. C. R., & Berlie, N. (2022). MAGEMin, an Efficient Gibbs Energy Minimizer: Application to Igneous Systems. Geochemistry, Geophysics, Geosystems, 23(7), e2022GC010427. https://doi.org/10.1029/2022gc010427
Google Scholar
Schmalholz, S. M., Moulas, E., Plümper, O., Myasnikov, A. V., & Podladchikov, Y. Y. (2020). 2D Hydro-Mechanical-Chemical Modeling of (De)hydration Reactions in Deforming Heterogeneous Rock: The Periclase-Brucite Model Reaction. Geochemistry, Geophysics, Geosystems, 21(11), e2020GC009351. https://doi.org/10.1029/2020gc009351
Google Scholar
Schmeling, H., Kruse, J. P., & Richard, G. (2012). Effective shear and bulk viscosity of partially molten rock based on elastic moduli theory of a fluid filled poroelastic medium. Geophysical Journal International, 190(3), 1571–1578. https://doi.org/10.1111/j.1365-246x.2012.05596.x
Google Scholar
Schmid, D. W., & Podladchikov, Y. Y. (2003). Analytical solutions for deformable elliptical inclusions in general shear. Geophysical Journal International, 155(1), 269–288. https://doi.org/10.1046/j.1365-246x.2003.02042.x
Google Scholar
Schultz, A. J., & Kofke, D. A. (2018). Comprehensive high-precision high-accuracy equation of state and coexistence properties for classical Lennard-Jones crystals and low-temperature fluid phases. The Journal of Chemical Physics, 149(20), 204508. https://doi.org/10.1063/1.5053714
Google Scholar
Schwerdtfeger, P., Burrows, A., & Smits, O. R. (2021). The Lennard-Jones Potential Revisited: Analytical Expressions for Vibrational Effects in Cubic and Hexagonal Close-Packed Lattices. The Journal of Physical Chemistry A, 125(14), 3037–3057. https://doi.org/10.1021/acs.jpca.1c00012
Google Scholar
Sekerka, R. F., & Cahn, J. W. (2004). Solid–liquid equilibrium for non-hydrostatic stress. Acta Materialia, 52(6), 1663–1668. https://doi.org/10.1016/j.actamat.2003.12.010
Google Scholar
Shinoda, W., Shiga, M., & Mikami, M. (2004). Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Physical Review B, 69(13), 134103. https://doi.org/10.1103/physrevb.69.134103
Google Scholar
Shukla, R. C., & Cowley, E. R. (1985). Helmholtz free energy of an anharmonic crystal to O(λ4). III. Equation of state for the Lennard-Jones solid. Physical Review B, 31(1), 372–378. https://doi.org/10.1103/physrevb.31.372
Google Scholar
Spear, F. S. (1993). Metamorphic phase equilibria and pressure– temperature–time path: Vol. Mineralogi.
Google Scholar
Speziale, S., Zha, C.-S., Duffy, T. S., Hemley, R. J., & Mao, H. K. (2001). Quasi-hydrostatic compression of magnesium oxide to 52 GPa: Implications for the pressure-volume-temperature equation of state. Journal of Geophysical Research: Solid Earth, 106(B1), 515–528. https://doi.org/10.1029/2000jb900318
Google Scholar
Spiegelman, M., & McKenzie, D. (1987). Simple 2-D models for melt extraction at mid-ocean ridges and island arcs. Earth and Planetary Science Letters, 83(1), 137–152. https://doi.org/10.1016/0012-821x(87)90057-4
Google Scholar
Spycher, N. F., Sonnenthal, E. L., & Apps, J. A. (2003). Fluid flow and reactive transport around potential nuclear waste emplacement tunnels at Yucca Mountain, Nevada. Journal of Contaminant Hydrology, 62–63, 653–673. https://doi.org/10.1016/s0169-7722(02)00183-3
Google Scholar
Steinhardt, P. J., Nelson, D. R., & Ronchetti, M. (1983). Bond-orientational order in liquids and glasses. Physical Review B, 28(2), 784–805. https://doi.org/10.1103/physrevb.28.784
Google Scholar
Stephan, S., Thol, M., Vrabec, J., & Hasse, H. (2019). Thermophysical Properties of the Lennard-Jones Fluid: Database and Data Assessment. Journal of Chemical Information and Modeling, 59(10), 4248–4265. https://doi.org/10.1021/acs.jcim.9b00620
Google Scholar
Stukowski, A. (2010). Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool. Modelling and Simulation in Materials Science and Engineering, 18(1), 015012. https://doi.org/10.1088/0965-0393/18/1/015012
Google Scholar
Tajčmanová, L., Vrijmoed, J., & Moulas, E. (2015). Grain-scale pressure variations in metamorphic rocks: Implications for the interpretation of petrographic observations. Lithos, 216–217, 338–351. https://doi.org/10.1016/j.lithos.2015.01.006
Google Scholar
Thompson, A. P., Aktulga, H. M., Berger, R., Bolintineanu, D. S., Brown, W. M., Crozier, P. S., in ’t Veld, P. J., Kohlmeyer, A., Moore, S. G., Nguyen, T. D., Shan, R., Stevens, M. J., Tranchida, J., Trott, C., & Plimpton, S. J. (2022). LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Computer Physics Communications, 271, 108171. https://doi.org/10.1016/j.cpc.2021.108171
Google Scholar
Thompson, A. P., Plimpton, S. J., & Mattson, W. (2009). General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions. The Journal of Chemical Physics, 131(15), 154107. https://doi.org/10.1063/1.3245303
Google Scholar
Toxvaerd, S., & Dyre, J. C. (2011). Communication: Shifted forces in molecular dynamics. The Journal of Chemical Physics, 134(8), 081102. https://doi.org/10.1063/1.3558787
Google Scholar
Van der Molen, I., & Van Roermund, H. L. M. (1986). The pressure path of solid inclusions in minerals: The retention of coesite inclusions during uplift. Lithos, 19(3–4), 317–324. https://doi.org/10.1016/0024-4937(86)90030-7
Google Scholar
Vaughan, P. J., Green, H. W., & Coe, R. S. (1984). Anisotropic growth in the olivine-spinel transformation of Mg2GeO4 under nonhydrostatic stress. Tectonophysics, 108(3–4), 299–322. https://doi.org/10.1016/0040-1951(84)90241-5
Google Scholar
Virieux, J. (1986). P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics, 51(4), 889–901. https://doi.org/10.1190/1.1442147
Google Scholar
Wheeler, J. (2014). Dramatic effects of stress on metamorphic reactions. Geology, 42(8), 647–650. https://doi.org/10.1130/g35718.1
Google Scholar
Yang, J., & Faccenda, M. (2020). Intraplate volcanism originating from upwelling hydrous mantle transition zone. Nature, 579(7797), 88–91. https://doi.org/10.1038/s41586-020-2045-y
Google Scholar
Yarushina, V. M., & Podladchikov, Y. Y. (2015). (De)compaction of porous viscoelastoplastic media: Model formulation. Journal of Geophysical Research: Solid Earth, 120(6), 4146–4170. https://doi.org/10.1002/2014jb011258
Google Scholar
Yoo, S., Xantheas, S. S., & Zeng, X. C. (2009). The melting temperature of bulk silicon from ab initio molecular dynamics simulations. Chemical Physics Letters, 481(1–3), 88–90. https://doi.org/10.1016/j.cplett.2009.09.075
Google Scholar
Yoo, S., Zeng, X. C., & Morris, J. R. (2004). The melting lines of model silicon calculated from coexisting solid–liquid phases. The Journal of Chemical Physics, 120(3), 1654–1656. https://doi.org/10.1063/1.1633754
Google Scholar
Yu, C., Ouyang, Y., & Chen, J. (2022). Enhancing thermal transport in multilayer structures: A molecular dynamics study on Lennard-Jones solids. Frontiers of Physics, 17(5), 53507. https://doi.org/10.1007/s11467-022-1170-5
Google Scholar
Zheng, X., Cordonnier, B., Zhu, W., Renard, F., & Jamtveit, B. (2018). Effects of Confinement on Reaction-Induced Fracturing During Hydration of Periclase. Geochemistry, Geophysics, Geosystems, 19(8), 2661–2672. https://doi.org/10.1029/2017gc007322
Google Scholar
Zhou, Y., He, C., Song, J., Ma, S., & Ma, J. (2005). An experiment study of quartz-coesite transition at differential stress. Chinese Science Bulletin, 50(5), 446–451. https://doi.org/10.1007/bf02897461
Google Scholar
Zhukov, V. P., & Korsakov, A. V. (2015). Evolution of host-inclusion systems: A visco-elastic model. Journal of Metamorphic Geology, 33(8), 815–828. https://doi.org/10.1111/jmg.12149
Google Scholar

APPENDICES

APPENDIX A
Modelling the distribution of stress in granular rocks

In our model we consider a rock that is composed of two-media (solid-fluid) that have an effective viscosity ratio ηs/ηf =104. We chose this ratio as an approximation since larger contrasts in effective viscosity do not affect the magnitude of pressure variations significantly (Moulas et al., 2014; Schmid & Podladchikov, 2003). Both the solid and the fluid are treated as isotropic and incompressible. The incompressible approximation leads to identical results to the compressible case for timescales that are ~ 10 times larger than the Maxwell viscoelastic timescale (Moulas et al., 2019). The Maxwell timescale is defined here as t=η/K, where η is the effective viscosity and K is the bulk modulus of the material. Finally, we would like to note that the viscous incompressible solution is mathematically equivalent to the elastic incompressible solution if one replaces viscosities by shear moduli and velocities by displacements (Goodier, 1934). Thus, we provide the viscous solution as a suitable approximation for the development of stresses over large timescales, but the stress pattern would be identical if one considers a purely elastic rheology.

To calculate the stress pattern of figure 2 we solve the conservation of momentum for the case of slow viscous flow:

Pxj+τijxi+ρgj=0

where P is the mean stress, τij are the components of the deviatoric stress tensor, ρ is the density and gj is the component of gravity acceleration in the jth direction. The components of the deviatoric stress tensor are found as:

τij=μ(Vjxi+Vixj)

where V represents the velocity. The incompressibility condition is given by equation (A3).

Vkxk=0

The system of equations is solved in an iterative manner by introducing an artificial compressibility in equation (A3) and by utilizing a staggered, finite-difference numerical grid (Chorin, 1997; Virieux, 1986) under free slip boundary conditions. Our model was solved in dimensionless form using the vertical lengthscale (L), the lithostatic pressure (ρsgL), and the solid’s effective viscosity (ηs) as independent scales. To complete the solution, we have considered a density ratio ρs/ρf=2700/1000.

Our system of equations was solved by using the pseudo-transient approach (Chorin, 1997; Duretz et al., 2019; Räss et al., 2022). However, we have modified the iteration strategy by imposing that the stress tensor within the fluid domains is given by the following relationships:

σij=Pδij+τij=ρfg(Ly)δij

τij=0

where σij are the components of the total stress tensor and δij is the Kronecker delta. The previous modification ensures that the fluid has zero deviatoric stress and its state of stress is hydrostatic as it would be the case if the porosity was interconnected.

APPENDIX B
Equilibrium between a fluid and an elastically isotropic solid under stress

Sekerka and Cahn (2004) obtained a relation equivalent to equations (13) and (14) for an elastically isotropic solid under stress in equilibrium with a fluid of the same component. Assuming that the shear stresses in the x1x2 plane are zero (σ12=0), they derived the following expression for the change of the equilibrium temperature of the system due to the application of a nonhydrostatic stress in the solid at constant pressure in the fluid (eq 20 of Sekerka & Cahn, 2004):

TTH= ωsH2 ΔsH E[(qs11)22νqs11qs22+(qs22)2]

where E and ν are the Young’s modulus and Poisson’s ratio of the solid in the hydrostatic state, respectively. In the same way, for an elastically isotropic solid, the change in equilibrium pressure of the fluid due to the application of a nonhydrostatic stress at constant temperature can be described as:

PfPfH= ωsH2 ΔωHE[(qs11)22νqs11qs22+(qs22)2]

For an elastically isotropic solid, the formulations of equations (13)–(B1) and (14)–(B2) are equivalent, as can be demonstrated by recognizing that for an elastically isotropic material the following equivalences hold true: D1111=D2222=1E, D1122= νE. However, for an elastically anisotropic solid without shear stresses (qs12=0), the use of isotropic relationships (eqs B1 and B2) would lead to a discrepancy compared to the corresponding anisotropic solutions (eqs 13 and 14, respectively). The relative discrepancy between the isotropic and the anisotropic formulations is controlled by the elastic anisotropy of the solid. For an elastically anisotropic solid under an axially-symmetric stress state (qS11=qs22, qs12= qs13=qs23=qs33=0) the relative discrepancy between the isotropic and anisotropic formulations becomes:

(TTH)iso (TTH)aniso 1=(PfPfH)iso (PfPfH)aniso 1=2(1