Faults as Volumetric Weak Zones in Reservoir-Scale Hydro-Mechanical Finite Element Models—A Comparison Based on Grid Geometry, Mesh Resolution and Fault Dip

: An appropriate representation of faults is fundamental for hydro-mechanical reservoir models to obtain robust quantitative insights into the spatial distribution of stress, strain and pore pressure. Using a generic model containing a reservoir layer displaced by a fault, we examine three issues which are typically encountered if faults have to be incorporated in reservoir-scale ﬁnite element simulations. These are (1) mesh resolution aspects honoring the scale di ﬀ erence between the typical cell size of the ﬁnite element (FE) reservoir model and the heterogeneity of a fault zone, (2) grid geometry relative to the fault geometry and (3) fault dip. Di ﬀ erent fault representations were implemented and compared regarding those on the modeling results. Remarkable di ﬀ erences in the calculated stress and strain patterns as well as the pore pressure ﬁeld are observed. The modeling results are used to infer some general recommendations concerning the implementation of faults in hydro-mechanical reservoir models regarding mesh resolution and grid geometry, taking into account model-scale and scope of interest. The goal is to gain more realistic simulations and, hence, more reliable results regarding fault representation in reservoir models to improve production, lower cost and reduce risk during subsurface operations.


Introduction
Since the beginning of the 21st century, the hydrocarbon industry has shifted towards faulted and structurally more complex conventional or unconventional reservoirs, which require a thorough understanding of both the hydraulic and the mechanical reservoir behavior [1][2][3]. Fluid flow simulations are a well-established tool in the industry [4,5] and geomechanical modeling has also turned out to be of tremendous help when trying to gain quantitative insights into the spatial distribution of stress and strain on the reservoir-scale [6][7][8][9][10]. Due to the interaction of fluid flow and mechanical behavior, fully-coupled hydro-mechanical simulations gain more and more importance.
Various numerical modeling techniques have been tried and tested, e.g., finite difference (FD), boundary element (BE), discrete element (DE) and hybrid methods [11][12][13][14][15], but the most commonly used approach for hydro-mechanical reservoir simulations is the finite element (FE) method which the present study focuses on. Such FE reservoir models typically have a lateral size between kilometers to tens of kilometers and applications can range from hydrocarbon and geothermal reservoirs to sites for underground gas storage [16][17][18]. In order to obtain a realistic subsurface representation for reliable stress and fracture predictions as well as fluid flow path analysis, the reservoir models have to take into accounts faults, i.e., discontinuities offsetting the strata, in addition to the lithostratigraphic horizons.
However, faults are not only characterized by local discrete deformation. Faults can also induce stress perturbations, i.e. local changes in stress orientation and stress magnitude, which differ distinctly from the regional trends [19]. Hydraulically, faults can either act as barriers or conduits for fluid flow between different reservoir compartments [20,21]. Considering the usually weaker mechanical properties of fault zones, if compared to the undeformed rock, they are more sensitive to pore pressure changes resulting from fluid injection into or fluid withdrawal from the reservoir. Those pore pressure changes can reactivate preexisting faults, leading to induced seismicity potentially causing critical situations like, fault seal breach, land subsidence and well collapse [6,21,22]. Consequently, a realistic representation of faults in a hydro-mechanical reservoir simulation proves to be invaluable for the reliability of the numerical model predictions. This has been an issue in various studies over the past years [8,10,21,23,24]. Several authors suggest that in order to honor the architecture and internal heterogeneity of faults, they should be treated as volumetric features in the numerical simulations, even if the fault zone thickness seems to be negligible in contrast to the overall model and element size, respectively [25,26]. This is considered of particular relevance for complex, intensively faulted reservoirs [4,27] and could improve projects where faults are part of the full-scale reservoir model e.g., [25,28,29].
In this paper, we analyze the representation of fault zones as volumetric weak zones in hydro-mechanical simulations regarding three commonly used assumptions in FE reservoir modeling. These are (1) mesh resolution aspects honoring the scale difference between the typical cell size of the FE reservoir model and the heterogeneity of a fault zone, (2) grid geometry relative to the fault geometry and (3) fault dip. Different combinations of these three parameters were investigated and compared regarding their impact on the modeling results. The primary goal is to provide guidelines for appropriate fault representations in numerical reservoir models regarding mesh resolution, grid geometry and scope of interest. Such guidelines will assist in building more realistic hydro-mechanical simulations of faulted reservoirs. Furthermore, the modeling techniques developed and experiences gained are also considered as a starting point for further studies to investigate the effect of fault zone heterogeneities and to develop refined upscaling techniques for hydro-mechanical fault zone properties.
General information about the state of the art of fault modeling in FE reservoir models can be found in Section 2 for both fluid flow and geomechanical simulations as well as coupled hydro-mechanical simulations. Section 3 describes the three commonly used assumptions to incorporate faults in reservoir models which we investigate in this study. The model setup as well as the constitutive laws for the simulation can be found in Section 4. Section 5 provides the modeling results, which are discussed in Section 6 and lead to the conclusions given in Section 7.

Faults in Fluid Flow Simulations
For fluid flow simulations, the use of transmissibility multipliers is widely used in reservoir engineering [30][31][32][33]. Transmissibility multipliers are assigned to the cells of the calculation grid intersected by the fault plane. Depending on the selected hydraulic properties which are averaged over the entire fault zone and adjacent rock, these cells then act as barriers or conduits for fluid flow. Different concepts have been suggested over the past years to estimate fault transmissibility multipliers [34][35][36]. More recently, a fault facies concept was developed to improve fault representation in reservoir-scale fluid flow simulations by representing faults as 3D structures with variable material properties for different parts of the fault zone. The fault facies concept tries to include the internal complexity of fault zones and their permeability structure by assigning various (hydraulic) material properties to the rock bodies influenced by faulting. It is assumed that the fault volume can be populated with a fault facies similar to how geological models are populated with a sedimentary facies concept [3,5,20,24,25,37]. However, the basic challenge remains how to properly upscale the complex internal structure and the heterogeneous permeability distribution which exist in real fault zones.

Faults in Geomechanical Simulations
In addition to hydraulic properties, fault zones are also highly heterogeneous with respect to mechanical properties. This heterogeneity together with the complex fault geometry and internal architecture will strongly control the mechanical behavior of a fault zone and how it affects the stress field in its vicinity [38][39][40]. There are different approaches to implement faults in geomechanical models, but for FE modeling, which is the most commonly used approach in hydro-mechanical reservoir simulations, there are: (1) Contact or interface elements [11,41,42] and (2) volumetric weak zones [43][44][45].
Contact or interface elements represent faults as discrete discontinuities [46,47], which is an expansion of the classical continuum approach used in FE simulations. Those elements enable differential movements between separate and individually meshed model parts. A contact surface consists of contact and their correspondent target elements located at opposing sides of the predefined fault. However, a contact-and target-element pair has the same coordinates. Hence, the fault is represented as a line in 2D and a surface in 3D models, respectively [9,47], instead of a volumetric feature.
Infinitely high normal stiffness values would be needed to enforce compatibility between adjacent fault surfaces and to create a zero thickness of the contact elements. In order to prevent numerical instabilities, a certain mesh penetration at the interface has to be accepted. Thus, stiffness values close to Young's moduli of the adjacent rocks are commonly used for the contact elements. Shear and normal stress can be transmitted through the contact elements and material properties like cohesion and a friction coefficient can be used to describe frictional sliding. Therefore, once the stress state at the contact elements violates the defined failure criterion [42], relative displacement between corresponding nodal points will occur.
In contrast, the concept of volumetric weak zones describes faults in a FE model by assigning (usually weaker) fault rock material properties to the grid cells intersected by faults. This concept can be applied to a rectangular grid, resulting in a stair-stepped fault representation (Figure 1 Left) [48,49]. Alternatively, the FE mesh can be adapted to the fault zone geometry, resulting in a curvilinear grid (Figure 1 Right) [12,44,50]. Both grid geometries are further explained in Section 3.2.
Energies 2020, 13, 2673 4 of 28 Figure 1. Images (not to scale) showing the basic differences between the different grid geometries applied. Left: The rectangular approach (R) uses a regular, rectangular grid in which fault zone properties are assigned to those cells intersected by the fault plane. Due to the underlying grid geometry, a staircase-shaped fault representation results. Right: The curvilinear approach (C) uses an irregular grid geometry which is adapted to the fault geometry. This allows the representation of the true shape of the fault zone in the FE model.

Faults in Coupled Hydro-Mechanical (HM) Simulations
Hydraulic and mechanical processes are closely related [19,21,38]. Hence, pore pressure changes can lead to mechanical reactivation of a fault zone and fault reactivation can lead to increasing permeability, fluid flow and fault seal break [1,53]. Coupled hydro-mechanical simulations can reproduce those relationships. It is the volumetric representations of faults that can map both weaker mechanical properties and fault-specific hydraulic parameters onto fault-intersecting elements [1,54].

Some General Aspects of Representing Faults in Finite Element Reservoir Models
The goal of this study is to show the sensitivity of the fault representation as a volumetric weak zone regarding three common assumptions used to incorporate fault zones into coupled hydromechanical reservoir simulations. showing the basic differences between the different grid geometries applied. Left: The rectangular approach (R) uses a regular, rectangular grid in which fault zone properties are assigned to those cells intersected by the fault plane. Due to the underlying grid geometry, a staircase-shaped fault representation results. Right: The curvilinear approach (C) uses an irregular grid geometry which is adapted to the fault geometry. This allows the representation of the true shape of the fault zone in the FE model. For reservoir models with a typical cell size distinctly larger than the thickness of the fault zone and its subunits, no details of the fault architecture and material heterogeneity can hardly be considered explicitly. Instead, mechanical properties, which represent the joint effects of host rock, damage zone and fault core, i.e., following a homogenized continuum approach, are used. Optionally, anisotropic properties can be assigned to the fault cells [4,51]. However, at least to a certain extent, local grid refinement techniques or detailed submodeling of specific faults in reservoir models allow for incorporation of further details in the fault zone architecture and mechanical property distribution [3,52].

Faults in Coupled Hydro-Mechanical (HM) Simulations
Hydraulic and mechanical processes are closely related [19,21,38]. Hence, pore pressure changes can lead to mechanical reactivation of a fault zone and fault reactivation can lead to increasing permeability, fluid flow and fault seal break [1,53]. Coupled hydro-mechanical simulations can reproduce those relationships. It is the volumetric representations of faults that can map both weaker mechanical properties and fault-specific hydraulic parameters onto fault-intersecting elements [1,54].

Some General Aspects of Representing Faults in Finite Element Reservoir Models
The goal of this study is to show the sensitivity of the fault representation as a volumetric weak zone regarding three common assumptions used to incorporate fault zones into coupled hydro-mechanical reservoir simulations.

Mesh Resolution
FE analysis is affected by both the mesh quality and element quantity [55,56], so meshing is a key issue in FE modeling. While it is established that a finer grid resolution leads to better modeling results, the grid size used in FE reservoir models is often several orders of magnitude larger than the internal architecture of a fault zone and its material heterogeneities [57,58]. Next to upscaling techniques for both hydraulic and mechanical properties, which include the internal features and material heterogeneity of the fault zone [4], detailed rendering of the fault zone architecture through high-resolution local grid refinements around fault zones can lead to simulation results closer to reality [27]. However, more elements result in an exponential increase in the computing time for the simulations [55,56]. Since longer run-times of simulations mean higher costs, it is of importance to find a proper mesh resolution and simultaneously minimize the computational costs. Consequently, a thorough understanding of the influence the mesh resolution has on numerical fault simulations is crucial.

Grid Geometry
We apply two different basic approaches used in previous studies to implement faults as volumetric weak zones in FE hydro-mechanical models ( Figure 1). These approaches differ regarding the basic grid geometry and whether the grid is adapted to the fault geometry or not.

Rectangular Approach
Rectangular grid geometry is a frequently used approach for the representation of faults in FE reservoir models (Figure 1 Left) [48,49]. Such a rectangular element grid often is already available from preceding flow simulations or property modeling. A volumetric weak zone is generated by intersecting this grid with a structural model, e.g., a fault interpretation derived from 3D seismic interpretation and assigning fault rock properties to those cells cut by the fault planes [48,49,59,60].
The fault rock properties assigned to the corresponding grid cells normally use a homogenized continuum, optionally with anisotropic properties, to incorporate the joint effects of fault core and fractured damage zone [4,51] since the cell size is typically much larger than the material heterogeneities inside a fault zone. Neighboring cells maintain the intact rock properties. Such a Energies 2020, 13, 2673 5 of 27 rectangular grid is fast and easy to generate, and therefore, saves time in model generation. However, the staircase-shaped fault geometry can suppress deformation within the fault zone [45].

Curvilinear Approach
Compared to the rectangular grid with a staircase-shaped fault geometry, a mesh adapted to the fault zone geometry allows a representation of the fault which is closer to reality, e.g., a lenticular or elliptical fault shape, which is common for normal faults in sedimentary environments [61,62]. If the FE grid geometry is adjusted to the fault geometry (Figure 1 Right) [43], an irregular, curvilinear grid following the shape of the fault zone results. The approach permits to represent the fault zone as a continuous feature in contrast to the stair-stepped geometry inevitably connected to the rectangular approach. In most cases, the cell size representing the fault zone will be too large to reproduce details of the internal fault zone architecture. Thus, similar to the rectangular grid, the actual heterogeneity of the fault zone has to be upscaled to the size of the elements representing the fault, e.g., by using a homogenized continuum approach (optionally with anisotropic properties). Adjusting the grid to the fault geometry takes significantly more time to generate, especially for large reservoir models including multiple fault zones. Thus, understanding the influence the fault shape has on the modeling results can allow to decide how accurate the shape has to be regarding to the aim of the hydro-mechanical study.

Fault Dip
Faults, besides other subsurface structures, are commonly detected through geophysical methods. Due to the ability of elastic waves to sense a fault in the subsurface, active seismic is the tool mostly used to detect faults [62][63][64][65]. However, there are certain limitations and dependencies, like e.g., the wavelength of the seismic source or the rock density [66,67], influencing the seismic resolution and, hence, the ability to detect details of the fault zone. The vertical seismic resolution is dependent on the wavelength and depth, but for reservoirs, it is usually in the range of tens of meters. Further difficulties exist when trying to detect small structures, such as small-scale faults or particular fault zone features, like fault core, fault throw or damage zone [68]. Identification of steeply inclined fault zones is complex [65,69,70] and although different migration techniques exist [67,71,72], finding the correct fault dip remains difficult for the seismic interpreter. Thus, a frequently used simplification in reservoir modeling is to assume vertical faults, i.e., a fault dip of 90 • . We compare this value with the typical dip angle of 60 • for normal faulting [73]. The aim is to show the impact that improper fault dip estimations can have on the mechanical response of the fault zone.

Model Setup
In order to compare the different options for volumetric fault representations in the FE models outlined above, we use as the basic model setup a high-permeability reservoir layer embedded in low permeability over-/underburden and offset by a normal fault. The various scenarios studied differ regarding (1) the basic grid geometry (rectangular vs. curvilinear), (2) the mesh resolution of the fault zone (1, 3 and 9-elements width) and (3) the fault dip (60 • vs. 90 • ). The material properties for each of the three model units (fault zone, reservoir rock and over-/underburden) remain the same in all scenarios. Table 1 and Figure 2 give an overview of the characteristics of the seven cases studied.
The geological rationale for the model setup, i.e., a lenticular fault zone displacing a reservoir horizon, is presented in Figure 3a. The resulting model domain, the dimensions as well as all initial and boundary conditions are shown Figure 3b. For the fully coupled hydro-mechanical simulations, the FE software Ansys 19.2 is used [74].

Model Geometry
The 3D model represents a slice through the center of a normal fault offsetting a reservoir horizon (Figure 3a). The overall geological setting is a normal faulting regime, including an elliptically shaped fault zone offsetting a reservoir horizon into upper and lower compartments. The model dimensions of the 3D-slice are 1 × 1 km² and 1-element thickness (Figure 3b). The total height

Model Geometry
The 3D model represents a slice through the center of a normal fault offsetting a reservoir horizon ( Figure 3a). The overall geological setting is a normal faulting regime, including an elliptically shaped fault zone offsetting a reservoir horizon into upper and lower compartments. The model dimensions of the 3D-slice are 1 × 1 km 2 and 1-element thickness ( Figure 3b). The total height of the fault is 500 m, whereas the offset of the 25 m thick reservoir horizon is 47 m. Since the fault zone is represented as one material rather than explicitly subdividing it into fault core and damage zone, scaling relationships [75] for both features are combined and the total width of the fault zone is assumed to be 12 m. The fault dip is either 60 • or 90 • , depending on the scenario studied.
Depending on the mesh resolution inside the fault zone (1-, 3-or 9-element width) the whole model has between 126,000 and 1,134,000 elements for the rectangular and between 45,684 and 312,164 elements for the curvilinear approach. The differences arise from the different mesh resolution towards the model boundaries. The fault zone itself is represented by 119 (1-element width), 358 (3-element width) or 3226 (9-element width) elements for the rectangular and 564 (3-element width) or 4004 (9-element width) elements for the curvilinear models. The actual size of the elements describing the fault geometry are 12 × 12 × 9 m (1-element width), 4 × 4 × 3 m (3-element width) and 1.3 × 1.3 × 1 m (9-element width).

Boundary Conditions
The initial as well as the mechanical and hydraulic boundary conditions used for the simulations are shown in Figure 3b. No displacements orthogonal to model boundaries ('roller boundary conditions') are allowed at the basal, left, right, front and back side of the model. A lithostatic pressure boundary condition equivalent to the weight of the overburden holds for the model top. As the top of the model is assumed to be 1000 m beneath the earth's surface, the corresponding pressure acting on the model can be calculated according to [76]: where ρ r is the average density of the overburden rock (assumed to be 2300 kg/m 3 ) [21], g is gravitational acceleration (9.8 m/s 2 ) and d is depth (1000 m). Since no further tectonic (horizontal) stress components are considered in the simulations, the initial stress field results entirely from gravitational loading. Consequently, the model is located in a normal faulting regime, i.e., outside the area affected by the fault the vertical stress S v corresponds to the maximum principal stress S 1 . Hydraulic boundary conditions assume impermeable boundaries at the top and base of the model as well as no horizontal flow through the vertical model sides. Initially, hydrostatic pore pressure conditions are assumed throughout the model. The corresponding pore fluid pressures P f can be calculated according to: where ρ f is the fluid density (1000 kg/m 3 ; assumed average density of the pore fill in the overburden) and d n is the depth of the corresponding element node, i.e., between 1000 m (model top) and 2000 m (model base). After the initial load step, i.e., after mechanical and hydraulic equilibrium in response to the boundary conditions has been achieved, successively higher pore pressures are applied to the boundary nodes of the reservoir layer in the hanging wall of the fault (lower right in Figure 3b). This represents fluid injection in the downthrown block and is continued for 5 years at a rate of 0.75 MPa every 3 months, i.e., in 20 time steps, until the maximum injection pressure of 15 MPa is reached. Therefore, pore pressure (p f ) for the nodes at the boundary of the reservoir layer is increased according to: where p i is the pressure increment (0.75 MPa) and t is the time step (1 to 20).

Constitutive Laws
Poroelasto-plastic material behavior and fluid flow through a porous medium is assumed for this hydro-mechanical simulation. In a saturated porous medium, stresses are split between the solid Energies 2020, 13, 2673 9 of 27 phase and the fluid phase, so the total stresses are reduced by the pore pressure in the rock volume. This relationship can be described by: where the effective stresses (σ ij ) are derived from the total stress tensor (σ ij ) by subtracting the pore pressure (p), weighted by the Biot coefficient (α) and Kronecker's delta (δ ij ) [74][75][76].
The mechanical material behavior in the poroelastic domain is described by the following stress-strain relationship: It links strain (ε ij ) and effective stress (σ ij ) through the material properties of Young's modulus (E) and Poisson's ratio (υ) and the sum of the effective principal stresses (σ' kk ) [77][78][79]. The failure criterion delimiting the poroelastic range and initiating plastic material behavior in a hydro-mechanical analysis can be defined by a variant of the Mohr-Coulomb law according to: where τ crit is the critical shear stress at failure, σ n is the normal stress, ϕ' is the effective angle of internal friction and c' is the effective cohesion [76]. Injection leads to an increase in pore pressure, which reduces the effective stresses and the effective normal stress. In the σ n -τ diagram, this pore pressure increase shifts the Mohr Circle towards the left until it ultimately reaches the failure line and shear failure occurs. Tensile failure is incorporated via a tension cut-off delimiting the Mohr-Coulomb failure line for tensile stresses (Figure 4).
where the effective stresses (σ'ij) are derived from the total stress tensor (σij) by subtracting the pore pressure (p), weighted by the Biot coefficient ( ) and Kronecker's delta (δij) [74][75][76]. The mechanical material behavior in the poroelastic domain is described by the following stressstrain relationship: It links strain (εij) and effective stress (σ'ij) through the material properties of Young's modulus (E) and Poisson's ratio ( ) and the sum of the effective principal stresses (σ'kk) [77][78][79]. The failure criterion delimiting the poroelastic range and initiating plastic material behavior in a hydromechanical analysis can be defined by a variant of the Mohr-Coulomb law according to: where Ʈcrit is the critical shear stress at failure, σn is the normal stress, φ' is the effective angle of internal friction and c' is the effective cohesion [76]. Injection leads to an increase in pore pressure, which reduces the effective stresses and the effective normal stress. In the σn-Ʈ diagram, this pore pressure increase shifts the Mohr Circle towards the left until it ultimately reaches the failure line and shear failure occurs. Tensile failure is incorporated via a tension cut-off delimiting the Mohr-Coulomb failure line for tensile stresses (Figure 4). Besides the hydrostatic pressure, also the injection pressure affects the model and the increase in pore pressure propagates through the model. This is described by Darcy's law as fluid flow (q) is determined according to: where k is the intrinsic permeability of the porous medium, n the fluid viscosity, p the pore pressure, f the fluid viscosity, S the specific storage (as a function of porosity) and t the time [77]. The fully coupled hydro-mechanical simulation accounts for 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 [11,76,80].  Besides the hydrostatic pressure, also the injection pressure affects the model and the increase in pore pressure propagates through the model. This is described by Darcy's law as fluid flow (q) is determined according to: where k is the intrinsic permeability of the porous medium, n the fluid viscosity, p the pore pressure, ρ f the fluid viscosity, S the specific storage (as a function of porosity) and t the time [77].
The fully coupled hydro-mechanical simulation accounts for 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 [11,76,80].

Material Input Parameters
The hydraulic and mechanical material properties assigned to the various model units are listed in Table 2. For the mechanical properties, only Young's modulus and the frictional properties are varied between the fault zone on one side and the reservoir rock and over-/underburden on the other. For ease of comparison, Poisson's ratio and density are kept the same for all three units. A friction angle of 10 • was used for the fault zone, which is substantially lower than the typical values for most intact rocks [38]. However, various studies have pointed out that processes like cataclasis and formation of clay minerals, among others, lower the friction angle and reduce the strength of fault rocks, thus leading to mechanically weak fault zones [81][82][83]. In this study, we focus on mechanical weak fault zones, where plastic deformation and therefore fault reactivation are most likely to occur. In order to study fluid flow along the fault zone and how it alters its mechanical strength by changes in effective stress, we decided to choose a permeable fault zone. However, impermeable faults frequently also occur in sandstone-shale regimes [82,84,85].

Results
The following section presents the simulation results for the various scenarios. First, the stress and strain patterns as well as the pore pressure distribution of a reference model are shown, which provides a baseline for the following comparison between different model setups. The modeling results are visualized by means of contour and vector plots as well as detailed sections of the fault zone. Please refer to Table 1 for relating the model name to the details of the different model setups.

Reference Model (C3-60)
The model setup C3-60 describing the fault zone as a 3-element wide unit embedded in a curvilinear grid is used as a reference model. Applying the material properties and boundary conditions outlined in Sections 4.2 and 4.4 (see also Figure 3) leads to the initial model state, i.e., before fluid injection into the reservoir horizon starts, shown in Figure 5a-c. As was defined, the hydrostatic pore pressure increases linearly from 9.81 MPa at the model top (1000 m depth) to 19.62 MPa at the base (2000 m depth) (Figure 5a). The effective maximum principal stress (S 1,eff ) ranges from 15 MPa at the top to 37 MPa at the bottom of the model (Figure 5b). Due to its higher Biot coefficient, the fault zone exhibits lower effective stress magnitudes than the surrounding rock. Some minor perturbations in the stress pattern occur in the immediate vicinity of the fault. For the initial and boundary conditions selected, the fault zone is entirely in an elastic state prior to fluid injection. This ensures, that initial loading does not already cause plastic deformation and fault reactivation (Figure 5c).
Due to the assumed normal faulting regime, the regional maximum principal stress orientation is vertical. Only the fault cells properly exhibit a slight rotation of about 5° relative to regional ( Figure  6b).
Examining the model state after 5 years of injection, both the minimum and maximum peaks at the upper fault tip in the close surrounding of the fault zone are around 15 to 20 MPa higher than the magnitudes for S1,eff (Figure 6c). The orientation of S1 near the fault tip has rotated up to 30° clockwise at the downthrown side. The rotation decreases towards the right until it reaches vertical again. On the left side of the fault tip, the rotation is lower, i.e., about 5° counterclockwise near the fault ( Figure  6d).    (Figure 5d). The pore pressure increase in the reservoir as well as in the fault zone is also indicated in the S 1,eff pattern, as these model parts show 5 to 10 MPa lower stress magnitudes than the surrounding rock mass and for the initial state (Figure 5e). While the reservoir rock as well as the over-/underburden remain in an elastic state, pore pressure increase results in plastic straining of the fault zone elements. The corresponding maximum effective stress pattern shows an increase at the upper left and lower right end of the fault and a decrease over the entire fault zone. Plastic straining occurs in the whole fault zone with deformation ranging between 0.5% within the fault zone and about 3% at the fault tips (Figure 5f).
In Figure 6, the magnitude of the maximum total principal stress (S 1 ) is displayed together with a detailed view of the orientation of S 1 at the fault tips for both the initial (Figure 6a,b) and the final calculation step (Figure 6c,d) of the base model approach (C3-60). The basic pattern is similar to the magnitudes for S 1,eff , albeit at higher values as pore pressure effects are not considered (Figure 6a). Due to the assumed normal faulting regime, the regional maximum principal stress orientation is vertical. Only the fault cells properly exhibit a slight rotation of about 5 • relative to regional (Figure 6b). Energies 2020, 13, 2673 13 of 28

Mesh Resolution
The impact of the mesh resolution of the fault zone is compared by implementing different element sizes for the same total width of the fault zone, i.e., the number of elements used to resolve the fault zone varies. Figures 7 and 8 display the results after the last load step, i.e., after 5 years of injection, for the magnitude of S1,eff and the orientation of S1 for different mesh resolutions. Thereby, the fault is represented either as 1-, 3-or 9-element wide zone in a rectangular grid. In general, significant differences between the 1-element (R1-60) on one side and the 3-(R3-60) and 9-element (R9-60) width on the other can be observed.  (Figure 6d).

Mesh Resolution
The impact of the mesh resolution of the fault zone is compared by implementing different element sizes for the same total width of the fault zone, i.e., the number of elements used to resolve the fault zone varies. Figures 7 and 8 display the results after the last load step, i.e., after 5 years of injection, for the magnitude of S 1,eff and the orientation of S 1 for different mesh resolutions. Thereby, the fault is represented either as 1-, 3-or 9-element wide zone in a rectangular grid. In general, significant  After the last load step, approach R1-60 displays a similar stress pattern, with subhorizontal contour lines like the initial step of the base model. Only the reservoir horizon in the downthrown block shows a reduction of the effective stress because of injection and the resulting increase in pore pressure. The upper left reservoir horizon, however, has not experienced any pore pressure increase (Figure 7a).
In contrast to the 1-element wide fault models, both, the 3-and 9-element wide fault zones clearly undergo a pore pressure increase inside the fault zone. In both models, the entire fault zone is precisely visible by the corresponding decrease in the magnitude of the effective maximum principal stress (S1,eff). Peaks of S1,eff occur at the upper and lower fault tips after 5 years of injection (Figure 7b,c). Those peaks are missing in the results for the 1-element width fault zone.
Comparing both higher mesh density models R3-60 and R9-60, only minor differences are found for the S1,eff values. Somewhat higher and, respectively, lower values are produced at the fault tips for the higher mesh resolution. Mesh resolution R3-60 has maximum values of 30 MPa for the upper and 50 MPa for the lower fault tip, while for mesh resolution R9-60 the corresponding values are 40 and 60 MPa. Although these small areas up to 10 m around the fault tips have different maximum values for S1,eff, the overall appearance of the S1,eff-pattern is the same for the multiple element row approaches (Figure 7b,c).
Likewise, approach R1-60 does not show rotation of the vertical orientation of S1 except a very slight one (less than 5°) for the actual fault cells for the last load step (Figure 8a). The R3-60 and R9-60 approaches exhibit rotations of S1 similar to the base model. Both display up to 30° clockwise rotations on the right side of the upper fault tip and between 5° and 10° counter clockwise rotation  After the last load step, approach R1-60 displays a similar stress pattern, with subhorizontal contour lines like the initial step of the base model. Only the reservoir horizon in the downthrown block shows a reduction of the effective stress because of injection and the resulting increase in pore pressure. The upper left reservoir horizon, however, has not experienced any pore pressure increase (Figure 7a).
In contrast to the 1-element wide fault models, both, the 3-and 9-element wide fault zones clearly undergo a pore pressure increase inside the fault zone. In both models, the entire fault zone is precisely visible by the corresponding decrease in the magnitude of the effective maximum principal stress (S1,eff). Peaks of S1,eff occur at the upper and lower fault tips after 5 years of injection (Figure 7b,c). Those peaks are missing in the results for the 1-element width fault zone.
Comparing both higher mesh density models R3-60 and R9-60, only minor differences are found for the S1,eff values. Somewhat higher and, respectively, lower values are produced at the fault tips for the higher mesh resolution. Mesh resolution R3-60 has maximum values of 30 MPa for the upper and 50 MPa for the lower fault tip, while for mesh resolution R9-60 the corresponding values are 40 and 60 MPa. Although these small areas up to 10 m around the fault tips have different maximum values for S1,eff, the overall appearance of the S1,eff-pattern is the same for the multiple element row approaches (Figure 7b,c).
Likewise, approach R1-60 does not show rotation of the vertical orientation of S1 except a very slight one (less than 5°) for the actual fault cells for the last load step (Figure 8a). The R3-60 and R9-60 approaches exhibit rotations of S1 similar to the base model. Both display up to 30° clockwise rotations on the right side of the upper fault tip and between 5° and 10° counter clockwise rotation After the last load step, approach R1-60 displays a similar stress pattern, with subhorizontal contour lines like the initial step of the base model. Only the reservoir horizon in the downthrown block shows a reduction of the effective stress because of injection and the resulting increase in pore pressure. The upper left reservoir horizon, however, has not experienced any pore pressure increase (Figure 7a).
In contrast to the 1-element wide fault models, both, the 3-and 9-element wide fault zones clearly undergo a pore pressure increase inside the fault zone. In both models, the entire fault zone is precisely visible by the corresponding decrease in the magnitude of the effective maximum principal stress (S 1,eff ). Peaks of S 1,eff occur at the upper and lower fault tips after 5 years of injection (Figure 7b,c). Those peaks are missing in the results for the 1-element width fault zone.
Comparing both higher mesh density models R3-60 and R9-60, only minor differences are found for the S 1,eff values. Somewhat higher and, respectively, lower values are produced at the fault tips for the higher mesh resolution. Mesh resolution R3-60 has maximum values of 30 MPa for the upper and 50 MPa for the lower fault tip, while for mesh resolution R9-60 the corresponding values are 40 and 60 MPa. Although these small areas up to 10 m around the fault tips have different maximum values for S 1,eff , the overall appearance of the S 1,eff -pattern is the same for the multiple element row approaches (Figure 7b,c).
Likewise, approach R1-60 does not show rotation of the vertical orientation of S1 except a very slight one (less than 5 • ) for the actual fault cells for the last load step (Figure 8a). The R3-60 and R9-60 approaches exhibit rotations of S 1 similar to the base model. Both display up to 30 • clockwise rotations Energies 2020, 13, 2673 14 of 27 on the right side of the upper fault tip and between 5 • and 10 • counter clockwise rotation on the left side, respectively. The orientation rotates back to the vertical orientation with increasing distance from the fault (Figure 8b,c). Figure 9 shows the results for the von Mises plastic strain after the last load step for the three mesh resolutions examined. The 1-element wide fault zone of model R1-60 does not show any plastic straining at all and, hence, is not reactivated (Figure 9a). This contrasts significantly with the strain patterns calculated for the 3-and 9-element fault zone models. Thereby, the finer meshed fault zone shows a larger spatial extent of higher plastic strain values reaching up to 1.3% inside the fault zone (Figure 9c). Apart from this difference, the overall results for the plastic strain are rather similar, i.e., both approaches show larger plastic straining in the fault center, which is reducing towards the fault tips (Figure 9b,c). on the left side, respectively. The orientation rotates back to the vertical orientation with increasing distance from the fault (Figure 8b,c). Figure 9 shows the results for the von Mises plastic strain after the last load step for the three mesh resolutions examined. The 1-element wide fault zone of model R1-60 does not show any plastic straining at all and, hence, is not reactivated (Figure 9a). This contrasts significantly with the strain patterns calculated for the 3-and 9-element fault zone models. Thereby, the finer meshed fault zone shows a larger spatial extent of higher plastic strain values reaching up to 1.3% inside the fault zone (Figure 9c). Apart from this difference, the overall results for the plastic strain are rather similar, i.e., both approaches show larger plastic straining in the fault center, which is reducing towards the fault tips (Figure 9b,c).

Pore Pressure and Effective Stress Magnitude
In general, the increase in pore pressure by injection into the downthrown reservoir section propagates through the fault zone in both models. However, there is a slight difference between the rectangular and curvilinear approach in the absolute value reached in the upper reservoir section. While approach C9-60 shows pore pressures up to 25 MPa in this part of the reservoir, they remain below 20 MPa in approach R9-60 stays. Within the fault zone proper, the pore pressure increase expands more towards the fault tips in model C9-60 e.g., the pore pressure for C9-60 is between 27.5 to 32.  (Figure 10a,c).
The stress magnitudes for S1,eff display the same overall behavior, but due to the higher pore pressures, S1,eff is lower in the central part of the fault zone for the curvilinear approach (C9-60). There, magnitudes of less than 10 MPa are observed and the spatial extent is up to 12 m wide to both sides of the fault. The magnitudes achieved with the rectangular approach (R9-60) only gain around 10 to 15 MPa in the middle of the fault and the spatial extent is only about 8 m. Considering the magnitudes of S1,eff at the fault tips for both grid geometries, the approach C9-60 exhibits higher peaks with 68 MPa for the lower and 47 MPa for the upper fault tip. The values observed for the rectangular approach are around 2 MPa lower with 66 and 45 MPa respectively (Figure 10b,d).
A detailed view of S1,eff in the vicinity of the lower fault tip is presented in Figure 11. The spatial extent for both the minimum and maximum peak values is smaller for approach C9-60. The lower values (< 10 MPa) for S1,eff taper towards the fault tip. In contrast, R9-60 has a wider zone of low values (< 10 MPa) which seems to expand towards the fault tip. At distances of more than 10 m from the fault zone, both approaches exhibit the same pattern.

Pore Pressure and Effective Stress Magnitude
In general, the increase in pore pressure by injection into the downthrown reservoir section propagates through the fault zone in both models. However, there is a slight difference between the rectangular and curvilinear approach in the absolute value reached in the upper reservoir section. While approach C9-60 shows pore pressures up to 25 MPa in this part of the reservoir, they remain below 20 MPa in approach R9-60 stays. Within the fault zone proper, the pore pressure increase expands more towards the fault tips in model C9-60 e.g., the pore pressure for C9-60 is between 27.  (Figure 10a,c).
The stress magnitudes for S 1,eff display the same overall behavior, but due to the higher pore pressures, S 1,eff is lower in the central part of the fault zone for the curvilinear approach (C9-60). There, magnitudes of less than 10 MPa are observed and the spatial extent is up to 12 m wide to both sides of the fault. The magnitudes achieved with the rectangular approach (R9-60) only gain around 10 to 15 MPa in the middle of the fault and the spatial extent is only about 8 m. Considering the magnitudes of S 1,eff at the fault tips for both grid geometries, the approach C9-60 exhibits higher peaks with 68 MPa for the lower and 47 MPa for the upper fault tip. The values observed for the rectangular approach are around 2 MPa lower with 66 and 45 MPa respectively (Figure 10b,d).
A detailed view of S 1,eff in the vicinity of the lower fault tip is presented in Figure 11. The spatial extent for both the minimum and maximum peak values is smaller for approach C9-60. The lower

Stress Orientation
The magnitude and orientation of the maximum total principal stress (S1) is shown in Figure 12  achieved for the curvilinear approach. From the fault to the sides of the model, the orientation rotates back towards the vertical orientation of S 1 again (Figure 12b,d).
Energies 2020, 13, 2673 17 of 28 Figure 11. Detailed view on the results after the last load step for a 9-element wide fault zone with 60° dip embedded with different grids. Rectangular grid: (a) magnitude S1,eff (b) von Mises plastic strain. Curvilinear grid: (c) magnitude S1,eff (d) von Mises plastic strain.

Stress Orientation
The magnitude and orientation of the maximum total principal stress (S1) is shown in Figure 12 for the last load step. The S1 magnitude distribution shows similar patterns for both grid types. The values range from about 25 MPa at the model top to about 45 MPa at the model bottom. Both models exhibit similar S1 peaks at the fault tips, which reach 56 and 58 MPa at the upper fault tip and 77 and 78 MPa at the lower fault tip, respectively (Figure 12a,c). Likewise, the orientation of S1 in the vicinity of the fault tips is very similar. However, there are slight differences in the actual numbers. For the rectangular grid the rotation to the vertical orientation is about 5° counter clockwise on the left and 20°-25° clockwise on the right side of the fault. While left of the fault tip, the rotation is also 5° counter clockwise, the reorientation differs on the right side, where between 25° and 30° clockwise are achieved for the curvilinear approach. From the fault to the sides of the model, the orientation rotates back towards the vertical orientation of S1 again (Figure 12b,d).

Elastic and Plastic Strain
After the last load step, the von Mises elastic strain shows the highest values at the fault tips for both fault grid geometries. At the lower tip, the elastic strain is up to 0.13%. Comparing the elastic strain throughout the fault zone, it is obvious that R9-60 achieves higher values towards the middle of the fault. In contrast, elastic strain in the curvilinear approach decreases strongly from the lower tip towards the middle of the fault and to a lesser extent from the upper tip towards the fault center (Figure 13a,c).
The main difference between both approaches is observed for the von Mises plastic strain. After the last load step, the plastic strain occurs mostly at the fault tips for C9-60. The plastic strain for the curvilinear approach reaches up to 3.2% at the lower fault tip. Towards the middle of the fault, plastic strain decreases to about 0.3%. For R9-60, the plastic strain behaves the other way around. It decreases from the middle of the fault towards the fault tips. The rectangular approach has its highest values of about 1.5% in the middle of the fault and decreases to about 0.4% at the fault tips. Outside the fault zone no plastic straining occurs (Figure 13b,d). A more detailed view of these plastic strain patterns in the vicinity of the lower fault tip is pictured in Figure 11b,d. This figure shows the von Mises plastic strain for both approaches within 50 m from the fault tip.
Energies 2020, 13, 2673 18 of 28 Figure 12. Results after the last load step for a 9-element wide fault zone with 60° dip embedded with different grids. Rectangular grid: (a) magnitude of the maximum principal stress (S1). (b) detailed view of the orientation of S1. Curvilinear grid: (c) magnitude of S1. (d) detailed view of the orientation of S1.

Elastic and Plastic Strain
After the last load step, the von Mises elastic strain shows the highest values at the fault tips for both fault grid geometries. At the lower tip, the elastic strain is up to 0.13%. Comparing the elastic strain throughout the fault zone, it is obvious that R9-60 achieves higher values towards the middle of the fault. In contrast, elastic strain in the curvilinear approach decreases strongly from the lower tip towards the middle of the fault and to a lesser extent from the upper tip towards the fault center (Figure 13a,c).
The main difference between both approaches is observed for the von Mises plastic strain. After the last load step, the plastic strain occurs mostly at the fault tips for C9-60. The plastic strain for the curvilinear approach reaches up to 3.2% at the lower fault tip. Towards the middle of the fault, plastic strain decreases to about 0.3%. For R9-60, the plastic strain behaves the other way around. It decreases from the middle of the fault towards the fault tips. The rectangular approach has its highest values of about 1.5% in the middle of the fault and decreases to about 0.4% at the fault tips. Outside the fault zone no plastic straining occurs (Figure 13b,d). A more detailed view of these plastic strain patterns

Fault Dip
The effect of different fault dips is investigated using the curvilinear approach with a 3-element wide fault zone with 60 • (C3-60) and 90 • fault dip (C3-90), respectively. Figure 14 shows the results for pore pressure, S 1,eff and total plastic strain after the last load step. Similar results are achieved for the pore pressure distribution, albeit somewhat higher pore pressures are observed in the upper reservoir section for the vertical fault (Figure 14a,d). For S 1,eff , approach C3-60 shows peak values of up to 60 MPa for the lower and 50 MPa for the upper fault tip. Those local stress peaks do not occur in case of a vertical fault, i.e., C3-90. Otherwise, the effective stress magnitudes within the fault zone of about 5 MPa are quite similar for both fault dips (Figure 14b,e).
Energies 2020, 13, 2673 19 of 28 in the vicinity of the lower fault tip is pictured in Figure 11b,d. This figure shows the von Mises plastic strain for both approaches within 50 m from the fault tip. Figure 13. Results after the last load step for a 9-element wide fault zone with 60° dip embedded with different grids. Rectangular grid: (a) total von Mises elastic strain (S1). (b) total von Mises plastic strain. Curvilinear grid: (c) total von Mises elastic strain (S1). (d) total von Mises plastic strain.

Fault Dip
The effect of different fault dips is investigated using the curvilinear approach with a 3-element wide fault zone with 60° (C3-60) and 90° fault dip (C3-90), respectively. Figure 14 shows the results for pore pressure, S1,eff and total plastic strain after the last load step. Similar results are achieved for the pore pressure distribution, albeit somewhat higher pore pressures are observed in the upper reservoir section for the vertical fault (Figure 14a

Discussion
For the initial load step, the expected model behavior is observed. The hydrostatic pressure and the loading result in vertical gradients of both pore pressure and the maximum principal stress (S1 = Sv). Outside the area affected by the predefined fault zone, S1 is vertical and the stress regime displays normal faulting [73]. No plastic straining is observed after the initial load step which proves that the fault is not already reactivated by the initial and boundary conditions selected. The results of the last load step, i.e., after 5 years of injection, show higher pore pressure in both reservoir horizons, and Substantial differences are again observed for plastic straining, as no plastic straining occurs for a vertical fault. In contrast, von Mises plastic strain in case of a fault zone dipping at 60 • is 0.6% in the middle of the fault and up to 3.2% at both fault tips (Figure 14c,f).

Discussion
For the initial load step, the expected model behavior is observed. The hydrostatic pressure and the loading result in vertical gradients of both pore pressure and the maximum principal stress (S 1 = S v ). Outside the area affected by the predefined fault zone, S 1 is vertical and the stress regime displays normal faulting [73]. No plastic straining is observed after the initial load step which proves that the fault is not already reactivated by the initial and boundary conditions selected. The results of the last load step, i.e., after 5 years of injection, show higher pore pressure in both reservoir horizons, and therefore, indicates fluid migration through the fault zone. Plastic straining occurs throughout the entire fault zone. Hence, fault reactivation is achieved from pore pressure increase due to fluid injection into the reservoir.
The calculated stress pattern can be compared to other numerical simulations which incorporate faults [1,50,86]. The corresponding stress perturbations derived from such numerical models resemble the results of the base model regarding both spatial extent and magnitude of the stress changes induced by the fault.

Mesh Resolution
A general recommendation is to use a finer mesh in those areas where stress and strain gradients are large. To identify the regions where greater mesh density or local grid refinement is required, preliminary simulations with a coarse mesh appear useful [53,[87][88][89].
In our study, the representation of a fault as 1-element row in a rectangular grid appears to be a special case, but they are not unusual for reservoir simulations. In the last several years, different authors have used this kind of fault zone representation to investigate a multitude of reservoir-related tasks like fault reactivation [49,60], CO2-Storage [59] and fully-coupled reservoir simulations [29,49]. However, our modeling results indicate substantial differences compared to fault representations with multiple element rows. For the mechanical part, the stress and strain patterns as well as magnitudes are different and for the hydraulic part, the fault seems to act as barrier for fluid flow in the case of a 1-element representation.
Only minor differences in the resulting stress and strain patterns as well as stress magnitudes exist if the fault is represented as a 3-element or 9-element wide zone. The higher mesh resolution inside the fault zone only increases the maximum values for both stress and strain at the fault tips, but the spatial extent of these differences remains less than 10 m near the fault tips.
Several studies [87,90,91] have evaluated the effects of mesh density on finite element analysis. There are two issues which can be responsible for the punctually higher stress and strain values in the simulations using the higher mesh resolution. First, stress and strain fields have higher gradients in the localization zone, such as the tip of a fault zone, for higher mesh densities [92]. This implies that the finer mesh resolution leads to a more punctually concentrated stress, which in turn, is higher [93,94]. The higher stress values ultimately cause higher plastic strain, which is also observed in the finer mesh resolution. The second issue is that models with a lower mesh resolution, i.e., with less elements, appear to be somewhat stiffer, while increasing the number of elements softens the model slightly and improves the accuracy of the stiffness integration [89,95]. Therefore, the slightly higher stress and strain values for the models using a higher resolution may originate from the lower effective stiffness at the fault tips.
However, the spatial extent of this variations is very small and the overall stress and strain pattern as well as the hydraulic behavior for the whole model is the same for 3-and 9-element wide fault zones. In order to determine the required mesh resolution inside the fault zone, it is crucial to evaluate the required precision as well as the target location regarding the specific aim of the particular study. The differences between the one element (R1-60) approach and both multiple row approaches (R3-60 and R9-60) stem from the arrangement of the elements to which the fault zone properties have been assigned. Regarding R1-60, the rectangular grid geometry leads to the fault zone elements being connected by only one common node, i.e., all neighboring elements which share joint edges exhibit the non-permeable and mechanically stronger host rock properties. This subdues the formation of both a through-going fluid pathway and deformation zone. In contrast, the fault zone with a higher mesh resolution and more element rows has numerous common element faces and, hence, a continuous fluid pathway through the entire fault zone can develop. Likewise, strain can accumulate and form a continuous zone of deformation. The larger amount of fluid flow through the fault zone increases the pore pressure inside the fault zone and leads to reduced effective stresses, which in turn, decreases the shear strength again.

Grid Geometry
Comparing the two fundamental grid geometries studied, the different results in both hydraulic and mechanical simulations can be explained by the special features of the rectangular grid for fluid flow as well as stress and strain propagation.
On the hydraulic part, the rectangular grid induces a spatial restriction in the fluid pathway since some of the elements forming the fault zone are connected by only one node. The reduced fluid migration through the rectangular fault zone also explains the lower pore pressures observed in the upper reservoir section for the rectangular grid geometry (see Figure 10a).
In a coupled simulation the hydraulic behavior directly effects the mechanical response. According to the Mohr-Coulomb criteria (Figure 4), a pore pressure increase causes a decrease in the effective stresses, which shifts the Mohr-Coulomb circle towards the shear failure line [21,38,76]. Thus, the higher pore pressures in the curvilinear approach lead to a higher reduction of the effective stresses [19,78,96], which explains the difference between both gird geometries regarding the magnitudes and the spatial extent of S 1,eff (Figure 10b,d).
For the rectangular grid, the increase in pore pressure is mainly in the center of the fault which is also the reason for the location of the highest plastic strain values there. However, this does not fully explain the different plastic strain distribution for both grid geometries. Fault geometry may also be partly responsible for this different mechanical response. While the curvilinear approach offers a linear, smooth boundary between the undeformed reservoir rock and the fault zone, the rectangular approach has a stair-stepped geometry ( Figure 15). This results in some kind of interlocking with the stronger rock properties outside the fault zone, and therefore, constrains strain accumulation compared to the curvilinear approach.
In addition, the grid geometry seems to be responsible for the different results of the plastic strain occurring at the fault tips (Figure 11b,d and Figure 13b,d). The rectangular fault zone ends with a rectangular block of elements. This favors the plastic strain to disperse over multiple elements, which in combination with interlocking of the host rock cells, reduces the plastic strain values at the fault tips even further ( Figure 15R). In contrast, the fault in the curvilinear fault converges towards the fault tip ( Figure 15C). This causes the strain to accumulate at the fault tips, and thus, exhibits larger plastic strain there.
Particularly the plastic strain pattern of the curvilinear approach seems to be closer to nature. During the development and the growth of a fault, stress concentrations at the fault tips lead to strain localization [97][98][99]. Small discontinuities (i.e., Griffith cracks) in the so-called "intact rock" can be stimulated and coalesced by those strain localizations [100][101][102] similar to the ones observed for the curvilinear approach (see Figure 13d). This macroscale shear failure of the intact host rock at the fault tips can expend through ongoing strain to extend a continuous fault surface [103]. Although the FE method used in our approach cannot model fracturing and consequently, fault propagation, the plastic strain accumulations at the fault tips for curvilinear fault representations seem to mimic the behavior in nature. stimulated and coalesced by those strain localizations [100][101][102] similar to the ones observed for the curvilinear approach (see Figure 13d). This macroscale shear failure of the intact host rock at the fault tips can expend through ongoing strain to extend a continuous fault surface [103]. Although the FE method used in our approach cannot model fracturing and consequently, fault propagation, the plastic strain accumulations at the fault tips for curvilinear fault representations seem to mimic the behavior in nature.

Fault Dip
Different fault dip angles do not significantly influence fluid flow. Therefore, the pore pressure distribution is similar for the 60° and 90° dip scenarios. However, the dip of the fault zone has a strong influence on the mechanical behavior. Both stress and strain patterns are reduced for the 90° fault dip compared to the 60° fault dip. Particularly the difference in the plastic strain is striking. There is no plastic straining and, consequently, no fault reactivation for the 90° case. This behavior can be explained according to the Mohr-Coulomb failure criterion [21,38]. The stresses acting on the fault plane are divided into the normal stress (σn) acting orthogonal to the fault and the shear stress (τ) acting parallel on the fault, which can be calculated according to:

Fault Dip
Different fault dip angles do not significantly influence fluid flow. Therefore, the pore pressure distribution is similar for the 60 • and 90 • dip scenarios. However, the dip of the fault zone has a strong influence on the mechanical behavior. Both stress and strain patterns are reduced for the 90 • fault dip compared to the 60 • fault dip. Particularly the difference in the plastic strain is striking. There is no plastic straining and, consequently, no fault reactivation for the 90 • case. This behavior can be explained according to the Mohr-Coulomb failure criterion [21,38]. The stresses acting on the fault plane are divided into the normal stress (σ n ) acting orthogonal to the fault and the shear stress (τ) acting parallel on the fault, which can be calculated according to: with maximum (σ 1 ) and minimum (σ 3 ) principal stresses and the fault dip (α) [73,94]. According to the model setup, the initial stress field is the same for both, the 60 • -and 90 • -fault dip. If σ 1 and σ 3 acting on the fault are the same, the only different variable for the different fault dips in Formulas (8) and (9) is the fault dip (α) itself. The shear stress, calculated by formula (8), is zero when the sinus-term is zero, which is the case for a 90 • fault dip. Since the shear stress acting on the fault is zero, the shear failure line is not exceeded, and the rock remains intact. Hence, the deformation takes only part of the elastic domain, no plastic strain occurs for the models with a 90 • fault dip.

Practical Aspects of Model Building
Between the various approaches there are huge differences due to the amount of work and time needed, i.e., the costs required to incorporate faults into a hydro-mechanical reservoir model. Computing time increases exponentially with the number of elements, so finer mesh resolutions within and near the fault zone significantly increase the runtime of the simulations. Even if it is not readily available from preceding property modeling or a flow simulation, generating the grid geometry for a regular rectangular model and intersecting it with a fault interpretation is quite straightforward and by far the most rapid technique. Implementing the curvilinear approach for fault representations requires a special grid adapted to the fault geometry. This step can be a rather time-consuming and labor-intensive task, depending on the complexity of the geometry of the fault network being modelled. As shown above, a rectangular grid with the same fault zone resolution can lead to similar results, at least at larger distances from the fault zone.
The rectangular approach can therefore be suitable for a first-order evaluation before more complex models with a curvilinear grid are performed. In addition, for reservoir simulations focusing on the central parts of fault compartments rather than the fault zone itself, the rectangular approach can be reasonable. In contrast, for reservoir simulations focusing on fault zones and their immediate vicinity, better results can be expected by using the curvilinear approach.

Conclusions
Using a simple generic model setup, different scenarios for the incorporation of faults as volumetric weak zones into hydro-mechanical reservoir models are analyzed. The various scenarios differ regarding the mesh resolution of the fault zone (1-, 3-or 9-element width), the grid geometry (rectangular vs. curvilinear) and the fault dip (60 • and 90 • ). Significant differences in the stress and strain patterns are indicated in the results, which are induced by the fault depending on its incorporation in the numerical model. Based on the numerical simulation results, five general recommendations can be given on how to represent faults in FE reservoir models: 1.
The mesh resolution has to be considered very carefully, since it can-combined with a rectangular grid-lead to serious errors. It needs to be ensured that the fault cells do not interlock with the surrounding, stronger and less permeable host rock as this effects both fluid flow through and straining of the fault zone. This interlocking effect mainly occurs for 1-element width fault zones. For grids with multiple element wide fault zones, the only differences are observed in the vicinity of the fault tips. Three element wide fault zones appear to be appropriate for most reservoir-scale models. Only if the aim of the study is within 10 m of the fault zone, a finer resolution should be considered.

2.
If the aim of the study is to model the fault zone properly, a curvilinear representation is recommended. In addition to somewhat better fluid migration throughout the whole fault zone, this approach shows higher plastic strain at the fault tips, which appears to be closer to reality.

3.
If the aim of the study is at a distance of more than 10 m from the fault, both grid geometries are interchangeable. They show similar stress and strain patterns. In this case, the advantage of the rectangular grid is that it generally takes much less time to generate.

4.
Different fault dips produce different mechanical results, i.e., stress and strain patterns. Therefore, care should be taken to consider a realistic fault dip, e.g., from interpretation of depth-converted seismic sections, rather than using a vertical fault dip for simplification.

5.
Different fault dips produce similar hydraulic results. If the aim of the study is primarily on hydraulic issues, vertical faults can be an acceptable simplification.
Regarding the upscaling of the material heterogeneity and complex fault zone architecture to a volumetric weak material representing the fault zone, challenges remain. For a more detailed analysis, including architectural and material heterogeneity, a volumetric fault zone description with a grid adapted to the fault zone geometry seems more appropriate. The possibility of incorporating further details can be presented by local mesh refinement rather than upscaling them to the size of one element.
This research investigates a simple generic fault model in a siliciclastic succession, but the methodology and the findings can be transferred to structurally more complex models as well as other lithologies, e.g., carbonate reservoirs. However, while the general recommendations can be used, the reservoir structure has to be modelled according to the specific geological setting. Funding: This research received no external funding. The research was financed from the unrestricted-use budget of the chair of Engineering Geology from TU Darmstadt. It did not receive funding from agencies in the public, commercial, or not-for-profit sectors.