Post-entrapment modification of residual inclusion pressure and its implications for Raman elastic thermobarometry

Residual pressure can be preserved in mineral inclusions, e.g. quartz-in-garnet, after exhumation due to differential expansion between inclusion and host crystals. Raman spectroscopy has been applied to infer the residual pressure and provides information on the entrapment temperature and pressure conditions. However, the amount of residual pressure relaxation cannot be directly measured. An underestimation or overestimation of residual pressure may lead to significant errors between calculated and actual entrapment pressure. This study focuses on three mechanisms responsible for the residual pressure modification: (1) viscous creep; (2) plastic yield; (3) proximity of inclusion to the thin-section surface. Criteria are provided to quantify how much of the expected residual pressure is modified due to these three mechanisms. An analytical solution is introduced to demonstrate the effect of inclusion depth on the residual pressure field when the inclusion is close to the thin-section surface. It is shown that for a quartz-in-garnet system, the distance between the thin-section surface and inclusion centre needs to be at least 3 times the inclusion radius to avoid pressure release. In terms of viscous creep, representative case studies on a quartz-in-garnet system show that viscous relaxation may occur from temperatures as low as 600–700 ^∘C depending on the particular pressure–temperature (P–T) path and various garnet compositions. For quartz entrapped along the prograde P–T path and subject to viscous relaxation at peak T above 600–700 ^∘C, its residual pressure after exhumation may be higher than predicted from its true entrapment conditions. Moreover, such a viscous resetting effect may introduce apparent overstepping of garnet nucleation that is not related to reaction affinity.

1 Introduction

During metamorphism, the growth of porphyroblasts often results in the entrapment of solid or fluid inclusions, which then provide important information about the rock's history (e.g. Farber et al., 2014; Yardley and Bodnar, 2014; Ferrero and Angel, 2018). Due to the differences in the elastic compressibility and thermal expansion coefficient between the inclusion and host, residual inclusion pressures may be preserved after exhumation (e.g. Rosenfeld and Chase, 1961; Gillet et al., 1984; Zhang, 1998; Angel et al., 2015). The residual pressure can be inferred by Raman shift based on experimental calibrations, e.g. quartz inclusions (Liu and Mernagh, 1992; Schmidt and Ziemann, 2000). This allows the application of Raman thermobarometry to infer the entrapment pressure and temperature (P–T) conditions (e.g. Ashley et al., 2014; Bayet et al., 2018; Enami et al., 2007; Izraeli et al., 1999; Kohn, 2014; Spear et al., 2014; Taguchi et al., 2019a; Zhong et al., 2019a). Existing models that link residual pressure and entrapment P–T conditions are based on elastic rheology and often assume an infinite host radius (Rosenfeld and Chase, 1961; Van Der Molen and Van Roermund, 1986; Guiraud and Powell, 2006; Angel et al., 2017b). Despite these simplifications, recent experimental works have been successfully performed to compare the measured residual pressure with modelled residual pressure under well-controlled P–T conditions for synthetic samples with a quartz-in-garnet system (Thomas and Spear, 2018; Bonazzi et al., 2019).

Although many studies using Raman spectroscopy reported maximal residual pressure close to the predictions from elastic models (e.g. Ashley et al., 2014; Enami et al., 2007; Zhong et al., 2019a), a large amount of inclusion pressure estimates are lower than theoretically predicted by the elastic model (Korsakov et al., 2009; Kouketsu et al., 2016; Yamamoto et al., 2002). The modification of inclusion pressure can be due to various reasons and a systematic investigation is critical to better apply Raman thermobarometry to natural samples. Meanwhile, Raman thermobarometry has been employed to investigate the amount of overstepping for garnet growth by comparing the P–T constraints from phase equilibria and elastic thermobarometry (Spear et al., 2014; Castro and Spear, 2017; Wolfe and Spear, 2017). Particularly, when comparing the determined paleopressures based on phase equilibria and elastic barometry using a quartz-in-garnet system, careful examination of the amount of residual quartz pressure modification due to the creep of garnet host becomes critical.

When a mineral inclusion maintains residual pressure, differential stress is concomitantly developed around the inclusion on the host side to maintain mechanical equilibrium (e.g. Zhang, 1998; Tajčmanová et al., 2014). The host mineral may experience viscous creep, which is manifested by the dislocation structures (Chen et al., 1996; Yamamoto et al., 2002; Taguchi et al., 2019b). Furthermore, the host mineral may also experience rate-independent plastic yield when the differential stress exceeds the yield criterion (e.g. Zhang, 1998). In the mechanics literature, plastic deformation is commonly considered as any inelastic deformation (time-dependent and time-independent) (e.g. Kachanov, 1971). In this work, we distinguish between viscous creep, i.e. the rate-dependent inelastic deformation, and the rate-independent plastic flow. Mechanical models show that both viscous creep (dislocation or diffusion creep of host) and plastic flow during decompression and cooling can cause a significant inclusion pressure drop (Dabrowski et al., 2015; Zhang, 1998). This would lead to an underestimate of residual inclusion pressure (Zhong et al., 2018) (Fig. 1). Meanwhile, during the thin-section preparation, mineral inclusions are positioned into proximity towards the thin-section surface (Fig. 1). The thin-section surface is stress-free and may elastically release the residual pressure (Mindlin and Cheng, 1950; Seo and Mura, 1979; Zhong et al., 2019b). It is of petrological interest to study how deep the inclusion needs to be in order to preserve the residual pressure. Experimental works and numerical simulations with the finite-element method have been performed to test the safe inclusion depth (inclusion radius less than one-third of host radius) so that the residual inclusion pressure can be preserved for the application of Raman barometry (Campomenosi et al., 2018; Mazzucchelli et al., 2018).

Figure 1 Schematic illustration for the residual pressure. The grey and black curves are retrograde P–T paths for host and inclusion, respectively. Pressure drop is possibly due to the following reasons: (1) viscous relaxation preferentially occurs at high-temperature conditions; (2) plastic yield commonly occurs at low confining pressures where residual pressure is high; (3) thin-section preparation that drives inclusion close to the thin-section surface. Note that this illustration refers to systems where the inclusion is elastically softer than its host (e.g. quartz in garnet).

In this contribution, we systematically investigate the following mechanisms for residual inclusion pressure modification: (1) viscous creep of the host materials, (2) plastic yield within the host, and (3) pressure release due to the proximity of inclusion towards the thin-section surface. For the first and second purposes, a 1-D visco-elasto-plastic mechanical model is developed in a radially symmetric spherical coordinate frame. The derived system of equations is non-dimensionalized to extract the key parameters that control the amount of viscous relaxation and plastic yield. For the third mechanism, a simple analytical solution for the residual inclusion pressure field close to the thin-section surface is introduced based on the existing work of Seo and Mura (1979). The analytical solution demonstrates the effect of the inclusion depth that controls the amount of pressure release. This solution applies to the case where the inclusion possesses the same elastic moduli as the host. The inclusion is initially subject to an arbitrary hydrostatic pressure in an infinite host and its pressure is released as it moves towards a stress-free thin-section surface. In comparison, for a natural quartz-in-garnet system, numerical solutions are applied to investigate the safe distance that causes negligible pressure release due to the presence of the thin-section surface (stress-free boundary). In this study, both inclusion and host are treated as elastically isotropic as an assumption to put focus on the effect of these three mechanisms on preserved residual pressure. The effects of elastic anisotropy for commonly encountered quartz inclusion have been studied experimentally and theoretically by, e.g., Murri et al. (2018) and Campomenosi et al. (2018) and are discussed at the end.

2 1-D mechanical model with visco-elasto-plastic rheology

2.1 Governing equations

We develop a 1-D mechanical model with spherical symmetry that is based on the conservation of mass and momentum. In 1-D radially symmetric spherical coordinate frame, mechanical equilibrium is expressed as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial {\mathit{\tau }}_{rr}}{\partial r}}+{\displaystyle \frac{\mathrm{3}{\mathit{\tau }}_{rr}}{r}}-{\displaystyle \frac{\partial P}{\partial r}}=\mathrm{0},\end{array}$$

where τ_rr is the radial component of deviatoric stress (Pa), P is pressure (Pa), and r is radial coordinate (m). We apply the Maxwell visco-elasto-plastic rheology as follows:

$$\begin{array}{}\text{(2)}& {\dot{e}}_{rr}={\dot{e}}_{rr}^e+{\dot{e}}_{rr}^v+{\dot{e}}_{rr}^p,\end{array}$$

where ${\dot{e}}_{rr}$ is the radial components of the deviatoric strain rate (s⁻¹) composed of elastic, viscous (rate-dependent) and plastic (rate-independent) parts. The elastic and viscous strain rates are expressed as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\dot{e}}_{rr}^e={\displaystyle \frac{{\dot{\mathit{\tau }}}_{rr}}{\mathrm{2}G}},\\ & {\dot{e}}_{rr}^v={\displaystyle \frac{{\mathit{\tau }}_{rr}}{\mathrm{2}\mathit{\eta }}},\end{array}\end{array}$$

where the dot above ${\dot{\mathit{\tau }}}_{rr}$ denotes time derivative, G is shear modulus (Pa), and η is viscosity (Pa s). The non-Newtonian (effective) viscosity is expressed as follows:

$$\begin{array}{}\text{(4)}& \mathit{\eta }=A\mathrm{|}{\mathit{\tau }}_{rr}{\mathrm{|}}^{\mathrm{1}-n},\end{array}$$

where A is the temperature-dependent pre-factor and n is the stress exponent (e.g. Dabrowski et al., 2015, Eq. 10). The plastic strain rate is obtained by using the Tresca yield criterion (e.g. Kachanov, 1971):

$$\begin{array}{}\text{(5)}& F=\left|{\mathit{\tau }}_{rr}-{\mathit{\tau }}_{tt}\right|-C,\end{array}$$

where C is plastic yield strength (Pa) and τ_tt is the tangential component of deviatoric stress. Due to spherical symmetry, we have ${\mathit{\tau }}_{tt}=-\mathrm{1}/\mathrm{2}{\mathit{\tau }}_{rr}$. Applying the plastic flow law (e.g. Vermeer and De Borst, 1984), we get

$$\begin{array}{}\text{(6)}& {\dot{e}}_{rr}^p=\mathit{\lambda }{\displaystyle \frac{\partial F}{\partial {\mathit{\tau }}_{rr}}}=\mathit{\lambda }\mathrm{sgn}\left({\mathit{\tau }}_{rr}\right)=\mathit{\lambda }\mathit{\delta },\left\{\begin{array}{l}\mathit{\lambda }=\mathrm{0}\mathrm{for}F\le \mathrm{0}\\ \mathit{\lambda }\ne \mathrm{0}\mathrm{for}F>\mathrm{0}\end{array}\right.,\end{array}$$

where λ is the plastic multiplier (s⁻¹), which provides the amount of plastic strain (rate) that guarantees the yield criterion is not exceeded, and δ is the sign of τ_rr. For isotropic materials, the pressure (negative mean stress) can be expressed as a function of volume and temperature via the equation of state (EoS), and its time derivative is as follows:

$$\begin{array}{}\text{(7)}& \dot{P}=-{\dot{\mathit{\epsilon }}}_{kk}/\mathit{\beta }+\mathit{\alpha }\dot{T}/\mathit{\beta },\end{array}$$

where β is compressibility (Pa⁻¹), α is the thermal expansion coefficient (K⁻¹), and $\dot{T}$ is the rate of temperature change (K s⁻¹). Temperature is treated as homogeneous within the inclusion–host system. Einstein summation is used here for the volumetric strain rate (${\dot{\mathit{\epsilon }}}_{kk}={\dot{\mathit{\epsilon }}}_{rr}+\mathrm{2}{\dot{\mathit{\epsilon }}}_{tt}$, where due to spherical symmetry the two tangential strain rates are equal). No viscous or plastic volumetric strain is considered. This assumption is a good approximation for non-porous, crystalline materials (e.g. Moulas et al., 2019).

Substituting Eqs. (3) and (6) into Eq. (2) and applying first-order finite difference in time to Eqs. (2) and (7) (i.e. ${\dot{\mathit{\tau }}}_{rr}=\frac{{\mathit{\tau }}_{rr}-{\mathit{\tau }}_{rr}^o}{\mathrm{\Delta }t}$ and $\dot{P}=\frac{P-P^o}{\mathrm{\Delta }t}$), we can explicitly express τ_rr and P as

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{\tau }}_{rr}=\mathrm{2}\mathit{\eta }Z{\dot{e}}_{rr}+\left(\mathrm{1}-Z\right){\mathit{\tau }}_{rr}^o-\mathrm{2}\mathit{\eta }Z\mathit{\lambda }\mathit{\delta },\text{(9)}& {\displaystyle }P=P^o-\mathrm{\Delta }t{\dot{\mathit{\epsilon }}}_{kk}/\mathit{\beta }+\mathit{\alpha }\mathrm{\Delta }t\dot{T}/\mathit{\beta },\end{array}$$

where $Z=\frac{\mathrm{\Delta }tG}{\mathrm{\Delta }tG+\mathit{\eta }}$ is the viscoelastic coefficient, Δt is the time increment, ${\mathit{\tau }}_{rr}^o$ is the radial component of deviatoric stress in the previous time step, P^o is the pressure in the previous time step. If the yield criterion in Eq. (5) is exceeded (F>0), the plastic multiplier must be adjusted to drive F to zero. This can be achieved by substituting the deviatoric stress Eq. (8) into Eq. (5), and we let F=0. Therefore, we obtain λ as follows:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathit{\lambda }& =\mathit{\delta }{\dot{e}}_{rr}+{\displaystyle \frac{\left(\mathrm{1}-Z\right)\mathit{\delta }}{\mathrm{2}\mathit{\eta }Z}}{\mathit{\tau }}_{rr}^o-{\displaystyle \frac{\mathit{\delta }C}{\mathrm{3}\mathit{\eta }Z}},\\ & \mathrm{if}F>\mathrm{0}(\mathrm{otherwise}\mathit{\lambda }=\mathrm{0}).\end{array}\end{array}$$

2.2 Non-dimensionalization

The variables in the above equations can be scaled to derive non-dimensional parameters for better understanding the behaviour of the inclusion–host system. This is done by choosing the following six parameters to non-dimensionalize the system of equations: the temperature drop in the host-inclusion system ΔT, the inclusion radius R, the time of the P–T path t^∗, the host's viscosity pre-factor A_h, the host's plastic yield strength C_h, and the expected pressure perturbation P_exp that is given as follows:

$$\begin{array}{}\text{(11)}& P_{\exp }={\displaystyle \frac{\mathrm{\Delta }P({\mathit{\beta }}_{\mathrm{i}}-{\mathit{\beta }}_{\mathrm{h}})-\mathrm{\Delta }T({\mathit{\alpha }}_{\mathrm{i}}-{\mathit{\alpha }}_{\mathrm{h}})}{{\mathit{\beta }}_{\mathrm{i}}+\mathrm{3}/\mathrm{4}{G}_{\mathrm{h}}}},\end{array}$$

where ΔP and ΔT are the confining pressure and temperature drops from entrapment to the Earth's surface, β_i and β_h are the compressibility of inclusion and host, α_i and α_h are the thermal expansion coefficients of inclusion and host, and G_h is the shear modulus of host. The number P_exp is the expected amount of residual inclusion pressure after exhumation assuming linear thermo-elasticity and an infinite host (Zhang, 1998). It is noted that this is not the actual final residual inclusion pressure but merely a scale to non-dimensionalize the stress (and pressure). By choosing P_exp as the stress scale, the inclusion residual pressure is expected to be between 0 and 1 for a case of cooling and decompression. This pressure scale allows convenient quantification for the amount of pressure modification due to viscous creep and plastic yield. The involved physical variables are scaled as shown below:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& r=R\stackrel{\mathrm{‾}}{r},\\ & \mathit{\beta }={\displaystyle \frac{\mathrm{1}}{P_{\exp }}}\stackrel{\mathrm{‾}}{\mathit{\beta }},\\ & G=P_{\exp }\stackrel{\mathrm{‾}}{G},\\ & \mathit{\alpha }={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }T}}\stackrel{\mathrm{‾}}{\mathit{\alpha }},\\ & P=P_{\exp }\stackrel{\mathrm{‾}}{P},\\ & \dot{T}={\displaystyle \frac{\mathrm{\Delta }T}{t^{\ast }}}\stackrel{\mathrm{‾}}{\dot{T}},\\ & {\mathit{\tau }}_{rr}=P_{\exp }\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}},\\ & C=C_{\mathrm{h}}\stackrel{\mathrm{‾}}{C},\\ & \mathit{\eta }=P_{\exp }{t}^{\ast }\stackrel{\mathrm{‾}}{\mathit{\eta }},\\ & F=C_{\mathrm{h}}\stackrel{\mathrm{‾}}{F},\\ & \mathrm{\Delta }t=t^{\ast }\stackrel{\mathrm{‾}}{\mathrm{\Delta }t},\\ & A=A_{\mathrm{h}}\stackrel{\mathrm{‾}}{A},\\ & \mathit{\lambda }={\displaystyle \frac{\mathrm{1}}{t^{\ast }}}\stackrel{\mathrm{‾}}{\mathit{\lambda }},\\ & v_r={\displaystyle \frac{R}{t^{\ast }}}\stackrel{\mathrm{‾}}{v_r},\end{array}\end{array}$$

where the overhead bars indicate dimensionless properties. Substituting these scaling equations into Eqs. (1), (8), and (9), we get

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}}}{\partial \stackrel{\mathrm{‾}}{r}}}+{\displaystyle \frac{\mathrm{3}\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}}}{\stackrel{\mathrm{‾}}{r}}}-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{P}}{\partial \stackrel{\mathrm{‾}}{r}}}=\mathrm{0},\text{(14)}& {\displaystyle }\stackrel{\mathrm{‾}}{P}=\stackrel{\mathrm{‾}}{P^o}+{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathit{\beta }}}}\left[-\stackrel{\mathrm{‾}}{\mathrm{\Delta }t}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{r}}^{\mathrm{2}}\stackrel{\mathrm{‾}}{v_r}}{{\stackrel{\mathrm{‾}}{r}}^{\mathrm{2}}\partial \stackrel{\mathrm{‾}}{r}}}+\stackrel{\mathrm{‾}}{\mathit{\alpha }}\stackrel{\mathrm{‾}}{\dot{T}}\right],\text{(15)}& {\displaystyle }\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}}={\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\stackrel{\mathrm{‾}}{\mathit{\eta }}\stackrel{\mathrm{‾}}{Z}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{v_r}}{\partial \stackrel{\mathrm{‾}}{r}}}-{\displaystyle \frac{\stackrel{\mathrm{‾}}{v_r}}{\stackrel{\mathrm{‾}}{r}}}\right)+\left(\mathrm{1}-\stackrel{\mathrm{‾}}{Z}\right)\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}^o}-\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\eta }}\stackrel{\mathrm{‾}}{\mathit{\lambda }}\mathit{\delta }\stackrel{\mathrm{‾}}{Z},\end{array}$$

where dimensionless viscosity, viscoelastic coefficient, and plastic multiplier are expressed as

$$\begin{array}{}\text{(16)}& {\displaystyle }\stackrel{\mathrm{‾}}{\mathit{\eta }}=\mathit{\text{De}}\cdot \stackrel{\mathrm{‾}}{A}\stackrel{\mathrm{‾}}{\mathrm{|}{\mathit{\tau }}_{rr}}{\mathrm{|}}^{\mathrm{1}-n},\text{(17)}& {\displaystyle }\stackrel{\mathrm{‾}}{Z}={\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathrm{\Delta }t}\stackrel{\mathrm{‾}}{G}}{\stackrel{\mathrm{‾}}{\mathrm{\Delta }t}\stackrel{\mathrm{‾}}{G}+\stackrel{\mathrm{‾}}{\mathit{\eta }}}},\text{(18)}& {\displaystyle }\begin{array}{rl}\stackrel{\mathrm{‾}}{\mathit{\lambda }}& ={\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\mathit{\delta }\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{v_r}}{\partial \stackrel{\mathrm{‾}}{r}}}-{\displaystyle \frac{\stackrel{\mathrm{‾}}{v_r}}{\stackrel{\mathrm{‾}}{r}}}\right)+{\displaystyle \frac{\left(\mathrm{1}-\stackrel{\mathrm{‾}}{Z}\right)\mathit{\delta }}{\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\eta }}\stackrel{\mathrm{‾}}{Z}}}\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}^o}\\ & -C^{\ast }{\displaystyle \frac{\stackrel{\mathrm{‾}}{C}}{\mathrm{3}\stackrel{\mathrm{‾}}{\mathit{\eta }}\stackrel{\mathrm{‾}}{Z}}},\mathrm{if}\stackrel{\mathrm{‾}}{F}>\mathrm{0}.\end{array}\end{array}$$

Two dominant dimensionless numbers emerge after non-dimensionalization. They are the Deborah number De and dimensionless yield strength C^∗ defined as follows:

$$\begin{array}{}\text{(19)}& {\displaystyle }\mathit{\text{De}}={\displaystyle \frac{A_{\mathrm{h}}/P_{\exp }^n}{t^{\ast }}},\text{(20)}& {\displaystyle }{C}^{\ast }={\displaystyle \frac{C_{\mathrm{h}}}{P_{\exp }}}.\end{array}$$

The Deborah number (De) is the ratio between the characteristic viscous relaxation time ($A_{\mathrm{h}}/P_{\exp }^n$) and model duration (t^∗) (Reiner, 1964). If De>1, the system behaves in an elastic manner, and if $\mathit{\text{De}}<\mathrm{1},$ viscous creep becomes important. The pre-factor of viscosity is temperature-dependent. By choosing the pre-factor A_h at peak temperature, one can directly use De to estimate the maximal amount of viscous relaxation. This is especially suitable for the process of isothermal decompression in many high-pressure (HP) rocks.

The dimensionless yield strength C^∗ characterizes the ability of a host mineral to plastically yield and a high C^∗ implies that the material is less prone to plastic yield given the amount of residual inclusion pressure P_exp. The viscosity of different mineral phases may vary by orders of magnitude, and the plastic yield strength of different minerals may also vary by several factors. Therefore, these two dimensionless numbers have a dominant effect on the amount of inclusion pressure modification due to viscous relaxation and plastic yield.

2.3 The numerical approach for the visco-elasto-plastic model

The dimensionless viscosity (Eq. 16), viscoelastic coefficient (Eq. 17), and plastic multiplier (Eq. 18) can be substituted into pressure equation (Eq. 14) and deviatoric stress equation (Eq. 15). Together with mechanical equilibrium equation (Eq. 13), they form a system of three equations with three unknowns, namely $\stackrel{\mathrm{‾}}{v_r}$, $\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}}$, and $\stackrel{\mathrm{‾}}{P}$. The numerical model is based on a finite-difference scheme over a 1-D staggered grid (on the numerical stencil, see, e.g., Gerya, 2010, chap. 7). The initial pressure $\stackrel{\mathrm{‾}}{P}$ is set at the beginning of the numerical model. If upon entrapment, the inclusion and host possess the same hydrostatic pressure, the deviatoric stress $\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}}$ is zero in the inclusion and host. If pressure heterogeneity exists upon entrapment, the deviatoric stress of the host ($\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}}$) needs to be precalculated with the elastic model $\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}}=-(P_{\mathrm{inc}}-P_{\mathrm{host}})/{\stackrel{\mathrm{‾}}{r}}^{\mathrm{3}}$ to ensure that mechanical equilibrium is satisfied at the beginning of the model (P_inc is the initial inclusion pressure and P_host is the initial host pressure).

For predefined P–T path, the inclusion–host system is loaded by the increment of confining pressure and temperature. Both temperature and far-field pressure are obtained directly from the P–T–t path as prescribed. Temperature is treated as homogeneous in the model and the new pressure is set as the outer boundary value. Because viscosity, the viscoelastic coefficient, and the plastic multiplier are functions of deviatoric stress (viscosity is also a function of temperature as prescribed by the P–T path), the system of the mechanical equations is nonlinear. We solve for the three variables ($\stackrel{\mathrm{‾}}{v_r}$, $\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{rr}}$, and $\stackrel{\mathrm{‾}}{P}$) using an iterative (Picard) method. Within the iteration loop, an elastic test stress is first evaluated by letting $\stackrel{\mathrm{‾}}{\mathit{\lambda }}=\mathrm{0}$ so that the prediction for the yield function $\stackrel{\mathrm{‾}}{F}$ is obtained. If $\stackrel{\mathrm{‾}}{F}<\mathrm{0}$, no plastic yield occurs and $\stackrel{\mathrm{‾}}{\mathit{\lambda }}$ remains zero. Otherwise the prediction of the yield function is positive and $\stackrel{\mathrm{‾}}{\mathit{\lambda }}$ is computed based on Eq. (18). The calculated $\stackrel{\mathrm{‾}}{\mathit{\lambda }}$ is then substituted back into Eq. (15) to adjust the amount of plastic strain rate. This will drive $\stackrel{\mathrm{‾}}{F}$ to zero (on the plastic yield surface). After the Picard iteration loop, the residuals of the three Eqs. (13), (14), and (15) are minimized below ca. 10⁻¹².

The elastic moduli are updated based on pressure and temperature conditions from tabulated lookup tables within the iteration. The lookup tables are pre-computed based on EoS. We used the EoS for quartz crystal from Angel et al. (2017a) and the EoS for pyrope, grossular, and almandine crystals from Milani et al. (2015). The detailed expressions of EoS can be found in the EoSFit7c software documentation (Angel et al., 2014). The EoS for spessartine is from Gréaux and Yamada (2014). The compressibility and thermal expansion coefficient for garnet are averaged based on the molar percentage of garnet endmembers. The shear moduli of garnet endmembers are from Bass (1995). The numerical model has been benchmarked using the analytical solution with elastic, non-Newtonian viscous rheology in Zhong et al. (2018). The numerical benchmark for elasto-plastic rheology is performed by using the analytical solution of Zhang (1998) (see Supplement).

3 Inclusion pressure modification due to visco-plastic deformation of host

3.1 Residual pressure affected by viscous or plastic flow

The solutions of the system of equation (Eqs. 13, 14, 15) are obtained using the elasticity of a quartz-in-garnet system. The host radius is set to be 50 times the inclusion radius to make boundary effects negligible. Temperature is treated as homogeneous in space. At the beginning of the model, a pressure perturbation within the inclusion is prescribed, and the far-field host maintains zero confining pressure. The pre-factor of viscosity is fixed as temperature does not vary in this case. The amount of inclusion pressure relaxation is systematically investigated for the two inelastic deformation mechanisms (i.e. viscous creep and plastic yield) as a function of De and C^∗. The results are shown in Fig. 2 with the purpose of systematically demonstrating how much the initially prescribed residual pressure can be reduced due to viscous creep and plastic yield. This diagram may assist petrological investigations because De and C^∗ can be evaluated based on experimental rock deformation data for different minerals. The Deborah number can be evaluated using the experimental flow law of a single crystal, e.g. garnet (Karato et al., 1995; Wang and Ji, 1999), as given in the next section. The plastic yield strength is evaluated using microhardness test data (see Discussion below for details). The thickness of the plastic yield region is plotted as contours. The thick grey contour represents the onset of plastic yield starting from the inclusion–host interface and propagating towards the host side (Fig. 2). Based on the amount of inclusion pressure relaxation, three regimes are distinguished.

Figure 2 Inclusion pressure as a function of the Deborah number and dimensionless yield strength given different stress exponents. The contours denote the radius of plastic yield region R_yield scaled by inclusion radius. The thick grey contour represents the onset of plastic yield. Three regimes are labelled: (1) elastic (De>1, $C^{\ast }>\mathrm{1}$); (2) viscous (De<1 and C^∗ is above the onset of plastic yield); (3) plastic ($C^{\ast }<\mathrm{1}$, De is above the onset of plastic yield). To obtain the results, a residual pressure is prescribed at the beginning and the confining pressure and temperature are fixed, i.e. no temporal variations in P–T conditions.

The elastic regime takes place when De and C^∗ are higher than 1. Under these circumstances, no viscous relaxation and plastic yield occurs. The residual inclusion pressure is close to the expected residual pressure (P_inc≈P_exp). This regime is the most suitable for the application of Raman thermobarometry.

The viscous regime dominates when De is lower than 1, and C^∗ is above the plastic onset shown by the thick grey contour. In this case, the main mechanism responsible for the inclusion pressure relaxation is the viscous creep. The effect of stress exponent on the amount of viscous relaxation is also significant. In general, a higher stress exponent delays pressure relaxation (cf. Dabrowski et al., 2015). As the viscosity of natural minerals is low at high-temperature conditions, the viscous regime may be reached at a high temperature, which leads to the relaxation of residual pressure.

The plastic regime prevails when C^∗ is lower than 1 and De is above 1. Under this circumstance, the residual pressure is not significantly relaxed by viscous creep, but by plastic yield. In general, the radius of the plastic yield region is positively related to the amount of residual pressure release.

It is noted here that although viscous relaxation and plastic yield of the host mineral have the same effect in reducing the residual inclusion pressure after exhumation, there is a fundamental difference between them. Viscous relaxation is time-dependent (De includes time), which means that the residual pressure will vanish given an infinite amount of time. Plastic yield refers to a time-independent process, and it will only limit the amount of deviatoric stress supported by the host mineral. If the yield criterion is reached, plastic strain (rate) in the host immediately occurs, which leads to the inclusion pressure release. Both viscous relaxation and plastic yield are irreversible; as a consequence, if the inclusion-host system is placed back in the original entrapment condition, the stress state would be different.

3.2 Viscous relaxation of garnet host

Assuming that the thin-section surface is sufficiently far away from a quartz inclusion and no plastic yield occurs around the quartz inclusion, only viscous creep may contribute to the modification of residual pressure. Here, we show the effect of viscous relaxation, particularly influenced by the temperature, on the preserved residual pressure. Using De as a criterion to estimate the amount of viscous relaxation, we show the relationship between temperature, inclusion pressure, and relaxation time given De=1 (see Eq. 19) in Fig. 3. The flow law of garnet from Wang and Ji (1999) is applied. The pre-factor A of the effective viscosity (Eq. 4) is as follows:

$$\begin{array}{}\text{(21)}& A={\displaystyle \frac{G^n}{\mathrm{2}B}}\exp \left({\displaystyle \frac{g\cdot T_{\mathrm{m}}}{T}}\right),\end{array}$$

where B=exp (40.1) in s⁻¹ and g=32. The stress exponent n=3. Geometric correction based on experimental set-up (simple or pure shear) was not applied. The melting point T_m of pyrope-rich garnet, grossular, and spessartine is from the Table 1 in Karato et al. (1995). For almost pure almandine, the garnet melting point is found to be 1588 K from Mohawk Garnet Inc., which is slightly higher than 1570 K for almandine-rich (Alm_0.68Prp_0.20Grs_0.12) garnet in Karato et al. (1995).

Figure 3 Viscous relaxation time (in years) of different garnet hosts as functions of temperature and inclusion overpressure. The viscous relaxation time is calculated based on the expression of the Deborah number (De=1) in Eq. (19). The viscosity's pre-factor is T-dependent and is obtained using the flow law from Wang and Ji (1999). The melting temperature is from Karato et al. (1995) (the melting temperature of almost pure almandine is taken from the data of Mohawk Garnet Inc. to be 1588 K). Shear modulus is from Bass (1995). The flow law is given in the main text (Wang and Ji, 1999). The four garnet endmembers are almandine (Alm), grossular (Grs), pyrope (Prp), and spessartine (Sps).

Table 1 Averaged plastic strength from microhardness tests for some minerals at room conditions. Strength is converted from microhardness based on $C_{\mathrm{h}}=H/C_{\mathrm{g}}$, where the geometry constant C_g is taken as 3. Raw data are dependent on crystallographic orientation, composition, and applied load that are examined in some of the involved references.

¹ Data reported in Whitney et al. (2007). ² Data reported in Wong and Bradt (1992). The reported data for calcite and dolomite are averaged from the applied load and azimuthal angle from [$\mathrm{10}\stackrel{\mathrm{‾}}{\mathrm{1}}\stackrel{\mathrm{‾}}{\mathrm{1}}$]. ³ Data reported in Smedskjaer et al. (2008). ⁴ Data reported in Yuan et al. (2017) ⁵ Data reported in Dekker and Rieck (1974). The reported data are averaged from the applied load at [110] and [100].

As an example (Fig. 3), for a quartz inclusion possessing 0.5 GPa residual pressure maintained at 650 ^∘C, significant viscous relaxation is expected during 1 Myr for almandine-rich garnet host based on the applied flow law. This temperature becomes higher (700 ^∘C) for pyrope-rich garnet. If the residual pressure is used to recover the entrapment pressure given that the temperature experienced by the garnet-host system was higher than 650–700 ^∘C, an underestimate of the entrapment pressure may potentially occur.

Figure 4 (a) Synthetic retrograde P–T paths from eclogite-facies metamorphic conditions. The quartz inclusions are entrapped within almandine at different peak P–T conditions along the same isomeke; thus a purely-elastic model would predict the same value for the residual inclusion P. Due to viscous relaxation, the residual P is lower than the pressure predicted by an elastic model. In (b), the apparent entrapment P is calculated based on the relaxed residual inclusion pressure given different entrapment T along the elastic isomeke that is given in (a). Pressure relaxation is manifested by lower values of apparent entrapment P and it becomes more significant if the host experiences high temperatures with time. Panels (c) and (d) are the same plots for amphibolite-facies entrapment conditions. The amount of viscous relaxation is less compared to eclogite facies due to the lower magnitude of inclusion overpressure and the stress-dependent viscosity of garnet host. Pure almandine garnet is used as host, and its flow law is from Wang and Ji (1999).

In Fig. 4, synthetic retrograde P–T paths from eclogite and amphibolite-facies metamorphic conditions are prescribed with different peak temperatures. The entrapment P–T conditions for the three synthetic P–T paths are along an elastic isomeke, which is the isopleth where no relative elastic interaction exists between inclusion and host. Thus, the residual inclusion pressure would be the same if the inclusions were entrapped along the same isomeke and they were purely elastic. By involving the viscous rheology of the garnet host, different residual inclusion pressures are predicted. For the P–T path starting at 800 ^∘C, 2 GPa, the quartz inclusion pressure is predicted to be less than 0.2 GPa. This reduced value of the residual pressure is then used to determine the apparent entrapment pressure (Fig. 4b). In Fig. 4b, it is shown that for the entrapment pressure within eclogite-facies conditions at 700 ^∘C, and by using a purely-elastic model, a value of entrapment pressure is inferred that is approximately 10 % less than the actual value. The amount of underestimation of entrapment pressure increases to 30 % when the entrapment temperature reaches 800 ^∘C. These values are conservative estimates since the total exhumation time is set to 1 Myr. Longer residence at high-temperature conditions would result in larger modifications of the residual pressure.

For amphibolite-facies entrapment conditions, the residual pressure that is preserved in the quartz inclusion is significantly lower compared to the case where the entrapment occurred at eclogites-facies conditions. As shown in Fig. 4d, ca. 5 % and 20 % underestimates of true entrapment pressure are predicted depending on whether the entrapment occurred at 700 ^∘C or 800 ^∘C, respectively. Similarly, the amount of underestimation will be larger if the duration of exhumation is longer than 1 Myr.

3.3 Pressure relaxation along the prograde P–T path and apparent overstepping

The pressure relaxation problem becomes more complicated when the quartz inclusion is entrapped not at the peak P–T conditions but along the prograde P–T path. In this case, viscous relaxation occurs also along the prograde P–T path and the pressure difference between host and inclusion will relax with time and increasing temperature. This effect starts before the rock reaches the peak P–T conditions. Two synthetic P–T paths are illustrated in Fig. 5. In Fig. 5a, the quartz is entrapped in the almandine-garnet host at 400 ^∘C, 1 GPa, and further experiences eclogites-facies P–T conditions. During the prograde path, the quartz inclusion will develop underpressure (e.g. Angel et al., 2015, Fig. 5), which will also be subject to viscous relaxation over geological time. The quartz pressure starts to converge towards the garnet host pressure at T > 600 ^∘C. Nearly complete viscous resetting is observed when the system is brought up to 800 ^∘C. The prograde time is set to 1 Myr or 10 Myr to compare the amount of viscous relaxation as a function of time in Fig. 5.

Figure 5 The prograde P–T path for inclusion (dashed curve) and host (solid curve). Panel (a) is for rocks that experienced eclogite-facies peak P–T conditions. The quartz inclusion is entrapped at 400 ^∘C and 1 GPa. Along the given prograde P–T path, viscous relaxation becomes significant at > 600 ^∘C. The duration of the prograde P–T path is illustrated with different colours (1 and 10 Myr; see legend). At 800 ^∘C, the quartz inclusion pressure is reset to the confining pressure (host). For rocks that experienced amphibolite-facies peak P–T conditions (b), viscous relaxation becomes significant at ca. 650–700 ^∘C and the quartz inclusion pressure is partially reset at 700 ^∘C. Pure almandine garnet is used as a host, and its flow law is from Wang and Ji (1999).

An alternative scenario is considered where the rock may also stay at the peak P–T conditions before decompression occurs. A synthetic clockwise P–T path reaching eclogite-facies metamorphic conditions is constructed as shown in Fig. 6. The quartz inclusion is entrapped in the garnet host at 400 ^∘C, 0.6 GPa, which is considered to be along the entrance of the garnet stability field. Subsequently the system is brought to 700–750 ^∘C, 1.8–1.9 GPa, conditions and stays there for 5 Myr. Afterwards, the retrograde P–T path takes 10 Myr. Two different P–T paths of quartz inclusions are constructed based on the implemented elastic and visco-elastic rheologies. Interestingly, the residual pressure of the inclusion that was subject to viscous relaxation is significantly higher (by 0.2 GPa) than the prediction of the pure elastic model as shown by the black dashed curve (0.14 GPa). The apparent entrapment pressure is calculated using the predicted residual pressure for the inclusion whose host experienced viscous relaxation. A large discrepancy exists between the apparent entrapment pressure (ca. 1 GPa at the entrapment T 400 ^∘C) and the true entrapment pressure (0.6 GPa). The overall overestimate of true entrapment pressure (0.6 GPa) is about 0.3–0.4 GPa, which may potentially be interpreted as overstepping the garnet growth or nucleation.

4 Inclusion pressure modification due to proximity to the thin-section surface

Despite the importance of viscous or plastic relaxation in the post-entrapment modification of pressure, residual pressure measurements may be different when the inclusions are closer to the thin-section surface (Enami et al., 2007). When a pressurized mineral inclusion in an infinite host under mechanical equilibrium is moved towards the thin-section surface, its pressure is released and the pressure distribution within the inclusion becomes heterogeneous. Mindlin and Cheng (1950) provided a closed-form analytical solution of stress field inside and outside a spherical inclusion with thermal strain in a semi-infinite host. The analytical solution has been generalized to ellipsoidal inclusion (Seo and Mura, 1979). Substantial mathematical investigations have also been done in deriving the analytical solution of the elastic field for inclusion in half-space (e.g. Tsuchida and Nakahara, 1970; Aderogba, 1976; Jasiuk et al., 1991). Although the analytical formulations for individual stress components of inclusion are non-trivial, here, we show that the formula for the pressure distribution of a pressurized inclusion can be significantly simplified (detailed derivations are provided in the Appendix):

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{P_{\mathrm{inc}}}{P_{\mathrm{ini}}}}=\mathrm{1}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\displaystyle \frac{\mathrm{1}+\mathit{\nu }}{{\stackrel{\mathrm{‾}}{R}}_{\mathrm{2}}^{\mathrm{3}}}}({\displaystyle \frac{\mathrm{3}{\left(\stackrel{\mathrm{‾}}{z}+\stackrel{\mathrm{‾}}{L}\right)}^{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{R}}_{\mathrm{2}}^{\mathrm{2}}}}-\mathrm{1}),\end{array}$$

where P_ini is the initial inclusion pressure in an infinite host under mechanical equilibrium, $\stackrel{\mathrm{‾}}{L}$ is the scaled inclusion depth ($\stackrel{\mathrm{‾}}{L}=L/R$) and ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{2}}=\sqrt{x_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}^{\mathrm{2}}+(x_{\mathrm{3}}+L{)}^{\mathrm{2}}}/R$ is a function of position in Cartesian coordinate system (Fig. 7), and ν is the Poisson ratio of the inclusion and host. It is emphasized that in this situation the inclusion and host possess the same elastic moduli.

The released inclusion pressure due to proximity to the thin-section surface is plotted in Fig. 7b and c using Eq. (22). Pressure release is concentrated at the top of the inclusion, while the bottom of the inclusion is subject to minimal pressure releases (< 10 %). Interestingly, the top of the inclusion is subject to negative pressure (expansion) when the inclusion is very close to the thin-section surface (see, e.g., the case of $\stackrel{\mathrm{‾}}{L}=\mathrm{1.1}$). Based on the analytical solution, the safe inclusion depth to avoid residual pressure release is ca. $\stackrel{\mathrm{‾}}{L}=\mathrm{2.5}$ (the amount of pressure release is less than 3 % within the entire inclusion). Here, the simple analytical solution in Eq. (22) can precisely model the inclusion's residual pressure due to stress release at the thin-section surface, where the inclusion possesses the same elastic moduli as the host. In a natural mineral inclusion–host system, the inclusion and host possess different elastic properties. As a case study, the stress fields of quartz-in-almandine and almandine-in-quartz systems are numerically modelled using a finite-difference (FD) thermo-elastic model (model benchmarks are provided in the Supplement). These examples are chosen to investigate two endmembers: elastically stiffer host (quartz-in-almandine in Fig. 8a) and softer host (almandine-in-quartz in Fig. 8b). Pressures at three points within the inclusion (top, centre, and bottom) are contoured as a function of $\stackrel{\mathrm{‾}}{L}$. The pressures evaluated at these three localities based on the analytical solution in Eq. (22) are also shown by the dashed curves for comparison with numerical solutions. With decreasing distance to the thin-section surface, the heterogeneity of the pressure field increases. It is shown that the pressure release is less significant in elastically stiffer host (garnet) than in elastically softer host (quartz).

Figure 6 (a) Clockwise P–T path of inclusion (dashed curve) and host (solid curve). The dashed black curve shows the inclusion P–T path based on a pure elastic model, and the blue dashed curve is based on a visco-elastic model. The quartz inclusion is entrapped in almandine garnet at 400 ^∘C, 0.6 GPa. The prograde P–T path lasts 5 Myr, and the rock stays at peak P for 5 Myr before the retrograde P–T path, which lasts 10 Myr. The residual pressure preserved by the quartz inclusion that was subject to viscous relaxation is in fact higher than the elastic limit. Therefore, the apparent entrapment pressure, calculated using elastic isomeke, is higher than the actual entrapment pressure as shown in (a). This may lead to a ca. 3–4 kbar apparent overstepping effect. The almandine flow law is from Wang and Ji (1999).

Figure 7 (a) Analytical model configuration of a mineral inclusion close to the thin-section surface. The distance between the surface to the inclusion centre is denoted by L. (b) Pressure distribution on an x–z plane. The pressure is scaled by the initial inclusion pressure (P_ini). The initial inclusion pressure is under force equilibrium in an infinite host. The analytical model describes the amount of pressure release when the inclusion approaches the thin-section surface. (c) Pressure at three localities (inclusion top, centre, and bottom) as a function of dimensionless depth L∕R. The analytical solution of Eq. (22) is used for the pressure plot.

It is shown that the difference between the analytical and the numerical solution due to the difference in elastic moduli becomes significant when the inclusion depth is shallow. The analytical and solutions are similar when evaluated at the bottom point at any depth (Fig. 8). For a quartz-in-garnet system, the analytical solution overestimates the amount of pressure release (Fig. 8a). Assuming 3 % pressure release as acceptable for the application of Raman barometry, the analytical solution yields a safe distance of ca. $\stackrel{\mathrm{‾}}{L}=\mathrm{2.0}$ for the bottom and centre point, while the numerical solution yields ca. $\stackrel{\mathrm{‾}}{L}=\mathrm{1.5}$. For the top point, the safe distance of ca. $\stackrel{\mathrm{‾}}{L}=\mathrm{2.5}$ based on the analytical solution is again higher than the prediction of ca. $\stackrel{\mathrm{‾}}{L}=\mathrm{2.0}$ based on a numerical solution. The difference in safe distance between analytical and numerical solution is due to the presence of elastically stiffer garnet host.

Differential stress (|σ₁−σ₃|) is also shown in Fig. 8b. High differential stress at the host appears when the inclusion is close to the thin-section surface. Differential stress may also exist inside the inclusion, but it is in general smaller than that of the host. For a quartz-in-garnet system, the differential stress forms a “ring”-shaped pattern with a peak at the surface. The differential stress may reach up to 3 times the expected residual pressure. This may potentially trigger plastic failure at the thin-section surface. However, for the garnet-in-quartz system, such pattern is not observed even if the inclusion depth is shallow.

Figure 8 Dimensionless pressure and differential stress plotted on an x–z plane or as a function of dimensionless depth. (a) Quartz-in-almandine system; (b) almandine-in-quartz system. For the profiles, pressure and differential stress are measured at different locations denoted by the coloured dots. In the top panel, the dashed curves in the pressure plot are based on the analytical solution in Eq. (22) considering the same elastic moduli between inclusion and host, while the solid curves are based on finite-difference results. The discrepancy between the solid (numerical solution) and dashed (analytical solution) curves in (a) is due to the fact that the host elasticity is different from the inclusion.

5 Discussion

5.1 What may cause the residual pressure modification?

The mechanisms investigated here, i.e. visco-plastic flow of the host and proximity of the inclusion to the thin-section surface, can all be responsible for the modification of the residual inclusion pressure. The amount of inclusion pressure change due to these mechanisms is controlled by the Deborah number (De), dimensionless yield strength (C^∗), and dimensionless depth ($\stackrel{\mathrm{‾}}{L}$). In the examples of quartz-in-garnet systems, the residual pressure is considered to be sealed in a perfectly elastic garnet host. Based on our study, the presence (radius) of a plastic yield region and preserved residual inclusion pressure are dominated by dimensionless yield strength $(C^{\ast }=C_{\mathrm{h}}/P_{\exp }$) as shown in Fig. 2. Strength C_h can be converted from hardness test data using the formula below (e.g. Evans and Goetze, 1979):

$$\begin{array}{}\text{(23)}& C_{\mathrm{h}}=H/C_{\mathrm{g}}\end{array}$$

where H is the measured microhardness and C_g is a constant accounting for the indenter's geometry in the experiment. Taking C_g=3 (Evans and Goetze, 1979), the yield strength of garnet host is between 4.4 and 5 GPa at room conditions (Whitney et al., 2007), which leads to $C^{\ast }=\mathrm{4.4}\sim \mathrm{5}$ given residual inclusion pressure P_exp=1 GPa. This suggests that plastic yield does not occur in an idealized scenario of isotropic, spherical quartz inclusion entrapped in an infinite garnet host. However, such an ideal scenario is highly improbable in natural samples. Localized plastic yield may still occur due to the following reasons: (1) elevated differential stress when the inclusion is close to the thin-section surface (“ring”-shaped pattern in Fig. 8); (2) stress concentration at the corners of the quartz inclusion (Whitney et al., 2000); (3) anisotropic elastic deformation of the quartz inclusion (e.g. Murri et al., 2018); (4) pre-fractures or weakness in garnet host before the entrapment of quartz inclusions that leads to the localization of dislocations. Although our model does not predict exact conditions for plastic yield due to the above possibilities, it gives a lower bound for the strength and provides information on what type of host mineral phase cannot be used for Raman barometry. Plastic yield strength of some common rock-forming minerals measured in hardness tests are compiled and provided in Table 1. As an example, given P_exp=1 GPa, the dimensionless yield strength of calcite host is ca. 0.6 and dolomite is ca. 1.5 (Wong and Bradt, 1992). This implies that plastic flow will most likely affect the residual pressure P_exp in the calcite host but not in dolomite host.

After thin-section preparation, the inclusion pressure may be (partially) released. The dimensionless depth can be evaluated by performing depth-step scan analysis with Raman spectroscopy in order to observe if the pressure gradually decreases towards the thin-section surface (Enami et al., 2007). For a quartz-in-garnet system, to avoid significant pressure release (> 3 %) in the bottom half of the inclusion, the dimensionless depth needs to be above at least 1.5 (Fig. 8). To avoid significant pressure release in the entire quartz inclusion, the dimensionless depth needs to be above ∼2. Therefore, we recommend a safe dimensionless depth of 2–2.5 (from the surface to the centre of the inclusion) for quartz-in-garnet Raman barometry (see also Mazzucchelli et al., 2018). For a 30 µm thick thin section, the maximal radius of an inclusion is ca. 6 µm (12 µm in diameter) located at the mid-point of the thin section. In practice, it is difficult to precisely measure the depth of an inclusion, and it is uncommon that an inclusion is located right in the middle of a thin section. Therefore, it is ideal to choose smaller inclusions or prepare a thicker thin section for measurement (Campomenosi et al., 2018; Mazzucchelli et al., 2018).

For commonly used quartz-in-garnet Raman barometry, our results show that below 550–600 ^∘C, the effect of viscous relaxation can be negligible. Above ca. 650–750 ^∘C, the effect of viscous relaxation needs to be taken into account depending on a particular P–T path, garnet composition, and timescale (Figs. 3, 4). This is similar to the empirical estimate ca. 750 ^∘C in Walters and Kohn (2014). It is also shown that the preserved residual pressure may even increase due to viscous relaxation if viscous resetting occurs at peak P condition (Fig. 6). This is simply because viscous creep does not only relax the overpressure in the quartz inclusion but also the underpressure that develops along the prograde P–T path. Meanwhile, the amount of viscous relaxation is time-dependent (De is a function of the operating time of viscous relaxation). Thus, the above temperature criterion for Raman barometry applies only for exhumation lasting at a million years' timescale. A higher temperature criterion for Raman barometry (e.g. ∼1000 ^∘C for garnet host at high pressure close to the coesite–quartz transition) is applicable for more rapid exhumation, e.g. xenolith ascent carried by magma (Zhong et al., 2018) or garnet synthesis experiments that lasts hours or days (Thomas and Spear, 2018; Bonazzi et al., 2019).

5.2 Implications to garnet overstepping

Quartz-in-garnet Raman barometry has been used to determine the entrapment pressure, i.e. garnet nucleation or growth conditions, and compared it to the P–T conditions determined based on phase equilibria or classical chemical thermobarometry (Castro and Spear, 2017; Spear et al., 2014). As has been shown in Fig. 6, viscous resetting occurs when the inclusion–host system is brought to high temperature (> 600–700 ^∘C). Even if the quartz inclusion is entrapped at lower P–T conditions, e.g. the garnet entrance conditions, the preserved residual inclusion pressure may still be significantly higher than predicted from the actual entrapment P–T conditions using a pure elastic model. In this case, erroneous results may emerge if one uses the relaxed residual quartz inclusion pressure to determine the entrapment pressure. In case of significant viscous resetting at peak T conditions followed by decompression, as in the case of some HP rocks, apparent garnet growth overstepping will be inferred (see Fig. 6b). In that case, care must thus be taken to interpret the discrepancy between the results of quartz-in-garnet Raman barometry and phase equilibria. As shown in the example with a synthetic clockwise P–T path (Fig. 6), ca. 3–4 kbar apparent overstepping is predicted by considering viscous resetting at peak T condition. The amount of apparent overstepping will be even larger if the exhumation process happens faster (current model assumes 10 Myr decompression time).

6 Conclusions

We first presented a 1-D visco-elasto-plastic model to study the inclusion–host system undergoing a prograde or retrograde P–T path. The non-dimensionalization of the governing equations yields two controlling parameters: the Deborah number (De) and dimensionless yield strength (C^∗), which control the amount of pressure drop due to viscous relaxation and plastic yield. Both De and C^∗ must be higher than 1 to avoid underestimating the residual pressure. Subsequently, a simplified analytical solution for inclusion pressure (Eq. 22) close to a stress-free thin-section surface is presented based on the existing analytical solution from Seo and Mura (1979). It is suggested that the distance between the thin-section surface and inclusion must be higher than 2–3 times the inclusion radius to avoid stress release.

The relevance of our presented visco-elasto-plastic model to quartz-in-garnet elastic barometry has been systematically studied. Although plastic yield is not expected for garnet host due to its high yield strength, the residual inclusion pressure preserved in quartz inclusion can be partially modified at high temperature due to time-dependent viscous creep. It is shown that above 650–700 ^∘C over a Myr timescale, viscous creep of garnet host may partially reset the quartz pressure. This may have important implications for the determination of entrapment pressure of quartz inclusion. Additionally, this may also cause apparent overstepping of garnet growth; thus care must be taken when applying quartz-in-garnet barometry at rocks which experienced high temperatures (> 600–700 ^∘C).

Appendix A

Here, we introduce a simplified formula for the pressure distribution of an initially pressurized inclusion in an infinite host moved toward a stress-free surface based on the existing analytical solution of Seo and Mura (1979). A Cartesian coordinate system is employed as shown in Fig. 7. The full stress tensor σ_ij of inclusion loaded with eigenstrains is represented as follows (Seo and Mura, 1979):

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\mathit{\sigma }}_{ij}& ={\displaystyle \frac{{\mathit{\epsilon }}^{\ast }\left(\mathrm{1}+\mathit{\nu }\right)G}{\mathrm{2}\mathit{\pi }\left(\mathrm{1}-\mathit{\nu }\right)}}\left[-\mathrm{4}\mathit{\pi }{\mathit{\delta }}_{ij}-{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial x_{\mathrm{i}}{x}_j}}+\mathrm{4}\mathit{\nu }{\mathit{\delta }}_{ij}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{3}}}}\right.\\ & +\left(\mathrm{3}-\mathrm{4}\mathit{\nu }\right)\left({\mathit{\delta }}_{\mathrm{3}i}+{\mathit{\delta }}_{\mathrm{3}j}-\mathrm{1}\right){\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{i}}{x}_j}}\\ & -\left.\left({\mathit{\delta }}_{\mathrm{3}i}+{\mathit{\delta }}_{\mathrm{3}j}\right){\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{i}}{x}_j}}-\mathrm{2}{x}_{\mathrm{3}}{\displaystyle \frac{{\partial }^{\mathrm{3}}\mathit{\varphi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{i}}{x}_j}}\right].\end{array}\end{array}$$

While for the host, stresses are given below:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{\mathit{\sigma }}_{ij}& ={\displaystyle \frac{{\mathit{\epsilon }}^{\ast }\left(\mathrm{1}+\mathit{\nu }\right)G}{\mathrm{2}\mathit{\pi }\left(\mathrm{1}-\mathit{\nu }\right)}}\left[-{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial x_{\mathrm{i}}{x}_j}}+\mathrm{4}\mathit{\nu }{\mathit{\delta }}_{ij}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{3}}}}\right.\\ & +\left(\mathrm{3}-\mathrm{4}\mathit{\nu }\right)\left({\mathit{\delta }}_{\mathrm{3}i}+{\mathit{\delta }}_{\mathrm{3}j}-\mathrm{1}\right){\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{i}}{x}_j}}\\ & \left.-\left({\mathit{\delta }}_{\mathrm{3}i}+{\mathit{\delta }}_{\mathrm{3}j}\right){\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{i}}{x}_j}}-\mathrm{2}{x}_{\mathrm{3}}{\displaystyle \frac{{\partial }^{\mathrm{3}}\mathit{\varphi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{i}}{x}_j}}\right],\end{array}\end{array}$$

where the indices of x_i ($i=\mathrm{1},\mathrm{2},\mathrm{3}$) are in Cartesian coordinate frame following the order of x, y, and z(see Fig. 7) and ε^∗ is the isotropic eigenstrain that is expressed as the difference of volumetric strain between inclusion and host assuming that they are not bounded by each other. The elliptic integrals ψ and ϕ are expressed below:

$$\begin{array}{}\text{(A3)}& \mathit{\psi }=\mathit{\pi }{R}^{\mathrm{3}}{\int }_{\mathit{\lambda }}^{\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\frac{R_{\mathrm{1}}^{\mathrm{2}}}{R^{\mathrm{2}}+s}}{{\left(R^{\mathrm{2}}+s\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}}}\mathrm{d}s,\end{array}$$

where $\mathit{\lambda }=R_{\mathrm{1}}^{\mathrm{2}}-R^{\mathrm{2}}$ for host, λ=0 for inclusion, and $R_{\mathrm{1}}=\sqrt{x_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}^{\mathrm{2}}+(x_{\mathrm{3}}-L{)}^{\mathrm{2}}}$.

$$\begin{array}{}\text{(A4)}& \mathit{\varphi }=\mathit{\pi }{R}^{\mathrm{3}}{\int }_{\mathit{\lambda }}^{\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}-\frac{R_{\mathrm{2}}^{\mathrm{2}}}{R^{\mathrm{2}}+s}}{{\left(R^{\mathrm{2}}+s\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}}}\mathrm{d}s,\end{array}$$

where $\mathit{\lambda }=R_{\mathrm{2}}^{\mathrm{2}}-R^{\mathrm{2}}$ for both host and inclusion and $R_{\mathrm{2}}=\sqrt{x_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}^{\mathrm{2}}+(x_{\mathrm{3}}+L{)}^{\mathrm{2}}}$. Here, we focus on the stress experienced by the inclusion and derive a simplified form for the pressure of inclusion. For the inclusion, the elliptic integrals are derived:

$$\begin{array}{}\text{(A5)}& {\displaystyle }\mathit{\psi }=\mathrm{2}\mathit{\pi }\left(R^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{R}_{\mathrm{1}}^{\mathrm{2}}\right),\text{(A6)}& {\displaystyle }\mathit{\varphi }={\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\mathit{\pi }{R}^{\mathrm{3}}{R}_{\mathrm{2}}^{-\mathrm{1}}.\end{array}$$

The normal stresses in the inclusion are

$$\begin{array}{}\text{(A7)}& {\displaystyle }\begin{array}{rl}{\mathit{\sigma }}_{\mathrm{11}}& ={\displaystyle \frac{{\mathit{\epsilon }}^{\ast }\left(\mathrm{1}+\mathit{\nu }\right)G}{\mathrm{2}\mathit{\pi }\left(\mathrm{1}-\mathit{\nu }\right)}}\left[-\mathrm{4}\mathit{\pi }-{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial x_{\mathrm{1}}{x}_{\mathrm{1}}}}+\mathrm{4}\mathit{\nu }{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{3}}}}\right.\\ & \left.-\left(\mathrm{3}-\mathrm{4}\mathit{\nu }\right){\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{1}}{x}_{\mathrm{1}}}}-\mathrm{2}{x}_{\mathrm{3}}{\displaystyle \frac{{\partial }^{\mathrm{3}}\mathit{\varphi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{1}}{x}_{\mathrm{1}}}}\right],\end{array}\text{(A8)}& {\displaystyle }\begin{array}{rl}{\mathit{\sigma }}_{\mathrm{22}}& ={\displaystyle \frac{{\mathit{\epsilon }}^{\ast }\left(\mathrm{1}+\mathit{\nu }\right)G}{\mathrm{2}\mathit{\pi }\left(\mathrm{1}-\mathit{\nu }\right)}}\left[-\mathrm{4}\mathit{\pi }-{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial x_{\mathrm{2}}{x}_{\mathrm{2}}}}+\mathrm{4}\mathit{\nu }{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{3}}}}\right.\\ & \left.-\left(\mathrm{3}-\mathrm{4}\mathit{\nu }\right){\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{2}}{x}_{\mathrm{2}}}}-\mathrm{2}{x}_{\mathrm{3}}{\displaystyle \frac{{\partial }^{\mathrm{3}}\mathit{\varphi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{2}}{x}_{\mathrm{2}}}}\right],\end{array}\text{(A9)}& {\displaystyle }\begin{array}{rl}{\mathit{\sigma }}_{\mathrm{33}}& ={\displaystyle \frac{{\mathit{\epsilon }}^{\ast }\left(\mathrm{1}+\mathit{\nu }\right)G}{\mathrm{2}\mathit{\pi }\left(\mathrm{1}-\mathit{\nu }\right)}}\left[-\mathrm{4}\mathit{\pi }-{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{3}}}}+\mathrm{4}\mathit{\nu }{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{3}}}}\right.\\ & \left.+\left(\mathrm{3}-\mathrm{4}\mathit{\nu }\right){\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{3}}{x}_{\mathrm{3}}}}-\mathrm{2}{x}_{\mathrm{3}}{\displaystyle \frac{{\partial }^{\mathrm{3}}\mathit{\varphi }}{\partial x_{\mathrm{3}}^{\mathrm{3}}}}-\mathrm{2}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial x_{\mathrm{3}}^{\mathrm{2}}}}\right].\end{array}\end{array}$$

By substituting ψ and ϕ into the equations above, the normal stresses can be obtained. In deriving the pressure, i.e. negative mean stress, many terms in Eqs. (A7)–(A9) can be cancelled out. A simplified form is obtained as follows:

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{P}_{\mathrm{inc}}& ={\displaystyle \frac{\mathrm{4}{\mathit{\epsilon }}^{\ast }\left(\mathrm{1}+\mathit{\nu }\right)G}{\mathrm{3}\left(\mathrm{1}-\mathit{\nu }\right)}}\\ & \left[\mathrm{1}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\displaystyle \frac{R^{\mathrm{3}}}{R_{\mathrm{2}}^{\mathrm{3}}}}\left(\mathrm{1}+\mathit{\nu }\right)\left({\displaystyle \frac{\mathrm{3}{\left(z+L\right)}^{\mathrm{2}}}{R_{\mathrm{2}}^{\mathrm{2}}}}-\mathrm{1}\right)\right].\end{array}\end{array}$$

The pre-factor $\frac{\mathrm{4}{\mathit{\epsilon }}^{\ast }\left(\mathrm{1}+\mathit{\nu }\right)G}{\mathrm{3}\left(\mathrm{1}-\mathit{\nu }\right)}$ is in fact the initial pressure of the inclusion in an infinite host loaded by the eigenstrain ε^∗ under mechanical equilibrium. Therefore, we may simplify Eq. (A10) as follows:

$$\begin{array}{}\text{(A11)}& P_{\mathrm{inc}}=P_{\mathrm{ini}}\left[\mathrm{1}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\displaystyle \frac{R^{\mathrm{3}}}{R_{\mathrm{2}}^{\mathrm{3}}}}\left(\mathrm{1}+\mathit{\nu }\right)\left({\displaystyle \frac{\mathrm{3}{\left(z+L\right)}^{\mathrm{2}}}{R_{\mathrm{2}}^{\mathrm{2}}}}-\mathrm{1}\right)\right].\end{array}$$

where P_ini is the inclusion pressure in an infinite host loaded by eigenstrain ε^∗ under mechanical equilibrium before moving it close to the thin-section surface. The equation can be non-dimensionalized by using R as a length scale shown below:

$$\begin{array}{}\text{(A12)}& {\displaystyle \frac{P_{\mathrm{inc}}}{P_{\mathrm{ini}}}}=\mathrm{1}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\displaystyle \frac{\mathrm{1}+\mathit{\nu }}{{\stackrel{\mathrm{‾}}{R}}_{\mathrm{2}}^{\mathrm{3}}}}\left({\displaystyle \frac{\mathrm{3}{\left(\stackrel{\mathrm{‾}}{z}+\stackrel{\mathrm{‾}}{L}\right)}^{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{R}}_{\mathrm{2}}^{\mathrm{2}}}}-\mathrm{1}\right).\end{array}$$

The analytical solution for pressure in the mineral inclusion subject to an initial residual pressure P_ini is obtained. When the inclusion is far from the thin-section surface ($\stackrel{\mathrm{‾}}{L}\to +\mathrm{\infty }$, and ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{2}}\to +\mathrm{\infty }$), the actual residual pressure approaches the expected residual pressure based on a classical elastic model (P_inc→P_ini). Another potential application of the solution in Eq. (A12) is for benchmarking numerical solutions. Due to the simplicity of the pressure expression, it is particularly suitable for quick validation of numerical models, e.g. the finite-difference model presented in the Supplement.