Elastic geothermobarometry: Corrections for the geometry of the host-inclusion system

Application of conventional geothermobarometry is extremely challenging in many rock types due to alteration processes, chemical reequilibration and diffusion, and kinetic limitations. Elastic geothermobarometry on host-inclusion systems is a complementary method independent of chemical equilibrium. An inclusion completely entrapped inside a host mineral will develop a residual stress upon exhumation because of the contrast in elastic properties (Rosenfeld and Chase, 1961). If the host does not undergo plastic deformation or brittle failure after trapping the inclusion, the entrapment pressure (P_trap) can be calculated from the measured residual pressure on the inclusion (or remnant pressure, P_inc), provided that the elastic properties (equations of state, EoS) for the host and inclusion are known (e.g., Zhang, 1998; Angel et al., 2014). Elastic geothermobarometry is increasingly applied to metamorphic rocks, where measurements of Raman shifts on quartz inclusions trapped in garnet (quartz-in-garnet inclusion barometry, QuiG) give information on the residual stresses that can be used to infer growth conditions (e.g., Kouketsu et al., 2016) and the degree of overstepping of garnet isograds (e.g., Spear et al., 2014).

The validity of elastic geobarometric methods was discussed by Ashley et al. (2016), who reported that the P_trap inferred from measured P_inc of quartz inclusions in garnets do not match those obtained by conventional geobarometry on the same rocks. However, the calculation of P_trap currently assumes that the minerals are elastically isotropic with ideal geometry where the inclusion is spherical and isolated at the center of the host (Goodier, 1933; Eshelby, 1957; Van der Molen and Van Roermund, 1986). None of these conditions apply in natural systems; neither quartz nor garnet are elastically isotropic, inclusions are often close to grain boundaries or other inclusions, and they are often not spherical. The resulting changes in P_inc can only be quantified using numerical approaches.

In this paper, we use finite element (FE) models of elastically isotropic host-inclusion systems with non-ideal geometries to determine the magnitude of the geometric effects on P_inc, and in turn on the calculated P_trap. We show that the discrepancies reported by Ashley et al. (2016) are only partly due to the geometry of their samples. We provide guidelines as to which geometries of host-inclusion systems lead to deviations smaller than the typical experimental uncertainties in inclusion pressures obtained from conventional μ-Raman measurements, and can therefore be safely used for geobarometry without any correction.

The final stress state of an inclusion is path independent, and it is convenient to split the pressure-temperature (P-T) change from the entrapment conditions (P_trap, T_trap; see Fig. 1) to the final pressure and temperature (P_end, T_end) into two parts (see Angel et al., 2014). Figure 1 illustrates the stepwise procedure used to calculate the residual pressure from known entrapment conditions. During step 1 the temperature is reduced from T_trap to T_end along the isomeke (Rosenfeld and Chase, 1961; Adams et al., 1975), thus preserving the reciprocal mechanical equilibrium between the host and the inclusion. The change in external T and P required to maintain the pressure in the inclusion equal to the external P can be calculated directly from the thermodynamic properties of the minerals without any influence of the geometry of the system. In step 2, the isothermal decompression from P_foot, T_end to the final P_end, T_end (P_footT_end and P_end as defined by Angel et al. [2014]) causes a mechanical disequilibrium between the host and the inclusion. Consequently, the stresses are readjusted through the relaxation process. Because the relaxation depends on force balance at the interface between host and inclusion, in this step the geometry becomes important. The exact amount of relaxation in step 2 can only be calculated if the geometry of the system is ideal; for all other cases a numerical approach is required.

In our study we used two commercially available engineering packages (MARC Mentat by MSC Software, http://www.mscsoftware.com/product/marc/, and Abaqus by Dassault Systèmes, https://www.3ds.com/products-services/simulia/products/abaqus/) to create and solve two-dimensional (2-D) axisymmetric and 3-D models using FE numerical simulations. Always using isotropic elastic properties, we explored the effects of several deviations from ideal geometry, including the size of the inclusion relative to the host and its proximity to external surfaces. To evaluate the effects of non-spherical shapes we modeled ellipsoids of revolution with aspect ratios 1:1:1, 2:1:1, 1:2:2, 5:1:1, and 1:5:5. The effects of edges and corners were then determined by comparing the results against cylindrical and prismatic models (with quadrilateral cross sections) with the same aspect ratios.

To simulate the effects of external pressure, edge loads (for 2-D models) or face loads (for 3-D models) were applied to the external boundaries of the models. Stationary boundary conditions were placed on the relevant edges and faces to prevent rigid body rotations and translations. An example of a model mesh and the elastic properties used in the models are given in the GSA Data Repository¹; more details are in Burnley and Davis (2004), Burnley and Schmidt (2006), and Abaqus (2016). For each model, we performed calculations using different elastic isotropic properties for the host and the inclusion to probe possible scaling laws.

This corrected P_inc can then be used to calculate the P_trap using isotropic elastic geobarometry models (e.g., Angel et al., 2014, 2017b).

Our FE models have been validated against the analytical exact solution by modeling an ideal infinite spherical system. In practice, the host can be considered infinite when the simulation results do not change upon further increase in the size of the host (Fig. 2A). Our FE models then reproduce the analytical solution for the pressure inside a spherical inclusion well within the expected numerical precision (i.e., 0.2%).

The stress in the region of the host close to the inclusion is always deviatoric (e.g., Zhang, 1998). Therefore, when a large inclusion is surrounded by a thin layer of host crystal, the deviatoric stress extends throughout the volume of the thin host layer, causing the outer boundary of the host to deform. The host is thus no longer able to shield the inclusion from the external pressure. Consequently, the P_inc will be partially released. For a spherical inclusion at the center of the host, the pressure release is a function of the size and the properties of the inclusion with respect to the host (Fig. 2A). Hosts much stiffer than the inclusion (e.g., quartz in garnet) can preserve a larger P_inc. Our results indicate that, if the radius of the host is at least four times that of the inclusion, both the P_inc^non-ideal and the P_trap^non-ideal are within 1% of the value expected for an infinite host. For the same reason, the capacity of the host to act as a pressure vessel for the inclusion is also reduced when a small inclusion is close to the external surface of the host (Fig. 2B). Stiffer hosts preserve more residual pressure than softer hosts. Regardless of the contrast in elastic properties, if the inclusion is at least 3 radii from the external surface of the host the effect on P_inc is <1%. If a spherical and isotropic inclusion is close to the external surface of the host, the normal stresses in the inclusion are not homogeneous, and the domains of the inclusion closer to the external surface record stresses lower than those toward the center of the host. For a quartz inclusion in pyrope the variation of the pressure across the inclusion can reach 8% when the distance to the surface is half the radius of the inclusion (inset in Fig. 2B). Note that these conclusions do not depend on the absolute size of the inclusion, but upon the relative sizes of the inclusion and host.

For fluid inclusions the aspect ratio and the presence of corners and edges are two major influences on the pressures of isolated inclusions (e.g., Burnley and Davis, 2004; Burnley and Schmidt, 2006). In our models of solid inclusions, we find that the aspect ratio of the inclusion gives rise to deviations in P_inc >7% for soft platy inclusions (aspect ratio 1:5:5) in stiff hosts (e.g., quartz in pyrope; see Fig. 3). The presence of edges and corners further enhances the deviations (≈9%). For non-spherical shapes with edges and corners, the stress in the inclusion is neither homogeneous nor hydrostatic. The pressure varies from the center of the inclusion toward its external surface, by different amounts in different directions. For a quartz inclusion with aspect ratio 1:5:5 in pyrope, the pressure variation along the longer axes of the inclusion is ∼5%, while it is <1% along the shortest axis (see Figs. DR2 and DR3 in the Data Repository). For a residual pressure at the center of the inclusion of 0.3 GPa, the differential stress (σ_max – σ_min) within the inclusion reaches 0.28 GPa. For a stiff inclusion in a soft host with the same shape, the pressure variation within the inclusion is typically much larger (e.g., 22% for diamond in pyrope) and of the opposite sign (see Fig. DR2).

The exact effect of inclusion shape on P_inc^non-ideal is a complex interplay between the bulk and shear moduli for both host and inclusion (see Fig. 3). In general, the influence of non-ideal shapes becomes greater when the bulk modulus of the host and the inclusion are similar, provided there is a significant contrast in shear moduli. For a soft inclusion in a stiffer host (quartz in garnet, or pyrope in diamond) with aspect ratios less than 1:2:2, the deviations induced by the shape are typically <5% (Fig. 3).

Ashley et al. (2016) used Raman spectroscopy to determine the remnant pressure in quartz inclusions in garnets while they were heated up to 500 °C. As the T increases, the P_inc^exp increases because of thermal pressure effects, but the P_trap for a single inclusion should always be a unique value independent of the temperature (Tend) at which the P_inc^exp is measured. However, Ashley et al. (2016) reported large variations on P_trap for the same inclusion calculated from the various P_inc^exp measured at different T_end and none of the calculated P_trap agreed with the results from conventional geobarometry. Ashley et al. (2016) ascribed this unphysical behavior to the use of unrealistic EoS for quartz close to the α-β structural phase transition. We chose this example to assess if the shape of the inclusion could explain these discrepancies.

We consider the case of sample MT 09-09, where several quartz inclusions are entrapped in an almandine-rich garnet (Ashley et al., 2015, 2016). The P_trap values at 540 °C were recalculated from the experimental P_inc^exp of 0.300 GPa and 0.491 GPa at the minimum and maximum T_end (31 °C and 500 °C) using a more reliable EoS for quartz (Angel et al., 2017a) that explicitly includes the α-β transition. Assuming ideal geometry for the quartz inclusion, the discrepancy between the two P_trap values is 0.186 GPa (Table 1), similar to that reported by Ashley et al. (2016), demonstrating that the differences cannot be ascribed to errors in the EoS. To eliminate the discrepancies in P_trap values the volume thermal expansion of almandine must be increased by more than 30% to α_298K ∼2.76 × 10⁻⁵ K^–1; this is unrealistic given that this value is much greater than those of any silicate garnet end member.

Because the shapes of the inclusions measured by Ashley et al. (2016) were not reported, we then overestimated the shape effects by modeling the inclusion as a platy prism (aspect ratio 1:5:5). At room temperature the correction factor is Γ = −0.094, similar to that for quartz in pyrope (Fig. 3), but decreases to Γ = −0.078 at 500 °C (Table 1) due to the elastic softening of quartz as it approaches the phase transition (Lakshtanov et al., 2007). The inclusion pressures corrected for shape, P_inc^corrected, are then 0.331 GPa (at 31 °C) and 0.532 GPa (at 500 °C), and result in a small but insignificant reduction of 0.02 GPa in the differences in P_trap calculated from the two measurements. Furthermore, even with the geometrical correction, the P_trap values (1.091 and 0.929 GPa) are not in agreement with those obtained from the conventional methods (0.82 GPa; Ashley et al., 2015).

Thus neither the EoS nor the shape of the inclusion can explain the discrepancies found by Ashley et al. (2016), and other factors not yet included in the current models must be responsible for the discrepancies in the P_trap. One factor is that quartz inclusions in a garnet host will be subject to isotropic strain (leaving aside further perturbations arising from the elastic relaxation and the geometry) because garnet is cubic. As quartz is elastically anisotropic, the isotropic strain will result in a non-hydrostatic stress in the inclusion. The effect of this deviatoric stress on the Raman spectrum of quartz is not known in detail, but both theory (Key, 1967) and experiments (Briggs and Ramdas, 1977) show that Raman peak shifts will be different from those predicted from hydrostatic calibrations used by Ashley et al. (2016) to convert measured Raman shifts into pressures. Therefore, the mismatch in the P_trap is probably due to the combination of an inappropriate Raman stress calibration and the assumption of elastic isotropy in the geobarometric models.

Current elastic geobarometric models assume isotropic elastic properties for the host and the inclusion, and that the inclusions are isolated and spherical. These conditions do not commonly occur in natural rocks. Regardless of the relative stiffness of host and inclusion, for a big inclusion in a small host and for an inclusion close to the external surface of the host, the P_inc^exp is reduced relative to the ideal case, but a simple correction factor cannot be defined and Γ should be evaluated on a case-by-case basis with finite element method (FEM) analysis carried out on realistic digital models of the inclusions.

For isotropic elasticity, our FEM results show that for an inclusion at least at 3 radii from external surfaces or other inclusions, the geometric effects on P_inc^exp are <1% (Fig. 2B). Under these conditions the shape effects then dominate the geometric corrections to the measured P_inc^exp (Fig. 3). For soft inclusions in a stiff host (e.g., quartz in garnet), non-spherical inclusions (Γ< 0), will exhibit a lower pressure than spherical inclusions. Correction of the measured P_inc^exp for the shape effects will therefore result in P_inc^corrected > P_inc^exp and thus an increase in the calculated P_trap. By contrast, for stiff inclusions in soft hosts (Γ > 0), the correction will lead to P_inc^corrected < P_inc^exp and therefore a reduction in P_trap. Experimental uncertainties on P_inc^exp are typically <5% when measured by Raman spectroscopy (e.g., Ashley et al., 2016, Kouketsu et al., 2016). For P_inc <1 GPa, the uncertainties propagated into the P_trap are smaller than those on the P_inc. Therefore, pressures from inclusions for which the geometrical effects on P_inc^exp are <5% will provide reliable estimates of P_inc, and hence P_trap without the need for correction. For soft inclusions in stiff hosts, such as quartz in garnet, this means the following:

• (1) The radius of the inclusion must be smaller than one-half of that of the host.

• (2) The distance from the external surface is larger than one-half the radius of the inclusion.

• (3) The inclusion aspect ratio is lower than 1:3:3, with few sharp edges and corners.

These guidelines do not apply to inclusions stiffer than the host (e.g., diamond in garnet) that require much larger corrections of opposite sign (Fig. 3).

This project received funding from the European Research Council under the European Union’s Horizon 2020 research and innovation program grant agreements 714936 to Alvaro and 307322 to Nestola. Alvaro has also been supported by the MIUR-SIR (Ministry of Education, University and Research–Scientific Independence of Young Researchers, Italy) grant MILE DEEp (Mineral Inclusion Elasticity for a New Deep Subduction Geobarometer; RBSI140351). Support for P. Burnley and computational resources was provided by the U.S. National Nuclear Security Administration under the Stewardship Science Academic Alliances program through U.S. Department of Energy Cooperative Agreement DE-NA0001982. We thank Matthew Steele-MacInnis and two anonymous reviewers for very constructive comments on the manuscript.