Articles | Volume 35, issue 4
Research article
11 Jul 2023
Research article |  | 11 Jul 2023

Elasticity of mixtures and implications for piezobarometry of mixed-phase inclusions

Ross J. Angel, Mattia L. Mazzucchelli, Kira A. Musiyachenko, Fabrizio Nestola, and Matteo Alvaro

Elastic thermobarometry (or piezobarometry) is the process of determining the P (pressure) and T (temperature) of entrapment of inclusions from their pressure, stress or strain measured when their host mineral is at room conditions. The methods and software used for piezobarometry are currently restricted to inclusions consisting of single phases. In this contribution we describe the theory of the elasticity of mixtures of different phases and combine it with the existing isotropic analysis of the elastic interactions between single-phase inclusions and their hosts to calculate the inclusion pressures of mixed-phase inclusions. The analysis shows that the reliability of calculated entrapment conditions for mixed-phase inclusions, including those containing fluid plus minerals, depends in a complex way upon the contrasts between the elastic properties of the host and the phases in the inclusion. The methods to calculate the entrapment conditions of mixed-phase inclusions have been incorporated into the EosFit7c program (version 7.6) that is available as freeware from

1 Introduction

Elastic thermobarometry (or piezobarometry) is the process of determining the P (pressure) and T (temperature) of entrapment of inclusions from the pressure measured in them when their host mineral is at room conditions. It has recently provided new constraints on the PT conditions attained during metamorphism (e.g. Zhong et al., 2018; Gonzalez et al., 2019; Korsakov et al., 2020; Kosminska et al., 2020) and on the environments in which lithospheric and mantle diamonds grow (e.g. Izraeli et al., 1999; Sobolev et al., 2000; Nestola et al., 2011; Anzolini et al., 2018, 2019). Current applications of elastic thermobarometry rely on several fundamental assumptions. At the time of entrapment of the inclusion it is assumed that the host and inclusion both experience the same temperature and the same homogeneous isotropic stress (i.e. hydrostatic pressure) and that, once sealed, there is no void space in the inclusion. The standard analysis (e.g. Angel et al., 2017b; Mazzucchelli et al., 2021; Angel et al., 2022a; Cisneros and Befus, 2020) is based on the assumption that the host and inclusion behave entirely elastically from entrapment until the time of measurement of the residual pressure of the inclusion, Pinc, in the laboratory. Partial or complete resetting of Pinc during exhumation by plastic flow or cracking of the host mineral is therefore not considered in the calculations, although plastic flow in the host mineral during its post-entrapment exhumation history (e.g. Ferrero and Angel, 2018; Anzolini et al., 2019; Campomenosi et al., 2021) can reset the inclusion stress state so that the measured Pinc may indicate the conditions of resetting rather than entrapment (Zhong et al., 2018; Campomenosi et al., 2023). The standard analysis also assumes that the inclusion is spherical and is elastically isolated within the host and that the host and inclusion are both elastically isotropic. With these assumptions the measured Pinc can be used with the equations of state (EoSs) of the inclusion and host phases to calculate an entrapment isomeke, a line in PT space that represents possible entrapment conditions (Adams et al., 1975; Rosenfeld and Chase, 1961). Several programs are available that perform this calculation (e.g. Kohn, 2014; Angel et al., 2017b; Mazzucchelli et al., 2021; Cisneros and Befus, 2020).

Minerals are not elastically isotropic, which means that even if they are trapped in cubic mineral hosts they will be under non-hydrostatic stress when measured in the laboratory (e.g. Bonazzi et al., 2019; Murri et al., 2022). But it has been shown for common inclusion minerals such as zircon and quartz that the biggest influence of this deviatoric stress is on the frequencies of Raman bands used to determine the inclusion “pressure” (Bonazzi et al., 2019). If the mode Grüneisen tensor (Angel et al., 2019b) of the inclusion mineral is used to determine its strain state from the measured Raman frequencies and then its stress state (Gilio et al., 2021a), the average normal stress can then be used as a good approximation of Pinc from which a reliable entrapment isomeke can be calculated (Bonazzi et al., 2019). Corrections to the measured Pinc of non-spherical inclusions can be calculated so that the model can still be used to calculate their entrapment conditions (Mazzucchelli et al., 2018). A full analysis of systems in which both the host and the inclusion are elastically anisotropic is currently under development (Mazzucchelli et al., 2019; Gonzalez et al., 2021).

All of this analysis has been focused on single-phase inclusions. But many inclusions contain more than one phase, some of the most obvious being the multiphase “nanogranitoid” inclusions derived from the crystallisation of melt inclusions (Cesare et al., 2015) and those containing two phases together as a result of partial back transformation of a high-pressure polymorph (e.g. quartz and coesite; Chopin, 1984; Ye et al., 2001). Fluid films are often found around silicate inclusions in diamonds (Nimis et al., 2016) and around diamond inclusions in silicates along with graphite (e.g. Perraki et al., 2006; Kotková et al., 2021). Diamonds are also found in fluid-dominated inclusions (Frezzotti et al., 2011). Rims of amorphous material occur around diamond inclusions in silicate hosts (e.g. Kotková et al., 2021; Jakubová et al., 2022) and oxide inclusions in silicates (Alifirova et al., 2022). In this paper we therefore describe the theory of the elasticity of mixtures of different phases to calculate the elastic properties of mixed-phase inclusions. We start with an analysis of the elastic properties of multiphase inclusions that contain sufficient fluid to ensure hydrostatic pressure conditions throughout the inclusion and use this in the isotropic analysis of host–inclusion systems, as described above, to calculate the effect of the fluid on inclusion pressures. We then extend this analysis to multiphase inclusions without fluids. This allows us to determine when the Pinc of various types of multiphase inclusions can be used to reliably calculate their entrapment conditions and when they cannot.

2 Methods

In the context of this work we use the term “fluid” to indicate any material that cannot support significant shear stress because it has small or zero shear moduli. Thus fluid can indicate a gas, liquid or melt, with or without dissolved species. If there is sufficient fluid in an inclusion, then there will be no elastic interactions between the solid mineral grains, and all of the minerals within the fluid and the fluid itself will be under the same hydrostatic pressure. In such situations (e.g. Fig. 1a, b) the fluid is said to “percolate”, which means that it occupies a single continuous connected volume. The temperature is also assumed to be uniform across the inclusion and the host mineral at all times. All other situations, especially those in which the fluid is not a continuous phase in the inclusion or the mineral grains form a continuous network (e.g. Fig. 1c), lead to stress gradients within the solid phases, and no simple thermodynamic analysis with EoS can be made. For the purposes of making the presentation of the thermodynamics straightforward, we further assume that there are no reactions between the minerals and the fluid or between the minerals themselves, there is no precipitation of solid from the fluid and no dissolution of the minerals or the host mineral into the fluid, and there is no loss of volatiles (e.g. Ni et al., 2017) from the inclusion. Thus, we assume that the chemical composition and amount of each mineral phase and that the fluid remains constant. Under these conditions the PVT relation of each mineral phase in the inclusion is defined by its own hydrostatic EoS. Note that this approach can incorporate the effects of isochemical phase transitions in the solid mineral phases (e.g. the αβ transition in quartz) as well as in the fluid. Under these constraints and assumptions the PVT relationship of a composite inclusion can be derived from the individual EoS of the solid and fluid phases in the inclusion. The expressions that we derive here also apply in general to any mixture of phases in a closed system under uniform P and T. For extension to closed systems in which there are chemical reactions between the phases, see Stixrude and Lithgow-Bertelloni (2022).

Figure 1Various configurations of fluid (blue) and solid (white and grey) in composite inclusions. In (a) and (b) the fluid occupies a single continuous volume, and the solid does not form a continuous network across the inclusion. The pressure of the fluid will be hydrostatic and the same as the solid phase or phases. This situation is exactly analogous to experiments in a diamond anvil pressure cell with a hydrostatic pressure medium. (c) Examples in which the solid phases form a continuous network so there will be stress gradients in the solid and the stress state in the inclusion will not be hydrostatic.


2.1 EoS of mixtures

The derivation of the elastic properties and EoS of a composite inclusion is easiest to follow if we use the molar volumes of the phases, Vmi, and define the number of moles of each phase i in the inclusion as mi (see Appendix). The total volume Vi of each phase in the inclusion is then miVmi, where mi remains constant in a given inclusion, and Vmi at any P and T is given by the EoS of the phase i. The total volume, V, of the mixture of the phases in the inclusion at any P and T is the sum of the volumes of the individual phases:

(1) V = i = 1 n phase m i V m i .

The volume fraction fi of any phase in the inclusion is given by

(2) f i = m i V m i V .

It is important to note that while the number of moles of each phase mi remains constant because the inclusion is a closed system that does not undergo chemical exchange between the phases, the volume fractions fi will vary with P and T because the phases have different elastic properties. We will show below how the volume fractions fi vary with pressure.

The volume thermal expansion coefficient α of the mixture of the phases in the inclusion is obtained by taking the temperature derivative of Eq. (1):

(3) α = 1 V V T = 1 V i = 1 n phase m i V m i T .

The thermal expansion coefficient of a single phase is αi=1VmiVmiT, so substitution into Eq. (3) shows that the thermal expansion of the mixture is the sum of the individual thermal expansion coefficients of the phases weighted by their volume fractions but not by the molar fractions of the phases:

(4) α = i = 1 n phase m i V m i V α i = i = 1 n phase f i α i .

The volume compressibility β=1VVP of the inclusion can be obtained in the same way by taking the derivative of Eq. (1) with respect to pressure:

(5) β = i = 1 n phase m i V m i V β i = i = 1 n phase f i β i .

Equations of state are normally expressed not in terms of the compressibility but in its inverse, the bulk modulus K=-VPV=-1β. From Eq. (5) the bulk modulus of the mixture of phases in the inclusion is thus

(6) 1 K = i = 1 n phase f i K i .

This is the well-established result that the Reuss bound on the elastic properties of a mixture of phases all under the same hydrostatic pressure is determined, in general, as a linear sum of the compliances or compressibilities of the components and not as a linear sum of the moduli (e.g. Hill, 1963). The latter is the Voigt bound, which is appropriate for situations when all of the components of the mixture are subject to the same isotropic strain but not the same pressure. Although the Voigt bound is appropriate for the composite properties of a host and inclusion in the fictive “unrelaxed state” (Angel et al., 2014b) where the inclusion is subject to the volume strain imposed by the host mineral, the Voigt bounds are not relevant for the properties of a mixture including a fluid and are therefore not discussed further in this work.

The rate of change of the bulk modulus with pressure K=K/P is also an important parameter of EoS. It can be obtained for the mixture of phases in the inclusion by taking the pressure derivative of Eq. (5); thus,

(7) β P = i = 1 n phase f i β i P = i = 1 n phase f i β i + β i f i P ,

in which we use βi=βi/P. Equation (7) shows that the change in the compressibility of a mixture with pressure is more complicated than might be expected. It is made up of two contributions. The first, fiβi, is the sum over the compressibility derivatives of the individual phases, weighted by their volume fractions, and it represents the change in the total compressibility with pressure if the volume fractions of the phases fi do not change with pressure. The second term represents the contribution from the changing volume fractions of the phases as pressure is changed and can lead to some unexpected properties. It can be evaluated by taking the pressure derivative of Eq. (2):

(8) f i P = m i V V m i P - V m i V V P = m i V m i V 1 V m i V m i P - 1 V V P = f i β i - β .

This shows that the volume fractions of the phases in mixtures remain constant only if their compressibilities (or bulk moduli) are the same as the compressibility of the mixture. If there are only two phases in the mixture, this would require the two phases to have identical compressibilities and thus EoSs. Substituting the result of Eq. (8) back into Eq. (7) and using Eq. (5) gives the pressure derivative of the Reuss bound on the compressibility:

(9) β P = - β 2 + i = 1 n phase f i β i + i = 1 n phase f i β i 2 ,

from which the pressure derivative of the Reuss bulk modulus of the mixture follows as

(10) K = K P = K 2 β P = 1 β 2 β P .

The form of these relationships between the properties of the individual phases and those of the mixture shows that, while the molar proportions of the phases in a mixture remain constant in a closed system such as an inclusion, the elastic properties of the mixture depend (Eqs. 4, 6, 10) on the volume fractions that change (Eq. 8) with P and T. This has the important consequence that it is not possible to define the parameters of an EoS for a mixture by simply taking an average of the EoS parameters of the phases at reference conditions weighted by either volume or molar fractions. Instead, the properties of a mixture can only be defined as the appropriate sums, given above, of the properties of the individual phases calculated at P and T from their individual EoS.

2.2 Mixtures in inclusions

Apart from the calculation of the elastic properties of the mixed-phase inclusion outlined in the previous section, host–inclusion calculations for multiphase inclusions can be based on the standard idealised model of an isotropic isolated spherical inclusion embedded within an isotropic homogeneous host phase (Angel et al., 2017b). The effect of the shapes of non-spherical multiphase inclusions on measured Pinc values can be corrected by a shape factor (Mazzucchelli et al., 2018). The final measured inclusion pressure Pinc can then be considered to arise from two contributions, Pthermo and ΔPrelax:

(11) P inc = P thermo + Δ P relax .

The first term, Pthermo, is the pressure generated in the inclusion due to the volume change imposed upon it by the volume change of the host on exhumation from entrapment to room conditions (Angel et al., 2017b). This is easily calculated from the EoSs of the host and the inclusion phases (e.g. Angel et al., 2014b, 2017b). However, it is a fictive pressure that is never actually measurable because an inclusion at Pthermo cannot be in mechanical equilibrium with its host at room pressure. If Pthermo is greater than room pressure, the inclusion therefore expands within the host and the inclusion pressure decreases. This has been described as the mutual elastic relaxation of the system and is denoted ΔPrelax (e.g. Angel et al., 2014b, 2017b; Zhang, 1998). In general, when Pthermo is more positive than the external final pressure, then ΔPrelax<0, and vice versa. This means that the final measured inclusion pressure Pinc always lies between Pthermo and the external pressure.

The magnitude of ΔPrelax can be calculated in several ways by relying on various assumptions. One approach uses the concept of the entrapment isomeke (Rosenfeld and Chase, 1961; Adams et al., 1975; Angel et al., 2014b), which is the line through entrapment conditions along which the volume ratios V/Vtrap of a free host crystal and a free inclusion crystal would be equal. The pressure of the isomeke can be determined at any T from the EoS of the host and inclusion by finding the pressure at which their two V/Vtrap curves cross (Fig. 2). The pressure of the entrapment isomeke at the temperature at which Pinc is measured is denoted Pfoot. The entrapment isomeke is a thermodynamic concept whose position depends only upon the entrapment P and T and the EoSs of the host and inclusion. The change in inclusion volume along the entrapment isomeke is the same as for a free crystal under the same change in P and T. Therefore the inclusion pressure remains equal to the external pressure on the entrapment isomeke. Hence, the entrapment isomeke is useful because it separates the areas of P and T space where the inclusion is at higher pressures than the external pressure applied to the host from areas where the inclusion pressure is lower than the external pressure (e.g. Rosenfeld and Chase, 1961; Angel et al., 2015b).

Figure 2Plot of V/Vtrap against pressure at 25 C for quartz and zircon for entrapment at 1.5 GPa and 900 C. The crossing point of these two curves corresponds to the pressure on the entrapment isomeke at 25 C, denoted as Pfoot. EoS from Angel et al. (2017a) and Ehlers et al. (2022).


The slope of an isomeke between the host and inclusion can be calculated from their elastic properties:

(12) P T isomeke = α I - α H β I - β H ,

in which the subscripts “H” and “I” indicate the properties of the host and inclusion phases respectively. Because the inclusion pressure remains equal to the external pressure on the entrapment isomeke, the inclusion is in mechanical equilibrium with its host, and there is no relaxation. The isomeke therefore provides a reference state of uniform stress Pfoot at room T from which ΔPrelax can be calculated isothermally (Angel et al., 2014b, 2017b). The effects of mixtures on ΔPrelax can be understood by using a linearised approximation of the approach of Angel et al. (2014b):

(13) Δ P relax K I K H K 21 P foot ,

in which K21=KI-KHKI-43GH is a function of the bulk moduli as well as the shear modulus GH of the host.

2.3 Implementation in EosFit

The equations given in Sect. 2.1 have been coded into the EoS module of the CrysFML Fortran subroutine library (Rodriguez-Carvajal and Gonzalez-Platas, 2003) that underlies the EosFit suite of programs (Angel et al., 2014a; Gonzalez-Platas et al., 2016) for EoS calculations. While the properties of a mixture can be calculated at any P and T directly through these equations from the EoS of the individual phases, it is not possible to invert Eq. (1) in order to directly calculate the P at a known V and T. To overcome this problem the calculation of P is done in two steps within the EoS module of the CrysFML. First, a simple search in P at constant T is performed until the calculated volume Vcal is approximately equal to the specified V. Then VcalV is calculated for a series of P values, and the array of values is fitted with a spline to determine the pressure at which Vcal-V=0. The use of a spline interpolation means that the precision of the resulting pressure value is not significantly limited by numerical precision.

Calculations at any P and T of the elastic properties (K, K and α), volume and pressure of mixtures of phases from their individual EoS can be performed in a dedicated utility, named mphase (for multiphase calculations), within the EosFit7c console program (Angel et al., 2014a). The program is freely available at in compiled form for Windows™, Linux™ and MacOS™ operating systems. The architecture of the EosFit7c program also allows it to be called from other software such as MATLAB™ to perform EoS calculations, including those in the mphase utility, without the need to cross-compile code directly with the CrysFML library.

The EoS parameters of individual phases can be loaded into the mphase utility from the standard “.eos” files used throughout the EosFit program suite. For convenience, the file format has been extended to allow the EoSs of several phases to be stored together in a single file. The mphase utility also includes the common EosFit7c commands to input the EoS directly or to modify EoS parameters and associated metadata. The EoS of fluid phases can be imported directly from tables of V at P and T values in text files. The elastic properties of fluid EoSs provided in this way are calculated by numerical differentiation and interpolation on the table values. In order for calculations in mphase to be meaningful, the pressure units for the bulk moduli of the various phases must all be in the same units (e.g. kbar or GPa) and the reference volume given for the EoS of each phase must be the molar volume in consistent units (e.g. cm3 mol−1). This is checked by the mphase utility before calculations can be made, and inconsistencies are clearly reported to the user.

Within mphase the molar quantities of each phase, denoted mi above, are termed the “stoichiometries” of the phases. Properties can be calculated and listed for both a fixed stoichiometry of a mixture over ranges of P and T or for ranges of stoichiometry and volume fractions at a single P and T. Isochors of mixtures can also be calculated. The use of negative stoichiometries for some phases enables the variation of ΔV of reactions to be calculated as a function of P and T. In studies of multiphase inclusions the stoichiometries of the phases within the inclusion are typically not known, but the volume fractions of the phases when the host is at room conditions can be measured or estimated. Therefore, mphase includes a command to calculate the stoichiometries from the volume fractions of the phases specified by the user.

Apart from the calculation of the elastic properties of the mixed-phase inclusion, all of the calculations for host and inclusion systems use the same methods as the EosFit-Pinc program (Angel et al., 2017b). Specifically, mphase allows the calculation of isomekes between a multiphase inclusion and its host, the entrapment isomeke given an observed final inclusion pressure Pinc and the Pinc from known entrapment conditions Ptrap and Ttrap. The value of ΔPrelax (Eq. 13) is calculated by the iterative method described in Angel et al. (2017b) which uses the elastic properties of the host at the final room conditions and the bulk modulus of the inclusion at the final Pinc. In the cases in which the calculation of Pfoot fails, for reasons we discuss below, the relaxation is approximated by ΔPrelax-3KI4GHPinc (from Eq. 21 in Zhang, 1998), which is a linear-elastic approximation to our Eq. (13) (Appendix in Angel et al., 2014b). As we will show below, the cases in which Pfoot cannot be calculated normally correspond to relatively small values of Pthermo and Pinc, and this approximation does not introduce any significant discontinuity in the trends of Pinc with inclusion composition.

3 Mixed fluid–solid inclusions

3.1 Properties of fluid plus solid mixtures

The properties of fluid and solid mixtures are most easily demonstrated with a simple mixture of one fluid and one solid mineral. We express the properties in terms of the volume fraction ffluid in the mixture, with the volume fraction of the solid therefore being (1−ffluid). Equations (3) and (5) clearly show that the thermal expansion coefficient and the compressibility of this mixture are linear in the volume fraction of fluid. However, EoSs are normally parameterised in terms of the bulk modulus, which is the inverse of the compressibility. As a consequence, the fact that fluids are typically much softer than solids means that the variation in the bulk modulus is very non-linear (Fig. 3a) and that the bulk modulus decreases very rapidly as small amounts of fluid are added to the solid; the properties of a mixture with small amounts of fluid are therefore very sensitive to the exact volume fraction of the fluid. At larger volume fractions of fluid the bulk modulus becomes close to that of the pure fluid. Thus, adding a small amount of solid to a fluid has no significant effect on the bulk modulus (Fig. 3a). The magnitude of the effects of adding fluid depends on the contrast in the bulk moduli of the fluid and solid. As Fig. 3a shows, the reduction in the bulk modulus upon mixing CO2 with olivine is much more rapid than adding water, because CO2 at room conditions is a gas and thus its bulk modulus is several orders of magnitude less than that of water. For practical applications to inclusions, the effects of less than a few volume percent of any fluid will not be those shown in Fig. 3 because the solid will remain as a connected network of crystals (e.g. Fig. 1c). Thus, the properties of the inclusion will remain between those of the minerals and the fluid until the solid phases no longer form a continuous network, at which point the properties will rapidly change to those shown in Fig. 3.

Figure 3The variation of the (a) bulk modulus and (b) its pressure derivative K of solid–fluid mixtures. Properties are shown for olivine with water at 0 GPa and 25 C (black lines) and 0.5 GPa and 25 C (blue lines). Properties for olivine plus CO2 mixtures are shown as red lines. The EoSs used to calculate these properties are the BM3 isothermal model (Angel et al., 2018) for mantle composition (Fo90-92) olivine, the water EoS of Zhang and Duan (2005), and the CO2 EoS of Pitzer and Sterner (1994).


The variation of the pressure derivative of the compressibility βP (Eq. 9) is almost linear in ffluid because only the first term (in β2) is non-linear. But the behaviour of the pressure derivative K of the bulk modulus is more complex in mixtures (Fig. 3b) and deserves some detailed explanation. It also arises from the much lower bulk modulus of the fluid compared to that of the solid. This means that there is a rapid decrease in the volume fraction ffluid as pressure is increased, and this contributes to an increase in the bulk modulus of the mixture (Eq. 6, Fig. 3a) in addition to the normal stiffening of the individual phases with increasing pressure. This behaviour of K can be simply expressed as

(14) K = K f fluid × f fluid P .

Both derivatives in Eq. (14) are negative, the value of Kffluid being the slope of the lines shown in Fig. 3a and ffluidP because the compressibility of the fluid is larger and thus more negative than the mixture of fluid plus solid (Eq. 8). Thus, the value of K increases rapidly with the first addition of the fluid to the solid because Kffluid is very large and negative for small amounts of fluid. At much larger fluid fractions, corresponding to mixtures dominated by the fluid, Kffluid is smaller, so the value of K is smaller, although it always remains higher than the K of the stiffer solid (Fig. 3b). This significant non-linearity (Fig. 3b) can also be considered to be arising from the factor of 1β2 (Eq. 10) relating K to βP. Differentiation with respect to ffluid of the expression for β2 (from Eq. 5) shows that β2 has a minimum value at

(15) f fluid max = β solid β fluid - β solid ,

at which point K has a maximum value (Fig. 3b). The value of β2 can never become zero because the bulk modulus must always be positive, so the maximum value of K always remains finite. As pressure is increased the compressibility of the fluid increases faster than that of the solid, so βfluidβsolid (Eq. 15) becomes smaller, and the maximum in K therefore shifts towards larger fluid fractions as shown in Fig. 3b for olivine plus water at 0.5 GPa. At the same time, the increase in ffluidmax means that the maximum value of K is reduced (Eq. 14) by increasing pressure. Conversely, if the contrast in the compressibilities between the solid and fluid is increased, for example on changing the fluid from water to CO2, the maximum value of K is increased and it occurs at smaller fluid fractions (Fig. 3b). This pattern of variation in elastic properties for a simple mixture of a fluid with a single solid will also occur for mixtures of several solid phases with a single fluid phase provided the fluid remains much softer than the solid minerals.

3.2 Effect of fluids on inclusions in a stiff host

In the absence of fluid, and assuming similar thermal expansion coefficients for the host and the inclusion, an inclusion that has a lower bulk modulus than its host would exhibit a positive Pinc when measured in the laboratory at room conditions. We will use olivine in diamond as an example because olivine inclusions in gem-quality diamonds are surrounded by a thin film of hydrous silicic fluid (Nimis et al., 2016). Such films seem to be almost ubiquitous around silicate and oxide inclusions in lithospheric diamonds (e.g. Nestola et al., 2018; Nimis et al., 2019), and their existence may explain the hydrostatic stress state of some of the inclusions (Angel et al., 2022a). Although there is evidence that the fluid has precipitated material after entrapment and that the fluid is not pure water, we will use the properties of water (Zhang and Duan, 2005) to illustrate the possible effects on Pinc of the presence of fluids in these inclusions.

The final pressures Pinc of inclusions can be understood in terms of how their isochors compare to their host and their isomekes with their host (e.g. Rosenfeld and Chase, 1961; Adams et al., 1975; Angel et al., 2014b). The slopes of isochors of a single phase i are given by P/TVi=αiKi. Although water has a much lower bulk modulus than olivine (Fig. 3a), it also has a much larger thermal expansion coefficient in this range of P and T, and, as a consequence, the isochors of both phases are very similar (Fig. 4).

Figure 4The isochors (solid lines) of diamond (Angel et al., 2015a), olivine (Angel et al., 2018) and water (Zhang and Duan, 2005) that all pass through 5.3 GPa and 1100 C. Dashed lines are the isomekes of water and olivine with diamond and represent the P and T required for the inclusion phase to have the same pressure as the external pressure applied to the diamond host after entrapment at 5.3 GPa and 1100 C.


The isochor of diamond through P=5.3 GPa and T=1100C also passes through P=0.0 GPa and T=25C. Thus, the volume of diamond would not change from entrapment at these conditions to the examination at room conditions. In this special case the volume of an inclusion phase trapped in the diamond at P=5.3 GPa and T=1100C would not change on exhumation (apart from the mutual elastic relaxation), and therefore the Pthermo of the inclusion will be defined by the isochor of the inclusion phase through entrapment conditions. The similarity of the isochor pressures of water and olivine at 25 C (1.365 and 1.249 GPa respectively) therefore qualitatively suggests that the effect of water on the inclusion pressures of olivine may be smaller than might be expected from their great contrast in bulk moduli. The contrast between the elastic properties of diamond and water is extreme with αHαI and βHβI, so the slope of their isomeke is, from Eq. (12), very similar to that of the isochor of water (Fig. 4):

(16) P T isomeke α I β I = P T V inc = α I K I .

This is why the isochors of fluid inclusions in minerals can be used to estimate their conditions of entrapment (e.g. Roedder and Bodnar, 1980). By contrast, the much stiffer compressibility of olivine compared to that of water means that the isomeke of olivine with diamond for entrapment at P=5.3 GPa and T=1100C lies at pressures significantly above those of its isochor (Fig. 4). The consequence is that, as water is added to an inclusion of olivine trapped at these conditions, the Pfoot decreases quite rapidly. This contributes to a decrease in ΔPrelax (Eq. 13). The bulk modulus of the composite inclusion also falls rapidly (Fig. 3a), so KIKH also decreases and the amount of relaxation becomes smaller at larger fluid fractions (Fig. 5) and almost zero for a pure fluid inclusion. The result, for these entrapment conditions, is that Pinc increases with fluid content (Fig. 5). For lower entrapment pressures, such as P=4 GPa and T=1100C, the Pfoot of olivine is lower than that of pure water inclusions, so Pfoot also increases with increasing fluid content, but the value of ΔPrelax is less than the previous entrapment example (Fig. 5). Quite surprisingly, this means the presence of pure water will increase the measured inclusion pressures of olivine inclusions in diamond when they have been trapped in the lithosphere or upper mantle.

Figure 5The variation of pressures of olivine plus water inclusions in diamond as a function of the fluid volume fraction for two different entrapment conditions. The value of ΔPrelax=Pinc-Pthermo is indicated by the arrow. Note that although lower entrapment pressures lead to Pfoot for olivine inclusions being less than for water, the trend of increasing inclusion pressure Pinc applies for all entrapment conditions.


Figures 5 and 6 show that the difference in Pthermo of pure solid inclusion compared to that of the pure fluid inclusion trapped at the same conditions is the primary control on whether fluid increases or decreases Pinc. For solid inclusions in diamond Pthermo tends to decrease as the bulk modulus of the solid increases, so inclusions of garnet (K0 170 GPa) and spinels (K0 190–200 GPa) have lower Pthermo and Pinc than olivines (K0 125 GPa) trapped at the same conditions, and thus the presence of fluid also raises their inclusion pressures. Coesite (K0=101 GPa; Angel et al., 2001) and orthoenstatite (K0=105 GPa; Angel and Jackson, 2002) are sufficiently soft that Pthermo and Pinc are significantly higher than water trapped at the same conditions, and thus the presence of fluid will, in contrast to the case of olivine, garnets and spinels, reduce the inclusion pressures (Fig. 6). The effect is particularly strong for coesite whose isomekes with diamond are approximately horizontal in PT space (Angel et al., 2022a) as a result of its thermal expansion coefficient being very small (Kulik et al., 2018) and similar to that of diamond (see Eq. 12).

Figure 6The variation of pressures of coesite plus water (dashed lines) and orthoenstatite plus water (solid lines) inclusions in diamond as a function of the fluid volume fraction for the same entrapment conditions as Fig. 4 (P=5.3 GPa and T=1100C). The value of Pthermo for both minerals is higher than for pure water, so addition of fluid decreases the measured inclusion pressure Pinc in contrast to stiffer inclusion minerals such as olivine (Fig. 5).


3.3 Effect of fluids on inclusions in a soft host

The effect of the addition of the fluid to a solid inclusion phase that is stiffer than the host (the fluid always remaining softer than the host mineral) follows the same principles illustrated in Sect. 3.2. For zircon trapped in pyrope at metamorphic conditions, the effect of fluid is relatively minor, but this system serves to show how the effect of fluid can be reversed as entrapment conditions change. For entrapment at 3 GPa and 800 C the Pinc of pure water and pure zircon inclusions differ by less than 0.1 GPa, so the effect of adding water to a zircon inclusion is minimal (Fig. 7). For lower entrapment pressures (e.g. 2 GPa; Fig. 7) and/or higher entrapment temperatures, the addition of water lowers the Pinc, while Pinc is increased by the presence of water at higher pressures and lower temperatures of entrapment.

Figure 7The variation of pressures of zircon plus water inclusions in pyrope as a function of the fluid volume fraction for two different entrapment conditions. Note that for entrapment at 3 GPa and 800 C the effect of water on the Pinc is small. Lower entrapment pressures (or higher temperatures) lead to water decreasing the Pinc.


Also note in Fig. 7 the divergence of Pfoot to extreme values for low fluid contents in this system in which the solid inclusion phase is stiffer than the host. This arises because the addition of fluid to the solid inclusion reduces the bulk modulus of the inclusion and it becomes very similar to that of the host over a range of P, T and ffluid. The resulting problems for host–inclusion calculations are difficult to illustrate with zircon in garnet because of the similarity of the bulk moduli of garnet and zircon (K0 170 and 225 GPa respectively), which means that the amount of fluid required to make the bulk modulus of the inclusion equal to that of the garnet is very small, around 1 vol %. Calculations at such low water contents are very sensitive to numerical round-off errors. We can better illustrate the concepts and difficulties with the system of zircon plus water in a quartz host, in which the much larger contrast in the bulk moduli moves the critical compositions to  6.5 vol % water at room conditions and  20 vol % at 1.5 GPa and 900 C (Fig. 8a).

Figure 8(a) The variation of the bulk modulus of zircon plus water mixtures as a function of fluid content (solid lines) for room conditions (black) and 1.5 GPa and 900 C (red). The values of the bulk modulus of quartz at the same conditions of P and T are indicated by the horizontal dashed lines. (b) Calculated inclusion pressures of zircon plus water inclusions in quartz trapped at 1.5 GPa and 900 C. The grey areas indicate where Pfoot is not defined, corresponding to inclusion compositions whose bulk moduli are close to that of the quartz host at room conditions (left) and entrapment conditions (right). Note that the x axis only extends to 50 % water by volume. Pressures for more water-rich inclusions are very similar to those of an inclusion with 50 vol % water; see also the isomekes in Fig. 9.


Fluid-dominated inclusions in the system of zircon plus water in quartz with ffluid> 0.3 (30 vol % fluid) have thermoelastic properties of the inclusion dominated by the properties of the fluid (Fig. 8a), so the entrapment isomeke closely follows that for pure water inclusions (Fig. 9). The resulting Pfoot and Pinc are therefore very similar to those of pure water (Fig. 8b). As the ffluid is reduced towards 0.2, the bulk moduli of the inclusion and the host become very similar near to the entrapment conditions, so the entrapment isomeke becomes very steep in PT (Eq. 12) and does not extend to room T at reasonable or finite pressures. Therefore Pfoot cannot be calculated, and ΔPrelax has to be approximated by the alternative expression of Zhang (1998). This does not introduce significant error because for these amounts of fluid Pthermo remains small and therefore ΔPrelax is small. At slightly lower ffluid (in this case 0.08 to 0.18) the bulk moduli are sufficiently different at entrapment (Fig. 8a) for the entrapment isomeke (Fig. 9) to be well-defined. However, with decreasing P and T along the entrapment isomeke, the difference between the bulk moduli of the host and inclusion decreases and the isomeke steepens. In this example, this steepening is additionally modified by the softening of quartz associated with its αβ transition (Angel et al., 2017a), and the isomeke terminates. For other hosts without phase transitions, the isomeke would become vertical (Eq. 12). In both cases the entrapment isomekes do not extend continuously to room T. Nonetheless, the entrapment isomeke can still exist at lower T and still represents conditions under which the fractional volume change of the host and the inclusion from entrapment are equal (Fig. 9). At these lower temperatures it is not unusual to find the isomeke so strongly curved that there are two isomeke pressures for each T as shown for ffluid=0.16 and 0.10 in Fig. 9. This arises from the composite inclusion being softer than the host while its K is larger than the host due to the rapid decrease in fluid volume fraction with increasing pressure (Fig. 3b, Eqs. 7–10, 14). Consequently, the PV line of the inclusion is strongly curved compared to that of the normal EoS of the host, so the curves of V/Vtrap for the host and inclusion cross twice for a single pressure (see 10 vol % water curve in Fig. 10). Both crossing points represent the entrapment isomeke and thus possible values of Pfoot. But the relevant isomeke pressure can be identified by interpolation or extrapolation from adjacent inclusion compositions with more normal isomekes and single values of Pfoot from the same entrapment conditions. For this example, it is obvious that the correct Pfoot to be used in the range ffluid= 0.08–0.18 is the lower value that is very similar to that of the more water-rich inclusions. This can also be verified by comparing the final Pinc calculated from the two Pfoot values to the Zhang (1998) approximation. At even lower fluid contents (ffluid=0.02 to 0.06) the bulk modulus of the inclusion differs significantly from the host at entrapment conditions but instead becomes very similar to that of the host at room conditions (Fig. 8a), and, as a result, the curves of V/Vtrap for the host and inclusion against P can become parallel, or not cross at all (see 4 vol % water curve in Fig. 10), so that Pfoot is not defined. For these compositions, relaxation must be calculated with the approximation of Zhang (1998). Finally, at very low ffluid (<0.02 for this example) the properties of the composite inclusion and thus the entrapment isomeke, Pfoot and Pinc approach those of a pure zircon inclusion (Figs. 2, 9).

Figure 9Entrapment isomekes for zircon plus water inclusions in quartz, entrapped at 1.5 GPa and 900 C. Inclusion compositions in volume percent fluid at entrapment conditions are indicated in the same colours as the isomekes. The red dots show that the entrapment isomeke for 2 % water is very similar to that of pure zircon, and the blue dots for 50 % water show that the entrapment isomekes for water-rich inclusions are indistinguishable from those of pure water. Note that the entrapment isomekes for 10 % and 16 % water are terminated below entrapment by the αβ phase transition in quartz (dark grey line) but exist as separate line segments at lower temperatures. The phase boundaries for water and quartz are shown in light grey. The transformation from α to β quartz is responsible for the deflection of the isomeke for pure water at ca. 1000 C.


Figure 10Plot of V/Vtrap against pressure at 25 C for quartz and zircon plus water for entrapment at 1.5 GPa and 900 C to illustrate potential problems in calculating Pfoot for systems in which the solid inclusion phase is stiffer than the host. Normally, the host and inclusion EoSs define a single crossing point and a unique value for Pfoot as for pure zircon (black line). However, the curve for 10 vol % water plus zircon (at entrapment) crosses the curve for quartz twice, defining two values for Pfoot. This is caused by the severe curvature of the V/Vtrap line of the mixture due to its large value of K (Eqs. 7–10, 14) compared to that for pure quartz. At lower water contents (e.g. 4 vol %) the large value of K means that its V/Vtrap line does not intersect the line for quartz, so Pfoot is not defined.


So far, we have only considered the effects of a fluid for systems and entrapment conditions that would lead to the pure solid having a positive Pinc. When Pinc for a solid inclusion falls below P=0, this corresponds to a stretched state of the inclusion being kept in tension by bonding to the surrounding host mineral that has expanded more from entrapment than a free crystal of the inclusion phase would. Measured examples include quartz trapped in metamorphic garnet at high T and low P (Gilio et al., 2021b), apatite in olivine (Cisneros and Befus, 2020), and synthetic zircon in quartz (Gonzalez et al., 2022). If the Pinc of a pure fluid inclusion trapped at the same conditions would be positive, the addition of small amounts of fluid to the solid-phase inclusion would result in the calculated Pinc first becoming less negative and then positive. However, while EoSs of solids at negative pressures in the extensional regime are valid, at least for small negative pressures (Angel et al., 2019a), they are not valid for fluids or gases. Instead, Pinc will be buffered to the fluid–vapour equilibrium curve by the formation of a vapour bubble in the fluid phase, which for water lies below 0.02 GPa at temperatures below the critical point 374 C (near-horizontal grey line in Fig. 9). In these circumstances the entrapment isomeke should not be calculated from the measured Pinc. Instead, the P and T of homogenisation of the fluid component in the inclusion will define a point on the entrapment isomeke (Fig. 9).

3.4 Implications of fluids in inclusions

The examples discussed above show that the practical consequences of fluid contents of inclusions for the calculation of their entrapment conditions strongly depend on the elastic contrast between the solid phases within the inclusion and the fluid phase and between both of them and the host phase, as well as the entrapment conditions themselves. Together these control the dependence of Pinc on fluid content (e.g. Figs. 5, 6, 7 and 8). As an example, the influence of water on the pressures of olivine inclusions in diamonds is not as serious as might be expected from the contrast in elastic properties shown in Fig. 3. A 1 µm thick rim of water on an olivine-containing inclusion of 100 µm diameter represents  6 % of the total volume of the inclusion. The measured Pinc of such a water-containing inclusion trapped at 4 GPa and 1100 C will be almost 0.3 GPa higher than a pure dry olivine trapped at the same conditions (Fig. 5). If the Ptrap is then calculated from the measured Pinc by ignoring the presence of fluid and using only the EoS for olivine, the Ptrap will be 0.4 GPa too high at 1100 C. This error decreases with increasing depth of entrapment and is approximately zero for entrapment above 6.5–7 GPa. Therefore, failure to allow for the fluid content of olivine inclusions will result in estimates of depth of their entrapment that are in error by less than 12 km. This is probably significantly less than the errors introduced by not modelling the effects on Pinc of either the precipitation of solids from the fluid in these inclusions after their entrapment (Nimis et al., 2016) or the escape of fluids into cracks around the inclusions (e.g. Angel et al., 2022a). Similar levels of error, but in the opposite direction (Fig. 6), will arise from ignoring fluids in pyroxene inclusions. In contrast, the much stronger variation in Pinc with fluid content in coesite inclusions (Fig. 6) means that ignoring 6 vol % water leads to a larger underestimate of entrapment pressures, of the order of 1 GPa. This serves to illustrate that the potential errors due to low fluid contents must be carefully evaluated for each host–inclusion system in turn over the range of expected entrapment conditions. In addition, the contributions from uncertainties in the EoS of the fluid should also be evaluated by calculation with alternative EoSs when they are available.

4 Multiphase solid inclusions

4.1 Properties of solid mixtures

In the absence of a percolating fluid completely surrounding the solid crystals, the stress (and pressure) distribution in a mixture of solid phases, or indeed in a polycrystalline aggregate of one anisotropic phase, is inhomogeneous (e.g. Chap. 6 in Pollard and Fletcher, 2006) because the stress state at any point depends on how the crystal grains interact with one another (e.g. Fig. 1c). In this case the local distribution of the stress must be calculated from full mechanical calculations (e.g. Schmalholz et al., 2020). The average elastic response of an aggregate of a very large number of randomly oriented crystals falls between the Reuss and Voigt bounds. The latter corresponds to a state of uniform strain applied to each of the constituent crystals and leads to unequal pressures or stresses in each of them. The Reuss bound corresponds to what we have already described for fluid-dominated mixtures, in which all of the grains of solid are subject to the same uniform hydrostatic pressure, and it also implies that the boundaries between the grains are completely free to slip so that there are no coherency-induced strains when P or T are changed. In an inclusion the number of crystals is very limited, and it may be that neither bound applies. In principle, the stress state of inclusions containing multiple phases can be determined by measuring the stress or strain states of the individual component crystals by diffraction or spectroscopy. If the crystals are all at the same pressure, the elastic properties will be determined by the Reuss bound and the equations presented in Sect. 2.2, and they will be independent of the internal microstructure of the inclusion.

Under these assumptions, the calculations for multiphase solid inclusions are the same as for those involving a fluid, with the difference that the contrast in elastic properties will, in general, be much less because no mineral has elastic properties close to those of fluids such as water. Thus, the variation in the bulk modulus with volume fraction of solid phases is less curved than for fluids (e.g. compare Fig. 11a with Figs. 3 and 8), which means that the value of K shows a much weaker maximum than found for mixtures with fluids in agreement with Eq. (14). This also means that the volume fractions of phases in a purely solid mixture change much less with P and T than those of a solid plus fluid. Similarly, the variation in thermal expansion coefficient in a mixture of solids will be much less than for a mixture of solid plus liquid (Eq. 4). Consequently, for inclusions containing phases that are either all more stiff or all less stiff than the host there are no difficulties in calculating inclusion pressures which will lie between those expected for the individual phases. Problems arise, as for fluids, when the bulk modulus of the host lies between those of the individual phases within the inclusion.

4.2 Hosts of intermediate stiffness

This case can be illustrated (Fig. 11a) by composite rutile plus amphibole (cummingtonite) inclusions that have been found in eclogitic garnets (e.g. Musiyachenko et al., 2021). Rutile is stiffer than garnet, and therefore pure rutile inclusions trapped in the eclogite field are calculated to have Pinc<0, while amphibole is considerably softer than garnet (Fig. 11a), so pure amphibole inclusions should therefore show positive Pinc. With the addition of amphibole to rutile inclusions the calculated Pinc will therefore initially become less negative, then zero, and then positive as shown in Fig. 11b. This figure also demonstrates a significant difference in the pressures of mixed-phase solid inclusions from those including fluids. The trend of Pinc for a mixture of solids is far closer to linear (Fig. 11b) than the trends seen with fluids (e.g. Figs. 5, 6, 7, 8), so the addition of a small amount of the stiffer solid component to the softer solid has an effect on Pinc that is similar in magnitude (but opposite in sign) to that of adding the soft component to the stiffer mineral (Fig. 11b).

Figure 11(a) The variation of the bulk modulus and K of rutile–amphibole mixtures at room conditions as a function of the volume fraction of amphibole. (b) Calculated pressures of rutile plus amphibole inclusions trapped in pyrope garnet at 3 GPa and 800 C. The grey area indicates the compositions for which Pfoot is either not defined or exceeds ±10 GPa. Calculations with the EoSs of pyrope (Angel et al., 2022b), rutile (Angel et al., 2020) and cummingtonite (Holland and Powell, 2011).


Figure 12Entrapment isomekes for rutile plus amphibole inclusions trapped in pyrope garnet at 3 GPa and 800 C. Isomekes are labelled with the volume percent of amphibole in each inclusion, and the grey arrow indicates how the entrapment isomekes rotate with increasing amphibole content. The intersection of the entrapment isomekes with 25 C (vertical dashed line) correspond to the Pfoot values shown in Fig. 11; note that the entrapment isomekes for inclusion with 9 vol %–12 vol % amphibole do not arrive to 25 C and therefore do not have a Pfoot. Calculated with the EoSs of pyrope (Angel et al., 2022b), rutile (Angel et al., 2020) and cummingtonite (Holland and Powell, 2011).


The absence of Pfoot values at intermediate compositions (Fig. 11b) can, as for the example of zircon plus water in quartz, be explained with calculations of the entrapment isomekes (Fig. 12) and the elastic properties of the phases. At entrapment conditions of 3 GPa and 800 C both rutile and amphibole have lower thermal expansion coefficients than pyrope, so any mixture of the two in an inclusion will similarly have a lower α (Eq. 4). Rutile is stiffer than pyrope over the entire PT range of interest, which means that βpyrope>βrutile, and thus the entrapment isomeke of pure rutile in pyrope has a positive slope at these entrapment conditions. However, at low pressures at ambient T, αpyrope is significantly less than that of rutile and amphibole, so the slope of the isomeke must become negative, leading to a minimum in the isomeke pressure at some intermediate T. As amphibole is added to a rutile inclusion, the bulk modulus of the inclusion at entrapment conditions is reduced towards that of the pyrope host, and the entrapment isomekes become steeper. The thermal expansion coefficient of rutile plus amphibole mixtures at any given P and T does not vary significantly with volume fraction, so the addition of amphibole to rutile leads to both a depression of the minimum of the entrapment isomeke and a decrease in Pfoot, as shown for 7 vol % and 8 vol % amphibole in Fig. 12. Further small increases in amphibole content (e.g. 9 vol % and 10 vol %) continue to reduce the bulk modulus of the inclusion towards that of the host pyrope, leading to even steeper isomekes that become undefined at lower pressures because the V/Vtrap curves of the host and inclusion no longer have a crossing point, a situation similar to that shown in Fig. 10 for the inclusion containing 4 vol % water. For such compositions, the entrapment isomeke also exists at temperatures far below room T but not at room T, and thus a Pfoot at room T does not exist. At ca. 11 vol % amphibole, the isomeke becomes vertical at entrapment conditions. At pressures just below entrapment the bulk modulus of the host is greater than that of the inclusion, so βpyrope<βinc, leading to an isomeke with a negative slope. But the larger K of the inclusion mixture (Fig. 11a) causes the bulk modulus of the inclusion to increase more rapidly than that of the pyrope host with increasing pressure, with the result that above the entrapment pressure the situation is reversed with Kpyrope<Kinc and βpyrope>βinc, so the entrapment isomeke has a positive slope. Thus, the entrapment isomeke is curved towards higher temperatures (Fig. 12), with the PV curves at any temperature resembling those of the inclusion of 10 vol % water shown in Fig. 10. With a further increase in amphibole content the bulk modulus of the inclusion mixture falls below that of the host, so the isomeke slope at entrapment has a negative slope. By 20 vol % the slope is sufficiently shallow that the higher K of the inclusion no longer has the effect of curving the isomeke back on itself, so Pfoot at room T exists but is larger than 10 GPa. Further increase in amphibole content decreases the contrast between the bulk moduli of the inclusion and host (Fig. 11a), making the isomekes shallower until that for 100 % amphibole is almost flat because of the similarity of the thermal expansion coefficients of amphibole and pyrope at 3 GPa for all temperatures. Note that this anti-clockwise rotation of the isomekes as amphibole is added to rutile means that the isomekes for mixtures of the two phases never lie between those of the pure phases (Fig. 12).

This example also shows that the exact behaviour of entrapment isomekes is very sensitive to the details of the EoSs of both the host and inclusion phases and the relative change in their elastic properties with P and T, especially at compositions for which the thermal expansion and/or compressibility of the inclusion are very similar to those of the host mineral. As one example, if the EoS of Milani et al. (2015) for pyrope is used instead of that of Angel et al. (2022b), then the isomekes for 10 vol %–12 vol % amphibole have the opposite curvature to those shown in Fig. 12 and instead curve back to yield two values of Pfoot at room temperature. Pinc for these entrapment conditions is calculated to be zero not at  11 vol % but at  14 vol % amphibole (Musiyachenko et al., 2021). Nonetheless, the comparison of Fig. 11 with Fig. 9 of Musiyachenko et al. (2021) shows that while relatively small changes in mineral elastic properties can significantly change the isomekes for compositions whose properties are close to those of the host, the values of Pthermo and Pinc vary much less and are therefore more robust. Conversely, the calculation of entrapment isomekes from measured values of Pinc becomes less reliable as the properties of the inclusion approach that of the host, whether at the final conditions or at conditions around entrapment.

The example of rutile inclusions in garnet also serves to illustrate one further point. Measurements of pure rutile inclusions generally show them to be at room pressure (Musiyachenko et al., 2021; Cisneros and Befus, 2020) instead of having negative Pinc expected from calculations. This does not indicate that the calculations described here are incorrect but shows that instead of being “stretched” into tension by bonding to the surrounding host garnet, the rutile inclusion crystals have been released by the host. The inclusions must therefore include either a void space due to contraction of the rutile back to the volume of a free rutile crystal at room conditions or an interface region of a lower-density material (Alifirova et al., 2022). In such cases we then expect Pinc=0 for amphibole–rutile mixtures for which the calculated Pinc<0 and then increasing positive values of Pinc as the amphibole content is increased further. The measured Pinc=0 for inclusions that have undergone contraction cannot be used to infer entrapment conditions using the methods described in this paper.

5 Summary

Mixtures of mineral phases are, of course, common as rocks. The thermodynamic properties of a mixed-phase inclusion, which are the same as those of the same mixture in a rock closed to chemical exchange, are not simple sums of the properties of the individual phases. This means that it is not possible to define the parameters of a single EoS for a mixture by simply taking an average of the EoS parameters of the phases at reference conditions weighted by either volume or molar fractions. Thermodynamic calculations of phase equilibria in rocks normally also consider that the P and T are externally imposed and that the volume of the system is free to respond to changes in these two intensive variables and is not modified by elastic interaction between the phases. Multiphase inclusions obviously differ from mixtures of phases in rocks in two important respects; the volume of the system is imposed by the host mineral (in response to the external P and T) which also interacts elastically with the inclusion. When the host–inclusion system is subject to uniform temperature, the inclusion pressure then becomes the dependent variable and can be calculated as the sum of two contributions (Eq. 11). The first, Pthermo, arises from the imposition of the change in the volume of the host mineral on the closed inclusion and is calculated without ambiguity from the EoS of the host and that of the phases comprising the composite inclusion. The second, ΔPrelax, arises from the elastic interaction between the inclusion and the host (e.g. Angel et al., 2014b, 2017b; Zhang, 1998). Calculations for multiphase inclusions have been implemented in the EoS module of the CrysFML Fortran subroutine library (Rodriguez-Carvajal and Gonzalez-Platas, 2003) and made available to users in the EosFit7c program (Angel et al., 2014a). If it is assumed that the properties of a solid solution are ideal and the coefficients of volume compressibility and thermal expansion vary linearly with molar proportions of the end-members (but see Myhill, 2018), then the equations given here and the mphase utility of EosFit7c can also be used to calculate the EoS of any composition in the solid solution.

We have shown that, as for single-phase inclusions, it is the contrast in the bulk moduli and the thermal expansivities between the host and the inclusion that determines Pthermo, which is quite often more important than the relaxation in determining the final inclusion pressure. Problems and uncertainties arise when the properties of the inclusion become similar to those of the host. This will always occur for certain volume fractions in two-phase inclusions when one phase is stiffer than the host and the other softer, whether the softer phase is a fluid or a solid. At these inclusion compositions the entrapment isomekes can become vertical in PT space and can become undefined, as illustrated by a zircon plus 4 vol % water inclusion in quartz (Fig. 10). Often associated with such behaviour is the occurrence of entrapment isomekes that are so strongly curved that two Pfoot are defined for a given temperature (Fig. 10), creating difficulties for the calculation of ΔPrelax. A third possibility that we have not illustrated is when the curves of V/Vtrap for both the host and inclusion become identical over a pressure range because both K and K of host and inclusion are the same. For this last case there is no unique Pfoot. These situations only affect the ability to calculate the entrapment isomeke and Pfoot, while Pthermo remains well-defined in all cases and the relaxation can be approximated by other methods not involving Pfoot. This can be understood in another way by considering a small portion of a single mineral grain as an inclusion with exactly the same composition and EoS as its surrounding host. Then the “inclusion”, which is really part of the host crystal, will always be at the same P and T as the host. Therefore, the host and inclusion have the same fractional volume change from entrapment to any other P and T, the isomeke is no longer a unique line in PT space, and Pfoot cannot be defined. Nonetheless, in this special case Pthermo=Pinc are always the same as the external pressure, and thus Pthermo can be calculated from the EoS of the mineral. In the cases where Pfoot cannot be determined, an alternative method of calculating ΔPrelax must be employed. We have found that a linear-elastic approximation to Eq. (13) provides a good estimate of ΔPrelax as -3KI4GHPinc, but there are other alternative methods (e.g. Gonzalez et al., 2021; Zhong et al., 2021; Morganti et al., 2020) that should converge to very similar values of Pinc.

The calculated inclusion pressures depend on the volume fractions of the phases in the inclusion, which can be measured or estimated by a variety of methods including 3D Raman mapping (e.g. Musiyachenko et al., 2021), diffraction or optical imaging of the entrapped inclusion, or mapping of an inclusion exposed for further analysis, for example by microprobe. The sensitivity, and hence reliability, of calculated entrapment pressures to the volume fractions of the constituent phases depends on the contrast in the elastic properties of the inclusion phases themselves. For inclusions containing two or more solid phases the calculated Pinc values for the same entrapment conditions are approximately linear with the volume fractions of the phases (e.g. Fig. 11), and the same applies to Ptrap values calculated from measured Pinc. The uncertainty in Ptrap for mixed solid inclusions therefore depends on the contrast in elastic properties, and especially bulk moduli, of the phases involved. If this contrast is small, then Ptrap values will be insensitive to estimates of the volume fractions of the phases in the inclusion. Conversely, although it might be thought that the large contrast between the bulk moduli of solids and fluids would always lead to very strong variations in Pinc values with the amount of fluid in the inclusion, this is not always the case because there can be a compensating effect of the contrast in thermal expansion coefficients (e.g. Fig. 4). For example, the variation in Pinc values with fluid fraction in silicate inclusions in diamond differs considerably between olivine, coesite and orthopyroxene (Figs. 5, 6). Therefore, no general guidelines can be given for the sensitivity of Ptrap values to volume fractions of fluids, but we have provided the tools for calculations to be made for any inclusion system.

The examples that we have discussed illustrate that the measured Pinc of mixed-phase inclusions must be interpreted with care. In addition to considering the uncertainties that may be introduced into host–inclusion calculations by uncertainties in phase fractions in the inclusion, and when the properties of the two become similar, one also should consider whether the inclusion has undergone reactions or transformations of any sort between the phases. Such changes would affect both the volumes of the phases due to changes in composition and their elastic properties and cannot be easily generalised. They are therefore excluded from our current analysis because they should be analysed on a case-by-case basis. One can, however, make the general statement that any reactions or transformations in the inclusion that would lead to an increase in the volume of an unconstrained system will result in the inclusion pressure being buffered to a higher value than expected had the reaction or transformation not occurred (Gillet et al., 1984; Perrillat et al., 2003; Anzolini et al., 2016; Arlt and Angel, 2000).

Lastly, careful consideration must be given to the microstructure of the inclusion before applying the analysis described in this paper. Our analysis strictly only applies to inclusions in which there are no significant elastic interactions between the solid phases, a situation which occurs when the crystals are bathed in sufficient fluid (Fig. 1); one can say that the fluid percolates across the inclusion but the solids do not. As the amount of fluid is decreased there will be a point at which the solids also percolate and form a continuous network (e.g. Fig. 1c). This will lead to stress gradients within the solid grains, and the methods presented here should be used with caution unless it can be demonstrated that the mineral grains are, on average, all at the same pressure and that the stress gradients are limited either in magnitude or extent within the grains. Otherwise, the stress state of multiphase inclusions must be analysed explicitly to account for the elastic interaction between the grains, especially for microstructures in which the elastic interactions are expected to be significant and the stress state is unlikely to be isotropic and hydrostatic. Examples would include inclusions within inclusions or grains containing coherently exsolved phases.

Appendix A: Definition of symbols used in this paper
Equation of state parameters, single phases
Vmi Molar volume of phase i
αi Thermal expansion coefficient of phase i
βi Volume compressibility of phase i
Ki Bulk modulus (inverse of volume compressibility) of phase i
Ki Pressure derivative of the bulk modulus of phase i;
a subscript 0, as in V0 or K0, indicates the property at reference conditions
Parameters for mixtures of phases
mi Number of moles of phase i in a mixture (constant in a given mixture)
fi Volume fraction of phase i in a mixture (varies with P and T)
V,K, K, α, β Equation of state parameters of a mixture
Inclusion pressures
Ptrap The pressure of entrapment of the inclusion
Ttrap The temperature of entrapment of the inclusion
Pinc The final measured pressure of the inclusion, normally when the host mineral is at room conditions
Pfoot The pressure of the entrapment isomeke at the temperature of measurement of the inclusion pressure Pinc
Pthermo The pressure the inclusion would have if there was no elastic relaxation; it is the pressure of the inclusion subject to the volume change in the host
ΔPrelax The pressure change in the inclusion due to the mutual elastic relaxation: ΔPrelax=Pinc-Pthermo
KH,αH Bulk modulus and thermal expansion coefficient of the host mineral
GH Shear modulus of the host mineral
KI, αI Bulk modulus and thermal expansion coefficient of the inclusion
Code availability

The routines used for calculation are available as the cfml_eos module (Angel et al., 2014a) which is part of the CrysFML Fortran subroutine library (Rodriguez-Carvajal and Gonzalez-Platas, 2003) from

Data availability

There are no original data reported in this paper. Citations are given in the text for all of the equations of state used in calculations.

Author contributions

Conceptualisation and initial development of the thermodynamic model by RJA, MA and FN. Writing of code in EosFit by RJA. Development and testing of algorithms for host–inclusion calculations by RJA, MLM and KAM.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Probing the Earth: Melt and solid inclusions as probes to understand nature”. It is not associated with a conference.


This project was originally developed with funding from the European Research Council to Fabrizio Nestola and completed under the European Union's Horizon 2020 research and innovation program grant agreement 714936 TRUE DEPTHS to Matteo Alvaro. We thank Joe Gonzalez and Mattia Gilio (Pavia) and Hugo W. van Schrojenstein Lantman (Oslo) for discussions on various aspects of this work and Yuuki Hagiwara (JAMSTEC) for developing and testing the table for the EoS of CO2. We thank Miguel Cisneros, Frank Spear and one anonymous reviewer for their comments that helped improve the manuscript.

Financial support

This research has been supported by the H2020 European Research Council (grant no. 714936).

Review statement

This paper was edited by Silvio Ferrero and reviewed by Miguel Cisneros, Frank Spear and one anonymous referee.


Adams, H. G., Cohen, L. H., and Rosenfeld, J. L.: Solid inclusion piezothermometry I: comparison dilatometry, Am. Mineral., 60, 574–583, 1975. 

Alifirova, T., Daneu, N., Kohn, K., Griffiths, T., Rečnik, A., and Habler, G.: Atomic scale structure of garnet-rutile interfaces in host-inclusion systems, 23rd General Meeting of the International Mineralogical Association, Lyon, France, 2022. 

Angel, R. J., Mosenfelder, J. L., and Shaw, C. S. J.: Anomalous compression and equation of state of coesite, Phys. Earth Planet. In., 124, 71–79, 2001. 

Angel, R. J. and Jackson, J. M.: Elasticity and equation of state of orthoenstatite, MgSiO3, Am. Mineral., 87, 558–561, 2002. 

Angel, R. J., Gonzalez-Platas, J., and Alvaro, M.: EosFit7c and a Fortran module (library) for equation of state calculations, Z. Kristallogr., 229, 405–419, 2014a. 

Angel, R. J., Mazzucchelli, M. L., Alvaro, M., Nimis, P., and Nestola, F.: Geobarometry from host-inclusion systems: the role of elastic relaxation, Am. Mineral., 99, 2146–2149, 2014b. 

Angel, R. J., Alvaro, M., Nestola, F., and Mazzucchelli, M. L.: Diamond thermoelastic properties and implications for determining the pressure of formation of diamond-inclusion systems, Russ. Geol. Geophys., 56, 211–220,, 2015a. 

Angel, R. J., Nimis, P., Mazzucchelli, M. L., Alvaro, M., and Nestola, F.: How large are departures from lithostatic pressure? Constraints from host-inclusion elasticity, J. Metamorph. Geol., 33, 801–803, 2015b. 

Angel, R. J., Alvaro, M., Miletich, R., and Nestola, F.: A simple and generalised P-T-V EoS for continuous phase transitions, implemented in EosFit and applied to quartz, Contrib. Mineral. Petr., 172, 29,, 2017a. 

Angel, R. J., Mazzucchelli, M. L., Alvaro, M., and Nestola, F.: EosFit-Pinc: A simple GUI for host-inclusion elastic thermobarometry, Am. Mineral., 102, 1957–1960,, 2017b. 

Angel, R. J., Alvaro, M., and Nestola, F.: 40 years of mineral elasticity: a critical review and a new parameterisation of equations of state for mantle olivines and diamond inclusions, Phys. Chem. Miner., 45, 95–113, 2018. 

Angel, R. J., Miozzi, F., and Alvaro, M.: Limits to the validity of thermal-pressure equations of state, Minerals, 9, 562,, 2019a. 

Angel, R. J., Murri, M., Mihailova, B., and Alvaro, M.: Stress, strain and Raman shifts, Z. Kristallogr., 234, 129–140,, 2019b. 

Angel, R. J., Alvaro, M., Schmid-Beurmann, P., and Kroll, H.: Commentary on “Constraints on the Equations of State of stiff anisotropic minerals: rutile, and the implications for rutile elastic barometry”, Mineral. Mag., 84, 339–347,, 2020. 

Angel, R. J., Alvaro, M., and Nestola, F.: Crystallographic methods for non-destructive characterization of mineral inclusions in diamonds, in: Natural Diamonds: Their Mineralogy, Geochemistry and Genesis, edited by: Pearson, G. and Smit, K., Reviews in Mineralogy and Geochemistry, Mineralogical Society of America, Washington DC, 257–306,, 2022a. 

Angel, R. J., Gilio, M., Mazzucchelli, M. L., and Alvaro, M.: Garnet EoS: a critical review and synthesis, Contrib. Mineral. Petr., 177, 54,, 2022b. 

Anzolini, C., Angel, R. J., Merlini, M., Derzsi, M., Tokár, K., Milani, S., Krebs, M. Y., Brenker, F. E., Nestola, F., and Harris, J. W.: Depth of formation of CaSiO3-walstromite included in super-deep diamonds, Lithos, 265, 138–147, 10.1016/j.lithos.2016.09.025, 2016. 

Anzolini, C., Prencipe, M., Alvaro, M., Romano, C., Vona, A., Lorenzon, S., Smith, E. M., Brenker, F. E., and Nestola, F.: Depth of formation of super-deep diamonds: Raman barometry of CaSiO3-walstromite inclusions, Am. Mineral., 103, 69–74,, 2018. 

Anzolini, C., Nestola, F., Mazzucchelli, M. L., Alvaro, M., Nimis, P., Gianese, A., Morganti, S., Marone, F., Campione, M., Hutchison, M. T., and Harris, J. W.: Depth of diamond formation obtained from single periclase inclusions, Geology, 47, 219–222,, 2019. 

Arlt, T. and Angel, R. J.: Pressure buffering in a diamond anvil cell, Mineral. Mag., 64, 241–245, 2000. 

Bonazzi, M., Tumiati, S., Thomas, J., Angel, R. J., and Alvaro, M.: Assessment of the reliability of elastic geobarometry with quartz inclusions, Lithos, 350–351, 105201,, 2019. 

Campomenosi, N., Scambelluri, M., Angel, R. J., Hermann, J., Mazzucchelli, M. L., Mihailova, B., Piccoli, F., and Alvaro, M.: Using the elastic properties of zircon-garnet host-inclusion pairs for thermobarometry of the ultrahigh-pressure Dora-Maira whiteschists: problems and perspectives, Contrib. Mineral. Petr., 176, 36,, 2021. 

Campomenosi, N., Angel, R. J., Alvaro, M., and Mihailova, B.: Resetting of zircon inclusions in garnet: implications for elastic thermobarometry, Geology, 51, 23–27,, 2023. 

Cesare, B., Acosta-Vigil, A., Bartoli, O., and Ferrero, S.: What can we learn from melt inclusions in migmatites and granulites?, Lithos, 239, 186–216,, 2015. 

Chopin, C.: Coesite and pure pyrope in high-grade blueschists of the Western Alps: a first record and some consequences, Contrib. Mineral. Petr., 86, 107–118, 1984. 

Cisneros, M. and Befus, K. S.: Applications and limitations of elastic thermobarometry: insights from elastic modeling of inclusion-host pairs and example case studies, Geochem. Geophy. Geosy., 21, e2020GC009231,, 2020. 

Ehlers, A. M., Zaffiro, G., Angel, R. J., Boffa-Ballaran, T., Carpenter, M. A., Alvaro, M., and Ross, N. L.: Thermoelastic properties of zircon: implications for geothermobarometry, Am. Mineral., 107, 74–81,, 2022. 

Ferrero, S. and Angel, R. J.: Micropetrology: are inclusions grains of truth?, J. Petrol., 59, 1671–1700,, 2018. 

Frezzotti, M. L., Selverstone, J., Sharp, Z. D., and Compagnoni, R.: Carbonate dissolution during subduction revealed by diamond-bearing rocks from the Alps, Nat. Geosci., 4, 703–706, 2011. 

Gilio, M., Angel, R. J., and Alvaro, M.: Elastic geobarometry: how to work with residual inclusion strains and pressures, Am. Mineral., 106, 1530–1533,, 2021a. 

Gilio, M., Campomenosi, N., Musiyachenko, K. A., Angel, R. J., Cesare, B., and Alvaro, M.: Elastic geobarometry of quartz inclusions in garnet at high temperature, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-1508,, 2021b. 

Gillet, P., Ingrin, J., and Chopin, C.: Coesite in subducted continental crust: P-T history deduced from an elastic model, Earth Planet. Sc. Lett., 70, 426–436, 1984. 

Gonzalez, J. P., Thomas, J. B., Baldwin, S. L., and Alvaro, M.: Quartz-in-garnet and Ti-in-quartz thermobarometry: Methodology and first application to a quartzofeldspathic gneiss from eastern Papua New Guinea, J. Metamorph. Geol., 37, 1193–1208,, 2019. 

Gonzalez, J. P., Mazzucchelli, M. L., Angel, R. J., and Alvaro, M.: Elastic geobarometry for anisotropic inclusions in anisotropic host minerals: quartz-in-zircon, J. Geophys. Res.-Sol. Ea., 126, e2021JB022080,, 2021. 

Gonzalez, J. P., Thomas, J. B., Mazzucchelli, M. L., Angel, R. J., and Alvaro, M.: Experimental evaluation of anisotropic elastic thermobarometry applied to stiff mineral inclusions in soft hosts: First evaluation of zircon inclusions in quartz, Goldschmidt, 10–15 July 2022, Hawai'i, 2022. 

Gonzalez-Platas, J., Alvaro, M., Nestola, F., and Angel, R. J.: EosFit7-GUI: A new GUI tool for equation of state calculations, analyses, and teaching, J. Appl. Crystallogr., 49, 1377–1382,, 2016. 

Hill, R.: Elastic properties of reinforced solids: some theoretical principles, J. Mech. Phys. Solids, 11, 357–372, 1963. 

Holland, T. J. B. and Powell, R.: An improved and extended internally consistent thermodynamic dataset for phases of petrological interest, involving a new equation of state for solids, J. Metamorph. Geol., 29, 333–383,, 2011. 

Izraeli, E., Harris, J., and Navon, O.: Raman barometry of diamond formation, Earth Planet. Sc. Lett., 173, 351–360, 1999. 

Jakubová, P., Kotková, J., Wirth, R., Škoda, R., and Haifler, J.: Morphology and Raman spectral parameters of Bohemian microdiamonds: implications to elastic geothermobarometry, J. Geosci., 67, 253–271,, 2022. 

Kohn, M. J.: “Thermoba-Raman-try”: Calibration of spectroscopic barometers and thermometers for mineral inclusions, Earth Planet. Sc. Lett., 388, 187–196, 2014. 

Korsakov, A. V., Kohn, M. J., and Perraki, M.: Applications of Raman Spectroscopy in Metamorphic Petrology and Tectonics, Elements, 16, 105–110,, 2020. 

Kosminska, K., Spear, F. S., Majka, J., Faehnrich, K., Manecki, M., Piepjohn, K., and Dallmann, W. K.: Deciphering late Devonian-early Carboniferous P-T-t path of mylonitized garnet-mica schists from Prins Karls Forland, Svalbard, J. Metamorph. Geol., 38, 471–493,, 2020. 

Kotková, J., Fedortchouk, Y., Wirth, R., and Whitehouse, M. J.: Metamorphic microdiamond formation is controlled by water activity, phase transitions and temperature, Sci. Rep., 11, 7694,, 2021. 

Kulik, E., Murzin, V., Kawaguchi, S., Nishiyama, N., and Katsura, T.: Thermal expansion of coesite determined by synchrotron powder X-ray diffraction, Phys. Chem. Miner., 45, 873–881,, 2018. 

Mazzucchelli, M. L., Burnley, P., Angel, R. J., Morganti, S., Domeneghetti, M. C., Nestola, F., and Alvaro, M.: Elastic geothermobarometry: corrections for the geometry of the host-inclusion system, Geology, 46, 231–234, 2018. 

Mazzucchelli, M. L., Reali, A., Morganti, S., Angel, R. J., and Alvaro, M.: Elastic geobarometry for anisotropic inclusions in cubic hosts, Lithos, 350–351, 105218,, 2019. 

Mazzucchelli, M. L., Angel, R. J., and Alvaro, M.: EntraPT: an online application for elastic geothermobarometry, Am. Mineral., 106, 829–836,, 2021. 

Milani, S., Nestola, F., Alvaro, M., Pasqual, D., Mazzucchelli, M. L., Domeneghetti, M. C., and Geiger, C.: Diamond–garnet geobarometry: The role of garnet compressibility and expansivity, Lithos, 227, 140–147, 2015. 

Morganti, S., Mazzucchelli, M., Alvaro, M., and Reali, A.: A numerical application of the Eshelby theory for geobarometry of non-ideal host-inclusion systems, Meccanica, 55, 751–764, 2020. 

Murri, M., Gonzalez, J. P., Mazzucchelli, M. L., Prencipe, M., Mihailova, B., Angel, R. J., and Alvaro, M.: The role of symmetry-breaking strains on quartz inclusions in anisotropic hosts: implications for Raman elastic geobarometry, Lithos, 422–423, 106716,, 2022. 

Musiyachenko, K. A., Murri, M., Prencipe, M., Angel, R. J., and Alvaro, M.: A new Grüneisen tensor for rutile and its application to host-inclusion systems, Am. Mineral., 106, 1586,, 2021. 

Myhill, R.: The elastic solid solution model for minerals at high pressures and temperatures, Contrib. Mineral. Petr., 173, 12,, 2018. 

Nestola, F., Nimis, P., Ziberna, L., Longo, M., Marzoli, A., Harris, J. W., Manghnani, M. H., and Fedortchouk, Y.: First crystal-structure determination of olivine in diamond: composition and implications for provenance in the Earth's mantle, Earth Planet. Sc. Lett., 305, 249–255,, 2011. 

Nestola, F., Prencipe, M., Nimis, P., Sgreva, N., Peritt, S. H., Chinn, I. L., and Zaffirio, G.: Toward a robust elastic geobarometry of kyanite inclusions in eclogitic diamonds, J. Geophys. Res.-Sol. Ea., 123, 6411–6423,, 2018. 

Ni, P., Zhang, Y., and Guan, Y.: Volatile loss during homogenization of lunar melt inclusions, Earth Planet. Sc. Lett., 478, 214–224,, 2017. 

Nimis, P., Alvaro, M., Nestola, F., Angel, R. J., Marquardt, K., Rustioni, G., Harris, J. W., and Marone, F.: First evidence of hydrous silicic fluid films around solid inclusions in gem-quality diamonds, Lithos, 260, 384–389,, 2016. 

Nimis, P., Angel, R. J., Alvaro, M., Nestola, F., Harris, J. W., Casati, N., and Marone, F.: Crystallographic orientations of magnesiochromite inclusions in diamonds: what do they tell us?, Contrib. Mineral. Petr., 174, 29,, 2019. 

Perraki, M., Proyer, A., Mposkos, E., Kaindl, R., and Hoinkes, G.: Raman micro-spectroscopy on diamond, graphite and other carbon polymorphs from the ultrahigh-pressure metamorphic Kimi Complex of the Rhodope Metamorphic Province, NE Greece, Earth Planet. Sc. Lett., 241, 672–685,, 2006. 

Perrillat, J. P., Daniel, I., Lardeaux, J. M., and Cardon, H.: Kinetics of the quartz-coesite transition: application to the exhumation of ultrahigh-pressure rocks, J. Petrol., 44, 773–788, 2003. 

Pitzer, K. S. and Sterner, S. M.: Equations of state valid continuously from zero to extreme pressures for H2O and CO2, J. Chem. Phys., 101, 3111–3116,, 1994. 

Pollard, D. D. and Fletcher, R. C.: Fundamentals of Structural Geology, Cambridge University Press, Cambridge, UK, ISBN 978-0521839273, 2006.  

Rodriguez-Carvajal, J. and Gonzalez-Platas, J.: Crystallographic Fortran 90 Modules Library (CrysFML): a simple toolbox for crystallographic computing programs, IUCr Computing Commission Newsletter, 1, 50–58, 2003 (code available at:, last access: 20 June 2023). 

Roedder, E. and Bodnar, R. J.: Geologic pressure determinations from fluid inclusion studies, Annu. Rev. Earth Planet. Sc., 8, 263–301, 1980. 

Rosenfeld, J. L. and Chase, A. B.: Pressure and temperature of crystallization from elastic effects around solid inclusion minerals?, Am. J. Sci., 259, 519–541, 1961. 

Schmalholz, S. M., Moulas, E., Plümper, O., Myasnikov, A. V., and Podladchikov, Y. Y.: 2D hydro-mechanical-chemical modeling of (de)hydration reactions in deforming heterogeneous rock: the periclase-brucite model reaction, Geochem. Geophy. Geosy., 21, e2020GC009351,, 2020. 

Sobolev, N. V., Fursenko, B. A., Goryainov, S. V., Shu, J. F., Hemley, R. J., Mao, H. K., and Boyd, F. R.: Fossilized high pressure from the Earth's deep interior: The coesite-in-diamond barometer, P. Natl. Acad. Sci. USA, 97, 11875–11879, 2000. 

Stixrude, L. and Lithgow-Bertelloni, C.: Thermal expansivity, heat capacity and bulk modulus of the mantle, Geophys. J. Int., 228, 1119–1149,, 2022. 

Ye, K., Liou, J. B., Cong, B. L., and Maruyama, S.: Overpressures induced by coesite-quartz transition in zircon, Am. Mineral., 86, 1151–1155, 2001. 

Zhang, Y.: Mechanical and phase equilibria in inclusion–host systems, Earth Planet. Sc. Lett., 157, 209–222,, 1998. 

Zhang, Z. and Duan, Z.: Prediction of the PVT properties of water over wide range of temperatures and pressures from molecular dynamics simulation, Phys. Earth Planet. In., 149, 335–354, 2005. 

Zhong, X., Moulas, E., and Tajčmanová, L.: Tiny timekeepers witnessing high-rate exhumation processes, Sci. Rep., 8, 2234,, 2018. 

Zhong, X., Dabrowski, M., and Jamtveit, B.: Analytical solution for residual stress and strain preserved in anisotropic inclusion entrapped in an isotropic host, Solid Earth, 12, 817–833,, 2021. 

Short summary
We have developed the thermodynamic theory of the properties of inclusions consisting of more than one phase, including inclusions containing solids plus a fluid. We present a software utility that enables for the first time the entrapment conditions of multiphase inclusions to be determined from the measurement of their internal pressure when that is measured in a laboratory.