Rotation of the Stress Tensor in a Westerly Granite Sample During the Triaxial Compression Test

We simulated the spatiotemporal modelling of 3D stress and strain distributions during the triaxial compression laboratory test on a westerly granite sample using finite-difference numerical modelling implemented with FLAC3D software. The modelling was performed using a ubiquitous joint constitutive law with strain softening. The applied procedure is capable of reproducing the macroscopic stress and strain evolution in the sample during triaxial deformation until a failure process occurs. In addition, we calculated focal mechanisms of acoustic emission (AE) events and resolved local stress field orientations. This detailed stress information was compared with that from numerical modelling. The comparison was made based on the 3D rotation angle between the cardinal axes of the two stress tensors. To infer the differences in rotation, we applied ANOVA. We identified the two time levels as the plastic deformation phase and the after-failure phase. Additionally, we introduced the bin factor, which describes the location of the rotation scores in the rock sample. The p values of the test statistics F for the bin and phase effects are statistically significant. However, the interaction between them is insignificant. We can, therefore, conclude that there was a significant difference in the time between the rotation means in the particular bins, and we ran post hoc tests to obtain more information where the differences between the groups lie. The largest rotation of the stress field provided by the focal mechanisms of AE events from the numerically calculated stress field is observed in the edge bins, which do not frame the damage zone of the sample.

approaches have been proposed to invert earthquake focal mechanisms for stress orientation (Maury et al. 2013), but the most popular are those developed by Michael (1984), Gephart and Forsyth (1984) and Angelier (2002) with extensions proposed by Lund and Slunga (1999), Hardebeck and Michael (2006), Arnold and Townend (2007), Maury et al. (2013), Vavrycuk, (2015 and others. Usually, the more widely used methods obtain similar results for similar data sets. However, the stress inversion results are sensitive to the number of earthquakes, focal mechanism uncertainties, and fault plane variability (e.g., Hardebeck and Hauksson 2001;Bohnhoff et al. 2004;Vavrycuk 2015). To guarantee the high resolution of the stress inversion results, a high number of seismic events with well-constrained focal mechanisms is required in a spatially determined area over a certain time period. These conditions are well satisfied by laboratory experiments, which provide a convenient framework to test the performance of stress inversion techniques in determining the details of stress field variations during loading. Triaxial compression (TC) tests on rock samples performed to analyse rupture mechanics are frequently accompanied by monitoring acoustic emission (AE) activity. The analysis of AE activity allows the spatiotemporal evolution of damage in the sample to be tracked and provides detailed information on fracturing and frictional processes in rock samples subjected to loading (Stanchits et al. 2006). The seismic moment tensor (MT) can be inverted from AE data and then decomposed to describe volumetric (ISO), double-couple (shear, DC) and compensated linear vector dipole (CLVD) strain components. This helps improve the understanding of physical processes taking place within seismic sources and the sample, such as rupture dynamics, fault complexity or the radiation of the seismic energy related to the damage accumulation (Ben-Zion and Ampuero 2009;Castro and Ben-Zion 2013). Based on the MT solutions, the focal mechanisms of AE events are estimated.
The stress inversion method provides information about the directions of the three principal stress axes and a measure of the size of the intermediate principal stress, r 2 , relative to the maximum, r 1 , and minimum, r 3 , principal stresses, called the stress ratio R. The stress ratio is very often used to determine temporal local rotations of the stress tensor and to determine the potential processes responsible for these variations.
Systematic temporal stress rotations have been observed in reservoirs in relation to fluid injections (Martínez-Garzón et al. 2013, 2014aSchoenball et al. 2014). These stress variations appear as a response to the pore pressure changes and decrease in in situ temperatures by the cold fluid [Jeanne et al. 2015;Yoon 2015]. More recently, stress rotation has been identified in subduction zones in relation to slow slip events (Warren-Smith et al. 2019). These stress changes were interpreted as the accumulation and release of fluid pressure within the subducting oceanic crust, impacting the timing of slow slip event occurrence. Spatiotemporal local stress tensor rotations were also identified before and after large tectonic earthquakes [Hardebeck and Hauksson 2001;Hardebeck 2012;Ickrath 2015], indicating that earthquakes are capable of causing significant stress partitioning along their rupture [e.g., Bohnhoff et al. 2006].
Local rotations of the stress field at a fault are very difficult to detect. They are of a slightly higher order than the error of the stress field estimation. Theoretically stress axes may rotate up to 45°due to an earthquake (Hardebeck and Hauksson 2001), and temporal rotations [ 20°have been observed (e.g. for the Tohoku earthquake, Hasegawa et al. 2011). Different spatial zones may have completely different stress states, so even larger spatial rotations are possible.
Therefore, it is important to recognize whether the possible range of stress axis rotations along the fault is statistically significant. Here, we exploit laboratory tests (Kwiatek et al. 2014) to develop a study to characterize the evolution of the stress field along the rupture in a sample during the whole laboratory experiment. We develop a numerical modelling method to reproduce the evolution of stress and strain changes in the classical triaxial fracture experiment performed on westerly granite (WG) samples. We did not intend to reproduce the global stress and stress evolution in the sample but to reproduce the peculiarities of the local stress and strain field. Numerical modelling delivers a more detailed picture of the spatiotemporal evolution of stress and strain in the rock sample, which is typically only available from associated macroscopic measurements (axial load, axial displacement) or point measurements on the sample surface (strain meters). Then, we discuss the implications of our findings in the context of stress rotation monitoring, spatiotemporal stress partitioning along the rupture, and the capability of numerical modelling to capture these stress variabilities. This work contributes to an improved understanding of the physical mechanics underlying the rupturing process and assesses the possibilities to monitor stress heterogeneities in the rupture region. The study aims to provide highlights of answers to the following questions: What happens to the stress field tensor during the experiment? When and which stress rotation pattern may be observed during the preparatory rupture process, coseismic phase and postseismic phase, what is the magnitude of this rotation and is it statistically significant? Does numerical modelling give a reliable representation of local spatiotemporal changes in the stress field orientation? To answer these questions, the numerical modelling results of the spatial stress field orientation are compared with local stress field data inverted from seismic data. The AEderived focal mechanisms (e.g., Kwiatek et al. 2014) are used to invert the local stress field orientation using stress tensor inversion algorithms (Martínez-Garzón et al. 2014a, b).
The paper is organized as follows: first, we present the experimental procedure of the triaxial compression test, discuss the observed seismicity and then describe the MT and stress inversion procedure. Next, the numerical simulation based on Itasca FLAC3D, a three-dimensional finite difference method (FDM) software, is presented; this section involves the calibration of the experimental and synthetic stressstrain curves. Then, we compare local stress field orientations originating from stress tensor inversion of AE events with those derived in the model. To compare the observed and synthetic stress orientations, rotation angles (Kagan 2007) were used, and analysis of variance (ANOVA), a statistical tool, was used to quantify the differences.

Experimental Procedure and Seismicity
The triaxial test was performed on a cylindrical, intact sample of WG. The sample size was 40 9 107 mm, and it was loaded with a constant strain rate of 3 9 10 -6 s -1 . The oven-drained sample was notched 2.5 cm deep at 30°to the cylinder axis. The laboratory experiment as well as AE measurement procedure are described in detail in Kwiatek et al. (2014).
International Society for Rock Mechanics (ISRM) suggests that the specimen diameter should not be less than 54 mm for the cylinder specimen. Authors would like to add that the other ISRM recommendations have been met: • the test specimens shall be right circular cylinders, • height to diameter ratio of 2.5:3.0 (2.68), • the diameter of the specimen should be related to the size of the largest grain in rock by the ratio of at least 10:1 (largest grain around 300 lm; The grain size of Westerly granite ranges from about 0.05-2.2 mm (Moore et al. 1987), so even taking into the consideration the biggest values mentioned, the diameter of analysed specimen is over 18 times bigger).
First, the sample was confined to 75 MPa. Then, deviatoric loading was applied at a constant displacement rate of 3 9 10 -6 s -1 until the sample fractured (Fig. 1a). The generated fracture was complex at the top part and was composed of two major subsurfaces, as presented by the spatial distribution of the AE events (Fig. 1b). The displacements (strains) in the sample were monitored by two strain meters located directly on the sample. The sample failed at a maximum vertical stress equal to 574 MPa (Fig. 1a), when the strains in the sample reached a maximum value of 1.19 (Fig. 1a).

Acoustic Emission Monitoring and MT Inversion
The AE activity was monitored during the experimental test with 16 AE sensors glued to the sample surface, ensuring almost full azimuthal coverage of the seismic events. AE monitoring resulted in 14,583 events detected and located within the sample (Figs. 1b, 2-blue bars). In Fig. 1b hypocenters of AE activity recorded during the loading were presented. All dots correspond to seismic events. Additionally, moving with 50 s window, we have presented seismic events occurred at the maximum loading equal to 574 MPa with yellow color. That procedure allows for highlighting the fault plane that appeared due the applied loading.
Additionally, for 6561 events, full MT inversion (FMTI) was performed (Fig. 2, light bars) following the procedure described in Kwiatek et al. (2014). For  Histogram of seismic activity during triaxial compression tests (dark bars-events with the hypocentre location, light bars-events with the hypocentre location and seismic moment tensor); distribution of AE magnitudes with the occurrence time of the events (stem plot); seismic events FMTI, 14 first P wave amplitudes were used. The most seismically active period started just before sample failure. The largest AE events occurred close to the moment of fracture initiation (Fig. 2, stem plot). The blue background represents the period when the sample was in the postfailure phase after reaching the maximum compressive strength. The increasing number of seismic events was connected with increasing compression stress that was applied to the bottom plate of the rock sample, with the peak seismic activity occurring just after the sample failed. The most seismically populated volume is located close to the sample centre (Fig. 3, dark bars); the highest magnitude events are also located centrally in the sample and in the upper half (Fig. 3, stem plot).
The moment tensor can represent different seismic sources and to identify its type, the MT is usually diagonalized and decomposed into elementary parts: deviatoric and isotropic. The isotropic part describes the volumetric strain component. By definition, a positive isotropic tensor is related to volumetric expansion. The deviatoric part can be additionally divided into double-couple or CLVD components; MT components are widely used in seismology for physical interpretations. Seismic source components have already been described in many papers (see, e.g., Vavrycuk 2015). The DC source is generally connected with shear-faulting. However, a CLVD source has no verified geological meaning and is often interpreted as a source component representing the residual radiation from the best DC source (Dahm and Krueger 2014). In seismology, especially in anthropogenically induced seismicity, seismic sources can have different MT compositions and can be complex (e.g., Orlecka-Sikora et al. 2014;Lizurek et al. 2015;Rudzinski et al. 2016Rudzinski et al. , 2017Lasocki et al. 2017;Lasocki and Orlecka-Sikora 2020).
The vast majority of the analysed AE events have a double-couple component of the source mechanism, indicating that shear slip occurs. These are faultparallel AE events with dip directions in accordance with macroscopic slip (Kwiatek et al 2014). The authors also noticed that larger AEs contain fewer ISO components, whereas small events contain more ISO components. There are 3 times more DC events than CLVD events and almost 9 times more DC events than ISO events (Fig. 4).

Stress Tensor Inversion
Stress tensors were inverted from AE data using two inversion methods: MSATSI and BRTM. MSATSI (Martínez-Garzón et al. 2014a, b) provides a framework for calculating the deviatoric stress tensor together with its uncertainties using the bootstrap resampling method. The BRTM method uses the right tetrahedral method and Bayesian approach for the determination of the stress field from focal mechanism datasets and provides a probability function over the focal sphere for both the r 1 and r 3 principal stress directions (Massa et al. 2016).
The WG sample was divided into nine spatial bins in which stress inversion was performed from available AE mechanisms (Fig. 5a). We considered two Fig. 3 Histogram of the seismic activity (bars of histogram) and distribution of the AE magnitude (stem plot) as a function of the sample height phases of sample deformation under the triaxial compression test: plastic and postfailure phases. The elastic phase was omitted due to the very limited number of seismic events during this phase. We consider the plastic phase as a period from the beginning of the experiment until the maximum stress of r 1 (574 MPa) is reached. The post-failure phase is the consecutive one. Generally, at least 24 event windows were chosen for stress inversion calculations. Only for one case, namely, bin 7, during the failure phase was a 17-event window chosen for the same reason as for the rejection of the elastic part.
Such event windows ensure the largest number of cases (220 stress inversions) where the reliable limit of the seismic events for stress inversion was fulfilled. The number of seismic events and stress inversions performed in particular bins is presented in Table 1. Italic cells correspond to the plastic phase of the deformation of the WG sample, and the failure phase is highlighted in bold italics.  green-r 1 , yellow-r 3 ) and MSATSI (marked with triangles; dark-r 1 , light-r 3 ) software in all nine bins; c variability of the dip angles of r 1 and r 3 within all bins obtained in BRTM; dark lines correspond to r 1 , whereas light lines correspond to r 3 The comparison of the trends and plunges of the maximum (r 1 ) and minimum (r 3 ) principal stresses obtained from stress inversion performed in BRTM and MSATSI software for all nine bins (the bin indicator is marked in the right corner of each plot), and just one calculation window is presented in Fig. 5b. Both inversion methods give similar results, except for bins 4 and 6. Figure 5c presents the variability of r 1 and r 3 within all bins for the whole process of triaxial loading obtained by BRTM; dark lines represent r 1 , whereas r 3 is represented with light lines. The median dip angles of the principal stresses are equal to 72.6°and 9.7°for r 1 and r 3, respectively.

Numerical Modelling
To answer the questions raised in the introduction, the FLAC 3D model of WG rock samples under triaxial compressive pressure was developed to model the spatial and temporal evolution of stresses and strains within the sample. Thereafter, the model was calibrated using the data from the actual laboratory experiment. This is to justify the stress changes originating from the seismic observations gathered from AE sensors and the seismic moment tensor calculations from the laboratory triaxial experiment. The calibration process in numerical modelling involves explicit fine-tuning of the geomechanical parameters in the numerical model to achieve a significant level of agreement with the experimental data (Zeh-Zon Lee 2014).

Geometry and Boundary Conditions
Cylindrical specimens with a circular base with a radius of 20 mm and a height of 107 mm were modelled using hexahedral-shaped zones. The rock specimen model consisted of 42,400 elements and 43,254 grid points.
The boundary conditions applied to the model consist of fixed vertical displacements at the top of the sample and at the surface contacting the loading plate. The displacements in the horizontal directions are allowed, ipso facto, implementing zero friction. Generally, special low friction materials are used in practice to reduce the base and loading plate friction, which justifies the application of the previously mentioned boundary conditions. Additionally, friction added to the sample top could result in the increased strength of the sample, whose origin would be difficult to distinguish later: this could be influenced by the applied friction or by the influence of intermediate principal stress (Senent et al. 2013).

Constitutive Model
Rock responses to compressive loading with different phenomena: elastic deformation, microfracture initiation, plastic deformation and finally failure. The stages before the failure occurrence are characterized by the elastic response, while the post-failure of the brittle rocks, such as the granite, is characterized by massive crack propagation leading to shear bands. The mentioned behaviours are strongly influenced by the amount of confining pressure (Tan et al. 2015).
There are several constitutive models since the late 1950s, when researchers started to use plastic theory in constitutive models for rock specimens. The most commonly used Mohr-Coulomb failure criterion is still evolving, combining ever newer behaviours such as cohesion softening-friction hardening and adding different parameters to control the strength degradation behaviour of the rock (Tan et al. 2015).
The typical axial stress-strain curves, depending on the degradation behaviour, are presented in Fig. 6. The intact strength of the rock and the presence of joints strongly govern the strength and deformation of a jointed specimen. The constitutive model used in this study is the ubiquitin joint model with strain softening (SUBI). Strain softening was reached using cohesion softening-friction hardening behaviour. It accounts for the presence of orientation of weakness (weak plane) in a Mohr-Coulomb model. In the SUBI model, yield can occur in either the solid or along the weak plane or both (FLAC3D, Version 5.0, User's Guide). The presence of the joint is accounted for in the plastic corrections but has no effect on the elastic behaviour, and the model is restricted to one set of joints (Dehkordi 2008;Ismael and Konietzky 2017). When general failure is detected within Itasca's FLAC software and plastic corrections are applied, then new stresses are analysed for failure on the weak plane and updated accordingly (FLAC3D, Version 5.0, User's Guide).
The criterion for failure on the plane, whose orientation is given, consists of a composite Mohr-Coulomb envelope with tension cut-off. The position of a stress point on the latter envelope is controlled by a nonassociated flow rule for shear failure and an associated rule for tension failure. The SUBI constitutive model was chosen due to the occurrence of the two notches created before the loading was applied to the rock sample and not modelled explicitly in the sample. The SUBI model allows for specifying both matrix and joint properties.

Geomechanical Parameters
Although the strain softening constitutive model allows for changes in parameter values with plastic strain, the stiffness moduli are independent of plastic strain. The majority of the published values, such as the elastic, shear and bulk moduli, density and Poisson's ratio, can be used as-is in the constitutive model, but the cohesion, friction and tensile strengths are more difficult to directly relate between the laboratory test values and modelling values. The hardening and softening parameters for this case were back-calculated from the results of the laboratory triaxial test (Tables 2, 3).

Simulation of Triaxial Loading
To model the triaxial loading applied to the sample, the following modelling steps were applied to the model: Step 1 Boundary conditions were applied to the model geometry. Then, the hydrostatic stress state was initialized at 75 MPa.
Step 2 Confining pressure was applied to the outer faces circumfluent to the cylinder.
Step 3 A constant strain rate was applied in the vertical direction at the bottom end of the cylinder, leading to an increase in axial stress until failure occurred. The model was loaded by applying a constant grid velocity of 3.21 -8 m/step (determined based on the constant strain rate 3 9 10 -6 s -1 that was obtained during experiment in the lab and the sample height).
Additionally, servo controlling loading was involved in the numerical calculations. When the failure mechanism is initiated, the stress state in the sample becomes nonuniform. To better control the deformation in the system and to reduce numerical errors in the modelling, the magnitude of the strain rate can be monitored and subsequently adapted as a The FLAC3D software allows for calculating the stress distribution within the model with the accuracy specified by the user. Accuracy is provided to the model by means of densification of the mesh grid. Stress information is stored in each zone of the numerical model. The principal stress values were monitored at nine points within each of the nine bins of the WG sample: eight monitoring points were located in the bin corners, and one was centrally located.

Global Stress Relations
A low (0.036) RMSE was achieved for the fitting of the stress-strain curve obtained from the triaxial compression test made in the laboratory, and the curve resulting from numerical modelling (Fig. 7).
The experimentally observed and numerically derived peak strengths of the sample reach 574 MPa and 584 MPa, respectively (with a difference of 1.62%). In both the experimental and numerical tests, sample failure occurred at similar axial strains, which were equal to 1.19 and 1.15, respectively (difference of 3.36%). In the numerical simulation, slight strength hardening was observed after failure. Analysing the appearance time of yielded zones during numerical loading, these zones start to become visible during the strength-hardening period, so we can assume that the failure appears at slightly greater strains. The FLAC3D software also allows the development of localization of the shear bands to be tracked by plotting zones in which plastic yielding occurred. In Fig. 8, a comparison between the failure pattern created during the experiment and numerical modelling is shown. The fracture plane generated during laboratory experiments is visibly complex, revealing a broad damage zone composed of two subsurfaces with a rough fracture surface. The numerical modelling shows the yielded zones in shear (Fig. 8b). The created shear band is similar to the failure pattern obtained from the laboratory test. Figure 8b shows the vertical displacement modelled on the fault plane reaching the maximum value of 2.53 mm. The piston displacement observed during the experimental trial reached a value of 2.52 mm at the moment of sample failure, which shows a significant convergence with the modelled value. The distributions of other parameters within the WG sample, such as the maximum, intermediate and minimum principal stresses along with maximum shear strain rate or vertical, horizontal Z and horizontal X displacement, are presented in the auxiliary materials (Appendix 1).

Results and Discussion-Local Stress Relations
We compared the stress orientations derived from AE data and modelled them with FLAC3D software. The comparison was made based on the 3D rotation angle (Kagan 2007) between the cardinal axes of the two stress tensors. We calculated the angles at which one principal stress coordinate system (e.g., from AE stress inversion) has to be rotated to obtain the remaining one (numerically modelled) to find where two stress fields deviate from each other the most.
To infer the differences in the rotation of the stress field derived from the stress inversion based on focal mechanisms of AE data and the stress field calculated numerically, we apply analysis of variance (ANOVA). ANOVA compares mean scores with each other to detect any overall differences between the results. First, we identified the two phases as the plastic deformation phase and the after-failure phase to check if the rotations were higher before or after failure. Additionally, we introduced the already mentioned bin factor, which describes the location of the rotation scores in the rock sample. We tested for the null Fig. 8 Post-mortem cross-section of the sample. The damage zone is filled with epoxy; a The broad damage zone was composed of two sub-surfaces, b the distribution of the vertical displacement along the yielded zones, and the c not-yielded zones to present the created fracture hypothesis that there are no differences in the mean rotation values with respect to time phases or spatial bins. The alternative hypothesis states that the related mean values of rotation are not equal; at least one mean is different from another mean. We calculated the F statistics, which is the ratio of two variances, between-groups and within-groups, that are expected to be approximately the same value when the null hypothesis is true, which yields F-statistics near 1. The F-statistics calculated based on our sample are compared to the critical F-value of the F-distribution for a population where the null hypothesis is true at a certain significance level (usually the 0.05 level). We evaluated whether our sample F-value was so rare that it justified rejecting the null hypothesis for the entire population. The probability of obtaining an F-value that is at least as high as our study's value is the p value, which indicates the level of statistical significance. A low probability, low p value, indicates that our sample data are unlikely when the null hypothesis is true, and the null hypothesis is rejected. Table 4 shows the results of the ANOVA for the within-subjects effects. The p values of the test statistics F for the main effects-bin and phase-are 0.0283 and 0.0008, respectively, and both effects are statistically significant. However, the interaction between them is insignificant, since F is 1.11 and p is 0.3519. We can, therefore, conclude that there was a significant difference in the time between the mean values of rotation in the particular bins, and we ran post hoc tests to obtain more information about where the differences between the groups lie.
The Tukey honest significant difference (HSD) test was used to determine if the interaction between two sets of data was statistically significant. The value of the Tukey test is given by dividing the absolute value of the difference between pairs of means and by the standard error of the mean (SE). The post hoc Tukey's test confirmed that the rotations of the 2nd bin are significantly lower than the rotations observed in other bins (except the 1st and 5th-7th bins). When we look at the phase of deformation, we see that the rotations during the failure phase become significantly higher in comparison to the plastic phase. The analysis of the interactions between the main effects provides insights into the behaviour of the stress field in the particular bins during the experiment. It is visible that the rotation of the stress field is affected by the phase of deformation for most of the bins. However, for the failure phase for the 5th bin, the stress field is significantly less rotated compared to that of the modelled one. For the last phase of deformation, when the damage of the sample appears, the rotation does not increase, contrary to the rest of the bins.
From Fig. 9, it can be concluded that the largest rotation of the stress field provided by the focal mechanisms of AE events from the numerically calculated stress field is observed at the edge of bins 1 and 3 and bin 7, which do not frame the damage zone of the sample. Relatively large distances were present to the failure plane and additionally, for bins 1 and 7, the lowest number of inversions performed during both the plastic and failure phases, could result in such rotations. The lowest rotations are observed in the central bins of the sample, numbered 2 and 5, and the rotation in the upper right corner, namely, bin 9. Such rotations for bins 5 and 9 can possibly be explained by the largest number of inversions performed within these bins (more reliable results) as well as its location regarding the newly created fractures. However, not only the location to the fracture plane but also the rotations can be lower due to a more stable stress field, which is related to the offset from the model boundaries and hence the boundary conditions. The variation in the trends and plunges of r 1 in all the sample bins from the AE in comparison to the numerically modelled bins is presented in Fig. 10. Coloured circles present strikes and dips for the particular bin (the bin number is presented inside the circles) retrieved from AE data, whereas the numerically modelled data are presented with black crosses. The latest mentioned could be presented with just one point due to the same plunge of r 1 for all bins, which is additionally equal to 90°. In this case, different values of the r 1 strike do not influence the cross position. Vivid colours correspond to the plastic phase of deformation, while weathered colours present a state of failure. During the plastic deformation for the two bins, namely, the 4th and 7th bins, trends of r 1 are in the range of the 2nd (3 bins) and 3rd (1 bin) quarters of the spheroid, whereas the rest of the trends (5 bins) lay in the 1st and 2nd quarters. In the failure state, the r 1 strike for two bins lay in the 1st quarter now, with four bins in the 2nd quarter, two in the 3rd and one in the 4th.

Conclusions
This paper presents a detailed approach for calibrating a strain-hardening/softening ubiquitin-joint model (based on the Mohr-Coulomb model) for simulating the stress-strain behaviour of WG samples. The procedure was presented by calibrating the model for triaxial testing data. Successful calibration of global stress-strain behaviour of the numerically modelled sample with the experimental data allowed us to attempt to compare local stress tensors coming from two numerical modelling and stress tensor inversion of AE events.
The calibration of the specimen shows the following results: • The bilinear ubiquitin-joint model with friction hardening/cohesion softening is capable of reproducing Westerly Granite behaviour during the triaxial fracture experiment.  Additionally, the distribution of the vertical displacement along the yielded zones creating the failure plane within the Westerly Granite sample is presented in c Fig. 11 Distributions of the a vertical, b horizontal Z and c horizontal X displacements within the sample • We identify the two time levels as the plastic deformation phase and the after-failure phase. The p values of the test statistics F for the main effects-bin and phase-are statistically significant. However, the interaction between them is insignificant. We can, therefore, conclude that there was a significant difference in the time between the rotation means in the particular bins, and we ran post hoc tests to obtain more information about where the differences between the groups lie.
• The largest rotation of the stress field provided by the focal mechanisms of AE events from the numerically calculated stress field is observed in the edge bins, which do not frame the damage zone of the sample. The lowest rotations are observed in the central bins of the sample. • The rotation of the stress field is affected by the phase of deformation for most bins. Rotations during the failure phase become significantly higher in comparison to the plastic phase. The post hoc test confirmed abovementioned observations. Funding This work is partially funded by Science4CleanEnergy (S4CE), a European research consortium funded by European Union's Horizon 2020 research and innovation programme under grant Agreement No. 764810.
Availability of Data and Material The datasets generated numerically during the current study are available from the corresponding author on reasonable request. The authors would like to thank Grzegorz Kwiatek from Deutsche GeoForschungsZentrum (GFZ) for providing triaxial experimental data as well as for his invaluable comments and advice on this paper.

Declarations
Conflict of interest All authors have participated in (a) conception and design, or analysis and interpretation of the data; (b) drafting the article or revising it critically for important intellectual content; and (c) approval of the final version. This manuscript has not been submitted to, nor is under review at, another journal or other publishing venue. The authors have no affiliation with any organization with a direct or indirect financial interest in the subject matter discussed in the manuscript.
Informed Consent Informed consent was obtained from all individual participants included in the study. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References
Angelier J (2002) Inversion of earthquake focal mechanisms to obtain the seismotectonic stress IV -a new method free of choice among nodal planes. Geophys J Int 150:588-609 Fig. 15 Differences in the principal stress values monitored at 9 points within each of the 9 bins within the WG sample. The colours of the curves correspond to the bin which is being compared to the 5th bin and each curve combines recordings from all 9 points within a particular bin