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

Elastic geothermobarometry on inclusions is a method to determine pressure-temperature conditions of mineral growth independent of chemical equilibrium. Because of the difference in their elastic properties, an inclusion completely entrapped inside a host mineral will develop a residual stress upon exhumation, from which one can back-calculate the entrapment pressure. Current elastic geobarometric models assume that both host and inclusion are elastically isotropic and have an ideal geometry (the inclusion is spherical and isolated at the center of an infinite host). These conditions do not commonly occur in natural rocks, and the consequences for inclusion pressures can only be quantified with numerical approaches. In this paper, we report the results of numerical simulations of inclusions with the finite element method on elastically isotropic systems. We define and determine a geometrical factor (Γ) that allows measured residual pressures to be corrected for the effects of non-ideal geometry. We provide simple guidelines as to which geometries can safely be used for elastic geobarometry without correcting for the geometry. We also show that the discrepancies between elastic and conventional geobarometry reported in literature are not due to geometrical effects, and therefore result from other factors not yet included in current models.


INTRODUCTION
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.

METHODS
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 foot T 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 1 ; 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.
For each geometry we calculated the actual inclusion pressure P inc non-ideal by performing FE simulations upon isothermal decompression from P foot to P end . We define a geometrical factor (hereafter Γ) as the normalized deviation of the actual inclusion pressure from that expected for an ideal isolated spherical inclusion, P inc ideal , for the same decompression: (1) The value of Γ is obtained using the linear elastic approximation, so it is independent of the magnitude of P foot -P end . Because pressures in natural inclusions are typically <1 GPa (e.g., Ashley et al., 2016), this linear approximation is not significant for most inclusions, and the Γ 1 GSA Data Repository item 2018057, elastic geothermobarometry (computational details, elastic properties used for all calculations, and some additional examples), is available online at http://www.geosociety.org/datarepository/2018/ or on request from editing@geosociety.org. parameter can be used to correct experimentally determined inclusion pressures (P inc exp = P inc non-ideal ) for geometric effects: (2) This corrected P inc can then be used to calculate the P trap using isotropic elastic geobarometry models (e.g., Angel et al., 2014Angel et al., , 2017b.

Insights from Finite Element Models
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 Figure 1. Calculation procedure for two host inclusion pairs; one with ideal geometry (red dot inclusion) and one with non-ideal geometry (green dashed rectangle), that are trapped at the same pressure and temperature (P trap , T trap ) conditions.
Step 1: Along the isomeke the host and inclusion are in reciprocal mechanical equilibrium. Therefore, the pressure on the isomeke at the final T end (i.e., P foot ) will be the same for any geometry of the system.
Step 2: The host is decompressed to the final pressure P end . The relaxation of the inclusion is geometry dependent and therefore the final P inc will be different for the two systems (P inc non-ideal ≠ P inc ideal ). The geometrical factor Γ is a measure of this discrepancy. The inset illustrates how to apply Γ to correct the experimental P inc exp measured on natural rocks with non-ideal geometry. The corrected P inc corrected can then be used to back-calculate the P trap using currently available elastic geobarometry models.
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 (T end ) 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.

Calculation of Entrapment Pressures
We consider the case of sample MT 09-09, where several quartz inclusions are entrapped in an almandine-rich garnet (Ashley et al., 2015(Ashley et al., , 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 −5 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  Note: P trap is entrapment pressure. ∆P trap is calculated as the difference between P trap from the pressure on the inclusion (P inc ) at the end temperature T end = 31 °C, and that from the P inc at 500 °C (exp is experimental). The equations of state used for quartz and almandine are reported in the GSA Data Repository (see text footnote 1). Figure 3. Geometrical factor Γ for several shapes plotted versus the normalized aspect ratio. The latter is calculated with the unique axis as the denominator (e.g., aspect ratio 2:1:1 becomes ½ = 0.5). For a soft inclusion in a stiffer host (e.g., quartz in garnet), Γ < 0 and therefore P inc non-ideal < P inc ideal (as in Fig. 1). The opposite occurs for a stiff inclusion in a softer host (e.g., diamond in pyrope). Note that Γ values greater than zero are plotted with a compressed vertical scale. Exp.-experimental. 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.

CONCLUSIONS
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).