Geomechanics in Depleted Faulted Reservoirs

This paper examines the impact of the effective stresses that develop during depletion of a faulted reservoir. The study is based on finite element modeling using 2D plane strain deformation analysis with pore pressure and elastoplastic deformation of the reservoir and sealing shale layers governed by the Drucker–Prager plasticity model. The mechanical properties and response of the rock formations were derived from triaxial test data for the sandstone reservoirs and correlation functions for the shale layers. A normal fault model and a reverse fault model were built using seismic data and interpretation of field data. The estimated tectonic in-situ stress field was transformed to the plane of the modeled geometry. Sensitivity studies were performed for uncertainties on the values of the initial horizontal stress and for the friction of the fault surfaces. It was found that the stress path during depletion is mainly controlled by the initial lateral stress ratio (LSR). The developed effective stresses with depletion are influenced by the fault geometry of the compartmentalized blocks. Plastic deformation develops for low LSR whereas for high values the system tends to remain in the elastic region. When plastic deformation takes place, it affects mainly the region near the fault. The reservoir deformation is dominated by vertical displacement which is higher near the fault region and nearly uniform in the remote area. The volumetric strain is dominated by compaction. More volatile conditions in relation to change of the friction coefficient and LSR were found for the normal fault geometry.


Introduction
Faults are shear fractures or extension fractures on which shearing displacement has occurred. Faults can enhance the flow communication through rock, or they can form barriers to fluid flow, depending on their thickness and its opening and the material texture and composition within the fault zone. The operators need to assess early in the producing life of a reservoir whether production is affected by the presence of naturally occurring fractures [1] and whether communication between fault block compartments exists. This can have an economic and operational impact to the field. Due to the great importance of faults in compartmentalized reservoirs, a specialized discipline focusing on a fault-seal analysis is being developed [1].
Reservoir depletion increases the stress carried by the load-bearing grain frame of the reservoir rock [2]. The existence of the fault affects the deformation and stress during depletion resulting in nonuniform distribution, expressed by a local stress path. Geomechanical analysis can be extended to identify rock failure conditions that can lead to the creation of new faulted systems in the subsurface formations. Geomechanics play an important role in identifying the stress conditions in a faulted reservoir system and the potential of slip activation of an existing fault. The stress path affects the geomechanical behavior of the reservoir in terms of compaction and associated surface subsidence. Furthermore, the stress paths in the reservoir and the overburden affect the stability of boreholes during drilling and hydrocarbon production. A similar study [3] emphasized that the stress path inside the reservoir and in the overburden has a significant influence on the ability to drill wells in reservoirs that undergo depletion. Additionally, in the reservoir zones the horizontal stress drops as pore pressure decreases which can cause mud losses risk during the drilling. Furthermore, the increase of the effective stresses often causes the problem of sand production. Therefore, the stress changes during depletion is a key issue for well planning after start-up of production as the decline of the reservoir pressure during reservoir depletion increases the effective stresses. Understanding fault reactivation during the exploitation of a reservoir is also important for preventing well damage [4].
Production-induced fault reactivation cannot be explained by the direct local pore pressure change which by itself, tends to stabilize faults, but it requires a global stress equilibrium which takes into account the geometry and the stiffnesses of the reservoir and surrounding rocks [4][5][6][7][8][9]. A similar study [10] found that the stress path may deviate when nonlinear elasticity is considered and that the assumption of elasticity may result in underestimation of stress and strain changes in and around the reservoir. Geomechanical safety for underground gas storage in compartmentalized reservoirs due to fault reactivation is presented by [11]. Furthermore, production optimization in faulted and fractured reservoirs considering geomechanics is studied in [12].
This study presents results for two faults, a normal fault (NF) and a reverse fault (RF), from the Aphrodite field in East Mediterranean, which are considered as representative cases for the other faults in the field. Previous related work that identified geometrical dimensions of the faults was presented in [13]. The applied stress field, magnitude, and orientation, of the maximum horizontal stress (S H ) is based on [14]. This study examines the generated by depletion effective stresses and stress paths, the rock yield tendency and fault slipping potential, and the displacement field and volumetric strain in the reservoir and surrounding zones. Sensitivity analysis of the results is performed for the assumed friction coefficient of the fault surfaces and lateral stress ratios (LSR). The main objective of this study is to show what consequences an elastoplastic analysis may have on the prediction of effective stress changes and deformation in a faulted reservoir and surrounding rocks.

Faults
A fractured system consists of faults and joints that potentially formed the thrust, with the joins created in recent time. Joints develop where the stress regime is in compression. Joints can also develop in areas that rocks are being folded. The main characteristics in faulted systems are (a) the geometry in terms of dip-angle and block compatibility movement, (b) the applied stress field magnitude and its orientation, (c) the fault-friction coefficient, and (d) the rock properties. Faults are characterized by the magnitude of principal stresses S v , S H , and S h which define the applied regime to a faulting system. At normal faulting conditions S v > S H and at reverse faulting conditions S h > S v . Furthermore, strike-slipping is developing at the intermediate state. Additionally, horizontal stress anisotropy is determined by the ratio of S H /S h [15], which in this study is considered to be close to unity as probably an equilibrium of the horizontal stresses in the area is achieved.
In this study the faults-joints system is simplified by considering two major faults. The geometry of the considered two faults is shown in Figure 1 where the seismic profiles are illustrated for (a) the normal fault in NW-SE section and (b) the reverse fault in SW-NE section. The reservoir zone is illustrated within the yellow marked horizons and the fault geometry is shown by black lines. Faults are generally not straight uniform surfaces and commonly range from planar to listric or curved geometries, or conjugate faults. In this case study, shales encase the sandstone reservoir. The reservoir consists of sandstones interbedded by shales but in the modeling is simplified to a uniform zone. The 3D sketch of Figure 2 shows the representative volume which defines the fault contacts between the blocks. The face normal to SH corresponds to the normal fault shown in Figure 1a and the face normal to Sh corresponds to reverse fault of Figure 1b. From the geological understanding, one can assume that the reverse fault developed first followed by the formation of compressional fractures which formed the normal fault perpendicular to the reserve fault plane. Currently, the field stress is recognized to be at normal fault conditions (SH < Sv), that it had passed from the intermediate strike-slip state (SH > Sv). The analysis considers two simplified structural shapes that reduce the computational cost and complexity without losing the qualitative characteristics of the overall process. Each structural model consists of two fault-blocks in a normal fault and reverse fault patterns. The grid dimensions have a length (x-axis) of 4000 m and a total thickness (yaxis) of 900 m. Based on the seismic data, the fault dip for normal fault geometry is 57° and for the reverse fault geometry is 66°. The fault parameters are presented in Table 1. The vertical displacement fault throw is 100 m, hence, the right-side block has been lowered by 100 m from its original depth in the NF case and raised by 100 m in the RF case. The 3D sketch of Figure 2 shows the representative volume which defines the fault contacts between the blocks. The face normal to S H corresponds to the normal fault shown in Figure 1a and the face normal to S h corresponds to reverse fault of Figure 1b. From the geological understanding, one can assume that the reverse fault developed first followed by the formation of compressional fractures which formed the normal fault perpendicular to the reserve fault plane. Currently, the field stress is recognized to be at normal fault conditions (S H < S v ), that it had passed from the intermediate strike-slip state (S H > S v ). The 3D sketch of Figure 2 shows the representative volume which defines the fault contacts between the blocks. The face normal to SH corresponds to the normal fault shown in Figure 1a and the face normal to Sh corresponds to reverse fault of Figure 1b. From the geological understanding, one can assume that the reverse fault developed first followed by the formation of compressional fractures which formed the normal fault perpendicular to the reserve fault plane. Currently, the field stress is recognized to be at normal fault conditions (SH < Sv), that it had passed from the intermediate strike-slip state (SH > Sv). The analysis considers two simplified structural shapes that reduce the computational cost and complexity without losing the qualitative characteristics of the overall process. Each structural model consists of two fault-blocks in a normal fault and reverse fault patterns. The grid dimensions have a length (x-axis) of 4000 m and a total thickness (yaxis) of 900 m. Based on the seismic data, the fault dip for normal fault geometry is 57° and for the reverse fault geometry is 66°. The fault parameters are presented in Table 1. The vertical displacement fault throw is 100 m, hence, the right-side block has been lowered by 100 m from its original depth in the NF case and raised by 100 m in the RF case. The analysis considers two simplified structural shapes that reduce the computational cost and complexity without losing the qualitative characteristics of the overall process. Each structural model consists of two fault-blocks in a normal fault and reverse fault patterns. The grid dimensions have a length (x-axis) of 4000 m and a total thickness (y-axis) of 900 m. Based on the seismic data, the fault dip for normal fault geometry is 57 • and for the reverse fault geometry is 66 • . The fault parameters are presented in Table 1. The vertical displacement fault throw is 100 m, hence, the right-side block has been lowered by 100 m from its original depth in the NF case and raised by 100 m in the RF case. The top of the upper grid surface is at the depth of 4800 m and the top reservoir is at 5100 m. The sealing layers amount of to two shale bodies enclosing the sandstone reservoir. The  Figure 3. The boundaries conditions include vertical rolling for the sides and horizontal rolling for the base of the model. The fault blocks are connected with frictional contact surfaces. The vertical stress, σ 1 is applied as an input on the top of the model along the vertical z-axis. The selected sampling locations are carried out at defined pseudo wells and at the fault slip surface as depicted in Figure 3. Table 1. Geometrical parameters of the examined faults based on Markou and Papanastasiou [13]. The top of the upper grid surface is at the depth of 4800 m and the top reservoir is at 5100 m. The sealing layers amount of to two shale bodies enclosing the sandstone reservoir. The applied loads and boundary conditions in each model geometry are illustrated in Figure 3. The boundaries conditions include vertical rolling for the sides and horizontal rolling for the base of the model. The fault blocks are connected with frictional contact surfaces. The vertical stress, σ1 is applied as an input on the top of the model along the vertical z-axis. The selected sampling locations are carried out at defined pseudo wells and at the fault slip surface as depicted in Figure 3. Pore pressure elements are considered at the given pressure conditions in each depletion scenario. The main assumption is that the formation bodies are in pressure equilibrium in each fault block. For simplifying the current process, the applied pore pressure is controlled as an input in each fault block. The model is based on the Terzaghi effective stress relation for the pore pressure and stresses. The principal stresses are common in both situations. The fault plane dip and strikes are different for NF and RF geometries. The applied transformed stresses perpendicular to the strike orientation present a stress related equilibrium to the current geometrical dimensions. Figure 4 shows the orientation of the in-situ stress field that acts relative to the fault geometry.

Parameter
The lateral confining stresses play an important role in reservoir geomechanics and associated problems such as wellbore stability, sanding, and hydraulic fracturing. Therefore, knowledge of the present-day tectonic regime at depth, can provide a quantitative assessment of the horizontal stress magnitude. A first estimate of the horizontal stresses can be obtained assuming a uniaxial compression during deposition according to the theory of elasticity. In this case the ratio of the horizontal effective stresses to the vertical effective stress is expressed by Ko, as a function of the Poisson ratio ν [16], Other horizontal stress values can be defined by assuming a direct value for Ko or by Ko provided by other theories and assumptions, such as plasticity and limit equilibrium according to Jacky equation [17], Pore pressure elements are considered at the given pressure conditions in each depletion scenario. The main assumption is that the formation bodies are in pressure equilibrium in each fault block. For simplifying the current process, the applied pore pressure is controlled as an input in each fault block. The model is based on the Terzaghi effective stress relation for the pore pressure and stresses.
The principal stresses are common in both situations. The fault plane dip and strikes are different for NF and RF geometries. The applied transformed stresses perpendicular to the strike orientation present a stress related equilibrium to the current geometrical dimensions. Figure 4 shows the orientation of the in-situ stress field that acts relative to the fault geometry. where φ is the effective internal friction angle. Alternatively, excess horizontal stresses can be obtained by simulating the tectonic movement in FE analysis [18]. In this study, Ko is further defined in terms of changes of the effective stresses and presented in the stress path section.

Figure 4.
Applied stress field and transformation to 2D plane models. The transformation angle for the NF is 42° and for the RF is 51°.
Due to the uncertainty in the magnitude of the minimum horizontal stress a sensitivity analysis will be carried out for three different values of the lateral stress ratio (LSR) which is defined as the ratio of the total horizontal stress to the total vertical stress, LSR = SH/Sv. The base case scenario is the most likely one. In this case, the in-situ stresses are inserted directly in the model as residual stresses. Additional uncertainty on the results arises from the value of the friction coefficient of the fault surfaces. All the examined scenarios are shown in Table 2. The lateral confining stresses play an important role in reservoir geomechanics and associated problems such as wellbore stability, sanding, and hydraulic fracturing. Therefore, knowledge of the present-day tectonic regime at depth, can provide a quantitative Energies 2021, 14, 60 5 of 17 assessment of the horizontal stress magnitude. A first estimate of the horizontal stresses can be obtained assuming a uniaxial compression during deposition according to the theory of elasticity. In this case the ratio of the horizontal effective stresses to the vertical effective stress is expressed by Ko, as a function of the Poisson ratio ν [16], Other horizontal stress values can be defined by assuming a direct value for Ko or by Ko provided by other theories and assumptions, such as plasticity and limit equilibrium according to Jacky equation [17], where ϕ is the effective internal friction angle. Alternatively, excess horizontal stresses can be obtained by simulating the tectonic movement in FE analysis [18]. In this study, Ko is further defined in terms of changes of the effective stresses and presented in the stress path section. Due to the uncertainty in the magnitude of the minimum horizontal stress a sensitivity analysis will be carried out for three different values of the lateral stress ratio (LSR) which is defined as the ratio of the total horizontal stress to the total vertical stress, LSR = S H /S v . The base case scenario is the most likely one. In this case, the in-situ stresses are inserted directly in the model as residual stresses. Additional uncertainty on the results arises from the value of the friction coefficient of the fault surfaces. All the examined scenarios are shown in Table 2. Corresponding to LSR scenarios, the effective stresses are shown in Figure 5 where σ v is shown by a continuous line and the value of σ 3 by dashed lines. The dotted line represents the elastic solution for σ 3 in the case that no tectonic stresses affected the field.  6 of 17 Figure 6 shows how the LSR evolves from the initial in-situ conditions, for LSR = 0.76, 0.84, and 1.0., with reservoir depletion, assuming different stress paths, for Ko = 0, 0.3, and 0.5. From this graph one can realize how the total horizontal stresses decrease during the reservoir depletion. These results are based on the assumption that the stress state and pore pressure are uniform, ignoring the formation geometry, parameters, and stiffnesses and the stress equilibrium in the horizontal direction. As it will be shown later when these assumptions are relaxed in the FE analysis the stress path is not uniform and it varies across the layers.

Elastoplasticity
The rock stress-deformation with reservoir depletion is performed with finite element analysis (FEA) that incorporates an elastoplasticity theory for a nonassociated hardening Drucker-Prager material. The use of plasticity provides more accurate estimation of the stress and strain changes in the reservoir and overburden. As it is mentioned earlier, a pure elastic model can lead to significant deviation of the results [10]. The material parameters and functions were derived from calibration of triaxial test data for the sandstone reservoirs and from acoustic velocities and correlation functions for the shale layers.
The experimental triaxial test data for the reservoir sandstone are shown by the stress-strain curves of Figure 7 in terms of the axial and radial strains as a function of the applied differential stress. The confinement stress which refers to σ′3 represents the minimum effective principal stress. A triaxial test at low confining pressure, typically shows a concave stress-strain behavior during initial loading, due to closure of microcracks, with stiffening as the stress increases [15].
The relationships for the Poisson ratio, ν, is given in Equation (3) and the Young modulus, E is calculated from Equation (4) Figure 6. LSR evolution with reservoir depletion assuming different stress paths.

Elastoplasticity
The rock stress-deformation with reservoir depletion is performed with finite element analysis (FEA) that incorporates an elastoplasticity theory for a nonassociated hardening Drucker-Prager material. The use of plasticity provides more accurate estimation of the stress and strain changes in the reservoir and overburden. As it is mentioned earlier, a pure elastic model can lead to significant deviation of the results [10]. The material parameters and functions were derived from calibration of triaxial test data for the sandstone reservoirs and from acoustic velocities and correlation functions for the shale layers.
The experimental triaxial test data for the reservoir sandstone are shown by the stressstrain curves of Figure 7 in terms of the axial and radial strains as a function of the applied differential stress. The confinement stress which refers to σ 3 represents the minimum effective principal stress. A triaxial test at low confining pressure, typically shows a concave stress-strain behavior during initial loading, due to closure of microcracks, with stiffening as the stress increases [15]. Using the three confinement tests of σ′3 = 9.0, 17.9, and 26.9 MPa that reflect the expected range of stress conditions from the initial to depleted states, the average elastic Young's modulus is found to be E = 14.5 GPa and the average Poisson ratio ν = 0.24 derived from the experimental data as the best estimates [14]. The relationships for the Poisson ratio, ν, is given in Equation (3) and the Young modulus, E is calculated from Equation (4) as ∆σ 2 = ∆σ 3 = 0, Using the three confinement tests of σ 3 = 9.0, 17.9, and 26.9 MPa that reflect the expected range of stress conditions from the initial to depleted states, the average elastic Young's modulus is found to be E = 14.5 GPa and the average Poisson ratio ν = 0.24 derived from the experimental data as the best estimates [14].
Plastic yielding takes place after the stress state reaches the yield surface, defined here by the linear Drucker-Prager criterion and where p is the mean pressure and t is a function of the Mises equivalent stress given by An extended D-P model uses the third invariant of deviatoric stress where S is the stress deviator defined as and a parameter K (0.8-1) in order to deform the D-P cone according to the stress path to match better the experimental results in Abaqus [19]. The D-P material parameters are related to the well know Mohr-Coulomb strength parameters via the following relationships for the internal cone: Thus, the calculated D-P angle, β for the sandstone is β = 40.6 • . The evolution of the yield surface with plastic deformation is described by the equivalent stress as a function of the equivalent plastic strain ε pl , which is expressed by the D-P shear stress (cohesion), d = 7.4 MPa for initial yield (ε pl = 0), 17.8 MPa for ultimate strength (ε pl = 3%), and 6.4 MPa for softening (ε pl = 7%). Plastic strains can be obtained from an associated flow rule of the Drucker-Prager model. When dilation, ψ = β, the model can provides unrealistic volume increases, so, we used a nonassociative model defined by dilation, ψ given by which for the Tamar sandstone with ϕ = 30 • , results in ψ = 18.2 • . For the shale layers, the rock parameters were determined from acoustic data and correlation functions. For the compressional velocity measurements of Vp = 3.5 km/s [15], the correlation in Equation (15) yields the uniaxial compressive strength [20]. (15) which gives C 0 = 30.2 MPa. For a representative friction angle of shale, ϕ = 20 • , the corresponding value of the cohesion c is derived from the Mohr-Coulomb relation via which gives c = 10.5 MPa [18]. These values correspond to the D-P material parameters β = 31.6 • , d = 17.7 MPa, and ψ = 0. No plastic hardening was assumed for the shale layers. The rock properties are summarized below in Table 3.

Results
The parametric studies present results for the base case scenario which corresponds to the most likely parameters and further sensitivity analysis for the uncertain parameters. The results are presented and compared for the two reservoir fault geometries. Reservoir depletion was simulated by decreasing the reservoir pressure to a final equilibrium state. Results were obtained and presented for the effective stresses, volumetric strain, yield tendency, slipping potential, and displacement field in the reservoir and surrounding zones.

Yield Factor and Plastic Deformation
A yield factor shows how close a stress state is to a yield condition even in the elastic region. The yield factor ranges between 0 and 1 in the elastic region, with the value of 1 attained when the material has reached the yielding limit. The yield factor for the Drucker-Prager criterion is written in the form, Yield factor = − tan (17) Figure 9 presents the yield factor after the reservoir is depleted down to pore pressure value of 40 MPa. The results show that the higher the coefficient of friction, the higher is the yield tendency of the rock formation (Figure 9a). The response to the surface friction is more sensitive in the normal fault geometry, especially close to the fault. Figure 9b shows that the lower the LSR, the higher the yield tendency. The normal fault geometry presents more volatile conditions in relation to change of the friction coefficient μ and even more to the LSR magnitude. The reverse fault geometry shows less variation of the yield impact between different conditions. For the reverse fault, it is also noticeable that the HW is at higher yielding risk than the FW.

Yield Factor and Plastic Deformation
A yield factor shows how close a stress state is to a yield condition even in the elastic region. The yield factor ranges between 0 and 1 in the elastic region, with the value of 1 attained when the material has reached the yielding limit. The yield factor for the Drucker-Prager criterion is written in the form, Yield factor = t − p tan β d (17) Figure 9 presents the yield factor after the reservoir is depleted down to pore pressure value of 40 MPa. The results show that the higher the coefficient of friction, the higher is the yield tendency of the rock formation (Figure 9a). The response to the surface friction is more sensitive in the normal fault geometry, especially close to the fault. Figure 9b shows that the lower the LSR, the higher the yield tendency. The normal fault geometry presents more volatile conditions in relation to change of the friction coefficient µ and even more to the LSR magnitude. The reverse fault geometry shows less variation of the yield impact between different conditions. For the reverse fault, it is also noticeable that the HW is at higher yielding risk than the FW.

Yield Factor and Plastic Deformation
A yield factor shows how close a stress state is to a yield condition even in the elastic region. The yield factor ranges between 0 and 1 in the elastic region, with the value of 1 attained when the material has reached the yielding limit. The yield factor for the Drucker-Prager criterion is written in the form, Yield factor = − tan (17) Figure 9 presents the yield factor after the reservoir is depleted down to pore pressure value of 40 MPa. The results show that the higher the coefficient of friction, the higher is the yield tendency of the rock formation (Figure 9a). The response to the surface friction is more sensitive in the normal fault geometry, especially close to the fault. Figure 9b shows that the lower the LSR, the higher the yield tendency. The normal fault geometry presents more volatile conditions in relation to change of the friction coefficient μ and even more to the LSR magnitude. The reverse fault geometry shows less variation of the yield impact between different conditions. For the reverse fault, it is also noticeable that the HW is at higher yielding risk than the FW. As it is expected, plastic deformation increases as the LSR decreases due to the higher deviatoric stress. In this analysis for low LSR, it appears that bands of plastic deformations with values of 0.1% develop in the shale formations whereas in the sandstone zone the plastic distribution is more distributed in the range of 0.03-0.06% (Figure 10). For the mid and high LSR cases no significant signs of plastic deformation developed at the depleted condition of Pp = 40 MPa. As it is expected, plastic deformation increases as the LSR decreases due to the higher deviatoric stress. In this analysis for low LSR, it appears that bands of plastic deformations with values of 0.1% develop in the shale formations whereas in the sandstone zone the plastic distribution is more distributed in the range of 0.03-0.06% (Figure 10). For the mid and high LSR cases no significant signs of plastic deformation developed at the depleted condition of Pp = 40 MPa.

Stress Path
The stress path is traditionally calculated by the coefficient of earth pressure, defined as the ratio of the change of the minimum effective horizontal stress (Δσ′3) to the change of the effective vertical stress (Δσ′1) An alternative, common parameter to describe the stress path during production is defined as the ratio of the change in horizontal stress Δσ3 to the change in pore pressure Δp [21] = ∆ ∆ (19) Similar relations hold for the other directions defining the maximum horizontal and vertical coefficients, γH and γv, respectively. The γ is also called the depletion coefficient and the γv = Δσv/Δp is called arching coefficient. A related study [22] found that a reservoir follows a different stress path between depletion and injection. Additional studies [3] and [23] also examined the effect of faults on stress path evolution during reservoir pressurization. The nonuniform stress paths during the pore pressure depletion or injection is caused mainly by the reservoir geometry, the rock layer stiffness, rock heterogeneity, and plastic deformation. Figure 11 shows the stress path from the in-situ condition to the reservoir depleted conditions (shaded area). The analytical predictions of the stress path based on different values of Ko are also shown. The results show that a low stress path can lead to rock plastic yielding and shear failure during depletion. On the other hand, a high stress path is more likely to cause reservoir compaction. Sensitivity analysis for the surface friction coefficient does not show significant deviation from the baseline whereas changes in the LSR have a greater impact.

Stress Path
The stress path is traditionally calculated by the coefficient of earth pressure, defined as the ratio of the change of the minimum effective horizontal stress (∆σ 3 ) to the change of the effective vertical stress (∆σ 1 ) An alternative, common parameter to describe the stress path during production is defined as the ratio of the change in horizontal stress ∆σ 3 to the change in pore pressure ∆p [21] Similar relations hold for the other directions defining the maximum horizontal and vertical coefficients, γ H and γ v , respectively. The γ is also called the depletion coefficient and the γ v = ∆σ v /∆p is called arching coefficient. A related study [22] found that a reservoir follows a different stress path between depletion and injection. Additional studies [3] and [23] also examined the effect of faults on stress path evolution during reservoir pressurization. The nonuniform stress paths during the pore pressure depletion or injection is caused mainly by the reservoir geometry, the rock layer stiffness, rock heterogeneity, and plastic deformation. Figure 11 shows the stress path from the in-situ condition to the reservoir depleted conditions (shaded area). The analytical predictions of the stress path based on different values of Ko are also shown. The results show that a low stress path can lead to rock plastic yielding and shear failure during depletion. On the other hand, a high stress path is more likely to cause reservoir compaction. Sensitivity analysis for the surface friction coefficient does not show significant deviation from the baseline whereas changes in the LSR have a greater impact.  Figure 12 shows the final value of the stress path along the reservoir middle horizontal line in both FW and HW blocks. The results show the sensitivities with respect to LSR and fault friction μ. The LSR scenarios were carried out for the mid value of surface friction coefficient and the friction coefficient sensitivity for the mid LSR scenario. The friction coefficient influences the stress path only near the fault. The results are more sensitive to LSR for both fault geometries. Apparently, a low LSR initial condition results in higher stress paths to reach a stress equilibrium whereas for high initial LSR equilibrium is achieved at lower stress paths. For example, for low initial LSR = 0.76 equilibrium achieved with stress path K = 0.45 and the final LSR is 0.63 whereas for high initial LSR = 1.0 geostatic equilibrium was achieved with stress path K = 0.3 and the final LSR is 0.83. The stress path decreases near the normal fault whereas it increases near the reverse fault. Figure 13 shows contours of the final stress path after depletion for low and high LSR values in NF and RF geometries. Here, again the low LSR case results in higher values of the stress path in the reservoir after depletion and in a more variable distribution and localized zones elsewhere. In all cases, the value of initial LSR is the main parameter that controls the stress path. The fault geometry and surface friction play a role only near the fault.   Figure 12 shows the final value of the stress path along the reservoir middle horizontal line in both FW and HW blocks. The results show the sensitivities with respect to LSR and fault friction µ. The LSR scenarios were carried out for the mid value of surface friction coefficient and the friction coefficient sensitivity for the mid LSR scenario. The friction coefficient influences the stress path only near the fault. The results are more sensitive to LSR for both fault geometries. Apparently, a low LSR initial condition results in higher stress paths to reach a stress equilibrium whereas for high initial LSR equilibrium is achieved at lower stress paths. For example, for low initial LSR = 0.76 equilibrium achieved with stress path K = 0.45 and the final LSR is 0.63 whereas for high initial LSR = 1.0 geostatic equilibrium was achieved with stress path K = 0.3 and the final LSR is 0.83. The stress path decreases near the normal fault whereas it increases near the reverse fault. Figure 11. Stress paths in t-p plane for different Ko and LSR condition of the mid reservoir regions in normal and reverse fault geometries. Figure 12 shows the final value of the stress path along the reservoir middle horizontal line in both FW and HW blocks. The results show the sensitivities with respect to LSR and fault friction μ. The LSR scenarios were carried out for the mid value of surface friction coefficient and the friction coefficient sensitivity for the mid LSR scenario. The friction coefficient influences the stress path only near the fault. The results are more sensitive to LSR for both fault geometries. Apparently, a low LSR initial condition results in higher stress paths to reach a stress equilibrium whereas for high initial LSR equilibrium is achieved at lower stress paths. For example, for low initial LSR = 0.76 equilibrium achieved with stress path K = 0.45 and the final LSR is 0.63 whereas for high initial LSR = 1.0 geostatic equilibrium was achieved with stress path K = 0.3 and the final LSR is 0.83. The stress path decreases near the normal fault whereas it increases near the reverse fault. Figure 13 shows contours of the final stress path after depletion for low and high LSR values in NF and RF geometries. Here, again the low LSR case results in higher values of the stress path in the reservoir after depletion and in a more variable distribution and localized zones elsewhere. In all cases, the value of initial LSR is the main parameter that controls the stress path. The fault geometry and surface friction play a role only near the fault.   Figure 13 shows contours of the final stress path after depletion for low and high LSR values in NF and RF geometries. Here, again the low LSR case results in higher values of the stress path in the reservoir after depletion and in a more variable distribution and localized zones elsewhere. In all cases, the value of initial LSR is the main parameter that controls the stress path. The fault geometry and surface friction play a role only near the fault.  Figure 14 shows contours of the volumetric strain after depletion for the cases of low and high LSR for normal and reverse faults geometries. The negative values of the volumetric plastic strain show predominantly volume decrease in the elastic region. The elastic strains still dominate over any dilatant behavior with plastic shearing during depletion. The volume changes are less homogeneous near the fault area in the case of low LSR which triggers more plastic deformation.    Figure 14 shows contours of the volumetric strain after depletion for the cases of low and high LSR for normal and reverse faults geometries. The negative values of the volumetric plastic strain show predominantly volume decrease in the elastic region. The elastic strains still dominate over any dilatant behavior with plastic shearing during depletion. The volume changes are less homogeneous near the fault area in the case of low LSR which triggers more plastic deformation.

Slip Factor
The stress equilibrium conditions across a fault surface results in continuous normal and shear stress but the normal stresses parallel to the discontinuity are not necessarily equal ( Figure 15).

Slip Factor
The stress equilibrium conditions across a fault surface results in continuous normal and shear stress but the normal stresses parallel to the discontinuity are not necessarily equal ( Figure 15). The fracture slip factor (FSF) is determined from the shear and normal stresses that are acting on the fault plane according to [16] and [13].
where σ′n and τ are determined in each side by the cartesian components of the stress tensor according to where θ, is the angle between the direction of and the direction of the normal to the inclined plane.
The FSF must be smaller than or equal to the surface friction coefficient. It is assumed here that the fault surface friction is the same as of the intact material. For the sandstone μ = 0.58 which corresponds to friction angle, φ = 30° and for the shale is μ = 0.36 for φ = 20°.
The impact of surface friction and LSR on the fracture slip factor are presented in Figure 16 where the calculations were carried out on the near fault surface of the HW block. Figure 16a,c shows the slip factor for a low and high surface friction in the normal fault and reverse fault, respectively. In the case of a low surface friction, the determined slip factor reaches the friction limit of 0.36 in the upper overburden of both NF and RF geometries, leading further to stress redistribution. In the case of a high surface friction the stresses remain locked as they are below the slip condition of 0.75, everywhere along the surface. Figure 16b,d shows the slip factor for different LSR ratios. For these calculations the friction coefficient is equal with the mid value of 0.58. It is clear that the high LSR increases the normal stress on the surface mitigating the risk of slip. For low LSR the slip factor moves towards or reaches the limit of surface friction resistance in the overburden. The fracture slip factor (FSF) is determined from the shear and normal stresses that are acting on the fault plane according to [13,16].
where σ n and τ are determined in each side by the cartesian components of the stress tensor according to where θ, is the angle between the direction of σ x and the direction of the normal to the inclined plane. The FSF must be smaller than or equal to the surface friction coefficient. It is assumed here that the fault surface friction is the same as of the intact material. For the sandstone µ = 0.58 which corresponds to friction angle, ϕ = 30 • and for the shale is µ = 0.36 for ϕ = 20 • .
The impact of surface friction and LSR on the fracture slip factor are presented in Figure 16 where the calculations were carried out on the near fault surface of the HW block. Figure 16a,c shows the slip factor for a low and high surface friction in the normal fault and reverse fault, respectively. In the case of a low surface friction, the determined slip factor reaches the friction limit of 0.36 in the upper overburden of both NF and RF geometries, leading further to stress redistribution. In the case of a high surface friction the stresses remain locked as they are below the slip condition of 0.75, everywhere along the surface. Figure 16b,d shows the slip factor for different LSR ratios. For these calculations the friction coefficient is equal with the mid value of 0.58. It is clear that the high LSR increases the normal stress on the surface mitigating the risk of slip. For low LSR the slip factor moves towards or reaches the limit of surface friction resistance in the overburden.  Figure 17 shows the developed with reservoir depletion horizontal (top) and vertical (bottom) displacements along the middle line of the reservoir. The fault is located at the horizontal distance of 2000 m where the displacement fields are discontinuous as the line crosses from the FW to the HW. The results are shown for both normal fault and reverse fault geometries. It is clear that in all cases the displacement field is dominated by the vertical components. The horizontal displacements are mainly due to fault-slip and are nearly zero except in the case of a low friction in the normal fault geometry. In this case, the horizontal positive movement at the HW and negative movement at the FW can exceed 0.06 m whereas in all other cases the displacement range is from −0.01 to 0.01 m.

Displacement
The vertical displacement is dominated by the reservoir depletion due to the increase of the vertical effective stress. The influence of the fault movement is shown by the difference in the vertical displacement between the FW and HW which is greater and discontinuous near the fault. The vertical downward displacement in the FW can exceed 0.4 m whereas in the HW there is a separation of the curves for NF and RF cases. Overall, the higher displacements develop in the near fault region. In the remote from the fault area the vertical displacement is nearly constant as it is clearly governed by the reservoir depletion.  The vertical displacement is dominated by the reservoir depletion due to the increase of the vertical effective stress. The influence of the fault movement is shown by the difference in the vertical displacement between the FW and HW which is greater and discontinuous near the fault. The vertical downward displacement in the FW can exceed 0.4 m whereas in the HW there is a separation of the curves for NF and RF cases. Overall, the higher displacements develop in the near fault region. In the remote from the fault area the vertical displacement is nearly constant as it is clearly governed by the reservoir depletion.

Conclusions
This paper examined the impact of the effective stresses that develop with depletion of a faulted reservoir. The study was carried out with finite element modeling using 2D plane strain elements with pore pressure and elastoplastic deformation of the reservoir sandstone and sealing shale layers governed by the Drucker-Prager model. The mechanical properties and functions of the rock formations were derived from model calibration on triaxial test data for the reservoir sandstones and correlation functions for the surrounding shale layers. The faulted model geometry was built using seismic data and interpretation of the field data. The main findings and conclusions are summarized as follows: The stress evolution during reservoir depletion depends on the past and current tectonic stress activity of the region. The stress path during depletion is mainly controlled by the initial lateral stress ratio (LSR).


The developed effective stresses with the depletion of the field are influenced by the fault geometry of the compartmentalized blocks.  Plastic deformation develops for low LSR and increases as the reservoir is depleted affecting mainly the region near the fault. Plastic regions are most pronounced in the shale formation, rather than in the sandstone. More volatile conditions in relation to change of the friction coefficient and LSR were found for the normal fault geometry.  Lower values of the LSR increases the slip tendency at the fault surface.  The reservoir deformation is dominated by vertical displacement. The horizontal displacements are relatively small and only in the case of the normal fault it exceeds 0.07

Conclusions
This paper examined the impact of the effective stresses that develop with depletion of a faulted reservoir. The study was carried out with finite element modeling using 2D plane strain elements with pore pressure and elastoplastic deformation of the reservoir sandstone and sealing shale layers governed by the Drucker-Prager model. The mechanical properties and functions of the rock formations were derived from model calibration on triaxial test data for the reservoir sandstones and correlation functions for the surrounding shale layers. The faulted model geometry was built using seismic data and interpretation of the field data. The main findings and conclusions are summarized as follows: The stress evolution during reservoir depletion depends on the past and current tectonic stress activity of the region. The stress path during depletion is mainly controlled by the initial lateral stress ratio (LSR).

•
The developed effective stresses with the depletion of the field are influenced by the fault geometry of the compartmentalized blocks. • Plastic deformation develops for low LSR and increases as the reservoir is depleted affecting mainly the region near the fault. Plastic regions are most pronounced in the shale formation, rather than in the sandstone. More volatile conditions in relation to change of the friction coefficient and LSR were found for the normal fault geometry. In the remote fault area, the vertical displacement is clearly governed by the reservoir depletion and it is nearly uniform.

•
The volumetric strain predominantly indicates that compaction takes place. The higher volume decreases developed for low LSR. Funding: This research received no external funding.