Petroleum geomechanics modelling in the Eastern Mediterranean basin: analysis and application of fault stress mechanics Dynamics of sedimentary basins and underlying lithosphere at plate boundaries: The Eastern Mediterranean

. A fault stress analysis of a typical gas ﬁeld 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 Mohr-Coulomb 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 ﬁeld 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.


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 disk-shaped 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 Mohr-Coulomb 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 Mid-Jurassic passive rifting tectonic episode, which followed by further regional extension and spreading, with periodic eustatic sea-level 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 strike-slip 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 Pliocene-Pleistocene 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). 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 high-quality 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.

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 geo-pressure developed at depth. The model methodology is based on the Mohr-Coulomb 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.
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.
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 three-dimensional 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, fault-cutoff relationships, and small-displacement 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 geo-body, 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 three-dimensions (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

Fault and stress regimes
The main fault regimes, as per Anderson's classification system, are defined as normal faulting, strike-slip, 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. 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) strike-slip 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 strike-slip 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 strike-slip.

Fault plane stress determination
Mohr-Coulomb frictional faulting theory is applied to estimate the normal and shear stresses for each fault plane orientation in the 3D stress space. The two-dimensional (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). Three-dimensional Mohr diagram is also used to explain mechanism of faulting and reactivation of pre-existing 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 h and its strike direction is oriented at an angle a in relation to the maximum horizontal stress, S Hmax .
The simple form of the governing equations are based on the frictional faulting theory for normal (r 0 n ) and shear (s) 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 y-axis is parallel to the major principal stress and the x-axis 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 r 0 n and s are acting on an inclined plane with normal vector forming a direction b relative to the minor principal stress r 3 , which is parallel to the x-axis (Fjaer et al., 2008).
The values of r 0 n and s 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 (r 1 À r 3 )/2 and its centre is at the point (r 1 + r 3 )/2 on the r-axis. The coefficient of friction, l equals to tanu, where u is the friction angle. The failure plane forms an angle h = 45°+ u/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 ENE-WSW dextral strike-slip faults, onshore ENE-WSW latitudinal dextral strike-slip faults, and NE-SW sinistral (left-handed) strike-slip 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 ENE-WSW dextral strike-slip 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 NW-SE normal faults in the deep basin might not be caused by these regional geodynamics but rather to a local stress-field fluctuation affecting only the Oligocene-Miocene units (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).
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 strike-slip faulting. Figure 10 illustrates the location, depth map and a cross-section of the field (Christensen and Powers, 2013).

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.
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 h, and the triangular h-opposite angle, b =90°À h. Also, the relative to S Hmax orientation angle a, is calculated as, a = strike azimuth -S Hmax azimuth.
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 Miocene-Oligocene sandstones zone. The Bowers method is an effective stress approach that employs the virgin and unloading curve relations to account for both under-compaction 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  . 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). 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 r 0 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: where, U is a third parameter introducing the unloading effect, and The r 0 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: 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.
Fracture pressure is approximated using two equations based on the frictional faulting theory and coefficient of friction, l = 0.6 for (a) normal faulting regime and (b) the frictional equation for strike-slip regime which is derived as an intermediate stress state between the normal  faulting and vertical stress. These equations are given, respectively by 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 r 3 , the value of r 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 r 1 and r 3 that corresponds to critical stressed limit at l = 0.6 can satisfy the ratio r 1 /r 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 strike-slip 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, r 0 , is calculated from the total principal stress minus the pore pressure according to the Terzaghi effective stress theory.

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 strike-slip faulting (SS) regimes were tested. The aim of these multiple calculations is first to understand better the trend-line 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 strike-slipping. This information is used as an assumption in this analysis.

Parametric studies
The parametric calculations evaluated first the stress state of three common fault plane geometries: normal fault, strike-slip 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 r 3 and r 1 . However,      it is worth-mentioning that during the tectonic evolution of fault mechanics, the intermediate stress state between the normal faulting and reverse faulting is the strike-slip. 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, strike-slip, 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 a = 0°and h = 60°, for strike-slip is a = 30°and h = 90°and for reserve faulting is a = 90°and h = 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.
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 post-failure 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 (r 1 -r 3 ), the strike-slip is located in the circle (r 2 -r 3 ), and the reverse faulting in the circle (r 1 -r 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 (r 2 -r 3 ), the strike-slip moves to the upper part of large circle (r 1 -r 3 ), and the reverse faulting remains on the circle (r 1 -r 2 ). Stress conditions that approach or exceed the shear failure line as shown in the plot (for l = 0.6), can be considered as potential active events, or currently on a slippage mode.
Processing the calculated results we constructed the 2D Mohr stress nomograph as a function of the horizontal deviation a from S Hmax and of the vertical inclination b (the complimentary angle to dipping h). Figure 17 shows the calculated slip factor for the (a) normal faulting and (b) strike-slip regimes. The slip factor is calculated as the ratio of shear stress to normal effective stress (s/r 0 ).The upper limit of this slip factor is l = tanu = 0.6. Figure 18 shows the contours of the slip-factor on the (a,h) angle space.

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 strike-slip regime.
The results are plotted on 3D Mohr diagrams for a total of 15 realizations in each regime: three strike fault angle scenarios of a = 0°, 30°and 90°in relation to maximum principal horizontal stress, S Hmax and five fault dipping cases for each scenario, with dip h = 90°, 80°, 60°, 45°and 0°( b = h À 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 strike-slip regime, the stress failure risk is increased as the fault strike direction approaches the angle a = 30°relative to S Hmax . In terms of dipping angle, in normal faulting regime the failure line is approached at dip, h = 60° (Fig. 19a), while in strike-slip regime the failure line approaches the dip, h = 90° (Fig. 20b).

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 fluid-flow with geomechanics more important.
Fault sensitivity analysis to pore-pressure changes is performed in a series of calculations and the results are plotted on the Mohr-Coulomb diagram for shear failure, defined by the frictional coefficient l = 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.
The results are summarized in Table 4 for normal faulting regime and in Table 5 for strike-slip 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 sub-critical, 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 sub-critical. The rest values suggest stable conditions.
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 strike-slipping 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 trend-line 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 non-uniform compaction, the normal effective stress increases and shifting the Mohr-circle towards the shear slip line.
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 leak-off testing), and borehole image logs analysis for determining the orientation of the horizontal stresses.

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 trend-line 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. s y 0 z 0 ¼ À 1 2 r 1 À r 2 ð Þsin h sin 2k ð22Þ s x 0 z 0 ¼ 1 2 r 1 cos 2 k þ r 2 sin 2 k À r 3 Â Ã sin 2h: 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 drill-stem 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: where: P p : hydrostatic pressure in Pa; q: 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: where q (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 (1995), 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 , where: V i : interval velocity in ft s À1 ; r 0 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: where, U is a third parameter that accounts for the unloading effect, and where, r 0 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: