An integrated approach to discretized 3D modeling of geomechanical properties for unconventional mature field appraisal in the western Canadian sedimentary basin ORIGINAL PAPER - EXPLORATION GEOLOGY An integrated approach to discretized 3D modeling of geomechanical properties for unconventional mature field appraisal in the western Canadian sedimentary basin Kalyan Saikia1 • Clara E. Ikuku2 • Bhabesh C. Sarkar3 Received: 14 September 2016 / Accepted: 20 November 2017 / Published online: 30 November 2017 � The Author(s) 2017. This article is an open access publication Abstract In mature field appraisal and development, dis- cretized geomechanical property models play a vital role in providing information on in situ stress regime as a guide for placement of directional wells. Laboratory methods of measuring these properties, in most cases, take only small samples from consolidated rocks. These isolated samples may not be representative of the entire elastic regime existing in the reservoir owing to sample size. In general, geomechanical studies are performed on a well-by-well basis and then these measurements are used as calibration points to convert 3D seismic data (if available) to geome- chanical models. However, elastic properties measured this way are restricted to the well location and interpolation across the reservoir may not be always appropriate. To overcome these challenges, this paper describes an inte- grated approach for deriving 3D geomechanical models of the reservoir by combining results of 3D geocellular models and basin models. The basin model reconstructs the geologic history (i.e., burial history) of the reservoir by back-stripping it to the original depositional thickness. Through this reconstruction, the mechanical compaction, pore pressures, effective stress, and porosity versus depth relationships are established. Next, these mechanical properties are discretized into 3D geocellular grid using empirical formulas via lithofacies model even if no 3D seismic data are available for the reservoir. The dis- cretization of elastic properties into 3D grids results in a better understanding of the prevailing stress regimes and helping in design of hydraulic fracturing operations with minimal risks and costs. This approach provides an inno- vative way of determining effective horizontal stress for the entire reservoir through 3D distribution of elastic properties. Keywords Reservoir � Geomechanics � Basin model � Unconventional � Geostatistics Introduction In unconventional tight reservoirs, where hydraulic frac- turing is the key technique for enhanced oil production, a thorough knowledge of distribution of geomechanical properties helps to maximize return on investment. Hori- zontal wells with multistage hydraulic stimulation are the primary production strategies during development of the unconventional tight reservoirs. It is therefore important to understand in situ stress conditions such as effective min- imum horizontal stress as it guides the stimulation design which will in turn control the fracture propagation. Accu- racy of effective minimum horizontal stress model is pri- marily dependent on knowledge of the spatial distribution of reservoir’s elastic property. In general, for upstream reservoir modeling processes, reservoir geomechanical properties are determined either by (1) analysis of seismic data (S-wave and P-wave velocities) or (2) by laboratory measurement of core plugs or (3) a combination of both techniques. These are con- ventional laboratory methods which may not be repre- sentative of the entire elastic regime existing in the & Kalyan Saikia kalyansaikia@gmail.com 1 Halliburton, 24335 Piazza Drive, Richmond, TX 77406, USA 2 Clenik Petrotech, 104 Tremblant Way SW, Calgary, AB T3H0A7, Canada 3 Department of Applied Geology, Indian Institute of Technology (Indian School of Mines), Dhanbad 826004, India 123 J Petrol Explor Prod Technol (2018) 8:417–429 https://doi.org/10.1007/s13202-017-0406-3 http://orcid.org/0000-0002-8500-9690 http://crossmark.crossref.org/dialog/?doi=10.1007/s13202-017-0406-3&domain=pdf http://crossmark.crossref.org/dialog/?doi=10.1007/s13202-017-0406-3&domain=pdf https://doi.org/10.1007/s13202-017-0406-3 reservoir as in most cases these small samples are taken from consolidated rocks and used to measure geome- chanical properties. Furthermore, these measurements are performed on a well-by-well basis on core samples. Measurements taken at the core samples are often used in geomechanical analysis as calibration points to convert the 3D seismic data (if available) to 3D volumes of geome- chanical properties. Previous studies on the gross effect of mechanical properties on productivity of the unconventional reservoir clearly elaborated that the use of only seismic velocities and elastic properties measured on samples are not suffi- cient enough to describe spatial distribution pattern of elastic properties for entire reservoir. Havens (2012) car- ried out similar approach for Bakken petroleum system to establish the spatial behavior elastic properties and con- cluded that conventional approach of consideration of small core samples plugs in formulating 3D reservoir geomechanical model is not a practical solution. During fracture modeling of unconventional tight gas reservoirs, Gonzalez et al. (2014) and Deng et al. (2011) demonstrated that incorporation of 3D geomechanical model gives a more realistic representation of the orientation and geom- etry of hydraulic fractures compared to traditional way of semi-analytical calculations of fracture geomechanical properties. Other published studies on similar topics in the western Canadian sedimentary basin reveal that spatial distribution of rock mechanical properties is highly heterogeneous in nature and is result of complex basin- forming (depositional) process (Ferdous et al. 2015; Lavoie and Séjourné 2016; Tong et al. 2005). To overcome these challenges, this paper presents a practical approach for deriving and discretizing geome- chanical and other elastic properties in unconventional tight reservoir by integrating results of 3D geocellular and basin models. This technique provides a relationship between petrophysical properties and elastic rock proper- ties of the reservoir rock, explaining their response to mechanical stresses at reservoir condition. Such relation- ship is very important for defining trends in the reservoir description. The current approach also incorporates wire- line logs representing the entire sediment sequence which can be used as input for effective minimum horizontal stress modeling for the entire reservoir. The case study demonstrated here assumes that elastic properties measured conventionally are restricted to the well locations and cannot be interpolated across the reservoir if no 3D seismic data are available. Integrated modeling approach The integrated approach of deriving discretized geome- chanical model of the reservoir can be divided into three main steps. While the first two steps are highly dependent on the input data and selection of appropriate modeling algorithms, the third step is an integration of results of the first two steps. Step 1 This step involves creating a high-resolution 3D geocel- lular reservoir model of the petrophysical properties and depositional environment of the reservoir. Initially existing geological, geophysical, and well data are interpreted and incorporated to create a structural framework and empty geocellular grid. Next, facies curves derived from the well logs are used to generate cell-wise distribution of lithofa- cies within the reservoir layers employing advanced geo- statistical algorithms. Lastly, spatial distribution model of petrophysical property in terms of porosity, permeability, saturation etc., is generated and constrained to the litho- facies model. This helps in capturing geological uncer- tainty while creating 3D petrophysical property distribution model for the reservoir. While carrying out 3D geocellular modeling, mainly two geostatistical simulation algorithms, viz. sequential Gaussian simulation (SGS) and Plurigaussian simulation (PGS) were employed. Sequential Gaussian simulation (SGS) (Deutsch and Journel 1998) is one of the popular geostatistical algorithms used for interpolation of continu- ous variable in reservoir modeling. It was originally introduced as a solution to the smoothing problem of interpolation by kriging. Sequential Gaussian simulation algorithm reproduces a globally structured two-point statistics using variogram model, whereas kriging provides a best estimate at each location with minimum error vari- ance and ignores estimates made at other locations. Thus, with regards to modeling of hydrocarbon reservoirs, SGS algorithms provide more relevant 3D reservoir models honoring global spatial variations in sample points. The implementation of sequential simulation consists of reproducing the desired spatial properties through the sequential use of conditional distributions (Arpat 2005; Dimitrakopoulos and Luo 2004). Consider a set of N ran- dom variables Z(ua), a = 1,…, N defined at N locations ua. The aim is to generate L joint realizations {z(l)(ua),_ = 1, …, N} with l = 1, …, L of the N, conditional to n available data and reproducing the properties of a given multivariate distribution. To achieve this goal, the N-point multivariate distribution is decomposed into a set of N univariate con- ditional distributions: 418 J Petrol Explor Prod Technol (2018) 8:417–429 123 F(u1, …, uN; z1, …, zN|(n)) = F (uN; zN|(n ? N - 1)) 9 F (uN - 1; zN - 1|(n ? N - 2)) 9 ��� 9 F(u2; z2|(n ? 1)) 9 F(u1; z1|(n)), Where F (uN; zN|(n ? N - 1)) = Prob {Z(uN) B zN|(n ? N - 1)} is the conditional cumulative distribution function of Z(uN) given the set of n original data values and the (N - 1) realizations z(l)(ua), a = 1, …, N - 1 of the previously simulated values. This decomposition allows generating a realization by sequentially visiting each node on the simulation grid. Multivariate Gaussian distribution model is the only model for which the above-mentioned decomposition is analyti- cally available (Deutsch 2002). The conditional probability distribution in a Gaussian model is determined using sim- ple kriging of each unknown value at any node uj given the original data and previously simulated values (n ? j - 1) (Deutsch and Journel 1998). This concept is very conve- niently used in SGS making it one of the popular algo- rithms in reservoir modeling. The whole process of SGS used in current study can be summarized in five steps: 1. Generate a random path through the grid nodes. 2. Visit the first node in that path and use kriging to estimate a mean and standard deviation at that node based on surrounding data (i.e., local conditional probability distribution). 3. Select at random a value from the distribution and set the node value to that number. 4. Include the newly simulated value as part of the conditioning data. 5. Repeat steps 1–4 until all grid nodes have a simulated value. Truncated Gaussian simulation (TGS) (Matheron et al. 1987) and Plurigaussian simulation (PGS) (Galli et al. 1994) technique are two popular methods of reservoir modeling workflow for simulating lithotype and facies of sedimentary rocks. TGS technique is best suitable for reservoirs where the lithotypes occur in a sequential order and it helps in defining the geometry and internal archi- tecture of the reservoir rocks (Armstrong et al. 2003). PGS is a natural extension of TGS and capable of handling complex geological situations when the lithotypes occur in complicated relationship rather than simple sequence. The underlying principle of TGS and PGS simulations is to set up one or more simulations of Gaussian random fields at every grid point and then use some rules to convert these Gaussian values into lithotype indicators (Allard et al. 2012). In case of PGS simulation, generally two Gaussian random fields are used to define the lithofacies structure. Armstrong et al. (2003) distinguishes four steps in a PGS approach as: 1. Choosing the rock-type rules depending on the types of relations among lithotypes. 2. Estimating the parameter values. Two key factors control PGS including the thresholds at which the different Gaussians are truncated and the variogram model of the underlying Gaussian variable. The proportion of each facies, the ‘‘rock-type’’ rule, and the correlation between the two underlying Gaussian random functions determines the thresholds. Once the thresholds have been worked out, the variograms and cross-variograms of the indicators are calculated experimentally. 3. Generating Gaussian values at the conditioning points. 4. Simulating Gaussian values at each grid node with a usual Gaussian simulation algorithm. Step 2 In this step, a basin modeling study is initiated covering the reservoir area. The basin model reconstructs the geologic history (i.e., burial history) by back-stripping the reservoir to its original depositional condition. Through this recon- struction, the mechanical compaction, pore pressures, effective stress, and porosity versus depth relationships are established for the entire reservoir. Well data in the form of log curves constitute the main input for basin modeling process in association with other basin-related information such as regional temperature and pressure profile. Basin modeling is a technique which allows to under- stand effect of physical and chemical processes in a sedi- mentary basin to generate hydrocarbon. Basin modeling techniques help to reconstruct the burial and temperature history of the basin through time and to understand source rock maturation and subsequently hydrocarbon expulsion and migration (Hanstschel and Kauerauf 2009; Peters 2009). The basin modeling process dynamically models changes in rock properties through geologic time by numerically solving coupled partial differential equations with moving boundaries on discretized temporal and spa- tial grids. This computation results in many outputs of different rock properties including properties such as vit- rinite reflectance, temperature, effective stress, and poros- ity. In the final stage, the basin models need to be calibrated with existing data such as vitrinite reflectance, temperature, and porosity (AlKawai 2014) to make it more realistic. Overall implementation of basin modeling process contains a wide range of mathematical algorithms and methods each of them appropriate for each ‘‘sub-model,’’ and detailed discussion of these goes beyond the scope of this paper. Additional information and overview of these numerical methods are provided by Press et al. (2002) and J Petrol Explor Prod Technol (2018) 8:417–429 419 123 Hantschel and Kauerauf (2009). In basin modeling, numerical solution of differential equations is fundamental requirement and most challenging, complex and costly process to accomplish. In most of the scenarios, tempera- ture and pressure are modeled with parabolic diffusion equations. Multi-dimensional differential equations with complicated geometries are usually solved with the finite element method, whereas finite differences are often used in special cases of approximately one-dimensional prob- lems, such as simplified crustal layer models (Hantschel and Kauerauf 2009). Step 3 It is the integration step where results from basin modeling study (i.e., distribution of mechanical properties of reser- voir rocks over geological time) are correlated with each cell of geocellular model using empirical relationships (Eberhert-Phillips et al. 1989; Ingram and Urai 1999; Horsrud 2001; Yarus and Carruthers 2014). Use of these empirical relationships provided a practical solution in discretizing geomechanical properties in the current study area. Other published evidences of use of these relation- ships in basin modeling process include Alkawai (2014), Szydlik et al. (2015), and Schneider et al. (1996). The output of this process is a discretized geomechanical model of the reservoir along with other petrophysical properties. Thus, the methodology adopted in this study is a true integration of results and parameters of both geocellular and basin models. Diagrammatic representation of the integrated approach is provided in Fig. 1. A step-by-step description of the implementation of the entire modeling process for appraisal of unconventional mature field in western Canadian sedimentary basin is explained in the subsequent section of the paper. Geological background of the study area The study area is located at the western boundary of Sas- katchewan, Canada, close to the Alberta border in Hoosier field (Fig. 2). Previous studies by White (1969) identified two oil pools in the field. The study area falls in the North Hoosier oil pool at the west central Saskatchewan. In this field, the Bakken shale formation is the important forma- tion due to its high hydrocarbon potential. The Bakken formation was deposited during the geological age of Late Devonian and Early Mississippian. It lies unconformably over the Big Valley formation (Upper Devonian) and is conformably overlain by the Madison group as shown in Fig. 3 (Zhang et al. 2016). Deposition of Bakken formation occurred through a series of onlap–offlap cycles during the Tamaroa sequence (Wheeler 1963). Time-equivalent shale units of Bakken shale include the Exshaw/Banff in the Alberta Basin, the Woodford shale in the Anadarko Basin, the Chattanooga in the Southern Appalachian Basin, and the Antrim in the Michigan Basin (Meissner 1978). The Madison Group which overlies the Bakken formation is a thick carbonate sequence of Mississippian age. The rocks of Madison were deposited in a generally shallow marine setting as indicated by richly fossiliferous rocks and favor for accumulation of petroleum resources (Porter 1955). It has been divided into three formations, viz. Lodgepole, Mission Canyon, and Charles. The oldest Lodgepole con- sisting limestone and dolomites conformably overlies the Bakken formation is a major hydrocarbon-producing horizon (Heck 1979). Abnormally high pore pressures in the Bakken shales have been attributed to hydrocarbon generation associated with the thermal anomaly in some parts of the basin (Price et al. 1984). Hydrocarbon pro- duction in the Bakken depends on accurate placement of horizontal wells with multistage fracture stimulations. The effective minimum horizontal stress is a primary control- ling factor for fracture growth. Thus, knowledge of distri- bution of elastic properties of the Bakken shales is very important to accurately determine the effective minimum horizontal stress (Havens 2012). In terms of wireline log data interpretation, Bakken shale is easily recognizable because of the strong contrast in lithology (Havens 2012). The upper and lower shales have unusually high gamma ray readings and high resistivity, while the middle member has a signature similar to clastic and carbonate rocks. Data collection and conditioning For the current study, geologic data on Hoosier field were sourced and collected from public domain sources (Al- berta Geological Survey; Government of Saskatchewan). The data gathered for the current study mainly comprised of well location information, interpreted formation tops demarcating geological boundaries, wireline logs, geo- logical maps, etc. No seismic data are available covering the study area to enable understanding of lateral conti- nuity of lithological properties within reservoir layers. Thus, facies logs interpreted and derived from available wireline log data were heavily relied upon for establishing spatial continuity of depositional model of the reservoir. Three main geological markers (formation tops) repre- senting three most important geological units, i.e., Man- nville, Detrital, and Bakken, were identified from the wireline log interpretations of 27 well data. Figure 3 shows a stratigraphic column with the formations of interest. 420 J Petrol Explor Prod Technol (2018) 8:417–429 123 Fig. 1 Diagrammatic representation of geocellular model and basin model integration Fig. 2 Location Map of study area showing the Hoosier field. Modified after O’Connell et al. (2000) J Petrol Explor Prod Technol (2018) 8:417–429 421 123 3D geocellular modeling The interpreted formation tops namely Mannville, Detrital, and Bakken were first used to generate depth structure grids (i.e., maps). These depth structure grids (Fig. 4a) serve as input to create the structural framework of the reservoir with top and base of the framework as Mannville and Bakken formations. Internal stratigraphic unit thick- nesses and their relationship to each other helped in establishing the layering schemes in the structural frame- work. The structural framework constitutes the basic grid skeleton for discretization of reservoir properties. The 3D structural framework is then divided into blocks of size 150 m 9 150 m with vertical layer thickness of 0.5 m (Fig. 4b, d). A detailed multi-well lithostratigraphic correlation was carried out to establish the continuity of subsurface lithology. Multivariate statistical analysis technique is employed to group the log curves into electrofacies groups. Each of the electrofacies classes has been assigned a rock class based on the sedimentological interpretation of the wells. Lithological classes thus assigned are correlated with each of the wells covering the entire study area. The assigned lithology classes formed the percentage of rock values (Fig. 4c) that would be used to fill in the grid during facies modeling. The well-wise rock interpretation is next combined with petrophysical information to establish its effect on porosity, permeability, and saturation distribution controlling dynamic fluid behavior. Log measurement data from 27 wells in terms of petrophysical properties are converted to point sets and then blocked to the grid. A point set is collections of generic points that has 3D spatial location information, i.e., easting (x), northing (y), and depth (z) and rock property measurement values such a porosity and permeability corresponding to each point locations. While converting a log curve to point set, first measured depth (MD) (z) values and corresponding rock property measurements from log curves (such as porosity) are captured in columns in a data file. Next, easting (x) and northing (y) information is cap- tured from well head and added to each data points as two separate columns. Thus, the resultant point set data file has series of columns starting with x, y, z and followed by rock property values. Generally, each of these data points is taken at the scale of actual log measurement, i.e., 0.5 feet (approximately 15 cm). After creation of the point set, the next step is to block the point set to the grid resolution (Shepherd 2009). Usually, the vertical dimensions of geo- cellular grids are larger than the vertical sampling interval of log curves (or the point set). As a result, each geocells passing through well locations carry multiple points. In order to use the point information for modeling and inter- polation, the multiple points within each cell need to convert to one point per grid cells using an averaging technique such as arithmetic, geometric, and harmonic for continuous properties. In case of discrete properties (such as facies or lithology), the averaging approach is ‘‘most of,’’ i.e., most commonly occurring (mode) lithology within the points occurring in a single geocells (Cannon 2015). This process is referred as blocking of point sets to the geocellular grid. Lithotype assignment and vertical proportion curves in terms of three main lithotypes, i.e., dolomitic sandstone, siltstone, and shale are generated. These lithotype proportion curves resulted into lithotype proportion maps (Fig. 5) which are then incorporated in generation of facies distribution model for the reservoir. Variogram models have been constructed to establish data stationarity and spatial correlation of input data points over the study area. Establishing stationarity is expedient for an implicit assumption of a geostatistical calculation in the current study. The absence of a nugget effect in the variograms showed that there is no randomness in the sampled well population. The spatial position of well locations produced a large spatial correlation distance for each of the variables. This correlation allows to progress Fig. 3 Stratigraphic column of study area 422 J Petrol Explor Prod Technol (2018) 8:417–429 123 with a geostatistical approach to extrapolate the measured data for the entire grid cells. Plurigaussian simulation (PGS) algorithm (Armstrong et al. 2003) was employed with two variogram models for creating 3D lithofacies distribution model. In PGS simu- lation, the controlling component is lithotype rule, which defines the depositional environment (Armstrong et al. 2003). In addition to lithotype rules, the lithotype proportion maps resulted from lithotype curves provided a background guideline for cell-wise facies distribution. PGS is one of the best techniques in geologically complex depositional environment modeling. Next, sequential Gaussian simulation (SGS) (Deutsch and Journel 1998; Deutsch 2002) was used for creating 3D spatial distribution model of petrophysical properties (such as porosity and permeability) employing variograms created for each of Fig. 4 a Depth structural map from tops, b 3D grid layers, c Location of wells in the study area and d 3D block configuration Fig. 5 Vertical proportion curve and lithotype proportion map showing vertical and areal distribution of rock types J Petrol Explor Prod Technol (2018) 8:417–429 423 123 these properties. While creating the porosity model, pre- viously created facies model was used to constrain or bias reducing geological uncertainty of modeling of petro- physical properties. Distribution of facies and associated porosity in the reservoir is shown in Fig. 6. Building 1D basin model The concept of basin modeling in general refers to inte- gration of geological, physical, and geochemical processes in a sedimentary basin to simulate sediment burial, com- paction, heat transfer, hydrocarbon generation, and migration and entrapment processes within it. Basin mod- eling reconstructs the pressure regime of the basin during sediment deposition and predicts zone of overpressure. This process also delineates hydrocarbon migration path- ways and identifies structural and stratigraphic traps. Basin modeling can be performed in 1D (one-dimensional), 2D (two-dimensional) or 3D (three-dimensional) to simulate the burial and thermal histories of the basin with increasing complexities. If the area is relatively small, such that thermal, pressure, and other effects are uniform over the study area, then a 1D basin model might be sufficient (Yarus and Carruthers 2014). If the basin development processes vary significantly over the area, a 3D study is warranted. In the current study, the necessity of 2D/3D basin model is not felt due to relatively uniform thermal and pressure effect on the basin-forming process in the study area. Comparison 1D basin model results from few places within the study area also revealed that important parameters controlling the basin-forming process do not vary much. Thus, the current study was restricted to con- struction of 1D basin model and uses the results for derivation of geomechanical properties. The 1D basin model uses rock and fluid information to evaluate source rocks and thermal histories of basin along a single direction, i.e., vertical. In this process, layers in the 1D model are sequentially back-striped to its original depositional thickness resulting into effective stress– porosity relationships for all lithologic layers. The lithologies are generally assigned to the model layers by zones using petrophysical logs and other geological infor- mation. These lithologies provide thermal conductivity and heat capacity for each layer. In the current study, the input parameters used in developing 1D basin model of the reservoir include: • Geologic framework in terms of layer thicknesses and ages • Rock properties (i.e., lithology) for each framework layers • Well data information such as temperature and pressure • Identification of source layers with assignment of kerogen types • Regional geothermal gradient of the basin Figure 7 shows the output of basin modeling process in terms of variation of elastic properties and pressure regim while depositing the sediments over geological periods. Modeling of the parameters is governed by the rates of pressure generation (e.g., compaction disequilibrium) and pressure dissipation (e.g., as controlled by permeability) during basin-forming process. Another important output of 1D basin modeling process is burial history or geohistory plot of the reservoir as shown in Fig. 8. Geohistory plots indicate how sediments were deposited over geological time in the basin providing thickness variation information of layers. In other words, it shows the time versus thickness relationship for all the layers including the layers that have been eroded or periods of non-deposition (Hiatus). Based on the lithology information of the 1D model inputs, depth–porosity relationship is also calculated while estab- lishing burial history. The depth–porosity relationship is then converted to calculate the effective stress regime and temperature profile of each layer through the geological time. The effective stress thus derived for each layer of the basin provides crucial information for calculation of geomechanical properties of the reservoir rock. Fig. 6 Distribution of facies (rock types) and porosity (facies-constrained) 424 J Petrol Explor Prod Technol (2018) 8:417–429 123 Fig. 7 Output of 1D basin modeling process in terms of reservoir elastic properties and pressure regim Fig. 8 Geohistory plot (burial history)—output from 1D basin modeling process J Petrol Explor Prod Technol (2018) 8:417–429 425 123 Model integration The final stage of the current integrated approach is to combine the result of classical 1D basin model (burial history) with geostatistical 3D geocellular model to gen- erate discretized geomechanical model of the study area. Lithology identifier of each cell from 3D facies model provides the link of this integration. Each element of burial history model has also been assigned a lithological iden- tifier to define (1) mechanical compaction, i.e., porosity versus depth relationships, (2) pore pressure relationships, i.e., porosity versus effective stress, (3) fluid transport properties, i.e., porosity versus permeability and capillary threshold pressure, and (4) matrix properties such as ther- mal conductivity. The lithology identifiers available in both the models enable populating layer-wise distribution of effective stress for each cell in the 3D geocellular grid (Fig. 9). The effective stress estimation is obtained during the process of burial history modeling of the sedimentary basin. The burial (or geo) history model provides infor- mation on rate of sedimentation including depth and age of deposition. While calculating effective stress regime of the basin, we assume that sedimentation and subsidence is the primary cause of pressure and stress; stresses generally increase with depth and rock stress generally interacts negatively with porosity and compaction (Hantschel and Kauerauf 2009). In 1D basin modeling approach, Terza- ghi’s (1923) effective stress concept is widely used as the fundamental relationship for effective stress modeling (Hubert and Rubey 1959). Effective stress ðr0Þ ¼ Total stress ðrÞ�Pore pressure pð Þ: ð1Þ In this relationship, lithostatic pressure (equaling the sediment weight) is considered as total stress and additional stresses due to compaction or extensional forces are neglected (Hantschel and Kauerauf 2009). Biot (1941) established the theory of poro-elasticity for rocks and introduced Biot’s coefficient in the definition of the effective stress (Schneider et al. 1996). r0 ¼ r � a p; where a is Biot’s coefficient: ð2Þ In these models, pore pressure of formation is related to incomplete mechanical compaction and a fixed relation between porosity reduction and sediment compaction is assumed (Hantschel and Kauerauf 2009). Studies by various other researchers established a number of empirical relationships of porosity versus depth versus compaction. These relationships show that changes in porosity with increasing depth is a function of lithology of the sediments. Thus, in basin modeling process these lithology-wise empirical relationships are used to satisfy the depth–porosity function and pore pressure calculation. The complete mathematical explanations of these relationships and numerical solutions to complex equations can be obtained from Hantschel and Kauerauf (2009) and Schneider et al. (1996). Once the effective stress distribution model is generated for the entire reservoir geocellular grid, rock mechanical properties such as Vp, Vs, UCS, BRI, OCR, Young’s modulus, Poison’s ratio, and bulk modulus are calculated at each cell locations of geocellular grid using empirical formula as provided below (Eberhert-Phillips et al. 1989; Ingram and Urai 1999; Horsrud 2001; Yarus and Car- ruthers 2014): Vp ¼ 5:77 � 6:94/ � 1:73 p C þ 0:446 Pe � e�16:7Pe � � ; ð3Þ Vs ¼ 3:70 � 4:94/ � 1:57 p C þ 0:361 Pe � e�16:7Pe � � ; ð4Þ Log UCS ¼ �6:36 þ 2:45 log 0:86 Vp � 1172ð Þ; ð5Þ OCR ¼ Max Pe=Pe; ð6Þ BRIT ¼ UCS=UCSNC; ð7Þ E ¼ qV2s 3V 2 p � 4V 2 s � � = V2p � V 2 s � �h i ; ð8Þ m ¼ 0:5 V2p � 2V 2 s � � = V2p � V 2 s � �h i ; ð9Þ K ¼ q V2p � 4=3V 2 s � � ; ð10Þ where Vp = P-wave velocity, Vs = S-wave velocity, / = Porosity, C = Clay content, Pe = Effective stress, UCS = Unconfined compressive strength, OCR = Over consolidation ratio, BRIT = Brittleness, NC = Uncon- fined, q = Density, E = Young’s modulus, m = Poisson’s ratio, K = Bulk modulus. For each of the calculations lithology, porosity and effective stress serve as main input parameter. Figure 9 displays a representative view of the geocellular grid with discretized geomechanical property derived from the inte- gration of result of basin and geocellular modeling. Results and discussion Well coverage in the study area is not adequate as revealed in Fig. 4c. The wells are mostly clustered in the central and north-eastern part. Furthermore, in the study area conven- tional cores were collected in only five wells out of 27 wells made available. Thus, physical measurement and interpolation of geomechanical properties of these core samples were not considered as reliable source of infor- mation for future drilling and investment decision. Addi- tionally, due to non-availability of 3D seismic data, 426 J Petrol Explor Prod Technol (2018) 8:417–429 123 interpolation of well-based geomechanical properties beyond well control could not be validated through seismic attribute correlation and integration as secondary input. As a result, attempts to use the conventional approach of using well core sample analysis-based geomechanical property interpolation in the study area provided inappropriate information for full field appraisal and development plan- ning. However, use of current geostatistics-based modeling approach explained in this study helped to nullify these restrictions and provided a full field geomechanical prop- erty distribution model with the integration of parameters from basin modeling approach. Such integrated models are highly flexible and appropriate for full field development with reliable distribution of geomechanical property in 3D space of unconventional reservoir. The model developed in the current study clearly indi- cates that elastic moduli decrease with increasing porosity; therefore, such models would aid in the understanding of the distribution of brittle versus ductile zones in the reservoir. Statistical analysis in terms of correlation coefficient of Young’s modulus and porosity reveals a value of - 0.86. Figure 10 displays a plot of Young’s modulus and porosity at K-level 23 of 3D reservoir grid. Fig. 9 Distribution of geomechanical properties in the final integrated reservoir grid Fig. 10 Plot showing correlation of Young’s modulus and porosity for K-level 23 of 3D reservoir grid J Petrol Explor Prod Technol (2018) 8:417–429 427 123 Based on the spatial distribution pattern of geomechanical properties, the lower part of the Bakken formation in the study area has been identified as zone of prime importance with favorable values of geomechanical properties to maximize hydrocarbon production. In the absence of 3D seismic coverage in the study area, such information aids operators tremendously for effective design of depletion strategies with minimal geological risks and associated costs. Furthermore, for a variety of other subsurface operations, including flow simulation, history matching, fluid typing, and hydrocarbon production, these models serve as main input. Understanding the burial history of the reservoir is a step forward in modeling unconventional resources (Yarus and Carruthers 2014). As part of full field development project, a realistic picture of the depositional history of the reser- voir through geological time provides a more solid foun- dation to build a knowledge base for economic forecasting and maximum return of investment. In the unconventional upstream oil and gas industry, spatial distribution modeling of rock mechanical properties of matured reservoir has always been a major challenge. However, 3D spatial dis- tribution of model of geomechanical properties created in the current study provides an easy way of knowing pre- vailing rock elastic properties and stress regimes of the reservoir rocks. Conclusions Field development planning, in the current study area, involves the use of horizontal wells with multistage hydraulic stimulation in shale sweet spots as the primary production strategies. Knowledge of in situ stress condi- tions and facies is therefore very important for these operations as it guides the stimulation design which will in turn control the fracture propagation. Earlier studies on Bakken formation clearly elaborated the use of velocities and strains measured on samples from the core which are not sufficient enough to optimize field development strat- egy using hydraulic fracturing. This is because small sample plugs or slabs may not be always sufficient to describe elastic properties for the entire reservoir. For instance, Vp measured from one core sample may not be very diagnostic of lithology as it cannot be used to map shales and sandstones across a reservoir with accuracy. The current integrated approach incorporates wireline logs in the form of hard data covering entire sediment sequence. This provides a true representation of the sub- surface geological condition of the reservoir even if there are no 3D seismic data available. In this approach, use of geostatistical techniques results into highly accurate 3D facies model resembling real geological condition of deposition of sediments. The geostatistical techniques via spatial covariance models of sample values helped in quantifying and mapping of geological heterogeneity in sediment deposition in unconventional reservoirs. Addi- tionally, conditional simulation algorithms used in geo- statistical modeling approach honor the sample values at location of measurements and produce minimum biasness of interpolation of rock properties from well data. Thus, results generated from final integrated model are not only reliable, but also capable of providing accurate local error of estimation (i.e., heterogeneity) within them. Similar studies made elsewhere (Yarus et al. 2016; Yarus and Yarus 2014; Yarus and Carruthers 2014) indicated how combining geostatistics with basin modeling techniques accurately identifies sweet spots in unconventional reser- voirs. Verification and comparison of results obtained from such integrated models with actual and blind wells pro- vided high accuracy and reliability score with minimum deviation of results. Such convincing results are also true validation of consistency of result of current integrated technique over conventional approach. Traditional modeling methods for identifying good and poor reservoir quality in shales have not been successful as reported by many authors. It is mainly due to lack of geospatial understanding of shales properties and difficulty of quantifying their heterogeneity. However, the current integrated geostatistics-based geomechanical models with higher degree of confidence show how the velocity and elastic moduli are decreasing with increasing porosity throughout the reservoir. As a result, in the current study area the concerned field team was able to identify target locations or sweet spots with high degree of accuracy within the shale for fracturing and completion. Further- more, these models can be used alongside other reservoir characterization methods in shale play workflows to aid in the planning of unconventional well bores for enhanced production. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. References Alberta Geological Survey. http://ags.aer.ca/. Accessed July 2016 AlKawai W (2014) Integrated basin modeling with seismic technol- ogy and rock physics. M.Sc. Thesis, Stanford University, Stanford Allard D, D’Or D, Biver P, Froidevaux R (2012) Non-parametric diagrams for plurigaussian simulations of lithofacies. In: 9th international geostatistical congress, Oslo, Norway, pp 11–15 428 J Petrol Explor Prod Technol (2018) 8:417–429 123 http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ http://ags.aer.ca/ Armstrong A, Galli AG, Loch GL, Geoffroy F, Eschard R (2003) Plurigaussian simulations in geosciences. Springer, Berlin, p 144 Arpat GB (2005) Sequential simulation with patterns. PhD Thesis. Stanford University, Stanford Biot MA (1941) General theory of three-dimensional consolidation. J Appl Phys 12:155–164 Cannon S (2015) Petrophysics: a practical guide. Wiley, London, pp 16–161 Deng H, Aguilera R, Settari T (2011) An integrated workflow for reservoir modeling and flow simulation of the nikanassin tight gas reservoir in the western Canada sedimentary basin. SPE Annual Technical Conference and Exhibition, Denver, Colorado, USA. SPE # 146953-MS. doi:10.2118/146953-MS Deutsch CV (2002) Geostatistical reservoir modeling. Oxford University Press, Oxford Deutsch CV, Journel AG (1998) GSLIB: geostatistical software library and user’s guide, 2nd edn. Oxford University Press, Oxford Dimitrakopoulos R, Luo X (2004) Generalized sequential gaussian simulation on group size m and screen-effect approximations for large field simulations. Math Geol 36(5):567–591 Eberhert-Phillips D, Han D-H, Zoback MD (1989) Empirical relationships among seismic velocity, effective pressure, poros- ity, and clay content in sandstone. J Geophys 54(1):82–89 Ferdous H, Bergman K, Qing H (2015) Viking lowstand deposits in west central Saskatchewan: depositional model for the reservoir units in Dodsland-Hoosier Area, Saskatchewan, Canada. AAPG Search and Discovery Article #51076 Galli A, Beucher H, Le Loc’h G, Doligez B (1994) The pros and cons of the truncated gaussian method. In: Armstrong M, Dowd PA (eds) Geostatistical Simulations. Kluwer, Dordrecht, pp 217–233 Gonzalez L, Nasreldin G, Rivero J, Welsh P, Aguilera R (2014) 3D modeling of multistage hydraulic fractures and two-way-cou- pling geomechanics/fluid-flow simulation of a horizontal well in the nikanassin tight gas formation, Western Canada sedimentary basin. SPE Reserv Eval Eng 17(2):257–270 Government of Saskatchewan. http://www.ir.gov.sk.ca/ Accessed July 2016 Hantschel T, Kauerauf AI (2009) Fundamentals of basin and petroleum systems modeling. Springer, Berlin Havens J (2012) Mechanical properties of the Bakken Formation. M.S. Thesis, Colorado School of Mines, Golden Heck TJ (1979) Depositional environments and diagenesis of the Mississippian Bottineau interval (Lodgepole) in North Dakota: University of North Dakota Unpublished M.S. Thesis Horsrud P (2001) Estimating mechanical properties of shale from empirical correlations. SPE Drilling and Completion conference, pp 68–73 (SPE # 00056017) Hubert MK, Rubey WM (1959) Role of fluid pressure in mechanics of overthrust faulting. Bull Geol Soc Am 70:115–166 Ingram GM, Urai JL (1999) Top-seal leakage through faults and fractures: the role of mudrock properties. Geol Soc Lond Spec Publ 158:125–135 Lavoie D, Séjourné S (2016) Geomechanical properties of the upper ordovician Macasty shale and its caprock, Anticosti Island: a regional evaluation for a promising tight oil play. AAPG Search and Discovery Article #80542 Matheron G, Beucher H, de Fouquet C, Galli A, Guerulot D, Ravenne C (1987) Conditional simulation of the geometry of fluvio- deltaic reservoirs. SPE annual technical conference and exhibi- tion, Dallas, Texas. SPE# 16753 Meissner FF (1978) Petroleum geology of the Bakken Formation, Williston Basin, North Dakota and Montana. In: Proceedings of 1978 Williston basin symposium, Montana Geological Society, pp 207–227 O’Connell SC, Wood JM, Robinson SP (2000) An ancient tidal sand ridge complex in the Bakken Formation (Fammenian-Tour- naisian) of southwestern Saskatchewan. GeoCanada 2000 core workshop, Calgary, Alberta. http://webapp1.dlib.indiana.edu/ virtual_disk_library/index.cgi/2870166/FID3366/PDF/1228.PDF Peters KE (2009) Getting started in Basin and petroleum system modeling. AAPG CD-ROM #16, AAPG Datapages Porter JW (1955) Madison complex of Southeastern Saskatchewan and Southwestern Manitoba. J Alberta Soc Petrol Geol 3(8):126–130 Press WH, Teukolsky SA, Vetterling WT, Flannery BP (2002) Numerical recipes in C??, 2nd edn. Cambridge University Press, Cambridge Price LC, Ging T, Daws T, Love A, Pawlewicz M, Anders D (1984) Organic metamorphism in the Mississippian-Devonian Bakken shale, North Dakota portion of the Williston Basin. In: Wood- ward J, Meisner FF, Clayton JL (eds) Hydrocarbon source rocks of the greater Rocky Mountain region. Rocky Mountains Association-Geologists, Denver, pp 83–133 Schneider F, Potdevin JL, Wolf S, Faille I (1996) Mechanical and chemical compaction model for sedimentary basin simulators. Tectonophysics 263:307–317 Shepherd M (2009) Oil field production geology: AAPG Memoir 91, pp 182–183 Szydlik T, Helgesen HK, Brevik I, Prisco GD, Clark SA, Leirfall OK, Duffaut K, Stadtler C, Cogan M (2015) Geophysical basin modeling: methodology and application in deepwater Gulf of Mexico. Interpretation 3(3):SZ49–SZ58 Terzaghi K (1923) Die Berechnung der durchlssigkeitsziffer des tones aus dem verlauf des hydrodynamishen spannungserscheinungen. Szgber Akad Wiss Vieu Math Naturwiss Klasse Ila 132(3–4):125–138 Tong A, Chi G, Pedersen PK (2005) Sequence stratigraphy and preliminary diagenetic study of the lower cretaceous viking formation, Hoosier Area, West-Central Saskatchewan. Summ Investig Sask Geol Surv 1:2005 Wheeler HE (1963) Post-Sauk and Pre-Absaroka Paleozoic strati- graphic patterns in North America. AAPG Bull 46:1497–1526 White WI (1969) Geology and petroleum accumulations of the North Hoosier area, west central Saskatchewan. Department of Mineral Resources, Geological Science Branch, Sedimentary Geology Division, Report No. 133 Yarus JM, Carruthers D (2014) Modeling unconventional resources with geostatistics and basin modeling techniques: characterizing sweet-spots using petrophysical, mechanical, and static fluid properties. AAPG Search and Discovery Article #41362 Yarus JM, Yarus JM (2014) Modeling with geostatistics, basin modeling techniques, E & P Magazine. http://www.epmag.com/ modeling-geostatistics-basin-modeling-techniques-758646#p= full Yarus JM, Mohsenian E, Yarus JM (2016) Combing geostatistical reservoir modeling with classical basin modeling techniques to identify sweet-spots in unconventional reservoirs. AAPG Data- pages/Search and Discovery Article #90254 Zhang L, Buatois LA, Knaust D (2016) Sedimentology, ichnology and sequence stratigraphy of the upper devonian-lower Missis- sippian Bakken formation in eastern Saskatchewan. Bull Can Pet Geol 64(3):415–437 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. J Petrol Explor Prod Technol (2018) 8:417–429 429 123 https://doi.org/10.2118/146953-MS http://www.ir.gov.sk.ca/ http://webapp1.dlib.indiana.edu/virtual_disk_library/index.cgi/2870166/FID3366/PDF/1228.PDF http://webapp1.dlib.indiana.edu/virtual_disk_library/index.cgi/2870166/FID3366/PDF/1228.PDF http://www.epmag.com/modeling-geostatistics-basin-modeling-techniques-758646%23p%3dfull http://www.epmag.com/modeling-geostatistics-basin-modeling-techniques-758646%23p%3dfull http://www.epmag.com/modeling-geostatistics-basin-modeling-techniques-758646%23p%3dfull An integrated approach to discretized 3D modeling of geomechanical properties for unconventional mature field appraisal in the western Canadian sedimentary basin Abstract Introduction Integrated modeling approach Step 1 Step 2 Step 3 Geological background of the study area Data collection and conditioning 3D geocellular modeling Building 1D basin model Model integration Results and discussion Conclusions Open Access References