Abstract
In this research, an integrated workflow from geomechanics to reservoir simulation is suggested to accurately estimate performances of a shale gas reservoir. Rather than manipulating values of hydraulic fracturing such as fracture geometry and transmissibility, the workflow tries to update model parameters to derive reliable hydraulic fracturing results. A mechanical earth model (MEM) is built from seismic attribute and drilling and diagnostic fracture injection test results. Then, the MEM is calibrated with microseismic measurements obtained in a field. Leakoff coefficient and horizontal stress anisotropy are sensitive parameters of the MEM that influence the propagation of the fracture network and gas productions. Various combinations of calibration parameters from a single-well simulation are evaluated. Then, an appropriate combination is chosen from the whole simulation results of a pad to reduce the uncertainty. Finally, production estimations of the four wells which have slightly different fracture design are compared with seven-year production history. Their results are reasonably matched with actual data having 8% of global error due to successful development of the reservoir model with geomechanical parameters.
1. Introduction
Advancements in horizontal well drilling and multistage hydraulic fracturing have enabled economically viable gas productions from shale formations. Predicting productions of a fractured shale reservoir is essential for field assessments and establishing an optimum development plan. In the last decade, there have been a lot of researches to predict performances of shale gas reservoirs using a numerical simulation. It incorporates geological, geomechanical, and petrophysical properties and examines their effects on gas productions. Therefore, it will be more robust than analytical or empirical methods.
Numerical simulations in shale reservoirs could be categorized into two groups. One method is to use only a reservoir simulation model, and hydraulic fracture (HF) parameters such as fracture height and half-length are considered uncertain parameters. Nejadi et al. [1] presented a novel approach for characterization and history matching of hydrocarbon production from a hydraulic fractured shale gas. They used an assisted history matching algorithm to compare production performances with actual responses and updated HF and discrete fracture network (DFN) model parameters. Their approach could characterize fracture parameters. However, it is limited on the specific fracture design which was done already.
Yu et al. [2] performed a multistage fractured horizontal well numerical simulation on the basis of reservoir properties and completion parameters from an actual horizontal exploration well in China. They analyzed the effect on multistage fractured horizontal well’s performances for matrix permeability, permeability of stimulated reservoir volume (SRV), HF conductivity and half-length, SRV size, and bottomhole pressure. Their approach has some limitations since the method does not consider the relationship between fracture designs and the resulting fracture parameters.
Chang and Zhang [3] proposed a method for characterizing the SRV of shale gas reservoirs using production data. They treated major-fracture properties, the spatial extent of the SRV, and the properties of a dual-porosity, dual-permeability model as uncertain parameters. Using an iterative ensemble smoother, they performed the history matching and suggested the updated values of uncertain parameters. They tried to estimate shale gas productions in a stochastic approach. However, it is based on a synthetic model without considering actual hydraulic fracturing results.
The other method is to integrate a HF simulation and the reservoir simulation model. This method tries to describe HF propagations and their impacts on gas productions. Ajisafe et al. [4] applied a multidisciplinary integrated workflow to a horizontal well to model complex HFs and gas productions. They input the DFN and geomechanical properties into the unconventional fracture model (UFM) and calibrated them with the microseismic data and production history. Their approach was limited to the single-well application.
Cipolla et al. [5] illustrated an application of two complex fracture modelling techniques in conjunction with microseismic mapping to characterize fracture complexity and evaluated completion performances. The authors tried to match the fracture propagation results with measured microseismic mapping based on the various combinations of stress anisotropies and natural fracture spacing. They only focused on the fracture shape (fracture length and complexity) without considering the production performances.
Izadi et al. [6] carried out an integrated subsurface study in a tight gas field to evaluate the stimulation process in horizontal wells. They utilized all available data to build a 3D geomechanical model and tried to quantify heterogeneous in situ stress effects on HF propagation and stimulation efficiency using a 3D fully coupled hydraulic fracturing simulator. However, this method neglected calibration processes such as matching the results with microseismic or production data to reduce the uncertainties of the estimated geomechanical properties.
Lee et al. [7] constructed a numerical reservoir model based on an integrated workflow for modelling an unconventional reservoir. They investigated both fracture closure and proppant placement effects. Then, they concluded that the both parameters need to be considered to estimate gas productions in a shale gas well. They adopted the integrated workflow from a static modelling and carried out fracture simulation and production estimation in the same platform. However, the proposed method is still limited to a single-well analysis rather than incorporating multiple wells in a pad for comprehensive analysis.
Pankaj et al. [8] used an integrated modelling workflow to determine well spacing and well-to-well interference for the Marcellus Shale. Based on a geomechanical model, they simulated HF propagations. After that, they calibrated them with microseismic and well production data by varying parameters such as horizontal stress anisotropy, fracturing fluid leakoff, and the DFN. After that, they did well spacing sensitivity to reveal an optimum distance that the wells need to be spaced to maximize recovery and the number of wells per section. However, they only used interpreted vertical log data to construct the geomechanical model so they did not consider spatial heterogeneity of mechanical properties.
In this study, a novel workflow of geomechanics, hydraulic fracturing, and reservoir simulation is suggested for a shale gas reservoir. The workflow focuses on deriving hydraulic fracturing results which are similar to microseismic data in terms of their geometries and also matched with gas production data. It gives much consistent results than manipulating values of hydraulic fracturing to match measured data. The proposed workflow is verified in an actual pad having four wells of slightly different fracture designs.
A geomechanical model is constructed using drilling results, diagnostic fracture injection test (DFIT) analysis, and seismic attribute data. Rock mechanical properties are calculated from seismic attribute. Reservoir pressure is distributed based on the given borehole information: drilling mud densities and the result of DFIT analysis. After that, closure stress is computed using the above parameters, and it is verified with DFIT results. DFIT analysis also gives information of leakoff coefficients, and these values are used in the calibration of the model. During the calibration, uncertain parameters being difficult to obtain are derived and validated. Finally, gas productions of a whole well pad in the shale reservoir are estimated using the developed workflow and compared with actual production data for validation.
2. Materials and Methods
2.1. Diagnostic Fracture Injection Test
The diagnostic fracture injection test is useful to infer geomechanical properties including fracture closure stress, reservoir pressure, leakoff coefficient, and permeability. A small volume of fluid is injected into the formation, and it creates a hydraulic fracture. After the injection, the pressures in the wellbore are monitored for several days and these data are analyzed using special plots. A typical DFIT sequence is shown in Figure 1.
The analysis of DFIT data is conducted in two ways: before-closure analysis and after-closure analysis. Before-closure analysis focuses on the early pressure falloff period and provides fracture closure stress, instantaneous shut-in pressure (ISIP), and leakoff mechanism and coefficient. When the surface injection is stopped, the friction decreases rapidly and ISIP could be derived which is the sum of closure stress and net pressure. ISIP is estimated by placing a straight line on the pressure falloff plot.
Nolte [10] suggested a -function, which is a dimensionless measure of time to investigate the pressure decline behaviour (as shown in the equations below). where is the pumping time, is the actual time, is the dimensionless shut-in time, is the average decline rate function, is the initial decline rate, and is the dimensionless pressure difference function.
According to the leakoff mechanisms, pressure derivative () and -function derivative () curves show characteristic shapes (Figure 2).
(a) Normal leakoff
(b) Pressure-dependent leakoff
(c) Fracture height recession
(d) Fracture tip extension
In the normal leakoff, the fracture surface area for leaking fluid is constant so the -function derivative curve is a straight line before closure. The fracture closure is identified when the -function derivative starts to deviate downward from the straight line. The pressure-dependent leakoff (PDL) indicates the existence of natural fracture and shows a “hump” trend in the -function derivative. The leakoff rate of the fracture height recession or transverse storage is relatively smaller than the normal leakoff, and it has a “belly” trend below the linearity of the -function derivative. Fracture tip extension means that a fracture continues propagating after the shut-in caused by ultralow permeability of the formation.
After-closure analysis is performed to identify the flow regime and reservoir properties. Firstly, data points of pseudolinear or pseudoradial flow are derived using the slope of the derivative curve on a log-log plot. Then, the falloff pressure data () is plotted versus a linear time function () and radial time function () (as shown in the equations below). where is the fracture closure time.
A Cartesian plot of observed pressure during the falloff versus linear and radial time functions should yield a straight line with an intercept equal to initial reservoir pressure () and with slopes of and (shown below).
Permeability () can be determined from the slope, fracture closure time, and volume injected during the test (shown below): where is the permeability in md, is the net pay thickness in m, is the fluid viscosity in cp, is the injected fluid volume in m3, and the slope is in MPa.
2.2. Unconventional Fracture Model
An unconventional fracture model has been developed [12], which describes fracture propagation, rock deformation, fluid flow, and proppant transport simultaneously during hydraulic fracturing treatments. Compared to conventional planar fracture models, UFM could simulate the interaction of HFs with preexisting natural fractures and consider the interaction among HF branches by computing the “stress shadow” effect on each fracture exerted by the adjacent fractures.
Figure 3 shows a basic workflow for building a UFM. At each time step, flow and elasticity equations are solved to derive new pressure and flow distributions in the fracture networks. Injected fluids are modelled as a power law model, and fluid flow along fracture branches is described according to laminar or turbulent flow.
Fracture width is calculated using fluid pressure and elastic properties of the rock such as Young’s modulus and Poisson’s ratio through the elasticity equation. More detailed description of the governing equations can be found in Weng et al. [12]. Using UFM, the fluid flow equation is solved for each component of the fluids and proppant pumped, so proppant placement in the fracture networks can be described after hydraulic fracturing treatment as shown in Figure 4.
3. Results and Discussion
3.1. DFIT Analysis
In the field of interest, there are eight offset wells where DFIT was done at the toe of the horizontal wells. Figure 5 shows the recorded pressure data of one well, and it could be seen that there are two pressure peaks in the plot. Because a pressure-activated rupture disc valve is installed in the toe stage, pressure data after the 2nd peak should be used to analyze.
Figure 6 shows an example of DFIT analysis. At first, ISIP is derived. Then, fracture closure stress and leakoff mechanisms are estimated using a -function plot. Barree et al. [13] suggest that the transverse storage leakoff mechanism could mask the PDL leakoff characteristics so conducting a leakoff analysis is important. Although it looks like the “belly” trend in the -function plot, its slope has a positive slope which means the PDL leakoff is dominant (Table 1). From the log-log plot, flow regimes are identified and finally reservoir properties are derived by plotting the falloff pressure data versus time functions.
(a) ISIP plot
(b) -function plot
(c) Log-log plot
(d) Radial plot (after closure)
Table 1 shows the results of DFIT data analysis, and these data are used to construct a mechanical earth model. A MEM is defined as a model which describes rock elastic and strength properties, in situ stresses, and reservoir pressure based on a stratigraphic column [14].
A calculated PDL coefficient of leakoff is utilized to derive a leakoff multiplier using Equation (9), which is used as an input of the MEM. where is the altered leakoff coefficient, is the original matrix leakoff coefficient, is the leakoff multiplier, is the PDL coefficient of leakoff, and is the net pressure.
3.2. Mechanical Earth Modelling
Overburden stress () is integrated from the rock density () along the depth. Young’s modulus () and Poisson’s ratio () are derived from seismic attribute data. These seismic attribute data are the results of seismic inversion analysis. During the inversion analysis, well log data are used to calibrate the inversion results with measured data. where and are the -wave and the -wave velocities, respectively.
Reservoir pressure () is distributed based on the given borehole information: drilling mud densities and the result of DFIT analysis. Finally, closure stress () equivalent to minimum horizontal stress is computed from Equation (11) using the above parameters. where is the Biot constant and is the tectonic stress.
In Figure 7(a), calculated closure stress (color bar) was compared with the results of DFIT analysis (black dot). DFIT was performed only in the toe stage of the horizontal wells, so there is one measured point per well. Computed stress is well matched with the DFIT results among the offset wells, and its distribution is shown in Figure 7(b). In addition, the direction of the maximum horizontal stress was estimated as N45E from a world stress map and microseismic results.
(a) Calculated stress (color bar) vs. DFIT results (black dot)
(b) Closure stress distribution
However, there are still some uncertain parameters of the MEM such as leakoff coefficient, horizontal stress anisotropy, and natural fracture property. By comparing HF and reservoir simulation results with field measured data, these uncertain parameters are calibrated and the updated MEM is obtained.
3.3. Mechanical Earth Model Calibration
The calibration process was conducted for a single well, and then, it was expanded to a whole pad scale including the four wells. Firstly, four leakoff multiplier cases were constructed from the DFIT analysis. Using Equation (9), leakoff multipliers of offset wells were estimated as shown in Table 1. During the calculation, net pressure was assumed as 4.68 MPa which is the average value from hydraulic fracturing simulation. Case 1 assumed that there was not any PDL leakoff effect, and other cases were constructed from the range of leakoff multipliers. Hydraulic fracturing simulations were performed and compared with recorded microseismic data in the field. Microseismic data will give a rough implication of the generated fracture geometry. In case 3 of Figure 8, simulated fracture geometry is well matched with microseismic data and it is also identified in the history matching result in Figure 9.
(a) Leakoff case 1 ()
(b) Leakoff case 2 ()
(c) Leakoff case 3 ()
(d) Leakoff case 4 ()
Horizontal stress anisotropy means the difference between the maximum and the minimum horizontal stresses. Since measurement of the maximum horizontal stress is not directly possible, three different anisotropies were assumed in this research: the same with the minimum horizontal stress (case 1), 5% higher (case 2), and 10% higher (case 3). These anisotropy ranges are estimated from a correlation between a breakdown pressure and the maximum horizontal stress using DFIT data. The anisotropy affects both fracture geometry and gas production rates. Figure 10 shows HF simulation results of the three cases, and the shapes of case 2 and case 3 are similar to the microseismic data which aligned along the maximum horizontal stress direction. However, case 2 shows much lower gas production rates compared to the actual production data (Figure 11).
(a) Stress anisotropy case 1 ()
(b) Stress anisotropy case 2 ()
(c) Stress anisotropy case 3 (.1)
One of the most important factors controlling the development of a HF network is preexisting natural fracture (NF) networks, represented by DFN [5]. Sometimes, seismic attribute analysis is used to extract information of NF networks such as orientation, intensity (spacing), and stiffness [15]. However, there is not enough information to construct a DFN model in this field. Therefore, a statistical approach is used to investigate the sensitivities of the created fracture geometry to various DFN. Four models are constructed having different NF spacing: 75 m (case 1), 100 m (case 2), 150 m (case 3), and 200 m (case 4). During the hydraulic fracturing process, injected fluid could leak off into the NFs so it leads to creating a complex fracture network having shorter half-length.
Figure 12 shows HF simulation results of the four NF spacing cases. Blue patches in the background represent the DFN model. Using a focal mechanism analysis of microseismic data, the orientations of the NF network were determined as N40E and N71W. From the hydraulic fracturing simulation, case 2 and case 3 are similar to microseismic data obtained. Although it has influence on the fracture geometry, it does not have much effect on the gas production rates as seen in Figure 13. The narrower NF spacing is, the much interaction of the HF with the NFs exists. However, it interrupts to develop the main fracture longer. This affects the SRV and the resulting gas productions.
(a) NF spacing case 1 (75 m)
(b) NF spacing case 2 (100 m)
(c) NF spacing case 3 (150 m)
(d) NF spacing case 4 (200 m)
3.4. Production Estimation
The HF simulation of a pad was conducted using the calibrated MEM parameters of the single well. For NF spacing, both cases 2 and 3 were tested, and case 3 was appropriate for the whole pad. Figure 14 shows the fracture simulation results of the pad. It is found that overall created fracture geometry is similar to that of microseismic data, so the MEM is updated using the calibration factors.
(a) Top view
(b) Side view
Finally, production estimation of the pad is performed to validate the integrated workflow in the shale gas reservoir. An unstructured grid for the fracture network is generated consisting of 15 layers with 3 million cells. In addition to the reservoir grid, a fluid model, relative permeability, and pressure-dependent permeability are put into the reservoir simulator. Figure 15 represents pressure-dependent permeability multipliers of each region. As the reservoir pressure is depleting, permeability of each region is decreasing. Proppant placement helps to support against fracture closure, so permeability is maintained during the depletion.
About seven years, gas production data of the four wells are used for comparison and tubing head pressures of each well are utilized as a constraint. Figures 16 and 17 present simulated results of tubing head pressures and gas productions, respectively. The black dots indicate field measured data, and solid lines represent simulation results. Although fracture designs of four wells are slightly different, the model could estimate seven-year production history of the pad reasonably due to proper use of parameters calibrated. The global error shows 8% for the whole pad.
(a) Wells 1 and 2
(b) Wells 3 and 4
(a) Wells 1 and 2
(b) Wells 3 and 4
Therefore, the integrated approach proposed in this study will be useful to predict reliable future gas productions and to optimize HF design or operation conditions in the field. This study was conducted based on a deterministic single mechanical earth model, but stochastic approaches could also be taken to evaluate the uncertainties of rock mechanical, natural fracture, and stress-related parameters.
4. Conclusions
In this research, an integrated workflow is suggested and validated with the actual field data. It consists of building a mechanical earth model, calibrating the model, and estimating gas productions. The MEM is constructed using seismic attribute and drilling and DFIT analysis results. Based on the model, hydraulic fracture propagation of a single well in a pad is simulated and compared with microseismic data. In addition, gas production rates are compared for each of uncertain parameters for selecting proper values.
From the DFIT analyses, pressure-dependent leakoff mechanism is dominant in this field and leakoff multipliers are calculated. High leakoff coefficient causes most of injected fluid coming into the matrix and natural fracture. Therefore, it leads to generating shorter fracture half-length and less gas productions. The horizontal stress anisotropy affects the shape of the created fracture network. As the anisotropy increases, a planar-type fracture network is generated and the opposite case results in a complex-type fracture network. The natural fracture spacing is related to the interaction of the HF with NFs. However, it does not have much impact on the gas production rates.
The MEM is updated using calibrated parameters from a single-well case, and the pad scale simulation is conducted to validate the model parameters selected. Simulation results of the four wells, which have slightly different fracture designs, are well matched with microseismic data. After the hydraulic fractures are modelled and calibrated, an unstructured grid that contains fracture geometry and proppant distribution is built. Production estimations of the pad are compared with the seven-year production history, and their results are reasonably similar to the actual data (global error: 8%) due to proper use of geomechanical parameters calibrated.
The developed workflow has some availability for further studies. It will be useful to predict future gas productions and applicable to other shale gas reservoirs. It could be also used to decide an optimal HF design using the calibrated MEM.
Nomenclature
| : | PDL coefficient of leakoff (1/kPa |
| : | Original matrix leakoff coefficient (m/√min) |
| : | Altered leakoff coefficient (m/√min) |
| : | Young’s modulus (1/GPa) |
| : | Linear flow time function (dimensionless) |
| : | Radial flow time function (dimensionless) |
| : | Average decline rate function (dimensionless) |
| : | Initial decline rate (dimensionless) |
| : | -function (dimensionless) |
| : | Net pay thickness (m) |
| : | Permeability (md) |
| : | Slope of data on pseudolinear flow graph (MPa) |
| : | Slope of data on pseudoradial flow graph (MPa) |
| : | Initial reservoir pressure (MPa) |
| : | Reservoir pressure (MPa) |
| : | Falloff pressure (MPa) |
| : | Net pressure (MPa) |
| : | Actual time (minutes) |
| : | Closure time (minutes) |
| : | Pumping time (minutes) |
| : | Dimensionless shut-in time (dimensionless) |
| : | -wave velocity (km/s) |
| : | -wave velocity (km/s) |
| : | Biot constant (dimensionless) |
| : | Viscosity (cp) |
| : | Poisson’s ratio (dimensionless) |
| : | Density (kg/m3) |
| : | Closure stress (MPa) |
| : | Minimum horizontal stress (MPa) |
| : | Maximum horizontal stress (MPa) |
| : | Tectonic stress (MPa) |
| : | Vertical stress (MPa). |
Data Availability
The data used to support the finding of this study are included within the article.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The authors would like to thank Korea Gas Corporation and KOGAS Canada Limited for the permission to publish this article. The corresponding author also likes to thank MOTIE projects (20204010600250, 20162010201980) and the Institute of Engineering Research at Seoul National University, Korea.