1 Introduction

Geomechanical characteristics are increasingly incorporated into conventional reservoir models producing advanced dynamic models for the exploration and production of hydrocarbons. The significance of geomechanics lies in its ability to provide solutions for the challenges of reservoir production and well stability. The benefits of applying geomechanics encounter in the efficient coupling of rock mechanics with reservoir fluid flow simulators to produce dynamic earth representation models (Markou and Papanastasiou 2018). The integration of rock physics with porous fluids can enhance the accuracy of reservoir models. Geomechanical models can address several issues important for well drilling planning, pore pressure prediction, well design, casing integrity, and reservoir compaction. Hydrocarbon field characterization with geomechanics provides expanded insights on integrated geological, mechanical and flow domains in oil and gas reservoir studies. Stress analysis assists in understanding the existing structural geometry of subsurface features about the applied and past tectonics. Calculation of the effective stress changes during reservoir production simulation may influence various aspects of a field development. Geertsma (1973) presented an analytical solution to predict reservoir compaction. Related numerical simulations and analyses on subsidence can be found in Baade et al. (1988), Bruno and Bovberg (1992), Chin et al. (1993), Rattia and Ali (1981).

Compartmentalized blocks in reservoirs create complicated stress–strain conditions, which become more apparent near the fault surfaces. The effective stress changes during depletion can affected by the geometry of faults and generate different stress paths in local regions at the reservoir. Reservoir depletion increases the effective stress which is transferred to the rock grain matrix of the formation (Schutjens et al. 2004). Stress analysis can be applied to identify rock failure conditions and characterize the stability of fault conditions. During reservoir depletion or injection, due to changes in the effective stresses, the faults may cause different risks related to casing shearing and wellbore collapse.

Hence, reservoir geomechanics becomes more important in reservoirs highly impacted by faults and overpressure zones, they also facilitate the optimization of reservoir management planning in terms of compaction-induced changes to permeability, well drilling and completion design, and depletion strategy (Schutjens et al. 2004). An operator has to define, as accurately as possible, the fault structural geometry of the field and assess early before production commences whether production will be affected by the presence of faults and fractures (Wayne et al. 2006), though some of those cannot be identified even at the early stages of a field production.

Plumb et al. (1998) presented a methodology to estimate the in situ stresses in a tectonic active environment and fault geometry, using finite element analysis (FEA) and well data. Reservoir depletion-induced problems which impact rock deformation and compaction have been studied in relation to effective stress changes as one of the results of reservoir depletion (Hettema et al. 2002; Papamichos et al. 2001). The application of FEA models in geological structures and reservoir scales is widely used to assess the rock behavior under stress. They can affect the integrity of overlying cap rocks, the sealing faults, and can impose constraints on the drilling mud window (Addis et al. 1996, 2001; Carbognin et al. 2000; Fjaer et al. 2008; Hettema et al. 2000). Couples and Lewis (2016) in the geomechanical models acknowledge the importance of the mean stress on volumetric strains and localization behavior of geomaterials. The deformation of geomaterials in shear zones of faults is developed in complex patterns and drives the volumetric strain changes in dilation and compaction behaviors (Papamichos 2010; Papanastasiou and Vardoulakis 1992). Stress–strain interactions and their impact on rock deformation and shearing are difficult to be examined analytically; hence, the need for FEA simulation is recommended.

Bostrøm and Skomedal (2004) investigated possible shear or tensile fracturing of formations and faults, including crushing and reactivation of the faults. They also examined the sealing capacity of the interbedded shales in a faulted system as the reservoir is depleted. A study on cap rock breach was performed by Couples and Lewis (2016) to assess when a failure can occur during a shear offset increase. It was found that the initial fault-lying seal, across the middle of the model, becomes deformed (dilates in some places) and permeability increases. Amiri et al. (2019) and Taghipour et al. (2021) analyzed the fault reactivation risk using analytical and numerical methods. Gheibi et al. (2017) examined the impacts of faults on stress path changes in CO2 underground storage formations under changing reservoir pressure conditions and evaluated the criteria in tensile and shear modes in terms of fault propagation (Gheibi et al. 2018). In addition, a linear poroelastic model was presented by Xu et al. (2012) for geological sequestration and Teatini et al. (2020) examined the safety considerations in faulted reservoirs for CO2 storage.

Popular computational frameworks in modeling reservoir deformation are based on linear elasticity or plasticity models within the framework of small deformation (Gai et al. 2003; Liu and Liu 2017; Settari and Mounts 1998; Settari and Walters 2001). The accuracy of the geomechanical models is enhanced with the valid input of rock properties, where the importance of the quality of core samples in calibrating and validating simulation models was investigated by Holt et al. (2000). FEA enables the realization of a complete state of stress and pore pressure in a complex geometry where one can assess the fault displacement, the developed strain magnitudes, and the stress and pore pressure conditions. Such modeling processes help to better understand and characterize the reservoir geometry and its evolution, while can provide additional information related to matrix quality, the migration of the fluids through active faults, and potential structural trapping under the applied mechanical stress regimes.

Papanastasiou et al. (2017) presented a geomechanical model for the Eastern Mediterranean to assess and constrain the in situ stresses developed in a tectonically active basin. Markou and Papanastasiou (2020a, b) examined an optimum depletion strategy of a compartmentalized field in the Eastern Mediterranean, using a 2D FEA model based on the plane strain deformation method, incorporating the pore pressure and elastoplastic behavior of reservoir and overburden formations. In this study, we extend the analysis to 3D model simulation to identify the broader distribution of field stresses and the impact of rock deformation in the faulted reservoir zone of the field. Using 3D finite element geomechanical models, we reconstruct the historical rock and fault blocks displacement and transformation and assess the geo-dynamics that can impact the field geological structure. In particular with such tools and methods, we can examine the equilibrium forces applied in terms of magnitude and direction to a geological structure and improve the accuracy and quality of prediction performance of the reservoir models. Furthermore, we can explore potential regions of new faults to be generated or existing faults that are not currently visible from seismic imaging.

This modeling process examines a series of depletion scenarios to account for uniform and non-uniform depletion strategies. The 3D analysis presents results for the developed effective stresses, the displacements distribution, and deformations on the formations. They also examined the development of plastic zones, the permeability distribution, and slip tendency along the fault surfaces. This work provides new insights in geomechanical impacts on reservoirs located in a deepwater and high-stress environment, accommodated by a complex reservoir structural geometry.

2 Methodology

2.1 Faults and Slip Surface

Rock formations subjected to excessive stresses lead to failure and form fractured systems by creating faults and joints. Joints develop in both tensile and compressive stress regimes. Structural geometry of faults is constructed using the interpretation of available seismic data while the 3D representation of fault blocks and their mechanical reactions can create technical modeling challenges in meshing and shear failure (Healy et al. 2015; Markou and Papanastasiou 2020b) Fault displacement and slip tendency can support the estimation of stress state in a faulted system and Alpak et al. (2019) presented a model for structural connectivity of faults in a reservoir. Markou and Papanastasiou (2018) supported that the existing structural geometry of the fault blocks in reservoirs can indicate that the deformations occurred in the reservoir through the geological time driven by paleostresses. The 3D geometry includes interaction between multiple fault blocks providing improved model visibility, so one can appreciate the evolved stress–strain impacts in a dynamic stress environment.

Faults are characterized by the coefficient of friction of the contact surface, the geometry of faults, and the applied magnitude and orientation of the principal stresses Sv, SH, and Sh. The geometry of faults is related to the dip angle of the fault slip surface, which is affected by the applied principal stresses relationship. Markou and Papanastasiou (2018) examined how the stress field related to the faulting regimes are related to field stresses, while Haimson and Rudnicki (2010) presented that the fault dip angle increases with intermediate stress. The faulting regimes are primarily classified as normal faults, strike-slipping, and reserve faults; however, in actual conditions, a combination of the primary classes exists. Also, the horizontal stress anisotropy is defined by the ratio of SH/Sh and is examined by Hettema (2020) and Markou and Papanastasiou (2020b).

The reservoir geology consists of the Early Miocene–Oligocene sandstones covered by Mid–Early Miocene mudstones and shales. Figure 1a and b shows two seismic sections from the study area, where a normal fault and a reverse fault of the Aphrodite field in the Eastern Mediterranean are illustrated (Markou and Papanastasiou 2020b). The reservoir zone is defined between the horizons in yellow and the faults are delineated with the black lines. The reservoir zone consists of interbedded sands, and it is bounded by shale layers. A geological section is shown in Fig. 1c, as part of the 1D Mechanical Earth Model (MEM) constructed by Markou and Papanastasiou (2018) for the Aphrodite field, where it presents the complexity of the faulted blocks of the compartmentalized reservoir.

Fig. 1
figure 1

Seismic profiles illustrate two different fault geometries from the Aphrodite field where a normal fault is shown in (a), and a reverse fault is shown in (b) from Markou and Papanastasiou (2020b). Reservoir section of the Aphrodite field (c) after Markou and Papanastasiou (2018)

Full size image

Faults by definition are formed when the applied shear stress exceeds the rock shear strength and causes a relative displacement along the plane of the shear. From the geostatic condition, failure can occur if approximately σv/σh = 3, where 3 is related to the friction coefficient by (1 + sinφ) /(1-sinφ) that is based on Coulomb friction law, μ = τ/σ′n (Zoback 2010). The sandstone friction coefficient μ = 0.58 for an angle of friction φ = 30°. The shale μ = 0.36 corresponds to φ = 20°. In this FEA, the fault contact properties are modeled by a tangential behavior of penalty Coulomb friction formulation with a surface-to-surface contact interaction. The numerical model considers the varying normal deformation of contact surfaces in each layer. We use the frictional behavior of penalty Coulomb with a modified allowable elastic slip method, as a generally acceptable level of model accuracy for this study. The fault contact surface is modeled as a discretized inclined plane subject to effective stresses at depth. The main components include the faults geometry, the applied reservoir pore pressure, and the estimation of principal stresses, Sv from density logs, Sh from leak-off tests, and SH from borehole data. Thus, the estimation of the normal and shear stresses on each inclined fault plane is an output of stress equilibrium.

2.2 Regional Stress Field and Pore Pressure

A comprehensive source of information on the magnitude and direction of in situ stresses is the World Stress Map which was introduced by Heidbach et al. (2016). Working on these available data and integrating the regional structural geology knowledge, we present in Fig. 2 an integrated regional structural and stress map for the Eastern Mediterranean region that illustrates the main faults and incorporates the stress indicators associated with the direction of the maximum horizontal principal stress. The main stress indicators include measurements from focal mechanisms, borehole breakouts, drilling-induced fractures, and hydraulic fractures. Together with the magnitude values can provide information on the stress orientation and regime. Our study reservoir model is located in the deep Levantine sedimentary basin, approaching the southeastern flanks of Eratosthenes Seamount. The Early Miocene–Oligocene reservoir consists of multiple high-quality sandstone intervals. As few offset calibration data were provided in the area (Markou and Papanastasiou 2018), we can consider that it lies between the strike-slip regime at the south toward the Nile delta cone and the normal faulting regime at the north, from expected detachment blocks present at the flanks of Eratosthenes.

Fig. 2
figure 2

Regional structural and stress map for the Eastern Mediterranean region illustrating stress indicators and faulting regimes (Sources integration after Heidbach et al. 2016; Montadert et al. 2014; Tari et al. 2012; Tassy et al. 2015; Tugend et al. 2019)

Full size image

The geological structural setting of the Eastern Mediterranean basin starts from the Triassic to Mid-Jurassic with passive rifting tectonic episodes creating the North African margin. Further regional extension and spreading followed until the Early Cretaceous. Then after the end of seafloor spreading in the Neotethys Ocean, compressive stresses led to regional inversion episodes and structural deformation during the Late Cretaceous where the African plate initiated the convergence to the Eurasian plate. During that period, the evolution of the Cyprus Arc system and the Syrian Arc fold belt was initiated. A period of tectonic quiescence was retained following extensive clastic depositions. Then re-compressive stresses deformed further the region during the Early Miocene, including the evolution of the Dead Sea transform fault as a result of the lateral movement of the Arabian plate to the African plate. In the Late Miocene, the Mediterranean Sea isolated from the Atlantic Ocean and rapid deposition of thick Messinian Evaporites covered the basin. In the Early Pliocene, the subduction of the Eratosthenes Continental Block below the Anatolian plate along the Cyprus Arc resulted in the uplift of the island of Cyprus. Also, Pliocene–Pleistocene siliciclastic sediments that were deposited above the Messinian Evaporites were deformed by the regional tectonic activity and salt movements. (Markou and Papanastasiou 2018; Montadert et al. 2014; Nader 2011; Symeou et al. 2018).

The regional tectonic stresses have played an important role in the Eastern Mediterranean area in forming the current complex structural regime and supporting the evolution of local uplifts and extended sedimentary basins. Also, tectonics impact resulted in the generation of structural folds and faults in the area. The studied area is located in an uplifted and faulted geological structure. The model input parameters are based on the field data from the work by Markou and Papanastasiou (2018) where calibration was performed for the reservoir pore pressure, in situ stresses, LSR, and rock properties. The reservoir fluid is dry gas and the estimated pore pressure is Pp = 61.1 MPa. The depletion pressure is modeled to drop at Pdepl = 20.0 MPa. The estimated vertical stress is Sv = 88.6 MPa at the top of the reservoir. In sedimentary rocks, the lateral constraints generate horizontal stresses different from the vertical stress (Fjaer et al. 2008) and the horizontal stress determination can be expressed by the lateral stress ratio (LSR) as the ratio of horizontal stress to vertical stress. In this case, the ratio of the horizontal effective stresses to the vertical effective stress is expressed by K0, as a function of the Poisson ratio, K0 = ν/(1−ν) (Markou and Papanastasiou 2020b). Alternatively, if rock is assumed to be in limit, equilibrium K0 can be determined from the Jacky equation, K0 = 1−sinφ. In our study, the LSR = 0.83 for the Sh and LSR = 0.85 for the SH.

The orientation of applied in situ stresses is based on the borehole imaging presented in the reference work by Markou and Papanastasiou (2020a) using FMI resistivity and Stereonet plotting. As the three-dimensional model creates a blocky geometry, the original in situ stresses Sh and SH were transformed to stresses that correspond to the current geometrical model axis by introducing the Shmodel and SHmodel, respectively. In this process, a better geometrical definition is maintained that is related to the model geometry. Thus, one can recognize that the transformation of a nearly isotropic stress field to the model axis by 35 degrees clockwise rotation slightly changes the stress magnitude (Fig. 3).

Fig. 3
figure 3

Model geometrical transformation of horizontal stress components

Full size image

2.3 Elastoplasticity

FEA models use the elastoplastic theory to provide the mechanical behavior of the rock formations subject to stress and strain impacts. The elastoplastic models are efficiently used in reservoir depletion studies to assess the rock failure state at changing reservoir pore pressure conditions. The poromechanical theoretical basis was presented in earlier studies (Biot and Willis 1957; Brown and Korringa 1975; Detournay and Cheng 1993; Zimmerman 1991). Like the Mohr–Coulomb (M–C) model, the extended von-Mises yield and failure criterion provides the Drucker–Prager (D–P) model which is appropriate for frictional materials. The typical granular-like rocks that exhibit pressure-dependent yield and failure states become stronger as the mean pressure increases. The M–C and D–P models allow for a material with shearing to harden or soften isotropically for volume changes derived from the flow rule. The inelastic volume changes are determined from the dilation angle which is the ratio of volumetric change rate to the shear strain rate. In this analysis, the linear D–P model is used and is based on the three stress invariants.

The triaxial tests were performed on reservoir sandstone core samples at different levels of confinement and the data were used to calibrate the initial yield, the ultimate strength, and softening surfaces to the failure state of samples (Markou and Papanastasiou 2020a). The M–C material parameters model are obtained by plotting the Mohr circle for states of stress at initial yield and peak strength in the normal stress–shear stress diagram. The initial yield and peak strength envelopes are determined from the tangents to the Mohr circles. Though, in this study, the D–P plasticity model is used, we chose to calibrate first a M–C material model, which is widely used, and its parameters have a straight-forward physical meaning from the friction law of Coulomb. Thus, we derived the material parameters for the M–C and then we used mathematical relations to convert them to the D–P plasticity model (Durban and Papanastasiou 1997). The shale parameters were derived from seismic velocities and correlation functions. The rock physical and mechanical properties of the sand and shales used are presented in Table 1.

Table 1 Rock properties for the model design (Markou and Papanastasiou 2022)
Full size table

2.4 Flow, Consolidation, and Compaction in Porous Media

The consolidation theory was developed in soil mechanics for predicting subsidence that develops with time when a low permeability layer such as saturated clay or sandstone is loaded by surface loads due to structures. With the load application, the initial generated overpressure dissipates with time transferring the loads on the soil skeleton increasing the effective stresses which cause volume deformation that appears as subsidence at the surface. In reservoir flow, the decrease of the pore pressure is due to hydrocarbon withdrawal with reservoir depletion. Pressure depletion causes also an increase in the effective stresses to maintain equilibrium with the overburden load. The increase of the effective stresses results in rock deformation and associated problems described in the introduction. The deformation depends on the material behavior. In the present study, any time-dependent character (creeping) of flow-deformation problem was ignored.

2.5 Model Design

The FEA model follows the basic sequence shown in Fig. 4 starting with the reservoir inputs, the finite elements model, the analysis, and the process output. The main challenges in dealing with reservoir modeling using finite elements are the mesh size, which is related to the accuracy of the model, the time consumed to run a simulation, and the complex geometry that can cause local stress anomalies in dealing with the surface contact interactions. In this work, we use the Abaqus FEA software with the support of Python programming subroutines to efficiently handle the model outputs.

Fig. 4
figure 4

Workflow diagram showing the basic steps to perform a finite elements analysis

Full size image

The model is constructed to simulate a faulted sandstone reservoir zone, encased vertically by shales. In a real environment, there is a high complexity between the permeable and non-permeable zones, especially in interbedded sand situations. Nevertheless, a reasonable simplicity for understanding the overall dynamic stress behavior of the formation is kept by considering uniform bodies with averaged rock properties, especially for the net-to-gross parameters, porosity, and permeability. Our model is calibrated based on the field pressure data (Markou and Papanastasiou 2018) and the triaxial tests conducted on reservoir sandstone core samples (Markou and Papanastasiou 2020a). The results were validated in the initial conditions from the in situ field data, and as the field is not a producer yet, the modeled depletion conditions were validated through the analytical method of Hooke’s law and the outputs of triaxial tests. Also, the results in Table 3 can validate the model simulation field stress conditions. The model shown in Fig. 5 considers the three fault blocks that include a natural gas reservoir formation. The vertical fault throw reaches 100 m on both normal and reverse fault surfaces. The dip angle for the normal fault is 57° and for the reverse fault is 66°. In the discretized model geometry, the formation zones are illustrated in grey for shales and in yellow for sandstone. The model calibration, verification, and quality control were performed at observation points in key locations in the wells to ensure that the obtained results were consistent with the expected results. Three scenarios assessing different depletion schemes are described below.

Fig. 5
figure 5

3D model geometry showing the finite element meshing (a). Section illustrating the applied loads and boundary conditions (b) and schematic diagram of tetrahedron element (c)

Full size image

In this model, we use continuum pore pressure tetrahedron elements of 10 nodes (Fig. 5c) coupled with pore fluid diffusion-stress analysis that are suitable to mesh complex geometries. Grid dimensions are 4000 × 4000 × 900 m. The discretized domain includes the coarse elements with a size of approximately 200 m and the fine mesh elements at the slip regions with a size of approximately 50 m to capture steep gradients on the edges. In the FEA model, the actual principal stresses act in the form of effective stresses σʹv, σʹH, and σʹh. The vertical stress is applied at the top of the model, while the lateral sides and basement are constrained with roller boundary conditions. The lateral stress is set at the initial conditions defined as tectonic stress or residual horizontal stress.

The selected sampling locations at defined pseudo wells A and B are shown in the discretized model in Fig. 5a including also the fault slip surfaces. A section horizontal profile is also constructed through wells A and B to examine the lateral stress impact along the model (Fig. 5b). The fault blocks maintain uniform pressure equilibrium in scenario S1, while each block changes to isolated pressure conditions at the non-uniform depletion scenarios S2 and S3. The deformation model is based on the Terzaghi effective stresses derived from the equilibrium total stresses and pore pressure (Markou and Papanastasiou 2020a). The analysis examines three depletion scenarios as shown in Table 2 and compares the uniform and non-uniform depletion options. For example, the second option refers to production from Block A and assesses the impacts locally and at Block B. Our simulation focuses on the stress–strain behavior at different depletion stages. The linear relation considered is a simplified form to simulate the pore pressure conditions in the system. It provides the same stress–strain results at each equilibrium stage. Figure 6 presents the pressure depletion scenario S2 where Well A is active and Well B serves as an observer. This scenario refers to a non-uniform depletion state with the depletion to take place in Block A only.

Table 2 Design depletion scenarios presenting the geostatic stress state at the reference depth in the center fault block location
Full size table
Fig. 6
figure 6

Pore pressure profile plot during production of Well A for scenario S2. Time is presented in geostatic steps

Full size image

3 Results and Discussion

Three depletion scenarios were conducted through the simulation models. In each scenario, the depletion impact was tested under different pore pressure conditions. The objective of these multiple simulations is to understand the progressive stress distribution in the 3D space along the depleted reservoir zones, to assess the fault reactivation risk, and to analyze how depletion can impact the deformation and plasticity development in the reservoir.

Table 3 presents the main output for the well locations A and B through the different depletion scenarios. The depth measurement reference is at the mid-reservoir zone for the fault Blocks A and B. The reservoir at fault Block B is 100 m lower than A in a normal fault geometry scheme.

Table 3 Results of finite elements analysis output for Well A and B locations. Comparison between the initial and end of simulation model depletion state
Full size table

3.1 Effective Stresses

The effective vertical stress, σʹv, profile is a key parameter in pore pressure prediction studies and wellbore stability design. High-stress differences at the fault interface can lead to fault reactivation and as a consequence to wellbore instability issues.

Figure 7 presents results for the effective vertical stress, σʹv parameter, and its distribution at year 10 for the three different depletion scenarios. The results present how the geometry and rock properties can affect the distribution of effective vertical stress across the fault plane. The regions near to faults are subject to stress deviations from the average stress field, and we can see the change in stress distribution in each depletion scenario.

Fig. 7
figure 7

Effective vertical stress, σʹv distribution for uniform depletion in S1 (a) and non-uniform depletion of fault Block A in S2 (b) and Block B in S3 (c). The stress magnitude scale represents compressive stresses with negative sign in Pascal units (color limits: 0–100 MPa)

Full size image

Figure 8a presents the stress path which is defined as the ratio Δσʹh/Δσʹv and its evolution through the depletion. We can see that in the uniform depletion scenario (S1), the stress path deviates from its initial trend and increases from 0.32 to 0.45. This deviation is expected to occur due to plastic development in the rock at time step 3. Figure 8b presents the mean effective stresses at Wells A and B for the non-uniform depletion scenario (S2). In this scenario, we can see the gradual increase of stress difference, in terms of the mean stress, between two blocks considered without pressure communication. In such a case, the impact on the effective stress is greater in the depleted block.

Fig. 8
figure 8

Stress path evolution in relation to a time depletion for scenarios S1, S2, and S3 and b effective mean stress at the non-uniform depletion scenario S2 for Wells A and B. The fault Block A is the depleted block

Full size image

In Fig. 9, we present the results of the horizontal stress profile in a section at a depth of 5300 m along the Blocks A and B, and we compare the uniform S1 (left) and non-uniform depletion S2 (right) scenarios in terms of effective stress distribution. The simulation time corresponds to the year 10 depletion time, where the blocks have been depleted to a pore pressure of 20 MPa. Scenario S3 is not presented in that figure, as it demonstrates a similar behavior with the reverse shape of S2. In Fig. 9a, stress anomaly was developed in the fault surface region, as fault contacts can create small regions of high-stress concentration. Such regions are prone to local stress redistribution and local stress field rotation following the fault geometry. Figure 9b shows a stress step exceeding 40 MPa for σʹv and corresponding steps for horizontal stresses exceeding 20 MPa. A stress transition zone extends to 500 m from the fault contact on each side.

Fig. 9
figure 9

Effective stresses distribution for σv, σΗ, and σh for a scenario S1 in uniform depletion and b scenario S2 in non-uniform depletion. The fault contact surface is located at a distance of 2000 m

Full size image

The mean effective stress σʹmean follows a different distribution between the uniform and non-uniform depletion conditions with the highest values in the depleted block (Fig. 10). Compaction driven by an increase in the mean effective stress will mainly take place in the depleted block. The stress anomaly is shown in the fault region for the uniform depletion scenario S1, and a stress step of 25 MPa is identified in the non-uniform depletion scenarios S2 and S3.

Fig. 10
figure 10

Comparison of mean effective stress profiles for scenarios S1, S2, and S3

Full size image

3.2 Plastic Yielding

The equivalent or effective plastic strain is a scalar variable used to represent the material inelastic deformation. Plastic deformation generates a permanent change in the shape of rock before leading to fracture and failure in advanced stages.

Plastic strains develop in the reservoir in wider plastic zones in the sandstone bottom layer and localized areas along the fault zones. Plasticity develops gradually along the reservoir zones in the uniform depletion scenario. A key finding in this analysis is that the plastic strains develop in the adjacent fault block, which is not or partially depleted, to the contrary of the depleted block which remains in elastic deformation. Plastic yielding mainly develops near the fault slip zones at a narrow extent whereas a low risk of plastic behavior appears in the main reservoir body.

Accounts of effective plastic strain (PEEQ) are shown in Fig. 11 for the uniform depletion scenario. The central reservoir block regions yield plastically after time step = 3 (year 6) when the pore pressure is depleted to 40 MPa (Fig. 11b). For non-uniform depletion scenarios, no plastic deformation develops at mid-reservoir zones. Instead, a plastic transfer impact is developing at the near-fault-adjusted regions of the non-producing blocks. Plastic deformation appears also at the upper shales cap zones, as illustrated in Fig. 12 which shows a section profile along the A/B surface of a normal fault geometry. These plastic zones in the non-depleted Block A are due to stress transfer impact caused by the depletion of a nearby pressure isolated fault Block B. The plastic region expands normal to the fault surface at a distance of 200 m. Figure 13 presents the plastic regions at the surface of Block C which is attached to corresponding Blocks A and B surfaces. The results clearly show the plastic distribution behavior in the 3D space governed by the geometry of faults and rock properties.

Fig. 11
figure 11

The development of effective plastic strain (PEEQ) for uniform depletion scenario S1 a in relation to depletion time and b in relation to pore pressure change

Full size image
Fig. 12
figure 12

Stress transfer impact from depleted Block B to Block A. The plastic region expands normal to the fault surface at a distance of 200 m. Plastic region develops also at the upper shale cap rock near the upper fault contact

Full size image
Fig. 13
figure 13

Plastic development along the fault Block C surface where the left side attaches to Block A and the right side to Block B

Full size image

Overall, the plastic yielding is larger in the case of uniform depletion scenario S1 in the sand reservoir zone at the end depletion time. In all scenarios, the plastic zone is restricted to the near-fault zone at approximately 200 m. The distribution of plasticity is slightly different for each non-uniform depletion scenarios S2 and S3, demonstrating the stress transfer impact (Fig. 14).

Fig. 14
figure 14

Plastic yielding profile a along Blocks A and B and b detailed zoom to fault contact surface showing the extent of 200 m of plastic impact to the near-fault zone

Full size image

3.3 Slip Factor

The fault slip factor (FSF) can draw the tendency to induce shear slip on fractures. FSF can be calculated right on the fault-inclined surface; hence; it depends on the block geometry. The FSF is estimated using the shear and normal stresses acting at the fault plane, as previously conducted the appropriate 3D stress transformation (Fjaer et al. 2008; Markou and Papanastasiou 2018, 2020a). The stress equilibrium conditions applied in a fault surface are presented in Fig. 15.

Fig. 15
figure 15

Stress equilibrium condition at the fault plane between footwall FW and hanging-wall HW (After Markou and Papanastasiou 2020b)

Full size image

The FSF is determined according to Haddad and Eichhubl (2020) and Hettema (2020) as:

$${\text{FSF}}=\frac{\tau }{{\mathrm{\sigma^{\prime}}}_{n}}$$
(1)

where σʹn and τ are calculated using the stress tensor components according to

$${\sigma {\prime}}_{n}=\frac{1}{2}\left({\sigma^ {\prime}}_{y}+{\sigma^ {\prime}}_{{\text{x}}}\right)-\frac{1}{2}\left({\sigma^ {\prime}}_{y}-{\sigma^ {\prime}}_{{\text{x}}}\right){\text{cos}}2\theta -{\sigma }_{{\text{yx}}} sin2\theta$$
(2)
$$\tau =\frac{1}{2}\left({\sigma^ {\prime}}_{y}-{\sigma^ {\prime}}_{{\text{x}}}\right){\text{sin}}2\theta + {\sigma }_{{\text{yx}}} cos2\theta$$
(3)

The angle, θ is calculated between σx and the normal to the inclined plane.

The coefficient of friction for sandstone is μ = 0.58 and corresponds to the friction angle, φ = 30°. The coefficient of friction for shale is μ = 0.36 and corresponds to φ = 20°.

An increase of the normal stress at the surface contact can reduce the risk of slip. Its limit value depends on the surface friction, and it leads to further stress redistribution, once it is reached.

In this study, we calculated the delta pressure at the fault surfaces that is required for the fault surface to slip. A fault surface can be split into smaller sections that cross different rock properties along the vertical distribution of layers. Thus, the average values for each section can be used to improve the fault surface modeling process (i.e., the average sand/shale properties give μ = 0.47).

The equation becomes:

$$\text{Delta pressure to slip}= \frac{\tau }{\mu }-{\sigma^ {\prime}}_{n}$$
(4)

The fault slip surfaces in Fig. 16 illustrate the bounding surface for Blocks A/B. The calculated values represent the near slip surface region, as the fine mesh cells in slip zones are approximately 50 m in size. Therefore, each near slip surface side in a fault presents different values, while at the exact slip contact, the values converge. As it is mentioned, the friction coefficient μ for sandstone is 0.58 and for shale is 0.36. This parameter indicates the limit for which the two attached surfaces that drives them to slip. However, in the real world, there are multiple and different rock type surface contacts, that make slip interaction more complicated. Such cases are the sand/shale interbedded reservoirs.

Fig. 16
figure 16

Slip factor assessment of A–B surfaces for depleted cases of S1, S2, and S3 scenarios in terms of delta pressure to slip. The corresponding horizontal distance for each surface is shown in the horizontal axis

Full size image

The results for blocks A/B slip surface in depletion year 10 indicate that the sandstone zone does not present some risk of slip as the calculated values are below the sandstone limit of 0.58, which corresponds to positive delta pressure to slip. In particular, regions in blue (σʹn > 20 MPa) present no risk for slip tendency, while regions in red may indicate conditions of approaching the slip limit. The shale zones indicate a tendency to exceed the slip limit of 0.36 especially in the uniform depletion conditions. Also, some tendency for slip movement is expected at the edges of surfaces in multiple fault contacts. From the results appears that in the normal fault geometry, the footwall block shows a higher slipping tendency than the hanging-wall block, where the slipping tendency is limited to the shale covered regions.

Overall, there is variation in slip tendency along the near-fault surface contacts, but in all cases, the FSF is below the critical stress equilibrium limits, suggesting a low risk of displacement in the A/B surface. However, caution should be taken on the thin interbedded regions of fault contacts that any potential local slip displacement can drive unexpected leaks of contained fluids from one fault block to another. The slip conditions at the fault contact surfaces become complex with local areas close to fault connections being more pronounced to slip creating more unstable conditions and leading to the creation of localized areas of smaller faulted zones. Slip factor determination can potentially predict faults not recognized by the seismic resolution during seismic imaging mapping.

3.4 Deformation

In Markou and Papanastasiou (2020a, b) studies, it was found that in isolated depletion block schemes, different stress and displacement conditions evolve in the reservoir. In particular, in normal fault geometry, the footwall depletion causes the rock movement perpendicular to the fault surface, suggesting that the risk of fault activation is reduced. On the contrary, in the case of hanging-wall depletion, there is a differential displacement parallel to the fault surface indicating a risk of fault activation movement (Markou and Papanastasiou 2020a). The most predictable displacement conditions can exist in the uniform depletion scenario.

Furthermore, 3D stress distribution can naturally vary along the reservoir bodies, where FEA studies can assist in predicting rock deformation tendency, or even identify potential sub-fault segment zones. During the depletion and compaction, the volumetric strain is reduced, and the rock shrinks including the pore throats.

The volumetric strain provides a measure of rock dilation or compaction. It can be calculated by the first invariant of the strain tensor, with negative values corresponding to shrinkage (compaction) and positive values to expansion (dilation),

$${\varepsilon }_{v}={\varepsilon }_{1}+{\varepsilon }_{2}+{\varepsilon }_{3}$$
(5)

where εv is the volumetric strain and εi is the three principal strain components. The volumetric strain can be connected to changes in porosity and hence to permeability. This volumetric property is not important in elasticity but is a strong function of the dilation coefficient in the M–C and D–P model.

The section profiles of Fig. 17 present the different volumetric strain distributions between the uniform depletion S1 scenario (top), and the non-uniform depletion in S2 and S3 scenarios (bottom). This is a depletion-driven parameter that results in volumetric strain decrease with pressure depletion. Local anomalies can be found close to the fault zones and in the adjacent pressure isolated regions.

Fig. 17
figure 17

Section profiles of volumetric strain for scenarios S1 (a), S2, and S3 (b and c). Color limits − 0.30 to 0.15% where rock shrinkage (compression) is illustrated as negative values in green to red. The local areas in blue represent marginal rock expansion with positive values

Full size image

Figure 17a shows that uniform depletion results in nearly uniform compaction with εv = − 0.25% with slight variations in the fault region. In the non-uniform depletion scenarios S2 and S3 (Fig. 17b and c), compaction reaches εv = − 0.28% in the depleted block. The results show greater volumetric shrinkage of rock in non-uniform depletion compared to uniform depletion conditions. However, in non-uniform depletion, some regions at the adjacent idle blocks may present marginal volumetric local expansion anomalies that can become positive by reaching εv = 0.01%, as a consequence of overall geobodies deformation disturbance.

Figure 18 presents a comparison of volumetric strain changes along the reservoir between the examined depletion scenarios.

Fig. 18
figure 18

Comparison of volumetric strain change for scenarios S1, S2, and S3

Full size image

3.5 Displacement

Figure 19 shows the displacement fields for uniform depletion conditions. As it is expected, the displacement vector is dominated by the vertical component. The horizontal displacements are mainly driven by the fault slip and are developed in the near-fault surfaces. The horizontal displacement components propagate parallel to the corresponding fault surfaces and tend to be directed perpendicular to that surface (i.e., U1 with reverse fault along the z-axis and U3 with normal fault along the x-axis). Displacement magnitudes are controlled by the structural boundary conditions and the geometrical shape of each fault block compartment.

Fig. 19
figure 19

Displacement distribution for U2 vertical component (a). U1 represents x-axis movement (b) and U3 represents z-axis movement (c) for uniform depletion scenario S1 (color scales − 2.5 to 0 m for U2 and − 0.05 to 0.05 m for U1 and U3)

Full size image

By investigating the relative horizontal displacement patterns, we can explore potential regions of new faults to be generated or existing faults that are not currently visible from seismic imaging.

Figure 20 shows the displacement components in section profiles for the different depletion scenarios for the vertical U2 (top) and horizontal U1 (bottom left) and U3 (bottom right) components. The average sampling depth is 5300 along the fault Blocks A and B, and Block C profiles separately. The largest vertical displacement, U2, was found for uniform depletion with higher values of approximately 65 cm in the hanging block. The vertical displacement U2 presents − 15 cm displacement difference between Blocks A and B in the uniform depletion scenario S1. In non-uniform depletion scenarios, local subduction of 50 cm occurs in the depleted blocks. Transition zones of transfer impact also are recognized in the displacement fields.

Fig. 20
figure 20

Section profiles of displacements in Blocks A and B for a vertical U2, b horizontal U1, and c U3 components

Full size image

The horizontal displacement components U1 and U3 are shown in Fig. 20 (bottom). In uniform depletion conditions, a uniform behavior at the mid-depth of the reservoir at 5350 m was found, with horizontal displacements less than ± 2 cm. On the other hand, in non-uniform conditions, the horizontal displacement U1 becomes greater to reach up to + 30 cm in the middle of the reservoir, whereas the displacement U3 reaches up to + 60 cm close to the fault contact zone.

The U3 direction suggests that the region close to 500 m from the fault surface A/B may present potential risk of unstable conditions. This suggests avoiding locating a well close to the fault and justifying the need for proper modeling of geological structures of field reservoirs.

Overall, higher displacements are developing in the areas close to the fault whereas in areas away from the faults, the vertical displacement is the dominant parameter which is clearly governed by the reservoir depletion. Overall, the reservoir fault blocks and geometry can impact the displacement components.

3.6 Permeability

In this expanded 3D analysis, permeability calculation is based on previous work where a relationship between the normalized permeability and volumetric strain was established (Markou and Papanastasiou 2020a, b). Volumetric strain is a key parameter to correlate permeability curves derived from experimental rock data. The results show that the reservoir is highly impacted by the elastic strains and dilatant behavior with plastic shearing is mainly presented during the reservoir depletion at the near-fault areas.

The derived permeability plots in Fig. 21 present the normalized permeability factor over the sandstone region given by a linear relationship (Markou and Papanastasiou 2020a):

$$k=59\times {\varepsilon }_{vol}+1$$
(6)

where εvol is the difference between the model in situ geostatic state and the depleted geostatic state of the simulation. It was derived using the applied mean effective stress, calibrated from experimental data. Also, we assume isotropic permeability, but anisotropy of the permeability tensor can be also developed from deformation (Bagheri and Settari 2008). Permeability anisotropy can impact the producibility of the reservoir by creating unwanted water channelized paths to the wellbore. Alternatively, proper prediction of permeability heterogeneity can help to improve the field production planning.

Fig. 21
figure 21

Relationship of normalized permeability to a volumetric strain change and b relationship to reservoir pore pressure

Full size image

In this study, we identified that the slope of Eq. 6 can be further modified in future work as shown in dashed lines in Fig. 21a to cover different types of sandstone rocks with higher permeability moving upwards and lower permeability moving downwards. The new modified relations can become:

$$k={\left(59\times {\varepsilon }_{vol}+1\right)}^{n}$$
(7)

With the introduced exponent n, we can calibrate the initial relationship using the proposed experimental values for low permeability n = 2.8 and high permeability n = 0.59.

Figure 21 presents (a) the permeability change to volumetric strain relationship and (b) the permeability change to pore pressure reduction during depletion at Wells A and B locations. The results indicate that normalized permeability can be reduced at the model well locations up to 0.85 at year 10. The results of the 3D model confirm our previous studies based on 2D plane strain analysis for the range of permeability change (Markou and Papanastasiou 2020a).

The permeability reduction impact is expected during the depletion of the reservoir. The rate of decrease depends on the rock properties and the applied stresses. It also depends on the depletion scheme of the reservoir. The section profiles in Fig. 22 present the calculated normalized permeability distribution for scenarios S1, S2, and S3 along the Blocks A and B where the regions of higher permeability reduction are shown in dark grey. In the uniform depleted scenario S1, permeability changes’ distribution is uniform in the central reservoir zone, whereas there is a decrease near the fault region by creating a further local anomaly along the near-fault surface area. In the non-uniform depletion scenarios S2 and S3 there is a non-uniform decrease of permeability with smaller variations near the fault regions.

Fig. 22
figure 22

Section profiles of normalized permeability for scenarios a S1, b S2, and c S3. Color scale 0.8 to 1.1, with low permeability values in red

Full size image

Figure 23 shows the results in a line diagram of permeability changes when the reservoir reaches 20 MPa. For the uniform depletion scenario S1, the permeability reduces to 0.85 of its initial capacity. In the non-uniform depletion scenarios S2 and S3, the permeability reduces to 0.84 in the depleted block. In the non-uniform depletion conditions, the analysis suggests that a minimal permeability increase in the order of 0.01 units can occur to the adjacent pressure idle block due to the stress–strain transfer impact from the active depleted block.

Fig. 23
figure 23

Normalized permeability along the Blocks A and B for scenarios S1, S2, and S3

Full size image

The rock permeability decreases in each reservoir fault block that is depleted and gradually affects the adjacent blocks. Overall, different depletion scenarios drive the heterogeneity distribution patterns for permeability along the faulted reservoir zones, while in our examined reservoir, as it is of high permeability sandstones, it is derived that estimated stress-related reduction of permeability will not significantly affect the producibility of reservoir rocks.

4 Conclusion

We examined the stress-related impacts on rock deformation and permeability changes in a faulted and compartmentalized reservoir using a three-dimensional (3D) FEA. The model of reservoir depletion considers geomechanical effects by incorporating soil consolidation, plasticity mechanisms, and fault surface contact behavior calibrated to triaxial test data and correlation data and functions.

An integrated regional structural and stress map for the Eastern Mediterranean is also incorporated in this work that illustrates the main faults and stress indicators associated with the maximum horizontal principal stress in the region.

The current analysis in local reservoir scale determines the distribution of effective stresses and the development of plastic yielding in the field under different production scenarios. The volumetric strain was used to determine the permeability changes’ distribution of the reservoir formation under the developing stress conditions. The impact of effective stress in reservoir rocks increases as the field is depleted and can vary for different well depletion strategies.

Displacement magnitudes are controlled by structural boundary conditions and the geometrical shape of each fault block compartment. The results show the transfer impact in the horizontal stress field to an idle block from the depletion of an adjacent block. Also, local stress anomalies along the fault are prone to stress redistribution and rotation. High-stress differences at the fault interface can lead to fault reactivation and as a consequence to wellbore instability issues.

Surface slip conditions are affected by the rock properties and fault geometry. There is variation in slip tendency along the near-fault surface contacts, but in all cases, the slip factor is below the critical stress limits, providing a low risk of displacement in the fault surface. However, caution should be taken on the thin interbedded regions of fault contacts that any potential local slip displacement can drive unexpected leaks of contained fluids from one fault block to another. We found that the slip risk is higher in shales compared to sandstone, as they have lower slip limits.

Rock deformation is a depletion-driven parameter resulting in volumetric strain decreasing with pressure depletion. Local anomalies can be found close to the fault zones and in adjacent pressure isolated regions due to the transfer impact effect. Different depletion scenarios drive the heterogeneity distribution patterns for permeability along the faulted reservoir zones, but the estimated stress-related reduction of permeability will not significantly affect the producibility of reservoir rocks.

The results indicate the added value gained in reservoir studies using geomechanical FEA for screening reservoir depletion strategies for increasing the productivity and reducing potential rock-related instability risks.