Abstract
Due to the ultralow permeability of shale gas reservoirs, stimulating the reservoir formation by using hydraulic fracturing technique and horizontal well is required to create the pathway of gas flow so that the shale gas can be recovered in an economically viable manner. The hydraulic fractured formations can be divided into two regions, stimulated reservoir volume (SRV) region and non-SRV region, and the produced shale gas may exist as free gas or adsorbed gas under the initial formation condition. Investigating the recovery factor of different types of shale gas in different region may assist us to make more reasonable development strategies. In this paper, we build a numerical simulation model, which has the ability to take the unique shale gas flow mechanisms into account, to quantitatively describe the gas production characteristics in each region based on the field data collected from a shale gas reservoir in Sichuan Basin in China. The contribution of the free gas and adsorbed gas to the total production is analyzed dynamically through the entire life of the shale gas production by adopting a component subdivision method. The effects of the key reservoir properties, such as shale matrix, secondary natural fracture network, and primary hydraulic fractures, on the recovery factor are also investigated.
1. Introduction
With the growing demand on cleaner energy to sustain the global economy, the exploration and development of natural gas have attracted wide attention in the recent years. Shale gas, as a newly developed unconventional energy source, is expected to play a more important role. Shale formations are widely distributed in the earth’s crust. The Global Shale Gas Initiative estimates that there are 688 shale formations in 142 basins worldwide. This provides vast potential space for shale gas to be accumulated in biogenic manner, thermogenic manner, or the combined biogenic-thermogenic manner [1]. The Energy Information Administration (EIA) in 2013 estimated approximately 7,300 tcf shale gas resources in 137 shale formations in 41 countries that are technically recoverable, and 32% of the total estimated natural gas resources are in shale formation [2]. The most distinctive characteristics of shale gas reservoir are its extremely low permeability and porosity, which had made the shale formation be treated as source rock and seal previously. However, with the application of advanced technologies, such as multistage hydraulic fracturing and horizontal wells, it is possible to develop the organic-rich shale reservoir in an economic viable manner.
The flow mechanisms during the shale gas production are much more complicated than those in conventional gas reservoir due to the complex multiscale gas transport mechanisms in the stimulated shale gas reservoir. The shale gas can be adsorbed on the organic matter, and thus the desorption may be crucial for the ultimate gas recovery in shale formation if the organic content in the shale matrix is high [3]. The desorbed gas combined with free gas in intergranular pores needs to transport through the dominant nanopores in the shale matrix. Since the mean free pat of the gas molecule can be smaller than the pore radius in the nanopores, the diffusion can greatly improve the gas transport speed in this process which is now characterized by Knudsen diffusion. Many corrections to obtain the apparent permeability of shale matrix have been proposed to take the non-Darcy flow behavior and gas slippage effect into account [4–6]. The gas can then be travelled through a complex fracture network induced by hydraulic fracturing stimulation process, which creates a high permeable region called stimulated reservoir volume (SRV) [7, 8]. Mainly, two types of fractures are involved in this SRV region: the widely distributed reactivated preexisting natural fractures with complex geometry and the newly created sparsely distributed hydraulic fractures with relatively less complex geometry. When these two types of fracture system are well connected, the SRV region with enhanced permeability can greatly improve the well performance in low-permeable reservoirs [9, 10].
Numerical simulation of shale gas production is challenging. Generally, there exist two common conceptual approaches in the literatures: the dual/multiple porosity/permeability model and discrete fracture model. The dual-porosity model proposed by Warren and Root [11] has been widely used in well-connected small-scale fracture systems, such as natural fracture system [12–16]. In this approach, the reservoir is characterized by two overlapping continua to represent fractures and matrix. The flow interaction between these continua is described by using a transfer function called the shape factor. Due to the simplification on the complex geometry of the fractures, this approach cannot account explicitly for the effect of the individual fractures on the flow characteristics. Discrete fracture model method (DFM) can overcome the limitation of dual porosity method and take the realistic fracture geometry into account [17]. It is common to use unstructured grid during the simulation of the fluid flow process to make the grid system consistent with the geometry of each fracture. Due to the computational cost of this approach, DFM is usually applied when a small number of large-scale fractures play dominant roles in fluid flow. Recently, the embedded discrete fracture method (EDFM) has gained more interests when numerically characterizing the fractured reservoir [18, 19]. The conception of EDFM was proposed by Lee et al. [20, 21] by taking advantage of both dual-porosity model and discrete fracture model. EDFM characterizes the medium matrix by using structured grid and embeds the discrete fractures into the matrix grid by handling the intersection between fractures and matrix; thus all the issues associated with the mathematical difficulties of unstructured grid vanish.
The United States leads the shale gas production worldwide. EIA estimates that 15 trillion cubic feet of natural gas were produced in 2016 from shale gas wells, which takes up 47% of total natural gas production. China is the third country that commercializes the shale gas production after US and Canada, and Fuling shale gas play has been the dominant shale gas reservoir in China by far since its beginning of shale gas production in 2012. The shale gas production in China is still at its early stage, and the influence of many key reservoir properties and hydraulic fracturing design parameters on the well performance remains unclear. In this paper, we built a dual-porosity-based numerical model that is capable of considering the above special flow mechanisms to investigate contributions of free and adsorbed gas to the total production and the effect of the matrix permeability, the number of hydraulic fracturing stages, the conductivity of primary hydraulic fractures, half-length of primary hydraulic fractures, and the conductivity and development level of secondary fracture network on the recovery factors in SRV and non-SRV regions according to the field data collected from a shale gas field in Sichuan Basin. All these investigations are crucial to the reasonable design of development plan for shale gas production. A component subdivision method implemented in the UNCONG simulator is used to explicitly characterize the free gas and adsorbed gas produced from SRV and non-SRV regions, which makes the sensitivity analysis presented here unique and more complete compared with the previous works based on commercial simulators.
The rest of the paper is arranged as follows: the mathematical and numerical model to describe the gas flow process in the specified shale gas reservoir and the component subdivision method to differentiate the effect of free gas and adsorbed gas are presented in Section 2; the contributions of free gas and adsorbed gas to the total gas production are analyzed by using “component subdivision method” and the influences of reservoir properties and hydraulic fracturing parameters on the well performance are studied in Section 3; and finally conclusions are drawn based on all the analysis results in Section 4.
2. Methodology
2.1. Mathematical Model of the Fluid Flow
The entire shale gas reservoir system is divided into two types of media based on the concept of dual-porosity model: matrix and fracture system. The gas flow model in each type of medium is built based on its distinctive flow characteristics. The existence of massive nanopores in shale reservoirs leads to a large amount of gas adsorbed on the pore walls in the shale matrix. The amount of adsorbed gas can vary from 35 to 58% in Barnett shale to 60–85% in Lewis shale [22]. In reality, when the desorbed gas combined with free gas flows in the shale reservoir, it is difficult to differentiate them and quantify their contributions to the total gas production in the wellhead separately. However, it is crucial to understand the dynamic production characteristics of adsorbed gas and free gas in order to design the reasonable strategy to develop shale gas. To explicitly characterize the contribution of these two types of shale gas, we adopted the “component subdivision” method proposed by Yang et al. [23], where the free gas and adsorbed gas are numerically treated as two different components and they can be marked separately but modeled simultaneously as a unified flow system. This method is based on the compositional reservoir simulation model. The basic principle to build the flow model is the mass conservation for each component in different medium.
In the shale matrix, for a specific component , the mass conservation law can be written aswhere and are the control volume and area, is the porosity, is the saturations of gas, are the mole densities of oil and gas, is the mole fraction of gas, is the mole flux of gas, and is mole flux of the adsorbed component entering the matrix pore, and it can be computed aswhere is the matrix density, is the time interval, is the initial adsorption concentration during this time interval, and is the equilibrium adsorption concentration. can be determined based on the Langmuir isothermal adsorption curve:where is the Langmuir volume representing the maximum adsorption volume of the shale matrix and is the Langmuir pressure representing the pressure at which 50% of the Langmuir volume can be adsorbed.
In the fracture system, for a specific component , the mass conservation equation can be written as
One main difference of flow mechanisms between matrix and fracture system is that there is no need to take the adsorption effect in the fracture system because the fractures essentially act as the main flow path, not a storage space as the matrix does. In addition, it is required to take the intersections between the fractures and wellbore into consideration to accurately characterize the fracture flow, and is the mole flux to describe this effect. The flow between matrix and fractures can be considered as the fluid exchange between two different control volumes and can be computed as the flow term . The flux between two adjacent cells can be estimated via Darcy’s law:where is relative permeability, is viscosity, and is the potential difference between cells and . The transmissibility is calculated by volume-weighted average.where and are the volumes of cells and , while and are the apparent permeabilities of cells and , respectively.
In dual-porosity-dual-permeability (DPDP) model, the transmissibility between a matrix cell and its corresponding fracture cell is defined aswhere is the grid volume. The fracture spacing along the , , and directions is represented by , , and and , , and are the matrix permeability of the three directions, respectively. , , and are dimensionless parameters related to the geometric information of fracture, which are equal to for sugar cube model.
The flow between hydraulic fracture and well is calculated by the following equation:where wi is well index which represents the conductivity between well and the hydraulic fracture and is the potential difference between well and cell . In an orthogonal grid, one option is to use Peaceman’s formula to compute well index.
If the well is parallel to direction,
If the well is parallel to direction,
If the well is parallel to direction,where is well radius. , , and are the projection lengths of the perforation interval onto the -axis,-axis, and -axis; is the skin factor; , , and are the size of the fracture grid block.
As mentioned before, the Knudsen diffusion needs to be considered into the gas transport in shale matrix. This effect can be described by adding a correction term to the apparent permeability under the same pressure gradient:where and are the apparent permeability and the intrinsic permeability, respectively, is the reservoir pressure, and is the correction coefficient of Knudsen diffusion. The correction coefficient can be determined through fitting the experiment data.
Non-Darcy flow may occur in high permeable area, such as hydraulic fractures. In reservoir simulation, the phenomenon of non-Darcy flow is characterized by Forchheimer equation:where is the unit converge unit converter and is Forchheimer coefficient. Forchheimer coefficient also can be measured through experiment data fitting.
2.2. Component Subdivision Method
In order to quantify the distributions of free and adsorbed gas in different regions, “component subdivision” method is used. The free and adsorbed gas from SRV and non-SRV region are treated as different components with the same PVT features in (1), namely, SRV free gas, SRV adsorbed gas, non-SRV free gas, and non-SRV adsorbed gas, respectively. The PVT features are calculated by flash calculation according to PR Equation of State in this study. According to Langmuir equation (3), one component’s adsorbed content is in equilibrium with its partial pressure in the vapor phase. In the component subdivision method used here, a specific type of shale gas (free or adsorbed gas) produced from a specific region (SRV or non-SRV) is treated as an independent “component”; it should be in equilibrium with the total methane’s partial pressure rather than its own partial pressure. The mole fraction of methane gas can be defined aswhere , , , and represent the mole fractions of free gas from SRV region, adsorbed gas from SRV region, free gas from non-SRV region, and adsorbed gas from non-SRV region. Langmuir equation (3) in component subdivided method is modified as
In the beginning of the shale gas production, only free gas flows into the wellbore and no adsorbed gas has been produced yet; thus the mole fractions of adsorbed components, and , need to be set as zero.
2.3. Numerical Model of the Shale Gas Reservoir
To implement the mathematical model developed above, we built a numerical model based on dual-porosity-dual-permeability (DPDP) concepts to investigate the effects of all parameters on the well performance. The reservoir properties data used in this numerical model were obtained from the field measurement data in a shale gas reservoir located in Sichuan basin. These productive marine shale formations are mainly lower Silurian Longmaxi and Ordovician Wufeng with the burial depth of 2,500 meters and formation temperature of 85°C. Almost all the wells are hydraulically fractured horizontal wells with multiple transverse fractures. Here, we built a single-well model by using the simulation code UNCONG developed by Li et al. [24] as shown in Figure 1, and the key shale reservoir properties are listed in Table 1 for the base model used in the following analyses. All the main formation properties used in this study are collected from well logging and lab experiments in the Fuling shale gas play. It is assumed that the length of the horizontal well is 1000 m with 30 hydraulic fracturing stages. The half-length of the hydraulic fractures is 150 m. The hydraulic fractures are depicted as black lines in Figure 1, which are characterized by using local grid refinement method in the numerical models. Due to the natural fracture reactivation after the hydraulic fracturing, there would be a secondary fracture network around the primary hydraulic fractures. Here, we used the spacing between two fracture layers to represent the complexity or the developed level of the secondary fracture network as done in Warren and Root [11]. The less the spacing is, the more complex or developed the secondary fracture network is. In Figure 1, the secondary fracture network is depicted in red color, and it also indicates the coverage area of the SRV region. The non-SRV region as depicted in Figure 1 in blue color has not experienced any stimulation and thus holds the properties of the shale matrix. According to the geological and well logging data, the original-gas-in-place (OGIP) is 4.76 × 108 m3 and the adsorbed gas takes up 40% of the OGIP, and the OGIP in the SRV region and non-SRV region is, respectively, 2.25 × 108 m3 and 2.52 × 108 m3.
We first validated our simulator UNCONG by comparing the simulation results with those obtained through commercial simulator CMG on the basis of the base case. Figure 2 shows the production rate and total production of gas from UNCONG and CMG. As can be observed in the Figure, the results of production rate obtained from UNCONG are very similar to those obtained from CMG. There are slight differences in the comparison results of total production, especially in the late stage of the total production curve, due to the accumulation effect of the total production with time, but the accuracy of the UNCONG simulator can be accepted. The main reason why we used UNCONG but not CMG or other commercial simulators is that they do not have the ability to quantify the contribution of different types of gas produced from different regions. This function has been developed in UNCONG and can aid us to implement the sensitivity analysis in this work.
3. Result Analysis
3.1. Contributions of the Adsorbed and Free Gas to the Production
The adsorbed gas and free gas production behaviors are investigated explicitly and separately by using the “component subdivision” method. In the base model, the horizontal well is producing shale gas at the prescribed production rate of 60,000 m3, and the minimum bottomhole pressure is set as 6 MPa, which is consistent with the working condition in Fuling gas field. Figure 3 depicts the characteristics of adsorbed, free, and total shale gas production in terms of both production rate and cumulative production. It indicates that, under this working schedule, this well can produce shale gas stably in the first 4 years and then the production starts to decline. The cumulative production after 30 years can be as much as 1.81 × 109 m3.
To clearly show the contributions of the adsorbed and free gas on the total production, the entire production life of this well is divided into three stages based on its production characteristics: the first stage starts from the beginning to the end of the stable production; the second stage starts from the beginning of the production decline and ends after 10 years of production; and the third stage covers the period when the shale gas production enters its late stable production.
As shown in Figure 4, in the first stage, the shale gas can be produced at relatively stable rate of roughly 6 × 104 m3, and the total gas production is dominated by free gas in the SRV region. With decreasing the reservoir pressure, the free gas production starts decreasing and the gas used to be adsorbed on the matrix begins to desorb and the adsorbed gas production climbs up. At the end of the first stage, the production rate of the adsorbed gas approaches to the maximum at the rate of 6.5 × 103 m3. Our finding is consistent with Wang [25] who states that the existence of adsorption gas has a more obvious effect on the production decline in the early stage and this is contrary to the common belief that the adsorption gas becomes important only when the average pressure in the reservoir drops to a certain level. In the second stage, the production rate starts declining. Due to the relatively easy depletion of the free gas, most of the free gas has been produced in the first stage, which leads the contribution of the adsorbed gas to become more obvious in the second stage. In the third stage, the production rate remains at a low level. The produced gas is dominated by the free gas in the non-SRV region. At the end of the 30 years, 73.76% of the total production rate is composed of the free gas in the non-SRV region.
Figure 5 indicates that, during the 30-year production period, the cumulative gas production is composed of the free gas in the SRV region with 67.57%, free gas in the non-SRV region with 20.81%, the adsorbed gas in the SRV region with 10.22%, and the adsorbed gas in the non-SRV region with 1.40%. The free gas still greatly dominates the overall shale gas production, even though the adsorbed gas takes up 40% of the OGIP.
3.2. The Effect of the Matrix Permeability on the Recovery Factor
After analyzing the production dynamics in the above section, we focused on the effect of the postfracturing reservoir properties on the recovery factor in different regions, that is, the ratio of recovered gas at a specific production time and original-gas-in-place in the SRV and non-SRV regions. For the purpose of result presentation, we divided the 30-year production period into three 10-year intervals so that the dynamic change of the recovery factor can be better studied.
To investigate the effect of matrix permeability, we selected 3 typical values in the range of the field observed matrix permeability in the shale gas reservoir: 1 × 10−4 mD, 1 × 10−5 mD, and 1 × 10−6 mD, respectively. Figure 6 shows the effect of the matrix permeability on the recovery factor in the SRV region. The result implies that the matrix permeability can severely affect the recovery factor in the SRV region, and the recovery factor in the SRV region decreases with decreasing the matrix permeability. Along the time line, the recovery factor decreases dramatically, which is understandable due to the rapid production decline in the ultra-low-permeable reservoir. In the non-SRV region, the effective recovery is largely determined by the matrix permeability, since the hydraulic fracturing stimulation has no effect on this part of the reservoir. Figure 7 shows the effect of the matrix permeability on the recovery factor in the non-SRV. In general, it shows the same tendency as that in the SRV region, that is, the lower the matrix permeability is, the smaller the recovery factor is. However, one noticeable point in this analysis is that the shale gas in the non-SRV region almost cannot be recovered at all when the matrix permeability is lower than 1 × 10−6 mD.
3.3. The Effect of the Secondary Fracture Network Properties on the Recovery Factor
In this section, the effect of the properties of the secondary fracture network on the recovery factor in different region of shale gas reservoir is investigated. Two main properties, the density and the permeability of the secondary fracture network, are selected here.
As mentioned before, the density of the fracture network can be described by using the spacing of the fracture layer in the dual-porosity model. Here, we used 5 different values of fracture layer spacing, 10, 15, 20, 25, and 30, respectively. Figure 8(a) shows the effect of the secondary fracture network density on the recovery factor in SRV region. It can be found that the effect of secondary fracture network density on the recovery factor in the SRV region is not obvious for the base case, where matrix permeability is 1 × 10−4 mD. For comparison purpose, a new simulation based on the matrix permeability of 1 × 10−5 mD is performed. When the matrix permeability is decreased, the effect of secondary fracture network density on the recovery factor becomes more profound. It indicates that it is not helpful to create a more complex secondary fracture network to improve the recovery factor in the SRV region if the matrix permeability in this region is already relatively high. Figure 8(b) shows the effect of the secondary fracture network density on the recovery factor in the non-SRV region. It can be found that the complex fracture network is always helpful to improve the recovery factor in the non-SRV region, since it increases the contact area with the non-SRV region.
(a)
(b)
To investigate the fracture network permeability on the recovery factor in different regions, 4 typical values are selected in this analysis. They are 0.001 mD, 0.005 mD, 0.01 mD, and 0.1 mD, respectively. Figure 9 shows the analysis results. It indicates a similar tendency to the effect of the fracture network complexity. Improving the permeability of the fracture network has a more profound effect on the non-SRV region than on the SRV region. If the hydraulic fracturing cannot create the primary fractures but only a secondary fracture network, then it is obvious that the secondary fracture network permeability can have a severe effect on the recovery factor in both SRV and non-SRV regions, as shown in Figure 10. However, it also can be observed that the recovery factor can be very small when the secondary fracture network permeability is low, which indicates that the creation of primary hydraulic fractures can alleviate the requirement on the secondary fracture network permeability. Therefore, in practice, it would be more favorable to create primary hydraulic fractures with moderate complex and permeable secondary fracture network than to just create a more uniform fracture network.
(a)
(b)
(a)
(b)
3.4. The Effect of the Primary Hydraulic Fracture Properties on the Recovery Factor
Three main properties of the hydraulic fractures on the recovery factor are investigated in this section. They are the number of hydraulic fracturing stages, fracture conductivity, and the half-length of the hydraulic fractures, respectively. For the stages number, simulation results based on three cases are compared, 10, 20, and 30. Figure 11 shows the effect of the stage numbers on the recovery factors in SRV and non-SRV regions. It can be found that the recovery factors in both SRV and non-SRV regions are not influenced when the number of the hydraulic fracturing stages is decreased from 30 in the base case to 20. However, when further reducing the stage numbers to 10, the recovery factors are severely affected. To search for the reason behind the phenomenon, the pressure distribution with the different stage numbers is plotted in Figure 12 (the pressure distribution with stage number of 20 is very similar to that with 30 and thus it is not plotted). It can be observed that the pressure drops uniformly along the horizontal well (it forms a large blue rectangle) when the stage number is 30, but the pressure drawdown when the stage number is 10 forms several sections (small blue rectangular) and each section is separated relative to the others. Recall that our horizontal well is 1000 m, and each hydraulic fracturing stage is associated with secondary fracture network with the width of 30 m on one side. The spacing between two hydraulic fractures is about 33.3 m if the stage number is 30. Therefore, the pressure system associated with each primary hydraulic fracture is well connected. However, the spacing between two hydraulic fractures is about 100 m if the stage number is 10. It means that there would be a 40 m wide space which is not stimulated by hydraulic fracturing process, and thus the pressure system associated with each primary hydraulic fracture is independent of the others. In the design of the hydraulic fracturing, it requires increasing the stage number to create a fracture system (including primary hydraulic fractures and the reactivated natural fractures) that can form a well-connected SRV region, and it is not worthy of further increasing the stage number when this requirement can be satisfied.
(a)
(b)
(a)
(b)
Four typical values of hydraulic fracture conductivity are selected to investigate the effect on shale gas production. They are 0.1, 0.5, 1, and 10 mDm, respectively. Figure 13 shows the effect of the hydraulic fracture conductivity on the recovery factor in different regions. It can be seen that in general higher hydraulic fracture conductivity is favorable to improve recovery factors in both SRV and non-SRV regions.
(a)
(b)
The half-length of the hydraulic fractures has an effect on the SRV area and thus on the estimation of the OGIP in SRV and non-SRV regions. Here, three cases of designed hydraulic fracture half-length are selected, that is, 150 m, 130 m, and 110 m., and the corresponding OGIPs in SRV region are 2.25 × 109 m3, 1.95 × 109 m3, and 1.65 × 109 m3. Figure 14 shows the effect of the half-length of the hydraulic fractures on the recovery factors in SRV and non-SRV regions. It can be observed that the half-length only has a moderate effect on the recovery factor in SRV region. But the effect of half-length on the recovery factor in non-SRV region is profound. With the increase of the half-length, the recovery factor in the non-SRV region increases.
(a)
(b)
4. Conclusion
This work built a numerical model to analyze the gas production characteristics in shale gas reservoir and used the field data collected from a shale gas play in Sichuan Basin in China to investigate the influence of various factors on the production characteristics. It leads to the following key findings:(1)The contribution of adsorbed gas and free gas from different regions to the total gas production was explicitly and separately studied by using a “component subdivision method.” Free gas in the SRV region dominates the initial gas production. The adsorbed gas production increases as the pressure of the shale reservoir gradually decreases and approaches the maximum value at the end of the initial stable production stage. When gas production declines, the production of adsorbed gas can alleviate the decline level by compensating the sharp decline of free gas production in the SRV region. When the gas production comes into the late stable production stage, the free gas in the non-SRV region dominates the gas production. In terms of the cumulative production, free gas is the main component in the total gas production.(2)The shale matrix permeability in the SRV region has a profound effect on recovery factor in the SRV region as matrix permeability in the non-SRV region does to recovery factor in the non-SRV region. Higher matrix permeability has a positive effect on the recovery factor.(3)The density and the permeability of the secondary natural fracture network have more obvious effect on the recovery factor in the non-SRV region than that in SRV region. The existence of the primary hydraulic fractures is necessary even if a well-connected secondary natural fracture network can be created by hydraulic fracturing.(4)The stage number of the hydraulic fracturing needs to be sufficiently large to create a unified pressure system along the horizontal well. The conductivity and half-length of the primary hydraulic fractures have positive effects on the recovery factor in both SRV and non-SRV regions.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is funded by the National Natural Science Foundation of China (Grant no. 41402199), National Science and Technology Major Project (Grants nos. 2016ZX05037003 and 2016ZX05060002), the Science Foundation of China University of Petroleum, Beijing (Grant no. 2462014YJRC038), China Postdoctoral Science Foundation (Grant no. 2016M591353), and the Science Foundation of Sinopec Group (Grant no. P16058).