Applications and Limitations of Elastic Thermobarometry: Insights From Elastic Modeling of Inclusion‐Host Pairs and Example Case Studies

Elastic thermobarometry can be used to constrain the pressure and temperature conditions of mineral crystallization by exploiting the difference in the elastic evolution of a mineral inclusion and its host during cooling and decompression. In this work we examine the pressure‐temperature sensitivity of >5,000 untested inclusion‐host pairs. Hosts such as diamond and zircon are ideal host minerals because their low compressibility makes them rigid containment vessels. Highly compressible inclusions such as albite, graphite, and quartz serve as the most reliable barometers. We provide three case studies of inclusion‐host pairs from different geologic settings to demonstrate the advantages and challenges associated with these mineral pairs. Apatite inclusions in olivine from Yellowstone caldera mostly record negative residual pressures (tension) and suggest magmatic crystallization at ~0.4 GPa. Rutile inclusions in garnet from Verpeneset eclogites record near ambient conditions and do not recover reasonable metamorphic conditions of rutile entrapment. These results suggest that stiff inclusions may have a tensile strain limit, a possible limitation of elastic thermobarometry. Albite inclusions in epidote from a blueschist (Syros, Greece) record geologically reasonable entrapment pressures, but a large range of residual pressures that may be caused by the complex anisotropy of both phases. Our theoretical and applied results indicate that elastic thermobarometry has the potential to be used to understand petrologic processes in diverse geologic environments, including mantle, metamorphic, and magmatic settings but that each elastic thermobarometer requires careful evaluation.


Introduction
Methods that can directly monitor geologic processes that occur in the deep subsurface do not exist. Active processes in magmatic, metamorphic, and mantle environments remain at the forefront of geoscience research because they control tectonism and crustal deformation, seismic and volcanic hazards, the geomorphic evolution of landscapes, and the generation of economic mineral resources. Much effort has been directed toward advancing methods to understand subsurface processes.
Recent work has improved the reliability and the applicability of elastic thermobarometry through some of the following advances: updated mineral thermodynamic properties (e.g., , by providing shape corrections for non-ideal inclusions (e.g., Mazzucchelli et al., 2018), by quantifying the effect of anisotropic strain on calculations (e.g., Murri et al., 2018;Stangarone et al., 2019;Zhong et al., 2019), by advancing the theory of the technique (e.g., Mazzucchelli et al., 2019;Zhong et al., 2018Zhong et al., , 2020Zhukov & Korsakov, 2015), by easing the implementation of elastic modeling for users (e.g., Angel, Mazzucchelli, et al., 2017), and by comparing entrapment conditions calculated from elastic thermobarometry with independent constraints or known experimental conditions of garnet synthesis (e.g., Bonazzi et al., 2019;Thomas & Spear, 2018). Elastic thermobarometry has primarily been applied to diamond-and garnet-hosted inclusions (e.g., Ashley et al., 2014;Cohen & Rosenfeld, 1979;Howell et al., 2010;Izraeli et al., 1999;Nestola et al., 2019;Thomas & Spear, 2018), because both of these host minerals are nearly isotropic and they should impose a near isotropic stress on the inclusions. However, recent work suggests that anisotropic inclusion-host pairs may similarly provide reliable P-T conditions (Cisneros et al., 2020), yet many inclusion-host pairs that could be used to extract P-T constraints have not been explored. Here, we evaluate the theoretical potential of >5,000 inclusion-host pairs as barometers or thermometers, building upon the work introduced by Kohn (2014) (48 inclusion-host pairs). Most of these inclusion-host pairs have never been tested, and significant additional work is needed to fully evaluate the inclusion-host pairs that we present. Our contribution solely presents the theoretical potential of a suite of diverse inclusion-host pairs.
We demonstrate three simple worked examples to illustrate advantages and challenges associated with these mineral pairs. Apatite inclusionsinolivine (ap-in-ol) from Yellowstone caldera and albite inclusions in epidote (ab-in-ep) from Syros, Greece, record a complex array of entrapment conditions; however, some calculated entrapment pressures agree with P-T constraints from previous studies. Residual inclusion pressures recorded by rutile inclusions in garnet (ru-in-grt) in an eclogite from Verpeneset, Norway, are far different than expected residual pressures based on established metamorphic conditions of garnet growth and may suggest that stiff inclusions in hosts have a tensile strain limit. We also provide a script with a graphical user interface (GUI) that allows for evaluation of our code. The script can be used to calculate entrapment pressures and temperatures and allows linear mixing of molar volumes and shear moduli to account for mineral solid solutions.

Elastic Modeling Calculations: An Overview
The difference in compressibility and thermal expansivity between a crystalline inclusion and host generates compression or expansion of the inclusion as P-T conditions change. At ambient temperatures, encapsulated inclusions preserve a residual strain, which can be converted to an inclusion pressure (P inc ). Residual strains can be quantified by multiple techniques, including measurement of birefringent halos surrounding inclusions (e.g., Howell et al., 2010;Rosenfeld & Chase, 1961), X-ray Diffraction (e.g., Harris & Munn, 1970;Nestola et al., 2011), and by the approach we implement in this study, Raman Spectroscopy (e.g., Ashley et al., 2014;Enami et al., 2007). The residual inclusion pressure can be combined with elastic modeling to calculate the initial entrapment pressure (P trap ) or temperature (T trap ) (e.g., Adams et al., 1975aAdams et al., , 1975bRosenfeld, 1969;Rosenfeld & Chase, 1961;Van der Molen & Van Roermund, 1986;Zhang, 1998).
To calculate entrapment pressures and temperatures, we perform one-dimensional elastic modeling to determine the elastic response of inclusion-host pairs upon return to ambient conditions. We implement the elastic model derived by Angel, Mazzucchelli, et al. (2017, equation A5): where (V/V trap ) inc and (V/V trap ) host are the molar volume ratios of the entrapped inclusion and the free host at measurement (V) and entrapment (V trap ) conditions, respectively. (V trap /V foot ) inc is the ratio of the inclusion molar volume at entrapment conditions and at ambient temperature along the inclusion-host pair isomeke (V foot ) (e.g., Angel, Mazzucchelli, et al., 2017). G host , and P ext are the shear modulus of the host and ambient pressure (0.1 MPa), respectively. P inc is the residual pressure of the inclusion at measurement conditions. To calculate entrapment conditions with Equation 1, a unique solution is identified by iteratively calculating V trap until both sides of the equation are equal (e.g., Ashley et al., 2016;Angel, Mazzucchelli, et al., 2017;Guiraud & Powell, 2006;Kohn, 2014). V trap is a function of both the pressure and temperature of entrapment. Therefore, for many inclusions-host pairs it is necessary to constrain the temperature or pressure of entrapment so that the second unknown variable can be accurately determined. Equation 1 can also be solved to determine how the residual P inc changes across P-T space. The points in P-T space that represent solutions to a constant P inc have been referred to as isomekes ("iso" = equal, meke = "length") (Adams et al., 1975a(Adams et al., , 1975b. Isomeke lines represent points in P-T space along which fractional volume (V ) changes of an inclusion and host are equal: The instantaneous slope of an isomeke line can be defined by the equation where V, α, and β, are the molar volume, thermal expansivity, and compressibility of the inclusion and host, respectively (Angel et al., 2015) ( Table 1). The instantaneous slope is useful for evaluating if inclusion-host pairs are suitable barometers or thermometers and can also be used to constrain the curvature of isomekes. Isomekes often occur as nearly straight lines at geologically relevant conditions but sometimes exhibit curvature due to nonlinear responses of some minerals to changing P-T conditions (e.g., minerals that undergo phase transitions); isomekes should be carefully evaluated when working with minerals that potentially undergo phase transitions (e.g., titanite, quartz, calcite, and ilmenite). For example, quartz-in-host isomekes become increasingly nonlinear as temperatures approach the quartz α-β transition. Furthermore, at absolute zero (T = 0 K) thermal expansivities are zero, and thus, the isomeke slope is zero ([∂P/∂T] isomeke = 0), but the slope is nonzero above this temperature (Angel et al., 2015;Rosenfeld & Chase, 1961).
Isomeke slope and spacing are important criteria for elastic thermobarometry (e.g., Angel et al., 2015;Kohn, 2014). Inclusion-host pairs with steep isomeke slopes make the best thermometers (high [∂P/∂ T] isomeke ; e.g., zircon inclusions in garnet), whereas those with shallow slopes (low [∂P/∂T] isomeke ; e.g., quartz inclusions in garnet) are the best barometers. Intermediate isomeke slopes allow either P trap or T trap to be determined but require a good independent constraint of one of the variables. Isomeke spacing indicates the sensitivity of a thermobarometer to a change in a residual inclusion pressure, that is, the amount that P trap or T trap changes with a change in P inc . To calculate the slope and the spacing between isomekes, we assume that isomeke lines are linear. To quantify if a linear approximation is appropriate at geologically relevant conditions, we estimate the curvature (κ) of isomekes by calculating the difference in the angle (ϴ) of an isomeke slope (dP/dT, relative to the x axis) across two T trap windows (T trap = 500-700°C; 700-900°C). All isomeke curvature, slope, and spacing results are provided in the supporting information (Tables S1 and S2).

Molar Volume Calculations and Thermodynamic Properties
To perform molar volume calculations, we used thermodynamic properties from various sources (sources provided in the supporting information). Most of the thermodynamic properties were derived from the Holland and Powell (2011) database. We use the modified Tait Equation of State (EoS) and thermal pressure term to solve for molar volumes at given P-T conditions (Holland & Powell, 2011). We implement Landau and Braggs-William theory to calculate excess molar volumes that are associated with spontaneous strain development of phases that undergo phase transitions (e.g., Holland & Powell, 1996a, 1996b, 1998. To model quartz volumes, we use the thermodynamic properties and the curved-boundary modeling approach from Angel, Alvoro, et al. (2017). Because we use apatite, epidote, garnet, and rutile in our case studies, we use updated thermodynamic properties to model the molar volumes of these phases. To model apatite and epidote molar volumes, we use the thermodynamic properties given in Ashley et al. (2017) and Cisneros et al. (2020), respectively, which are derived from previous P-V-T experiments (apatite: Brunet et al., 1999;Hovis et al., 2014Hovis et al., , 2015Schouwink et al., 2010;epidote: Gatta et al., 2011;Pawley et al., 1996;Qin et al., 2016) and the Tait EoS with a thermal pressure term. To model grossular garnet, we use the P-V-T experiments and thermodynamic properties of Milani et al. (2017) and a Tait EoS with a thermal pressure term. To calculate rutile molar volumes, we use the revised thermodynamic properties and modeling approach of Zaffiro et al. (2019), which implements a Birch-Murnaghan third-order EoS and a Kroll thermal model. We account for solid solutions of apatite (OH-Cl-F) inclusions and garnet, olivine (Fe-Mg), and epidote-clinozosoite (Al-Fe) hosts by implementing ideal (linear) mixing of end-member molar volumes. This approximation has been shown to be appropriate for apatite (Hovis et al., 2015), olivine (Speziale et al., 2004), and epidote (Cisneros et al., 2020;Franz & Liebscher, 2004). Variations in molar volume between almandine-pyrope garnet solid solutions exhibit ideal behavior (Milani et al., 2015). Variations in molar volume between almandine-grossular garnet solid solutions exhibit slight nonideal behavior Geochemistry, Geophysics, Geosystems (Cressey et al., 1978), and pyrope-grossular solid solutions exhibit greater nonideal behavior (Geiger, 2000); thus, our linear approximation may introduce additional errors in our final P trap calculations. However, the difference in molar volume between pure almandine and grossular end-members far exceeds the difference between ideal and nonideal molar volume estimations, and not accounting for mixing introduces greater uncertainties. A comparison of P trap calculated from EoSFit-Pinc  and our script is provided in the supporting information (Table S3); this accounts for reproducibility of molar volume and elastic modeling calculations.
For minerals that have experimentally determined shear moduli, we use those shear moduli in our calculations. Otherwise, the shear modulus of the host is calculated from mineral Poisson ratios and bulk moduli, following the equation where G host , K host , and ν host are the shear modulus, bulk modulus, and Poisson ratio of the host phase. Shear moduli and Poisson ratios were compiled from the literature and are provided in the supporting information. We provide further reference information on the host shear moduli that we use in our case studies. For fayalite and forsterite, we use the experimentally determined shear moduli from Speziale et al. (2004) (G fa = 51.2 GPa; X fa = 0.94) and Zha et al. (1996) Ryzhova et al., 1966); however, we are not aware of experimental shear modulus data for clinozoisite.
To maintain consistency, we calculate the shear moduli of epidote group minerals by using Equation 7, with the aggregate Poisson's ratio of 0.26, from Mao et al. (2007) and the regressed bulk moduli from Cisneros et al. (2020) (G ep = 65.7 GPa, G cz = 84.0 GPa, and G zo = 59.7 GPa). We implement ideal mixing of shear moduli to account for solid solutions of host phases. We limit our inclusion-host pair calculations to 72 phases with reported shear moduli or Poisson ratios (Table S1).

Analytical Techniques
For each case study, we petrographically identified hosts containing mineral inclusions. We identified apatite inclusions in olivine and rutile inclusions in garnet in grain mounts of handpicked crystals. Albite inclusions in epidote were identified in~80 μm-thick petrographic sections. Inclusions were isolated and two-tothree times the inclusion radial distance from the host exterior, fractures, cleavage, or other inclusions (e.g., Campomenosi et al., 2018;Mazzucchelli et al., 2018;Zhang, 1998;Zhong et al., 2020).
We performed analyses with Raman systems at Baylor University and Virginia Tech. Measurements at Baylor University (ap-in-ol and ru-in-grt) were carried out on a Thermoscientific DXR Raman microscope with 1,600 lines mm −1 grating. Measurements used a 532 nm laser with 10 mW laser power at the sample surface and a 100X objective, which allowed a spatial resolution of 1 μm. The 1,000 cm −1 band of polystyrene was used for standardization. Measurements at Virginia Tech (ab-in-ep) were carried out on a JY Horiba LabRam HR800 Raman Spectrometer with 1,800 lines mm −1 grating, and a 100X objective with a 0.9 numerical aperture. A 514 nm Ar laser was used for analyses, and the 520.3 cm −1 plasmaline was used to simultaneously correct feldspar measurements for drift by applying a linear drift correction (Table 2; Figure 3). Based on the long-term repeat measurements of the~291, 479, and 507 cm −1 bands of Amelia Albite, peak position uncertainty at Virginia Tech is ±0.1 cm −1 . No thermal corrections were applied to the measured Raman spectra. A linear background subtraction was simultaneously applied to the spectra during peak fitting of inclusion and host Raman bands, and the Ar plasmaline (only VT analyses). We used a Voigt Area model within PeakFit v4.12 (by Systat Software Inc) to fit the~964 cm −1 , and~291 cm −1 ,~479 cm −1 and 507 cm −1 bands of apatite, and albite, respectively, and a Pearson IV model to fit the~444 cm −1 and 609 cm −1 bands of rutile; a Gaussian area model was used to fit the 520.3 cm −1 Ar plasmaline ( Figure 3). Quartz inclusions in garnet (qtz-in-grt) were also measured by Raman spectroscopy to compare with results from the ab-in-ep barometer. Further information about peak fitting inclusion and host bands and details about qtz-in-grt measurements and calculations are provided in the supporting information (Tables S5  and S6). Raman uncertainties that are reported propagate errors of simultaneous background subtraction and peak fit statistics of inclusion-host Raman bands, Ar plasmaline fit statistics (only Virginia Tech analyses), and instrument uncertainty (±0.1 cm −1 at VT and ±0.3 cm −1 at Baylor). P inc errors propagate Raman uncertainties and errors of constants from polynomial fits to data from hydrostatic experiments ( Figure S4; Table S4).
Fully encapsulated inclusions preserve strain that causes the Raman active vibrational modes of inclusions to be shifted to higher or lower wave numbers relative to minerals that are unstrained (fully exposed). We calculate the Raman shift(s) of inclusions (ω inc ) relative to Raman shift(s) of individual separated grains (apatite and rutile, ω ref ), and an unencapsulated standard (amelia albite; ω ref ) at ambient conditions (Δω = ω inc -ω ref ). To calculate residual inclusion pressures (P inc ) of minerals, we use pressure-dependent Raman shift(s) (P-Δω) relationships that have been calibrated under hydrostatic stress conditions by using hydrothermal diamond anvil cell experiments (HDAC), for example, apatite (Comodi et al., 2001;Schouwink et al., 2010); feldspar: (Aliatis et al., 2017;Befus et al., 2018); anatase/rutile: (Liu & Mernagh, 1992;Samara & Peercy, 1973;Swamy et al., 2005); and quartz /coesite: (Asell & Nicol, 1968;Hemley, 2013;Schmidt & Ziemann, 2000). The pressure-dependent Raman shift(s) of vibrational modes of many other minerals remain unconstrained. Recent work has explored the effects of nonhydrostatic stress on Raman active vibrational modes of minerals by using density functional theory strain-Δω simulations to determine Grüneisen tensors that can be used to convert mineral Raman shifts (Δω) to strains (e.g., Murri et al., 2018;Nestola et al., 2018;Stangarone et al., 2019). However, we use hydrostatic P-Δω calibrations to calculate residual P inc because no Grüneisen tensors have been determined for apatite, rutile, and feldspar. Some work suggests that apatite inclusions in garnet develop minimal anisotropy when heated at ambient pressure ; therefore, not accounting for anisotropic strain of apatite inclusions may introduce minimal errors. This estimation may be erroneous for rutile and albite inclusions.

Inclusion-Host Pair Modeling
We test the pressure-temperature sensitivity of 5,142 inclusion-host pairs by carrying out entrapment P-T calculations to produce isomekes in P-T space ( Figure 1). We calculated the slope and spacing of the resulting isomekes by assuming they are linear and parallel. We acknowledge that our assumption may misrepresent isomekes that are nonlinear (e.g., calcite or quartz inclusions produce curved isomekes). However, the majority of inclusion-host pairs produce near linear isomekes with curvature (κ) less than 2°at geologically relevant conditions (500-900°C); a curvature ≤2°is equal to a change in slope of ≤0.035 MPa°C −1 . Of the >5,000 pairs we analyzed, 75% have isomekes that curve less than 2°, and 86% curve less than 5°across the 400°C window (Table S2). Treatment of these isomekes as linear functions is an appropriate approximation, albeit,    Tables S1 and S2). The subsequent results we present serve as a first-order guide, but users should carefully examine the specific isomeke slope and spacing of a target inclusion-host pair.
Isomekes were constructed by solving for P trap in Equation 1, using P inc = −200 to 800 MPa and T trap = 300°C to 1,000°C, with 100 MPa and 100°C intervals, respectively. Isomeke slopes vary from 0.04 to 10 5.1 MPa°C −1 , with a median value of 1.8 MPa°C −1 (Figures 1a,  1c, and 1e; Table S1). Pairs at the 16th and 84th percentiles have isomekes slopes that are ≤0.5 and 5.6 MPa°C −1 , respectively (Figure 1a). At the median, an isomeke slope of 1.8 MPa°C −1 indicates that if the T trap estimate changes by 100°C, the calculated final P trap changes by 0.18 GPa (e.g., quartz-in-almandine; Figure 1a).
We measured isomeke spacing perpendicular to isomeke lines (Equation S2); hence, the units may be considered to be between MPa (horizontal isomekes) and°C (vertical isomekes); we do not give an exact spacing unit. Isomeke spacing was calculated by measuring the distance between the 100 and 200 MPa isomekes at T trap = 500°C. Isomeke spacing ranges from~0 up to 10 4.3 , with a median of 193. The 16th and 84th percentiles of inclusion-host pairs having isomeke spacings that are ≤92 and 383, respectively ( Figure 1d). Isomeke spacing of a single inclusion-host pair can vary between different isomekes and in P-T space, but we use our spacing calculation as a representative estimate that allows us to compare spacing with other inclusion-host pairs. We refer to sensitivity in terms of isomeke spacing. The relative terms "high sensitivity" and "low sensitivity" refer to smaller and greater isomeke spacing, respectively. Small isomeke spacing (high sensitivity) indicates that a large change in P inc causes a small change in P trap or T trap (Figure 1b). Tightly spaced isomekes reduce uncertainty in the final P trap or T trap , because the error associated with the residual P inc produces a smaller change in P trap or T trap . We note that our description of sensitivity does not refer to the sensitivity of Δω to P (P-Δω relationships from hydrostatic calibrations) and is independent of this relationship.

Case Studies
We present introductory case studies for three inclusion-host pairs from different geologic environments and discuss the Raman bands that correspond to specific vibrational modes (ν x ) that we focused on analyzing to calculate residual P inc . We present few analyses for the studied inclusion-host pairs, and we emphasize that these calculations only serve as illustrations and do not quantify the applicability of the studied pairs.

Apatite-in-Olivine: An Igneous Rock Application
Volcanologists and igneous petrologists use the P-T conditions of magmatic rocks to better understand melt genesis, assimilation and fractional crystallization, eruption triggers, and melt transport. Magma composition controls the mineral phases that may crystallize, and the diversity of magma compositions permits many potential inclusion-host pairs to be considered as elastic thermobarometers.

10.1029/2020GC009231
We first examine apatite inclusions hosted by olivine phenocrysts (ap-in-ol) from the rhyolitic obsidian Solfatara Plateau Lava flow in Yellowstone caldera (Figures 2a, 2d, and 2g). Apatite is rarely used as an inclusion, but previous studies have shown that apatite inclusions in garnet can record reasonable conditions of garnet growth in skarn deposits (Barkoff et al., , 2019. Elastic modeling shows that for ap-in-ol, over the window of T trap = 500°C to 900°C and P inc = 200 MPa: κ = 0.7°, dP/dT ≈ 1.7 MPa°C −1 (Figures 1a and 1c), and between the −200 and −100 MPa P inc isomekes at T trap = 500°C, spacing ≈280 (Figure 1d). The isomeke ). Pairs that fall in between these end-members require independent constraints to calculate the unknown variable of interest. (b) Spacing between isomekes is a metric of sensitivity. Wide spacing between isomekes indicates that a small change in P inc will produce large changes in P trap or T trap . Small spacing between isomekes indicates that a small change in P inc will produce small changes in P trap or T trap . We note that spacing calculations are unit-dependent, calculations using GPa units would results in different absolute values, but we are only concerned with relative differences. (c-f ) Cumulative frequency and histogram plots indicating the proportion of inclusion-host pairs that serve as barometers or thermometers and their sensitivities. Most inclusion-host pairs serve as thermobarometers that require independent entrapment pressure or temperature estimates to derive the unknown variable (entrapment P or T).

10.1029/2020GC009231
slope indicates that the ap-in-ol thermobarometer is near the median of all tested pairs and that it could be appropriate as a barometer. The wide isomeke spacing indicates that this is a low sensitivity thermobarometer and that small changes in P inc will cause large changes in calculated P trap (Figure 1b).
For apatite in the P6 3 /m space group, group theory predicts the following Raman active vibrational modes (Kroumova et al., 2003): A total of 54 vibrational modes are Raman active, but only a few modes have sufficient intensity to be observed in the Raman spectra of randomly oriented apatite (Comodi et al., 2001). We focus on thẽ 964 cm −1 band of apatite inclusions because of the strong intensity of this band (Figures 2a, 3a, and Isomekes graphs for (a, d, and g) apatite-in-olivine (X ap-fluoro = 1, X fa = 1), (b, e, and h) rutile-in-garnet (X ru = 1, X gr = 1), and (c, f, and i) albite-inepidote (X ab = 1, X ep,cz = 0.5). (a) Apatite-in-olivine isomekes show that the thermobarometer has a moderate T dependence, wide isomeke spacing (changing P inc by 100 MPa changes P trap by~1 GPa), and apatite inclusions primarily preserve compression at geologic conditions (tension at low P trap ), (d) P trap has a significant X fa dependence at high P inc (T trap = 500°C), and (g) T trap has a moderate X fa dependence (P trap = 1 GPa). (b) Rutile-in-garnet isomekes show that the thermobarometer has a moderate T dependence, wide isomeke spacing (changing P inc by 100 MPa changes P trap by~1 GPa), and rutile inclusions should preserve tension at geologic conditions; (e) P trap has a significant X gr dependence at low P inc (T trap = 500°C); and (h) T trap has a moderate X gr dependence (P trap = 1 GPa). (c) Albite-in-epidote isomekes show that the barometer has a minimal T dependence, narrow isomeke spacing (changing P inc by 100 MPa changes P trap by~0.25 GPa) and albite inclusions should preserve compression at geologic conditions; (f) P trap has a moderate X ep dependence at high P inc , but its dependence on composition is minimal at low P trap (T trap = 500°C); and (i) because of the shallow isomekes of this barometer, T trap shows an extreme X ep dependence (P trap = 1 GPa). Black solid lines indicate where P inc > 0 (thin) and P inc = 0 (thick). Dashed lines indicate P inc < 0 isomekes. Gray areas highlight the P-T space where measured residual P inc should retain tension (P inc < 0). 4a). Raman intensities of other apatite bands were too low to peak fit. The prominent~964 cm −1 band represents the phosphate symmetric stretch (ν 1 [A g ]; Comodi et al., 2001). No overlap was observed between apatite and olivine Raman bands, and the Raman spectra of apatite inclusions required no deconvolution of shoulder bands near 964 cm −1 (Figure 3a). Neighboring apatite peaks at~950 and 1,030 cm −1 represent distortion of the symmetric stretching mode owing to radiation damage, and the asymmetric stretch of the ν 3b (E 2g ) mode, respectively (Comodi et al., 2001;J. Liu et al., 2008). We postulate that the measured apatite inclusions have negligible radiation damage and are appropriate to use for this study, because of the relatively young age (<~100 ka) of the Solfatara Plateau Lava flow and the low intensity of the~950 cm −1 band (Figure 3a). Separated apatite grains have mean peak positions at 962.3 ± 0.1 cm −1 (n = 4, s.d. of measurements), whereas apatite inclusions in olivine have Geochemistry, Geophysics, Geosystems mean peaks at 961.7 ± 0.4 cm −1 (n = 10) ( Table 2). The enclosed apatite inclusions primarily exist under tension, as demonstrated by the mean negative peak shift of apatite inclusions relative to apatite in the matrix (Δω = −0.6 ± 0.4 cm −1 ).
The apatite 964 cm −1 band exhibits a pressure-dependent Raman shift that has been experimentally calibrated under hydrostatic conditions (Comodi et al., 2001;Schouwink et al., 2010). We calculate residual pressures by using the P-Δω relationship for fluoroapatite given in Ashley et al. (2017). Additional error may be introduced by using a P-Δω relationship derived from Durango fluoroapatite, but the difference in composition between the Durango apatite used for hydrostatic experiments (X F = 0.99) and our measured apatites (X F = 0.91) is minimal (Schouwink et al., 2010). The calculated mean residual P inc and s.d. recorded by apatite is −141 ± 96 MPa (n = 10), suggesting that these inclusions preserve tension (Table 2).

Rutile-in-Garnet: A Metamorphic Application
Metamorphic petrology is fundamental for understanding crustal and mantle dynamics, with processes that include continent and oceanic lithosphere subduction, obduction and exhumation, crustal strength, fault mechanics and rheology, and the effects of mantle dynamics on the evolution of the crust. A key approach used for understanding these problems is constraining the P-T evolution of metamorphic rocks from these different environments. A large number of metamorphic rocks have suitable thermobarometers; nonetheless, P-T conditions are difficult to extract from some common rock types (e.g., retrogressed metamorphic rocks, eclogites, skarn deposits, and spinel-bearing peridotites). Elastic thermobarometry offers great promise for metamorphic rocks and can potentially be used to constrain P-T conditions of distinct mineral growth.
Elastic modeling shows that for ru-in-grt, over the window of T trap = 500°C to 900°C and P inc = 200 MPa: κ = 5°, dP/dT ≈ − 3.8 MPa°C −1 (Figures 1a and 1c), and between the −400 and −300 MPa P inc isomekes at T trap = 500°C, spacing ≈310 (Figure 1d). The isomeke slope indicates that the ru-in-grt thermobarometer is above the median dP/dT of all tested pairs and that it requires a good independent temperature constraint. Furthermore, the wide isomeke spacing indicates that this is a low sensitivity thermobarometer and that small changes in P inc will cause large changes in P trap (Figure 1b).
For rutile in the P4 2 /mnm space group, group theory predicts the following Raman active vibrational modes (Kroumova et al., 2003): The E g and A1 g vibrational modes can be well observed in the Raman spectra as broad bands at~444 and 609 cm −1 , respectively (Figures 3b; e.g., Balachandran & Eror, 1982). Additional Raman bands at~145 and~826 cm −1 have been assigned to B 1g and B 2g modes, respectively (Frank et al., 2012), and a broad

10.1029/2020GC009231
Geochemistry, Geophysics, Geosystems band observed at~235 cm −1 has not been predicted by theoretical calculations. We focus on the intensẽ 444 and~609 cm −1 bands (Figure 3b). No overlap was observed between rutile and garnet Raman bands, but the Raman spectra of rutile inclusions required deconvolution of~2-3 shoulder bands near both the 444 and 609 cm −1 bands, suggesting that these bands result from multiphonon modes (see supporting information). Separated rutile grains have mean peak positions at 443.3 ± 0.7 and 609.2 ± 0.3 cm −1 (n = 5), whereas fully encapsulated rutile inclusions have average peak positions at 444.2 ± 0.6 and 609.5 ± 0.3 cm −1 (n = 12) ( Table 2). We use the hydrostatic experiments of Liu and Mernagh (1992) that demonstrate that both rutile bands experience pressure-dependent Raman shifts to higher wave numbers and regress their rutile data to determine P-Δω relationships ( Figure S4; Table S4). The mean P inc and s.d. calculated from the~444 and~609 cm −1 rutile bands are 135 ± 92 and 59 ± 77 MPa (n = 12), respectively ( Figure 4).

Feldspar-in-Epidote: A Second Metamorphic Application
Our final case study examines the potential of albite inclusions in epidote (ab-in-ep). The high β (low K) of feldspar and its widespread occurrence in diverse geologic settings make it a promising inclusion that can be used as barometer in many host phases (Figures 1 and 2c). We focus on testing the ab-in-ep barometer because epidote is a "stiff" host and is abundant in many mafic and felsic metamorphic rocks. Albite and epidote are also common in a wide variety of hydrothermal systems, and the ab-in-ep barometer can be useful for constraining pressures of deep roots below porphyry deposits, where both minerals are commonly found as alteration products (e.g., Runyon et al., 2019).
We analyze a mafic blueschist from the high-P/low-T Cycladic Blueschist Unit (Kalamisia Beach: 37°25′ 12.2″N, 024°57′23.9″E) on Syros, Greece (Figures 2c, 3c, and 4c) to test the ab-in-ep barometer. The metamorphic complex exposed on Syros reached max P-T conditions of~1.5-2.2 GPa and~500-550°C at~50 Ma (e.g., Laurent et al., 2018;Schumacher et al., 2008). It was subsequently exhumed to the midcrust by the Miocene with episodes of deformation and retrogression during exhumation prior to and post-core complex capture (Bröcker & Enders, 1999;Putlitz et al., 2005;Soukis & Stockli, 2013). Our sample is composed of an assemblage of garnet + omphacite + glaucophane + epidote + phengite + albite + quartz + chlorite ( Figure S6). The primary external foliation is defined by omphacite, glaucophane, epidote, and white mica. Garnets are rotated by the primary foliation and often contain glaucophane growth within pressure shadows and brittle garnet fractures but do not have a well-developed internal foliation, suggesting that the primary foliation formed after garnet growth. Omphacite displays alteration and breakdown to glaucophane, and glaucophane and other inclusions within epidote are commonly oriented parallel to the external foliation. No omphacite is observed as inclusions within epidote. Therefore, our observations suggest that epidote crystallization occurred during retrograde metamorphism and garnet grew during earlier high-P metamorphism. Elastic modeling shows that for ab-in-ep, over the window of T trap = 500 to 900°C and P inc = 200 MPa: κ = 0.2°, dP/dT ≈ − 0.05 MPa°C −1 (Figures 1a and 1c), and between the 100 and 200 MPa P inc isomekes at T trap = 500°C, spacing ≈ 250 (Figure 1d). The isomeke slope indicates that the ab-in-ep barometer is well below the 16th percentile of all tested pairs (0.5 MPa°C −1 ) and that theoretically, it is an excellent barometer. The barometer has isomeke spacing above the median of all pairs (193) but is more sensitive than the ap-in-ol and ru-in-grt thermobarometers (Figures 1b and 2c).

Discussion
Our results suggest that many untested inclusion-host pairs from diverse geologic environments have the potential to be used in elastic thermobarometry (Table 3). To test some of the newly modeled inclusion-host pairs, we provide three case studies that use samples with reference magmatic and metamorphic P-T constraints. We show that mechanical, theoretical, and analytical limitations, and/or inaccurate thermodynamic properties, may sometimes restrict the utility of elastic thermobarometry.

Isomeke Slope and Spacing
Previous work has shown that isomeke slopes are related to the difference between the K and α of an inclusion-host pair (Equation 3; Angel et al., 2015;Kohn, 2014). We similarly find that inclusions with high β (e.g., low K phases such as quartz, feldspar, apatites, amphiboles, graphite, and micas) serve as the best barometers, and inclusions-host pairs with large K differences are the best barometers. Inclusion-host pairs that have similar K, but a significant difference in α, serve as the best thermometers (e.g., zircon-in-garnet, zircon-in-rutile, and garnet-in-olivine). Similar to elastic thermometers, isomeke spacing is related to the difference of K between an inclusion and host. Inclusion-host pairs with small and large K differences will have large and small isomeke spacing, respectively.
Isomeke slope and spacing demonstrate the P-T dependence of the inclusion-host pair, and the sensitivity of the pair to changes in P inc , respectively. The useful range of each variable will depend on the application of interest. For example, quartz-in-almandine isomekes have a slope of~1.8 MPa°C −1 , near the median slope of 1.9 MPa°C −1 , and far above the arbitrary 16th percentile slope threshold of 0.5 MPa°C −1 (Figure 1c); however, quartz-in-garnet is known to be a reliable barometer (e.g., Ashley et al., 2014;Bonazzi et al., 2019;Thomas & Spear, 2018;Wolfe & Spear, 2018). Large isomeke spacing (≥~300) may make using some inclusion-host pairs (e.g., rutile-in-garnet) to constrain P-T conditions of crustal processes challenging. Inclusion-host pairs can also be sensitive to mineral and host compositions, affecting their theromobarometry potential. Table 3 demonstrates how different solid-solution compositions can change the isomeke slope and spacing; thus, the compositional sensitivity of inclusion-host pairs should be evaluated for reasonable geologic compositions.

Current Limitations of Elastic Thermobarometry and Calculations in This Study
We emphasize that our theoretical calculations (isomeke slope and spacing) and associated inclusion-host pair case study calculations (sections 4.3-4.5) have several limitations. Isomeke slope and spacing should be carefully evaluated for each inclusion-host pair of interest. At geologically relevant temperatures, nonlinear isomekes may become increasingly important for minerals that undergo phase transitions, and this may also limit the applicability of our spacing calculations that assume linear isomekes. We describe two other major limitations of our case-study calculations in this work: strain resolvability and inclusion-host anisotropy. We refer to strain resolvability as uncertainties associated with P-V-T data, and the Δω precision required to provide reasonable P trap errors.
Elastic thermobarometry requires careful evaluation of mineral thermodynamic properties and may require refitting P-V-T data with appropriate EoS. Implementing thermodynamic properties from large thermodynamic databases (e.g., Holland & Powell, 2011) should be approached with caution. In this work, we implement refit P-V-T data for apatite, rutile, garnet, and epidote; however, olivine and albite thermodynamic properties are derived from the Holland and Powell (2011) database and may require further evaluation. As further discussed in sections 4.3 and 4.4, for inclusion-host pairs such as apatite-in-olivine and rutilein-garnet, Δω errors must be small because these inclusion-host pairs have large isomeke spacing. Δω errors will increase P inc errors (calculated using Grüneisen tensors or hydrostatic calibrations) and resultant P trap 10.1029/2020GC009231 Geochemistry, Geophysics, Geosystems

Note.
Ranges are provided for mineral pairs to encompass compositional endmembers. Results for specific solid solutions are shown for some realistic pairs. Isomeke slope and spacing provide insight into a thermobarometers potential and sensitivity. P inc range provides the theoretical range of pressures an inclusion should preserve at ambient conditions, if it records reasonable geologic P-T conditions. a Spacing is measured perpendicular to the slope of the isomeke, and measured as the difference between 100 MPa isomekes. b It is unclear if high magnitudes of tension can be preserved in inclusions. Values as high as these may be unlikely. c Apatite is a common inclusion, generally exists under tension at ambient conditions, and is typically a thermometer. Radioactivity may soften apatite and/or zircon. d See Izraeli et al. (1999) for further applications of inclusions in diamond.

10.1029/2020GC009231
errors. Δω errors can be reduced by using high spectral resolution Raman systems, appropriate spectral calibrations, and many repeat measurements.
The effect of inclusion and host mineral anisotropy (different compressibility and thermal expansivity along different crystallographic axes) on elastic thermobarometry calculations is an issue that has been addressed by recent studies (e.g., Angel et al., 2019;Murri et al., 2018). Recent work has shown that using hydrostatic stress calibrations to constrain quartz inclusion pressures may sometimes lead to inaccurate estimates of garnet growth conditions (3.0 GPa experiments); however, at lower pressures (2.5 GPa experiments), hydrostatic calibrations replicate known pressure conditions of garnet growth (Bonazzi et al., 2019). Furthermore, previous experimental studies suggest that calculating quartz P incl by using hydrostatic calibrations produces accurate P trap estimates of garnet growth (Thomas & Spear, 2018). The current limitations of using hydrostatic calibrations to calculate P trap remain unclear and certainly require further experimental testing. We reiterate that our case studies use hydrostatic calibrations to constrain residual pressures of apatite, rutile, and albite inclusions and that we do not account for anisotropy. Therefore, the inclusion-host pairs in our case studies require further evaluation before any pairs can be implemented for accurate P trap calculations. Furthermore, the accuracy of modeling anisotropic minerals with an isotropic elastic model remains unknown, though recent results suggest that an isotropic elastic model satisfactorily simulates the P inc evolution of an anisotropic inclusion-host pair during isobaric heating (quartz inclusions in epidote; Cisneros et al., 2020).

Apatite-in-Olivine: An Igneous Rock Application
Apatite inclusions from the Solfatara Plateau, Yellowstone caldera, preserve tension, with a mean residual P inc = −141 ± 96 MPa (n = 10, Table 2), calculated from the apatite 964 cm −1 band. The s. d. around the mean residual P inc is large and may result from inclusion and/or host anisotropy that we have not accounted for, apatite shape effects, or differences in compositions between inclusion and matrix grains. To calculate P trap , we estimate apatite compositions by using empirical volatile relationships to be X Cl = 0.04, X F = 0.91, and X OH = 0.05 (see the supporting information). The pre-eruptive magmatic temperature (T trap ) and olivine composition of Solfatara Plateau are 800 ± 50°C and Fa 92 ± 1 (Befus & Gardner, 2016;Vazquez et al., 2009).
We calculate that Solfatara Plateau apatite inclusions were entrapped by olivine at~0.44 ± 0.57 GPa using mean mineral compositions and our mean residual P inc (Figure 5a). Our best estimate P trap (~0.44 GPa) suggests crystallization in a deep crustal reservoir that is deeper than independent estimates for rhyolitic storage at Yellowstone caldera, but our uncertainty overlaps with reference conditions (Befus & Gardner, 2016;Gryger, 2017;Luttrell et al., 2013;Myers et al., 2016;Shamloo & Till, 2019). We attribute the discrepancy to significant limitations with the ap-in-ol thermobarometer. P trap estimates from the thermobarometer vary considerably as a function of olivine composition; this effect becomes more pronounced at higher residual P inc (e.g., 800 MPa; Figure 2d). At lower residual P inc (e.g., −100 MPa), the effect of olivine composition on the resultant P trap is smaller, but nonnegligible. Given our estimated isomeke slope of the ap-in-ol thermobarometer (dP/dT ≈ 1.7 MPa°C −1 ), changes in T trap should have an effect on the final Figure 5. Entrapment pressure estimates from two case studies. (a) Apatite-inolivine results (yellow star) suggest a higher P trap than previously published reference P-T conditions (gray box); however, uncertainties are large for these calculations. The thick red dashed line is our calculated residual P inc , and thin red dashed lines represent the 1σ P inc error. Apatite and olivine modeled as X ap-f luoro = 0.91, X ap-hydroxyl = 0.05, X ap-chloro = 0.04, and X fa = 0.92, X fo = 0.08, respectively. (b) Albite-in-epidote P trap isomeke (solid red) relative to the quartz-in-garnet P trap isomeke (solid blue), previously determined maximum P-T conditions (gray box), and the jadeite + quartz to albite reaction line (solid gray). Ab-in-ep and qtz-in-ep P trap isomekes are calculated from the weighted mean P inc of all analyses (Tables 2 and S6). Dashed red, blue, and gray lines bracket the 1σ error of P trap and the reaction line. Felspar and epidote modeled as X ab = 1, and X ep = 0.5, X cz = 0.5, respectively. Quartz and garnet modeled as X qtz = 1, and X alm = 0.7, X grs = 0.2, and X prp = 0.1, respectively. P trap , for example, changing T trap by 50°C changes P trap by ± 0.1 GPa. A further limitation associated with the ap-in-ol thermobarometer is the isomeke spacing (~280). The average P inc error of our individual inclusion measurements (±65 MPa) would change the final P trap by ± 0.4 GPa (Figures 2a and 5a). Raman system uncertainties of~0.1 cm −1 can reduce P inc errors to~34 MPa (P trap ± 0.2 GPa), but even ±0.2 GPa errors may not be appropriate for many applications when compounded with other uncertainties (e.g., T trap , mineral compositions, shape effects, softening of apatite from radiation damage, and anisotropy). In principle, the ap-in-ol thermobarometer can have wide-ranging applications for constraining pressures of magmatic systems but will require accurate constraints on all input variables, low instrumental uncertainties, or many repeat measurements.

Rutile-in-Garnet: A Metamorphic Application
The mean residual P inc calculated from the rutile 444 and 609 cm −1 bands are 135 ± 92 MPa and 59 ± 77 MPa (n = 12), respectively. Individual rutile inclusion analyses provide similar P inc from both rutile bands, and most inclusion analyses center near P inc = 0 (within error; Figure 4). Using our mean residual P inc , T trap = 676°C, and a Pyp 3 Alm 40 Grs 20 Sps 10 garnet core composition, we calculate that rutile inclusions from Verpeneset eclogites were entrapped by garnet at negative pressures. Cleary, these results are unrealistic and do not represent known geologic conditions of garnet growth and reference pressure estimates. The errant result can be sourced back to the calculated residual P inc from all analyses. Rutile inclusions should preserve negative pressures (tension) at Earth conditions, and high negative residual pressures at reference conditions (~−300 MPa; Figure 2b) .
Several possible processes can explain the low P trap that the ru-in-grt thermobarometer records. During exhumation, rutile inclusions may have in-elastically relaxed next to the garnet exterior or fractures, thus reducing residual pressures; however, in accordance with Mazzucchelli et al. (2018), the rutile inclusions were small relative to the garnet host, and sufficiently far away from the garnet exterior and fractures to avoid relaxation. Proximity to adjacent inclusions would produce higher inclusion pressures relative to those predicted by pure elastic relaxation; therefore, this possibility is also unlikely. Rutile inclusions are generally elongate and anisotropic crystals; therefore, elastic modeling of rutile inclusions as isotropic, spherical minerals may be inappropriate. Furthermore, radiation damage may lead to significant softening of rutile (decrease of K), and an increase of P inc , as we observe in our samples (lower P trap ; e.g., Beirau et al., 2016). Significant radiation damage is plausible, given the~400 Ma age of eclogite facies metamorphism in the Western Gneiss Region (Carswell et al., 2003). Additionally, tensile failure of the host may cause collapse around rutile inclusions, and an increase of P inc .
The preceding issues are plausible, especially for rutile inclusions that record positive residual P inc , but two of the primary issues that may limit ru-in-grt and other thermobarometers are the nature of rutile Raman bands, and that inclusions may not be able to retain ideal tensile strains at geologically relevant conditions. First, we find that peak positions of rutile Raman bands are extremely sensitive to the peak fitting approach that is implemented. The rutile Raman bands are broad, and the peak center position is challenging to constrain. Rutile Raman bands also exhibit varying degrees of asymmetry (especially the~444 cm −1 band), and we find that symmetric peak fitting functions did not adequately fit rutile bands ( Figure S2). Variable peak fitting approaches (e.g., symmetric vs. nonsymmetric distributions) can produce changes in peak center positions near~2 cm −1 . Second, isomekes for this thermobarometer are inverted; that is, lower residual pressures indicate higher pressures of entrapment (Figure 2b). If rutile inclusions are entrapped at high pressures, they may not be able to retain sufficient adhesion to the surrounding host cavity wall and therefore do not preserve an ideal tensional strain state. Similarly, apatite inclusionsinolivine may preserve a high P inc (~−140 MPa) relative to their ideal P inc (~−200 MPa) and may also suggest a tensile strain limit to elastic thermobarometry. For comparison, quartz inclusions in garnet have been shown to retain some tension (Ashley et al., 2015) and in some cases may preserve sufficient tension to give accurate pressures estimates (0.8 GPa experiments; P inc = −300 MPa; Thomas & Spear, 2018). The difference in result between previous studies and ours is problematic, given that the residual P inc required from both thermobarometers to replicate reference or experimental conditions is similar (P inc ≈ −300 MPa). The higher K of rutile relative to quartz (K ru = 205.1 GPa; K qtz = 64.3 GPa) may limit the amount of tension that rutile inclusions can accommodate and does not allow rutile to "stretch" to its ideal strain state. If so, this would imply that inclusions that should retain negative residual P inc at geologically relevant conditions but have a high K inlcusion may not efficiently attach to the host cavity and thus limits their applicability as elastic thermobarometers.
The ru-in-grt thermobarometer is further complicated by issues that are illustrated by associated isomekes. P trap estimates from the thermobarometer can vary considerably as a function of garnet composition, and this effect becomes more pronounced at lower residual P inc (e.g., −800 MPa; Figure 2e). Given our estimated isomeke slope of the ru-in-grt thermobarometer (dP/dT ≈ − 3.8 MPa°C −1 ), changes in T trap should affect the final P trap , that is, changing T trap by 50°C changes P trap by ±~0.2 GPa. An additional limitation with the rutile-in-garnet barometer is isomeke spacing (~310). Assuming a garnet composition of X gr = 1, T trap = 600°C, and P inc = −300 MPa (P trap = 1.72 GPa), the average P inc error of our individual inclusion measurements (±88 MPa) would change the final P trap by ±~1.1 GPa (Figure 2b). Cleary, future application of the rutile-in-garnet thermobarometer may be challenging, and successful application of this thermobarometer would require low uncertainty Raman measurements, and a better understanding of the tensile strain limit that inclusions can preserve.

Feldspar-in-Epidote: A Second Metamorphic Application
In comparison to the previous examples, we show an example that records a large, complicated range of residual P inc , but reasonable P trap conditions. The difference in P inc calculated from the~291 cm −1 (ν 14 ), 479 cm −1 (ν 22 ), and~507 cm −1 (ν 24 ) feldspar bands suggests that our measured inclusions record significant anisotropy ( Figure 4). P inc calculated from the albite 291 and 507 cm −1 bands generally record a higher P inc than P inc calculated from the albite 479 cm −1 band. The~291 cm −1 band of albite generally records the highest P inc . P inc calculated from the~479 cm −1 is generally about 1.5-2 times less than P inc calculated from thẽ 507 cm −1 band ( Figure 4; Table 2). The relationship may reflect the compression along two opposing feldspar axes, that is, tetrahedral ring compression in the ab-plane (ν 22 ) and compression in the c-plane (ν 24 ). To attempt to account for the significant P inc variability that feldspars record, we calculate a weighted mean P inc for each analysis. We assign a weight (w) of where σ i is the P inc error calculated for each albite Raman band and w i is the weight assigned to P inc calculated from each band. We calculate a mean P inc error based on the weighted variance of P inc calculated from each band, relative to the mean P inc of each analysis. This applies a lower mean P inc error to samples that have a lower P inc s.d., to apply more weight to albite analyses that exhibit less P inc variability ( Figure 4, closer to hydrostatic stress line). We calculate a weighted mean P inc for the albite analyses population based on our previous weighting, to calculate a mean P trap . We estimate the composition of the feldspar we analyzed to be pure albite (see the supporting information) and assume an epidote composition of X ep 0.5, X cz = 0.5.
To provide a comparison against P trap calculated from the ab-in-ep barometer, we carried our qtz-in-grt measurements from the same thick section. To calculate P trap from the qtz-in-grt barometer, we assume an ideal quartz and Alm 0.7 Grs 0.2 Pyp 0.1 garnet composition and account for quartz anisotropy by calculating our residual quartz P inc from strains Bonazzi et al., 2019;Murri et al., 2018). Further information about quartz P inc calculations is provided in the supporting information. Between T trap = 500°C and 550°C, using our mean P inc (574 ± 63 MPa), the qtz-in-grt barometer records P trap between 1.38 to 1.44 GPa (Figure 5b; Tables S5 and S6). At T trap = 450°C, P trap calculated from the ab-in-ep barometer is 1.3 ± 0.2 GPa (Table 2). Changing T trap from 400°C and 500°C only changes the mean P trap from 1.34 to 1.32 GPa (Figure 5b). The P trap we calculate from the qtz-in-grt and ab-in-ep barometers are consistent with high-P garnet growth near peak metamorphic conditions, and retrograde epidote growth at lower pressures. Our calculated ab-in-ep P trap is also in reasonable agreement with the reaction Jadeite + Quartz to Albite (Figure 5b) (Holland, 1980).
The higher isomeke spacing sensitivity and shallow isomeke slope of the ap-in-ep barometer relative to the ap-in-ol and ru-in-grt thermobarometers significantly reduces P trap errors; however, compounded 10.1029/2020GC009231 Geochemistry, Geophysics, Geosystems uncertainties may result from using albite P-Δω calibrations and an unstrained Amelia albite standard, if the inclusion feldspar composition does not closely approximate albite. P trap estimates also vary as a function of epidote composition, and this effect becomes more pronounced at higher residual P inc (e.g., 800 MPa; Figure 2f). Given the estimated isomeke slope of the ab-in-ep barometer (dP/dT ≈ −0.05 MPa°C −1 ), changes in T trap have a negligible effect on the final P trap . Changing T trap by up to 700°C changes P trap by less than 0.2 GPa. The ab-in-ep barometer has moderate isomeke spacing (~250). Assuming an epidote composition of X ep = 0.5, T trap = 450°C, and P inc = 400 MPa (P trap = 1.11 GPa), the population weighted mean P inc error (±76 MPa) would change the final P trap by ±~0.2 GPa (Figure 2c and 5b). This preliminary work suggests that the ab-in-ep barometer may record reasonable P trap conditions but that careful evaluation will be needed in the future to understand the effects of feldspar anisotropy and an anisotropic host. Synthesis experiments of albite inclusions that are entrapped by epidote at known conditions would be appropriate to test the accuracy of the barometer and hydrostatic calibrations (e.g., Bonazzi et al., 2019;Thomas & Spear, 2018). Furthermore, albite may be a more suitable inclusion in other host phases that exhibit less anisotropy than epidote (e.g., garnet).

Ongoing Directions
The purpose of this work is to present new inclusion-host pairs that have promising potential as elastic thermobarometers. Our contribution primarily presents the theoretical potential of inclusion-host pairs that can be used for elastic thermobarometry. We provide case studies to demonstrate inclusion-host pairs will present complexities that must be carefully considered.
We acknowledge that significant future work is needed to successfully implement the elastic thermobarometers that we present. Remaining work includes (but not limited to) better understanding the effects of anisotropy on residual P inc , tensile strain limits, the limitations of 1-dimension elastic modeling, visco-elastic effects, and the effects of nonideal geometries; however, ongoing research is making progress with many of these issues (e.g., Mazzucchelli et al., 2018;Murri et al., 2018;Stangarone et al., 2019;Zhong et al., 2020). When the elastic properties and postentrapment modifications of inclusion-host pairs are carefully considered, the technique has far-reaching petrologic potential and could allow geoscientists to constrain P-T conditions of processes in diverse mantle, metamorphic, magmatic, and extraterrestrial settings.