Elastic geothermobarometry is a method of determining metamorphic conditions from the excess pressures exhibited by mineral inclusions trapped inside host minerals. An exact solution to the problem of combining non-linear Equations of State (EoS) with the elastic relaxation problem for elastically isotropic spherical host-inclusion systems without any approximations of linear elasticity is presented. The solution is encoded into a Windows GUI program EosFit-Pinc. The program performs host-inclusion calculations for spherical inclusions in elastically isotropic systems with full P-V-T EoS for both phases, with a wide variety of EoS types. The EoS values of any minerals can be loaded into the program for calculations. EosFit-Pinc calculates the isomeke of possible entrapment conditions from the pressure of an inclusion measured when the host is at any external pressure and temperature (including room conditions), and it can calculate final inclusion pressures from known entrapment conditions. It also calculates isomekes and isochors of the two phases.

Geobarometry, inclusion, Equations of State, elasticity

The determination of the remnant pressures in inclusions, as measured by X-ray diffractometry, birefringence analysis, or Raman spectroscopy, provides an alternative and complementary method to conventional geothermobarometry by using elasticity theory. A remnant pressure in an inclusion is developed because the inclusion and the host have different thermal expansions and compressibilities, and therefore the inclusion does not expand in response to P and T as would a free crystal. Instead it is restricted by the host mineral, and this confinement can result in inclusions exhibiting over-pressures, or under-pressures, when the host is studied at room conditions. By measuring the remnant pressure the possible temperatures and pressures of entrapment can be calculated by using the elastic properties of the host and inclusion minerals. This basic concept has been known for a long time (Rosenfeld and Chase 1961). Difficulties arise because the classic solutions for the stress distribution in host-inclusion systems (e.g., Goodier 1933; Eshelby 1957) are derived for linear elasticity, which assumes that the stresses and strains are small, and that the elastic properties do not change with pressure or temperature. However, minerals are subject to large changes in pressure and temperature from formation to room conditions, so their elastic properties are not constant but are described by non-linear Equations of State (EoS).

Several approaches have been used to apply the classic host-inclusion elastic solutions to mineral systems. All of them assume that the two minerals are elastically isotropic, and that the inclusion is spherical and isolated from the host surface and any other inclusions or defects in the host mineral. The simplest approach has been to ignore the variation of the elastic properties of minerals with pressure and temperature (e.g., Zhang 1998). This leads to errors in inclusion pressures, especially when they are calculated for prograde metamorphic conditions following entrapment (e.g., Angel et al. 2014b). A second approach has been to calculate the evolution of the inclusion pressure in a series of small steps from entrapment conditions by adjusting the elastic properties of the host and inclusion at each step according to either a full or approximate EoS, and then using the linear solution at each step to calculate mechanical equilibrium (Gillet et al. 1984; van der Molen and van Roermund 1986; d'Arco and Wendt 1994). A third approach is to consider the “thermodynamic pressure”, Pthermo, in the inclusion when it is constrained to have the same volume change as the host crystal from entrapment Ptrap and Ttrap to the final external Pend and Tend (Fig. 1). Pthermo is different from the final external pressure on the host, and this drives a further mutual elastic relaxation that reduces the difference between the inclusion pressure and Pend. This relaxation must be calculated in a second step. The advantages of this approach are that the calculation of Pthermo can be exact by using appropriate non-linear EoS, and the only linear elasticity approximation is in the relaxation term. However, the correct solution for the pressure in the spherical inclusion requires that the relaxation is evaluated during isothermal decompression from a state of uniform stress (Goodier 1933), and not along any P-T path as often incorrectly assumed (e.g., Guiraud and Powell 2006). The first step is therefore to consider a temperature change from Ttrap to Tend and to calculate the change in external pressure required to induce an equal pressure change in the inclusion (Fig. 1). This thermodynamic path is an isomeke of the host and inclusion phases (Rosenfeld and Chase 1961; Adams et al. 1975). The pressure, Pfoot, on the entrapment isomeke at Tend can be determined from Ptrap and Ttrap and the EoS for the host and inclusion phases. The second step is to calculate the pressure change in the inclusion during an isothermal change in the external pressure from Pfoot to the final Pend. Angel et al. (2014b) calculated this pressure change in the inclusion with the solution from Goodier (1933) and the assumption that the pressure derivatives of the elastic moduli of the host and inclusion are the same. In this paper we now present the solution to the problem in a way that avoids any such approximations of linear elasticity. The solution is encoded into a GUI program, EosFit-Pinc, which calculates entrapment conditions from measured remnant inclusion pressures, and vice versa.

As in Angel et al. (2014b) the full EoS of the two minerals are first used to calculate the external pressure Pfoot on the entrapment isomeke at the final temperature (Fig. 1), where the host and inclusion are at the same pressure and temperature. For the second step of the calculation we start from Pfoot and we first consider a small change in pressure dPinc imposed internally in the inclusion and, simultaneously, a small pressure change dPext imposed externally on the host. We then calculate the radial displacement u of the host-inclusion boundary necessary to return the system to mechanical equilibrium. Displacements of points at a distance R from the center of a spherically symmetric system are governed by the equation u = AR + (B/R2). The parameters A and B are constants of integration determined by the particular boundary conditions. For a small finite inclusion of radius rinc much less than the host, the constants A and B can be obtained from the more general solution to the “pressurized hollow sphere problem” (e.g., Bower 2010, www.solidmechanics.org) as:
(1)
KH and GH are, respectively, the bulk and shear moduli of the host. The fractional volume change (strain) of the inclusion from the initial conditions to the final conditions is
(2)
If the displacement u is small compared to the radius of the inclusion, as for the relaxation of mineral inclusions, then the volume strain of the inclusion due to applying the inner and outer pressures is
(3)
(Bower 2010). Substituting for the parameters A and B yields:
(4)
Because the isothermal bulk modulus is defined as
(5)
the first term in Equation 4 is the volume strain, εhost, of the free host crystal unaffected by the inclusion, so:
(6)

The second term is therefore the fractional volume change of the inclusion due to elastic relaxation. The important point is that this solution has been derived from considerations of force balance, and we are only calculating the displacement of the host-inclusion boundary in terms of a constant pressure applied to it from within the inclusion, and the external pressure on the host. There are no assumptions about elastic properties of the material inside the inclusion, and there are no assumptions of linear elasticity (which would mean constant values of KH and GH). Because the force balance is calculated for the final conditions, the value of GH to be used is that at the final conditions with the host under the external pressure change dPext. The only assumption used here is that GH does not change as a result of the deviatoric stresses that develop in the host around the inclusion (e.g., Eshelby 1957; Zhang 1998; Howell and Nasdala 2008).

The result can also be written (by following Bower 2010) in terms of finite changes in pressure from Pfoot to the final Pinc and Pend as:
(7)

The second term for the fractional volume change of the inclusion upon relaxation has been derived previously (e.g., Zhang 1998). This expression is then used, sometimes with constant KI and GH for all P-T conditions, in several programs (e.g., Ashley et al. 2014; Kohn 2014). These programs follow Guiraud and Powell (2006) and perform a one-step calculation of the entrapment conditions from the measured Pinc. In doing so, they incorrectly mix in Equation 7 the volume strains from entrapment to final conditions (as εhost and εinc) with the volume strain of relaxation, 3(Pinc – Pend)/4GH, which is relative to Pfoot. The consequence is that entrapment conditions estimated in this way do not fall exactly on a thermodynamic isomeke (see Appendix online1).

By applying the hollow sphere solution to the isothermal decompression of the inclusion and host from the isomeke at Tend, we have shown that Zhang's (1998) widely used assumption of linear elasticity to determine the pressure change on relaxation was un-necessary. Instead, we can use the full non-linear EoS of host and inclusion to evaluate their mutual relaxation and the final inclusion pressure Pinc provided we start this step of the calculation from Pfoot on the entrapment isomeke. The final pressure in the inclusion cannot be calculated directly because εinc depends on Pinc, so an iterative approach has to be used, with the following steps:

  1. Calculate the Pfoot of the entrapment isomeke from the given entrapment conditions and the full P-V-T EoS of the two phases.

  2. Calculate the volume strain of the host εhost on changing pressure from Pfoot to Pend, from the EoS of the host. This value stays fixed.

  3. Estimate a Pinc.

  4. Calculate εinc=εhost+3(Pinc−Pend)4GH.

  5. Calculate a new Pinc from εinc and the full EoS of the inclusion.

  6. Use this new value of Pinc in step 4.

  7. Repeat steps 4-5-6 until the iteration converges to a final value of Pinc.

To calculate entrapment conditions from a known Pinc, the Pfoot that produces the observed Pinc is found by the same iterative cycle, and then the entrapment isomeke can be calculated using the EoS of the two phases.

EosFit-Pinc is a Windows GUI program that performs host-inclusion elasticity calculations for spherical inclusions in an effectively infinite matrix for elastically isotropic materials by implementing the solution to the host-inclusion relaxation problem described above. The EosFit-Pinc GUI consists of a series of tabs. When the program is started the calculation tabs are not active until valid EoS are loaded into the program. The Settings tab (Fig. 2) provides the user with options for the GUI, including setting the temperature scale (Kelvin or Celsius) for the display. The display mode has an “expert” setting in which more intermediate results of calculations, and the results from different elastic relaxation models (e.g., Zhang 1998; Angel et al. 2014b), are displayed when relevant. EoS parameters for the host and inclusion can be loaded into the program as .eos files, the standard file format for transferring EoS parameters between programs in the EosFit program suite (Angel et al. 2014a; Gonzalez-Platas et al. 2016). These files can be created by the other EosFit programs, which can be launched directly from the EosFit-Pinc GUI, or EoS parameters can be imported directly from version ds62 of the Thermocalc database (Holland and Powell 2011). Other import options can be added if required.

Figure 2 shows a calculation of entrapment conditions for a cubic ferropericlase inclusion within a diamond. The remnant inclusion pressure of 1.139 GPa was calculated from the unit-cell parameters of sample GU4A (Hutchison 1998) using the EoS for ferropericlase determined from the measurements of Mao et al. (2011). The user enters the remnant pressure and the external conditions under which it was measured into the GUI and sets the temperature range for the calculation of the entrapment isomeke. The P,T points along the entrapment isomeke representing possible entrapment conditions are then displayed in the lower results panel of the GUI. If the entrapment conditions are known or estimated, for example from the petrology and chemistry of the rock, then the CalcPinc tab of the GUI allows the remnant inclusion pressure to be calculated. Calculations of the isochors of the two phases separately, or their common isomekes can also be performed, and the results cut-and-pasted from the output window to other programs for plotting.

EosFit-Pinc is written in Fortran-95 using the CrysFML (Rodriguez-Carvajal and Gonzalez-Platas 2003) library, with the cfml_eos module (Angel et al. 2014a). The program is free for non-commercial use and does not require any commercial software or libraries other than those provided with the program. It is available for download from www.rossangel.net for Windows operating systems, together with .eos files for common minerals and complete help documentation. The same calculations can be performed with the EosFit7c console program, which also provides results for other relaxation models and is available for Mac, Linux, and Windows from the same web site.

The EosFit-Pinc program provides an easy-to-use GUI for rapid calculations of inclusion pressures of minerals. It allows more flexibility than recently published programs and methods (e.g., Ashley et al. 2014; Kohn 2014) because the EoS are not built-in, but are provided to the program by the user. For quartz in garnet with Pinc < 0.3 GPa, all three programs predict entrapment conditions within ~50 bars. When Pinc is 1 GPa, the use of linear elasticity and the incorrect path for the relaxation calculation in other programs (Ashley et al. 2014; Kohn 2014) becomes significant and results in the entrapment conditions being under-estimated by more than 500 bars (0.05 GPa) at 800 °C compared to EosFit-Pinc. EosFit-Pinc is not restricted to calculating inclusion pressures when the host is at room conditions. For example, it can be used to calculate that a quartz inclusion (EoS from Angel et al. 2017) trapped inside an almandine garnet (Milani et al. 2015) on a pro-grade metamorphic path at 400 °C and 0.7 GPa only experiences a pressure of 1.8 GPa when the garnet reaches the conditions of the quartz=coesite phase boundary at 750 °C, 2.7 GPa. And when the host reaches conditions in the diamond stability field at ~4.4 GPa and 800 °C, the inclusion will be at 2.6 GPa. By contrast, by importing the rutile EoS from ds62 of Thermocalc (Holland and Powell 2011) one can calculate with EosFit-Pinc that a rutile inclusion trapped at the same conditions as the quartz exhibits pressures within 0.2 GPa of the external conditions on the same metamorphic path, as a consequence of the smaller contrast in bulk moduli between rutile and garnet.

1Deposit item AM-17-96190, A full algebraic proof is provided in the Appendix. Deposit items are free to all readers and found on the MSA web site, via the specific issue's Table of Contents (go to http://www.minsocam.org/MSA/AmMin/TOC/2017/Sep2017_data/Sep2017_data.html).

Software development and analysis was supported by ERC starting grant 307322 to Fabrizio Nestola, and by the MIUR-SIR grant “MILE DEEp” (RBSI140351) to Matteo Alvaro. We thank Javier Gonzalez-Platas (La Laguna) for continuing collaboration and development of the CrysFML, and Frank Spear (RPI) and Kyle Ashley (Austin, Texas) for comparison calculations, detailed discussions, and helpful reviews.