Elastic and Frictional Properties of Fault Zones in Reservoir-Scale Hydro-Mechanical Models—A Sensitivity Study

: The proper representation of faults in coupled hydro-mechanical reservoir models is challenged, among others, by the di ﬀ erence between the small-scale heterogeneity of fault zones observed in nature and the large size of the calculation cells in numerical simulations. In the present study we use a generic ﬁnite element (FE) model with a volumetric fault zone description to examine what e ﬀ ect the corresponding upscaled material parameters have on pore pressures, stresses, and deformation within and surrounding the fault zone. Such a sensitivity study is important as the usually poor data base regarding speciﬁc hydro-mechanical fault properties as well as the upscaling process introduces uncertainties, whose impact on the modelling results is otherwise di ﬃ cult to assess. Altogether, 87 scenarios with di ﬀ erent elastic and plastic parameter combinations were studied. Numerical modelling results indicate that Young’s modulus and cohesion assigned to the fault zone have the strongest inﬂuence on the stress and strain perturbations, both in absolute numbers as well as regarding the spatial extent. Angle of internal friction has only a minor and Poisson’s ratio of the fault zone a negligible impact. Finally, some general recommendations concerning the choice of mechanical fault zone properties for reservoir-scale hydro-mechanical models are given. scenarios studied and the corresponding model names. The scenarios di ﬀ er regarding the elastic and the frictional properties assigned to the model unit representing the fault zone.


Introduction
Hydro-mechanical simulations have developed into a standard tool for various subsurface applications ranging from hydrocarbon and geothermal reservoirs to underground storage sites for CO 2 [1][2][3][4]. However, challenges exist regarding the proper implementation of faults into such numerical models. Faults not only have a profound impact on fluid flow, but also effect the stress field in their vicinity. In addition, pore pressure changes due to injection or production can induce slippage and fault reactivation, respectively [5][6][7][8]. This may cause induced seismicity, land subsidence and well collapse [9][10][11][12]. Fault reactivation may also breach the reservoir seal causing up-fault leakage and allowing fluid migration due to enhanced permeability inside the fault zone [13][14][15]. Thus, proper incorporation of faults into hydro-mechanical models is of crucial relevance for various reasons.
In reality, faults are characterized by a complex geometrical structure and a heterogeneous material distribution [16][17][18][19][20][21]. They should be regarded rather as a fault zone, i.e., a volumetric feature, than as a discrete fault surface [20,[22][23][24]. A fault zone typically contains a fault core accompanied by damage zones on either side, but may be structurally more complex including the appearance of multiple fault cores as well as damage zones of variable width and fracture intensity [16,21,[25][26][27]. Therefore, the appropriate representation of faults in numerical reservoir models is challenged by two related aspects: (1) the difference in scale between the heterogeneity of the fault zone (centimeters to meters) and the typical size of the calculation cells of the numerical grid (meters to tens of meters) and (2) the material properties assigned to the fault zone which stem-if available at all-from rock mechanical testing on core samples with a diameter of a few centimeters. Both aspects require upscaling, i.e., merging the heterogeneous material properties of the fault zone [28][29][30][31], in order to assign reasonable material properties to the fault zone elements of a reservoir-(kilometer)-scale hydro-mechanical model.
In the present study we use a generic finite element (FE) model with a volumetric fault zone description to examine what effect such upscaled parameters have on pore pressures, stresses, and deformation within and surrounding the fault zone. Such a sensitivity study is important as the usually poor data base regarding specific hydro-mechanical fault properties as well as the upscaling process introduces uncertainties, whose impact on the modeling results is otherwise difficult to assess. As the prime focus is on fault reactivation due to pore pressure changes, we explore a parameter range for Young's modulus, Poisson's ratio, cohesion, and angle of internal friction which is typical for fault zone rocks. Modeling results are expected to provide a framework to understand quantitatively how the elastic and frictional-plastic properties assigned to a fault zone in a hydro-mechanical reservoir simulation actually affect the modeling results.

Elastic and Frictional Fault Zone Properties
Faults are usually not a target for drilling operations. They are rather avoided to reduce well stability problems. Even if coring of a fault zone is attempted, intense fracturing and poor consolidation of fault zone rocks frequently results in very limited core recovery [32,33]. As a consequence, there usually is a lack of fault-specific material parameters to populate a hydro-mechanical model [34,35]. Instead, the material parameters have to be estimated from literature sources considering how faulting and fracturing change the material properties of the intact rock. In the following we review the possible range of mechanical fault rock parameters in the elastic and frictional-plastic domain, which also forms the basis for the subsequent numerical sensitivity study.

Elastic Material Properties
Deformation in the elastic domain can be described by Hooke's Law, which relates stress and strain via an elasticity (or stiffness) matrix [36,37]. In case of a linear-elastic and isotropic medium, only two material properties are required: Young's modulus and Poisson's ratio. The Young's modulus or modulus of elasticity of a fault zone can vary by several orders of magnitude between the value of the intact rock and the one of the fault-related rocks in the damage zones and fault core. The Young's modulus for common lithologies (intact rock) ranges from 1 GPa to 100 GPa [24,38]. It decreases with increasing fracture density, i.e., in the damage zone and tectonic breccias, respectively. Fault gauges which frequently make up the fault core can exhibit Young's modulus values as low as 0.01 GPa [24,39,40].
The Young's modulus required for numerical modelling is the static Young's modulus which-ideally-is derived from uniaxial or triaxial testing of cores in the laboratory. Such rock mechanical testing is typically carried out on samples with a diameter of a few centimeters. However, mechanical properties like Young's modulus and unconfined compressive strength (UCS) derived from lab tests are usually significantly larger (1.5×-10×) than the corresponding values of a larger rock mass [24,37]. This scale-dependence can be accounted for by using empirical correlations, but this inevitably adds uncertainties to the material properties assigned to the fault zone in the numerical model [38,41,42]. If no cores are available, an option is to calculate dynamic Young's modulus from well logs (p-and s-wave velocities, density). However, these dynamic Young's modulus values need to be converted to static ones based on empirical correlations which introduces yet another uncertainty to the upscaled properties.

•
The higher the deformation and the strain values during fault zone formation, the more likely the rock mass is highly deformed and powdered. Thus, cohesion decreases with increasing deformation during fault formation [15,48,51,55].

•
Hydrothermal circulation (depending on the chemical content of the circulating fluids) can overprint and cement the fault zone rocks leading to increased cohesion [15,26,56].
It should be noted that for sedimentary rocks the intercept of a linear Mohr-Coulomb criterion actually overestimates rock cement strength. In such cases, use of a nonlinear failure envelope would be more appropriate to describe real cohesion at low effective confining stresses.
The second frictional material parameter used in the MC criterion is the angle of internal friction. It describes the increase in shear strength of a rock with respect to an increase in normal stress [15,24,39]. Friction angles of common lithologies range between 5 • and 50 • [24,37]. A frequently used average value for rocks is 40.5 • equivalent to a coefficient of internal friction of 0.85 [50]. In general, the friction angle is related to grain size, so coarse-grained rocks tend to have a higher friction angle, while lower friction angles are typical for fine-grained rocks and rocks with a high clay content, respectively [57][58][59][60][61]. Besides grain size, the degree of fracturing also affects the angle of internal friction [62]. Typically, the friction angle decreases with increasing fracturing of the rock [63,64].
However, other studies [65,66] indicate that during rock degradation in shear and without further weakening due to mineral alteration and clay smearing, the friction angle remains largely unaltered. Thus, weakening in such cases results primarily from cohesion softening.
A fault-specific determination of cohesion and angle of internal friction again requires core samples for triaxial rock mechanical testing in the lab. If the material is not suitable for multistage testing [67][68][69], at least three samples are required for one triaxial test. Thus, it is very rare that proper rock samples for triaxial testing of fault zone rocks are available to populate a numerical model with specific frictional fault rock properties [15,24].

Modelling Concept
A finite element (FE) model was set up to compare the effect on hydro-mechanical simulation of different elastic and frictional-plastic properties assigned to a fault zone. It consists of over and underburden sections of low permeability with a high-permeability reservoir layer in between and offset by a normal fault. Figure 1a exhibits an elliptical fault zone displacing a reservoir horizon, the geological rationale for this model setup. The reservoir horizon embedded into the 1 km × 1 km × 0.003 km model frame is displaced by a 500 m long normal fault centered in the middle of the model. At the top of the model frame, loading conditions simulate 1 km of overlaying rock material. So the model bottom is adopted to be 2 km beneath the surface. Additionally, the weight of the overburden implemented as pressure acting on the model top, no further mechanical forces, or displacements are assigned to the model. On the hydraulic part, the initial model assumes a hydrostatic pore pressure field. Subsequently, pore pressure at the nodes on the right side of the reservoir unit is increased to simulate fluid injection into the downthrown block. For five years the pore pressure is constantly increased every three months by 0.75 MPa until it finally reaches 15 MPa above hydrostatic.

Modelling Concept
A finite element (FE) model was set up to compare the effect on hydro-mechanical simulation of different elastic and frictional-plastic properties assigned to a fault zone. It consists of over and underburden sections of low permeability with a high-permeability reservoir layer in between and offset by a normal fault. Figure 1a exhibits an elliptical fault zone displacing a reservoir horizon, the geological rationale for this model setup. The reservoir horizon embedded into the 1 km × 1 km × 0.003 km model frame is displaced by a 500 m long normal fault centered in the middle of the model. At the top of the model frame, loading conditions simulate 1 km of overlaying rock material. So the model bottom is adopted to be 2 km beneath the surface. Additionally, the weight of the overburden implemented as pressure acting on the model top, no further mechanical forces, or displacements are assigned to the model. On the hydraulic part, the initial model assumes a hydrostatic pore pressure field. Subsequently, pore pressure at the nodes on the right side of the reservoir unit is increased to simulate fluid injection into the downthrown block. For five years the pore pressure is constantly increased every three months by 0.75 MPa until it finally reaches 15 MPa above hydrostatic.  Starting from a base model which is used as reference for the subsequent sensitivity studies, Young's modulus, Poisson's ratio, cohesion, and angle of internal friction are varied within a range typical for fault zone rocks. At the same time, the model geometry and the material properties for the other model units (reservoir rock and over-/underburden) remain the same in all scenarios. Altogether, 87 FE models with different mechanical parameter combinations assigned to the fault zone have been studied. FE software Ansys 19.2 [70] is used for the fully coupled hydro-mechanical simulations.

Model Geometry
The 3D model represents a slice through the central part of a normal fault (Figure 1a). The resulting model domain, the dimensions as well as the initial and boundary conditions are shown in Figure 1b. The model dimensions are 1 km × 1 km × 0.001 km (equivalent to a model slice of one element thickness) and comprise a fault zone of 500 m height and 60 • dip. The fault offsets a 25 m thick reservoir section by 47 m. Using scaling relationships of Johri et al. [33], the total width of the fault zone including the fault core as well the damage zone is assumed to be 12 m. Treffeisen and Henk [71] have shown that in reservoir-scale models discretizing the width of a fault zone into three elements is sufficient to achieve modelling results very similar to a finer discretization. Hence, the size of the elements describing the fault geometry is 4 m × 4 m × 3 m, whereas the whole model domain comprises 45,684 elements.

Constitutive Laws
The hydro-mechanical simulation uses fluid flow through a porous medium and poroelastic-perfectly plastic material behavior. In the following the principle behind such fully coupled hydro-mechanical simulations as well as the assigned boundary conditions are briefly explained. The Modelling utilizes the commercial software program Ansys 19.2 and the reader is referred to the extensive theory and user manuals [70] for more detailed information.
The mechanical part of the fully coupled simulation is based on the total stresses are reduced by the pore pressure in the rock volume leading to effective stresses. This relationship can be described by [72][73][74]: where the effective stresses (σ ij ) are derived from the total stress tensor (σ ij ) by subtracting the pore pressure (p f ), weighted by the Biot coefficient (α) and Kronecker's delta (δ ij ). The hydraulic part of the coupling is defined in Ansys 19.2 with the following mass-balance equation [70]: where the flux due to Darcy's law q f (m/s), the volumetric strain ε v (-), the compressibility parameter Q* and the degree of fluid saturation S f (-). Q* is calculated from the bulk modulus of both the skeleton and the fluid. These two formulas briefly outline the interaction of the hydraulic and mechanical behavior as changes in pore pressure result in effective stress changes and related volumetric strain. This in turn alters porosity and permeability, which again affect the pore pressure field [1,37]. Figure 1b shows the mechanical and hydraulic boundary conditions, both for the initial state and the subsequent injection stage. For the mechanical calculations, displacements orthogonal to the model boundaries are not allowed at the base, left, right, front, and back sides of the model (roller boundary conditions). Since the model top is located 1000 m beneath the earth surface, the weight of the overburden is assigned through a pressure boundary condition with an average density of 2300 kg/m 3 . For the hydraulic calculations, all model boundaries are considered to be impermeable.

Initial and Boundary Conditions
During the initial load step, mechanical and hydraulic equilibrium in response to the boundary conditions is achieved. So the whole model domain experiences hydrostatic pore pressure conditions. Concurrently the initial stress field is established through gravitational loading only and no other tectonic (horizontal) stress components are considered in the simulations. Therefore, the tectonic regime is normal faulting with the vertical stress S v being the maximum principal stress S 1 .
Afterwards, fluid injection in the downthrown block is simulated by successively applying higher pore pressures for five years to the boundary nodes of the reservoir layer in the hanging wall of the fault (lower right in Figure 1b). The injection pressure increases for 20 time steps by 0.75 MPa, until the maximum injection pressure of 15 MPa above hydrostatic is reached.

Material Parameters
The hydraulic and mechanical material properties assigned to the fault zone, reservoir and host rock are listed in Table 1. Details of the changes in Young's modulus, Poisson's ratio, cohesion, and friction angle assigned to the fault zone in the various sensitivity studies as well as the acronyms dedicated to the corresponding scenarios are shown in Table 2. The coding used for each model variant reflects the differences to the base model (BM). For example, BM-c1 implies that the model varies only regarding cohesion of the fault zone rocks and that a value of 1 MPa (instead of 5 MPa as in model BM) is used for the calculations. All other material parameters remain the same.   The material properties assigned to the fault zone elements cover a range which is typical for fault zones in siliciclastic successions, i.e., Young's modulus was varied between 0.01 GPa and 10 GPa, Poisson's ratio between 0.1 and 0.25, cohesion between 0.01 MPa and 5 MPa and friction angle between 10.0 • and 25.0 • . For example, the combination of the maximum values from each category, e.g., a Young's modulus of 10 GPa, a Poisson's ratio of 0.25, a cohesion of 5 MPa and a friction angle of 25 • , would be typical for a fault zone containing slightly deformed sandstone and re-cemented fault gouge [15,24,39]. Combination of the minimum values, i.e., a Young's modulus of 0.01 GPa, a Poisson's ratio of 0.1, a cohesion of 0.01 MPa and a friction angle of 10 • would represent a fault rock containing mostly clay and siliciclastic rock flour [26]. Between these extremes, various other rock types and, hence, material properties can be imagined, depending on the fraction clay or rock and the fraction of lithified sandstone [43].
The permeability of fault zones in siliciclastic successions is usually controlled by the clay content of the offset rocks. Thus, impermeable faults are more likely in sandstone-shale regimes [16,59,75]. However, we decided to choose a semipermeable fault zone since we are particularly interested in how fluid flow effects the mechanical strength of the fault and how this potentially leads to fault reactivation [2,15]. Faults in sandstones with less than 15% clay content show already a slightly reduced permeability in comparison to the host rock [76]. For this sensitivity study, we choose a fault zone permeability about two orders of magnitude lower than the reservoir rock to account for some fault damage on hydraulic conductivity. However, modelling techniques are generally applicable and could also be applied to fault zones which act as impermeable barriers or high-permeability conduits for fluid flow.

Results
The following section presents some of the modeling results. Firstly, the results of the base model (BM) are presented, which provides the baseline for the following comparisons with some of the other 86 scenarios studied (see Table 2 for the various models). Afterwards stress and strain patterns for both the whole model area and the fault zone itself are compared between exemplary model setups.

Base Model
The base model (BM) comprises a fault zone with material properties typical for siliciclastic fault zones separating a sandstone reservoir, whereby the fault acts as conduit for fluid flow between different reservoir compartments. The corresponding mechanical material properties for a fault zone in such a setting are inferred from literature [24,40,43]   Results of the base model (BM) of the spatial variation in pore pressure (a,d), the magnitude of the effective maximum principal stress (S1,eff; b,e), and orientation of S1,eff (deviation from vertical; c,f).

Magnitude and Orientation of S1,eff
The effective maximum principal stress (S1,eff) displays 15 MPa at the top and up to 37 MPa at the bottom of the model (Figure 2b). This reflects the combined effects of boundary conditions and gravitational loading and considers poroelastic effects. The fault zone shows the lowest effective stress values due to its higher Biot coefficient in comparison to the surrounding rock mass. Figure 2c displays the orientation of S1,eff. For ease of visualization, this is presented in form of a contour plot which show the deviations (clockwise (CW) and counterclockwise (CCW)) from the vertical direction, i.e., the undisturbed orientation of S1,eff in a normal faulting regime. As a consequence, outside the model areas affected by the fault hardly any deviation occurs, i.e., S1,eff is vertical. However, inside the fault zone S1,eff has rotated up to 17° CW. In the immediate vicinity of the fault zone and near the fault tips, stress rotations between 5° CCW and 5° CW can be observed.
After five years of fluid injection the pore pressure inside the reservoir ( Figure 2d) and the fault zone has increased which directly effects the S1,eff pattern ( Figure 2e). The lower, right reservoir section stands out against the surrounding rock mass with S1,eff values about 5-10 MPa lower, while the effective stresses near the fault zone show an increase at the upper left and lower right fault tips as well as a corresponding decrease at the opposing sides. Inside the fault zone, S1,eff has decreased up to less than 10 MPa. Injection does not only modify effective stress magnitudes but also stress orientations ( Figure 2f). Significantly larger model areas up to 300 m distance from the fault zone are now showing S1,eff orientations differing from vertical. These stress perturbations follow some kind of a point symmetry with a maximum rotation of 21° CW to the left of the lower and to the right of the upper fault tip, respectively. On the opposite sides of the fault tips extending to the middle part of the fault zone, the direction of S1,eff deviates between 5° and 15° CW from vertical. Figure 3 shows the von Mises total strain as well as the corresponding elastic and plastic strain components. For the first load step, total strain ranges from 0 to 0.16% ( Figure 3a) and is dominated by elastic straining (compare to Figure 3b,c). After five years of injection, the total strain exhibits values between 0.0 and 1.6% ( Figure 3d) and is primarily controlled by plastic straining, which shows

Magnitude and Orientation of S 1,eff
The effective maximum principal stress (S 1,eff ) displays 15 MPa at the top and up to 37 MPa at the bottom of the model (Figure 2b). This reflects the combined effects of boundary conditions and gravitational loading and considers poroelastic effects. The fault zone shows the lowest effective stress values due to its higher Biot coefficient in comparison to the surrounding rock mass. Figure 2c displays the orientation of S 1,eff . For ease of visualization, this is presented in form of a contour plot which show the deviations (clockwise (CW) and counterclockwise (CCW)) from the vertical direction, i.e., the undisturbed orientation of S 1,eff in a normal faulting regime. As a consequence, outside the model areas affected by the fault hardly any deviation occurs, i.e., S 1,eff is vertical. However, inside the fault zone S 1,eff has rotated up to 17 • CW. In the immediate vicinity of the fault zone and near the fault tips, stress rotations between 5 • CCW and 5 • CW can be observed.
After five years of fluid injection the pore pressure inside the reservoir ( Figure 2d) and the fault zone has increased which directly effects the S 1,eff pattern (  Figure 3d) and is primarily controlled by plastic straining, which shows significantly higher values than the elastic strain component (compare to Figure 3e,f). Elastic straining only in relation to the initial stress field is shown in Figure 3b. In the vicinity of the fault zone von Mises elastic strain values range from 0.04% at the top to 0.1% at base. Somewhat larger strain values are observed for the weaker fault zone where values range between 0.12% and 0.16%. These lower elastic strains are confined to the fault zone and do not extend into the surrounding rock mass. No von Mises plastic strain is observed for the initial load step (Figure 3c). Thus, it can be excluded that initial loading already causes plastic deformation in and reactivation of the fault zone. For the final load step, i.e., after five years of injection, elastic strain in the fault zone has decreased by about 0.04% compared to the initial load step, ranging now between 0.08% and 0.12% (Figure 3e). The larger elastic strain values extend into the surrounding host rock at the fault tips. While elastic strain has decreased inside the fault zone, plastic strain has increased and occurs throughout the entire fault zone (Figure 3f). Maximum values are about 1.6% at the fault tips, but even within the fault zone plastic straining of at least 0.5% can be observed.
Energies 2020, 13, x FOR PEER REVIEW 10 of 28 significantly higher values than the elastic strain component (compare to Figure 3e,f). Elastic straining only in relation to the initial stress field is shown in Figure 3b. In the vicinity of the fault zone von Mises elastic strain values range from 0.04% at the top to 0.1% at base. Somewhat larger strain values are observed for the weaker fault zone where values range between 0.12% and 0.16%. These lower elastic strains are confined to the fault zone and do not extend into the surrounding rock mass. No von Mises plastic strain is observed for the initial load step (Figure 3c). Thus, it can be excluded that initial loading already causes plastic deformation in and reactivation of the fault zone. For the final load step, i.e., after five years of injection, elastic strain in the fault zone has decreased by about 0.04% compared to the initial load step, ranging now between 0.08% and 0.12% (Figure 3e). The larger elastic strain values extend into the surrounding host rock at the fault tips. While elastic strain has decreased inside the fault zone, plastic strain has increased and occurs throughout the entire fault zone ( Figure  3f). Maximum values are about 1.6% at the fault tips, but even within the fault zone plastic straining of at least 0.5% can be observed.

Elastic and Plastic Strain
The impact of different fault zone Young's modulus values on elastic and plastic strain patterns after five years of injection is shown in Figure 5. Model BM shows a maximum elastic strain of 0.12% in the fault zone elements (Figure 5a; see also Figure 3e which uses a different color scale to visualize details). Decreasing Young's modulus by a factor of 10 (BM-Y1), the maximum values for the elastic strain inside the fault zone is increased to 1.03% (Figure 5b). A further decrease in Young's modulus The results for the plastic strain appear to be the other way around (Figure 5d-f). Model BM seems to show plastic straining primarily in the lower part of the fault zone (Figure 5d). However, Figure 3f, which uses a different color contour scheme, shows, that plastic straining actually occurs throughout the entire fault zone with a peak value of 1.5% at the lower fault tip. In model BM-Y1, no through-going plastic strain zone has developed, but the strain peak at the lower fault tip is about twice as much with 3.06% (Figure 5e). For BM-Y0.1 the plastic strain is limited to only a very small area at the lower fault tip, where the maximum value is 1.02% (Figure 5f). The remaining part of the fault zone does not experience any plastic straining. The results for the plastic strain appear to be the other way around (Figure 5d-f). Model BM seems to show plastic straining primarily in the lower part of the fault zone (Figure 5d). However, Figure 3f, which uses a different color contour scheme, shows, that plastic straining actually occurs throughout the entire fault zone with a peak value of 1.5% at the lower fault tip. In model BM-Y1, no through-going plastic strain zone has developed, but the strain peak at the lower fault tip is about twice as much with 3.06% (Figure 5e). For BM-Y0.1 the plastic strain is limited to only a very small area at the lower fault tip, where the maximum value is 1.02% (Figure 5f). The remaining part of the fault zone does not experience any plastic straining.  The effect of different Poisson's ratios of the fault zone rocks on the magnitude of both the vertical (S1,eff) and horizontal (S2,eff) effective stress is shown in Figure 6. The pattern of S1,eff magnitudes after five years of fluid injection is essentially the same for all scenarios i.e., for Poisson's ratios of 0.25 (BM, Figure 6a

Magnitude and Orientation of S1 ,eff and S 2,eff
The effect of different Poisson's ratios of the fault zone rocks on the magnitude of both the vertical (S 1,eff ) and horizontal (S 2,eff ) effective stress is shown in Figure 6. The pattern of S 1,eff magnitudes after five years of fluid injection is essentially the same for all scenarios i.e., for Poisson's ratios of 0.25 (BM, Figure 6a  Regarding the orientation of S1,eff and S2,eff, the calculated patterns exhibit no substantial differences depending of the Poisson's ratio assigned to the fault zone unit.  Regarding the orientation of S 1,eff and S 2,eff , the calculated patterns exhibit no substantial differences depending of the Poisson's ratio assigned to the fault zone unit.      The corresponding plastic strain pattern after five years of injection is presented in Figure 9d-f. Due to the color contour scheme required to compare the various simulation results, model BM seems to show plastic strain only in a small area near the base of the fault zone (Figure 9d). However, the detailed analysis shown in Figure 3f reveals that plastic straining actually occurs in the entire fault zone. Plastic strain varies between 0.31% and 0.62% throughout the fault zone and reaches a maximum of 1.54% at the lower fault tip. For a decrease in cohesion to 0.5 MPa (BM-c0.5), the values inside the fault zone increase to 0.81% to 1.61% with a local maximum at the lower fault tip of 4.03% (Figure 9e). A further decrease in cohesion to 0.01 MPa (BM-c0.01) increases the maximum peak to 4.43%, but within the fault zone itself, the values only increase to between 0.89% and 1.77%. However, the basic plastic strain pattern remains the same (Figure 9f). The corresponding plastic strain pattern after five years of injection is presented in Figure 9d-f. Due to the color contour scheme required to compare the various simulation results, model BM seems to show plastic strain only in a small area near the base of the fault zone (Figure 9d). However, the detailed analysis shown in Figure 3f reveals that plastic straining actually occurs in the entire fault zone. Plastic strain varies between 0.31% and 0.62% throughout the fault zone and reaches a maximum of 1.54% at the lower fault tip. For a decrease in cohesion to 0.5 MPa (BM-c0.5), the values inside the fault zone increase to 0.81% to 1.61% with a local maximum at the lower fault tip of 4.03% (Figure 9e). A further decrease in cohesion to 0.01 MPa (BM-c0.01) increases the maximum peak to 4.43%, but within the fault zone itself, the values only increase to between 0.89% and 1.77%. However, the basic plastic strain pattern remains the same (Figure 9f).     Regarding the orientation of S1,eff, Figure 10d

Elastic and Plastic Strain
Regarding elastic strain, modeling shows rather similar results for the three different friction angles. There are only minor differences in the strain distribution within the fault zone and the maximum values reached. Model BM model exhibits values up to 0.09%-0.1% in large parts of the fault zone, whereas the peak value of 0.13% is confined to a small area near the lower fault tip ( Figure  11a). For model BM-FA17.5, the elastic strain observed inside the fault zone decreases to 0.08%-0.09% (Figure 11b), but the area affected by higher elastic strains in continuation of the fault both updip and downdip is enlarged. This trend continues for a friction angle to 10.0° (BM-FA10) for which elastic strain within the fault zone is reduced to 0.07%-0.08% (Figure 11c), while the area affected by elastic straining in the vicinity of the fault zone widens.

Elastic and Plastic Strain
Regarding elastic strain, modeling shows rather similar results for the three different friction angles. There are only minor differences in the strain distribution within the fault zone and the maximum values reached. Model BM model exhibits values up to 0.09%-0.1% in large parts of the fault zone, whereas the peak value of 0.13% is confined to a small area near the lower fault tip (Figure 11a). For model BM-FA17.5, the elastic strain observed inside the fault zone decreases to 0.08%-0.09% (Figure 11b), but the area affected by higher elastic strains in continuation of the fault both updip and downdip is enlarged. This trend continues for a friction angle to 10.0 • (BM-FA10) for which elastic strain within the fault zone is reduced to 0.07%-0.08% (Figure 11c), while the area affected by elastic straining in the vicinity of the fault zone widens.
Considering plastic strain, the trend is reverse. Instead of a decrease in fault zone straining with a decrease in friction angle as was the case for elastic strain, plastic strain within the fault zone increases with decreasing angle of internal friction of the fault rock. Model BM shows plastic strain values of less than 0.4% in most parts of the fault zone ( Figure 11d) and a maximum of 1.54% at the lower fault tip. BM-FA17.5 exhibits plastic strain within the fault zone of up to 0.8% (Figure 11e) and to 2.07% at the lower fault tip. Both values increase even further in BM-FA10 to 1.0% plastic strain within the fault zone and 2.61% (Figure 11f). The impact of different elastic fault zone material properties on the simulation results after five years of fluid injection are examined by analyzing twelve scenarios with different Young's modulus and Poisson's ratios, respectively. Figure 12 shows the maximum values calculated for each model regarding the magnitude of S1,eff, elastic and plastic strain plotted against the Young's modulus for different Poisson's ratios. All three plots indicate that the different curves, each representing a specific Poisson's ratio, are very close, especially for the magnitude of S1,eff and the elastic strain (Figure 12a,  b). Only for plastic strain they differ between 1.0% and 1.75% for a Young's modulus of 0.1 GPa.  The impact of different elastic fault zone material properties on the simulation results after five years of fluid injection are examined by analyzing twelve scenarios with different Young's modulus and Poisson's ratios, respectively. Figure 12 shows the maximum values calculated for each model regarding the magnitude of S 1,eff , elastic and plastic strain plotted against the Young's modulus for different Poisson's ratios. All three plots indicate that the different curves, each representing a specific Poisson's ratio, are very close, especially for the magnitude of S 1,eff and the elastic strain (Figure 12a,b). Only for plastic strain they differ between 1.0% and 1.75% for a Young's modulus of 0.1 GPa.
Energies 2020, 13, x FOR PEER REVIEW 19 of 28 A large influence on the modelling results can be observed for Young's modulus. Considering Figure 12a, the magnitude of S1,eff decreases by 20 to 25 MPa if the Young's modulus increases about one order of magnitude. In addition, elastic strain shows decreasing maximum values for increasing Young's moduli (Figure 12b). However, plastic strain results are different as illustrated in Figure 12 c. It shows the maximum elastic, plastic, and total strain values for different Young's moduli. At first, plastic strain increases with increasing Young's modulus until a maximum value of about 3.4% is A large influence on the modelling results can be observed for Young's modulus. Considering Figure 12a, the magnitude of S 1,eff decreases by 20 to 25 MPa if the Young's modulus increases about one order of magnitude. In addition, elastic strain shows decreasing maximum values for increasing Young's moduli (Figure 12b). However, plastic strain results are different as illustrated in Figure 12 c. It shows the maximum elastic, plastic, and total strain values for different Young's moduli. At first, plastic strain increases with increasing Young's modulus until a maximum value of about 3.4% is reached for a Young's modulus of 0.5 GPa. From that point on plastic strain decreases with further increasing fault zone Young's modulus until maximum plastic straining is 1.75% for a Young's modulus of 10 GPa (Figure 12c). At the same time, elastic straining reaches its minimum. However, even though the maximum values at the fault tips differ, the general pattern within the fault zone proper remains the same: increasing Young's modulus lead to larger plastic strain observed in the fault zone.  (Figure 13 b), where the values are 4.43% (BM-c0.01), 4.03% (BM-c0.5) and 3.74% (BM-c1) before decreasing to 1.54% (BM). In both graphs of Figure 13, the distances between the curves derived for various internal friction angles are the same. Thus, for a constant cohesion value, the maximum S 1,eff magnitude and the maximum plastic strain decrease linearly with an increase in friction angle. Although the maximum values partially differ significantly, the general modeling results are rather similar for variations in cohesion and internal friction angle. In Figure 14 the stress pattern for the magnitude of S1,eff is shown as an example and similar trends can be observed in orientation plots as well. The base model with its rock properties is displayed in Figure 14a. The rows display the different results achieved for a constant cohesion value by varying the internal friction angle from left (25°) to right (10°) and the columns show the results for one specific friction angle while varying the cohesion from top (5 MPa) to bottom (0.01 MPa). An increase in both the peak magnitude of S1,eff as well as the spatial extent of the stress perturbation can be observed for a decrease in friction angle Although the maximum values partially differ significantly, the general modeling results are rather similar for variations in cohesion and internal friction angle. In Figure 14 the stress pattern for the magnitude of S 1,eff is shown as an example and similar trends can be observed in orientation plots as well. The base model with its rock properties is displayed in Figure 14a. The rows display the different results achieved for a constant cohesion value by varying the internal friction angle from left (25 • ) to right (10 • ) and the columns show the results for one specific friction angle while varying the cohesion from top (5 MPa) to bottom (0.01 MPa). An increase in both the peak magnitude of S 1,eff as well as the spatial extent of the stress perturbation can be observed for a decrease in friction angle as well as a decrease in cohesion. Apart from the stress peaks at the fault tips, the effect of varying the internal friction angle on the overall stress pattern is not visible at all and the differences can be only observed at the peak values located at the fault tips.

Discussion
Firstly, the base model (BM) is evaluated to ensure proper model behavior and, building on that, the modelling results for different material properties inside the fault zone are discussed.

Base Model
The BM results for pore pressure (Figure 2a) and S1,eff magnitudes outside the areas affected by the fault zone show the expected increase with depth resulting from the different fluid and rock densities (Figure 2b). The S1,eff orientation in these areas is vertical and thus, in agreement with the normal faulting regime [77] induced by the prescribed boundary conditions (Figure 2c). The fault zone is not reactivated as a result of initial loading as no plastic straining (Figure 3c) and no differential displacements are observed. After five years of fluid injection into the lower reservoir section, fluid migration through the fault zone has also led to an increase in pore pressure in the

Discussion
Firstly, the base model (BM) is evaluated to ensure proper model behavior and, building on that, the modelling results for different material properties inside the fault zone are discussed.

Base Model
The BM results for pore pressure (Figure 2a) and S 1,eff magnitudes outside the areas affected by the fault zone show the expected increase with depth resulting from the different fluid and rock densities (Figure 2b). The S 1,eff orientation in these areas is vertical and thus, in agreement with the normal faulting regime [77] induced by the prescribed boundary conditions (Figure 2c). The fault zone is not reactivated as a result of initial loading as no plastic straining (Figure 3c) and no differential displacements are observed. After five years of fluid injection into the lower reservoir section, fluid migration through the fault zone has also led to an increase in pore pressure in the upper reservoir part (Figure 2d). More importantly, the pore pressure increase within the fault zone has reduced the effective stresses (Figure 2e), thus, leading to plastic straining ( Figure 3f) and reactivation of the fault as well as differential displacements along its entire length. This fault reactivation has also modified the stress magnitudes ( Figure 2e) and orientations (Figure 2f) in the vicinity of the fault zone. Comparable patterns regarding spatial extent and magnitude of the stress changes have also been reported from other numerical simulations of faults [5,14,78].

Variations in Elastic Rock Properties
Comparison of the scenarios which differ with respect to Young's modulus shows that a decrease in Young's modulus assigned to the fault zone leads to an increase in the magnitude and re-orientation of S 1,eff as well as the spatial extent of the stress perturbations ( Figure 4). Thus, this relationship follows a negative correlation. For example, the models evaluated in Section 4.2.1 show for BM, BM-Y1, and BM-Y0.1 that the magnitude and the rotation angle of S 1,eff increases if the Young's modulus of the fault zone is decreased. These trends for Young's modulus can also be observed if combined with other material properties like different Poisson's ratios, cohesions, and internal friction angles. Thus, stiffness contrasts between the fault represented as volumetric weak zone and the host rock lead to magnitude changes as well as stress rotations. Softer rock properties in the fault zone lead to a larger stiffness contrast with respect to the surrounding, which results in broader ranges of between maxima and minima for S 1,eff magnitudes and stress reorientations, respectively [79,80]. Similarly, the spatial extent of the stress perturbations increases with decreasing stiffness of the fault zone and increasing stiffness contrast, respectively.
The stiffness contrast between the fault zone and the surrounding rock mass also affects the strain values observed in the fault zone. A larger stiffness contrast and therefore a broader ranges of S 1,eff magnitudes results in larger values of total strain (Figure 12c). With increasing Young's modulus in the fault zone and, hence, a reduction in stiffness contrast to the country rock, the total strain decreases.
The different Young's moduli imply different slopes of the curves in the stress-strain diagram. Thus, for a given failure stress lower Young's moduli result in larger elastic deformation prior to failure or, in other words, for a given amount of total deformation the share of plastic deformation is less in case of lower Young's modulus. These trends are supported by the modelling results, i.e., peak S 1,eff magnitude, elastic strain as well as total strain decrease with increasing Young's modulus of the fault zone rocks. Likewise, the plastic strain should increase with increasing Young's modulus, which is the case for Young's modulus between 0.1 and 0.5 GPa. However, the maximum value is observed for a Young's modulus of 0.5 GPa and the maximum plastic strain decreases between a Young's modulus of 0.5 and 10 GPa (Figure 12c). Thus, the total strain decrease in response to a Young's modulus increase delimits the maximum plastic strain which is achieved in response to respective stress states. The elastic strain component reaches its maximum at a low fault zone stiffness, whereas at high fault zone stiffnesses total strain is dominated by plastic straining.
Contrary to the large influence Young's modulus has on the modelling results, the Poisson's ratio assigned to the fault zone appears to have a negligible impact. This is somewhat surprising as under gravitational loading only Poisson's ratio relates the horizontal stress S H (MPa) to the vertical stress S V (MPa) according to [37]: Although this relationship holds for the rock units surrounding the fault zone, changing Poisson ratio ν (-) inside the fault zone (while keeping it constant in the rest of the model) has hardly any impact onto the modelling results. The effect of different Poisson's ratios inside the fault zone on the stress and strain pattern is apparently suppressed by the influence the Poisson's ratio of the surrounding rock mass has on the overall stress and strain patterns.

Variations in Plastic Rock Properties
If elastic fault zone properties are kept constant and only frictional-plastic material parameters are varied, modeling results are rather similar as long as cohesion is less than 1 MPa (Figures 8 and 9). For larger cohesion values, a substantial drop in both the magnitude of S 1,eff and plastic strain is observed ( Figure 13). The increase in pore pressure and the corresponding decrease in effective stresses within the fault zone shifts the Mohr circle in the Mohr-Coulomb diagram towards the failure line and results, in case of lower cohesion values, in larger plastic straining. In turn, larger plastic straining and fault reactivation leads to a stronger loading of the vicinity of the fault and, hence, larger S 1,eff magnitudes. For the modeling scenario selected in this study with gravitational loading only, a cohesion of 1 MPa appears to be some kind of a threshold value at larger cohesion values substantially less fault reactivation occurs resulting in reduced stress perturbations concerning both magnitude and spatial extent.
If cohesion is kept constant, the internal friction angle shows a negative correlation for all result parameters. Thus, an increase in the internal friction angle decreases the magnitude and rotation angle of the effective maximum principal stress (S 1,eff ) as well as the values for maximum elastic and plastic strain. This, again, can be explained by the Mohr-Coulomb diagram: the larger the friction angle and, hence, the steeper the slope of the failure line is, the more differential stress the fault zone rock can stand prior to failure. Thus, more stress states are in the plastic domain if the friction angle is lower which in turn results in stronger stress perturbations again.
Differences in the overall pattern for both stress and strain can hardly be detected and the results only seem to differ at the fault tips ( Figure 10, Figure 11, and Figure 13). In reality however a fault usually does not end within the reservoir and therefore the area of interest in reservoir modelling, hence the effects of that only occur within a few meters at the fault tips and are not important in reservoir engineering.

Conclusions
Using a simple generic model setup of a normal fault offsetting a reservoir horizon, the impact of different mechanical fault properties on the stress and strain distribution within a fault zone and its surrounding is analyzed. Thereby, the fault is represented as a volumetric weak zone with uniform properties which inherently implies upscaling of the small-scale heterogeneity of faults observed in nature as well as of the limited size of lab samples available for mechanical testing to the large cell size of the reservoir-scale hydro-mechanical model. This study provides insights on how the elastic and plastic material properties assigned to fault zones of a reservoir-scale model actually affect the results of the numerical simulation. Therefore, it helps to choose proper upscaled fault zone properties by knowing the effect each mechanical parameter has.
Starting with a base model as reference, the mechanical properties assigned to the fault zone are varied within a range typical for faults in sandstone-shale successions, i.e., Young's modulus between 0.1 and 10 GPa, Poisson's ratio between 0.10 and 0.25, cohesion between 0.01 and 5 MPa and angle of internal friction between 10 • and 25 • . Altogether, 87 different scenarios were studied.
Modelling results indicate that these four material properties have substantially different influences on the stress and strain perturbations induced by the fault. The Young's modulus of the fault zone and, more specifically, the contrast in Young's modulus between fault zone and surrounding rocks has by far the strongest impact onto the modelling result. The larger the difference in Young's modulus is, the larger the stress and strain perturbations are, both in absolute numbers as well as regarding the spatial extent. Cohesion turns out to be the second most important material property. In comparison to these two parameters, the impact of the internal friction angle is minor and of the Poisson's ratio assigned to the fault zone essentially negligible.
Consequently, the following guidelines for modelling of faulted reservoirs can be derived from this sensitivity study: If the stress and strain patterns in the vicinity of faults, i.e., within a few hundred meters, are not in the focus of the particular modelling project, the fault zone properties can be determined based on literature data and/or rules of thumb as outlined in Chapter 2. This holds because modelling results indicate that the exact fault properties have no impact on the stress and strain distribution in the far-field of the fault. This, however, changes in the vicinity of a reservoir-scale fault zone, i.e., within a few hundred meters distance to the fault.
If this near-field of a fault is of interest for the modeling study, the reservoir simulation workflow has to start with a first guess of the fault zone parameters based on literature data and experience as usually only very limited or no fault-specific material parameters are available. The initial model run should then be followed by an iterative calibration process by comparing the simulation results to data actually observed. These could be stress orientations (e.g., from borehole breakouts and drilling-induced fractures observed in image logs) and/or stress magnitudes (e.g., from hydraulic tests) observed in wells in the vicinity of the fault. According to the results of the present sensitivity study, during the calibration process the modeler can focus on Young's modulus and cohesion as the most important parameters to achieve a satisfactory fit between model predictions and reality. These two properties are the crucial parameters to assess the reactivation potential of a fault as well as the stress and strain patterns in its vicinity. The angle of internal friction is of less relevance and for Poisson's ratio the standard value of 0.25 can be adopted without concern. This reduction in the number of variables to be tuned during calibration substantially reduces time and computational costs to achieve a validated fault zone description for the hydro-mechanical reservoir model. This sensitivity studies cannot provide quantitative relationships between the parameters since too many variables besides the four material properties studied influence fault response in hydro-mechanical models. For example, different elastic properties in over and underburden relative to the fault zone would lead to variable stiffness contrast in different sections of the fault zone. This, in turn, would affect the stress and strain patterns in the vicinity of the fault.
However, this sensitivity study provides a guideline which material properties should be of prime interest during upscaling and model calibration. This is also valid, if local mesh refinement in the reservoir-scale model is performed in order to capture the geometrical heterogeneity of the fault zone in greater detail and with smaller elements, respectively. Perspectives for future modelling work include highly detailed fault zone models to provide a sound physical base for upscaling of hydro-mechanical material parameters.