Next Article in Journal
Thermo-Mechanical Stress Comparison of a GaN and SiC MOSFET for Photovoltaic Applications
Next Article in Special Issue
Fractional Time Derivative Seismic Wave Equation Modeling for Natural Gas Hydrate
Previous Article in Journal
Link between Energy Efficiency and Sustainable Economic and Financial Development in OECD Countries
Previous Article in Special Issue
Adaptive Processing for EM Telemetry Signal Recovery: Field Data from Sichuan Province
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Analytical Model for Production Analysis of Hydraulically Fractured Shale Gas Reservoirs Considering Irregular Stimulated Regions

1
Department of Energy and Resources Engineering, College of Engineering, Peking University, Beijing 100871, China
2
Key Laboratory of Theory and Technology of Petroleum Exploration and Development in Hubei Province, Wuhan 430074, China
3
Key Laboratory of Tectonics and Petroleum Resources, Ministry of Education, Wuhan 430074, China
4
School of Earth Resources, China University of Geosciences, Wuhan 430074, China
*
Author to whom correspondence should be addressed.
Energies 2020, 13(22), 5899; https://doi.org/10.3390/en13225899
Submission received: 23 September 2020 / Revised: 5 November 2020 / Accepted: 6 November 2020 / Published: 12 November 2020
(This article belongs to the Special Issue Development of Unconventional Reservoirs 2020)

Abstract

:
Shale gas reservoirs are typically developed by multistage, propped hydraulic fractures. The induced fractures have a complex geometry and can be represented by a high permeability region near each fracture, also called stimulated region. In this paper, a new integrative analytical solution coupled with gas adsorption, non-Darcy flow effect is derived for shale gas reservoirs. The modified pseudo-pressure and pseudo-time are defined to linearize the nonlinear partial differential equations (PDEs) and thus the governing PDEs are transformed into ordinary differential equations (ODEs) by integration, instead of the Laplace transform. The rate vs. pseudo-time solution in real-time space can be obtained, instead of using the numerical inversion for Laplace transform. The analytical model is validated by comparison with the numerical model. According to the fitting results, the calculation accuracy of analytic solution is almost 99%. Besides the computational convenience, another advantage of the model is that it has been validated to be feasible to estimate the pore volume of hydraulic region, stimulated region, and matrix region, and even the shape of regions is irregular and asymmetrical for multifractured horizontal wells. The relative error between calculated volume and given volume is less than 10%, which meets the engineering requirements. The model is finally applied to field production data for history matching and forecasting.

Graphical Abstract

1. Introduction

Unconventional hydrocarbon resources such as tight and shale oil/gas store in tight formations with ultra-low permeability. With the development of hydraulic fracturing technologies, multifractured horizontal wells (MFHWs) have rapidly emerged as the primary means for exploiting this type of resource. Meanwhile, some technologies are utilized, such as foam injection and carbon dioxide injection [1,2] on recovery enhancement, and photo-Fenton and floatation on sustainable management of flow-back water after hydraulic fracturing [3]. In unconventional reservoirs, due to the propagation of fractures in different directions, branch fractures will be created around the main hydraulic fractures, which have a significant impact on the pressure and rate transient analysis for the fluid flow in the reservoirs.
In order to analyze production data and make long-term forecasts, analytical and numerical tools have been developed. Among them, a large number of numerical approaches [4,5,6,7,8], such as finite element method and boundary element method, are adopted to study the multiple flow mechanisms and fracture characteristics in shale gas reservoirs. Compared with the numerical models, the analytical models are simpler to set up and require fewer original data. Therefore, they are widely used in practice for field applications. Extensive studies have focused on analyzing production data of transient flow period for such complex systems [9,10,11,12,13,14]. Lee et al. [15] derived a tri-linear model to represent the transient flow behavior within a single fracture in an infinite homogeneous reservoir. Wattenbarger et al. [16] proposed an analytical solution accounting for transient linear flow contributions from each region for a fracture with infinite conductivity in a closed reservoir. El-Banbi [17] and Bello et al. [18] extended previous analytical models to develop a series of linear dual-porosity solutions for fluid flow in tight reservoirs, where constant-rate and constant-pressure inner boundaries were considered. Ozkan et al. [19] and Brown et al. [20] extended the concept of tri-linear flow and presented a solution to simulate the pressure-transient behavior in three flow regions, including the hydraulic fracture, the inner region between hydraulic fractures, and the outer region beyond the tip of hydraulic fractures. Stalgorova et al. [21] and Heidari et al. [22] introduced an enhanced-fracture-region (EFR) model to consider the cases where branch fractures cover the interfracture region. Mahdi et al. [23] firstly provided the transient shape factor among different regions in the triple-porosity model.
The models represent an excellent attempt to capture reservoir parameters and flow characteristics for the complex system associated with the MFHWs. However, most analytical solutions were derived by Laplace transformation [17,18,19,20,21,22,23] and Boltzmann’s transformation [24]. Numerical inversion algorithm [25] is needed to transform the Laplace space into real-time space. In order to avoid the inconvenience of Laplace transformation, several approximated analytical solutions have been presented. Shahamat et al. [26] introduced the concept of continuous succession of pseudo-steady states based on the capacitance-resistance methodology to analyze the production data in tight and shale reservoirs. In their work, the gas adsorption was ignored, which would underestimate the shale gas production. Ogunyomi [27,28] introduced a model with two compartments in which the first compartment represents the volume of the enhanced fracture region and the second one is the volume of the matrix. The approximate solution was obtained by converting governing partial-differential equations into a system of ordinary differential equations through integration. Then the system of ordinary differential equations was solved by calculating the eigenvalues and eigenvectors. Qiu et al. [29,30] extended this method to a triple-porosity model for production analysis. Nevertheless, they both did not take the non-Darcy flow into account and just derived the approximate solutions based on the oil phase rather than gas phase with the pressure dependent parameters. In this work, we extended the method to consider the gas flow and even account for non-Darcy flow in the hydraulic fractures, and gas adsorption and slippage in the matrix.
In this work, our goal is to derive a new integrative analytical solution for shale gas reservoirs bypassing the Laplace transformation, which can effectively account for non-Darcy flow in the hydraulic fractures, and gas adsorption and slippage in the matrix. The analytical model is created on the basis of a three-region model, which consists of the hydraulic fracture region, the inner region also called stimulated region connected to the hydraulic fracture, and the outer region representing the low permeability matrix. Then we introduce modified pseudo-variables [31,32], i.e., pseudo-pressure and pseudo-time, to eliminate the nonlinearity of governing equations. The whole derivation is based on the integration in real-time space bypassing Laplace transform and numerical inversion. The derived analytical solution is validated with numerical simulations. Further, several numerical cases are designed to verify the applicability in an irregular stimulated region using the new model. The proposed method is finally applied for shale gas production analysis and forecasting in the field.

2. Method

Production of shale gas requires horizontal wells with multistage hydraulic fracturing. The propagation of fractures is uncertain and branch fractures will be created around the main hydraulic fractures, resulting in a stimulated region around each main hydraulic fracture. In order to model such a complex system, the reservoir is simplified as a three-region system, where the first region is the hydraulic fracture region which is regarded as the sole connection to the well. The second one represents the stimulated region with aggregated volume of all the microfractures and the third one is the adjacent ultra-low-permeability matrix directly connected to the stimulated region. Figure 1a shows the schematic 3D model. Three regions are contained in the model: region 1 is the hydraulic fracture, region 2 (darker color) with higher permeability around each hydraulic fracture, and region 3 (lighter color) with lower permeability connected to the region 2. The arrows represent the flow directions. For this model, the flow directions are parallel, and the system is symmetrical with respect to the hydraulic fracture and horizontal-well. Thus, it is feasible to use one quarter of the reservoir shown in Figure 1b to replace the whole reservoir in order to simplify the derivation process. According to Anderson et al. [33], they verified that when the permeability of region 2 is less than or equal to 500 nD, the contribution of the region beyond fractures can be neglected for the 20-years production. For the case that the distance from the fracture face to permeability boundary ( x 1 ) is less than the half distance between fractures ( x 2 ), the contribution would be smaller. Meanwhile, Stalgorova et al. [21] also set numerical models to illustrate the negligible contribution of region beyond fractures, and they obtained that the difference of 20-years production is negligible after comparing the results of numerical simulation with and without region beyond fractures. In the work of Heidari et al. [22], they also did not take the region beyond fractures into account. Therefore, the contribution of the region beyond the tips of the hydraulic fractures is assumed to be negligible in this work.
In this work, our analytical model is derived under the following assumptions:
  • The reservoir is homogeneous, isopachous and isothermal in each region.
  • Flow process is 1-D linear in each region.
  • Flow is single gas phase.
  • The high-velocity non-Darcy flow in the hydraulic fracture is considered.
  • The bottom-hole pressure is constant.
  • The impact of gravity is neglected.
  • Gas desorption meets the Langmuir isotherm adsorption equation.

2.1. Gas Adsorption/Desorption Effect

In contrast to conventional gas reservoirs, gas adsorption is an important feature of shale gas reservoir. The Langmuir isotherm adsorption equation [34] is widely used to calculate the shale gas adsorption and its expression is as follows:
V = V L P P L + P
where V is the volume of the adsorbed gas, and P is pressure. V L and P L stand for Langmuir volume and pressure, respectively. When considering the gas adsorption, the effect of adsorption on compressibility of the reservoir is essential. According to Bumb et al. [35], the new compressibility equation can be expressed as,
C t = C f + C w S w + C g ( 1 S w ) + C g d
C g d = 0.031214 ρ m V L B g ¯ P L ϕ ( p ¯ + P L ) 2
where C f is rock compressibility, C w is water compressibility, C g is free gas compressibility and C g d is adsorbed gas compressibility, S w is water saturation, ρ m is matrix density, B g ¯ is gas reservoir volume factor, and ϕ is the porosity.
Another important change is that the compressibility factor is modified by King [36],
z * = z 1 S w + 0.031214 ρ m z T ϕ p s c T s c V L 1 p + p L
where P s c is the standard condition pressure, T s c is the standard condition temperature, z is compressibility factor, T is the reservoir temperature.

2.2. High-Velocity Non-Darcy Flow Effect

For gas flow in the hydraulic fracture, high-velocity non-Darcy effect is considered in this study. Forchheimer [37] proposed that Darcy’s law is inadequate to describe high-velocity gas flow without adding an inertial effect, which is proportional to the square of the flow velocity. To account for the non-Darcy flow effect, an inertial term must be included. The Forchheimer’s flow equation is given as,
p = μ k v + β ρ v v
In order to reduce the nonlinearity, the equivalent permeability is introduced to obtain the extended Darcy’s law [37],
p = μ k e q v
Substituting Equation (6) into Equation (5), the equivalent permeability of hydraulic fracture yields,
k e q F = μ k i F μ + β ρ v k i F
where
β = 4.1 × 10 11 ( k i F ) 1.5
among which, β is the non-Darcy flow coefficient, k i F is the hydraulic fracture permeability, k e q F is the equivalent permeability of hydraulic fracture, v is the flow velocity, and ρ is the gas density.

2.3. Non-Darcy Flow in the Matrix

For the nano-porous media in shale reservoir, Darcy’s law has difficulty describing the actual gas behavior. The gas flow can be classified as four flow conditions, such as continuum flow, slip flow, transition flow, and free-molecular flow. According to the previous publications [38], the four flow conditions can be qualified by the Knudsen number. However, the Knudsen number varies from 10−3 to 1 in most shale reservoirs. In order to represent the non-Darcy gas flow in matrix, the apparent permeability is presented by the following general form:
k a m = k m f K n
where k m is the intrinsic permeability of the porous medium, f ( K n ) is the correlation term expressed as a function of the Knudsen number, which is modeled as [38],
f K n = 1 + α K n 1 + 4 K n 1 + K n
in which,
α = 128 15 π 2 tan 1 4 K n 0.4
Meanwhile, for a capillary tube of radius, r, the intrinsic permeability can be derived [39],
k m = r 2 8
The Knudsen number K n is defined as the ratio of the molecular mean free path length and the pore radius in shale matrix,
K n = λ r
The mean free path can also be calculated,
λ = μ g P π R T 2 M
Substituting Equation (14) into Equation (13) on basis of the real gas condition, we can obtain:
K n = μ g z r P π R T 2 M
where μ g is the gas viscosity, z is compressibility factor, r is the pore radius, P is the reservoir pressure, R is the universal gas constant, T is the reservoir temperature, and M is the gas molecular weight.

2.4. Derivation of Linearized Gas Diffusivity Equation.

For the flow of shale gas, the gas diffusivity equation will be nonlinear which makes deriving the analytical solution difficult. On one hand, with the reduction in average reservoir pressure, the gas properties like the gas viscosity ( μ ), gas compressibility ( C t ), and gas compressibility factor ( z ) will change with the pressure. On the other hand, when incorporating the significant mechanisms in shale gas reservoir like gas adsorption and non-Darcy flow effect, the permeability is a variate with pressure rather than a constant.
In order to deal with this problem, the pseudo-pressure and pseudo-time instead of the pressure and time are adopted to linearize the equations. We set a general real gas diffusivity equation in three-dimensional Cartesian coordinate system as example. When the non-Darcy effect is coupled, the k ( p ) can be calculated by Equation (7) or Equation (9). Certainly, the c t p can be replaced by c t * p to represent the gas adsorption effect.
x k p p μ p z p p x + y k p p μ p z p p y + z k p p μ p z p p z = ϕ c t p p z p p t
where k ( p ) is the pressure-dependent permeability, μ ( p ) is the pressure-dependent viscosity, c t p is the pressure-dependent compressibility.
Considering the pressure dependence of the permeability, we define a general modified pseudo-pressure transformation. Surely, k ( p ) can be calculated by Equation (7) or Equation (9) to couple with the non-Darcy effect. z ( p ) can be replaced with z * ( p ) calculated by Equation (4) to consider the gas adsorption.
m p = 1 k i p b p k p p μ ( p ) z ( p ) d p
where m ( p ) is the pseudo-pressure, k i is the intrinsic permeability, z ( p ) is the pressure-dependent compressibility factor.
Substituting Equation (17) into Equation (16), the right side of the partial differential equation is still nonlinear.
2 m p x 2 + 2 m p y 2 + 2 m p z 2 = ϕ μ p c t p k p m p t
To linearize Equation (18), a general modified pseudo-time transformation is defined [24].
t a = 1 k i 0 t k p ( μ c t ) i μ ( p ) c t ( p ) d t
where t a is the pseudo-time, μ i is the initial viscosity, c t i is the initial compressibility.
After substituting the Equation (19) into Equation (18), the general linear partial differential equation is derived as,
2 m p x 2 + 2 m p y 2 + 2 m p z 2 = ϕ μ i c t i k i m p t a
Considering the permeability and compressibility in three regions are different, therefore, the pseudo-time in three regions can be shown as follows.
t a F = 1 k i F 0 t k e q F p ( μ c t ) i F μ F ( p ) c t F ( p ) d t
t a 1 = 1 k i 1 0 t k 1 p ( μ c t ) i 1 μ 1 ( p ) c t 1 ( p ) d t
t a 2 = 1 k i 2 0 t k a m p ( μ c t * ) i 2 μ 2 ( p ) c t 2 * ( p ) d t
where t a F , t a F , t a 2 are the pseudo-time in the hydraulic fracture region, region 1, and region 2 respectively.
Therefore, we can obtain the linearized gas diffusivity equation in three regions respectively. Firstly, we consider the gas adsorption and gas slippage in matrix and the high-velocity non-Darcy flow in hydraulic fracture region. Thus, we use a new total compressibility coupling the gas adsorption effect in diffusivity equation of region 2. Finally, through the linearization by the modified pseudo-pressure and pseudo-time, the non-Darcy flow effect for matrix and fracture are included in the governing linear equations in different regions.

2.5. Model Description in Different Regions.

2.5.1. Model Description in Matrix Region (Region 2)

The system of equations based on the conceptual model is presented as follows. For the low-permeability matrix region (region 2), the diffusivity equation for gas flow is derived as,
2 m 2 ( p ) x 2 + 2 m 2 ( p ) y 2 + 2 m 2 ( p ) z 2 = ( ϕ μ i c t i * ) 2 k i 2 m 2 ( p ) t a 2
where m 2 ( p ) is the pseudo-pressure in region 2, x, y, z are the Cartesian coordinates, t a 2 is the pseudo-time for gas flow, c t i 2 * and μ i 2 are the initial modified total compressibility and viscosity in region 2, k i 2 and ϕ i 2 are the permeability and porosity in region 2, respectively.
The initial condition for the region is that the pseudo-pressure of the region is equal to initial pseudo-pressure at t = 0 .
m 2 ( p ) ( x , y , z , 0 ) = m ( p i )
where m ( p i ) is the initial pseudo-pressure.
The boundary conditions are defined as no-flow at top and bottom of the reservoir. Additionally, the both ends of y-direction can also be regarded as no-flow boundaries.
m 2 ( p ) y y = 0 , x f = 0
m 2 ( p ) z z = z 0 , z e = 0
Due to the plane of symmetry between adjacent fractures, the location x = x 2 is also a no-flow boundary.
m 2 ( p ) x x = x 2 = 0
The continuity of flux and pressure across the boundaries between the regions is assumed. In Ogunyomi’s [27] and Qiu’s [29] work, only one is chosen as a boundary condition. Considering the average pseudo-pressure will be adopted, the continuity of flux is chosen as the last boundary condition which can be given by,
k i 2 m 2 ( p ) x x = x 1 = k i 1 m 1 ( p ) x x = x 1

2.5.2. Model Description in Stimulated Region Volume (Region 1)

For the stimulated region volume (region 1), the diffusivity equation for gas flow is also expressed as,
2 m 1 ( p ) x 2 + 2 m 1 ( p ) y 2 + 2 m 1 ( p ) z 2 = ( ϕ μ i c t i ) 1 k i 1 m 1 ( p ) t a 1
where m 1 ( p ) is the pseudo-pressure in region 1, x, y, z are the Cartesian coordinates, t a 1 is the pseudo-time for gas flow, c t i 1 and μ i 1 are the initial total compressibility and viscosity in region 1, k i 1 and ϕ 1 are the permeability and porosity in region 1, respectively.
For the whole model, the initial condition remains the same.
m 1 ( p ) ( x , y , z , 0 ) = m ( p i )
The flow in region 1 is also the 1D linear in x-direction, therefore, the outer boundary conditions are identical with those of region 2.
m 1 ( p ) y y = 0 , x f = 0
m 1 ( p ) z z = z 0 , z e = 0
According to the continuity hypothesis, the boundary condition is expressed as,
k i 1 m 1 ( p ) x x = x 1 = k i 2 m 2 ( p ) x x = x 1
k i 1 m 1 ( p ) x x = w 2 = k i F m F ( p ) x x = w 2

2.5.3. Model Description in Hydraulic Fracture Region

For the hydraulic fracture region, the diffusivity equation for gas flow is also expressed as,
2 m F ( p ) x 2 + 2 m F ( p ) y 2 + 2 m F ( p ) z 2 = ( ϕ μ i c t i ) F k i F m F ( p ) t a F
where m F ( p ) is the pseudo-pressure in hydraulic fracture region, x, y, z are the Cartesian coordinates, t a F is the pseudo-time for gas flow, c t i F and μ i F are the initial total compressibility and viscosity in hydraulic fracture region, k i F and ϕ F are the permeability and porosity in hydraulic fracture region, respectively.
For the whole model, the initial condition remains the same.
m F ( p ) ( x , y , z , 0 ) = m ( p i )
With the assumption of the constant bottom-hole pressure, at the location of y = 0 , the pressure is the same as the bottom-hole pressure at any time.
m F ( p ) ( x , y = 0 , z , t a ) = m ( p w f )
The flow in hydraulic fracture region is also the 1D linear in the y-direction, therefore, the outer boundary conditions is expressed as,
m F ( p ) x x = 0 = 0
m F ( p ) y y = x f = 0
m F ( p ) z z = z 0 , z e = 0
According to the continuity hypothesis, the boundary condition is expressed as,
k i F m F ( p ) x x = w 2 = k i 1 m 1 ( p ) x x = w 2

3. Results

3.1. Derivation of Analytical Solution.

The diffusivity equations in the mathematical model are all PDEs. In order to obtain the analytical solution, the first step is to transform the sets of equations into ODEs. In this work, we adopted the integral method other than the common Laplace Transform. By integrating over the spatial domain, the spatial dependence is eliminated, which is even feasible to irregular regions as will be demonstrated later. The pressure in different regions is assumed as average value in our model. Therefore, the pseudo-time is supposed to be independent with space. Integrating the equations with respect to spatial coordinates,
x 1 x 2 0 x f z o z e x ( m 2 ( p ) x ) d x d y d z + x 1 x 2 0 x f z o z e y ( m 2 ( p ) y ) d x d y d z + x 1 x 2 0 x f z o z e z ( m 2 ( p ) z ) d x d y d z = ( ϕ μ i c t i * ) 2 k i 2 t a 2 x 1 x 2 0 x f z o z e m 2 ( p ) d x d y d z
In Equation (43), we move the pseudo-time outside the spatial integral since the pseudo-time is independent of the spatial coordinates. In order to get a simplified equation, the average pseudo-pressure and the effective pore volume is defined as,
m ¯ ( p ) = m ( p ) d x d y d z = m ( p ) V b
V p = ϕ V b
where V b is the volume of the region and V p is the pore volume of the region.
Equation (43) can be rewritten as,
0 x f z 0 z e m 2 ( p ) x x 2 m 2 ( p ) x x 1 d y d z + x 1 x 2 z 0 z e m 2 ( p ) y x f m 2 ( p ) y 0 d x d z + x 1 x 2 0 x f m 2 ( p ) z z e m 2 ( p ) y z 0 d x d y = ( ϕ μ i c t i * ) 2 V b 2 k i 2 d m ¯ 2 ( p ) d t a 2
Using the initial condition and boundary conditions, Equation (46) is simplified to
0 x f z 0 z e k i 2 m 2 ( p ) y x 1 d y d z = ( V p μ i c t i * ) 2 d m ¯ 2 ( p ) d t a 2
According to the equivalent Darcy’s law,
q 2 = T s c Z s c 2 T P s c 0 x f z 0 z e k i 2 m 2 ( p ) y x 1 d y d z
Define
C 2 = T s c Z s c 2 T P s c ( μ i c t i * ) 2
Therefore, Equation (48) can be rewritten as,
C 2 V p 2 d m 2 ¯ ( p ) d t a 2 = q 2
where V p 2 is the pore volume of region 2, q 2 is the flow rate in region 2.
For region 1, we also use the integration method to deal with the equation,
0 x f z 0 z e m 1 ( p ) x x 1 m 1 ( p ) x w 2 d y d z + w / 2 x 1 z 0 z e m 1 ( p ) y x f m 1 ( p ) y 0 d x d z + w / 2 x 1 0 x f m 1 ( p ) z z e m 1 ( p ) y z 0 d x d y = ( ϕ μ i c t i ) 1 V b 1 k i 1 d m ¯ 1 ( p ) d t a 1
After applying the initial condition and boundary conditions, Equation (51) can be rewritten as,
0 x f z 0 z e k i 1 m 1 ( p ) x x 1 k i 1 m 1 ( p ) x x 0 d y d z = ( V p μ i c t i ) 1 d m ¯ 1 ( p ) d t a 1
Note that,
q 1 = T s c Z s c 2 T P s c 0 x f z 0 z e k i 1 m 1 ( p ) x x 0 d y d z
Defining
C 1 = T s c Z s c 2 T P s c ( μ i c t i ) 1
Substituting Equations (30), (45), (49), and (50) into Equation (48) results in,
C 1 V p 1 d m ¯ 1 ( p ) d t a 1 = q 1 + q 2
where V p 1 is the pore volume of region 1, q 1 is the flow rate in region 1.
For hydraulic fracture region, we also use the integration method to deal with the diffusion equation,
0 x f z 0 z e m F ( p ) x w 2 m F ( p ) x 0 d y d z + 0 w / 2 z 0 z e m F ( p ) y x f m F ( p ) y 0 d x d z + 0 w / 2 0 x f m F ( p ) z z e m F ( p ) y z 0 d x d y = ( ϕ μ i c t i ) F V b F k i F d m ¯ F ( p ) d t a F
After applying the initial condition and boundary conditions, Equation (56) can be rewritten as,
0 x f z 0 z e k i F m F ( p ) x w 2 d y d z + 0 w / 2 z 0 z e k i F m F ( p ) y 0 d x d z = ( ϕ μ i c t i ) F V b F d m ¯ F ( p ) d t a F
The flow rate in hydraulic fracture can be expressed as,
q F = T s c Z s c 2 T P s c 0 x f z 0 z e k i F m F ( p ) x 0 d y d z
Defining
C F = T s c Z s c 2 T P s c ( μ i c t i ) F
Substituting Equations (42), (48), (58), and (59) into Equation (57) results in,
C F V p F d m ¯ F ( p ) d t a F = q F + q 1
The next step is to substitute the average pseudo-pressure with the relationship between pressure and the flow rate. Since it is assumed that gas flows sequentially from region 2 to region 1 then to hydraulic fracture, a general analytical solution for one-dimensional linear gas flow is derived to solve the problem (details are shown in Appendix A), which is given by:
m ¯ ( p ) = m ( p w f ) + 4 π 2 [ m ( p i ) m ( p w f ) ] n = 1 q D n ( 2 n 1 ) 2
where m ¯ ( p ) is the average pseudo-pressure, q D n is the dimensionless production from the nth mode.
Combining the assumptions, the average pseudo-pressure in each region can be expressed as,
m ¯ F ( p ) = m ( p w f ) + 4 π 2 [ m ( p i ) m ( p w f ) ] n = 1 q D F n ( 2 n 1 ) 2
m ¯ 1 ( p ) = m ¯ F ( p ) + 4 π 2 [ m ( p i ) m ¯ F ( p ) ] n = 1 q D 1 n ( 2 n 1 ) 2
m ¯ 2 ( p ) = m ¯ 1 ( p ) + 4 π 2 [ m ( p i ) m ¯ 1 ( p ) ] n = 1 q D 2 n ( 2 n 1 ) 2
We define the productivity index ( J ) and transmissibility ( T r 1 F and T r 21 ) as,
J = π 2 4 q i F m ( p i ) m ( p w f )
T r 1 F = π 2 4 q i 1 m ( p i ) m ¯ 1 ( p )
T r 21 = π 2 4 q i 2 m ( p i ) m ¯ 2 ( p )
where m ¯ 1 ( p ) , m ¯ 2 ( p ) is the average pseudo-pressure in hydraulic fracture region, region 1, and region 2, respectively. q i F = T s c Z s c 2 T P s c k i F A F m ( p i ) L F is the initial production rate from the hydraulic fracture region. q i 1 = T s c Z s c 2 T P s c k i 1 A 1 m ¯ 1 ( p ) L 1 is the initial production rate from the region 1. q i 2 = T s c Z s c 2 T P s c k i 2 A 2 m ¯ 2 ( p ) L 2 is the initial production rate from the region 2.
Substituting Equations (62)–(64) into Equations (50), (55), and (60), there are six ordinary differential terms on the left side of the system of equations, such as d q F n d t a F , d q F n d t a 1 , d q F n d t a 2 , d q 1 n d t a 1 , d q 1 n d t a 2 , d q 2 n d t a 2 . In order to solve the system of ordinary differential equations, an approximate pseudo-time t ˜ a is introduced. Therefore, the Equations (50), (55), and (60) can be rewritten as follows.
n = 1 d q 2 n d t ˜ a = T r 21 C 2 V p 2 + T r 21 C 1 V p 1 n = 1 ( 2 n 1 ) 2 q 2 n + T r 21 C 1 V p 1 n = 1 ( 2 n 1 ) 2 q 1 n
n = 1 d q 1 n d t ˜ a = T r 1 F C 1 V p 1 n = 1 ( 2 n 1 ) 2 q 2 n T r 1 F C 1 V p 1 + T r 1 F C F V p F n = 1 ( 2 n 1 ) 2 q 1 n + T r 1 F C F V p F n = 1 ( 2 n 1 ) 2 q F n
n = 1 d q F n d t ˜ a = J C F V p F n = 1 ( 2 n 1 ) 2 q F n + J C F V p F n = 1 ( 2 n 1 ) 2 q 1 n
where q D F n , q D 1 n , and q D 2 n are the initial production rate from the nth mode in hydraulic fracture region, region 1, and region 2, respectively. t ˜ a is an approximate pseudo-time, which is defined as the average of the t a 1 , t a 2 , and t a F .
Then we define three time-constant parameters as,
τ F = C F V p F J
τ 1 = C 1 V p 1 T r 1 F
τ 2 = C 2 V p 2 T r 21
We can rewrite this set of ODEs in matrix form,
d q 2 n d t ˜ a d q 1 n d t ˜ a d q F n d t ˜ a = ( 2 n 1 ) 2 1 τ 2 + T r 21 τ 1 T r 1 F T r 21 τ 1 T r 1 F 0 1 τ 1 1 τ 1 + T r 1 F τ F J T r 1 F τ F J 0 1 τ F 1 τ F q 2 n q 1 n q F n
where the initial conditions for the system of equations are,
q F ( t a = 0 ) = q i F
q 1 ( t a = 0 ) = 0
q 2 ( t a = 0 ) = 0
After solving Equation (70), we can get the nth flow rate in combination with the initial conditions. The flow rate is the summation of them. By converting the summation to an indefinite integral, the analytical solution in real-time space can be derived as follows (details are shown in Appendix B),
q F = β 3 q i F a 3 e λ 1 t ˜ a β 2 q i F a 6 e λ 2 t ˜ a + β 1 q i F a 9 e λ 3 t ˜ a + β 3 a 3 q i F π 4 λ 1 t ˜ a e r f c 3 λ 1 t ˜ a β 2 a 6 q i F π 4 λ 2 t ˜ a e r f c 3 λ 2 t ˜ a + β 1 a 9 q i F π 4 λ 3 t ˜ a e r f c 3 λ 3 t ˜ a
In Equation (78), the defined parameters are,
β 1 = a 1 a 1 a 5 a 2 a 4 a 1 a 9 a 3 a 7 a 1 a 5 a 2 a 4 a 1 a 8 a 2 a 7 a 1 a 6 a 3 a 4
β 2 = a 1 a 8 a 2 a 7 a 1 a 5 a 2 a 4 β 1
β 3 = β 2 a 4 a 7 a 1 β 1
where λ 1 , λ 2 , and λ 3 are the eigenvalues of matrix in the Equation (74), a 1 ~ a 9 are all the values in the eigenvector of matrix in the Equation (75), β 1 , β 2 , and β 3 are all the coefficients.
In our model, shale gas in region 2 must firstly flow into region 1 and then into the hydraulic fracture region and finally into the wellbore. Therefore, the flow rate in hydraulic fracture is equal to that in the wellbore. From the final solution, we can find that the flow rate is related to six variables. Through fitting the production data, these variables can be obtained and then the solution can be used for production analysis and forecasting.

3.2. Model Validation with Numerical Models

In order to verify the derived analytical solution, a numerical model is built with the commercial Eclipse reservoir simulator for comparing with the previous physical model, which is one-quarter of the volume around one hydraulic fracture. The numerical model has 27 grid cells in the x-direction, 50 grid cells in the y-direction, and only one grid cell in the z-direction. In order to capture the transient flow towards the hydraulic fracture, local grid refinement in x-direction is constructed. The top view of the model is shown in Figure 2, where the first column of grids represents half-length of hydraulic fracture, and the horizontal well located in the first row of the grids along the x-direction. The blue region represents the region 2 in which the gas adsorption plays an important role, while the red one is the region 1. Table 1 summarizes the input parameters used in the numerical models, which include the reservoir conditions, hydraulic conductivity, gas adsorption, and non-Darcy flow parameters.
In deriving the new analytical solution in this study, the pressure dependence of gas properties is considered using pseudo-pressure and pseudo-time. However, the results from the numerical models are gas rate versus real time. Therefore, the necessary step is to transform the numerical simulation results of gas rate with time into pseudo-time and then fit with the new model. Figure 3 presents the result analysis. According to the previous definition, the relationship between the pseudo-time and real time is shown in Figure 3a. The plot of 1/q vs. square root pseudo-time for regime diagnosis is shown in Figure 3b. The comparison of production rates obtained from numerical simulation and our analytical model is presented in Figure 3c. The blue dotted line represents the relationship between gas rate and pseudo-time, whereas the red one indicates that from the derived analytical solution. It can be seen that the results from the simulator and the analytical solution agree well with each other. Due to the high velocity gas flow in hydraulic facture region, the time constant in hydraulic fracture (0.001 days) is too short to be observed. Therefore, four flow regimes are identified in Figure 4. Regime 1 exhibits a half-slope straight line on the log-log plot and represents the transient linear flow in region 1. The permeability in this region is higher and hence the time constant in this region is shorter and nearly 22 days. Then an exponential curve of regime 2 indicates that the boundary of region 1 is reached, which is called the boundary-dominated flow in region 1 or inner-boundary-dominated flow. After that, the pressure propagates into the region 2. As for regime 3, it still presents an expected straight line with a half-slope signature. In our model, the permeability of region 2 is low-permeability matrix and thus the time constant is relatively longer (209 days). Regime 4 is the outer-boundary-dominated flow, which is affected by the boundary of region 2. Table 2 summarizes the four model parameters after fitting.
Based on the output parameters, we can predict the values of physical parameters to further validate our model according to the following step-by-step procedure:
Step 1: Calculate the productivity index.
As defined above, we can calculate the productivity index J combined with the values of initial rate q i F initial pseudo-pressure m ( p i ) and pseudo-bottom-hole pressure m ( p w f ) .
Step 2: Calculate the transmissibility between fractures and region 1 and region 1 and region 2.
Among the output parameters, we can obtain the ratio of transmissibility ( T r 21 / T r 1 F , T r 1 F / J ). Considering the calculated value in Step 1, the transmissibility can be calculated ( T r 21 , T r 1 F ).
Step 3: Calculate the pore volume of hydraulic fracture region, region 1, and region 2.
By transforming the formula that we defined in Equations (71)–(73), the pore volumes of different regions ( V p F , V p 1 , V p 2 ) can be obtained by:
V p F = J τ F C F   V p 1 = τ 1 T r 1 F C 1   V p 2 = τ 2 T r 21 C 2
Following the above steps, we calculate the physical parameters in the numerical case as shown in Table 3 to compare with the given data. It shows that our model solution is correct within the accepted error bound. Among them, V g and V c express the given volume and calculated volume, respectively.

3.3. Irregular Stimulated Region with One Hydraulic Fracture

In this section, three numerical cases with only one hydraulic fracture are designed for investigating the applicability of the analytical model for irregular stimulated region in Figure 4. For the three numerical cases, the input parameters are the same with the Case 1 except for the model dimensions. These three cases are identical with 31 grid cells in the x-direction, 51 grid cells in the y-direction, and only one grid cell in the z-direction, which represent the volume of 214 × 521 × 10 ft3. For Case 2, there is a rectangular stimulated region (region 1). As for Case 3, the stimulated region is irregular but symmetrical. In order to better describe the real condition, the stimulated region in Case 4 is designed to be neither regular nor symmetrical. The D-factor in hydraulic fracture is set as 0.0012 to represent the high-velocity non-Darcy flow.
In general, the results from the simulator and the analytical solution fit very well, as shown in Figure 5. There are also four flow regimes in the three cases. Comparing with Case 1 and Case 2, there are deviations for the cases with irregular stimulated region. The deviations are caused by the irregular inner boundary. Due to the irregular shapes of region 1, the time when the pressure propagates to inner boundaries changes with the distance away from the wellbore in different parts of inner boundary. However, the analytical solution is derived based on the regular inner boundary conditions. For the same reason, the variable starting time in regime 2 results in the deviation in the early time of the third flow regime. The flow time in regime 3 is long enough, thus, the curves in this regime fit well with each other at late time. For the three cases, the outer boundaries are identical and therefore the curves in the fourth regime also agree well with each other. According to the results summarized in Table 4, we can observe that the estimated results shown in Figure 5 are acceptable within engineering accuracy. In conclusion, our analytical solution is feasible for both regular and irregular region conditions.

3.4. Irregular Stimulated Regions with Several Hydraulic Fractures

In order to further validate the applicability in irregular region of the new derived model, three more numerical cases of multifractured horizontal wells are designed, as shown in Figure 6. These three cases have the identical dimension with 163 grid cells in the x-direction, 63 grid cells in the y-direction, and only one grid cell of 10 ft in the z-direction. The length of the horizontal well is 762 ft with 10 hydraulic fractures equally spaced along the x-direction. For Case 5, the stimulated regions are all the same regular and symmetrical regions. As a comparison, the stimulated regions in Case 6 are all irregular but symmetrical regions. Meanwhile, there is no interference between fractures in Case 5 and Case 6 and the length of the hydraulic fractures is 352 ft in Case 5 and 464 ft in Case 6. As for Case 7, all the stimulated regions around each hydraulic fracture are different, irregular, and asymmetric. The length of 10 hydraulic fractures ranges from 288 to 432 ft and there is interference between fractures. The other input parameters are the same with Case 1.
Figure 7 shows the comparison of gas rate versus pseudo-time obtained from numerical simulations and the new analytical solutions. It can be seen that the results agree very well. Four regimes are also identified. The linear flow time in region 1 based on the fitting results is 9, 7, and 5 days for the three cases, respectively. Then it displays the boundary flow of region 1, which lasts until about the 100th day. The time constant in region 2 is 72, 79, and 91 days, respectively. Finally, the boundary flow of region 2 lasts for hundreds of days. We can obtain that the decline rate of gas production is faster during regimes 1 and 2 and thus the regimes 3 and 4 can last for a longer time. Therefore, for the shale gas reservoir, it is crucial to enhance the volume of stimulated regions to keep a longer high production period. Meanwhile, the contribution of gas adsorption mainly reflects in regimes 3 and 4. It helps to extend the stable production period. According to the output parameters after fitting with the numerical cases, we can calculate the volume of the hydraulic fracture region, region 1, and region 2, as shown in Table 5. Considering that the relative errors meet the requirement of engineering accuracy, our new model is also suitable for the multifractured horizontal well.

3.5. Application to Field Case

The previous section demonstrated the accuracy of the derived analytical solution for production analysis. In this section, we apply the method to history matching and forecasting of field data in a shale gas reservoir. The gas well is selected for its relatively long production history and availability of pressure data. The main work flow is as follows:
  • Apply the given parameters and gas material balance equation to transform the time into pseudo-time and pressure into pseudo-pressure.
  • Make a diagnostic plot of production rate vs. pseudo-time.
  • Analyze the diagnostic plot to identify flow regimes.
  • Fit Equation (75) to the production data to obtain the model parameters τF, τ1, τ2, Tr21/Tr1F, Tr1F/J, and qiF.
  • Output the model parameters until a satisfactory matching is obtained.
  • Following the step-by-step procedure to calculate the volume of hydraulic fracture region, region 1, and region 2.
  • Make forecast of production rate with the model parameters.
We chose a horizontal well called Well B-15 with multiple fracture stages from Barnett shale [40]. The general data of the well for analysis is listed in Table 6. The well was produced under constant pressure and even for short early variable p/q the constant pressure was also applicable [41]. Firstly, we transformed the time into pseudo-time and the relationship between pseudo-time and time is shown in Figure 8a. The plot of 1/q vs. square root pseudo-time for regime diagnosis is shown in Figure 8b. Then the log-log diagnostic plot of gas rate versus pseudo-time was obtained, which exhibits a half-slope straight line in Figure 8c. According to the previous analysis, the flow time in matrix is longer and the decline rate in region 2 is slower. Therefore, the half-slope linear flow was diagnosed as the regime 3. Following the work flow, the next step was to match our model to the production data. The results of history matching and forecasting is shown in Figure 8c. The red marks in Figure 8c represent the production data and the green ones represent the analytical solution for history matching. It indicates a good matching. The six model parameters obtained from history matching are summarized in Table 7. Based on the parameters, we forecast the gas rate represented by the black markers in Figure 8c. According to the step-by-step procedure, we calculated the volume of the hydraulic fracture region as 89 ft3, the volume of stimulated region as 84.2 MMscf, and the volume of region 2 as 1784.2 MMscf. It is essential for increasing high production period by extending the volume of the stimulated region and decreasing the decline rate in the stimulated region.

4. Discussion

The results of this study can be summarized as four types: (1) derivation of the approximate analytical solution; (2) validation of the solution against different numerical models; (3) introducing a step-by-step procedure to predict the values of physical parameters; (4) application of the analytical model in the field case.
Result 1 is one of the novelty of the article. The model contains three regions and effectively accounts for non-Darcy flow in the hydraulic fractures, and gas adsorption and slippage in the matrix. During the derivation process, the governing PDEs are transformed into ordinary differential equations (ODEs) by integration, instead of the Laplace transform. There is no doubt that Laplace transform works extremely well. However, there are still some problems: (1) the first step for Laplace transform is dimensionless transformation. For some dimensionless variables such as dimensionless time, dimensionless pressure, dimensionless production, and so on, more inputs are needed and many are even estimated, which will lead to calculation error. (2) The numerical inversion is an essential step for converting Laplace space into real-time space. Among them, Stefest numerical inversion algorithm is most commonly used. In this algorithm, the number of inversion items N is uncertain. The improper value will result in deviations in real time. Certainly, Result 1 also has some imperfections. For example, a dual-porosity model with Knudsen diffusion [6] would be more representative and the heterogeneity deserves further study.
Result 2 is the key section in the article. Seven numerical cases were set for verification. According to the fitting results, it shows that the analytical solution is feasible for the irregular and asymmetric stimulated regions in a multifractured horizontal well. Considering that the shapes of the enhanced fracture regions are unknown in real cases, we creatively set three types of stimulated regions: regular region, irregular region, and very irregular region. At least, the validation results show that our model is robust.
Result 3 is another novelty of our work. It introduced a step-by-step procedure to calculate inversion parameters. Comparing with the given volumes and calculated volumes from case 1 to case 7, the model is also verified to meet the engineering requirements. Considering the simplifying assumptions, the model may need to be improved and the accurate microseismic data would be required to make further verification.
Result 4 is the most important section. Considering the decline characteristics of shale gas, the constant rate case is not as important as the constant pressure case for long term performance of tight/shale formations. Therefore, our model is derived based on the constant bottom-hole pressure condition, and namely one of limitations for our model is that it cannot be applied in a constant rate case. The second limitation is the data continuity, because the prerequisite for the model to make accurate prediction is great fitting. As for the discontinuous data like data missing, shut-in, or pressure/production jump, our model is not applicable. Combined with the assumption of single-phase flow in this work, the multiphase flow problem cannot be solved by this analytical solution. Gas–water two-phase flow is the most common in shale gas reservoirs, so our future work is to derive the new analytical solution for gas–water two-phase flow.

5. Conclusions

In this paper, we presented a practical analytical model to study the performance of MFHWs from shale gas reservoirs. Numerical models have been used to validate the analytical solutions and an excellent agreement was obtained. The following conclusions are drawn from this work:
  • A simple rate versus pseudo-time relationship is presented to account for transient linear and boundary-dominated flow periods in shale gas formation.
  • Incorporating the effect of gas adsorption, non-Darcy flow, and slippage flow in the analytical model by defining the modified pseudo-pressure and pseudo-time, accuracy is improved in production forecast in shale gas reservoir.
  • Comparing to the Laplace-transform solution, our analytical model is derived in real-time space and it is unnecessary to undertake dimensionless transformation and numerical inversion. It is more applicable in field scale.
  • Through the model parameters obtained from history matching the field data, the production rate and cumulative production can be forecasted. In addition, the pore volume of different regions can also be calculated by step-by-step procedure, which was validated to be feasible for the irregular and asymmetric stimulated regions in multifractured horizontal wells. According to the results, the calculation accuracy is less than 10% and meets the engineering requirements.

Author Contributions

Conceptualization, K.Q. and H.L.; methodology and investigation K.Q.; writing—original draft preparation K.Q.; writing—review and editing, H.L.; supervision, H.L.; project administration, H.L.; funding acquisition, H.L. All authors have read and agreed to the published version of the manuscript.

Funding

National Science and Technology Major Project of China: 2017ZX05039-005.

Acknowledgments

This work is funded by the National Science and Technology Major Project of China (Grant no. 2017ZX05039-005).

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

V the volume of the adsorbed gas, ft3
P the reservoir pressure, Psi
V L Langmuir volume
P L Langmuir pressure
C f rock compressibility, Psi−1
C w water compressibility, Psi−1
C g free gas compressibility, Psi−1
C g d adsorbed gas compressibility, Psi−1
C t total compressibility, Psi−1
C t modified total compressibility, Psi−1
B g ¯ gas reservoir volume factor
z * modified compressibility factor
z compressibility factor
Z s c standard compressibility factor
T reservoir temperature, K
T s c standard condition temperature, K
P s c standard condition pressure, Psi
p ¯ average reservoir pressure, Psi
k i F hydraulic fracture permeability, mD
k e q F equivalent hydraulic fracture permeability, mD
k a m matrix permeability, mD
k m intrinsic permeability, mD
k i 1 region 1 permeability, mD
k i 2 region 2 permeability, mD
r pore radius, ft
R the universal gas constant
M the gas molecular weight
K n Knudsen number
V b volume of the region
V p pore volume of the region
q i F initial production rate from the hydraulic fracture region
q i 1 initial production rate from region 1
q i 2 initial production rate from region 2
t production time, days
t a pseudo-time, days
t a F pseudo-time in fracture region, days
t a F pseudo-time in region 1, days
t a 2 pseudo-time in region 2, days
t ˜ a approximate pseudo-time, days
m 2 ( p ) pseudo-pressure in region 2
m 1 ( p ) pseudo-pressure in region 1
m F ( p ) pseudo-pressure in hydraulic region
m ¯ 2 ( p ) average pseudo-pressure in region 2
m ¯ 1 ( p ) average pseudo-pressure in region 1
m ¯ F ( p ) average pseudo-pressure in fracture region
w / 2 half-width of hydraulic fracture, ft
x 1 region 1 impact distance, ft
x 2 half distance between fractures, ft
y e half-length of macro-fracture, ft
z e depth of top reservoir, ft
z 0 depth of bottom reservoir, ft
τ 1 region 1 constant time, days
τ 2 region 2 constant time, days
τ F hydraulic region constant time, days
T r 21 transmissibility between region 1 and region 2, STB/D/psi
T r 1 F transmissibility between microfracture and matrix, STB/D/psi
J hydraulic fracture region production index, STB/D/psi
V p F pore volumes of hydraulic region
V p 1 pore volumes of region 1
V p 2 pore volumes of region 2
α correlation parameter
β non-Darcy flow coefficient
ϕ porosity
ρ m matrix density, g/cm3
μ fluid viscosity, cp
v gas flow velocity

Appendix A. Derivation of the Average Pseudo-Pressure

The solution to the 1-D linear oil flow with constant-pressure boundary condition is given by [27]:
p D = 1 + 4 π n = 1 ( 1 ) n 2 n 1 e ( 2 n 1 ) π 2 2 t D cos ( 2 n 1 ) π 2 x D
For the shale gas case, the pressure and time should be replaced by pseudo-pressure and pseudo-time shown in Equations (17) and (19). Following the same step, we derive the general solution applied for the 1-D linear gas flow with constant-pressure boundary condition.
m D ( p ) = 1 + 4 π n = 1 ( 1 ) n 2 n 1 e ( 2 n 1 ) π 2 2 t a D cos ( 2 n 1 ) π 2 x D
where
m D = m ( p i ) m ( p ) m ( p i ) m ( p w f ) is the dimensionless reservoir pseudo-pressure.
t a D = k t a ϕ μ c t L is the dimensionless pseudo-time.
x D = x e x x e x w f is the dimensionless distance in the x-direction.
Similarly, the production rate can be obtained as,
q D = 2 e ( 2 n 1 ) π 2 2 t a D
According to the definition of volume-weighted average pseudo-pressure,
m ¯ ( p ) = V m ( p ) d V V d V
Therefore,
m ¯ D ( p ) = 1 x D 1 + 4 π n = 1 ( 1 ) n 2 n 1 e ( 2 n 1 ) π 2 2 t a D cos ( 2 n 1 ) π 2 x D d x D 1 x D d x D
Performing the integration results in,
m ¯ D ( p ) = 1 + 8 π 2 n = 1 ( 1 ) n ( 2 n 1 ) 2 e ( 2 n 1 ) 2 π 2 4 t a D sin ( 2 n 1 ) π 2
For odd values of n, sin ( 2 n 1 ) π 2 is positive unity and ( 1 ) n is negative unity. As a result, the product of these two coefficients is always negative. Equation (A6) can be rewritten as,
m ¯ D ( p ) = 1 8 π 2 n = 1 1 ( 2 n 1 ) 2 e ( 2 n 1 ) 2 π 2 4 t a D
Substituting the Equation (A3) into the Equation (A7), we can obtain that,
m ¯ D ( p ) = 1 4 π 2 n = 1 q D ( 2 n 1 ) 2
Transforming the Equation (A8) into dimensional form, Equation (A8) can be rewritten as,
m ¯ ( p ) = m ( p w f ) + 4 π 2 [ m ( p i ) m ( p w f ) ] n = 1 q D n ( 2 n 1 ) 2

Appendix B. Solution of the System of ODEs

To solve the systems of ODEs in Equation (71), we extract the coefficient matrix as,
A n = 1 τ 2 + T r 21 τ 1 T r 1 F T r 21 τ 1 T r 1 F 0 1 τ 1 1 τ 1 + T r 1 F τ F J T r 1 F τ F J 0 1 τ F 1 τ F 3 × 3
For the 3 × 3 matrix, we can obtain three eigenvalues λ 1 , λ 2 , λ 3 and three eigenvectors ξ 1 , ξ 2 , ξ 3 where
ξ 1 = a 1 a 2 a 3 ξ 2 = a 4 a 5 a 6 ξ 3 = a 7 a 8 a 9
Combining the eigenvalues and eigenvectors, we obtain the solution of the ODEs,
q 2 n q 1 n q F n = ( 2 n 1 ) 2 C 1 a 1 a 2 a 3 e ( 2 n 1 ) 2 λ 1 t a + C 2 a 4 a 5 a 6 e ( 2 n 1 ) 2 λ 2 t a + C 3 a 7 a 8 a 9 e ( 2 n 1 ) 2 λ 3 t a
where C 1 , C 2 , and C 3 are unknown constants.
According to the initial assumption in Equations (72)–(74), we can calculate the value of C1, C2, and
C 3 C 1 = β 3 q i F   C 2 = β 2 q i F   C 3 = β 1 q i F
Where
β 1 = a 1 a 1 a 5 a 2 a 4 a 1 a 9 a 3 a 7 a 1 a 5 a 2 a 4 a 1 a 8 a 2 a 7 a 1 a 6 a 3 a 4
β 2 = a 1 a 8 a 2 a 7 a 1 a 5 a 2 a 4 β 1
β 3 = β 2 a 4 a 7 a 1 β 1
Therefore, we obtain the function of the final rate,
q F n = ( 2 n 1 ) 2 β 3 a 3 q i F e ( 2 n 1 ) 2 λ 1 t a β 2 a 6 q i F e ( 2 n 1 ) 2 λ 2 t a + β 1 a 9 q i F e ( 2 n 1 ) 2 λ 3 t a
To obtain the total production rate, we summate the Equation (A14),
q F = n = 1 q F n = β 3 a 3 q i F e λ 1 t a β 2 a 6 q i F e λ 2 t a + β 1 a 9 q i F e λ 3 t a + n = 1 ( 2 n 1 ) 2 β 3 a 3 q i F e ( 2 n 1 ) 2 λ 1 t a β 2 a 6 q i F e ( 2 n 1 ) 2 λ 2 t a + β 1 a 9 q i F e ( 2 n 1 ) 2 λ 3 t a
Through converting the summation to an integral,
q F = β 3 a 3 q i F e λ 1 t a β 2 a 6 q i F e λ 2 t a + β 1 a 9 q i F e λ 3 t a + lim z 3 λ 1 t a z β 3 a 3 q i F e z 2 t a β 2 a 6 q i F e λ 2 λ 1 t a + β 1 a 9 q i F e λ 3 λ 1 t a d z 2 λ 1 t a
Adopting the complementary error function, the rate versus pseudo-time function can be
q F = β 3 a 3 q i F e λ 1 t a β 2 a 6 q i F e λ 2 t a + β 1 a 9 q i F e λ 3 t a + β 3 a 3 q i F π 4 λ 1 t a e r f c 3 λ 1 t a β 2 a 6 q i F π 4 λ 2 t a e r f c 3 λ 2 t a + β 1 a 9 q i F π 4 λ 3 t a e r f c 3 λ 3 t a

References

  1. Hu, X.Y.; Li, M.T.; Peng, C.G.; Davarpanah, A. Hybrid Thermal-Chemical Enhanced Oil Recovery Methods—An Experimental Study for Tight Reservoirs. Symmetry 2020, 12, 947. [Google Scholar] [CrossRef]
  2. Jin, Y.; Davarpanah, A. Using Photo-Fenton and Floatation Techniques for the Sustainable Management of Flow-Back Produced Water Reuse in Shale Reservoirs Exploration. Water Air Soil Pollut. 2020, 231, 1–8. [Google Scholar] [CrossRef]
  3. Hu, X.; Xie, J.; Cai, W.; Wang, R.; Davarpanah, A. Thermodynamic Effects of Cycling Carbon Dioxide Injectivity in Shale Reservoirs. J. Pet. Sci. Eng. 2020, 195, 107717. [Google Scholar] [CrossRef]
  4. Wu, M.L.; Ding, M.C.; Yao, J.; Li, C.; Huang, Z.; Xu, S. Production-Performance Analysis of Composite Shale-Gas reservoirs by the Boundary-Element Method. SPE Reserv. Eval. Eng. 2019, 22, 238–252. [Google Scholar] [CrossRef] [Green Version]
  5. Sun, S.H.; Zhou, M.H.; Lu, W.; Davarpanah, A. Application of Symmetry Law in Numerical Modeling of Hydraulic Fracturing by Finite Element Method. Symmetry 2020, 12, 1122. [Google Scholar] [CrossRef]
  6. Sun, H.; Chawathe, A.; Hoteit, H.; Shi, X.; Li, L. Understanding Shale Gas Flow Behavior Using Numerical Simulation. SPE J. 2015, 20, 142–154. [Google Scholar] [CrossRef]
  7. Zhang, M.; Sun, Q.; Ayala, L.F. Semi-Analytical Modeling of Multi-Mechanism Gas Transport in Shale Reservoirs with Complex Hydraulic-Fracture Geometries by the Boundary Element Method. In Proceedings of the SPE Annual Technical Conference and Exhibition, Calgary, AB, Canada, 30 September–2 October 2019. SPE Paper 195902-MS. [Google Scholar]
  8. Davarpanah, A.; Shirmohammadi, R.; Mirshekari, B.; Aslani, A. Analysis of hydraulic fracturing techniques: Hybrid fuzzy approaches. Arab. J. Geosci. 2019, 12, 402. [Google Scholar] [CrossRef]
  9. Nobakht, M.; Clarkson, C.R. A New Analytical Method for Analyzing Linear Flow in Tight/Shale Gas Reservoirs: Constant-Flowing-Pressure Boundary Condition. SPE Reserv. Eval. Eng. 2012, 15, 370–384. [Google Scholar] [CrossRef] [Green Version]
  10. Zhao, Y.L.; Zhang, L.H.; Luo, J.X.; Zhang, B.-N. Performance of fractured horizontal well with stimulated reservoir volume in unconventional gas reservoir. J. Hydrol. 2014, 512, 447–456. [Google Scholar] [CrossRef]
  11. Clarkson, C.R.; Qanbari, F. An Approximate Semianalytical Multiphase Forecasting Method for Multifractured Tight Light-Oil Wells With Complex Fracture Geometry. J. Can. Pet. Technol. 2015, 54, 489–508. [Google Scholar] [CrossRef]
  12. Wei, S.M.; Xia, Y.; Jin, Y. Quantitative study in shale gas behaviors using a coupled triple-continuum and discrete fracture model. J. Pet. Sci. Eng. 2019, 174, 49–69. [Google Scholar] [CrossRef]
  13. Wang, Q.; Wan, J.F.; Mu, L.F.; Shen, R.; Jurado, M.J.; Ye, Y. An Analytical Solution for Transient Productivity Prediction of Multi-Fractured Horizontal Wells in Tight Gas Reservoirs Considering Nonlinear Porous Flow Mechanisms. Energies 2020, 13, 1066. [Google Scholar] [CrossRef] [Green Version]
  14. Zhang, Q.; Jiang, S.; Wu, X.Y. Development and Calibration of a Semianalytic Model for Shale Wells with Nonuniform Distribution of Induced Fractures Based on ES-MDA Method. Energies 2020, 13, 3718. [Google Scholar] [CrossRef]
  15. Lee, S.T.; Brockenbrough, J. A New Analytic Solution for Finite Conductivity Vertical Fractures With Real Time and Laplace Space Parameter Estimation. In Proceedings of the SPE Annual Technical Conference and Exhibition, San Francisco, CA, USA, 5–8 October 1983. SPE Paper 12013-MS. [Google Scholar]
  16. Wattenbarger, R.A.; El-Banbi, A.H.; Villegas, M.E.; Maggard, J.B. Production Analysis of Linear Flow into Fractured Tight Gas Wells. In Proceedings of the SPE Rocky Mountain Regional/Low-Permeability Reservoirs Symposium, Denver, CO, USA, 5–8 April 1998. SPE Paper 39931-MS. [Google Scholar]
  17. El-Banbi, A.H.; Wattenbarger, R.A. Analysis of Linear Flow in Gas Well Production. In Proceedings of the SPE Gas Technology Symposium, Calgary, AB, Canada, 15–18 March 1998. SPE Paper 39972-MS. [Google Scholar]
  18. Bello, R.O.; Wattenbarger, R.A. Multi-Stage Hydraulically Fractured Horizontal Shale Gas Well Rate Transient Analysis. In Proceedings of the North Africa Technical Conference and Exhibition, Cairo, Egypt, 14–17 February 2010. SPE Paper 126754-MS. [Google Scholar]
  19. Ozkan, E.; Brown, M.L.; Raghavan, R.S.; Kazemi, H. Comparison of Fractured Horizontal-Well Performance in Conventional and Unconventional Reservoirs. In Proceedings of the SPE Western Regional Meeting, San Jose, CA, USA, 24–26 March 2009. SPE Paper 121290-MS. [Google Scholar]
  20. Brown, M.; Ozkan, E.; Raghavan, R.; Kazemi, H. Practical Solutions for Pressure-Transient Responses of Fractured Horizontal Wells in Unconventional Shale Reservoirs. SPE Reserv. Eval. Eng. 2011, 14, 663–676. [Google Scholar] [CrossRef]
  21. Stalgorova, E.; Mattar, L. Practical Analytical Model to Simulate Production of Horizontal Wells with Branch Fractures. In Proceedings of the SPE Canadian Unconventional Resources Conference, Calgary, AB, Canada, 30 October–1 November 2012. SPE Paper 162515-MS. [Google Scholar]
  22. Heidari, S.M.; Clarkson, C.R. An Analytical Model for Analyzing and Forecasting Production from Multi-fractured Horizontal Wells with Complex Branched-Fracture Geometry. SPE Reserv. Eval. Eng. 2015, 18, 356–374. [Google Scholar] [CrossRef]
  23. Abassi, M.; Sharifi, M.; Kazemi, A. Development of New Analytical Model for Series and Parallel Triple Porosity Models and Providing Transient Shape Factor between Different Regions. J. Hydrol. 2019, 574, 683–698. [Google Scholar] [CrossRef]
  24. Yao, F. Production Evaluation of Fracturing Horizontal Wells for Solution Gas Drive in Tight Oil Reservoirs; China University of Petroleum (Beijing): Beijing, China, 2018. [Google Scholar]
  25. Stehfest, H. Algorithm 368: Numerical inversion of Laplace transforms [D5]. Commun. ACM 1970, 13, 47–49. [Google Scholar] [CrossRef]
  26. Shahamat, M.S.; Mattar, L.; Aguilera, R. A Physics-Based Method for Production Data Analysis of Tight and Shale Petroleum Reservoirs Using Succession of Pseudo-Steady States. In Proceedings of the SPE/EAGE European Unconventional Resources Conference and Exhibition, Vienna, Austria, 25–27 February 2014. SPE Paper 167686-MS. [Google Scholar]
  27. Ogunyomi, B.A. Simple Mechanistic Modeling of Recovery from Unconventional Oil Reservoirs; The University of Texas at Austin: College Station, TX, USA, 2015. [Google Scholar]
  28. Ogunyomi, B.A.; Patzek, T.W.; Lake, L.W.; Kabir, C.S. History Matching and Rate Forecasting in Unconventional Oil Reservoirs With an Approximate Analytical Solution to the Double-Porosity Model. SPE Reserv. Eval. Eng. 2016, 19, 70–82. [Google Scholar] [CrossRef]
  29. Qiu, K.X.; Li, H. A New Analytical Solution of the Triple-Porosity Model for History Matching and Performance Forecasting in Unconventional Oil Reservoirs. SPE J. 2018, 23, 2060–2079. [Google Scholar] [CrossRef]
  30. Qiu, K.X.; Li, H. A New Analytical Model for Production Forecasting in Unconventional Reservior Considering the Simultaneous Matrix-fracture flow. In Proceedings of the Asia Pacific Unconventional Resource Technology Conference, Brisbane, Australia, 18–19 November 2019. SPE Paper URTEC-198266-MS. [Google Scholar]
  31. Tabatabaie, S.H. Unconventional Reservoirs: Mathematical Modeling of Some Non-linear Problems; University of Calgary: College Station, AB, Canada, 2014. [Google Scholar]
  32. Behmanesh, H.; Hamdi, H.; Clarkson, C.R.; Thompson, J.M.; Anderson, D.M. Analytical Modeling of Linear Flow in Single-Phase Tight Oil and Tight Gas Reservoirs. J. Pet. Sci. Eng. 2018, 171, 1084–1098. [Google Scholar] [CrossRef]
  33. Anderson, D.M.; Nobakht, M.; Moghadam, S.; Mattar, L. Analysis of Production Data from Fractured Shale Gas Wells. In Proceedings of the SPE Unconventional Gas Conference, Pittsburgh, PA, USA, 23–25 February 2010. SPE Paper 131787. [Google Scholar]
  34. Langmuir, I. The Constitution and Fundamental Properties of Solids and Liquids. J. Am. Chem. Soc. 1916, 38, 2221–2295. [Google Scholar] [CrossRef] [Green Version]
  35. Bumb, A.C.; McKee, C.R. Gas-Well Testing in the Presence of Desorption for Coal bed Methane and Devonian Shale. In Proceedings of the SPE Unconventional Gas Technology Symposium, Louisville, KY, USA, 18–21 May 1986. SPE Paper 15227-MS. [Google Scholar]
  36. King, G.R. Material Balance Techniques for Coal Seam and Devonian Shale Gas Reservoirs. In Proceedings of the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 23–26 September 1990. SPE Paper 20730-MS. [Google Scholar]
  37. Forchheimer, P. Wasserbewegung Durch Boden. Z. Ver. Deutsch Ing. 1901, 45, 1782–1788. [Google Scholar]
  38. Wang, J.; Lou, H.S.; Liu, H.Q. An Intergrative Model to Simulate Gas Transport and Production Coupled With Gas Adsorption, Non-Darcy Flow, Surface Diffusion, and Stress Dependence in Organic-Shale Reservoir. SPE J. 2017, 22, 244–264. [Google Scholar] [CrossRef]
  39. Wang, H.Y.; Porcu, M.M. Impact of Shale-Gas Apparent Permeability on Production: Combined Effects of Non-Darcy Flow/Gas Slippage, Desorption and Geomechanics. SPE Reserv. Eval. Eng. 2015, 18, 495–507. [Google Scholar] [CrossRef]
  40. Samandarli, O. A New Method for History Matching and Forecasting Shale Gas/Oil Reservoir Production Performance with Dual and Triple Porosity Models; Texas A&M University: College Station, TX, USA, 2011. [Google Scholar]
  41. Hamdi, H.; Behmanesh, H.; Clarkson, C.R. A Semi-Analytical Approach for Analysis of Wells Exhibiting Multi-Phase Transient Linear Flow: Application to Field Data. In Proceedings of the SPE Annual Technical Conference and Exhibition, Calgary, AB, Canada, 30 September–2 October 2019. SPE Paper 196164-MS. [Google Scholar]
Figure 1. Schematic of a multistage fractured horizontal well with two regions around the hydraulic fracture. (a) The 3D model. (b) The plan view.
Figure 1. Schematic of a multistage fractured horizontal well with two regions around the hydraulic fracture. (a) The 3D model. (b) The plan view.
Energies 13 05899 g001
Figure 2. Reservoir grid for the numerical model. The red grids are the high-permeability region and the blue ones are the low-permeability region.
Figure 2. Reservoir grid for the numerical model. The red grids are the high-permeability region and the blue ones are the low-permeability region.
Energies 13 05899 g002
Figure 3. The result analysis for Case 1. (a) Comparison of the relationship between pseudo-time and time. (b) The plot 1/q vs. square root pseudo-time for regime diagnosis. (c) Comparison of the relationship between gas rate and pseudo-time. (Regime 1: transient linear flow in region 1 with slope -1/2. Regime 2: boundary-dominated flow for region 1. Regime 3: transient linear flow in region 2 with slope -1/2. Regime 4: boundary-dominated flow for region 2).
Figure 3. The result analysis for Case 1. (a) Comparison of the relationship between pseudo-time and time. (b) The plot 1/q vs. square root pseudo-time for regime diagnosis. (c) Comparison of the relationship between gas rate and pseudo-time. (Regime 1: transient linear flow in region 1 with slope -1/2. Regime 2: boundary-dominated flow for region 1. Regime 3: transient linear flow in region 2 with slope -1/2. Regime 4: boundary-dominated flow for region 2).
Energies 13 05899 g003
Figure 4. Reservoir grid for the numerical cases showing different shapes of stimulated region. The red grids are the stimulated region and the blue ones are the low-permeability region.
Figure 4. Reservoir grid for the numerical cases showing different shapes of stimulated region. The red grids are the stimulated region and the blue ones are the low-permeability region.
Energies 13 05899 g004
Figure 5. Comparison of the relationship of gas rate vs. pseudo-time from the numerical cases and new analytical solution. (a) Comparison of the relationship of gas rate vs. pseudo-time from the Case 2 and new model. (b) Comparison of the relationship of gas rate vs. pseudo-time from the Case 3 and new model. (c) Comparison of the relationship of gas rate vs. pseudo-time from the Case 4 and new model.
Figure 5. Comparison of the relationship of gas rate vs. pseudo-time from the numerical cases and new analytical solution. (a) Comparison of the relationship of gas rate vs. pseudo-time from the Case 2 and new model. (b) Comparison of the relationship of gas rate vs. pseudo-time from the Case 3 and new model. (c) Comparison of the relationship of gas rate vs. pseudo-time from the Case 4 and new model.
Energies 13 05899 g005
Figure 6. Reservoir grid for multifractured horizontal well (MFHWs).
Figure 6. Reservoir grid for multifractured horizontal well (MFHWs).
Energies 13 05899 g006
Figure 7. Comparison of the relationship of gas rate vs. pseudo-time from the numerical cases and new analytical solution. (a) Comparison of the relationship of gas rate vs. pseudo-time from the Case 5 and new model. (b) Comparison of the relationship of gas rate vs. pseudo-time from the Case 6 and new model. (c) Comparison of the relationship of gas rate vs. pseudo-time from the Case 7 and new model.
Figure 7. Comparison of the relationship of gas rate vs. pseudo-time from the numerical cases and new analytical solution. (a) Comparison of the relationship of gas rate vs. pseudo-time from the Case 5 and new model. (b) Comparison of the relationship of gas rate vs. pseudo-time from the Case 6 and new model. (c) Comparison of the relationship of gas rate vs. pseudo-time from the Case 7 and new model.
Energies 13 05899 g007
Figure 8. Summary of production profiles for the field example. (a) The relationship between time and pseudo-time. (b) The plot 1/q vs. square root pseudo-time for regime diagnosis. (c) The results of history-matching and forecasting on log-log plot.
Figure 8. Summary of production profiles for the field example. (a) The relationship between time and pseudo-time. (b) The plot 1/q vs. square root pseudo-time for regime diagnosis. (c) The results of history-matching and forecasting on log-log plot.
Energies 13 05899 g008
Table 1. Summary of input parameters for the synthetic numerical model.
Table 1. Summary of input parameters for the synthetic numerical model.
ParameterValue
Model dimension (X × Y × Z) (ft)80 × 500 × 10
Initial pressure (Psi)2500
Bottom-hole pressure (Psi)500
Temperature (K)318.15
Langmuir volume (Mscf/ton)0.096
Langmuir pressure (Psi)650
Compressibility (10−6 Psi−1)7.5
Matrix porosity0.06
Fracture porosity0.08
Permeability of hydraulic fracture (mD)500
Permeability in region 2 (mD)0.0005
Permeability in region 1 (mD)0.01
Fracture width (ft)0.1
D-factor (Day/Mscf)0.0012
Table 2. Model parameters used in the analytical model for validation against numerical simulation.
Table 2. Model parameters used in the analytical model for validation against numerical simulation.
ParameterValue
τF (days)0.001
τ1 (days)22
τ2 (days)209
Tr21/Tr1F0.09
Tr1F/J1.18
qiF (Mscf/D)109
Table 3. Inferred parameter values according to the fitting results.
Table 3. Inferred parameter values according to the fitting results.
ParameterGiven DataCalculated ValueRelative Error
(Vg-Vc)/ Vg × 100%
VpF (ft3)4038.374.1%
Vp1 (ft3)597059540.3%
Vp2 (ft3)18,00017,4553.0%
Table 4. Inferred parameter values according to fitting results.
Table 4. Inferred parameter values according to fitting results.
CaseτFτ1τ2Tr21/Tr1FTr1F/JqiF
(Mscf/D)
Given Data
(ft3)
Calculated Data
(ft3)
Relative Error
Case 20.029970.041.6327VpF = 31VpF = 323%
Vp1 = 26,258Vp1 = 25,1374.3%
Vp2 = 37,512Vp2 = 35,3435.8%
Case 30.0191890.051.4252VpF = 31VpF = 296.5%
Vp1 = 20,084Vp1 = 21,1235.2%
Vp2 = 43,686Vp2 = 42,3873%
Case 40.02111060.051.3316VpF = 31VpF = 303.2%
Vp1 = 24,145Vp1 = 25,1554.2%
Vp2 = 39,625Vp2 = 40,1351.3%
Table 5. Inferred parameters value according to fitting results.
Table 5. Inferred parameters value according to fitting results.
CaseτFτ1τ2Tr21/Tr1FTr1F/JqiF
(Mscf/D)
Given Data
(ft3)
Calculated Data
(ft3)
Relative Error
Case 50.029720.031.6990VpF = 21VpF = 239.5%
Vp1 = 106,022Vp1 = 103,8202.1%
Vp2 = 105,281Vp2 = 98,0546.7%
Case 60.037790.032.4937VpF = 28VpF = 276.5%
Vp1 = 81,062Vp1 = 79,3022.2%
Vp2 = 130,234Vp2 = 125,1393.8%
Case 70.025910.051.9911VpF = 24VpF = 234.2%
Vp1 = 67,264Vp1 = 66,0971.7%
Vp2 = 144,036Vp2 = 151,9505.5%
Table 6. Summary of the reservoir parameters for analysis.
Table 6. Summary of the reservoir parameters for analysis.
ParameterValue
Initial pressure (Psi)2950
Bottom-hole pressure (Psi)480
Temperature (K)344.3
Langmuir volume (Scf/ton)96
Langmuir pressure (Psi)650
Compressibility (Psi−1)4 × 10−6
Porosity0.06
Table 7. Model parameters used in the analytical model for application in the field example.
Table 7. Model parameters used in the analytical model for application in the field example.
ParameterValue
τF (days)0.01
τ1 (days)11
τ2 (days)3105
Tr21 / Tr1F0.075
Tr1F / J0.86
qiF (Mscf/D)19,041
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Qiu, K.; Li, H. An Analytical Model for Production Analysis of Hydraulically Fractured Shale Gas Reservoirs Considering Irregular Stimulated Regions. Energies 2020, 13, 5899. https://doi.org/10.3390/en13225899

AMA Style

Qiu K, Li H. An Analytical Model for Production Analysis of Hydraulically Fractured Shale Gas Reservoirs Considering Irregular Stimulated Regions. Energies. 2020; 13(22):5899. https://doi.org/10.3390/en13225899

Chicago/Turabian Style

Qiu, Kaixuan, and Heng Li. 2020. "An Analytical Model for Production Analysis of Hydraulically Fractured Shale Gas Reservoirs Considering Irregular Stimulated Regions" Energies 13, no. 22: 5899. https://doi.org/10.3390/en13225899

APA Style

Qiu, K., & Li, H. (2020). An Analytical Model for Production Analysis of Hydraulically Fractured Shale Gas Reservoirs Considering Irregular Stimulated Regions. Energies, 13(22), 5899. https://doi.org/10.3390/en13225899

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop