Surface Drilling Data for Constrained Hydraulic Fracturing and Fast Reservoir Simulation of Unconventional Wells

Abstract

The objective is to present a new integrated workflow which leverages commonly available drilling data from multiple wells to build reservoir models to be used for designing and optimizing hydraulic fracture treatment and reservoir simulation. The use of surface drilling data provides valuable information along every wellbore. This information includes estimations of geomechanical logs, pore pressure, stresses, porosity and natural fractures. These rock properties may be used as a first approximation in a well-centric approach to geoengineer completions. Combining these logs from multiple wells into 3D reservoir models provides more value including using them in reservoir geomechanics, 3D planar hydraulic fracturing design and reservoir simulation. When using these 3D models and their results in a fast marching method simulator, the impact of the interference between wells can be estimated quickly while providing results like those derived with a classical reservoir simulator. Integrating surface drilling data with 3D reservoir models, hydraulic fracturing design and reservoir simulation into a single software platform results in a fast and constrained approach which allows for a more efficient management of unconventional wells.

1. Introduction

The recent findings of the Hydraulic Fracturing Test Site I industry consortium (HFTS 1) are well summarized in the September 2018 Journal of Petroleum Technology article titled “Real Fractured Rock is So Complex it’s Time for New Fracturing Models” [1]. 600 ft. of core taken in a hydraulically fractured Wolfcamp reservoir in the Permian Basin, USA showed a more complex reality than what is accounted for in most hydraulic fracturing design and analysis software. This includes the interaction between hydraulic and natural fractures, which is largely ignored or poorly accounted for in most software currently used to model hydraulic fractures and their resulting geometry. Rassenfos [1] emphasized in his summary of multiple recent publications describing HFTS 1 findings, that the “fracture height is overrated. While microseismic testing indicated that fractures grew up about a 1000ft, the height of the propped fractures- the fractures most likely to produce oil and gas was about 30ft.” Rassenfos [1] summary article discusses the important role played by the natural fractures but misses multiple other challenges facing the realistic modeling of hydraulic fracturing. Among the most noticeably and urgent challenges is the lack of data to better characterize the key inputs needed by any hydraulic fracturing modeling approach and the role of interfaces and their impact on vertical fracture growth.

The nature of the unconventional revolution, and ensuing extensive use of hydraulic fracturing, is the prevalent belief that the majority of the wells need to be drilled and stimulated in a “factory mode” where useful data such as wireline logs are not acquired at a statically significant rate. This philosophy, exhibited by many major unconventional players, has left a void of data in most major fields, greatly undermining optimization efforts necessary to economically produce said fields. A solution exists in the use of surface drilling data, which is available at every well. These drilling data are acquired by all drilling contractors around the world and are used qualitatively and quantitively during drilling operations. The surface drilling data include torque (T), rate of penetration (ROP), weight-on-bit (WOB). Since this drilling data is available at any old, current and future well and most operators dealing with unconventional reservoirs are not acquiring wireline logs at all the wells, the authors investigated the possibility to estimate pseudo-logs from surface drilling data. Using surface drilling data to infer rock properties, pore pressure, and stresses, comes with multiple challenges, which when overcome, open the door to improvements to the physics used in modeling hydraulic fracturing.

2. Newly discovered value in surface drilling data and mechanical specific energy (MSE)

Ouenes et al. [2] introduced the use of surface drilling data to simultaneously estimate the rock geomechanical properties, pore pressure, stresses, porosity and natural fractures needed to guide the steering of horizontal wells within the most frackable rock in real time, and additionally provide a completion design for optimal hydraulic fracturing when drilling is finished. The Mechanical Specific Energy (MSE) [3] computed from commonly available surface drilling data such as torque (T), rate of penetration (ROP), weight-on-bit (WOB) and bit diameter (D) has been widely used to improve drilling efficiency. The Mechanical Specific Energy is defined as

$$\mathrm{MSE}=4\left(\frac{\mathrm{WOB}}{\mathrm{\pi }{\mathrm{D}}^2}\right)+\left(\frac{480\mathrm{N}\times \mathrm{T}}{\mathrm{ROP}\times {\mathrm{D}}^2}\right)$$E1

All the components used in the MSE equation are commonly measured at the surface during drilling operations. Most horizontal shale wells currently being drilled use a drilling motor which requires the use of a different term for the Rotation Speed N*. When using a motor, the MSE requires the use of the formula given below in Eq. (2) where N is the rotational speed of the drill pipe, Kn is the mud motor speed to flow ratio and Q is the total mud flow rate.

$$N^{\ast }=N+K_n\times Q$$E2

Once the MSE is available, it can be used for multiple purposes including deriving geomechanical properties, pore pressure, stresses, porosity, and natural fracture indicators. To fully take advantage of the derived MSE, it should be further combined with other drilling information. Once MSE is available, the Unconfined Compressive Strength (UCS) can be derived from it in multiple ways. However, the UCS derived from the MSE needs to be corrected as the bit performance is also significantly influenced by the differential pressure.

2.1 Corrected MSE

Most of the recent MSE applications for completion optimization use surface drilling data which do not represent the MSE at the drill bit. The challenge posed by the use of surface drilling data consists of finding a way to eliminate costly and risky downhole equipment to measure the downhole MSE while ensuring accurate results. The solution for this challenge is to correct the surface drilling data by removing the frictional losses along the borehole. The Corrected Mechanical Specific Energy (CMSE), which is calculated in real time and uses surface drilling data, wellbore geometry, and drilling equipment parameters to estimate the friction losses along the drill string, was shown [4, 5, 6] to be a viable solution. This new technology uses advanced drilling and wellbore mechanics to estimate the multiple factors that create the frictional losses in real time and has been validated in multiple wells and basins.

Commonly available surface drilling data such as torque, weight on bit, rate of penetration and RPM are used as inputs to the model, along with the mud motor differential pressure and flow rates. Mud motor specifications, such as the maximum limits of differential pressure and flow rates, are also important inputs to normalize the computed Mechanical Specific Energy. When a Rotary Steerable Tool (RST) is used to steer the well, the Bottom Hole Assembly (BHA) information is a necessary input along with the wellbore survey to accurately perform the torque and drag analysis and to estimate the friction pressure losses along the wellbore. The friction pressure losses are then extracted from the surface MSE to compute the Corrected MSE.

The Drilling Efficiency (DE), which is the ratio of the energy required over energy spent in breaking a unit volume of the rock, is computed based on the CMSE and the Confined Compressive Strength (CCS) as shown below in Eq. (3)

$$\mathit{DE}=\frac{\mathrm{CCS}}{\mathbf{C}\mathrm{MSE}}$$E3

As shown in Eq. (4) CCS accounts for the typically increasing Unconfined Compressive Strength (UCS) rock strength with depth as well as the effects of the confining stresses ($∆p$) and angle of internal friction factor ($\theta$) applied on the rock. By correlating the MSE with CCS through the DE and by fitting a trendline on the computed DE data-set, the pore pressure can be estimated by accounting for the variations of DE data-set from the DE trendline. The fit of the DE trendline should be calibrated with pore pressure measurements from DFIT tests and the DE trendline should be updated accordingly.

$$\mathit{CCS}=\mathrm{UCS}+\mathrm{\Delta p}\left(\frac{1+sin\mathbf{\theta }}{\mathbf{1}\mathbf{-}\mathbf{\text{sin}}\mathbf{\text{θ}}}\right)$$E4

Once these friction losses are correctly estimated, they can be used to correct the MSE measured from surface drilling data which can be compared to measured downhole MSE. The principle of the predictive model is that torque and drag forces in a directional wellbore are primarily caused by sliding friction. Sliding friction force is calculated by multiplying the sidewall contact force with a friction coefficient. A lumped-parameter model provides the basis for the prediction of torque and drag. Both torque and drag are caused entirely by sliding friction forces that result from contact of the drill string with the wellbore. The frictional forces are subtracted from the surface MSE to accurately estimate the corrected MSE.

2.2 Comparing corrected MSE from surface drilling data to downhole measured MSE

Figure 1 shows an example from the Gulf of Mexico data set where Majidi et al. [7] emphasized the importance of using the downhole torque measurements as the torque term in the MSE often dominates the WOB term (see Eq. (1)). The discrepancy between MSE from uncorrected surface data and from the use of downhole measurements is clear in Figure 1. Using all the available drilling information, the frictional losses are estimated in order to compute the CMSE. The similarity between Majidi et al. [7] MSE from downhole measurements and the results computed from corrected surface drilling data is illustrated in Figure 1. Using the CMSE approach where surface drilling data are corrected for friction losses, results qualitatively comparable to measured downhole MSE can be derived. This opportunity leads to a reasonable estimation of mechanical rock properties at any well using commonly available surface drilling data thus circumventing the major problem of lack of well data in unconventional reservoirs.

Figure 1.

Comparison of CMSE (blue) derived from surface drilling data vs. downhole MSE (orange) measured downhole. Notice the difference of the CMSE values from the MSE derived from surface drilling data without correction (gray).

2.3 Estimating rock properties from CMSE

The next step is to leverage the estimation of CMSE and UCS to build a real-time wellbore geomechanical model. The CMSE is directly used as a proxy for UCS by finding a linear correlation of the CMSE to the average UCS values in the zone of interest. Velocity in rocks primarily depends on three factors namely porosity/effective pressure, saturating fluids and lithology/rock minerals. When focusing primarily on the lateral section, it is reasonable to assume that saturating fluids are fairly homogenous. Thus, the two contributing factors to acoustic and shear velocity become lithology and porosity/effective pressure which are used to estimate these velocities and the rock mechanical properties. For example, the Young’s Modulus (YM) can be derived from UCS using multiple available correlations based on different lithologies. Knowing the YM could lead to using other correlations to estimate the Poisson’s ratio (PR), Shear Modulus (G), Porosity (PHI), Fracture Index (FI), and rock brittleness (STRBRT). Majidi et al. [7] showed how the MSE could also be used to derive pore pressure. Using frictional faulting theory, with the UCS and the pore pressure, in-situ stresses can also be estimated. Figure 2(A) illustrates the common input surface drilling data and the resulting outputs that include stresses in Figure 2(B) and rock properties in Figure 2(C). Using these key rock properties as inputs, multiple other properties combining both rock properties and stresses can be derived and used in completion optimization.

Figure 2.

(A) Using the commonly found surface drilling data to estimate the (B) pore pressure and stresses and (C) key geomechanical logs, porosity and natural fractures along any wellbore.

3. Engineered completion using surface drilling-derived logs in every unconventional well

Production diagnostic tools have time and time again confirmed the variable performance of stimulated stages, prompting the need to geoengineer completions to adapt the treatment to the variable nature of the stress and rock properties along the wellbore. Fortunately, the use of surface drilling data leads to the estimation of key rock property logs, minimum stress and CMSE that could be used as a reference log to geoengineer the different hydraulic fracturing stages. The objective is to set the stages and the clusters in rock with similar fracture gradients, so they break at a common treating pressure. Figure 3 shows a well that has benefited from the use of the CMSE derived from surface drilling data and the variable reference log used to geoengineer the completion. When starting with a geometric design an initial cluster efficiency of 62% is computed along the wellbore. As shown in Figure 3 in the two red rectangles, minor changes in the position of the stages and their clusters could increase the cluster efficiency to 70%. Current applications of this technology to various wells have shown a consistent increase from 10 to 30% of cluster efficiency when using the logs derived from surface drilling data. Fiber optic measurements have confirmed the ability of these surface drilling derived logs to predict the cluster efficiency within a reasonable engineering error. While fiber optic measurements are extremely rare and very expensive, on the order of a million dollars, the logs derived from surface drilling are available at every well at an insignificant cost.

Figure 3.

Using a surface drilling derived reference log to geoengineer completions by moving stages to target similar rock, increasing cluster efficiency.

Given the density of unconventional wells, this newly available information provides the opportunity to go beyond the wellbore and start propagating this well-centric information into a 3D representation of the reservoir.

4. Well-based 3D modeling using geostatistics and artificial intelligence: propagating the surface drilling derived logs across the reservoir volume

In Sections 2 and 3, it is shown how the rock mechanical properties, pore pressure and stresses are extracted along the lateral section of the wellbore and are used as reference properties to engineer the completion to improve the overall efficiency. In this section, these properties are propagated in 3D space to accurately characterize the reservoir, so that these inputs can be fed into the hydraulic fracturing design and reservoir simulation workflows explained in the subsequent sections. The large number of wells drilled in unconventional assets combined with the estimation of critical logs at all the wells from surface drilling data provides the opportunity to propagate the well information into a 3D reservoir model. Since many companies do not have seismic on their acreage or for budgetary reasons do not plan to license the existing seismic, these multiple logs derived at all of the wells allow the construction of reliable 3D reservoir models. These 3D models could be estimated in a stratigraphic framework over a large area that encompasses many wells. In such cases, geostatistics could be used to estimate the distribution of gamma ray, porosity, Young’s Modulus, Poisson’s ratio and shear modulus. However, the pore pressure, minimum stresses and natural fracture are more complex continuous properties that need to be estimated with neural networks [8, 9] and other artificial intelligence tools able to capture the complex geology that control their variability.

One major reason for propagating these rock properties in 3D is to provide that information to a 3D planar hydraulic fracturing simulator as well as to geomechanical software. To achieve this goal, all the wells are used together in a large reservoir grid to create the 3D models from which smaller well grids (Figure 4A) will be extracted around a well or a pad. With this approach, all the available well data will be used to improve the 3D distribution of the key properties needed for the 3D planar hydraulic fracturing simulator. The other benefit of these derived 3D models is the estimation of the stress gradients resulting from the interaction between regional stresses and the three sources of stress perturbation created by the local geology: (1) variable geomechanical properties, (2) pore pressure and (3) natural fractures all available from the extrapolation in 3D of the logs derived from surface drilling data. Because of their importance in estimating these stress gradients, the modeling of the natural fractures requires some particular attention.

Figure 4.

(A) A large stratigraphic 3D geocellular grid is built from all the available wells to propagate 3D reservoir models. A smaller, higher resolution grid is extracted around the well that provides to the 3D planar hydraulic fracturing simulator: (B) Young’s Modulus, (C) Poisson’s ratio, and (D) unconfined compressive strength (UCS). (E) A cross section in the well grid of the minimum stress and (F) Poisson’s ratio honoring the lateral and vertical variability captured by the 3D models.”

5. The importance of 3D natural fracture models

Natural fractures have a significant impact in unconventional reservoirs, yet they are rarely accounted for in most physical modeling related to hydraulic fracturing. Natural fractures could have a positive impact as they create additional surface contact during hydraulic fracturing which is commonly referred as fracture complexity [10]. This can be predicted with geomechanical modeling and validated with microseismic response [11, 12]. The contribution of the natural fractures could also be negative by creating direct links to water bearing faults [13] or by creating frac hits [14, 15], through poroelastic effects, that will often damage the production from child and parent wells. Given their importance, a predictive model that provides the 3D distribution of these natural fractures is a critical input for any model trying to predict the outcome of hydraulic fracturing. However, finding a 3D distribution of the natural fractures has two major challenges: how to define the natural fractures at the wells given the rare occurrence of core or image logs in unconventional wells, and how to distribute the limited well data in the 3D reservoir to create a predictive model.

One of the motivations behind the use of surface drilling data is to be able to extract a fracture indicator that can be used to enrich the poor and limited statistics of natural fractures indicators found at wells and using these to build a 3D natural fracture model. Since the natural fractures are not a result of a depositional phenomenon only, their prediction is very different from mapping a more conventional property such as reservoir porosity. Having few limited wells with core or image logs will likely not provide the full statistics of these natural fractures. In other words, most statistical methods such as Discrete Fracture Networks (DFN) may not have the proper statistics from the wells to make any reasonable prediction of their 3D distribution. Attempts were made to reduce this problem by constraining the DFN with a continuous property [16] to guide the statistical distribution but other issues made the use of DFN in natural fracture modeling very challenging. Among these challenges include the dramatic variations found in the upscaled properties [17] needed for additional use of the DFN in engineering applications.

To avoid all the issues related to the use of a DFN, the Continuous Fracture Modeling (CFM) approach was developed to create validated predictive models of natural fractures [8, 9]. The CFM approach takes full advantage of the surface drilling derived fracture index or even the limited statistics that can be found in any natural fracture proxy, image log or core. The CFM approach honors structural geology concepts and focuses on the drivers that influence the presence of natural fractures. For example, the density of natural fractures at a given point in the reservoir does not depend on poorly sampled statistics of various fracture sets measured through limited wireline data, but on the volumetric distribution and interaction of lithology, structural settings and distance to faults, porosity, and many other reservoir properties that create the resulting natural fractures. These reservoir properties commonly called natural fracture drivers could all be estimated directly or indirectly through geologic modeling and seismic processes that involve seismic inversion, spectral decomposition and volumetric curvatures [9] when the data is available. Since the relationship between the natural fracture drivers and the limited natural fractures available at the wells is complex, artificial intelligence [8] is used not only to retrieve any existing and potential correlations that honors the limited statistics but also all the structural geology concepts. With this approach, extrapolation beyond the limited statistics is possible and has been successfully applied during the last three decades to various problems requiring an accurate description on where the natural fractures are. Among these, are problems in geomechanics such as interactions between hydraulic and natural fractures. This interaction could be better understood if studied in a decoupled way; where the natural fractures role in altering the regional tectonic stress and their impact on the lateral propagation of the hydraulic fracture is separated from the effects of natural fractures during the vertical propagation of the hydraulic fractures.

6. Constraining the hydraulic fracture propagation in the horizontal direction

A major shortcoming in most of the current hydraulic fracture simulators is the assumption that shale reservoirs at the scale of the wellbore are subjected to a homogeneous stress environment. Hence, hydraulic fracture stages were designed based on a constant stress field, in terms of magnitude and orientation. Unfortunately, lateral stress gradients, and their effects on microseismicity, have evidenced a more complex situation. The origins of these lateral stress gradients are numerous and include variation of the rock geomechanical properties, pressure depletion around existing wells, and proximity to faults and their associated natural fracture systems.

A decoupled approach using a plane strain framework to capture the lateral stress gradients was used by Aimene and Ouenes [11]. Their geomechanical modeling uses as input the three key factors affecting the lateral stress gradients: rock elastic properties, reservoir pressure, and natural fractures. The elastic properties and reservoir pressure models derived in previous sections from surface drilling data are used as inputs for the geomechanical model. The model uses explicit fractures to describe the distribution of the natural fractures then simulates the proper initial stress conditions resulting from the various sources of stress variability followed by the simulation of hydraulic fracturing in this heterogeneous stress medium. Since microseismic data is limited to only a few wells, the geomechanical approach used to capture the lateral stress gradients must be able to predict microseismicity rather than use it as calibration. The resulting geomechanical simulation predicts the differential stress, stress rotations and strain which serves as a reasonable proxy for microseismic events for validation as shown in Figure 5.

Figure 5.

Differential stress (A) and strain (C) validated with microseismicity (B) and the resulting geomechanically constrained hydraulic fractures (D).

This geomechanical approach combines the advantages of the particle-based numerical method, Material Point Method (MPM), a meshless numerical method, and the CFM approach to solve for a general continuum mechanics problem where all reservoir realities (natural fractures, variable rock properties and reservoir pressure heterogeneity), can be accounted for when estimating the stress field prior to stimulation and its subsequent perturbation during hydraulic fracturing. A direct benefit of this geomechanical approach is the quick estimation of the differential stress.

6.1 Estimation of differential stress: geomechanics vs. surface drilling data

The first result of the reservoir geomechanics approach [11] is the differential stress (Figure 6) which can be used as shown in Paryani et al. [18] to geoengineer completions. The advantage of using differential stress for geoengineering completions is the ability to consider the complex geology beyond the wellbore. In other words, well-centric approaches such as the one relying entirely on using a reference log derived from surface drilling data, are approximations that work only if the geology is not highly variable around the considered well. When the geology is variable with significant variability of the geomechanical properties, natural fractures, and pore pressure, then the best approach is to use the derived 3D models as input in the reservoir geomechanics approach [11] to estimate the differential stress.

Figure 6.

(A) Interpreted faults used as input in the reservoir geomechanics that estimates (B) the differential stress and the lateral stress gradients needed to geoengineer the stages.

Figure 6 shows a two well pad where faults (Figure 6A) were interpreted from multiple wells and used as input in the geomechanical approach [11] to compute the differential stress (Figure 6B). Since the analysis of surface drilling data was available in both wells, the differential stress was also estimated by using only that limited information. The computed differential stress from surface drilling data (Figure 7, right track) is compared to the one derived from the reservoir geomechanical simulation (Figure 7, left track) extracted along the wellbore from the differential stress distribution shown in Figure 6B. This comparison shows very strong similarities between the differential stress derived from full reservoir geomechanics (Figure 7, left track) with the one derived from the well centric approach based only on surface drilling data (Figure 7, right track). Both curves indicate the same zones for high differential stress zones where engineered completions are required to overcome earth resistance to hydraulic fracturing. The engineered completion could adjust the pumping parameters, stage length and number of clusters according to the derived differential stress with the objective of pumping bigger stage lengths in areas of low differential stress and vice versa. Areas of low differential stress will promote complex fracturing whereas areas of high differential stresses will result in planar hydraulic fractures with lower cluster efficiency as shown in Figure 8. The resulting engineered completion can be derived within a few hours of the well reaching Total Depth (TD) which illustrates the benefits of using surface drilling data if no other information is available.

Figure 7.

(Left) Differential stress derived from reservoir geomechanics results shown in Figure 6B. (Right) differential stress derived from a well centric approach using only surface drilling data.

Figure 8.

Engineered completion using the differential stress as a reference log to adjust stage length and number of clusters.

6.2 Estimation of stress rotation

The knowledge of differential stress is important to predict the ability of creating fracture complexity and increased surface contact. Areas with higher differential stress will produce highly anisotropic hydraulic fractures with reduced surface contact. This problem will be further complicated if the maximum horizontal stress S_Hmax direction is not locally perpendicular to the wellbore, leading to undesirable parallel hydraulic fractures. Hence the need to estimate both the differential stress and the local direction of the maximum stress around the wellbore.

A good example that can illustrate this critical issue is the Grisham fault in the Permian basin, USA where stress rotates by up to 90°. Figure 9 shows the public domain faults [19, 20] used as input in the geomechanical simulation [11] to estimate the stress orientation.

Figure 9.

Stress orientation around the Grisham fault showing distinct behaviors of the stress field orientation north and south of the fault, and also along strike of the fault.

When the Woodford faults are subjected to a dominant E-W tectonic stress, rotations in the maximum horizontal stress direction arise as shown in Figure 9. It is important to note that faults which are oriented parallel or perpendicular to in-situ maximum stress direction, such as the N-S trending fault, cause little perturbation and critically stressed faults (roughly 30–60° from local maximum stress direction) cause large perturbations. While much of the basin is still subject to a S_Hmax within 10° of the input orientation, several areas evidence large deflections from this input orientation. These deflections can be further refined if using seismic data to define the local faults. Given the variability of the differential stress and the direction of the maximum stress direction, both captured with the plane strain modeling the next step is to evaluate the actual strain resulting from the hydraulic fracturing.

6.3 Estimation of strain for laterally constrained 3D planar hydraulic fractures

The 2D plane strain MPM modeling provides valuable stimulated reservoir volume (SRV) information by modeling the effect of a large increase in the stress around the wellbore and its distribution throughout the reservoir volume and interaction with fractures and faults as well as accounting for any variable geomechanical properties and pore pressure of the rock. This is achieved numerically by applying a large pressure on a hydraulic fracture plane with a given length varying between 100 and 200 ft. which is used to model the effects of the pumping pressure in the reservoir. The real surface contact available to the fluid to apply its pressure and create a stress front is much larger than the numerical hydraulic fracture assumed to be around 150 ft. Thus, the pressure applied to this limited surface must be higher than the pumping pressure and is approximately, in most realistic unconventional wells, about 2.5 times the minimum stress value. Since this stress can be modeled with a dynamic simulation, the pressure applied to the hydraulic fractures can be applied sequentially, in parallel, or in a zipper mode. This ability to simulate the sequence of hydraulic fracturing allows the proper representation of stress shadow effects between stages as well as those seen between wells. These stress shadowing effects are considered along with the complex geology present between the stages and wells. For each hydraulic fracturing sequence, the resulting strain will be able to provide useful indication on the resulting SRV as shown in Figure 10.

Figure 10.

(A) Equivalent fracture model (EFM) derived by CFM using only surface drilling derived fracture indicators logs. The EFM is used as input in the reservoir geomechanics that provides the initial perturbed stress field and the subsequent (B) asymmetric strain resulting from the hydraulic fracturing of the wells and (C) comparison to microseismic events (note: the heel section of the wells was not monitored due to operations). (D) and (E) The envelope of strain provides the gross geomechanical half lengths which provide the lateral stress gradients needed to constrain the 3D planar hydraulic fracturing simulator.

One simple way to account for the lateral stress gradients captured by the geomechanical simulation, is to estimate the geomechanical half lengths (Figure 5C) from an interpreted envelope of the strain (Figure 10B) that could represent a proxy for the SRV. These interpreted asymmetric geomechanical half lengths are used at each cluster or stage as a constraint in a 3D planar hydraulic fracture design. It is important to reemphasize that the use of the planar representation of the hydraulic fractures is not an indication that the hydraulic fractures are indeed planar but a simple mathematical discretization of an SRV estimated by the full geomechanical simulation.

Having a constraint in the lateral direction is very helpful for a better estimation of the fracture height when using a 3D planar hydraulic fracturing approach. In this model, the vertical fracture growth occurs in the simplified world of perfect interfaces where debonding does not occur in a layered anisotropic rock. Unfortunately, the fracture growth does not depend only on the lateral stress but also on the geologic nature of the laminations and the characteristics of their interfaces which could be weak and could shear and consume hydraulic fracturing energy thus reducing the hydraulic fracture height. Since we have successfully estimated and validated with microseismic data the lateral stress gradients estimated with 2D plane strain MPM model, this information can be used as an input in a 2D vertical problem where we will focus on the geologic factors affecting the vertical fracture growth.

7. Constraining the hydraulic fracture propagation in the vertical direction in the presence of weak interfaces and natural fractures

Interfaces are among the geological features that are known to have an impact on the vertical propagation pattern and the final fracture height in unconventional reservoirs. Interfaces limit adjacent lithologies with similar or contrasting properties. It could be a material or just contact between two adjacent lithologies. Typical material thickness is between 1 and 500 mm for volcanic ash layers, and μm to mm thickness for mineralized veins, highly or partially mineral filled fractures, organic matters layers in the form of bitumen lubricating film or kerogen coating the surface of the interface [21] and bentonite layers [22] than could vary in thickness up to a few centimeters.

Interface mechanical properties can either be strong or weak and can be further weakened by tectonic deformations. They are a source of displacement discontinuities and delamination and are associated with fracture propagation behaviors like kinking, offsets, bifurcation, stepping over and termination. This suggests that displacement continuity hypothesis that are used in many hydraulic fracturing models where contact mechanics between layers with stick conditions (no sliding) is used to model the interface between layers [23, 24], is not valid for modeling real interface effects. In fact, using stick condition within a stratified structure could not account for sliding and de-cohesion between the layers due to hydraulic fracture pressurization. Thus, the need to use a proper interface model that (i) accounts for the displacement discontinuity between layers and (ii) allows the decohesion and interface delamination at the interface.

Explicitly modeling the interfaces presents some modeling challenges: (i) interfaces could be very thin layers that require the deployment of high-resolution model, (ii) interface mechanical properties are difficult to access given the well logs limitation in detecting them.

Aimene et al. [25] introduced the combined use of Anisotropic Damage Mechanics (ADaM) model and interface models in MPM to model the effects of interfaces in 2D and 3D hydraulic fractures problems. Multiple interface modeling tools were deployed [25] to achieve a better understanding of the impact of the interface properties on the fracture propagation initiation, growth and path. These tools include the Coulomb frictional contact and imperfect interface models. In the Coulomb frictional contact model, the contact between materials is modeled by setting the tangential traction S proportional to the normal force N at the interface by using the friction coefficient μ, i.e. S = min (μ N, S_stick), where S_stick is tangential traction required for the interfaces to stick with zero discontinuity. In other words, the materials stick until the tangential forces required by sticking exceeds the frictional term after which the interface is modeled by frictional sliding.

In the imperfect interface model, interfaces are represented implicitly by their overall phenomenological effects described with a much smaller number of interfacial parameters [26]. In this context of imperfect interface model, surface drilling derived logs provide once more the information needed to detect major interfaces that could be responsible for loss of energy during hydraulic fracturing. These major interfaces could be inferred from the resulting Young’s Modulus or/and Poisson’s ratio computed from the surface drilling data and the CMSE. For example, Figure 2C shows an example of a Poisson’s ratio pseudo-log estimated from surface drilling data where multiple spikes of high values indicate changes of lithology occurring where the well crosses a geologic interface. To illustrate the impact of interfaces and natural fractures on the fracture height, we will consider a 2D and a 3D case.

7.1 2D effect of weak interface on hydraulic fracture height

A specimen, under high vertical overburden stress, made of three layers alternating soft and stiff rock as estimated from the surface drilling estimation of Poisson’s ratio and Young’s Modulus is used to highlight the potential of a hydraulic fracture to develop a step-over behavior. We consider two cases: case 1 has perfect interfaces and case 2 has weak interfaces. To illustrate the step-over phenomenon, the two cases include a flaw (small fracture) located at the interface, close to the vertical path from the injection point (Figure 11).

Figure 11.

Fracture propagation in the presence of a natural fracture flaw (in red) close to the vertical path from the injection point in specimen-2. Stiff layer (dark blue) bounded by soft layers (light blue) in the presence of a (top) perfect interfaces (case1, top), and weak interfaces (case 2, bottom). Notice the step-over occurred in the weak interface but not in the perfect interface. The animated figure can be seen in Video 1 available from (can be viewed at) https://youtu.be/oqDx96YXSvQ

Figure 11 shows that for case 1 with a strong interface, the vertical fracture growth was insensitive to the presence of the flaw and the fracture propagated in the direction of the applied stress giving rise to symmetric fracture half-heights (Figure 11, top). However, for case 2 with the weak interface and a flaw near the injection point, the hydraulic fracture was first arrested by the weak interfaces, and then stepping over occurred when the fracture reached the interface (Figure 11, bottom). The fracture height is much smaller in the case 2 with weak interfaces and a strong asymmetric fracture height was developed as a result of high shear in the weaker interfaces.

The weak interfaces gave rise to an asymmetric fracture half-height, where the flaw promoted the propagation toward the direction of its location. These flaws, which are mainly bed-bounded natural fractures, could affect the propagation path of the hydraulic fracture and generate asymmetric hydraulic fracture half-heights or arrest the parent fracture and promote the secondary child fracture. It is this complex geologic reality that makes current hydraulic fracturing simulators inadequate to capture the propped frac height. This geologic complexity gets more complicated as we consider actual 3D situations.

7.2 3D effects of dipping fractures planes on fracture geometry

There are multiple field conditions that cannot be modeled with 2D approximations, and full 3D modeling is needed. This is the case in a strike-slip stress regime, when dealing with natural fracture planes that are not vertical, or when stimulating with helical perforation, etc. In these cases, full 3D modeling tools are essential to accurately reproduce the 3D fracture propagation mechanisms that will lead to the correct fracture height. To illustrate a full 3D analysis, a 3D laboratory test in [27] was considered. Briefly, the laboratory specimen was made with a real reservoir sandstone. The specimen, subjected to the strike-slip stress regime, contains 3 natural fractures with a 40° strike relative to S_Hmax and 75° dip. In the experiment, the pressurization was accomplished by eight perforations, four in each side, parallel to the strike of the fractures. The natural fractures are not mineralized, so they were modeled with the Coulomb frictional contact law with μ = 0.85 according to the laboratory experiment. Figure 12 shows the experimental results vs. the 3D MPM numerical model. The geomechanical modeling tool was able to reproduce the main features, especially the turning of the hydraulic fracture and arresting of the fracture by the nearest natural fractures. No crossing of natural fracture was observed neither in the numerical results nor in the experiment. This result highlights the ability of the combined use of the fully 3D ADaM model with the interface modeling tools in capturing the interaction between hydraulic and natural fractures in field conditions requiring full 3D geomechanical modeling tools.

Figure 12.

(a) Experimental setup and laser scan showing final 3D fracture geometry (green) arrest against 2 natural fracture planes in [27] laboratory test (b) 3D ADaM MPM result showing the fracture geometry and its arrest against the dipping fracture planes represented as weakly bonded interface. Animated version of the 3D MPM result can be seen in Video 2 available from (can be viewed at) https://youtu.be/jdfDAM2qi-8.

The results from these decoupled processes, that attempt to quantify the geomechanical impact of geologic characteristics, can be used to constrain a 3D grid based planar hydraulic fracturing simulator that will include proppant transport.

8. Constrained 3D grid-based planar hydraulic fracturing simulator

Various hydraulic fracturing design scenarios can be analyzed starting from the geometric and engineered designs extracted in Section 3. The 3D reservoir models generated in Section 4 are fed in the hydraulic fracturing design model along with the various stress gradients in the lateral and vertical direction extracted from the differential stress model in Section 6. In the presence of stress gradients created by natural fractures, variable geomechanical properties and depleted reservoirs, the characteristics of the formation would be different on either side of the wellbore. This asymmetry and its correct quantification are important for an optimal well spacing of unconventional reservoirs. Fischer et al. [28] proposed a hydraulic fracture model (Eq. (5)) that explains the relative length change in two opposite wings of the hydraulic fracture.

$${\mathit{a}}_{\mathbf{1}}\left(\mathit{t}\right)=\frac{{\mathit{p}}_{\mathbf{0}}^{\mathit{net}}}{\mathit{g}}\left(\mathbf{1}-\sqrt{\mathbf{1}-\mathbf{2}\left(\frac{\mathit{g}}{{\mathit{p}}_{\mathbf{0}}^{\mathit{net}}}\right){\mathit{a}}_{\mathbf{2}}\left(\mathit{t}\right)-{\left(\frac{\mathit{g}}{{\mathit{p}}_{\mathbf{0}}^{\mathit{net}}}\right)}^{\mathbf{2}}{\mathit{a}}_{\mathbf{2}}{\left(\mathit{t}\right)}^{\mathbf{2}}}\right)$$E5

The model is based on the lateral change in stress that would result in preferential growth of the hydraulic fracture in the direction of the decreasing confining stress. However, the model only explains the relativity of the asymmetric behavior in presence of the stress gradient (g) and does not specify the absolute length or width of the fracture. The relationship between the shorter wing a₁ (t) and longer wing a₂ (t) is given by Eq. (5) where ${\mathit{p}}_{\mathbf{0}}^{\mathit{net}}$ is the initial net pressure.

Using this concept, a 3D grid based planar hydraulic fracture model that includes proppant transport is developed. The unique feature of the hydraulic fracturing model is that it combines both analytical and numerical formulations. The effects of stress field change on the relative growth of the fracture is estimated by iterating for the optimum fracture height based on the amount of proppant available. The ability to develop a semi-analytical asymmetric fracture model that solves for the optimum fracture height and lengths is made possible by using the constraints of the geomechanical half lengths derived from the strain map and the estimated asymmetric fracture height derived from the ADaM geomechanical simulation.

The net pressure which is the difference between the fluid pressure and the minimum horizontal stress or the closure stress determines the initiation and propagation of a hydraulic fracture. The effect of stress gradient, along the fracture length is incorporated in the fundamental pressure balance equation at the fracture tip which determines the growth of the fracture. The fluid flow in the fracture is computed numerically. Using the relation of velocity of the fracture fluid and the fracture length, a time dependent solution is achieved which forms the basis of the semi-analytical model. All the rock properties and stresses are input in the hydraulic fracture simulator as 3D models as shown in Figure 4.

Using all these constraints as inputs in the 3D planar hydraulic fracture simulator, the pressure monitored during the actual fracturing treatment is easily matched (Figure 13A) by altering only the pipe and perf friction and the leak off coefficient which depends on the input porosity or natural fracture model. The resulting frac geometry at one stage (Figure 13C) or along the entire wellbore (Figure 13B) shows the major lateral and vertical dimension and conductivity variations owing to the variable nature of the rock properties captured by the surface drilling data. With this result at each well, we have all what is needed for the reservoir simulation which is needed to compute the Estimated Ultimate Recovery (EUR) and the resulting asymmetric depletion, a fundamental information needed for the planning of future wells.

Figure 13.

(A) Pressure match at a stage and resulting (B) complex fracture geometry and conductivity along the wellbore with major lateral and vertical variations due to the variable nature of the rock properties captured by the surface drilling data along the well and at a single stage (C).

9. Estimating well performances using reservoir simulation

For unconventional reservoirs undergoing hydraulic fracturing the long-term performance of a well could be represented by its Estimated Ultimate Recovery (EUR) which will also affect the extent of the depleted zone. The EUR and the size of the depleted zone will be major inputs needed for planning the development of an unconventional reservoir. This critical information requires the use of dynamic fluid flow reservoir simulation. In unconventional reservoir simulation, in addition to all the usual matrix properties required as input, the properties of the stimulated zone around the wellbore are necessary to capture the effects of the hydraulic fracturing and the development of the SRV. Multiple software tools and ways to provide the necessary input have been used. Despite major progress made recently in reservoir simulation, it remains a time-consuming task which prompted the need to develop a faster approach using the fast marching method (FMM). Both approaches are illustrated and compared with a field example.

9.1 Fast marching method (FMM)

Classical reservoir simulation using finite difference and finite volume have been widely used to simulate fluid flow in conventional oil and gas reservoirs. With the advent of unconventional reservoirs, these widely used classical reservoir simulation tools have been adapted to the particular nature of permeability generated through hydraulic fracturing. Current state of the art and the challenges associated with unconventional reservoirs are discussed in other books [29]. In this chapter we address the issue of computation time and the need for unconventional reservoir simulators to provide an estimate of the EUR and pressure depletion in the fastest possible way to enable the engineers and geoscientist to compare multiple development scenarios very quickly. In other words, given the nature of the unconventional process and its fast pace, how we can trade a reduction in accuracy for a much faster reservoir simulation?

Multiple efforts have been made to reduce computation time in reservoir simulators by replacing the flow equations with a proxy model based on neural networks [30] or response surfaces and experimental design [31]. Other efforts include the use of fast front-tracking techniques using streamlines [32] and, more recently, using a fast marching method or FMM [33]. We use the FMM in our approach of modeling unconventional reservoirs for the 3D estimation of pressure depletion.

The FMM is a front tracking algorithm. It has been applied to wave propagation, and medical imaging [34] problems. For the subsurface porous media flow, the pressure diffusivity equation can be simplified to an Eikonal equation and solved using the FMM to obtain the diffusive time of flight contours which is a proxy for the pressure depletion time [35]. The simulation is very fast compared to a finite difference simulator thus its valuable application to the fast-paced world of unconventional reservoirs. This speed is gained through a loss of accuracy which we could estimate by comparing the resulting pressure to those derived in a finite difference simulator.

9.2 Classical finite difference reservoir simulation

A dynamic model using a finite difference compositional reservoir simulator is built around a two well pad by using previously derived 3D properties (Figure 4) as well as the final hydraulic fracture geometry and its resulting conductivity (Figure 13B). Another derived key input is the inter-well permeability resulting from the interaction between the hydraulic and natural fractures and captured by the strain resulting from the geomechanical simulation using the MPM 2D plane strain framework. This strain is converted into an effective inter-well permeability (Figure 14B) using the approach shown in [36] where a calibration factor that relates strain to the matrix effective permeability is estimated through history matching. This calibration factor was very quickly estimated and a match was found for both wells A and B (Figure 15) for oil, gas and water rates using a bottom hole pressure. The major drawback of using a finite difference reservoir simulator is the intensive numerical computation required even in today’s new generation parallel reservoir simulators. For example, the case described in this section required a 6-hour run on a good workstation. For unconventional reservoirs, we could trade the accuracy for a faster run time.

Figure 14.

(A) fracture geometry and conductivity resulting from the stimulation of two wells and (B) interwell permeability resulting from the interaction between the hydraulic and natural fractures.

Figure 15.

History matching of oil (green), gas (red), and water (blue) by using the bottom hole pressure (BHP) as a constraint. Notice the good match of both well measurements A and B, achieved very easily and quickly by using one single history matching parameter.

9.3 Unconventional reservoir simulation using fast marching method (FMM)

The motivation and the unique features of the Fast-Marching Method (FMM) simulator needed for unconventional reservoirs were described in Ouenes et al. [6] and Paryani et al. [18]. The input 3D models (Figure 4) and hydraulic fracture geometry (Figure 14) were input in the FMM simulator along with the PVT (Pressure, Volume, Temperature) and other inputs. The resulting pressure depletion at the end of the simulation derived from the FMM simulator (Figure 16B) shows the same features as those seen in the pressure depletion estimated in the classical reservoir simulator (Figure 16A). The same conclusions can be seen when examining the pressure distribution in a cross-section view as shown in Figure 17.

Figure 16.

(A) Areal view of the pressure depletion from finite difference reservoir simulator compared to the (B) pressure derived from the fast marching method simulator.

Figure 17.

(A)–(D) Cross section view of the pressure depletion from a finite difference reservoir simulator compared to the (E)–(H) pressure derived from the fast marching method (FMM) simulator.

There is a major difference between the classical finite difference reservoir simulation run and the one using the FMM: using the same computer, the full-scale heterogeneous model using the compositional finite difference reservoir simulator requires six (6) hours run time due to the large number of components while the FMM multi-phase black oil simulator results were derived in less than 1 minute.

With such a rapid evaluation tool and robust workflow that leverages the multiple constraints derived from the use of surface drilling data, the complex balance between finding the optimal Net Present Value (NPV) per well or per section could be easily estimated in few days or even hours. Using the current industry tools to achieve the same objective will take many weeks if not months and will have a large uncertainty if no well logs or seismic are available as it is very usual the case in unconventional reservoirs where well data and seismic are sacrificed at the altar of cost cutting measures. Fortunately, the surface drilling data provides a reasonable alternative that enables the entire reservoir modeling and management workflow.

10. Conclusions

The use of surface drilling data provides a reasonable engineering solution to the lack of well data in unconventional reservoirs. The correction of the MSE by removing friction losses turns surface drilling data into a major source of information for unconventional well planning. This information includes an estimation of geomechanical logs, pore pressure, stresses, porosity and natural fracture. These rock properties could be used as a first approximation in a well-centric approach to geoengineer completions. Moreover, combining these various logs from different wells into 3D reservoir models provides even more opportunities including using them in reservoir geomechanics, 3D planar hydraulic fracture design and reservoir simulation. The use of MPM modeling tools with the 2D horizontal plane strain framework allows the characterization of the lateral stress gradients, differential stress and the orientation of the maximum stress, all of which are key inputs needed to plan the optimal position and orientation of unconventional wellbores. Using geomechanical logs derived from surface drilling data to identify the geologic interfaces and deploying MPM tools in vertical 2D sections of the reservoirs, allows for capturing the effects of laminations and the loss of hydraulic fracturing energy in weak interfaces. These 2D decoupled approaches provide useful information to constrain a 3D grid based planar hydraulic fracturing approach. These geomechanical constraints define both the lateral and vertical directions which enable the estimation of distribution of proppants with a higher degree of certainty. The resulting realistic fracture geometry and its conductivity can be used in a commercial finite difference reservoir simulator or alternatively in a simplified fast marching method simulator providing similar information as the one given by the finite difference reservoir simulator, but in a fraction of the time. When using all these 3D models and their results in a fast marching method simulator the impact of the interference between wells and other optimization challenges can be estimated quickly while providing similar results as those derived with a finite difference reservoir simulator. By integrating the surface drilling data with 3D reservoir models, hydraulic fracturing design and reservoir simulation into a single software platform, this fast and constrained approach allows for a better management of unconventional wells within a competitive calculation time.

Acknowledgments

The authors would like to thank FracGeo and GO GeoEngineering staff for their contribution to the development, testing and deployment of these technologies. The contribution of Ron Dirksen, John Nairn, Djamel Sia, Bhavina Mistry, Mohcene Laraba, Thibaud Alengry, Nicholas Umholtz, Radouan Smaoui, Srichand Poludasu, Moussa Bari, Ravi Khandelwal, Xiaopeng Li, Bilel Attia, Aishwarya Gogoi, Abassi Ahmed, Dania Almadhun, Sonia K. Sanchez Lohff, Peter O’Conor, and Mark Morford to this large multi-disciplinary effort made it all possible. The authors also thank the Computer Modeling Group for providing the finite difference reservoir simulator.

Nomenclature

$\mathrm{WOB}$

weight on bit,

D

hole diameter

$T$

torque

$N$

rotational speed

ROP

rate of penetration

N

rotational speed of the drill pipe

Kn

mud motor speed to flow ratio

Q

total mud flow rate

g

stress gradient

Po

net pressure which is the difference between the fluid pressure and the minimum or the closure stress

a_1 (t)

short wing in planar 3D fracture

a_2 (t)

long wing in planar 3D fracture

S_Hmax

maximum horizontal stress S_Hmax