Regular Article
Petroleum geomechanics modelling in the Eastern Mediterranean basin: analysis and application of fault stress mechanics
^{1}
Department of Civil and Environmental Engineering, University of Cyprus, Nicosia, Cyprus
^{2}
Hydrocarbons Service, Ministry of Energy, Commerce, Industry and Tourism, Nicosia, Cyprus
^{*} Corresponding author: nmarkou@mcit.gov.cy
Received:
9
April
2018
Accepted:
9
July
2018
A fault stress analysis of a typical gas field in the Eastern Mediterranean is presented. The objective of this study is to provide estimates of the in situ stresses and pore pressure for populating a regional Mechanical Earth Model and to characterize the stability of faults under current and changing reservoir conditions. The fault stability analysis is based on the MohrCoulomb frictional faulting theory. The vertical in situ stress is estimated using seismic and density data and the bounds of the horizontal stresses were determined for different fault regimes. The pore pressure for determining the effective in situ stresses is estimated using the Bowers pore pressure prediction method. Fault stress analysis is performed in a series of calculations and the results are plotted on Mohr diagrams for shear failure. The fault stress analysis is performed on a wide range of alternative azimuth orientations for S_{Hmax} in order to capture the uncertainty on the actual orientation. Sensitivity with respect to reservoir pore pressure change suggests that pressure reduction in the reservoir improves the fault stress stability, ignoring in the current analysis any stress arching effects. Pore pressure increase decreases the normal stress on the fault leading to increasing risk of shear failure of the critically stressed faults. The case study examines eight faults on the Aphrodite gas field with the objective to characterize if the faults are active or remain dormant under current stress conditions and how the stability may change in reservoir injection or depletion conditions.
© N. Markou and P. Papanastasiou, published by IFP Energies nouvelles, 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Hydrocarbon reservoirs are often affected by faults that form compartmentalized zones, which can be pressure isolated or continue to maintain interconnectivity. During reservoir depletion or injection the existence of faults may pose different risks such as wellbore shearing and collapse due to changes in the effective stresses. These stress changes in the presence of fault blocks are not uniform in the reservoir, where areasin the middle, flanks and near the faults may follow different stress paths. Full account of all these effects requires numerical coupling of the reservoir simulator with geomechanical stress analysis.
Geomechanics integration into basin geology and reservoir modelling is an emerging application. It is based on the long history of rock mechanics which is continuously enhanced with new measurements, methods and computational tools. These advances couple efficiently the reservoir fluid mechanics with geomechanics and produce detailed dynamic earth representation models for the petroleum industry. Petroleum geomechanics are more important in reservoirs that are often affected by fault mechanics that form compartmentalized zones, which can be pressure isolated or continue to maintain interconnectivity. Extensive accounts on the importance of reservoir geomechanics can be found in the classical books of Fjaer et al. (2008) and Zoback (2010).
In one of the first studies that accounts for geomechanical effects, Geertsma (1973) used an analytical solution for displacement and stress changes in a diskshaped reservoir assuming that the reservoir and its surrounding formation and overburden have the same elastic properties. Addis et al. (1996) proposed relationships for estimating the magnitude of horizontal stresses in active faulted sedimentary basins that dependent on the internal friction angle of the rock formation. These relationships assume that the rock mass is on limit equilibrium condition. Plumb et al. (1998) presented a methodology based on wellbore information and finite element analysis to estimate the in situ stresses in active tectonic settings characterized by a major fault and Papanastasiou et al. (2017) presented an attempt to constrain the in situ stresses in a tectonically active offshore basin in Eastern Mediterranean. More recently, Gheibi et al. (2017) considered the effects of faults on stress path evolutions during pressurization in the context of CO_{2} geological storage and examined the criteria for further fault propagation in tensile and shear modes (Gheibi et al., 2018).
In this study, a model application is established for characterizing the faults in a given reservoir field area in the Eastern Mediterranean that are inactive or critically stressed. Faults that tend to exceed the MohrCoulomb shear failure envelope can be considered as potential active events, or currently on a slippage mode. The Eastern Mediterranean basin evolved through the Triassic to MidJurassic passive rifting tectonic episode, which followed by further regional extension and spreading, with periodic eustatic sealevel lowering until Early Cretaceous. Then, compressive stresses affected the area since the Late Cretaceous with plate convergence of Africa to Eurasia. That compressional tectonic regime resulted in regional structural deformations and formation of the Cyprus Arc thrust belt and the Syrian Arc. This was followed by a period of tectonic quiescence and subsequent extensive clastic depositions. During the Early Miocene, the second Syrian Arc compressional episode was taken place which further deformed the area with folds, reverse faults and strikeslip faults. At Late Miocene, the basin was covered rapidly by a thick section of Messinian Evaporites as a result of the isolation of Mediterranean Sea from the Atlantic Ocean. The rest overlying interval since present day consists of PliocenePleistocene siliciclastic sediments deformed by salt movements and regional tectonics (Nader, 2011; Montadert et al., 2014). Figure 1 shows a field example from a seismic section which highlights the interpreted faults on a structural analogue in the Levantine Basin by Montadert et al. (2010).
Fig. 1 Interpreted faults on a structural analogue to Aphrodite and Tamar formations in the Levantine basin offshore Cyprus (Montadert et al., 2010). 
This case study analyses the fault stress mechanics for the Aphrodite gas reservoir in offshore Cyprus, which was discovered in December 2011 by Noble Energy. The discovery well was encountered approximately 94 m of net natural gas pay in multiple highquality sandstone intervals. The discovery well was drilled to a depth of 5860 m in water depth of about 1689 m (Noble Energy, 2011). The Aphrodite sands deposited in a deep marine setting where the sediments were transported by gravity assisted mass turbidity flows. Interbedded sands and shales forming the reservoir zones which are deposited potential on a continental slope channel or fan system.
Section 2 presents a methodology for determining the in situ stresses and pore pressure which are used in assessing the current stability of faults and in changing conditions, e.g. due to reservoir depletion or injection. Section 3 presents the theoretical background for the fault regime characterization related to stress regime and the mathematical tool for determining the acting in situ stresses on fault surfaces at different orientations. Section 4 presents existing regional information on the stress field in the East Mediterranean region and field interpretation data at the scale of the Aphrodite field. The fault stress analysis is presented in Section 5 for different stress regimes and with sensitivity analysis to the fault plane orientation relative to the in situ stresses and to pore pressure changes. The final conclusions are summarized in Section 6.
2 Methodology
The model application will provide identification of critically stressed faults over the interpreted fault pattern of a typical field in the Eastern Mediterranean. The case study will examine eight faults on the Aphrodite gas field with the objective to characterize if the faults are active or remain dormant under current stress conditions and how their status will change under different reservoir conditions (i.e. reservoir depletion or injection).
Each fault surface is treated as an inclined plane under stresses incorporating the formation geopressure developed at depth. The model methodology is based on the MohrCoulomb frictional faulting theory. The modelling workflow is shown in Figure 2 and includes the faults interpretation and orientation, pore pressure estimation, a possible estimation of S_{Hmax} from borehole data, and final calculation of normal and shear stresses at each inclined fault plane. The outcome is further analyzed by performing sensitivity analysis to assess the stress state at different orientations or pore pressure conditions.
Fig. 2 Model workflow in the fault stress analysis. 
The fault stress state is directly related to the horizontal minimum and maximum principal stresses, S_{hmin} and S_{Hmax,} respectively. The common methods for determining the stress field orientation is to use information around the wellbore such as existing drilling induced fractures and formation of well breakouts (Fig. 3). Borehole ellipticity often indicates the direction of maximum stress where the borehole eccentricity is either elastically deformed or enlarged by borehole breakouts. Breakouts are formed in the borehole area of the highest circumferential (hoop) compressive stress (Haimson and Herrick, 1986; Papanastasiou and Vardoulakis, 1992; Papanastasiou and Zervos, 2004; Papamichos, 2010). Thus, the elongation should indicate the direction of the minimum horizontal stress. Drilling induced fractures appear at the borehole wall area of the lowest hoop stress, which indicates the direction of the maximum horizontal stress.
Fig. 3 Borehole diagram illustrating the areas of breakouts and drilling induced fractures. 
Structural dimensions of faults are constructed through the seismic data interpretation, where high seismic resolution can provide better visibility on the fault structural shape. Skorstad and Tveranger (2007) provide details on fault facies modelling in seismic data, that can carry information about the size of the fault envelope and consequently on the volumetric extent of the fault zone. Healy et al. (2015) present the polymodal fault patterns, where three or more sets of faults forming and slipping simultaneously can generate threedimensional strains from truly triaxial stresses. Such fault patterns can present fundamental challenge of understanding of shear failure in rocks. Morris et al. (2014) show that bed dips, faultcutoff relationships, and smalldisplacement fault patterns in the adjacent rock volume can inform about the strain and paleostress estimates. Also, Morris et al. (2016) present a methodology of using fault displacement and slip tendency to estimate stress states.
The structural shape of the fault blocks in the reservoir can provide certain understanding on structural deformations of the geobody, but caution should be given to the fact that present applied stress field usually does not coincide with the paleostresses that originated the faults or fractures over their formation. In such a case, the uncertainty on the current stress field orientation remains important and sensitivity analysis that incorporates alternative scenarios should be examined in relation to S_{Hmax} orientation.
Model studies based on parameterized stress components in threedimensions (3D) were conducted to assess the fault stress conditions. The model was first tested for the three main fault regimes and then calculations were performed to study the Aphrodite field area, offshore Cyprus. The application aims to analyze the stress state impact on faults at different angles of orientations from the maximum horizontal principal stress and dipping angles.
3 Fault regimes and in situ stress interpretation
3.1 Fault and stress regimes
The main fault regimes, as per Anderson’s classification system, are defined as normal faulting, strikeslip, and thrust (reverse) faulting (Fig. 4). The three stress orientations and magnitudes S_{v}, S_{Hmax}, and S_{hmin} describe the stress state in a geodynamic system.
Fig. 4 Main fault regimes as per Anderson’s classification (Heidbach et al., 2016). 
As there is no often way to precisely measure or predict the magnitudes of the horizontal principal stresses, S_{hmin} and S_{Hmax}, we can draw a range of potential values based on the fault field conditions. Figure 5, presents bounds on stress magnitude defined by frictional faulting theory in (a) normal faulting, (b) strikeslip and (c) reverse faulting regimes, assuming hydrostatic pore pressure. In fields with overpressure condition the S_{hmin} and S_{Hmax} window is reduced, where, as long as pore pressure approaches lithostatic, then principal stresses approach the S_{v} as the effective stress is reduced to minimum (Zoback, 2010). In normal faulting, the range of stress magnitude for S_{hmin} and S_{Hmax} is laying below S_{v}, while in reverse faulting, that range exceeds S_{v}. The range in strikeslip faulting is an intermediate state where the S_{hmin} is S_{v} and S_{Hmax} is above S_{v}. Thus, based on existing data and geological understanding of the area, this analysis can reduce the range of the magnitude of the stress S_{hmin} between the normal faulting and strikeslip.
Fig. 5 Range of stress magnitude for the main faulting regimes in an offshore field: (a) in normal faulting the critical S_{hmin} is highlighted, (b) in strikeslipping the critical S_{Hmax} is highlighted and (c) in reverse faulting the S_{Hmax} is highlighted. Modified after Zoback (2010). 
3.2 Fault plane stress determination
MohrCoulomb frictional faulting theory is applied to estimate the normal and shear stresses for each fault plane orientation in the 3D stress space. The twodimensional (2D) Mohr diagram is widely used in the fields of structural geology, seismology, engineering geology, rock mechanics and soil mechanics (Sibson, 1985; Streit and Hillis, 2002). Threedimensional Mohr diagram is also used to explain mechanism of faulting and reactivation of preexisting fault due to changes in applied stresses (Yin and Ranalli, 1992; Jolly and Sanderson, 1997; McKeagney et al., 2004).
The key terminology in faulting theory is illustrated in Figure 6a, where the fault orientation is related to the strike direction angle in relation to North and the dip plane angle in relation to the horizontal plane. Figure 6b represents an inclined fault plane in the 3D stress space which is dipping at angle θ and its strike direction is oriented at an angle α in relation to the maximum horizontal stress, S_{Hmax}.
Fig. 6 Schematic diagram illustrating: (a) the fault strike direction and dipping angle, and (b) the inclined fault plane in the 3D stress field. 
The simple form of the governing equations are based on the frictional faulting theory for normal (σ′_{ n }) and shear (τ) stresses which are represented by points on the 2D Mohr circle. In thick deep sedimentary basins where the overburden stress has a significant magnitude, and particularly in normal faulting regimes, the coordinate system is reoriented such that the yaxis is parallel to the major principal stress and the xaxis is parallel to the minor principal stress. A schematic illustration for the 2D force equilibrium bodies under the frictional faulting theory is shown in Figure 7a. The stresses σ′_{ n } and τ are acting on an inclined plane with normal vector forming a direction β relative to the minor principal stress σ_{3}, which is parallel to the xaxis (Fjaer et al., 2008).(1) (2)
Fig. 7 (a) Frictional faulting notation applied on twodimensional force equilibrium body and (b) typical Mohr circle diagram. 
The values of σ′_{ n } and τ on any inclined plane are given by a point on the Mohr circle which was defined by the major and minor principal stresses (Fig. 7b). The radius of the circle is (σ_{1} − σ_{3})/2 and its centre is at the point (σ_{1} + σ_{3})/2 on the σaxis. The coefficient of friction, μ equals to tanφ, where φ is the friction angle. The failure plane forms an angle θ = 45° + φ/2 from the direction of the minor principal stress. The C_{0} parameter shown in the diagram is the cohesion, which is assumed to be equal to zero at the fault plane.
The stress in transformation in 3D is given in Appendix A. These relations have been used by Bott (1959) for defining tectonic regimes. These equations will be used next in the fault stress analysis to assess the conditions of existing faults in the Aphrodite gas field.
4 Existing regional stress information and stress field interpretation data
4.1 Regional stress regimes
Ghalayini et al. (2014) investigated the stress regime offshore Lebanon. Figure 8 summarizes the results of that study where the faults in red are the structures that show evidence of current activity, including offshore ENEWSW dextral strikeslip faults, onshore ENEWSW latitudinal dextral strikeslip faults, and NESW sinistral (lefthanded) strikeslip faults. In the first phase of the Levant Fracture System (LFS) movement, all structures in this map were active. During the second phase of LFS in the Pliocene to present day, only the ENEWSW dextral strikeslip faults were active and might be linked to block rotations caused by the continuous sinistral movement of the LFS in Lebanon. The occurrence of the NWSE normal faults in the deep basin might not be caused by these regional geodynamics but rather to a local stressfield fluctuation affecting only the OligoceneMiocene units (Ghalayini et al., 2014).
Fig. 8 Stress tectonic regime study in the Levant basin offshore Lebanon (modified after Ghalayini et al., 2014). 
The world stress map also provides a view on the region, where the symbol directions refer to the maximum horizontal principal stress orientation (Fig. 9). The world stress map is a source that incorporates the contemporary crustal stress information along the global tectonic provinces. The location of our study area is between the strike slip regime at the south and normal faulting regime at the north. However, as the need for reliable measurements is important for analyzing and interpreting stress patterns, regional stress patterns despite they provide important information to a wider understanding in geodynamic processes for each province, they should be handled carefully when they incorporated information in a local basin context, especially when there are few calibration offset data. This view is consistent with the study of Ghalayini et al. (2014).
Fig. 9 World stress map generated using the CASMO online interface illustrating a series of stress indicators in the Eastern Mediterranean region (Heidbach et al., 2016). The inset shows the study area location on the regional seabed map after EMODnet bathymetry (EMODnet, 2018). 
The regional stress regimes trends can also be supported by data from another regional field, the Tamar field in offshore Israel, in which a series of faults affected the entire reservoir zones. These faulted structural patterns can reflect conditions of strikeslip faulting. Figure 10 illustrates the location, depth map and a crosssection of the field (Christensen and Powers, 2013).
Fig. 10 Overview of the Tamar field showing (a) location map, (b) structure map of top reservoir, and (c) structural crosssection (NW to SE) through the Tamar1 discovery well T1 (Christensen and Powers, 2013). 
4.2 Field information and data
The seismic section in Figure 11a shows the main reservoir fault compartments along the reservoir sand units and Figure 11b illustrates faults pattern in Aphrodite field. The study area covers the fault blocks #2, #3, and #4, in which eight faults were delineated and characterized. The Well A is located in fault block #2, and the Well B is located in fault block #4. The yellow shaded area covers the reservoir sand units of the field, which are interbedded with shales.
Fig. 11 (a) Seismic section and (b) depth structural map illustrating the faults pattern in the Aphrodite field offshore Cyprus (Hydrocarbons Service – MECIT, Courtesy: Noble Energy International Ltd). 
The interpreted faults geometries and their calculated parameters are summarized in Table 1 which presents the strike azimuth relative to grid north, the averaged dipping angles for each fault plane denoted by θ, and the triangular θopposite angle, β =90° − θ. Also, the relative to S_{Hmax} orientation angle α, is calculated as, α = strike azimuth – S_{Hmax} azimuth.
Fault geometries including strike orientation and dipping angles determined in the modelling process.
The sedimentary stratigraphic section, used in this modelling process, is shown in Figure 12. The constructed onedimensional (1D) earth mechanical model at wellbore location Well A location incorporates six main rock property zones: (1) seawater, (2) PlioQuaternary siliciclastic sediments, (3) Messinian evaporates, (4) LateMid Miocene mudstones and shales, (5) MidEarly Miocene mudstones and shales, and (6) Early MioceneOligocene sandstones.
Fig. 12 Sedimentary stratigraphic section along the Aphrodite field in SWNE orientation (Hydrocarbons Service – MECIT; Courtesy: Noble Energy International Ltd). 
The pore pressure and the fracture pressure are calculated based on the pore pressure prediction module which is built to accommodate alternative pore pressure prediction methods.
In this study, the pore pressure calculation is based on Bowers (1995) method for the MioceneOligocene sandstones zone. The Bowers method is an effective stress approach that employs the virgin and unloading curve relations to account for both undercompaction and fluid expansion overpressure. The effective stresses outside the velocity reversal zones are computed from the virgin curve. For the inside velocity reversal, offset well data are needed to calibrate the estimated pore pressure curve, while in the absence of well data, nearby well analogues and sensitivities can be used to establish lower and upper bounds on pressure. An analytical reference on Bower’s original equation is given in Appendix B.
The original equation to calculate the velocity virgin (normal) curve for shale can be rewritten using SI units system as(3)where, V_{p(normal)} is the velocity at the virgin (normal) curve (m s^{−1}), V_{p(water)} is the velocity of seawater (m s^{−1}) and σ′_{w} = S_{v} − P_{w} (MPa) is the effective hydrostatic stress. A and B are Bowers parameters calibrated with offset velocity vs effective stress data.
The Bowers unloading curve, which is a modification of the original equation is defined by an empirical relation which can be rewritten using SI unis as:(4)where, U is a third parameter introducing the unloading effect, and(5)
The σ′_{max} and V_{p(max)} are estimates of the effective stress and velocity at the onset of unloading.
Therefore, the transformed Bowers equation with respect to pore pressure (P_{p}) incorporating the unloading effect can be rewritten using the SI units as:(6)
Figure 13 shows the compressional velocity log, V_{p}, and the estimated normalized curves for the interpreted maximum V_{p} on the onset of unloading. In this study two potential main regions of pressure systems have been identified labeled related to normalized curves 1 and 2. The dotted normalized curves determine the projected V_{p} normal pressure trendline for each region.
Fig. 13 Determination of V_{p} maximum and σ_{max} from V_{p} log at well A (Hydrocarbons Service – MECIT; Courtesy: Noble Energy International Ltd) 
Fracture pressure is approximated using two equations based on the frictional faulting theory and coefficient of friction, μ = 0.6 for (a) normal faulting regime and (b) the frictional equation for strikeslip regime which is derived as an intermediate stress state between the normal faulting and vertical stress. These equations are given, respectively by(7) (8)
The equations (7) and (8) are produced based on Zoback (2010) reference on the in situ stress limits from the frictional strength of faults. These limits can be drawn from a Mohr diagram where for a given value of σ_{3}, the value of σ_{1} is established by the frictional strength of preexisting fault, where the Mohr circle cannot exceed the maximum frictional strength. Also, Jaeger and Cook (1979) showed that σ_{1} and σ_{3} that corresponds to critical stressed limit at μ = 0.6 can satisfy the ratio σ_{1}/σ_{3} = 3.1, which appears in the denominator in equation (7).
The calculated pore pressure curve is calibrated to selected wireline pressure points tested the reservoir formation at the wellbore A. Figure 14 shows the profiles for pore pressure at depth in the study area and the corresponding effective vertical stress. At the approximate reference depth, the hydrostatic or normal pressure is estimated at 51.1 MPa, the vertical stress at 90.7 MPa and the pore pressure is estimated between those bounds at 62.2 MPa. The fracture pressure is calculated for normal faulting regime to be equal to 71.4 MPa and for strikeslip to be equal to 81.0 MPa. Table 2 summarizes the model calculations and assumptions for the field stresses at the reference depth. In Figure 14b and Table 2 the effective stress, σ′, is calculated from the total principal stress minus the pore pressure according to the Terzaghi effective stress theory.
Fig. 14 (a) Pore pressure prediction at depth and (b) the corresponding estimated effective stress profile. 
Model assumptions, corresponding field stresses, and pore pressure at the approximate reference depth of the study sandstones zone.
5 Results and discussion
Two sets of calculations were conducted in this study, (a) the parametric calculations and (b) the fault stress analysis. For each set, two stress conditions for the normal faulting (NF) and strikeslip faulting (SS) regimes were tested. The aim of these multiple calculations is first to understand better the trendline behavior of fault stress mechanics at different field conditions and to appreciate the risk of fault reactivation for each case, and second to analyze the fault stress mechanics in the Aphrodite field. Based on the tectonic and geological understanding of the region, the stress state is expected to be between normal faulting and strikeslipping. This information is used as an assumption in this analysis.
5.1 Parametric studies
The parametric calculations evaluated first the stress state of three common fault plane geometries: normal fault, strikeslip and reverse fault, under the examined field stress regimes. These faults orientations are considered that were formed at the time of shear rock failure at its original historic stress regime, while at present time, these faults have been translated and deviated from their original orientation as paleostress does not necessary coincide with the present applied stress field. For example, an old normal fault that was formed under a normal faulting stress regime, after a post tectonic episode of compression in the area, where reverse faulting regime conditions were developed, at the same original normal fault feature its fault blocks can slip in the reversal direction due to the change in field stresses and in particular the swap between σ_{3} and σ_{1}. However, it is worthmentioning that during the tectonic evolution of fault mechanics, the intermediate stress state between the normal faulting and reverse faulting is the strikeslip.
The three common fault geometries are illustrated in the Mohr diagram to assess the location that mostly reflect each fault regime in the 3D Mohr circles area. For this set of calculations the fault regimes normal faulting, strikeslip, and reverse faulting were first tested under the determined effective stress field conditions and assumptions of Table 2. The structural fault orientation, predefined by the typical parameters, for normal faulting is α = 0° and θ = 60°, for strikeslip is α = 30° and θ = 90° and for reserve faulting is α = 90° and θ = 30°. The assumptions used for these fault regimes are shown in Table 3. The diagram in Figure 15 illustrates the three fault examples from each regime which are delineated on a horizontal top view plane.
Fig. 15 Illustration of the three typical fault regimes delineated over a top view plane. 
Model assumptions for the three typical fault regimes plotted in a 3D Mohr diagram.
A primary assumption for the fault stress analysis is the fact that faults were formed in the past, and since then were reoriented in the 3D stress space due to postfailure and regional tectonic deformations. Thus, as mentioned earlier, the present day field stress orientations are expected to differ from the paleostress regime applied in the past, millions of years ago.
The three typical stress fault regimes as calculated are plotted in the 3D Mohr space diagrams of Figure 16. Each of the examples is located into the corresponding plot area which is controlled by its geometry and effective local stress field conditions. For the normal fault stress conditions (S_{v} ≥ S_{Hmax} ≥ S_{hmin}), the normal faulting plane is located on the upper part of the large Mohr circle (σ_{1}σ_{3}), the strikeslip is located in the circle (σ_{2}σ_{3}), and the reverse faulting in the circle (σ_{1}σ_{2}). In the case of a strike slip stress conditions (S_{Hmax} ≥ S_{v} ≥ S_{hmin}), the normal faulting plane moves to the Mohr circle (σ_{2}σ_{3}), the strikeslip moves to the upper part of large circle (σ_{1}σ_{3}), and the reverse faulting remains on the circle (σ_{1}σ_{2}). Stress conditions that approach or exceed the shear failure line as shown in the plot (for μ = 0.6), can be considered as potential active events, or currently on a slippage mode.
Fig. 16 Illustration of the three common fault stress states in a 3D Mohr diagram for (a) normal faulting and (b) strikeslip stress regimes. The dashed lines in grey colour indicate the stress transformation with different fault plane orientation and dipping. 
Processing the calculated results we constructed the 2D Mohr stress nomograph as a function of the horizontal deviation α from S_{Hmax} and of the vertical inclination β (the complimentary angle to dipping θ). Figure 17 shows the calculated slip factor for the (a) normal faulting and (b) strikeslip regimes. The slip factor is calculated as the ratio of shear stress to normal effective stress (τ/σ′).The upper limit of this slip factor is μ = tanφ = 0.6. Figure 18 shows the contours of the slipfactor on the (α, θ) angle space.
Fig. 17 Mohr circle plots showing the slip factor for the modelled (a) normal faulting and (b) strikeslip stress regimes. 
Fig. 18 Contour diagrams of the slip factor as a function of α, θ angles for: (a) normal faulting and (b) strikeslip stress regimes. 
5.2 Uncertainty on the stress field orientation relative to the fault plane
One of the uncertainties that influences the fault stress analysis is the relative orientation of the stress field relative to the fault plane direction. This uncertainty arises mainly due the unknown orientation of the horizontal stresses.
The next set of calculations assesses the global stress response at field conditions considering different horizontal orientations of the horizontal maximum principal stress, S_{Hmax} and the fault dipping angles. As before, the Mohr circles are plotted for the effective stresses given in Table 3, as the stability of faults is governed by the effective stresses. This assessment establishes a nomograph on the 3D Mohr plot with parameterized values for each potential case.
The model assumptions and obtained results are shown in Figure 19 for the normal faulting regime and in Figure 20 for the strikeslip regime.
Fig. 19 Fault stress plot in normal faulting regime. Model parameters definition for different strike directions relative to S_{Hmax} and dipping angles. 
Fig. 20 Fault stress plot in strikeslipping regime. Model parameters definition for different strike directions relative to S_{Hmax} and dipping angles. 
The results are plotted on 3D Mohr diagrams for a total of 15 realizations in each regime: three strike fault angle scenarios of α = 0°, 30° and 90° in relation to maximum principal horizontal stress, S_{Hmax} and five fault dipping cases for each scenario, with dip θ = 90°, 80°, 60°, 45° and 0° (β = θ − 90°).
In the normal faulting regime plots we can see how the shear stress failure risk is reduced as the fault strike direction deviates from the S_{Hmax} orientation. On the other hand, in the scenario of strikeslip regime, the stress failure risk is increased as the fault strike direction approaches the angle α = 30° relative to S_{Hmax}. In terms of dipping angle, in normal faulting regime the failure line is approached at dip, θ = 60° (Fig. 19a), while in strikeslip regime the failure line approaches the dip, θ = 90° (Fig. 20b).
5.3 Sensitivity of fault stress stability to pore pressure changes
Another parameter which influences the stability of the faults is the magnitude of pore pressure. This influence could be due to uncertainty in estimating its value or due to regional variations in fault blocks with flow barriers. Furthermore, the pore pressure may increase during reservoir injection or decrease in reservoir depletion conditions. Differences in the pore pressure change can take place locally, e.g. between the centre of the reservoir and its flanks due to stress arching. Pressure differentiations across the field can be addressed by considering the geomechanical coupling. Furthermore, the presence of the fault blocks makes the coupling of fluidflow with geomechanics more important.
Fault sensitivity analysis to porepressure changes is performed in a series of calculations and the results are plotted on the MohrCoulomb diagram for shear failure, defined by the frictional coefficient μ = 0.6. The estimation of delta pressure to rock failure is based on the principle that the pore pressure changes affect only the normal stress and not the shear. Therefore, the stress path will shift horizontally to the left if pore pressure increases until it touches the shear line, as shown in Figure 21. The analysis is performed focussing on a wide range of alternative azimuth orientations of S_{Hmax}, for 45°, 90°, and 135°, in order to capture the uncertainty on the actual applied stress field orientation. Shear failure occurs at the point where the Mohr circle touches the shear failure line.
Fig. 21 Schematic diagram showing the estimation of delta pressure to shear failure. 
The results are summarized in Table 4 for normal faulting regime and in Table 5 for strikeslip regime. The values in red colour represent the combination of faultplane and stress field orientation under critical faulting conditions. Values in orange colour are classified as subcritical, while the rest values represent dormant faults with the minimum risk of reactivation at the given conditions. Red colour represents critical stressed cases where pore pressure increase less than 5 MPa will lead to fault activation. Orange colour represents values between 5 and 10 MPa which are classified as subcritical. The rest values suggest stable conditions.
Normal faulting conditions: Fault stress analysis results for delta pressure to shear failure in MPa.
Strikeslip conditions: Fault stress analysis results for delta pressure to shear failure in MPa.
The results show that potential S_{Hmax} orientation at 90° is the most critical for causing reactivation of faults D, G, and H in both faulting regime hypotheses. Orientation at 45° suggests that faults C, and D can be critical stressed at normal faulting regime only, whereas orientation at 135° suggests the faults A, B, E, F, G and H can become critically stressed under normal faulting regime only.
Additional sensitivity analysis examines the stress state relationship to pore pressure change by ±10%. Sensitivity results suggest that pressure reduction in the reservoir improves the stress stability in the field, ignoring in the current analysis any potential impact of stress arching in the formation. Figure 22 shows the sensitivity in pore pressure change in the reservoir, considering the S_{Hmax} azimuth at 90°. The middle plots represent the baseline case scenarios for normal faulting and strikeslipping regimes. The plots to the left simulate the effective stresses change in a pressure depletion process whereas the plots on the right simulate the changes of the effective stresses in field injection conditions. The overall trendline shows that when pressure decreases, the effective normal stress increases and the Mohr circle is shifted to the right moving away from the risk of shear slip. On the other hand, in pressure increase conditions due to injection or local nonuniform compaction, the normal effective stress increases and shifting the Mohrcircle towards the shear slip line.
Fig. 22 Influence of the pore pressure changes on the faultplane stability conditions. 
Some factors that can further reduce the uncertainty in the characterization of critical stressed faults, are the optimization in pore pressure prediction, data incorporation of pressure measurements (i.e. pressure points, mudlogging and leakoff testing), and borehole image logs analysis for determining the orientation of the horizontal stresses.
6 Conclusion
We presented a methodology for identifying the critically stressed faults of a typical gas field in the Eastern Mediterranean. The case study examined eight faults on the Aphrodite gas field with the objective to characterize if the faults are active or remain dormant under current stress conditions and how the stability may change in reservoir injection or depletion conditions.
The vertical in situ stress field was estimated using seismic and density data and the bounds of the horizontal stresses were determined for different fault regimes. The pore pressure for determining the effective in situ stresses was determined using the Bowers pore pressure prediction method.
A series of parametric studies were performed to understand the trendline behavior of fault stress state mechanics at different field conditions and to appreciate the risk of fault reactivation for each case. Fault stress analysis performed in a series of calculations and sensitivities of the results were plotted on Mohr diagram for shear failure. The analysis performed on a wide range of alternative azimuth orientations for S_{Hmax} in order to capture the uncertainty on the actual applied orientation. Sensitivity with respect to reservoir pore pressure change suggests that pressure reduction in the reservoir improves the fault stress stability, ignoring in the current analysis any stress arching effects. Pore pressure increase decreases the normal stress on the fault leading to increasing risk of shear failure of critically stressed faults.
Acknowledgments
The authors are grateful to Cyprus Ministry of Energy, Commerce, Industry and Tourism, for supporting this research and to the Hydrocarbons Service and Noble Energy International Ltd for permitting data access and publication.
Appendix A
Stress transformation
The stress analysis method is based on the stress transformations in 3D presented in Jaeger et al. (2007). The basic idea of the problem is to calculate the applied 3D stresses on an inclined fault plane and plot the results in a 3D Mohr diagram. The primary stresses are the input data, together with the zenith and longitudinal transformation angles. The model outputs include the effective normal and shear stress calculations over the inclined fault plane.
If the stress components are known in a given coordinate system, the stress components in a second coordinate system that is rotated from the first one, can be found by matrix multiplication of 3 × 3 dimensions (Jaeger et al., 2007).
First, a coordinate system (x′, y, z′) is considered that is rotated by some angle in 3D space from the original system (x, y, z) (Fig. 23).
Fig. 23 Coordinate system that utilizes the zenith and longitudinal angles (Jaeger et al., 2007). 
The direction cosines of the unit vector e_{ x′} , e_{ y′, } e_{ z′} relative to the unprimed coordinate system are: e_{ x′} = I_{11}, I_{12}, I_{13}, e_{ y′} = I_{21}, I_{22}, I_{23}, e_{z′} = I_{31}, I_{32}, I_{33}.
The traction on the plane whose outward unit normal vector is e_{ x′,} is given by p(e_{ x′} ) = τe_{ x′}, where τ is the stress tensor in the original coordinate system. The components of this traction in the e_{ x′} direction is found by taking the inner product of p(e_{ x′} ) with e_{ x′} as:(9)
The same approach is performed and for the other eight components of the stress tensor, e.g. τ_{ x′y′} = e_{ y′} p(e_{ x′}), etc. The resulting equations can be written in matrix form as following (Jaeger et al., 2007)(10)and in compact form as:(11)where, the rows of the rotation matrix L are formed from the direction cosines of the three primed unit vectors.
Considering the rotation to the zaxis, which corresponds to a rotation matrix in which I_{13} = I_{23} = I_{31} = I_{32} = 0 and I_{33} = 1, the equation (10) can be reduced to the following form:(12) (13) (14) (15) (16) (17)
Using the longitude angle, λ and the zenith angle, θ, to specify a plane, rather than the direction cosines of that plane relative to a particular Cartesian coordinate system, the axes associated with the primed coordinated systems are Pz′, which is in the radial direction, Px′, which is in the plane OPz and is associated with angle θ, and Py′ which is chosen so as to complete the righthanded coordinate system with points in the direction of increasing λ. The components of these unit vectors, relative to the unprimed coordinate system e′_{ z }, e′_{ x }, and e′_{ y } are now expressed in terms of λ and θ, and are given by:(18) (19) (20)
If the axes (x, y, z) are the principal axes, then the stress components on the plane Px′y′ can be given by:(21) (22) (23)
Appendix B
Bower’s original equation
Estimating pore pressure at depth
Direct measurement of pore pressure in relatively permeable formations is straight forward using a variety of commercially available technologies conveyed either by wireline (samplers that isolate formation pressure from annular pressure in a small area at the wellbore wall) or pipe (packers and drillstem testing tools that isolate sections intervals of a formation). Similarly, mud weight is sometimes used to estimate pore pressure in permeable formations as drilling mud must be in excess of the pore pressure and in the case that is lower formation fluids flow into the well (Zoback, 2010).
Hydrostatic pressure is estimated as following:(24)where:

P_{p}: hydrostatic pressure in Pa;

ρ: density of water in Kg m^{−3};

g: gravity acceleration constant is 9.81 m s^{−2};

Z: depth of the water column.
Similar to the above equation, the vertical stress, S_{v} (or overburden pressure) is estimated as following:(25)where ρ_{(Z)} is the bulk density of depth interval dZ, and g is the gravity acceleration.
There are many different methods of estimating the pore pressure at depth in the literature, such as the Eaton (1975), Bowers (1994), Hottmann and Johnson (1965), and Athy (1930). In the current study, the Bowers method was selected to be used having three calibration parameters A, B and U that will be discussed below.
The Bowers (1994) proposed the following relationship for interval velocity, V_{i},(26)where:

V_{i}: interval velocity in ft s^{−1};

σ′_{v}: Effective stress in psi;

A and B: Parameters calibrated with offset velocity versus effective stress data.
In addition to the above equation, the unloading curve is defined by the empirical relation as per Bower study as:(27)where, U is a third parameter that accounts for the unloading effect, and(28)where, σ′_{max} and V_{p(max)} are estimates of the effective stress and velocity at the onset of unloading. The units are in the imperial system and value of “5000” represents the water velocity equals to 5000 ft s^{−1}.
Therefore, the transformed Bowers original equation in the imperial units system, to estimate the pore pressure (P_{p}) incorporating the unloading effect can be rewritten using the SI units as:(29)
References
 Addis M.A., Last N.C., Yassir N.A. (1996) Estimation of horizontal stresses at depth in faulted regions and their relationship to pore pressure variations. SPE Formation Evaluation, 11, 1, 11–18, doi: 10.2118/28140PA. [CrossRef] [Google Scholar]
 Athy L.F. (1930) Density, porosity and compaction of sedimentary rocks, AAPG Bull. 14, 1–24. [Google Scholar]
 Bott M. (1959) The mechanics of oblique slip faulting, Geological Magazine 96, 2, 109–117, https://doi.org/10.1017/S0016756800059987. [CrossRef] [Google Scholar]
 Bowers G.L. (1995) Pore pressure estimation from velocity data: Accounting for overpressure mechanisms besides undercompaction, SPE Drill. Completion 10, 2, 89–95. [CrossRef] [Google Scholar]
 Christensen C.J., Powers G. (2013) Formation Evaluation Challenges in Tamar Field, Offshore Israel, SPWLA 54th Annual Logging Symposium, Society of Petrophysicists and Well Log Analysts, New Orleans, Louisiana. [Google Scholar]
 Eaton B.A. (1975) The Equation for Geopressure Prediction from Well Logs. Fall Meeting of the Society of Petroleum Engineers of AIME, SPE 5544, Soc. Petrol. Eng., https://doi.org/10.2118/5544MS [Google Scholar]
 Fjaer E., Holt R.M., Raaen A.M. (2008) Petroleum related rock mechanics, 2nd edn., Elsevier Science, Oxford, UK. [Google Scholar]
 Geertsma J. (1973) Land subsidence above compacting oil and gas reservoirs, Journal of Petroleum Technology 25, 6, 734–744, doi: 10.2118/3730PA. [CrossRef] [Google Scholar]
 Ghalayini R., Daniel J.M., Homberg C., Nader F.H. (2014) Impact of Cenozoic strikeslip tectonics on the evolution of the northern Levant Basin (offshore Lebanon), Tectonics 33, 11, 2121–2142. [CrossRef] [Google Scholar]
 Gheibi S., Holt R.M., Vilarrasa V. (2017) Effect of faults on stress path evolution during reservoir, pressurization, Int. J. Greenh. Gas Con. 63, 412–430. [CrossRef] [Google Scholar]
 Gheibi S., Vilarrasa V., Holt R.M. (2018) Numerical analysis of mixedmode rupture propagation of faults in reservoircaprock system in CO_{2} storage, Int. J Greenh. Gas Con. 71, 46–61, doi: 10.1016/j.ijggc.2018.01.004. [Google Scholar]
 Haimson B.C., Herrick C.G. (1986) Borehole breakouts – A new tool for estimating in situ stress?, in: Stephansson O., Stephansson O. (eds), Rock Stresses and Rock Stress Measurement, Proc. Int. Symp., Centek Pub., pp. 271–281. [Google Scholar]
 Healy D., Blenkinsop T.G., Timms N.E., Meredith P.G., Mitchell T.M., Cooke M.L. (2015) Polymodal faulting: Time for a new angle on shear failure, J. Struct. Geol. 80, 57–71. [CrossRef] [Google Scholar]
 Heidbach O., Rajabi M., Reiter K., Ziegler M. WSM Team. (2016) World Stress Map Database Release 2016,, GFZ Data Services, https://doi.org/10.5880/WSM.2016.001. [Google Scholar]
 Hottmann C.E., Johnson R.K. (1965) Estimation of Formation Pressures from LogDerived Shale Properties, J. Petrol. Tech. 17, 6, 717–722, https://doi.org/10.2118/1110PA. [CrossRef] [Google Scholar]
 Jaeger J.C, Cook, N.G.W. (1979) Fundamentals of rock mechanics, 3nd edn., Chapman and Hall, London. [Google Scholar]
 Jaeger J., Cook N.G., Zimmerman R. (2007) Fundamentals of rock mechanics, 4th edn., Blackwell Publishing. [Google Scholar]
 Jolly R.J.H., Sanderson D.J. (1997) A Mohr circle reconstruction for the opening of a preexisting fracture, J. Struct. Geol. 19, 887–892. [CrossRef] [Google Scholar]
 McKeagney C.J., Boulter C.A., Jolly R.J.H., Foster R.P. (2004) 3D Mohr circle analysis of vein opening, Indaramalodegold deposit, Zimbabwe: implications for exploration, J. Struct. Geol. 26, 1275–1291. [CrossRef] [Google Scholar]
 Montadert L., Lie O., Kassinis S. (2010) New seismic may put offshore Cyprus hydrocarbon prospects in the spotlight, First Break 28, 4, 91–101. [Google Scholar]
 Montadert L., Nicolaides S., Semb P.H., Lie Ø. (2014) Petroleum systems offshore Cyprus, in: Marlow L., Kendall C., Yose L. (eds), AAPG Memoir 106: Petroleum Systems of the Tethyan Region, 106, 301–334. [Google Scholar]
 Morris A.P., Ferrill D.A., McGinnis R.N. (2016) Using fault displacement and slip tendency to estimate stress states, J. Struct. Geol. 83, 60–72. [CrossRef] [Google Scholar]
 Morris A.P., McGinnis R.N., Ferrill D.A. (2014) Fault displacement gradients on normal faults and associated deformation, AAPG Bull. 98, 6, 1161–1184. [CrossRef] [Google Scholar]
 Nader F.H. (2011) The petroleum prospectivity of Lebanon: an overview, J. Petrol. Geol. 34, 2, 135–156. [CrossRef] [Google Scholar]
 Noble Energy. (2011, December 28) Noble energy announces significant natural gas discovery offshore Republic of Cyprus [Press Release]. Retrieved from http://investors.nblenergy.com/newsreleases/newsreleasedetails/nobleenergyannouncessignificantnaturalgasdiscovery . [Google Scholar]
 Papamichos E. (2010) Borehole failure analysis in a sandstone under anisotropic stresses, Int. J. Num. Anal. Meth. Geomech. 34, 6, 581–603. [CrossRef] [Google Scholar]
 Papanastasiou P., Vardoulakis I. (1992) Numerical treatment of progressive localization in relation to borehole stability, Int. J. Num. Anal. Meth. Geomech. 16, 389–424. [CrossRef] [Google Scholar]
 Papanastasiou P., Zervos A. (2004) Wellbore Stability: from linear elasticity to postbifurcation modeling, ASCE Int. J. Geomech. 4, 1, 2–12. [CrossRef] [Google Scholar]
 Papanastasiou P., Sarris E., Kyriacou A. (2017) Constraining the insitu stresses in a tectonically active offshore basin in Eastern Mediterranean, J. Petrol. Sci. Eng. 149, 208–217. [CrossRef] [Google Scholar]
 Plumb R., Papanastasiou P., Last N. (1998) Constraining the state of stress in tectonically active settings, SPE 47240 Proc. SPE/ISRM Eurock’98, Trondheim, Norway, pp. 179–188. [Google Scholar]
 Sibson R.H. (1985) A note on fault reactivation, J. Struct. Geol. 7, 751–754. [CrossRef] [Google Scholar]
 Skorstad A., Tveranger J. (2007) Fault facies – a new concept for data integration in reservoir fault zones, Conference: Data Fusion: Combining Geological, Geophysical and Engineering Data An SEG/AAPG/SPE Joint workshop, Vancouver. [Google Scholar]
 Streit J.E., Hillis R.R. (2002) Estimating fluid pressures that can induce reservoir failure during hydrocarbon depletion, Int. Rock Mechanics Conference, Texas, Irving, Paper SPE 78226. [Google Scholar]
 Yin Z.M., Ranalli G. (1992) Critical stress difference, faultorientation and slip direction in anisotropic rocks undernonAndersonian stress systems, J. Struct. Geol. 14, 237–244. [CrossRef] [Google Scholar]
 Zoback M.D. (2010) Reservoir geomechanics, Cambridge University Press, Cambridge, UK. [Google Scholar]
All Tables
Fault geometries including strike orientation and dipping angles determined in the modelling process.
Model assumptions, corresponding field stresses, and pore pressure at the approximate reference depth of the study sandstones zone.
Model assumptions for the three typical fault regimes plotted in a 3D Mohr diagram.
Normal faulting conditions: Fault stress analysis results for delta pressure to shear failure in MPa.
Strikeslip conditions: Fault stress analysis results for delta pressure to shear failure in MPa.
All Figures
Fig. 1 Interpreted faults on a structural analogue to Aphrodite and Tamar formations in the Levantine basin offshore Cyprus (Montadert et al., 2010). 

In the text 
Fig. 2 Model workflow in the fault stress analysis. 

In the text 
Fig. 3 Borehole diagram illustrating the areas of breakouts and drilling induced fractures. 

In the text 
Fig. 4 Main fault regimes as per Anderson’s classification (Heidbach et al., 2016). 

In the text 
Fig. 5 Range of stress magnitude for the main faulting regimes in an offshore field: (a) in normal faulting the critical S_{hmin} is highlighted, (b) in strikeslipping the critical S_{Hmax} is highlighted and (c) in reverse faulting the S_{Hmax} is highlighted. Modified after Zoback (2010). 

In the text 
Fig. 6 Schematic diagram illustrating: (a) the fault strike direction and dipping angle, and (b) the inclined fault plane in the 3D stress field. 

In the text 
Fig. 7 (a) Frictional faulting notation applied on twodimensional force equilibrium body and (b) typical Mohr circle diagram. 

In the text 
Fig. 8 Stress tectonic regime study in the Levant basin offshore Lebanon (modified after Ghalayini et al., 2014). 

In the text 
Fig. 9 World stress map generated using the CASMO online interface illustrating a series of stress indicators in the Eastern Mediterranean region (Heidbach et al., 2016). The inset shows the study area location on the regional seabed map after EMODnet bathymetry (EMODnet, 2018). 

In the text 
Fig. 10 Overview of the Tamar field showing (a) location map, (b) structure map of top reservoir, and (c) structural crosssection (NW to SE) through the Tamar1 discovery well T1 (Christensen and Powers, 2013). 

In the text 
Fig. 11 (a) Seismic section and (b) depth structural map illustrating the faults pattern in the Aphrodite field offshore Cyprus (Hydrocarbons Service – MECIT, Courtesy: Noble Energy International Ltd). 

In the text 
Fig. 12 Sedimentary stratigraphic section along the Aphrodite field in SWNE orientation (Hydrocarbons Service – MECIT; Courtesy: Noble Energy International Ltd). 

In the text 
Fig. 13 Determination of V_{p} maximum and σ_{max} from V_{p} log at well A (Hydrocarbons Service – MECIT; Courtesy: Noble Energy International Ltd) 

In the text 
Fig. 14 (a) Pore pressure prediction at depth and (b) the corresponding estimated effective stress profile. 

In the text 
Fig. 15 Illustration of the three typical fault regimes delineated over a top view plane. 

In the text 
Fig. 16 Illustration of the three common fault stress states in a 3D Mohr diagram for (a) normal faulting and (b) strikeslip stress regimes. The dashed lines in grey colour indicate the stress transformation with different fault plane orientation and dipping. 

In the text 
Fig. 17 Mohr circle plots showing the slip factor for the modelled (a) normal faulting and (b) strikeslip stress regimes. 

In the text 
Fig. 18 Contour diagrams of the slip factor as a function of α, θ angles for: (a) normal faulting and (b) strikeslip stress regimes. 

In the text 
Fig. 19 Fault stress plot in normal faulting regime. Model parameters definition for different strike directions relative to S_{Hmax} and dipping angles. 

In the text 
Fig. 20 Fault stress plot in strikeslipping regime. Model parameters definition for different strike directions relative to S_{Hmax} and dipping angles. 

In the text 
Fig. 21 Schematic diagram showing the estimation of delta pressure to shear failure. 

In the text 
Fig. 22 Influence of the pore pressure changes on the faultplane stability conditions. 

In the text 
Fig. 23 Coordinate system that utilizes the zenith and longitudinal angles (Jaeger et al., 2007). 

In the text 