Analytical solution for residual stress and strain preserved in anisotropic inclusion entrapped in an isotropic host

Abstract. Raman elastic thermobarometry has recently been applied in many
petrological studies to recover the pressure and temperature (P–T) conditions of
mineral inclusion entrapment. Existing modelling methods in petrology either
adopt an assumption of a spherical, isotropic inclusion embedded in an
isotropic, infinite host or use numerical techniques such as the finite-element
method to simulate the residual stress and strain state preserved in the
non-spherical anisotropic inclusions. Here, we use the Eshelby solution to
develop an analytical framework for calculating the residual stress and
strain state of an elastically anisotropic, ellipsoidal inclusion in an
infinite, isotropic host. The analytical solution is applicable to any class
of inclusion symmetry and an arbitrary inclusion aspect ratio. Explicit
expressions are derived for some symmetry classes, including tetragonal,
hexagonal, and trigonal. The effect of changing the aspect ratio on residual stress is investigated,
including quartz, zircon, rutile, apatite, and diamond inclusions in garnet
host. Quartz is demonstrated to be the least affected, while rutile is the
most affected. For prolate quartz inclusion (c axis longer than a axis), the
effect of varying the aspect ratio on Raman shift is demonstrated to be
insignificant. When c/a=5, only ca. 0.3 cm−1 wavenumber variation is
induced as compared to the spherical inclusion shape. For oblate quartz
inclusions, the effect is more significant, when c/a=0.5, ca. 0.8 cm−1
wavenumber variation for the 464 cm−1 band is induced compared to the
reference spherical inclusion case. We also show that it is possible to fit
an effective ellipsoid to obtain a proxy for the averaged residual
stress or strain within a faceted inclusion. The difference between the
volumetrically averaged stress of a faceted inclusion and the analytically
calculated stress from the best-fitted effective ellipsoid is calculated to
obtain the root-mean-square deviation (RMSD) for quartz, zircon, rutile,
apatite, and diamond inclusions in garnet host. Based on the results of 500 randomly generated (a wide range of aspect ratio and random
crystallographic orientation) faceted inclusions, we show that the
volumetrically averaged stress serves as an excellent stress measure and the
associated RMSD is less than 2 %, except for diamond, which has a systematically
higher RMSD (ca. 8 %). This expands the applicability of the analytical
solution for any arbitrary inclusion shape in practical Raman measurements.


where e.g. 0 is the reference lattice parameter measured at room conditions and is the lattice parameter measured at 95 entrapment conditions. Note that for hexagonal and trigonal minerals (e.g. quartz), if the symmetry of lattice parameters is maintained at entrapment conditions (e.g. for quartz we keep a=b, = = 90 o , = 120 o ), Eq. 1 still holds. For lower symmetry systems with non-orthogonal crystallographic axes (triclinic and monoclinic systems), it's not possible to align all the crystallographic axes parallel to the Cartesian coordinates and the angles of , and may also change at entrapment condition compared to reference room condition. Therefore, transformation is needed to convert strains from the 100 crystallographic axes in a unit cell to the Cartesian coordinate system for modelling the mechanical interaction between the inclusion and host. This can be done by using existing software such as PASCal (Cliffe and Goodwin, 2012), Win_Strain (http://www.rossangel.com/home.htm), and STRAIN (Ohashi and Burnham, 1973). A self-written MATLAB code is provided also in the Appendix following the Ohashi's method (Ohashi and Burnham, 1973) to calculate the lattice strain based on the changes of lattice parameters. The results are the same with all existing software. For the case of an isotropic 105 host under hydrostatic entrapment stress, the initial (entrapment) strain is expected to be isotropic and the principal strain components are simply one third of the volumetric strain, which can be directly obtained from the PVT relationship.
To simulate the exhumation of the inclusion-host system to room P-T conditions, we first unload the system by applying the strain opposite to the initial host strain state, i.e. −ε i host (Fig. 1B), a procedure which leads to a stress-and strain-free host at room conditions. This is an intermediate step that ignores elastic interaction between the inclusion and host, and the 110 inclusion will possess a virtual strain ε i incl − ε i host , as the internal inclusion-host boundary experiences the unloading strain −ε i host . At this moment, the stress state of the inclusion can be readily expressed using the linear-elastic constitutive law as: ij incl (ε j incl − ε j host ), where ij incl is the elastic stiffness tensor of the inclusion at room T. Einstein summation is used. It is straightforward to note that mechanical equilibrium is not satisfied at this intermediate step because, in general, there is a stress jump between the stressed inclusion and the stress-free host. Using the proposed approach, solving the original 115 mechanical problem is reduced to superposing the homogeneous unloading strain field −ε i host with a non-uniform solution for an initially stressed and strained inclusion embedded in a stress and strain free host at room T. The latter one is practically an eigenstrain problem (Eshelby, 1957). The thermal effects on the elastic stiffness tensor ij incl have no influence on the superposed deformation field driven by the eigenstrain load due to the mismatch between the lattice strains of the inclusion and the host. The stress and strain of the inclusion that serve as driving force for elastic interaction are as follows: 120 https://doi.org/10.5194/se-2020-180 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License.
where ε i * are referred to as inclusion eigenstrains and i * are eigenstresses (Mura, 1987). The eigenstresses correspond to the stress that a soft inclusion would experience if it was perfectly confined by a stiff host, i.e. the host was not allowed to deform elastically. The eigenstrain and eigenstress are the functions of the lattice strain of inclusion and host at entrapment conditions (taking room conditions as reference state) and the stiffness tensor of the inclusion at room P-T condition.
Because mechanical equilibrium is not satisfied for the stressed inclusion embedded in a stress-free host, elastic deformation 125 will occur (stage shown in Fig. 1b to Fig. 1c). The amount of elastic deformation that affects the inclusion with a pre-strain ε i incl − ε i host in a stress-free host to mechanical equilibrium is denoted as ε i . The strain ε i is shown in Fig. 1b to 1c, and the pre-strained state with strain ε i incl − ε i host is taken as the reference state for this elastic deformation field. By adding the strain ε i to the inclusion, we obtain the final residual stress and strain state as follows: where ε i res and i res are the final residual strains and stresses of the inclusion, which are the true physical stress and strain the 130 inclusion experiences. This final stage is shown in Fig. 1c. Finding the strain ε i will solve the anisotropic inclusion problem.
This will be sought in the next section using the equivalent eigenstrain method and Eshelby's solution.

Solving the problem with Eshelby's solution
The Eshelby's solution treats a homogeneous, ellipsoidal, isotropic inclusion embedded in an infinite isotropic host (Eshelby, 1957). Following Eshelby (1957), we replace the inclusion-host system by an isotropic homogeneous space with elastic 135 tensor ij host and load the ellipsoidal inclusion region with an equivalent eigenstrain i * (to differentiate from the previously introduced eigenstrain term ε i * due to lattice strain mismatch). Without elastic interaction, the inclusion would experience a stress − ij host j * under perfect confinement (e.g. a positive eigenstrain (expansion) leads to compressive stress, which is negative). After elastic interaction, the inclusion strain and stress under mechanical equilibrium due to a constant eigenstrain (eigenstress) load applied to an ellipsoidal region of otherwise homogeneous elastic space, can be expressed as follows 140 (Eshelby, 1957;Mura, 1987 where S ij is the Eshelby's tensor, which gives the one-to-one mapping between the equivalent eigenstrain ( j * ) and a homogeneous residual strain ( i res ) within the inclusion region. The Eshelby's solution is manifested in this tensor, which is only a function of the inclusion shape and the Poisson ratio ( ) of the isotropic host. For a spherical inclusion, S ij is symmetric and can be significantly simplified as follows: 145 All the other components are zero. For general ellipsoidal inclusions, the S ij tensor is given in the Appendix and a MATLAB script for calculating the S ij tensor is provided in the supplementary materials (see Appendix for more details for using the script).
Following the equivalent eigenstrain method (Mura, 1987, chapter 4), one may appropriately choose the equivalent eigenstrain j * to let S ij j * in Eq. 4 equals the strain ε i in Eq. 3 that drives the pre-strained anisotropic inclusion into 150 mechanical equilibrium with stress-free isotropic host, i.e. we have: By doing so, the stresses in the original anisotropic heterogeneity and the equivalent isotropic inclusion will be equal. This is because the host is stressed (and strained) by the same amount following Eq. 6, which leads to the same inclusion stress because the traction is matched between inclusion and host. By replacing the strain ε i = S ij j * into i res in Eq. 3 and equating the stresses i res in Eq. 3 with Eq. 4, we obtain the following relation: 155 The equivalent eigenstrain k * can be solved from this equation. By substituting the obtained k * back into Eq. 7, we may concisely express the final result for the residual strain and stress of the anisotropic inclusion embedded in isotropic infinite host: where I ij is the identity matrix. The dimensionless matrix M kj can be expressed as follows: https://doi.org/10.5194/se-2020-180 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License.
The dimensionless matrix M ij depends on the elastic stiffness properties of the inclusion and the host as well as the aspect 160 ratio of the inclusion manifested by the Eshelby's tensor. The components of this matrix are close to zero for a stiff host or a soft inclusion (no elastic relaxation so that i res → ij incl ε j * ) and it approaches the identity matrix for an infinitely soft host (in this case, i res → 0). An extreme case is represented by gas/liquid inclusion whose ij incl is low compared to the host stiffness, thus this dimensionless matrix M ij approaches zero and the isochoric assumption for the gas/liquid inclusion is justified. The M ij matrix can be readily calculated by using the published elastic stiffness tensor at room P-T conditions (e.g. 165 Bass, 1995). A MATLAB script is given in the supplementary data to perform this task (details of using the code can be found in Appendix)

Back-calculating eigenstrain terms based on residual inclusion strain
The wavenumber shifts of Raman peaks are induced by lattice strain. By measuring wavenumber shift of the inclusion in a thin-section, it is possible to recover the residual strain preserved within the inclusion (Angel et al., 2019;Murri et al., 2018). 170 This can be done by using the Grüneisen tensor. Therefore, ε i res can be obtained with e.g. least-square fitting method (Murri et al., 2018) and the residual stress can be readily expressed as i res = ij incl ε j res .
By inverting the right-hand matrix in Eq. 8, the eigenstrain terms can be expressed as a function of residual strain ε i res : For tetragonal or hexagonal minerals, e.g. zircon, rutile and apatite, the stiffness tensor comprises six independent components: 11 incl , 12 incl , 13 incl , 33 incl , 44 incl , 66 incl . For trigonal symmetry such as in the case of -quartz, another independent 175 component 14 incl = − 24 incl is also present. For all these mineral inclusions, and for axially symmetric residual strains, ε i * can be significantly simplified as follows (ε x * = ε y * ≠ ε * ): where ̅ ij incl = ij incl / is the dimensionless inclusion stiffness tensor scaled by the host shear modulus. Interestingly, for trigonal -quartz inclusions ( 14 incl ≠ 0), the stiffness tensor component 14 incl is not present in the expression given by Eq. 11.
In fact, the terms in the brackets are simply the residual stress components: By substituting the stiffness tensor components and the measured residual strains, the eigenstrains can be directly calculated.
The equation above thus allows estimation of the entrapment (hydrostatic or non-hydrostatic) stress (or strain) conditions by known the residual stress and strain conditions of the inclusion.

Cross validation against finite element solution
We have validated our implementation of the proposed analytical framework against independent finite element (FE) 185 solutions. A self-written 3D FE code is used to validate the presented analytical solution (Dabrowski et al., 2008;Zhong et al., 2018). For validation purposes, we used spheroidal quartz inclusions in an almandine garnet host. Adaptive tetrahedral computational meshes, with the highest resolution within and around the inclusion, are generated with Tetgen software (Si, 2015). The anisotropic elastic properties of quartz inclusion at room T are based on Heyliger et al. (2003). The host garnet elasticity is first treated as isotropic based on Jiang et al. (2004). The model length is set as 10 times (denoted as *10 below) 190 the inclusion's diameter (for spheroidal inclusion case, the model domain is a box where the side lengths are proportional to the corresponding axes lengths of the inclusion). For the model validation, we adjust the eigenstrain term to generate precisely 1 GPa compressive residual stress for spherical inclusion in infinite garnet host. This is done by letting xx res = yy res = zz res = −1 GPa hydrostatically and substituting i res into Eq. 8 to back-calculate the eigenstrain ε j * . Then, the spheroidal inclusion is loaded by the calculated eigenstrain ε j * in the FE code, and the residual stress is compared to results 195 obtained by the presented analytical solution. The choice of eigenstrains (either loading the inclusion to 1 GPa pressure or any other residual stress value) is not influential for the validation purposes, as long as both the analytical and FE methods take the same eigenstrains. This and other successful, more general tests with arbitrary aspect ratio and eigenstrains, have been performed but are not reported here.
In Fig. 2a, the numerically and analytically obtained residual stresses are plotted together as a function of the aspect ratio of 200 the tested spheroidal inclusions. In Fig. 2b, the difference is plotted as a function of element count and boundary distance (*5, *10 and *20). It is clearly shown that the two sets of solutions converge with increasing the number of mesh elements and the computational box size. The success of this convergence test validates the correctness of our presented analytical model (also FE code) for an anisotropic ellipsoidal inclusion entrapped in an isotropic space.
In addition, we have also tested the effect of applying cubic elastic stiffness tensor of almandine from Jiang et al. (2004) and 205 compared the residual stress with the results obtained for an elastically isotropic garnet (blue dots in Fig. 2a). The difference in residual stresses obtained with FE method using anisotropic garnet host and the analytical solution (implicitly assuming isotropic host) is less than 0.001 GPa. This suggests that it is not necessary in practice to consider the anisotropy of garnet https://doi.org/10.5194/se-2020-180 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License.
host. This has been also reported in e.g. Mazzucchelli et al. (2019). It was suggested that the elastic anisotropy of cubic garnet has no significant impact on the result of elastic barometry. Thus, effective isotropic elastic properties of garnets may 210 be used to model the inclusion-host elastic interaction.

Effect of inclusion aspect ratio on residual stresses
In Eq. 8, the aspect ratio of the inclusion only affects the Eshelby tensor. Here, we choose some common inclusions in metamorphic garnets, including quartz, apatite, zircon and rutile as examples to test the effect of inclusion aspect ratio on 215 residual stresses. The data sources of elastic stiffness tensors of the studied minerals are listed in the caption of Fig. 3. Here, we first focus on spheroidal inclusions, and let the crystallographic c-axis coincide with the geometrical z-axis of the inclusion. The inclusion aspect ratio is controlled by varying the lengths of the principal z-and x-(y)-axes, and it is denoted by c/a ratio for simplicity in the following text (note it is not the ratio of the lattice parameters c and a).
To isolate the effect of the inclusion aspect ratio, the eigenstrains for various inclusion minerals are all set to create 1 GPa 220 compressive hydrostatic residual stress for the reference spherical inclusion embedded in an isotropic, infinite garnet host, which is the same approach as in the previous "cross-validation" section. Therefore, the stress variations shown in Fig. 3 only represent the mechanical effect purely due to varying the geometrical aspect ratio c/a, and they allow direct comparison among different common mineral inclusions in metamorphic garnets.
Among all the tested minerals, the residual stress in quartz inclusions is the least sensitive to variations in aspect ratio. For 225 prolate quartz inclusion (c/a>1), the change of x res due to shape variation is within 0.1 GPa and the change of z res is within 0.2 GPa. The effects of varying the aspect ratio of prolate zircon and apatite inclusions are similar, with variation of x res from the reference spherical case reaching up to ca. 0.12 GPa and res up to ca. 0.2 GPa.
The residual stress in rutile is the most sensitive to aspect ratio ariations. With increasing c/a ratio from 1, x res in prolate rutile inclusions increases up to ca. 0.2 GPa and res decreases by ca. 0.6 GPa. This is relevant for practical measurement as 230 rutile often forms needle-shaped crystals.
The pressure (negative mean stress) is significantly less sensitive to inclusion aspect ratio variations. For prolate inclusions of quartz, zircon and apatite, the residual pressure differs from the reference level (spherical inclusion shape) by only ca. 0.01 GPa when the aspect ratio c/a is stretched up to 4. For oblate inclusion (c/a<1), the pressure variation is typically below 0.1 GPa. 235 The residual stress in mineral inclusions can be easily converted into residual strain, which can be directly translated into Raman shifts (Angel et al., 2019;Murri et al., 2018). The effects of varying the aspect ratio of a quartz inclusion on Raman wavenumber shifts (see Fig. 4) are determined using the calculated residual strain components and the Grüneisen tensor (Murri et al., 2018). The same model settings are applied, with 1 GPa compressive hydrostatic residual stress characterizing https://doi.org/10.5194/se-2020-180 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License. the case of a spherical quartz inclusion. It is shown that for prolate inclusions, aspect ratio only introduces minor effects on 240 the Raman shifts. For example, varying the c/a ratio between 1 and 5 induces a wavenumber variation of less than 0.3 cm -1 for the 464 cm -1 Raman peak. This is in most cases insignificant from the viewpoint of practical Raman measurements. The b/a ratio variations are also shown to be insignificant for the spectral shifts, with changes less than 0.2 cm -1 for the 464 cm -1 peak. For oblate inclusion, the impact of inclusion shape is shown to be more significant. For the 464 cm -1 peak, the change of wavenumber shift can reach 0.8 cm -1 for strongly oblate inclusion (c/a=0.25), and it is ca. 0.3-0.4 cm -1 for less oblate 245 inclusion (c/a=0.5).
Our results show good consistency with the Raman data reported in Kouketsu et al. (2014), where quartz inclusions with different aspect ratio were measured and no significant variation on spectral shift was found.

Effect of inclusion crystallographic orientation with respect to long-axis
In nature, the crystallographic axes of an inclusion are not necessarily aligned parallel to its geometrical axes. In this section, 250 the effect of varying the crystallographic orientation with respect to the geometrical axes on the resulting Raman wavenumber shift is systematically studied using the proposed analytical model. Here, we reorient the crystallographic caxis of a spheroidal quartz inclusion from its long axis (Fig. 5) between 0 (crystallographic c-axis parallel to the geometrical long axis) and 90 degrees (c-axis perpendicular to the long axis). The same eigenstrain is applied for quartz as from previous section. The elastic stiffness tensor of quartz is from Heyliger et al., (2003) and of isotropic almandine garnet is from Milani 255 et al., (2015). The predicted Raman spectral shifts are calculated based on the Grüneisen tensor calibrated by Murri et al. (2018).
The results are shown in Fig. 5. For an aspect ratio of 2, the 464 cm -1 band varies ca. ±0.2 cm -1 when the orientation of the crystallographic c-axis is varied between the long and the short geometrical axis. The effect of crystallographic orientation (c-axis) on the Raman shifts increases towards higher aspect ratio. For an aspect ratio of 5, the 464 cm -1 band varies ca. ±0.4 260 cm -1 . Similarly, for 206 and 128 cm -1 bands, the maximal wavenumber variations compared to the reference case of a spherical inclusion are ca. 0.8 and 0.3 cm -1 , respectively. The results suggest that the crystallographic orientation of a quartz inclusion with respect to its geometrical axes has no significant impact on the predicted Raman spectral shift, as long as the geometrical aspect ratio is not higher than 2-3.

Effect of faceted inclusion shape 265
The analytical solution presented in this study is derived for an ellipsoidal anisotropic inclusion. However, the shape of a natural mineral inclusion may exhibit corners, edges and facets, which results in stress concentration effects and may have an impact on the overall level of the residual stress.
Here, we explore the possibility of using an effective ellipsoid to fit the shape of a faceted inclusion. We use the equivalent aspect ratio to calculate the residual stress/strain based on the presented analytical solution for ellipsoidal inclusions. Fitting an arbitrary, irregular shape using an ellipsoid in 3D (or an ellipse in 2D) is a common practice in image analysis (Chaudhuri and Samanta, 1991;Li et al., 1999). A pixelated 3D image is used to calculate the second-order moment of the object shape to minimize the mismatch between the 3D irregular inclusion shape and the effective ellipsoid. The method allows for obtaining the lengths and orientations of the major, minor and intermediate axes of the effective ellipsoid (method described in Appendix and a MATLAB code is provided in supplementary data to perform this task). 275 Similar to previous sections, we use the eigenstrain components that can load the reference spherical inclusion of any given mineral embedded in isotropic almandine garnet host into 1 GPa compressive residual stress. The tested inclusion shapes include: cylinder, tetrahedron, cuboid, octahedron, hexagonal prism, and icosahedron. To vary the aspect ratio, the inclusion shape is stretched in the z-axis direction, which is parallel to the crystallographic c-axis as shown in Fig. 6.
We study five inclusion minerals: quartz (elastic tensor from Heyliger et al., 2003), zircon (Bass, 1995), rutile (Wachtman et 280 al., 1962), fluorapatite (Sha et al., 1994), and diamond (Bass, 1995). Almandine garnet is taken as the host grain (Milani et al., 2015). For each FE mesh, the size of the computational box is set more than 10 times the inclusion size and adaptive mesh is generated with highest mesh resolution within and close to the inclusion. 10-node tetrahedron elements with quadratic (second-order) shape (interpolation) functions for the displacement field are used. In total, there are ca. 2 million elements per model (numerical error less than 0.0003 based on benchmark results in Fig. 2). 285 The results of numerical simulations are shown in Fig. 6. The effective aspect ratio for all different inclusion shapes together with the residual stress components are given in supplementary data. The residual stress in non-ellipsoidal inclusions based on FE model is heterogeneous and we monitor the stress state: 1) at the centroid (shorten as CT) point (red dots in Fig. 6), 2) as the volumetric average (VA) within the entire inclusion (orange dots) The root mean square deviation (RMSD) is calculated by comparing the residual stress from the finite element solutions 290 based on various stress evaluation scheme (CT/VA) and analytical model using the best-fitted effective aspect ratio (Table   1). It is clearly shown that the VA stresses of quartz, zircon, rutile and apatite inclusion are remarkably similar to the analytical results, with RMSD generally lower than 0.02-0.03 GPa (ca. 2%). From CT to VA, a significant improvement on RMSD of a factor of 2 to 3 is obtained.
The only exception among the studied minerals is diamond, where the RMSD is higher than in the case of other inclusions, 295 which are elastically softer. This is consistent with the high "geometrical correction factor" reported for diamond in Mazzucchelli et al. (2018). However, as an improvement from Mazzucchelli et al. (2018), where the geometrical correction factor must be applied for all inclusion phases to correct the residual stress due to shape effects, we have found that the stress variation due to varying inclusion shape for minerals such as quartz, zircon, apatite and rutile can be satisfactorily approximated by using the proposed approach of the equivalent ellipsoidal inclusion, with RMSD generally lower than 3-4% 300 for most of the studied inclusion shapes. To achieve this improved and satisfactory level of prediction: 1) we have used bestfit ellipsoids to better approximate inclusion shapes, instead of a crude measure of the length/width ratio of e.g. cuboidal or https://doi.org/10.5194/se-2020-180 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License. cylindrical inclusion; 2) we have considered not only the centroid point of the inclusion (which indeed yields a larger RMSD), but also the volumetric average (VA) for the residual stress state sampled during stress measurements, which interestingly provides a significantly better approximation for the residual stress/strain state of the tested mineral inclusions. 305 This is practical and useful in Raman measurement because it is possible to perform either 1) multiple-point averaging during Raman analysis within the entire inclusion, or 2) defocus the laser beam to take into account a larger volume for the inclusion strain heterogeneity. For tiny inclusions (size of ca. 1~3 μm) and for a typical in-plane laser beam diameter of ca.
1~2 μm, the stress/strain averaging is, in fact, implicitly performed during measurements,. Based on FE analysis, it is clearly shown that the volumetrically averaged stress within the inclusion, rather than centroid point measurements, may provide a 310 closer match compared to the stress predicted based on the presented analytical solution developed for the best-fitted ellipsoid. This effect becomes statistically more significant when the faceted inclusion shape and crystallographic orientation are independent as demonstrated in the next section.

Irregularly faceted shapes and random crystallographic orientation
In nature, the shape of mineral inclusions is not necessarily highly symmetric as treated in the previous section and the 315 crystallographic orientation can be generally random with respect to the principal geometric axes of the inclusion. Here, a MATLAB script is used to generate completely random 3D inclusion shapes by prescribing random vertices (non-coplanar 5 to 24 vertices) and connecting them to form a closed 3D shape. Delaunay triangulation is used to form 3D volumes enclosed by the triangular faces. "Tetgen" software is again used to generate unstructured tetrahedron computational meshes fitted to the inclusion surface. The effective aspect ratio (geometrical longest to shortest axis of the best-fitted effective ellipsoid) is 320 controlled to be within 6. In total, we have generated ca. 500 random 3D inclusion shapes and performed finite element simulation for the previously studied set of anisotropic inclusion minerals (quartz, zircon, rutile, apatite and diamond) to calculate the elastic stress field. The generated random shapes are plotted in supplementary data (see the .gif animation to illustrate 100 selected examples of 3D inclusion shapes). We further allow the crystallographic c-axis to be pointing along the orientation randomly chosen parallel to either the longest, intermediate or shortest geometrical axis (best-fitted using the 325 method of Chaudhuri and Samanta, 1991). The FE results are compared to the analytical results based on the effective ellipsoidal inclusion with the same crystallographic orientation. This Monte-Carlo type FE simulation allows us to investigate how much stress deviation can be generated for irregularly faceted inclusion shapes with random crystallographic orientation, and how accurate the analytical approach based on the best fitted ellipsoid is to describe the residual stress state in an irregularly shaped inclusion, depending on the stress sampling scheme (CT/VA). The results are plotted in Fig. 7 (raw 330 data of FE simulations can be found in supplementary data).
For centroid point (CT), quartz inclusions have the lowest RMSD of ca. 0.03 GPa for all the three normal stress components and diamond inclusions have the highest RMSD of ca. 0.11 GPa. In general, CT stress shows a systematically higher RMSD than VA stresses. When the stress is volumetrically averaged within the inclusion (VA), the RMSD dramatically drops to a nearly perfect match between the FE results for irregularly faceted inclusion and the analytical prediction based on the best fitted ellipsoid. The RMSD of volumetrically averaged residual stresses (VA) of quartz, zircon, rutile and apatite are all lower than ca. 0.02 GPa (2%) and it shows no obvious dependence on the effective aspect ratio even for the extremely elongated or flattened inclusions (see the near-perfect alignment of the orange dots and 1-to-1 ratio line in the middle and bottom panel of Fig. 7).
Thus the volumetric average of the residual stress within the inclusion provides a sufficiently reliable match between the 340 exact results for irregularly shaped inclusions and the approximate predictions based on the analytical solution. This shows that it is possible to approximate the stress/strain state of the inclusion using an effective ellipsoid shape for the tested inclusions including quartz, zircon, rutile and apatite. Only diamond has a notably higher level of RMSD, exhibiting a systematic discrepancy between the exact numerical and approximate analytical results. This indicates that using the proposed equivalent analytical model for diamond inclusion may lead to a potential overestimation of the residual stress by 345 ca. 0.07 GPa (7%). However, this RMSD may still be acceptable as a crude estimate or may serve as an upper limit for elastic thermobarometry.

Non-linear elasticity at room T
The presented model builds on a linear elastic constitutive law at room temperature, i.e. i res = ij incl j res . This assumption is appropriate when the residual stresses/strains of the inclusion are low, thus the application of a constant anisotropic stiffness 350 tensor ij incl determined at room P-T conditions introduces no significant errors. For highly stressed mineral inclusions, e.g.
inclusions in diamond host from mantle xenoliths where the residual inclusion pressure may reach several GPa, this approximation may lead to non-negligible deviations. To eliminate such error, the stiffness tensor ij incl needs to be treated as a function of either non-hydrostatic stresses or anisotropic strains, i.e. ij incl ( i res ) or ij incl ( i res ). In experimental studies, ij incl is often described as a function of hydrostatic pressure, e.g. Bass (1995). It is beyond the scope of this paper to develop 355 a method of fitting ij incl with respect to the individual stress tensor components. If the stiffness tensor ij incl can be parameterized by i res or simply as a function of pressure as a first-order approximation, the residual stresses/strains are readily calculated by iterating Eq. 8, while updating the ij incl tensor using the calculated inclusion stress or strain during the iteration loop. Thus, the developed analytical method based on the Eshelby's solution can be extended to the case of a nonlinear inclusion phase as long as ij incl can be parametrized in terms of stress or strain components, or their invariants. 360

Concluding remarks and petrological implications
In this study, we use the classical Eshelby solution combined with the equivalent eigenstrain method to calculate the residual strain and stress in an anisotropic, ellipsoidal mineral inclusion embedded in an elastically isotropic host. The residual stresses can be expressed by a linear operator (Eq. 8) acting on the eigenstrain. The linear operator depends on the anisotropic elastic stiffness tensor of the inclusion evaluated at room P-T conditions, the shape of the inclusion, and the shear modulus and Poisson ratio of the host. The studied mechanical problem is loaded by an eigenstrain term, which is given by the difference between the lattice strains of the inclusion and host at the P-T conditions of entrapment.
The effect of inclusion aspect ratio on the inclusion residual stress and strain has been investigated quantitatively. The residual stress in quartz inclusions exhibits the least sensitivity to aspect ratio changes and rutile shows the most pronounced variation. The popularly used quartz-in-garnet system is studied in more details. For prolate quartz inclusions, the residual 370 stress variations caused by varying inclusion shape are shown to be insignificant when the crystallographic c-axis is subparallel to the geometrical long axis. The Raman wavenumber variation is less than 0.4 cm -1 for the 464 cm -1 peak even for highly elongated inclusions with an aspect ratio of 5. For oblate quartz inclusions with an aspect ratio of ca. 0.5, the additional wavenumber shift may reach ca. 0.8 cm -1 . Therefore, it is useful in practice, although potentially technically difficult, to have an estimate of the crystallographic c-axis orientation when studying highly stretched or flattened quartz 375 inclusions. As long as the c-axis is sub-parallel to the geometrical long-axis, the additional wavenumber shifts due to the inclusion aspect ratio is minor.
Our proposed analytical procedures to model residual inclusion stress and strain state do not require pre-FE simulation to obtain the 6-by-6 "relaxation tensor" as proposed by Mazzucchelli et al. (2019). For application purposes, as long as the lattice strains of inclusion (ε i incl ) and host (ε i host ) at high P-T conditions are available, it is possible to calculate the 380 eigenstrain term by subtracting them following Eq. 2. Given the driving eigenstrain, the residual strain and stress preserved within an anisotropic, ellipsoidal inclusion in isotropic host can be easily calculated using Eq. 8. The proposed procedure can be inversely applied to retrieve the residual strain/stress of any natural mineral inclusions embedded in elastically isotropic hosts, such as garnets.
The presented model is only exact for perfectly ellipsoidal inclusions. In nature, inclusions often possess different shapes 385 with facets and edges. Finite element simulations on various faceted inclusion shapes showed that the residual stress is modified to a different degree as compared to the simple ellipsoidal inclusion case, depending on the relative elastic properties between the inclusion and the host grain. However, the proposed approach of using the analytical result for the best-fitted effective ellipsoids yields remarkably good approximation for all the tested inclusion shapes. including highly irregular. The RMSD comparing the FE numerical solution for faceted inclusion and the analytical solution based on 390 effective best-fitted ellipsoid is typically less than 2% for quartz, zircon, apatite and rutile inclusions. The only exception are the elastically stiff diamond inclusions, where the RMSD reaches 7%. This finding expands the applicability of the analytical framework to arbitrarily shaped inclusions, whose elastic stiffness is not signifficantly higher than host (such as quartz, rutile, zircon and apatite). One important petrological implication is that it is possible to take the volumetrically averaged stress/strain within the inclusion and use it as a proxy to represent the residual stress/strain state of the inclusion. Then the 395 proposed analytical framework may be used to recover the entrapment condition by back-calculating the eigenstrain using the volumetrically averaged residual stress/strain and the effective ellipsoid aspect ratio of the inclusion (Eq. 8). In fact, https://doi.org/10.5194/se-2020-180 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License. Tables   Table 1. Root mean square deviation (RMSD) of finite element solution of symmetrically shaped non-ellipsoidal inclusion in Fig. 6 compared to the analytical solution of equivalent spheroidal inclusions. Isotropic almandine garnet is used as host.

Figures and
For each inclusion mineral and inclusion shape, the aspect ratio varies from 0.2~5. Effective aspect ratio is calculated for each shape and used for the analytical solution to obtain the residual stress state. The inclusion is loaded by eigenstrain that 525 creates 1 GPa hydrostatic residual pressure for spherical inclusion in infinite host. Thus, any stress variation can only be caused by shape change. The calculated stress data for each individual FE run is given in supplementary data. Stress is obtained for 1) the centroid point (CT), and 2) volumetric average (VA) of the entire inclusion (see Fig. 6 for illustration). Inclusion-host at entrapment conditions. The stress is homogeneous as σ trap but strains are different as ε i incl and ε i host . (b) First, relax the inclusion and host by −ε i host to room P-T conditions. Without elastic interaction, the inclusion has strain ε i incl − ε i host and stress ij incl (ε j incl − ε j host ) so the system is not in mechanical equilibrium. (c) Elastic interaction occurs to reach equilibrium by adding strain ε i to the inclusion (host also deforms). The residual inclusion stress is ij incl (ε j incl − ε j host + ε j ). (d) Equivalent scenario where the inclusion and host are initially stress free at room P-T and they both have 540 isotropic elasticity of ij host . (e) Equivalent eigenstrains i * are loaded to the inclusion. Without elastic interaction, the inclusion has stress − ij host j * . Eshelby's method is applied to obtain the final strain state in isotropic inclusion as S ij j * and stress as ij host (S jk k * − j * ), where S ij is the Eshelby's tensor (Eshelby, 1957). Equivalent eigenstrain method states that by properly choosing i * , the relation ε i = S ij j * can be satisfied (Mura, 1987, chapter 4). The stress of isotropic inclusion (f) as ij host (S jk k * − j * ) equals the stress of the anisotropic inclusion (c) as ij incl (ε j incl − ε j host + S jk k * ) (see Eq. 7). By solving for 545 j * , we obtain the residual stress and strain of anisotropic inclusion in isotropic host in Eq. 8.  inclusion (c/a=1, b/a=1). The initial residual inclusion pressure is assumed to be hydrostatic 1 GPa for the reference spherical inclusion. The wavenumber shift variation are due to changing the geometrical aspect ratios of c/a and b/a axes.
The Raman shifts are calculated using the residual strain and the Grüneisen tensor in Murri et al. (2018). For c/a>1, the 570 inclusion is prolate and for c/a<1, the inclusion is oblate. The stiffness tensor of quartz at room P-T is from Heyliger et al., inclusion c/a=1 (in this case the crystallographic orientation does not matter). The horizontal axis represents the aspect ratio of the spheroidal inclusion, and the vertical axis shows the angle between the crystallographic c-axis and the geometrical long axis. In the plot, c-axis is allowed to shift from parallel to the geometrical long axis to parallel to geometrical short axis of the spheroidal inclusion. The driving eigenstrain is set to produce a hydrostatic residual pressure of 1GPa in the reference 580 spherical inclusion.
https://doi.org/10.5194/se-2020-180 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License. Fig. 6. Finite element stress of various inclusion shapes (symbols) compared to the stress of effective spheroidal inclusion 585 based on analytical method (black curves). The effective aspect ratio of inclusion shape is calculated based on the fitting method of Chaudhuri and Samanta, (1991) and Li et al. (1999) (see Appendix). The inclusion is loaded with eigenstrain that generates 1 GPa compressive hydrostatic residual stress for spherical shape. The variation of stress is only caused by the shape change. The c-axis coincides with the streching direction. The red dots correspond to the stress at the inclusion's centroid (CT), the orange dots correspond to the volumetric average (VA) of the entire inclusion. The anisotropic elastic 590 stiffness tensor are listed in the caption of Fig. 3. The root mean square deviation (RMSD) for each inclusion shape and inclusion mineral phase is given in Table 1. The raw FE stress data can be found in supplementary excel file.
To obtain the similar transformation matrix relating the deformed crystallographic axes at entrapment conditions to Cartesian coordinates can be easily done by replacing 0 , 0 0 , 0 , 0 , 0 measured at reference P-T state to , , , , , that are measured at entrapment condition from Eq. A1. This transformation matrix is denoted as 1 . We then calculate the displacement gradient tensor : where I is the identity matrix. Without considering the antisymmetric rotation tensor, the infinitesimal Lagrangian strain 620 tensor can be expressed as follows: A MATLAB code is provided to perform this calculation (Calculate_Strain.m). The input values are the reference lattice parameters measured at room P-T conditions and the deformed lattice parameters at the entrapment conditions. The outputs are both the infinitesimal and finite Lagrangian strain tensor reported in MATLAB commend window. The results are https://doi.org/10.5194/se-2020-180 Preprint. Discussion started: 26 October 2020 c Author(s) 2020. CC BY 4.0 License.
where 1 , 2 and 3 are the lengths of three principal axes of the inclusions and they follow the order of 1 > 2 > 3 . In case this order needs to be changes, a simple 90 degree rotation can be executed on S ij . The provided code in supplementary data automatically perform such rotation to adjust the axes into correct order. The required tensors and are evaluated as 635 follows: The integrals ( , ) and ( , ) are evaluated using the method of the arithmetic-geometric mean (Abramowitz and Stegenm, 1964, chapter 17). Once ( , ) and ( , ) are obtained, can be computed and substituted into the Eshelby's tensor.

Fit arbitrary inclusion shape with effective ellipsoid
The method (details see Chaudhuri and Samanta, 1991;Li et al., 1999) requires a 3D data/image of the inclusion consist of regular cubic voxels, which has volume Δ in each voxel and coordinate x, y, z at the center of each voxel. The inclusion is denoted as domain R. Its second-order moment matrix is calculated as follows: 645 = ∭ ( 2 + 2 ) ≈ ∑ (