Rock Properties and Modelled Stress State Uncertainties: A Study of Variability and Dependence

The safety and sustainability of subsurface applications requires a profound knowledge of the local stress state which is frequently assessed using 3D geomechanical-numerical models. Various factors lead to generally large uncertainties in these models. The variabilities in the rock properties as one of the sources of uncertainties and their influence on the modelled stress state is addressed herein. A generic 3D geomechanical-numerical model is used to investigate the influence of different distributions of variability and their effect on different stress states. The variability in rock properties clearly affects the uncertainties in the stress state in a positive correlation with differences that depend on the affected component of the stress tensor. The basic observation is that largest uncertainties are observed in the normal components of the stress tensor where the variabilities apparently are most effective. The same rock property variabilities affect the shear components uncertainties to a significantly lesser extent. Variabilities in the Young’s modulus and the Poisson’s ratio chiefly affect the uncertainties in σxx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{xx}$$\end{document} and σyy\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{yy}$$\end{document}. The density variability, however, leads to highest uncertainties in σzz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{zz}$$\end{document}. In general, variabilities in the Young’s modulus are most effective, followed by the Density and then the Poisson’s ratio. Furthermore, an influence of the tectonic stress regime on how the variability in the rock properties affects the stress state is observed. At the same time only a small effect is observed for different stress magnitudes. The eventual uncertainties in a modelled stress state depend not only on the uncertainties in the rock properties but also whether the uncertainties are found mainly in the Young’s modulus, the Poisson’s ratio or the Density. These findings indicate the importance to regard variabilities in rock properties as a source for significant uncertainties in geomechanical-numerical models. It is proposed to use the derived relations for an inexpensive quantification of uncertainties by means of a post-computation assignment of uncertainties to a stress model. The influence of rock property variability on modelled stress state uncertainties is investigated using generic 3D geomechanical-numerical modelling. The variability in the Young’s modulus has the largest effect on the uncertainties in the stress state. The uncertainties in the normal components of the stress tensor are primarily affected by rock property variability. The tectonic stress regime and the absolute values of rock properties affect the results. The influence of rock property variability on modelled stress state uncertainties is investigated using generic 3D geomechanical-numerical modelling. The variability in the Young’s modulus has the largest effect on the uncertainties in the stress state. The uncertainties in the normal components of the stress tensor are primarily affected by rock property variability. The tectonic stress regime and the absolute values of rock properties affect the results.


Introduction
Subsurface geotechnical applications require a reliable description of the 3D stress state to ensure safe and sustainable operations. This can prevent economic losses related to damages of the infrastructure or due to loss of support from stakeholders after the occurrence of induced seismicity.
Information on the stress state is provided by a wide range of stress indicators (Amadei and Stephansson 1997) that probe very different rock volumes (Ljunggren et al. 2003). The volume ranges from 10 −2 -10 2 m 3 for borehole related stress indicator up to 10 8 m 3 for fault slip data or earthquake focal mechanisms (Amadei and Stephansson 1997). However, stress measurement campaigns indicate that even within small rock volumes the stress state can be highly variable (Martin and Christiansson 1991;Cuisiat and Haimson 1992;Singha and Chatterjee 2015;Liu et al. 2020;Feng et al. 2020). Assuming that this is related to the variability of elastic rock properties, leads to a generally small representative elementary volume (REV), i.e. the volume in which the rock properties are assumed to be homogeneous (Amadei and Stephansson 1997).
To predict the stress state 3D geomechanical-numerical models are used (Henk 2009;Stephansson and Zang 2012;Singha and Chatterjee 2015;Liu et al. 2020). Such a model in its basic form requires a geological model with distinct rock units that are assigned appropriate rock properties (Van Wees et al. 2003;Ahlers et al. 2021;Andersen et al. 2022). Within the model volume the stress state variability is controlled by the overburden and differences in the lithology which are a result of a complex succession of rock units and/or heterogeneities in the rock properties of a single unit (Warpinski and Teufel 1991;Martin and Chandler 1993;Hergert et al. 2015). While the succesion of rock units is addressed by accordingly designed geological models, the heterogeneities require information on the rock properties, their variability and the according REV which are divided in individual levels of heterogeneity, one for each rock property (Amadei and Stephansson 1997). This information is usually not available to a sufficient degree of detail since only a few-if any at all-rock property measurements are available (Leijon 1989;Amadei and Stephansson 1997). As a result, the modelled stress state is usually subject to significant uncertainties due to the unknown variability of rock properties (Hergert et al. 2015;Ziegler et al. 2016). This can be addressed by an assignment of rock properties to REV based sub-volumes of a model instead of assignment of fixed rock properties for an entire unit. An established approach is to derive rock properties and their distribution from seismic observations of dynamic rock properties that are transferred to static ones (Bosch et al. 2010;Trudeng et al. 2014). This results in the desired heterogeneity in the rock properties. However, the sub-volumes sizes are either a result of the resolution of the seismic survey, the cells used during inversion, or the geomechanical models discretization. They are thus not true REVs based on the rock properties distribution. In addition to many assumptions during the conversion from dynamic to static rock properties which decreases significance (Mavko et al. 2009;Fjaer 2019), no uncertainties associated with the heterogeneities are quantified. Furthermore, often only a small number of rock property measurements are available, which inhibits both, assignment and a meaningful analysis. In order to compile any information on the uncertainties, unit wide changes in the rock properties are used as a proxy for measurement errors (Hergert et al. 2015;Ziegler et al. 2016). Even though this approach captures some uncertainties and indicates a rough estimation of the range in expected stress states, it is not an appropriate approach to assess the full uncertainties due to the missing of intra-unit heterogeneity. Thus, the current geomechanical-numerical modelling workflows disregard the full range of uncertainties in the stress state due to rock property variability. Any modelled stress state hence contains hidden uncertainties since neither the qualitative nor the quantitative influence of the rock properties variabilities on the stress state are known (Ziegler et al. 2016;Fraldi et al. 2019;Harper 2020). This challenge is frequently addressed by experiments backed by numerical models on the influence of heterogeneity on the failure process (Tang et al. 2000;Wong et al. 2006;Lei and Gao 2019) or by studying the heterogeneity of the stress state in local highdensity stress measurement campaigns which can be seen, amongst others, as a result of heterogeneous rock properties (e.g. Martin and Chandler 1993;Wileveau et al. 2007;Schoenball and Davatzes 2017;Feng et al. 2020).
Herein, the perspective is reversed by starting with the rock property variability in order to investigate the impact of intra-unit variabilities of the rock properties on the six independent components of the modelled 3D stress tensor. A generic 3D geomechanical-numerical model with a single lithological unit is used. In order to simulate the variability of rock properties, they are individually assigned to subvolumes that are assumed to be REVs. The variability is realised by definition of probability distributions for each rock property, herein the following three elastic properties: Young's modulus, Poisson's ratio, and Density. A normal distribution is assumed for all three rock properties to account for possible outliers. Then, the rock properties are drawn from the according distribution for each subvolume individually. This leads to random assignment to sub-volumes of the three rock properties that are assumed to be independent from each other. The result is a stress state which is variable in a way that corresponds to the rock properties' variability. The aim was to investigate the uncertainties in the stress state due to the variability in the rock properties. Therefore, sufficient models with different randomly assigned rock properties from the same distributions are computed in order to estimate the range of possible stress states. Thereby, the expected uncertainties in the six components of the modelled stress state are found. Furthermore, changes in the rock properties, their variability, the stress magnitudes, and the tectonic stress regime allow exploration of their respective influence on the resulting modelled stress state uncertainties. The results show a first-order linear dependence of the modelled stress state uncertainties on the rock property variability. The rock stiffness (Young's modulus) is the key driver and the stress tensor's normal components have significantly higher uncertainties in comparison to the shear components.

General Concept and Model Setup
A 3D geomechanical-numerical finite element model to estimate the influence of the variabilities in the following three rock properties regarded in mechanical models:Young's modulus, Poisson's ratio, and Density-on the modelled stress state is used. The applied generic model with a single lithological unit has an extent of 10 × 10 × 5 km 3 and is discretised with 3000 2nd order brick elements which for numerical reasons are distinguished into 375 sub-volumes (see "Appendix") of which each is assumed to be a REV. By assigning different and independent rock properties to each sub-volume, their variability is realised. The respective properties are randomly drawn from normal distributions defined by a mean value and a standard deviation, i.e. different and independent values for the Young's modulus (pictured in Fig. 1), the Poisson's ratio, and the density are assigned. The now realised distribution of rock properties throughout the model volume is only one possibility out of many different scenarios. As such, the assigned rock properties throughout an individual model scenario do not perfectly follow the according distribution. Instead, the full distribution is realised by generation of various model scenarios (see Sect. 2.3). The model describes the equilibrium of volume (gravity) and surface forces (tectonic stress). The latter is imposed by displacement boundary conditions which are varied to find the best-fit (Reiter and Heidbach 2014) to a synthetic stress state at a single point at a given depth in the model volume ( Fig. 1 and Sect. 2.3).

Distribution of Rock Properties
The variability experienced in the properties of the same lithology has different origins. Partly it is due to inherent inhomogeneities in the rock (Martin and Chandler 1993;Sone and Zoback 2013;Zhang et al. 2020) and partly due to uncertainties that stem from deficiencies in our ability to assess the rock properties at all and in particular throughout the entire volume of a rock unit (Amadei and Stephansson (1997) and references therein). The three mechanically relevant rock properties that are discussed herein-in particular the Young's modulus-implicitly contains variability due to additional rock properties that are not relevant for this approach of mechanical modelling, e.g. porosity or rock strength. Thus, the approach is not limited to a certain reservoir type but universally applicable.
In the following, the five main reasons for variability in rock properties are indicated. (1) Systematic lateral changes in the lithology due to fluent transitions may remain undetected due to limited measurements. In addition, 3D geomechanical-numerical modelling usually simplifies a variable lithology by a single set of rock properties. (2) A certain random heterogeneity in the lithology is expected in most cases which adds noise to the rock properties. (3) If a rock property measurement is available, e.g. from a core sample, it is unclear whether the core is representative for a significant part of the lithological unit or if it is an exceptional outlier. The variability then has to be estimated from the range of values found in various samples from different locations. If only one sample is available, expert elicitation comes into play to assign a likely variability which increases the uncertainties. (4) If no rock property measurements of the considered lithology are available, often reference datasets of the same rock type but from a different location are used. This bears the possibility for significant deviations which have to be met by the assignment of a significant variability. (5) Eventually, measurement errors and uncertainties due to necessary assumptions are added to the results of laboratory tests.  The three elastic properties of rocks that are discussed herein-in particular the Young's modulus-showcase a great variability that depends, amongst others, on the rock type, location, weathering, or burial history (Phoon and Kulhawy 1999;Cai 2011;Zhang et al. 2020). A possible correlation of the rock properties with or a dependency on each other is disregarded herein. No universally applicable correlation or dependency is known. As default average values used to define the center of the distribution, values are chosen that represent average rocks-a Young's modulus of 50 GPa, a Poisson's ratio of 0.28, and a density of 2700 kg/ m 3 (Aladejare and Wang 2016; Bär et al. 2020).
Rock properties and their variability are studied with regards to several issues. These range from literature studies that compile a very broad view on different rock typesigneous, sedimentary, and metamorphic-and their variability (Aladejare and Wang 2016;Bär et al. 2020) to detailed local studies that explore the rock property variability in an individual location and lithology (Pepe et al. 2016;Zhang et al. 2020). These studies results led to the assignment of the variabilities to the rock properties. In contrast to other studies, where a Weibull type distribution is found for the rock strength (Lobo-Guerrero and Vallejo 2006;Wong et al. 2006) herein, a normal distribution is applied, where a variability of 1 includes roughly 68% of the rock properties values. This is due to the comprehensive regard for reasons for variability addressed in this study which involves measurement errors. Furthermore, in most field applications the few data records on rock properties disallow a reliable estimation of a distribution. Due to the exclusion of extreme outliers by definition of a range the normal distribution is chosen to represent the rock properties. The range of 1 tested for the influence on the stress state is from 10% up to 60% of the average value, e.g. the average Young's modulus of 50 GPa is tested with a variability of 1 between 5 and 30 GPa.

Model Calibration and Uncertainty Estimation Due to Rock Property Variability
The presented generic model features an assumption on the stress state which is realistic according to various stress measurement campaigns and data compilations (Brudy et al. 1997;Hickman and Zoback 2004;Morawietz et al. 2020). In a depth of 2900 m a strike slip stress regime with xx = 100 MPa, yy = 48 MPa, and zz = 77 MPa is assumed. By varying the displacement boundary conditions the model is calibrated at a single sub-volume at that depth ( Fig. 1, Reiter and Heidbach 2014). The rock properties in this subvolume used for calibration are fixed to the mean values of rock properties as stated in Sect. 2.2. This model scenario with a spatial variability in rock properties provides a single quantification of the resulting spatial stress variability. However, the rock properties are assigned to the elements in an arbitrary way. Thus, this model scenario is only one possibility amongst many ( Fig. 2a, b). No exhaustive information on the range of other possibilities, i.e. the uncertainties of the stress state at a given location due to the variabilities in the rock properties is provided yet. An estimation of uncertainties requires a full scan of all possible distributions of the rock properties variability. This can be achieved by numerous model scenarios with the same model geometry and the same calibration point, but different rock properties from the same probability distribution. It is automatized using Python scripts that automatically assign the rock properties (Ziegler 2021), calibrate the desired stress state (Ziegler and Heidbach 2021), and evaluate the effect of the randomly assigned variability. Once a sufficient amount of different model scenarios have been evaluated the script is halted and the uncertainties of all components of the stress tensor are displayed.
In order to estimate the sufficient amount of model scenarios, a break condition is defined. It is tied to the tapering off of changes in the stress states standard deviation from all model scenarios that have been computed up to this point. The criterion is tested after each model scenario run by an estimation of the standard deviation of each stress tensor component in this and all previous model scenario runs by with the average stress where n is the number of model runs that have been computed. If the changes in the standard deviations compared to the ten previous condition checks n − 1, n − 2, … n − 10 are < 0.01 MPa for each of the six independent components of the stress tensor, the break condition is met and the script is halted.
It is then assumed that even with additional model scenarios no significant changes in the stress states standard deviation or in the average stress state will occur due to the randomly assigned rock properties. This implicitly assumes that the parameter space has been scanned sufficiently, i.e. the rock properties have been evenly distributed throughout the model volume in the sum of the model scenarios. The uncertainties in the stress state due to the rock properties variability can now be quantified. In order to prevent a premature break of the script this procedure is not limited to one location but repeated at 150 locations throughout the model volume.

Analysis Procedure
The result of the model workflow and analysis is the standard deviations of the six independent stress tensor components throughout the model volume. They are estimated from all model scenarios once the break criterion of the script was met. The rock properties for each element of the models only lithological unit are drawn from the same distribution of rock properties and thus a smooth average stress gradient is achieved (Fig. 2c) and the standard deviation is comparable throughout the model volume (Fig. 2d). The only deviation is observed at the location used for calibration of the model, where the rock properties are fixed on their average value in all model scenarios (Fig. 2d).
In order to test different probability distributions of rock properties, the described procedure is repeated in a different batch of model runs. Within each batch the calibrated stress state and the probability distributions of the rock properties are the same. The only variation between model scenarios is the randomly performed assignment of rock properties based on the specified probability distributions. Different batches, however, allow to change the calibrated stress state and the probability distributions of the rock properties in order to investigate their effect on the uncertainties of the stress state.
The bottom row of Fig. 2 only shows the resulting average stress state (Fig. 2c) and the associated uncertainties in the stress field (standard deviation, Fig. 2d) due to one set of rock properties and their variability of 1 = 30% in all three rock properties. Also, only two components of the stress tensor ( xx and zz ) are visualized. However, the aim is to characterize all six components of the stress tensor and their uncertainties due to differently large variability in the rock properties. Thus, the observed standard deviations from within each sub-volume are presented in a histogram ( Fig. 3a) for different batches with different variabilities of the rock properties. To allow for a better comparability, the same information is also displayed in as boxplot (Fig. 3b) which will be exclusively used in the later figures.

Results
In order to investigate the influence of the variability in the rock properties on the stress states uncertainties, different batches with differently large variabilities in the three rock properties in only one of them are tested while the others are fixed. Furthermore, the influence of different stress regime (i.e. changes of the stress magnitude ratios in the sub-volume with the a priori known stress state and fixed rock properties), and the influence of different rock types-different mean rock properties-on the uncertainties of the stress field are investigated. zz (horizontal plane) from two model scenarios with different randomly assigned rock properties (a, b). Bottom: Average stress (c, left) and standard deviation (d, right) in the xx (vertical profile) and zz components of the stress tensor in a model setup with a single unit. The sub-volume where the stress state is a priori known and where rock properties are fixed is clearly indicated on the vertical profile by the small standard deviation (d)

Influence of the Three Rock Properties
The resulting uncertainties in Fig. 3 are due to an amalgamation of the effects of the variabilities in the three rock properties Young's modulus, Poisson's ratio, and Density. To investigate the effect of the variability in each individual rock property, only one of the rock properties is assigned a variability while the others are fixed on their mean value. The results are shown in the three panels of Fig. 4, each panel for variabilities in one rock property.
As already indicated in Fig. 3, mainly affected by the variability in rock properties are the normal components of the stress tensor and in particular xx and yy (Figs. 3 and 4). The separate investigation of variabilities in only one of the Fig. 3 a Colour-coded standard deviations of the modelled stress states components (x-axis) throughout the model volume and its variability (y-axis) as a histogram. Different variabilities between 10% (back) and 50% (front) of all three rock properties is assumed. b The same data is displayed with the variability (x-axis) and the standard deviation in the modelled stress state (y-axis) in a box plot. The bold line indicates the range of 50% of the data which contains the median (dot). The thin line represents the whiskers, the range of values that are no outliers but outside the interquartile range. The maximum length of the whiskers is the 150% of the interquartile range. Individual outlier data points are indicated by a circle a b c Young's modulus Density P oisson's ratio three rock properties confirms that the uncertainties in the shear components xy , yz , and zx are generally small. The absolute magnitude of shear stresses in this model setup are expected to be small. However, only the variability of the standard deviations are interpreted which reduces the impact of absolute stress magnitudes. Nevertheless, a potential bias in the variability of the shear components due to the small absolute stress magnitudes is investigated in Sect. 3.2. It is indicated in Fig. 4a that the main contributor to uncertainties in the stress state is the Young's modulus. The uncertainties are mainly focused on the xx and yy components of the stress tensor. The zz component is next with a significantly higher uncertainty than the shear components but also clearly below the other two normal components. Conversely, variabilities in the density mainly affect the zz component of the stress tensor and only secondarily xx and yy components with the shear components again clearly negligible (Fig. 4b). In the upper kilometres, this can be easily explained by the overburden as the source for zz which depends on the density. The influence of variabilities in the Poisson's ratio is only significant for the xx and yy components while all other components including zz are negligible (Fig. 4c). In all three panels of Fig. 4 the one sub-volume that is used for calibration and has accordingly no variability is clearly indicated by single outliers in xx and yy .
For variabilities (up to 40% to 50%) the relationship between the variability in the rock properties and the uncertainties in the stress state follow a linear trend. For larger variabilities, this linear trend is flattening (Fig. 4). In particular for the Poisson's ratio (Fig. 4c) this is clearly visible. This bend in the curve is caused by the absence of negative rock properties. With increasing variability, the drawing of negative values from the distribution becomes exceedingly probable (Fig. 5a). To prevent this, in such cases the negative value is discarded and a new value is drawn until a positive value is found (Ziegler 2021). This imposes a biased distribution of assigned rock properties. The somewhat more pronounced bend in the Poisson's ratio curve (Fig. 4c) is due to the double limit of the Poisson's ratio. Neither negative values nor values larger than 0.5 are allowed (Fig. 5b).

Influence of the Stress Magnitudes
The stress state itself is variable in terms of the absolute magnitudes and the differential stresses. It is investigated whether stress magnitudes affect the stress states perceptibility to rock property variability in a way that the same rock property variability results in different uncertainties of the modelled stress state. This influence is investigated by comparing the uncertainties in the three depths of − 600 m, − 2300 m, and − 4000 m with accordingly different stress magnitudes while maintaining the same rock property variability.
Most prominently it is observed that uncertainties in zz are positively correlated with an increasing stress magnitude (Fig. 6). However, this is not uniquely due to the larger stress magnitude but rather due to the increase in the overburdens variability. The range of possible overburden increases with depth due to the variability in the density which adds up with every overlying sub-volume. Additional differences in uncertainties due to different stress magnitudes are observed but have only a limited significance compared to the three different rock properties themselves. The dominance of the uncertainties in xx over yy is established and becomes more pronounced with an increase in stress magnitude/depth.
Even though the absolute uncertainties in the shear stresses are negligible, it is interesting to note that xy is positively correlated with the absolute stress magnitudes/depth while yz and zx are negatively correlated. Still, the observed differences in how the the shear components are affected in different depths (i.e. by different absolute stress magnitudes) is small. This is in particular small in comparison to the normal components variabilities and also compared to how zz is affected. This indicates that the uncertainties in the shear components are not biased by their small absolute magnitudes.

Influence of the Stress Regime
In addition to the previously addressed different stress magnitudes, the orientation of the principal stress components is important for characterization of the stress state and define the three tectonic stress regimes normal, strike slip, and thrust faulting (Anderson 1905 (Ziegler and Heidbach 2021). The analysis of the results is focussed on the respective depths in order to assure the comparability of all stress magnitudes (including zz ) and to exclude any influence of the stress magnitudes. Figure 7 shows the results in a 3-by-3 matrix that juxtaposes the three different rock properties with the three stress regimes. The uncertainties in the modelled stress state due to variabilities in the Young's modulus show a clear dependency on the stress regime (Fig. 7, left column). The differences are chiefly concentrated on the horizontal normal components xx and yy which also have the largest uncertainties in general. In a normal faulting stress regime uncertainties in yy are primarily affected by the Young's modulus variability. Concurrently, in a thrust faulting regime uncertainties in xx are primarily affected by the Young's modulus variability.
Variabilities in the density have no significantly different effect on the uncertainties in the modelled stress state in the different tectonic stress regimes (Fig. 7, central column). Despite a slightly larger uncertainty in zz in normal faulting stress regimes, followed by strike slip and finally thrust faulting no changes are observed. These deviations are probably due to the decrease in absolute zz magnitude from normal via strike slip to thrust faulting but are so small that they are not significant.
Variability in the Poisson's ratios affects the uncertainties in the modelled stress state in a way that is clearly dependent on the tectonic stress regime (Fig. 7, right column). Contrary to the behaviour of the Young's modulus, in a normal faulting regime the uncertainties in xx are primarily affected and in a thrust faulting regime uncertainties in yy . Both stress tensor components are clearly less and almost equally affected by variabilities in the Poisson's ratio in a strike slip regime.

Absolute Rock Properties
Initially, this study aims at representation of average rock types with regards to their mean properties. However, they are by no means representative for the entire diversity of existing rocks. Thus, the influence of the mean rock property value on the uncertainties in the stress state is evaluated, i.e. different variabilities around different mean rock properties are tested. Young's moduli between 10 and 90 GPa (in steps of 20 GPa), Densities between 2100 and 3300 kg/m 3 (in steps of 300 kg/m 3 ), and Poisson's ratios between 0.20 and 0.36 (in steps of 0.04) are evaluated. Previously, the variability was altered based on percentage of the mean value. However, with changes in the mean value the definition of Even though the different mean rock properties increases or decreases the resulting modelled uncertainties, Fig. 8 illustrates that in contrast to changes in the tectonic stress regime, no relative changes between the uncertainties in the different components of the stress tensor occur. Considering variabilities in the Young's modulus, the modelled uncertainties in xx are always larger than those in yy , no matter the absolute rock properties values.
For a fixed variability (Fig. 8 top row) of the Young's modulus (± 15 GPa) and the Poisson's ratio (± 0.084) higher mean rock property result in a decrease of the uncertainties (Fig. 8a, c). In particular in case of the Young's modulus this decrease is dramatic while for the Poisson's ratio it is almost negligible. Independent of the average rock property, no changes are observed in case of the density (Fig. 8b). Apparently, the effect of the increase in average rock property is mitigated by the relatively decreasing variability which is constant at 810 kg/m 3 .
If a variability of ± 30% relative to the average rock property is assigned to the Young's modulus, changes in the absolute value of the rock property do not affect the uncertainties in the stress state (Fig. 8d). This is due to the  variabilities increase together with the average value of the rock property. The relative value of 1 remains identical and is thus reflected in the stress states uncertainties that remains constant. Conversely, for an increase in the density an increase in the stress states uncertainties is observed with a relative variability of 30% (Fig. 8e). The same is observed for the Poisson's ratio (Fig. 8f).
The Young's modulus variabilities influence on the uncertainties in the modelled stress state is in agreement with the intuitive expectations. The densities variabilities influence on the uncertainties of the modelled stress state, however, behaves counterintuitive. This is because the density is affected by an increase in absolute variability with depth, i.e. the range of possible total overburden increases (see Sect. 3 usually non-linear even without considering heterogeneities due to different rock units different densities (Warpinski and Teufel 1991). The overburden is responsible for the vertical normal stress according to with the density , the depth z, and the gravitational acceleration g. This formula indicates the possibility to estimate the uncertainties solely introduced by the density on the modelled stress state by an analytical approach as well. Furthermore, the other components of the stress state are affected by a variability in the density via the Poisson's ratio according to with the affected stress tensor lm , the unaffected stress tensor ij , and ij as the Kronecker delta. This explains the contrasting behaviour compared to the Young's modulus in all components of the stress tensor. The Poisson's ratio behaves peculiar in that the average values are negatively correlated with the uncertainties in the stress field when assuming a constant variability and positively correlated when assuming a relative variability. Even though the Poisson's ratio itself is not depth-dependent, the relationship with the overburden, described in the previous equations, makes the impact of its effects depth-dependent (as the density in Fig. 8e). Concurrently, the non-depthdependent ties to the Young's modulus lead to the similar behaviour (Fig. 8d).

General Discussion
Numerical modelling of the in-situ stress state is commonly performed for engineering, exploratory, or seismic hazard assessment purposes by means of 3D geomechanicalnumerical models (Goertz-Allmann and Wiemer 2013; Ziegler et al. 2015;Yoon et al. 2015;Jeanne et al. 2016;Blöcher et al. 2018). The uncertainties of such models are indicated to be large (Hergert et al. 2015;Ziegler et al. 2016) and recently approaches and methods to quantify these uncertainties became available (Lecampion and Lei 2010;Ziegler and Heidbach 2020). However, these approaches usually concern the uncertainties due to stress magnitude data records. Thus, the influence of uncertainties in the 3D geological model (Wellmann et al. 2010;Bond et al. 2015) is disregarded by assumption of perfectly known strata. In addition, perfectly known (no measurement error) and homogeneous (one set of rock properties for an entire unit) rock properties are assumed. However, rock properties are one of the key drivers of uncertainties (Hergert et al. 2015;Ziegler et al. 2016). Their influence on the stress field is indicated by heterogeneities found by extensive stress estimation campaigns (Warpinski and Teufel 1991;Zhao et al. 2015;Feng et al. 2020) and theoretically by means of modelling (Roche and van der Baan 2017; Lei and Gao 2019;Fraldi et al. 2019).
Herein, only the influence of rock property variability on the modelled stress state is regarded and quantified. This is done by means of a one-unit generic hexahedron model that does not include any geological boundaries, faults, or other structures that are known to somehow affect the stress state. A more complex geometry is likely to have an additional impact on the variability of the stress state due to rock properties. Furthermore, for simplicity, the expected uncertainties due to measurement errors in the stress data records are disregarded by assuming a perfectly known stress state for calibration. Furthermore it is assumed that the rock properties in the REV used for stress state calibration are perfectly know. The assessment of the influence of aforementioned additional drivers for aleatoric uncertainties as well as any epistemic uncertainties are outside the scope of this work. Nevertheless, they should be addressed in future research.
This study is focused on heterogeneous rock properties and their influence on mechanical models of the in-situ stress state. Such models only require definition of the three regarded rock properties Young's modulus, Poisson's ratio, and Density. The properties are drawn from representative distributions which do not include all possible outliers of rock types with extreme properties. Such materials should be investigated separately which may even require a different rheological behaviour. The chosen properties are assumed to be a representation of the rock mass behaviour within a model element. It is implicitly assumed that the heterogeneity in the three properties is also caused by variabilities in additional parameters such as rock strength or porosity. With respect to the element size significant disturbances of the rock (e.g. wide faults, natural caves, mined out areas) should be discretized and not parameterised. Kinematic modelling approaches (Hergert et al. 2011) or fracture propagation models (Yoon et al. 2015) will face additional variabilities that need to be explicitly regarded. This also holds for Hydro-Mechanical or Thermo-Hydro-Mechanical models which require a range of additional property information such as permeability or thermal conductivity (Altmann et al. 2010;Ziegler et al. 2017).

Impact of Two-layer Geometries with Different Rock Property Variability
The results are scale invariant due to the linearity of the problem formulation which allows a broader application to real world settings. However, contrary to the previously presented generic approach, usually stress models contain more than just one lithological unit. Different units accordingly have different rock properties and different associated variabilities. The influence of the rock properties variability in one unit on the modelled stress states uncertainties in another unit is tested. A two-unit setting is investigated in an advanced approach which has an extent of 10 × 10 × 5 km 3 (Fig. 9a) and is discretised with 3000 elements that are assigned to 375 sub-volumes. At a depth of 1500 m a horizontal boundary separates the sedimentary Top unit from the crystalline Bottom unit. For each of the two units, the rock properties for each individual element are drawn from the according probability distributions around the average rock properties (Fig. 9a). An estimation of the uncertainties in the stress state due to the initially specified variability is displayed in Fig. 9b which indicates the smaller uncertainties in the stress state of the upper unit due to the smaller variability in the rock properties.
In order to further test the behaviour, the variability in the upper unit is varied while the one in the bottom unit is This effect that even a reduction of the rock properties' variability only in one unit leads to a reduction of the stress states uncertainties in all units is beneficial for geomechanical modelling that aims at predicting the in-situ stress state. Often the rock properties knowledge for in-situ stress models is limited to a restricted area of interest. Well-known rock properties in this area hence potentially positively affect the uncertainties in the entire model. In addition, a decrease in variability in a certain unit, e.g. due to an extensive measurement campaign that reduces the measurement error, is far more feasible than a decrease of uncertainties in all present units. In particular, this holds for shallow units which are feasibly accessed from an economic and technical point of view.

Applications
This study is designed to address the effects of rock property uncertainties on the modelled 3D stress state in order to enable quantification and reduction of model uncertainties and improve interpretation of modelled stress tensors. A case study and application of the findings on a particular engineering project should be subject of future research. This requires detailed data on the according rock properties distributions. The intuitive approach is to apply the herein presented method to a full-scale 3D geomechanical-numerical model of a particular project. However, for large models this is currently not possible due to the required computation time that exceed any feasible amount of time by orders of magnitude. Thus, such an exhaustive approach to quantify the uncertainties is not an option.
Instead, it is proposed to establish an ad-hoc quantification of uncertainties. The relationship between rock property variabilities and the modelled stress state are investigated using generic models such as the presented. The results are used to generate lookup tables that contain the uncertainties in the modelled stress state due to different variability distributions in the rock properties in different settings. Therefore, sufficient combinations of rock properties and variabilities need to be assessed to set up a significant lookup table. The table is then used for post-computation assignment of uncertainties to the modelled stress state.
The presented results indicate important factors that influence a lookup table. The rock properties and their variability, the stress magnitudes and the tectonic stress regime have an impact on the models perceptibility to rock property variability. Accordingly, they need to be addressed in a lookup table. Further research is required to investigate whether there is a significant impact of the geometry, e.g. pinching or cropping out units, salt domes. In particular complex geometries other than horizontal layering may have an effect but an investigation is beyond the scope of this work. Furthermore, the applicability and limits of such a lookup table approach need to be addressed in future research.
Any application needs to address the possibility that site-specific or rock-specific correlation and/or dependencies between the rock properties exist. Their investigation requires robust information on rock properties from comprehensive measurement campaigns or databases such as Bär et al. (2020). If a significant correlation or dependency can be observed, the variabilities in all three rock properties can be reduced, i.e. the variability in the modelled stress state due to all three rock properties (Fig. 3) are reduced while the variability in the modelled stress state due to each individual rock property remains the same (Fig. 4). However, in terms of correlation or dependencies, an approach that does not consider them cannot underestimate the variability in the stress state. Any correlation or dependency would necessarily reduce the variability in the modelled stress state.

Conclusion
The presented generic model workflow quantifies systematically the impact of rock property variability on the uncertainties of the six independent components of the 3D stress tensor. In particular the model result for the one-layer problem show the following: • The variability in the Young's modulus has the largest effect on the uncertainties, followed by the density and the Poisson's ratio. • The normal components are significantly stronger affected by a factor of 2-6 in comparison to the shear components. • The tectonic stress regime changes the results, but the general findings are still valid. • These general findings are independent from the number of layers that are implemented. If the variability of rock properties in only one unit is reduced, a reduction of uncertainties is also observed in the other units stress state.
However, the application of this approach to real settings is in most cases not feasible due to extraordinary long computation time required. To mitigate this while also enabling a significant assessment of associated uncertainties, it is proposed to use the previously found relations between the variability in rock properties and the modelled stress states' uncertainties for post-computation assignment of uncertainties. This lookup table approach will significantly reduce the cost but further research is required to evaluate the applicability, provide rules for implementation, and generate accordingly refined lookup tables.

Appendix: Numerical Implementation
Using the finite element method requires continuity of the field variable at the nodes of the elements. The field variable in this case is the displacement since the problem of the equilibrium of forces is a second order partial differential equation formulated in displacements. The stress tensor is estimated at the integration points that are located at predefined positions within the elements. However, to visualize the stress state the values are extrapolated from the integration points to the nodes and thus, the stresses at the nodes are a mean of the stresses estimated at the nodes that are shared by several elements. Thus, if elements that share a node have different rock properties the effect of these is reduced. This effect is stronger in 1st order finite elements where the stress is constant within the elements and significantly lower in 2nd order elements since stress varies linearly in these. The generic model described herein is constructed entirely of brick elements and thus the following description refers to elements of the according geometry. For 8-nodeelements with linear basis functions (1st order elements) at each node seven other elements border (Fig. 11a). In case of 20-node elements with quadratic basis functions (2nd order elements) there are between three and seven adjacent elements (Fig. 11b). High contrasts in the rock properties between neighbouring elements are thus smoothed in case of using 1st and 2nd order elements.
To mitigate this smoothing effect, the assignment of rock properties is not done element-wise but to groups of eight elements each (Fig. 11c). Such a group features seven nodes whose stress values are independent from other element groups with different rock properties. In each direction in space three of these consecutive independent integration points are available (Fig. 11c). At these nodes the full impact of variabilities in the rock properties can be observed. Hence, only the stress state from these central nodes is used for displaying the results. At the borders to other elementboxes a smoothed stress state that does not represent the actually expected uncertainties is observed (Fig. 11d) due to the amalgamation of effective rock properties.
Acknowledgements The author thanks two anonymous reviewers whose comments improved the manuscript. The work leading to these results has received funding from the Initiative and Networking Fund of the Helmholtz Association through the project "Integrity of nuclear waste repository systems -Cross-scale system understanding and analysis (iCross)" project by the Federal Ministry of Education and Research (project number 02NUK053D), the Helmholtz Association (project number SO-093) and the Helmholtz Centre Potsdam -Deutsches GeoForschungsZentrum GFZ.
Funding Open Access funding enabled and organized by Projekt DEAL.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.