High-magnitude stresses induced by mineral-hydration reactions

Fluid-rock interactions play a critical role in Earth’s lithosphere and in engineered subsurface systems. In the absence of chemical mass transport, mineral-hydration reactions will be accompanied by a solid-volume increase that may induce differential stresses and associated reaction-induced deformation processes, such as dilatant fracturing to increase fluid permeability. However, the magnitudes of stresses that manifest in natural systems remain poorly constrained. Here we show that the simplest hydration reaction in nature MgO + H2O⇔ Mg(OH)2 can induce stresses of several hundred megapascals, with local stresses up to ∼1.5 GPa. We demonstrate that these stresses are dissipated not only by fracturing but also induce plastic deformation with dislocation densities (1015m−2) exceeding those typical of tectonically deformed rocks. If these reaction-induced stresses can be transmitted across larger length scales they may influence the bulk stress state of reacting regions. Moreover, the structural damage induced may be the first step towards catastrophic rock failure, triggering crustal seismicity.


Introduction
When fluids infiltrate rocks, resultant mineral reactions modify physicochemical rock properties, having a first-order effect on geodynamic and geochemical processes within the lithosphere. These fluid-driven mineral reactions may govern, e.g., the stability of the lithosphere (Jackson et al., 2004;Hilairet et al., 2007), the formation of ore deposits (e.g., Sillitoe, 2010) and environmental subsurface processes, such as the sequestration of anthropogenic CO 2 (e.g., Matter et al., 2016). For more than a century, scientists have been arguing whether fluid-rock interactions occur at constant solid volume or are accompanied by solid-volume changes ,∆V s ,with immediate implications for element mobility and the generation of differential stress (Lindgren, 1918;Carmichael, 1987;Fletcher and Merino, 2001;Wheeler, 2020).
While fluid-driven mineral-replacement reactions can occur without substantial positive ∆V s (Plümper et al., 2017), observations of natural rocks and laboratory experiments reveal that mineral reactions involving an increase in solid volume in a confined space can generate stresses on the confining boundaries (e.g., Jamtveit et al., 2009;Kelemen and Hirth, 2012;Plümper et al., 2012;Zheng et al., 2018). This phenomenon is known as the development of a 'force of crystallization' (Becker and Day, 1905). More recently, this process has received increased interest as it may explain the widespread hydration and carbonation of mantle peridotite (Rudge et al., 2010;Plümper et al., 2012), having implications for the deep water and carbon cycles and the sequestration of CO 2 (Kelemen and Matter, 2008). Moreover, variation in local stresses during mineral reactions may have consequences for the thermodynamics of the reaction (Wheeler, 2020).
Although efforts have been made to quantify the magnitudes of stresses that can be exerted during mineral growth, theoretical predictions of stresses exceeding 1 GPa (Kelemen and Hirth, 2012;Wolterbeek et al., 2018), experimental constraints of 20-160 MPa (Skarbek et al., 2018;Wolterbeek et al., 2018;Zheng et al., 2018) and observations of stresses 300 MPa in natural rocks (Kelemen and Hirth, 2012) vary widely. An argument that has been put forward that may explain the limitation of reaction-induced stresses is the healing of grain contacts (Wolterbeek et al., 2018) and the expulsion of the interfacial fluid film (Zheng et al., 2018;Guren et al., 2021) ceasing the reaction. Likewise, investigations of stress-relaxation mechanisms have largely focused on brittle fracturing (Jamtveit et al., 2009;Kelemen and Hirth, 2012;Plümper et al., 2012). However, mineral reactions could also induce a multitude of other elastic and inelastic deformation processes.
Here, we present an integrated observational and numerical approach to assess the magnitude of reaction-induced stresses generated during the simplest mineral-hydration reaction in nature MgO + H 2 O⇔ Mg(OH) 2 . (1) The results show that differential stresses locally reach gigapascal levels and that these stresses induce numerous elastic and inelastic deformation processes in the encompassing marble. Furthermore, pore-network modelling suggests that reaction-induced fracturing increases the rock's permeability by several orders of magnitude.

Geological setting and methods
The samples investigated here are part of the contact aureole of the Adamello massif, Italy (Müller et al., 2009, Supplemental Material, SM). The aureole developed within marbles where hydrothermal alteration caused dolomite decomposition into periclase and calcite at about 600°C. Subsequent hydration resulted in reaction (1) and an associated positive ∆V s of up to 119%. Balancing reaction (1) at constant volume requires the removal of 54% MgO (SM). In turn, Mg-transport away from the hydration sites would likely trigger the formation of new dolomite, massive brucite precipitation in cracks or considerable Mg-metasomatic effects within the vicinity of the contact aureole. None of these products were observed. The carbonate-silicate phase assemblages formed during contact metamorphism suggest a pressure of approximately 100 MPa (Müller et al., 2009). Estimates of maximum pressures do not exceed 250-350 MPa. The marbles were analyzed using multi-scale optical, electron, and X-ray microscopy (SM). Calcite domains around periclase-to-brucite hydration sites were specifically investigated using highangular resolution electron backscatter diffraction (HR-EBSD, Wallis et al., 2019). This image cross-correlation technique allows simultaneous mapping of lattice rotations and elastic-strain heterogeneities. The lattice rotations and elastic strains can be related to the densities of geometrically necessary dislocations (GNDs, ρ GND ) and residual stresses.

Microstructures
We define two microstructures within the samples; (1) a 'background microstructure (BM)' without any brucite (Fig. 1A) and (2) a 'reaction microstructure (RM)' where brucite domains are abundant on the millimeter scale ( Fig. 1B; SM). The BM is characterized by calcite grains with thin to tabular twins and an average twin density, N , of 20±15 mm -1 (Fig. 1A). Analysis of 20 thin sections shows that ∼60% of the calcite encompassing hydration sites exhibits an increase in tapered twinning with an average twin density of 295±41 mm -1 ( Fig. 1B; SM), indicative of deformation (Barber and Wenk, 1979). In ∼20% of all investigated hydration sites, optical birefringence is observed within the encompassing calcite (SM), interpreted to be the result of high local residual stresses. Also, carbonate grains encompassing hydration sites exhibit extensive fracturing normal to the brucite-carbonate interface . Comparison between backscattered electron (BSE) images (Fig. 1E) and Mg-element maps (Fig. 1F) reveals that cracks are only partially filled with brucite (arrowheads in Fig. 1E & F), suggesting that Mg-transport distances were limited. X-ray tomography  highlights that the cracks connect hydration sites spanning distances of up to hundreds of micrometers. Data from HR-EBSD and transmission electron microscopy (TEM) reveal that dislocation densities vary by orders of magnitude among different microstructural domains. Comparison between the BM and RM datasets shows an increase of two orders of magnitude in ρ GND to 10 15 ±0.4 mm -2 within the calcite enclosing the hydration sites ( Fig. 2A, B; SM). A difference in dislocation density, ρ, is also apparent from TEM observations made within the BM calcite ( Fig. 2C; ρ = 10 12 ±1 mm -2 ) and across the brucite-calcite interface ( Fig. 2D; ρ = 10 15 ±0.4 mm -2 ). Very high lattice distortions within the calcite did not allow HR-EBSD to determine ρ GND at the brucite-calcite interface (white areas; Fig. 2A $ B). However, ρ GND decays from 10 15 mm -2 to background levels of 10 12 mm -2 over distances of 50-100 µm away from the brucite-calcite interface, implying that dislocation generation is directly linked to processes that occurred at the interface.

Estimating stress magnitudes
Although HR-EBSD is able to measure residual stresses (Wallis et al., 2019), the stress fields of dislocations surrounding the hydration sites likely obscure stresses that potentially arose due to reaction (1), such that the residual stresses are dominated by stresses induced by the dislocations themselves (Wallis et al., 2021, SM). Nevertheless, the measurements of ρ GND within the calcite surrounding brucite do provide a lower bound on the differential stress, σ d , at the brucite-calcite interface via the dislocation-density piezometer (de Bresser, 1996, Fig. 2A & B second color axis) σ d = 10 −6.21(±0.86) ρ 0.62(±0.07) . (2) Correcting for background stress magnitudes (95 MPa; SM) results in a minimum σ d around the hydration sites of 120-220 MPa with local maxima of up to ∼1.5 GPa. These estimates of σ d are lower bounds because HR-EBSD does not detect either (1) orientation gradients in the direction normal to the specimen surface that may result from additional GNDs, nor (2) statistically stored dislocations that generate cancelling lattice curvature. Furthermore, if the deformation around the inclusions had not reached steady state, then ρ had not necessarily increased to the piezometric value corresponding to the applied stress. We may use the twin (Rybacki et al., 2013) and crack (Kelemen and Hirth, 2012) densities observed ( Fig. 1) to estimate the stress magnitudes via the relationships log (σ t ) = (1.29 ± 0.02) + (0.05 ± 0.05) log (N ) and where Υ Cc is the Young's modulus (70 GPa), γ Cc is the calcite surface energy (0.1-0.6 J m -1 ), and ω Cc is the crack spacing within calcite surrounding brucite. Measurements of N and ω Cc (SM) result in σ t = 280 +120 −30 MPa and σ c = 620±125 MPa (γ Cc = 0.1 J m -2 ) to σ c = 1500±300 MPa (γ Cc = 0.6 J m -2 ). The variation in σ c arises from the variability of γ Cc with varying degrees of surface coverage by aqueous fluids (Røyne et al., 2011).
The positive ∆V s during reaction (1) may not be the only cause of stress at the brucite-calcite interface. All 'host-inclusion' systems that experience pressure-temperature (P-T) changes develop pressure differences between the host and inclusion due to differences in elastic properties and thermal expansivities. The only exception is when the system follows the isomeke line, that is, a line in P-T space where the volumetric deformation of the inclusion is accommodated by the deformation of its host (Angel et al., 2015). As this path is unique, any other change in P-T will induce a pressure difference between brucite and calcite (Moulas et al., 2020). For the case of isobaric cooling during contact metamorphism the pressure of the host mineral can be approximated as constant, and the developing stresses are controlled by differential thermal expansion. Using a linear-elastic approximation (SM), we calculated the potential pressure difference due to cooling from 700°C (equilibrium periclase-to-brucite temperature at 350 MPa) to 25°C. The calculation indicates that an 'underpressure' of -300 MPa would develop. However, reaction-induced cracks emanate radially, whereas negative thermoelastic pressures within the brucite would result in concentric crack patterns (Van der Molen and Van Roermund, 1986).
Overall, we suggest that the stresses obtained via the dislocation-density piezometer, twin piezometer and crack analysis are indeed a result of the positive ∆V s during reaction (1) under confinement within the marble. To further constrain the reaction-induced stress magnitudes, we used a non-linear continuum mechanics model for viscoelastic pressure relaxation within a host-inclusion configuration (Dabrowski et al., 2015, SM). Stress magnitudes obtained via the dislocation-density piezometer (Fig. 2) show a best fit with the numerical model at 350 MPa maximum pressure (Fig. 3). The model also predicts an effective stress exponent n between 2 and 3, consistent with dislocation creep in the host matrix and, therefore, with the observed microstructures. As such, stress magnitudes imposed onto the brucite-calcite interface were at least 100 MPa, likely several hundred MPa, with local stresses up to ∼1.5 GPa. Moreover, these elevated stresses suggest that if fracturing occurred concomitantly with dislocation motion and twinning it did not fully dissipate the growing reactioninduced stresses.

Dissipation and storage of reaction-induced strain energy
Work W must be done on the environment to accommodate the periclase-to-brucite volume change in the host matrix. Based on the microstructural observations, W can be expressed as where E dislocation , E fracture , E twinning and E elastic are the energies associated with dislocation formation, reaction-induced fracturing, mechanical twinning, and residual elastic strain (SM). Assuming that the energy dissipated during dislocation formation is equal to the currently stored total strain energy around the dislocations (Anderson et al., 2017) and neglecting energy dissipation through heat and acoustic emission upon dislocation motion we estimate a mean E dislocation = 20±10 J/mol. We infer E fracture = 1.0±0.5 J/mol, E twinning = 16±2 J/mol and E elastic = 385±125 J/mol. These estimates demonstrate that dislocation formation and glide was the most important energy dissipation mechanism. However, most strain energy is stored elastically, albeit E elastic being likely an overestimation due to the contribution of local dislocation strain fields within the HR-EBSD maps. The increased dislocation density and elastic strain may influence the overall dissolution kinetics of calcite. This increase together with high, local reaction-induced stresses could promote pressure solution adding an additional dissipation mechanism. To estimate the magnitude of the potential increase in calcite dissolution rate we calculate the crystal activity of calcite ∆a Cc compared to the reference activity 1 via ∆a Cc = exp(E i /RT ) (Schott et al., 1989), where E i is either E dislocation or E elastic , R the gas constant and T is temperature (25°C). Using E dislocation results in ∆a Cc = 0.004 with an upper bound of 0.01. Hence, an increase in dislocation density has a minor effect on the dissolution rate in agreement with previous observations (Schott et al., 1989;Blum et al., 1990). On the contrary, using E elastic results in ∆a Cc = 0.4±0.1. This increase is significant and may contribute to enhancing calcite dissolution under stress. However, no direct evidence for pressure solution of the calcite grains, such as cut-off chemical zonation or indentations, was identified. This could be due to a severe reduction in interfacial fluid film thickness, hence decreased water diffusivity (Guren et al., 2021) along the calcite-brucite interface, or due to a strong chemical potential gradient towards the periclase-brucite reaction interface limiting the amount of fluid at the calcite-brucite interface. Nevertheless, complete periclase hydration is evident in all samples, implying that fluid penetrated to the periclase-brucite reaction front, at least over geological timescales.

Reaction-induced fracturing and permeability enhancement
Using pore network modelling, we estimated the fluid permeability through the reaction-induced crack network obtained via X-ray tomography ( Fig. 4; SM). Assuming that the network was transiently transmissive to aqueous fluids, we obtain a permeability of 0.5×10 -14 m 2 . This computed, core-scale permeability is five to six orders of magnitude higher than measurements on intact, corescale marbles (Fischer and Paterson, 1992) and is in agreement with the permeability-depth relationship of Ingebristen and Manning (2010) for metamorphically disturbed rocks (1×10 -13.3 to 1× -14.5 m 2 at the estimated depth of contact metamorphism).

Conclusions
Our observations of one of the simplest natural mineral-hydration reactions shows that reactioninduced stresses within rocks can locally reach GPa-level. Hence, volume-changing reactions are able to (1) exert high-magnitude differential stresses onto their surrounding matrix resulting in deviations from lithostatic pressure and (2) induce a multitude of irreversible deformation processes with, among others, implications for the transient fluid permeability of the lithosphere. If these reaction-induced stresses can be transmitted across larger length scales they may influence plate tectonic processes and the induced microscopic damage may be the first step towards catastrophic rock failure, triggering crustal seismicity.

Data availability
The data to this paper is freely available at the Utrecht University Yoda data repository under: https://public.yoda.uu.nl/geo/UU01/XVJYBS.html Wheeler, J., 2020, A unifying basis for the interplay of stress and chemical processes in the Earth: support from diverse experiments: Contributions to Mineralogy and Petrology, v. 175, no. 12, p. 1-27.

Geological setting and samples
The samples investigated here are part of the contact aureole of the Adamello massif, Italy (Fig.  S1) and were initially collected for the study by Müller et al. (2009). The aureole developed within marbles where hydrothermal alteration caused dolomite decomposition into periclase and calcite at about 600°CC. The carbonate-silicate phase assemblages formed during contact metamorphism suggest a pressure of approximately 100 MPa (Müller et al., 2009). Estimates of maximum pressures do not exceed 250-350 MPa (Pennacchioni et al., 2006).

Volume change and mass removal
Within the main text, we state that the reaction is accompanied with a positive volume increase unless considerable amounts of MgO are transported away from the site of reaction to preserve volume. Below we discuss and visualize the volume change of reaction (1) and potential constant-volume consideration in more detail.
The molar volumes at ambient conditions for MgO and Mg(OH) 2 are 11.25 and 24.63 cm 3 /mol, respectively. If the reaction occurs stoichiometrically balanced, as in reaction (1), the solid-volume change based on the given molar volumes is 119%. To balance reaction (1) at constant volume, it is convenient to rewrite the reaction to MgO + 0.46 H 2 O ⇔ 0.46 Mg(OH) 2 + 0.54 MgO (2) This notation shows that 54% MgO would need to be transported away from the site of reaction. As no Mg-metasomatism was observed at the outcrop-scale (compare to Fig. S1) nor any massive vein formation with Mg-bearing minerals the likelihood that substantial amounts of Mg were transported away from the reaction site to counterbalance the volume increase is unlikely. However, we cannot constrain whether all Mg released upon periclase hydration did take part in the immediate growth of brucite at the periclase-brucite reaction interface. Nonetheless, even if 10% of Mg was transported away from the site of reaction the positive volume change upon periclase hydration would still have been substantial. To visualize the overall volume change we have modified a classic figure published by Carmichael (1987) shown in Figure S2. 3 Sampled and analytical methods

X-ray micro-tomography
Cylinders, 8 mm in diameter, were cut out of thin-section billets from which areas of interest were determined using optical polarization microscopy. These cylinders were scanned using a TES-CAN UniTOM (TESCAN facility, Ghent, Belgium) and a Zeiss Versa Xradia 610 (Department of Earth Sciences, Utrecht University, the Netherlands). The collected three-dimensional volume was pre-processed within the Thermo Fisher Scientific Avizo 9.2 software package. A bilateral, edgepreserving filter was applied to reduce noise. Segmentation of reaction-induced cracks and brucite domains was executed using a supervised random-forest machine-learning algorithm implemented into the open-source ILASTIK software (Berg et al., 2019). Subsequently, Avizo was used to calculate the connected components and generate a skeleton of the component network. This skeleton serves as input for the pore network model (see below).

Electron backscatter diffraction
Electron backscatter diffraction (EBSD) data were acquired on a Philips XL30 field-emission gun scanning electron microscope (SEM). The SEM was equipped with an Oxford Instruments Nordlys Nano detector and Aztec acquisition software. Reference-frame conventions were validated following the approach of Britton et al. (2016). For all high-angular resolution EBSD (HR-EBSD) maps presented here, x1 is horizontal, x2 is vertical and x3 is out of the plane of the map. For HR-EBSD analysis, diffraction patterns were saved at a resolution of 1344 x 1024 pixels. Cross-correlation-based processing for HR-EBSD was performed using a restricted developmental code in MATLAB® that performs the analyses detailed by Wilkinson and coworkers (Wilkinson et al., 2006;Britton and Wilkinson, 2011;Britton and Wilkinson, 2012). This image cross-correlation technique allows simultaneous mapping of lattice rotations and elastic-strain heterogeneities with precision on the order of 0.01°and 10 −4 , respectively. The same analyses can be performed using the proprietary software CrossCourt4 distributed by BLG Vantage (http://www.hrebsd.com/wp/ crosscourt/). The cross-correlation procedure maps the shifts in 100 regions of interest within the diffraction patterns relative to their positions in a reference pattern within each grain to overdetermine the displacement-gradient tensor describing lattice rotations and elastic strains. Densities of geometrically necessary dislocations (GND) were calculated from the rotation fields following the method by Wallis et al. (2019). We fit densities of the dislocation types associated with the slip systems summarized by de Bresser and Spiers (1997) to the measured lattice curvatures but, for simplicity, we present the sum of these densities as the total GND density. Stress heterogeneities were calculated from the elastic-strain heterogeneities using Hooke's law and the elastic constants for calcite at room temperature and confining pressure of 0.1 MPa. Measured elastic strains and residual stresses are relative to the strain and stress state at the reference point in each grain. Reference points were chosen in areas of good pattern quality away from grain boundaries.

Focused ion beam scanning electron and transmission electron microscopy
Using locations identified in the HR-EBSD maps, electron-transparent thin foils were prepared for (scanning) transmission electron microscopy ((S)TEM) using a FEI Helios Nanolab G3 focused ion beam-scanning electron microscope (FIB-SEM). These foils were investigated in a FEI Talos F200X (S)TEM equipped with four energy-dispersive X-ray detectors (Super-X EDX). Bright-field (BF), dark-field (DF) and high-angle annular dark-field (HAADF) images were acquired simultaneously in STEM mode. All SEM, FIB-SEM and TEM analyses were carried out at the Microscopy Center, Utrecht University. To determine dislocation densities within the TEM images, we used the lineintercept method (Ham, 1961). Here, the dislocation density, ρ, is given as a relationship between the intersection count of a given grid and the dislocations in a TEM image with the thickness of the TEM foil via the relationship where n is the intersection count, l is the total length of the grid, and t is the thickness of the foil.

Twin-density piezometer
Although σ t agrees with values obtained through the ρ GND analysis, it is worth noting that the calcite twin piezometer has not been specifically calibrated for heterogeneous stress fields around grain boundaries (Rybacki et al., 2013). Nonetheless, the widespread increases in twin density around many brucite grains, and around many faces of those grains ( Figures S3 and S4), suggest that twin densities averaged over large numbers of grains (and hence large ranges of Schmid factors for twinning) likely provide a reasonable broad estimate of approximate average stresses within those domains.

Energy calculation
Work W must be done on the environment to accommodate the periclase-to-brucite volume change in the host matrix. Based on the microstructural observations, W can be expressed as the sum of the different energies E i related to the induced deformation processes where E dislocation , E fracture , E twinning and E elastic are the energies associated with dislocation formation, reaction-induced fracturing, mechanical twinning, and residual elastic strain. Below we outline how we calculated the different energy contributions based on the HR-EBSD results. The dislocation energy E dislocation was calculated using an approach discussed by (Blum et al., 1990), where with Here r is the distance from the dislocation, µ is the isotropic shear modulus, b the magnitude of the Burgers vector, V m the molar volume, ∆H f the enthalpy of fusion, K is a constant related to the Burgers vector orientation, where K = 1 for a screw dislocation, and K = 1 − ν for an edge dislocation with ν as the Poisson's ratio. For calcite we use µ = 32 GPa (de Bresser, 1991), the Burgers vectors b are taken from Table S1, V m = 36.9×10-6 m 3 /mol, ∆H f = 5740 J/mol and ν = 0.317.
The fracture energy E fracture was calculated via where A is the fracture area and G the energy release rate with σ the applied stress (50 MPa), a the crack length and Υ the Young's modulus (70 GPa). The twin energy E twinning is calculated via where γ e is the twin energy for e-twin {1018}⟨4041⟩ boundaries of 0.183 J/m 2 (Bruno et al., 2010) and L the length of twins identified in optical microscopy images of the samples studied. Calcite e-twins were identified to be the majority of twins around periclase-to-brucite hydration sites based on EBSD.
The elastic strain-energy density was estimated assuming that all measured strain is linearly proportional to stress via where C ijkl is the elasticity tensor of calcite after Chen et al. (2001), ϵ kl and ϵ ij are the different strain components obtained from HR-EBSD. Processed three-dimensional X-ray tomography data was imported to the PoreFlow simulation platform to compute fluid permeability via pore-network modelling (Raoof et al., 2013). Modelling parameters were chosen as follows; fluid density = 1000 kg/m 3 , fluid viscosity = 0.001 Pa·s, inlet pressure = 1 Pa, and outlet pressure = 0 Pa. The module uses an adapted version of the Hagen-Poiseuille equation, which describes the total volumetric fluid flow through the pore throat.

Calculation of elastic residual stress
For the calculation of the elastic residual pressure, we utilized the linear-elastic solution for a spherical, isotropic inclusion within a spherical host (Zhang, 1998). Furthermore, we assumed that the initial pressure of the host and of the inclusion were equal where P is the pressure, T is temperature, R is the radius of the sphere, α is the thermal expansion, and K,G are the bulk and shear moduli. The superscripts ini and fin represent the initial and final states while the subscripts h and i represent the properties for the host and the inclusion, respectively. It is implicitly assumed that host and inclusion are always in thermal equilibrium.
To obtain the values shown in the main text, we used the bulk moduli from Holland and Powell (1998), that is 48.5 GPa for brucite and 76 GPa for calcite, respectively. For the thermal expansion of brucite, we used the value 6.69×10 -5 K -1 after Fukui et al. (2003), and for the thermal expansion of calcite we used the value 4.25×10 -5 K -1 that represents thermal expansion at reference conditions (Holland and Powell, 1998). The shear modulus of calcite is taken as 32 GPa (Bass, 1995). 6.3 Estimation of stress in the host based on the model by Dabrowski et al. (2015) We consider a spherical and isotropic mineral with a power-law viscous flow law as host material. To calculate the differential stress, σ d , within this host we used the analytical solution from Dabrowski et al. (2015) where dP is the pressure difference between the inclusion and the far-field region of the host, n is the apparent stress exponent in the power-law viscous flow law and r is the radial distance from the center of the inclusion. The solution applies for values of ω d in the host, that is for r Ri > 1 (with R i being the radius of the inclusion) and for situations in which the host is thick enough so that differential stress in the host can decrease to essentially zero away from the inclusion, that is R h >> R i (with R h being the radius of the host). The strong decrease of σ d radially away from the inclusion also implies that there is no need for a strict confinement in the host to maintain the high values of σ d close to the inclusion and the pressure difference between host and inclusion. The host must simply be thick enough so that σ d decreases to essentially zero and does, hence, not cause any deformation anymore in the far-field of the host. The fact that a strict confinement is not required is important for geological applications because in natural rock units a strict confinement never exists, for example due to the Earth's free surface.
The misfit between the analytical solution and the inferred differential stresses can be quantified using a cost function of the form where the superscripts m and o represent modeled and observed values. By optimizing the function Ω we can infer the parameters dP and n. As an input to the model, we have used the differential stress depicted in the upper righthand corner of the visualization (based on HR-EBSD measurements) shown in Figure 2A of the main text as well as shown below (Fig. S9). Figure S9. Differential stress inferred from geometrically necessary dislocations obtained via HR-EBSD. Black dashed square shows the area from which the data has been taken as an input to compare with the analytical solution by Dabrowski et al. (2015). The area marked by the black dashed square is also where a FIB-SEM foil was excavated that is shown in Fig. 2D in the main text.