Determination of stress state based on well logging data and laboratory measurements – a CBM well in the southeastern part of the Upper Silesian Coal Basin (Poland)

The main objective of this study is to present calculation methods of horizontal stress profiles, taking into account the stress boundaries model, poro-elastic horizontal strain model and the effective stress ratio approach, using calibration with wellbore failure. The mechanical earth model (MEM) parameters from log measurements and well testing data were estimated for a well located in the southeastern part of the Upper Silesian Coal Basin. Log-derived horizontal stresses of the well are commonly treated as the final product of geomechanical modeling in oil and gas practices. A less popular method for estimating horizontal stresses is based on Kirsch equations juxtaposed with compressional and tensile failure observed on a micro-imager or six-arm caliper. Using this approach, horizontal stresses are determined based on the fact that when hoop stresses exceed the formation’s tensile strength, tensile fractures are created, and when those stresses exceed the compressive strength of the formation, breakouts can be identified. The advantage of this method is that it can be run without in situ stress measurements. The presented workflow is recommended every time there is an image log and dipole sonic measurement in the available dataset, both being necessary to observe the failure zones and MEM.


INTRODUCTION
Geomechanical modeling has become a fundamental element of the hydrocarbon exploration process, supporting the reduction of risks associated with wellbore instability (Słota-Valim 2014) and the designing of hydraulic fracturing treatment. Understanding the mechanical properties of the reservoir will enable its subsequent optimal development. In the case of structurally complex reservoirs, such knowledge is indispensable for building static and dynamic models (Dudek et al. 2020). Thus, improving methods of predicting stress state and failure in hydrocarbon wells is of unceasing concern for the industry.
In this study, the mechanical earth model (Zoback 2007, Jędrzejowska-Tyczkowska & Słota-Valim 2012) was built for the vertical well X (Fig. 1) in six steps (Fig. 2) by using log measurements and well-testing data from well X and information obtained from offset wells X-1, X-2H, X-3K and X-4H located at a distance of 1.5 km within the southeastern part of the Upper Silesian Coal Basin. The MEM was created in Schlumberger Techlog 2019.1. Several findings based on the drilling reports provide constraints for the MEM. Neither mud losses nor kicks were reported in well X. No hole instability (shear failure, i.e. cavings) was reported either. Apart from the measurements from hydraulic fracturing treatments in X-2H and X-4H, no strong calibration for stress magnitude and direction, or rock properties (rock strength, elastic properties) is available in the dataset from the offset wells. The horizontal stress profiles of the well under study were determined with the use of the presented workflow, and their outputs were checked against the compressive and tensile wellbore failures interpreted in the study.

GEOLOGICAL SETTING
The Upper Silesian block is a part of the Brunovistulia Terrane (Kotas 1982, Buła 2000, Finger et al. 2000, Żelaźniewicz et al. 2009), the northern part of which extends from the area around Vienna-Brno-Wrocław in the west, to the Kraków-Lubliniec Fault Zone in the northeast (Buła & Żaba 2005. The southern part of this unit lies beneath the Carpathians and reaches at least the Peri-Pieniny Fault Zone (Kalvoda et al. 2003). The studied area is located in the Upper Silesian Coal Basin (USCB) which was formed during the last phases of development of the Moravo-Silesian Paleozoic Basin (Kotas 1995, Zdanowski & Żakowa 1995, Różkowski 2004, Aleksandrowski et al. 2011. USCB stretches within the foreland of the Variscan belt, in the NE of the Brunovistulicum basement dating back to Cadomian age (Kotas 1982, Buła 2000. In Devonian times, siliciclastic sedimentation took place, with subsequent sedimentation of carbonates which continued in the Tournaisian (Kotas 1982). In the Early Carboniferous, the deposition of siliciclastic sediments of deep-sea flysch (Culm) began, gradually transforming into the sedimentation of a paralytic and terrestrial carbonate molasse (Kotas 1982). The Carboniferous molasse formations cover the geological time range from the Namurian A to the Westphalian D (Zdanowski & Żakowa 1995). The lithological profile of USCB, whose thickness reaches 8000 m in the center, is composed of two main parts comprised of four lithostratigraphic units. The lower part is the Paralic Series, characterized by an extensive range of clastic deposits, including marine, terrestrial-marine, as well as thin coal seams. The portion of sandstones in the profile of this series is 20-50% and of coals 3-4% (Zdanowski & Żakowa 1995). The upper part of the coal-bearing sediments consists of continental sediments divided into three units: the Upper Silesian Sandstone Series, Mudstone Series and Cracow Sandstone Series (Dębowski 1973, Kotas 1982, Jureczka & Kotas 1995, Słoczyński & Drozd 2018. The thickness of all series varies and decreases toward the east (Kotas & Porzycki 1984). The Upper Silesian Sandstone Series is dominated by sandstones and conglomerates, accompanied by inserts of clay and silt rocks, as well as a few coal seams. In turn, the Mudstone Series is dominated by mud and clay sediments, sometimes accompanied by sandstone layers and numerous coal seams. The profile of the Upper Carboniferous ends with the Cracow Sandstone Series, characterized by the predominance of sandstones and conglomerates (>70%) and by the presence of a few coal seams (Zdanowski & Żakowa 1995). The boundary with younger rocks is unconformable and erosional paleosurface is irregular. The Carboniferous is covered by Permian rocks in the northeastern part of the USCB, by Triassic clastic and calcareous rocks in the northern part and by Jurassic rocks in the eastern part (Słoczyński & Drozd 2018). The lowest part is the Miocene sequence covering a large part of the basin, and consisting of clays, claystone and conglomerate (Słoczyński & Drozd 2018).
The USCB comprises a variable geological structure in terms of both stratigraphy and tectonics. Within the basin, fold tectonics, block tectonics, and fault-and-block tectonics, all related to Variscan orogenesis can be recognized. Recent tectonic stresses development is associated with Alpine orogenesis (Kotas 1985, Dubiński et al. 2019. The stress state of the studied area has been influenced by crustal tectonic processes together with gravitational stress caused by overburden, and stress perturbation caused by mining activity (Dubiński et al. 2019).

GEOMECHANICAL MODELING Overburden stress (S v )
The magnitude of S v can be calculated by the integration of all rock densities from the surface down to the depth of interest z (Zoback et al. 2003): Due to issues with the caliper and lack of density measurements in the uppermost interval to derive the overburden (S v ), pseudo-density data from the Gardner equation (where a = 0.23 and b = 0.25) was used (Gardner et al. 1974), calibrated with bulk density logs (Fig. 3): (2) where: To derive a depth-continuous vertical stress (overburden) profile (S v ), bulk density logs from the investigated well were integrated with the derived exponential trend line and pseudo density. The exponential trend line, generalized for all lithological types, was determined by adjusting the reference points to match the density log over the depth interval for which the density data were available. Density is quite low at the shallow intervals and increases with depth. Hence, the manner in which measured values are extrapolated to the surface is important (Zoback et al. 2003). The resulting overburden density (mud-weight equivalent EMW) at the depth of 1000 mMD is 2.3 g/cm 3 .

Pore pressure (P p )
According to the interpretation of the test data in the well X-2H made by the Reservoir Engineering Division (PGNiG), there is no overpressure in the wells X-2H. In the course of further analysis, hydrostatic pore pressure was assumed due to the lack of evidence for any overpressure (no inflows during the drilling).

Rock properties
Where shear velocity logs have not been acquired, the empirical Castagna et al. (1985) relationship was used. Shear velocity modeling with the lithology-dependent original Greenberg-Castagna VP-VS relations (VP -compressional velocity [km/s], VS -shear velocity [km/s]) is used for both quality control of the shear slowness measurements and the shear slowness prediction. Due to the lack of data covering the whole well, the simplified version of the equation was used. The parameters thus obtained were calibrated with dipol sonic measurements ( Fig. 4): where: Δt shear_Castagna -shear slowness calculated using the Castagna relationship, Δt comp -compressional slowness, a, b -the coefficients used, −0.002223 and 0.748408, respectively, based on calibration with measured shear slowness.
The mechanical properties of the rocks determine the susceptibility to failure and deformation, thus affecting the stress field (Jędrzejowska-Tyczkowska & Słota-Valim 2012). The dynamic elastic properties were determined directly from the sonic and density data (Gassmann 1951).

Fig. 4. Calibration of shear slowness from the Castagnia relationship in well X: DT -compressional slowness, DTS -measured shear slowness, DTS (Castagna) -calculated shear slowness
The mechanical properties derived from the above equations are dynamic properties. For geomechanical applications, they need to be converted into static properties. The conversion from the dynamic to the static Young Modulus was carried out using the correlation with static mechanical properties obtained from laboratory core measurements performed at the University of Warsaw -Geomechanics Department, Faculty of Geology (Łukaszewski et al. 2019) (Fig. 5): Based on the empirical relationship between wireline acoustic logs and rock strength, the Unconfined Compressive Strength (UCS) profiles for each lithological unit (coal, sandstone and shale) were derived. For each lithotype, triples of plugs from approximately the same depth were available. They were treated with different confining pressure, making it possible to determine the value of the unconfined compressive strength (Tab. 1). Laboratory core measurements were used for calibration. The UCS was based on the McNally (1987), Horsrud (2001) and Kumar et al. (2015) correlations. Within the coal seam, the UCS was iteratively fitted during wellbore stability modeling. The final UCS results from multiplying the combined UCS by a 1.5 gain to calibrate to the laboratory measurements: coal (Kumar et al. 2015, modified): where UCS is in MPa and V p in m/s, shale (Horsrud 2001): where Δt comp is in μs/ft, As the tensile strength of rock is commonly in the order of 1/12 to 1/8 of its UCS, the tensile strength was calculated from UCS with 0.1 gain. Internal Friction (FANG) was calculated based on the empirical correlation with V p log (Lal 1999). Laboratory core measurements performed at the University of Warsaw were used for calibration: where V p is in m/s.

Maximum horizontal stress azimuth
The interaction of wellbore with regard to lithology and the recent stress field leads to the creation of shear and tensile failures. In order to determine with confidence, the maximum horizontal stress (S Hmax ) in the wells under study, the interpretation of the breakouts (BO) is necessary. Such structures form perpendicularly to S Hmax (Bell & Gough 1979). (Heidbach et al. 2008, modified) with location of X-2H and orientation of S Hmax marked accordingly (B) A B To detect and analyze breakouts, and thus to determine the direction of S Hmax , typically micro--imagers and six-arm calipers are used. Since, according to image log within target zone, breakouts and DITFs do not exist in this particular case, the direction S Hmax = 150° ±10° was assumed, based on the results of microseismic monitoring (distribution of events) during hydraulic fracturing in the X-2H well (Fig. 6A) and the World Stress Map (Heidbach et al. 2016) (Fig. 6B). This is in agreement with previously published data from the region (Jarosiński 2005(Jarosiński , 2006.

Horizontal stresses and wellbore stability (WBS)
For well X, the Mohr-Coulomb stress boundaries model, the poro-elastic horizontal strain model and the effective stress ratio approach were used.
The Mohr-Coulomb stress boundaries model (Zoback 2007) is a failure model which describes the relationship between two stresses (on the assumption that there is no impact of intermediate principal stress upon failure), presuming that the formation is at failure. Supposing that the vertical stress is a principal stress, the limits of horizontal stresses in the stress domain are the lower limit of minimum horizontal stress and the upper limit of maximum horizontal stress. Both are obtained from the Mohr-Coulomb stress boundaries model. This model does not give the correct stress profiles for a particular case, but it provides minimum and maximum stress limits which cannot be crossed without any failure of rock.
On the basis of the poro-elastic horizontal strain model, the most commonly used method for horizontal stress calculation, and taking into account horizontal tectonic stress, the magnitudes of S Hmax , S Hmin were calculated. This equation is also known as the evolution of Eaton's equation. Assuming a flat-layered poro-elastic deformation in the formation rock, particular constant strains S Hmin and S Hmax are applied to the formation in the directions of minimum and maximum stress, respectively. The poro-elastic horizontal strain model can be expressed using the static Young's modulus, Poisson's ratio, Biot's constant, overburden stress and pore pressure. S Hmin and S Hmax cannot be directly measured by adjusting these strains, yet the calculated stresses can be calibrated with the measured horizontal stresses at depth. In this analysis, 0.00025 and 0.0004 were used respectively, developed by matching modeling results to observation of the borehole. Reasoning in the category of a sedimentary basin with a span of, e.g. 100 km, this will give an absolute deformation of 25 m and 40 m, respectively.
For a fluid-saturated porous material which is assumed to be linear, elastic and isotropic, considering anisotropic tectonic strain, the horizontal stresses (S hmin and S Hmax ) are equal to (Blanton & Olson 1999): where: S hmin -minimum horizontal stresses, S Hmax -maximal horizontal stresses, S v -vertical stress, a -Biot's coefficient, P p -pore pressure, ν -Poisson's ratio, E -Young's modulus, ε x , ε y -strains in the minimum horizontal stress and maximum horizontal stress direction, respectively.
The equations take into account the phenomenon of a stronger rock (with a higher Young's modulus) supporting higher horizontal stress in a tectonically active area. This model can account for situations where sandstones are under higher horizontal stress than the neighboring shales. When the strains are close to zero, this model could be considered a uniaxial strain model.
The model does not deliver reliable and realistic results due to incorrect assumptions. From a geological perspective over time, the earth is not elastic (well-cemented reservoirs are elastic-plastic -with stress at a frictional limit; poorly cemented reservoirs are viscoplastic -with stresses relaxing over time) and stress is not instantaneous. The existence of consistent directions of principal stresses over broad regions is an obvious manifestation of anisotropic magnitudes of the horizontal stress. Moreover, tectonic sources of stress often result in one (or both) of the horizontal stresses exceeding the vertical stress, as required in areas of strike-slip or reverse faulting. Attempts to make corrections for existence of tectonic stress using arbitrary coefficients only make the situation worse by adding more empirically determined parameters (Zoback 2007). Finally, the effective stress ratio approach was implemented. At depths where reliable failure observation and stress calibration were available, the effective stress was calculated as follows: where: ESR Shmin -effective stress ratio for S hmin , S v -overburden stress, S hmin -ISIP stress measured from frac in G-2H, P p -pore pressure.
The resulting points were then fitted with a curve from which a continuous S hmin trend was calculated. The calibration was based on the interpreted frac gradient from the Mini-Frac in X-2H in the interval 2046.5-1973.1 m -0.236 bar/m and from the DataFRAC in X-4H in the interval 2143.4-2116.7 m -0.243 bar/m (horizontal sections).
An ESR Shmin of 0.98 was used. In this study, ISIP is taken as a reasonable estimate of the lowest stress. However, the best estimate comes from the fracture closure pressure (FCP) observed after ISIP. The FCP should be lower than the ISIP, and in the well X-2H this was not the case, which is why the ISIP was used. It is possible that the step-down affected the decline associated with the closing of the fracture. The ISIP can be accepted as an upper bound of S hmin . During the DataFRAC, the FCP was not observed either.
The S Hmax magnitude cannot be measured; observations of wellbore failure along with the other geomechanical model parameters are used in order to determine the S Hmax magnitude. In other words, using Techlog software, back-calculation is performed to establish the S Hmax magnitude necessary to create the observed failure, with the other parameters for the given hole orientation. To estimate the stress redistribution around the wellbore drilled in an elastic medium with applied stresses, the Kirsch equations (Kirsch 1898) can be used: where σ ΔT -the thermal stresses caused by the difference between the mud temperature and formation temperature.
BO form when the hoop stress (σ θθ ) is the largest and exceeds the UCS, whereas DITF form when the hoop stress (σ θθ ) is below zero and exceeds the T 0 of the near-wellbore rock (Aadnoy 1990).
To satisfy the regional experience (Carpathians) (Jarosiński 2005), a ESR SHmax of 1.2 was used as datum point. Several other values were also tested.
In order to verify the geomechanical model parameters, the failure was predicted (Fig. 7A) taking into account the well trajectory and the mud weight (the lowest mud weight which was experienced by the particular hole interval during operations). These predicted breakouts were compared to observations of wellbore failure from the imager. The predicted failure generally matched the failure observed on the imager. There may have been slightly too much failure predicted as a result of uncertainties regarding the model parameters. Due to the fact that in the well under study very few breakouts were observed (at the depth of 1060.0-1060.5 m) and that the drilling-induced tensile fractures (DITFs) are not clear, there is huge uncertainty as to the model. The features interpreted by Geofizyka Toruń SA (Fig. 8)

(Final Evaluation
Report of well X) as DITFs are in most cases probably drilling-enhanced natural fractures. Where the borehole hoop stress is tensile, either from in-situ stresses or drilling/coring practices, the intersection of a natural fracture or foliation plane with the tensile region of the borehole may be preferentially opened in tension. Drilling-enhanced natural fractures can be easily mistaken for inclined or J-hook tensile wellbore failures, resulting in serious errors in geomechanical modeling.
The stress polygon (Fig. 9) demonstrates the range of possible values for horizontal stresses (assuming a particular depth, P p and the coefficient of friction), for normal, strike-slip and reverse faulting regimes, respectively, according to Anderson's stress and faulting classification system, and the relationship described by Jaeger & Cook (1979). They showed that the values of σ 1 and σ 3 correspond to the situation where a critically oriented fault is at the frictional limit (Equation (23) where: P p -pore pressure, μ -coefficient of frictional sliding. To predict the limiting stress differences at depth, Anderson's faulting theory (Anderson 1951) has to be applied to determine which of the stresses S Hmax , S hmin , and S v correspond to principal stresses S 1 , S 2 , and S 3 . This will depend upon whether it is a normal, strike-slip, or reverse faulting environment: normal faulting: The crosscut of the horizontal stresses, independent of the stress state, must be placed within the stress polygon. The pair of horizontal stresses which falls on the outer margin of the polygon results in fault reactivation. The stress polygon often allows a wide range of stress values at depth of each analysis. In order to limit this range, leak-off tests or hydraulic fracturing for S hmin estimation and the observation of BO and DITF for S Hmax are necessary (Zoback et al. 2003, Zoback 2007).
An analysis of stress polygons (using the values from the previously built MEM), assuming the existence of DITFs, indicates there is a huge anisotropy of horizontal stress, i.e. an ESR SHmax of approximately 2.9 at ESR Shmin of 0.98. Taking into account the calculated elastic properties and rock strength, such a stress state would cause the well to be very hard to drill (and to be almost without a safe mud window) (Fig. 7B) which is not the case either in the described well X or in offset wells, resulting in no reported instability during drilling. Moreover, the S Hmax thus obtained exceeds the upper limit of the S Hmax calculated from the Mohr-Coulomb stress boundaries model (SHMAX_MC_UB) (Fig. 10).

Fig. 10. S Hmax and S hmin calculated from the Mohr-Coulomb stress boundaries model (MC), poro-elastic horizontal strain model (PHS) and effective stress ratio
To reproduce the wellbore condition in well X, the modified Lade failure criterion (Lade 1977) with its unconfined compressive strength (UCS) and the angle of internal friction (φ) to assess failure was used. It is well known that the Mohr-Coulomb failure model is a conservative model, meaning that it shows results which are worse for drillers (more failure) than actual, as it ignores the effect of intermediate principal stress (Colmenares & Zoback 2002). Lade (1977) developed a three-dimensional failure criterion for frictional materials without effective cohesion. Ewy (1999) modified the original Lade criterion to calculate the linear shear strength increase with increasing mean stress, I 1 /3. While considering materials with non-zero cohesion, Ewy (1999) incorporated P p as a required parameter and entered additional material constants S related to the cohesion of the rock and η (internal friction) (Zoback 2007).

CONCLUSIONS AND DISCUSSION
In this study, the stress boundaries model, poro-elastic horizontal strain model and effective stress ratio approach combined with an analysis of wellbore failure were used to derive continuous stress profiles. In this approach, the stresses obtained through the application of different methods could be compared, leading to a reliable assessment of the stress state which is a crucial issue for future projects in this region. The post-drill geomechanical model for well X was calibrated with well data from the wells G-1, G-2H, G-3K and G-4H. Although calibration is still not altogether satisfactory, the model verification process shows that the MEM parameters work well within the uncertainty limits.
All in all, the conclusions from the application of the geomechanical model in this study can be described as follows: -The strike slip stress regime: S hmin < S v < S Hmax has been modelled, being consistent with the published data from this region (Jarosiński 2006). However, the difference between values of S v and S hmin recorded during fracturing treatment in both G-2H and G-4H is very small, which has never been indicated by previous authors regarding this part of Poland (Jarosiński 2006). Since the relation is confirmed in both offset wells, this should be treated with high confidence. Such a state requires a careful overburden calculation. The source of tectonic compression most likely is Alcapa's push transmitted by Outer Carpathians (Jarosiński 2006). Nevertheless, it has to be confirmed by further treatments expected during following CBM's exploration and production. Owing to the presence of a strike-slip regime, propagation of the vertical fracture during hydraulic fracturing is expected in direction of S Hmax , which should be around 150° ±10 o , according to microseismic events recorded during fracturing in X-2H, observations from image data, and published data (Jarosiński 2005). This direction is associated with the fore-Carpathian stress domain (Jarosiński 2005). It is worth noticing that the shadowing effect can potentially cause the flipping of the orientation of hydrofrac to a horizontal position because of stress regime change. This effect is related to the process of stress state disturbance and a subsequent gradual increase of minimum horizontal stress during following stages, due to the injection of frac fluids into the formation. Another option is growth and propagation upwards, where the stress contrast is lower (Dohmen et al. 2015). Such a process was observed in the G-4 well. The impact of mining activity on the recent stress state should be examined in the course of further analysis due to the presence of Brzeszcze coal mine in the vicinity. Local stresses could be affected as well because of Jawiszowice fault located to the south of the investigated area. -Hydrostatic pore pressure was assumed due to the lack of evidence for any overpressure (no inflows during drilling), i.e. the pore pressure gradient does not exceed 0.013 MPa/m. This assumption seems to be valid since no effective overpressuring mechanism was identified on the basis of the available data. Regarding microscale analysis, small variation of pore pressure should be considered. However, such an approach will be challenging due to limited measurements. Further exploitation of CBM and accompanying water production, resulting in pressure decrease, will impede the determination of the initial pore pressure; Determination of stress state based on well logging data and laboratory measurements...
-The overburden gradient at the depth of 1000 mMD is estimated to be 2.3 g/cm 3 . An overburden gradient up to 328 mMD was determined using pseudo-density from acoustic log, due to the lack of density logs from the uppermost interval. -The magnitudes of S hmin and S Hmax are modeled using both the evolution of Eaton's equation and the effective stress ratio method. The S hmin magnitude was constrained using Mini-Frac projected from well X-2H. ESR SHmax was applied to satisfy observed well performance (failure of borehole). The analysis of stress polygons combined with the observation of wellbore failure evidence ruled out huge stress anisotropy, initially suggested by the interpretation of DITFs. Despite the range of horizontal stress anisotropy at the regional scale, the discussion should be continued with the acquisition of new data. There are uncertainties arising from the available dataset: -Overburden: the density log was not recorded along the whole wellbore. -Pore pressure: uncertainty exists as there are no calibration points in the overburden. -Minimum horizontal stress magnitude (S hmin ): uncertainty exists as to the calibration of S hmin . Extended leak-off tests with pressure-time curves in shales in future wells will reduce the uncertainty margin in S hmin and also help in the modeling of the S Hmax magnitude. -Maximum horizontal stress magnitude (S Hmax ): it is important to understand that uncertainties related to issues such as differentiating borehole failure due to mechanical instability and time-dependent instability can have a huge impact on S Hmax calibration. -Rock strength: in the next wells to be drilled, it may prove useful to perform more rock strength (uniaxial, triaxial or multistage) tests with a spread across profiles.
The author wishes to thank PGNiG SA for its permission to publish this paper. He is also indebted to the reviewers for their constructive and helpful comments and corrections. Considerable thanks are extended to K. Kępińska for proofreading this translation.