An equivalent method to assess the production performance of horizontal wells with complicated hydraulic fracture network in shale oil reservoirs

Journal of Natural Gas Science and Engineering 71 (2019) 102975

Contents lists available at ScienceDirect

Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse

An equivalent method to assess the production performance of horizontal wells with complicated hydraulic fracture network in shale oil reservoirs

T

Yan Donga, Peichao Lib,∗∗, Wei Tianc, Yuxi Xiana, Detang Lua,∗ a

Department of Modern Mechanics, University of Science and Technology of China, Hefei, 230027, China School of Mechanical and Automotive Engineering, Shanghai University of Engineering Science, Shanghai, 201620, China c Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen, 518055, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Complex fracture network Equivalent well model Production performance XFEM

The production performance of horizontal well with complex fracture networks remains a challenging issue. To more accurately evaluate the stimulation effect of fracture network on the reservoirs, an equivalent method is presented based on the extended finite element method (XFEM) and the Peaceman well model. We discretize the governing equations with finite volume method and solve the linearized equations with fully implicit scheme. The reliability of the equivalent method is validated against the available numerical method in the literature. After that, two cases of different fracture network structures are investigated. Results show that the equivalent method is feasible to assess the effect of fracture network on horizontal well production. And more broadly distributed fracture network has better enhancing effect. Meanwhile, fracture tips have little stimulation effect on the production. The findings are of benefit to provide us a simplified way to investigate the flow in reservoirs containing complex fracture network.

1. Introduction Unconventional shale oil resource is of broad prospect for its massive reserve and the lower environmental impact. The technology of multiple fractured horizontal wells (MFHWs), which may develop complex fracture network, is widely applied to enhance the producing efficiency in shale reservoirs (Shen et al., 2017; Jiao et al., 2018). With the developing of hydraulic fracturing techniques, forecasting the production performance of horizontal wells with complex fracture network and optimizing treatment design have become significant issues. Based on the bi-wing fracture model, many researchers have developed analytical models of MFHWs in shale reservoirs. With the hypothesis of infinite-conductivity fractures, Wang (2014) developed a semi-analytical well testing model of MFHWs in consideration of flow diffusion, viscosity, desorption, and stress-sensitivity of natural fracture permeability. Ozkan et al. (2011) introduced a trilinear-flow model considering linear flow performance of fractures to assess the performance of MFHWs and to guide the hydraulic fracturing design of shale reservoirs. Ren and Guo (2015) presented a finite-conductivity MFHWs model to investigate the influence of different angles of hydraulic fractures and fracture asymmetry on the performance of MFHWs. Zhao

∗

et al. (2016) presented a mathematical gas flow model in nano-pores to describe a horizontal well with multiple fractures taking into account the slippage and diffusion effect. Compared to analytical methods, the advantage of numerical methods is the capacity of handling more complex conditions (Shen et al., 2016). Most of the numerical models adopt constant well-pressure condition to analyze the cumulative production of reservoir. Yu et al. (2014) used the simulator CMG-IMEX to study the effect of biwing fractures half-length on the production performance. Zhao et al. (2017) established a mathematical model of an arbitrarily off-centered, infinite-conductivity fractured well in gas reservoir and implemented the boundary element method (BEM) to study the well transient pressure response. Except for bi-wing fracture, the dendritic-like fracture geometry was presented based on micro-seismic mapping technique. Xu et al. (2017a) investigated the influence of dendritic-like fracture shape on the horizontal wells production performance using discrete fracture method (DFM) based on unstructured Voronoi grids (PEBI). DFM believes that properties of fractures are independent and each fracture has individual contribution to the well performance. AlTwaijri et al. (2018) investigated the influence of complex fracture network on the production performance of shale gas horizontal well based on the embedded discrete fracture method (EDFM). However, this model ignored the

Corresponding author. Corresponding author. E-mail addresses: [email protected] (P. Li), [email protected] (D. Lu).

∗∗

https://doi.org/10.1016/j.jngse.2019.102975 Received 16 May 2019; Received in revised form 28 July 2019; Accepted 17 August 2019 Available online 22 August 2019 1875-5100/ © 2019 Published by Elsevier B.V.

Journal of Natural Gas Science and Engineering 71 (2019) 102975

Y. Dong, et al.

interaction between hydraulic fractures and natural fractures, which is discrepant with the actual situation. Zhang et al. (2019b) introduced the construction of hydraulic fractures on the basis of micro-seismic monitoring data, and then proposed a production prediction model of MFHWs with complex fracture network combining with multiple interacting continua theory and DFM. Actually in shale reservoirs, hydraulic fractures are not just straight cracks. Natural fractures, which may affect the hydraulic fracture propagation, are often observed within the reservoir formation (Gale et al., 2007). Numerous researchers believed that it’s possible to generate highly complex fracture network during the hydraulic fracturing process. Some researchers established the propagation model of hydraulic fracture network based on the discrete element method (DEM) (Zou et al., 2016; Zhang et al., 2017, 2019a) and the displacement discontinuity method (DDM) (Ren et al., 2016; Zeng and Yao, 2016; Kresse and Weng, 2018). Wang et al. (2016) and Zeng et al. (2018b) investigated the propagation of induced fracture in orthotropic formations and found that the induced fracture tends to change to the direction of material axis with larger modulus based on the XFEM. Zeng et al. (2018a) used the embedded discrete fracture model combining with the XFEM to examine the flow-deformation coupling effect upon the propagation of multiple fractures. However, lots of production evaluating models have been presented without considering the interaction between hydraulic fractures and natural fractures. Considering that the branch fractures of complex fracture network are related to each other, it's still difficult to accurately describe the fluid flow and the flow distribution in fracture network (Wu et al., 2014). Some methods only carry out the stimulation of the fracture network by increasing the permeability of a certain matrix area, while ignoring the influence of the specific fracture structure. Meanwhile, as for the commercial software, the production performance prediction for discrete fracture network model remains a challenging issue because of the complicated structure of the discrete fracture network. That is, the production performance evaluation method still needs further study. Analytical well test method usually involves well model. Peaceman (1983) presented a commonly applied well model when the well trajectory locates in the center of grid. The main purposes of this work are to present an equivalent method combining the shape of complex fracture network with the Peaceman well model and more accurately investigate the effect of arbitrary crack network on production performance in unconventional shale oil reservoirs.

Fig. 1. Schematic illustration of fracture propagation model.

vector, and ai and biα are the enriched degrees of freedom vector of different enriched elements. For crack-crossed elements, H(x) is Heaviside step function

1 , dist (x ) > 0 H (x ) = ⎧ ⎨ ⎩− 1 , dist (x ) < 0

where dist(x) is the shortest distance from x to crack. In order to determine it's positive or negative, n is defined as the unit normal vector to the crack. The distance is positive if x locates in the side where n points, while it is negative if x locates in the other side (Zhuang and Cheng, 2011). For embedded elements, φα(r, θ) could be denoted as

φα (r , θ) =

In this section, we present the fracture network model based on the XFEM, including deformation governing equation, interaction process between hydraulic fractures and natural fractures, and ultimately the construction of fracture network.

∫Ω σ: ε(∗) dΩ=∫Ω ⋅∗dΩ+∫Γ ⋅∗dΩ+∫Γ p⋅w∗ds t

(3)

c

(4)

where σ is the stress of the domain, ε is the strain vector, b is the body force, u* is the virtual displacement, and w is crack aperture. By substituting the displacement interpolating function Equation (1) into Equation (4), we obtain the discrete equilibrium equation as

2.1. Deformation governing equation The discontinuity of fracture is characterized by the enriched shape functions in XFEM. That allows the XFEM to overcome the limits intrinsic in meshing and to simulate fracture with arbitrary pathways. In 2D space, the node displacement approximation is as follows (Ren et al., 2009):

⎡ ⎤ 4 H (x ) ai + ∑ φα biα ⎥ ∑ Ni (x ) ⎢ui +   ⎢ α = i∈N 1⎥ i ∈ NH  i ∈ Nemb ⎦ ⎣

θ θ θ θ r ⎡sin , cos , sin sin θ , cos sin θ ⎤ 2 2 2 2 ⎣ ⎦

where r and θ are polar coordinates of crack-tip. As shown in Fig. 1, consider an isotropic domain Ω with a boundary Γ, which consists of Γu, Γt, and Γc. The prescribed displacement u is constrained on Γu, traction t is imposed on Γτ. Γc is fracture plane which is distinguished by Γc+ and Γc−, and n+ and n− are the normal vectors of Γc+ and Γc−. Water pressure p is imposed on Γc in the direction of n+ and n−. By virtual work principle, the weak form of equilibrium equation is obtained as

2. Fracture network model

uh (x ) =

(2)

K ⋅U = F

(5)

where K, U, and F are global stiffness matrix, global displacement vector, and external force vector, respectively. They are obtained by assembling the element stiffness matrix Kije, element displacement vector Ue, and element external force vector fie. Here, Ue, Kije, and fie are defined as:

U e = {ui, aj , bk1, bk2 , bk3, bk4} (i ∈ N , j ∈ NH , k ∈ Nemb) (1)

where N is the set of all nodes, NH is the set of nodes in crack-crossed elements, and Nemb is the set of nodes in embedded elements. Ni (x) is the standard shape function, ui is the continuous degrees of freedom

Kije

2

uu ua ub ⎡ K ij K ij K ij ⎤ ⎢ au ⎥ aa = ⎢ K ij K ij K ijab ⎥ ⎢ K bu K ba K bb ⎥ ij ij ⎦ ⎣ ij

(6)

(7)

Journal of Natural Gas Science and Engineering 71 (2019) 102975

Y. Dong, et al.

f ie = {f iu , f ia , f ib1 , f ib2 , f ib3 , f ib4 }

T

determined to obtain the route of the induced fracture in next time step.

(8)

2.3. Construction of fracture network 2.2. Interaction process Fig. 3(a) is the schematic of MFHW. The green region represents the horizontal well which is perforated for multiple clusters, and the blue solid lines represent natural fractures. The horizontal well is parallel to the direction of minimum horizontal in-situ stress σhmin, and the induced fractures initial direction is perpendicular to the wellbore. Assuming that initial induced fractures are far enough away to ignore the stress interference effect, we select a typical research region marked with red dashed line. The complex fracture model is shown in Fig. 3(b). The left boundary is fixed in the x direction, and other boundaries are free. The black solid line represents the initial induced fracture. According to the theories mentioned above, we can obtain the ultimate fracture network shape as shown in Fig. 4. Arbitrary azimuth angles and the locations of natural fractures result in the formation of various intersection patterns during the hydraulic fracturing process, which makes the complex hydraulic network eventually appear. Crack aperture is related to the normal stress across the fracture. It can be described as

We adopt the maximum circumferential stress criterion as the fracture propagation criterion. The generation of fracture network will involve the intersection between hydraulic fracture and natural fracture. When a hydraulic fracture intersects a natural fracture with the approaching angle of α, the crossing criterion and opening criterion are described as follows:

⎧ σT > T0 ⎨ ⎩ p > σn + Tcoh

crossing criterion opening criterion

(9)

where σT is the normal stress component parallel to natural fracture, σn is the normal stress of natural fracture, T0 and Tcoh are the tensile strength and cohesive strength of natural fracture respectively. If the natural fracture is opened, the induced fracture may re-initiate on the opposite. Blanton (1986) presented the criterion of hydraulic fracture re-initiation as below:

σh max − σh min −1 > T0 cos 2 α − b sin 2 α

(10)

w=

Here, b is the offset coefficient, which is introduced in detail in Blanton (1986). The interaction process is illustrated in Fig. 2. First, we focus on the opening criterion. The induced fracture will open the natural fracture if the opening criterion is satisfied. And then, the stress on the re-initiation point x0 is calculated to identify whether the fracture will re-initiate based on Equation (10). If re-initiation occurs, the induced fracture will offset. If not, simultaneity-mode will appear when the crossing required pressure is equal to this open pressure. If the opening criterion is dissatisfied, then the crossing criterion will determine whether the induced fracture can cross the natural fracture. If not, the hydraulic fracture will be arrested by natural fracture. That is, there are five patterns of interaction between natural fracture and hydraulic fracture when the induced fracture intersects the natural fracture, namely penetration, arrest, opening, offset, and simultaneity-mode. Consequently, the intersection pattern can be

2(p − σn ) hf (1 − ν 2) E

(11)

where w is the crack aperture, p is the water pressure, and hf is the fracture height. Considering water pressure p is constant, the bigger normal stress means the smaller aperture. This indicates that crack aperture is varying in different fracture segments. Thus, we need to take into account that the different fracture widths have different stimulation effects on the formation. 3. Equivalent method In this section, the flow governing equation is introduced first. Then, we present an equivalent well model to simplify the flow in the shale reservoir containing complex fracture network. To develop the equivalent method for production performance evaluation, the following assumptions are made: (1) Single phase oil flows through isothermal, incompressible and single layer reservoir, and the gravity effect is ignored. (2) The fluid (oil) is weakly compressible, and the migration of fluid from the matrix to the fracture conforms to the Darcy's law. (3) Before the production stage, the reservoir is under steady condition. The initial pore pressure of the reservoir is assumed to be equal in the production stage. That is, the impact of fracturing on formation pore pressure is not considered. (4) Enlightened by analytical well test method, this equivalent method is based on constant well flux condition. 3.1. Flow equation The flow equation in matrix is given as follows:

K ∂ φ ⎛ ⎞+q ∇⋅⎡ (∇p − ρg∇Z ) ⎤ = ⎢ ⎥ ∂t ⎝ B ⎠ ⎣ μB ⎦

(12)

where p is the pore pressure, K is the matrix permeability, μ is the fluid viscosity, φ is the matrix porosity, B is the fluid volume factor, and q is the source/sink term. Hydraulic fractures enhance the horizontal well productivity by creating complex fracture network system with higher conductivity. An accurate well model is necessary to quantify the enhancement in numerical simulation process. Peaceman (1983) presented a commonly used well model when the well trajectory locates in the center of grid. The flux of the well is given as follows:

Fig. 2. Schematic of intersection criterion. 3

Journal of Natural Gas Science and Engineering 71 (2019) 102975

Y. Dong, et al.

Fig. 3. Geometry of multiple fractured horizontal well.

Fig. 5. Schematic of equivalent scheme (a) actual element crossed by crack; (b) equivalent element.

normal direction of the crack plane. Considering the aperture of same crack segment has little change, the lengths of two crack planes are approximately equal, which are denoted by li. The flow rate of crack segment in element i is

Qfi = 2li h

Fig. 4. Ultimate fracture network shape.

q=

2πKh (p − pwf ) μB ln(re/ rw ) i

K ∂p μ ∂n

(14)

Fig. 5(b) shows that the equivalent well locates on element center. According to Darcy law for radial flow, the flow rate of equivalent well in element i is

(13)

where q is the production of this well, h is the formation thickness, re is the equivalent radius of a well grid, rw is the wellbore radius, and B is fluid formation volume factor. pi is the numerically calculated average pressure for the well grid. pwf is the bottom hole pressure. re is defined as that radius at which the steady-state flowing pressure for this well is equal to pi.

qwi =

2πK e hrw ∂p ∂r μ

= Qfi r = rw

(15)

where Ke is equivalent permeability. The flow rate of the crack segment and the equivalent well should be equal. Based on the assumption that the normal pressure gradient is equal to the radial pressure gradient, we can obtain the equivalent permeability is

3.2. Equivalent well model

Ke = As shown in Fig. 5, the fracture segment in an element is transformed into an equivalent well model located in the grid center. The following assumptions are made:

li K rw

(16)

By substituting the above equivalent permeability into the Peaceman well model, the equivalent well model can be expressed as

2πK e h e (p − pwf ) μB ln(re/ rw ) i

(1) The equivalent wellbore flow rate/flux is assigned proportionately according to the length and aperture of crack segment. (2) Steady flow is considered near the crack. (3) The element is sufficiently small, and the normal pressure gradient in Fig. 5(a) is equal to the radial pressure gradient in Fig. 5 (b).

qwi =

Fig. 5 (a) shows the element crossed by a single crack whose aperture is wi. The solid blue lines represent the crack planes. According to the Darcy law, fluid migrates from matrix to the crack along the

where Bref is reference formation volume factor, Cf is compressibility coefficient and pref is reference pressure. Considering well produces at constant flux, the flow rate is assigned

where

B=

4

pewf

(17)

is the equivalent wellbore pressure. B is defined as

Bref [1 + Cf (p − pref )]

(18)

Journal of Natural Gas Science and Engineering 71 (2019) 102975

Y. Dong, et al.

Fig. 6. Three typical element types containing multi-crack segments.

equations

Table 1 Basic parameters of verification model.

∫Ω f dΩ = Ax − B = 0

Parameter

Value

Unit

Initial reservoir pressure Matrix porosity Formation thickness Matrix permeability Equivalent wellbore radius Production Reservoir dimensions Fracture length Production time Coefficient of oil compressibility

35 0.1 30 1 0.001 14 100 × 100 50 30 0.001

MPa fraction m md m m3/day m m day 1/MPa

where f is the residual of the equivalent well equations and flow equations in block Ω. And x = [Δp]T. Coupling the change of equivalent wellbore pressure, the matrix form of Equation (21) is:

⎡ Amm Amw ⎤ ⎡ Δpm ⎤ = ⎡ Bm ⎤ ⎢ Awm Aww ⎥ ⎢ Δpw ⎥ ⎢ B w ⎥ ⎣ ⎦⎣ ⎦ ⎣ ⎦

(22)

where m is the number of all elements, w is the number of crack-crossed elements. The linearized equation system is solved with newton iteration method through the matrix solver GMRES (Saad and Schultz, 1986).

proportionately according to the length and aperture of crack segment, which is expressed as below

qwi = Qfi = Q

(21)

4.2. Verification

li wi ∑j ∈ Eenr l j wj

The advantage of the presented equivalent method is that it can simplify the flow between matrix and fracture and has the ability to conduct the flow of complex fracture network. In order to verify the presented method, a numerical case with single straight fracture is studied using a fully implicit numerical simulation method based on PEBI grid (Li et al., 2016). Basic parameters are listed in Table 1. Based on XFEM, we simulate and obtain a hydraulic fracture (solid red line) with length of 50 m as can be seen in Fig. 7. Then we investigate the changes of formation pore pressure along two lines as illustrated in Fig. 7. Validation results are represented in Fig. 8 and Fig. 9. A good agreement of pore pressure along Line 1 and Line 2 is observed between our model and Li's model from Figs. 8 and 9. The difference between our results and Li's results along Line 2 is less than

(19)

where Q is the total production, and wi is the average aperture of crack segment in element i. Eenr is the set of all elements crossed by the fracture network. In practice, the generation of complex hydraulic fracture network inevitably involve fracture intersections (Xu et al., 2017b). At this point, there may be multiple fracture segments in one element, which need to be considered separately. There are three typical branching situations: two separate cracks, a junction, and two intersecting cracks, as shown in Fig. 6. Because of the stress shadow effect, the crack aperture exhibits a much sharper drop at the intersection point. Cracks are divided into several segments in these cases, and the aperture of different segments is unequal. In case (a), there are two segments: AB and CD. And in case (b), there are 3 segments: AO, BO and CO. Meanwhile, there are 4 segments in case (c): AO, BO, CO and DO. Thus, the flow rates of these elements are expressed as 2

−

⎧ q i = Qi = Q ∑k = 1 lk wk− , f ⎪ w ∑j ∈ E l j wj enr ⎪ − 3 ⎪ i lk wk ∑ qw = Qfi = Q k = 1 − , ∑j ∈ E l j wj ⎨ enr − ⎪ ∑k4 = 1 lk wk i i ⎪q = Q = Q −, f ⎪ w l w ∑j ∈ E enr j j ⎩

in case(a) in case(b) in case(c) (20)

4. Numerical solution and verification 4.1. Numerical solution The partial differential equations is solved with finite volume method. And the fully implicit scheme and upstream weighted are applied to discrete the equations. After the linearization of the nonlinear terms, the volume integral can be expressed as a system of linear

Fig. 7. Schematic of the boundary straight fracture model. 5

Journal of Natural Gas Science and Engineering 71 (2019) 102975

Y. Dong, et al.

Fig. 8. Pore pressure along Line 1. Fig. 11. Schematic of initial fractures distribution.

0.1 MPa when 0 < x < 60, and the relative error is less than 1%. However, the difference increases when x > 60, which is caused by the distinction of fracture aperture in the two models. In our fracture network model, the aperture of straight fracture is approximately equal except crack tip region. This causes the pore pressure to drop more on the right side. Nevertheless, the biggest difference of two model is less than 0.2 MPa and the relative error is less than 2%. In addition, we compare the results of our model and Li's model at three different coefficients of oil compressibility. The changes of formation pore pressure along line 1 are represented in Fig. 10. As we can see, as the coefficients of oil compressibility increases, the changes of formation pore pressure decrease. And the relative error under the three conditions is less than 2%. Consequently, the accuracy of our model can be validated. 5. Results and discussion Fig. 9. Pore pressure along Line 2.

In this section, two examples are presented to analyze the effect of fracture network on production performance. The model is 100m × 100m × 30 m in size. The left boundary is fixed in the x direction, and other boundaries are free. The block contains one hydraulic fracture and a natural fracture system with arbitrary azimuths as shown in Fig. 11. The black solid line between red points represents the initial hydraulic fracture, and the blue solid lines represent natural fractures. A linear elastic and isotropic block is considered with Young's modulus of 30 GPa and Poisson's ratio of 0.22. Initial in-situ stress is distributed as σhmax = 4 MPa and σhmin = 3 MPa, and the tensile strength T0 and the cohesive strength Tcoh are 1.5 MPa and 0.1 MPa. Water pressures in the hydraulic fracture is 10 MPa. With seeding 200 nodes each side, the model is divided into 39601 elements. Then, the shape of fracture network can be simulated and obtained based on XFEM. Fig. 12 shows the construction of fracture network in different natural fracture distribution conditions, as marked with (a) and (b). As we can see in Fig. 12 (a), the fracture network contains many crack branches, and it extends a long way in the x direction. Meanwhile, the fracture network in Fig. 12 (b) tends to expand in the y direction, and several branches are connected at the initial stage of fracturing. All the crack branches in Fig. 12 (b) are in the left region of the formation. In comparison, it can be stated that the fracture network in Fig. 12 (a) is more broadly distributed. It follows from Fig. 11 that the locations and azimuths of natural fractures are important factors affecting the structure of fracture network. The hydraulic fracture in case (b) intersects many natural fractures in the early stage and several

Fig. 10. Pore pressure along Line 1 for three different coefficients of oil compressibility.

6

Journal of Natural Gas Science and Engineering 71 (2019) 102975

Y. Dong, et al.

Fig. 12. Shape of different fracture networks.

branches interconnect. Consequently, the azimuths of natural fractures in case (b) enable the crack branches to be inclined to propagate along the y direction. Suppose the production is 12.96 m3/day, and the other physical parameters are the same as the validation example. Before the production stage, the reservoir is regarded as steady condition. The initial pore pressure throughout the reservoir at the production stage is assumed to be 35 MPa.The pore pressure field of case (a) and (b) are visually depicted in Figs. 13 and 14, respectively. It can be seen that pressure drawdown takes place in the shape of fracture network outline. Comparing with the validation example, much larger stimulated reservoir volume (SRV) coverage is generated. In other words, fracture network practically stimulates low permeability formations with higher conductivity. Another interesting phenomenon from Figs. 13 and 14 is that the ends of cracks have little stimulation effect on the production capacity. That is, the fracture tips basically have no influence on the pressure field. The reason is that the aperture of crack tip is much smaller than that of anterior fractures. In general, the pressure drawdown phenomenon is more obvious near the injection point. In other words, the induced cracks generated in the early stage have greater Fig. 14. Pressure distribution of case (b) after 30 days of production.

stimulation effect on the reservoir, which is consistent with the Weng et al. (2014)'s results qualitatively. Considering the fracture network in case (b) is mainly in the left region, we can see that the SRV of fracture network in case (a) is larger. Furthermore, the biggest formation pressure difference of case (a) is 0.91 MPa while that of case (b) is 1.34 MPa. It is noteworthy that the biggest formation pressure difference of single straight fracture is 2.37 MPa. This indicates that fracture network can enhance the production capacity. And the stimulation effect will be enhanced if the fracture network is more broadly distributed as shown in Fig. 12 (a). 6. Conclusions In this paper, we present an equivalent method combining the XFEM with analytical well test method to assess the production performance of horizontal wells with hydraulic fracture network. The key feature of this method is to transform the fracture segments into an equivalent well model at the center of grid. It can provide a simplified way to more accurately investigate the effect of arbitrary crack network on production performance in unconventional shale reservoirs without

Fig. 13. Pressure distribution of case (a) after 30 days of production. 7

Journal of Natural Gas Science and Engineering 71 (2019) 102975

Y. Dong, et al.

the need for re-meshing. The validation of this new method is performed by comparison with a fully implicit numerical simulation based on PEBI grid. The influence of fracture network on the pore pressure is investigated. The main conclusions are summarized as follows:

permeability reservoir with pseudo threshold pressure gradient. J. Pet. Sci. Eng. 147, 308–316. Ozkan, E., Brown, M.L., Raghavan, R., Kazemi, H., 2011. Comparison of fractured-horizontal-well performance in tight sand and shale reservoirs. SPE Reserv. Eval. Eng. 14 (02), 248–259. Peaceman, D.W., 1983. Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability. SPE J. 23 (03), 531–543. Ren, J., Guo, P., 2015. A novel semi-analytical model for finite-conductivity multiple fractured horizontal wells in shale gas reservoirs. J. Nat. Gas Sci. Eng. 24, 35–51. Ren, L., Su, Y., Zhan, S., Hao, Y., Meng, F., Sheng, G., 2016. Modeling and simulation of complex fracture network propagation with SRV fracturing in unconventional shale reservoirs. J. Nat. Gas Sci. Eng. 28, 132–141. Ren, Q., Dong, Y., Yu, T., 2009. Numerical modeling of concrete hydraulic fracturing with extended finite element method. Sci. China Ser. E Technol. Sci. 52 (3), 559–565. Saad, Y., Schultz, M.H., 1986. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 7 (3), 856–869. Shen, W., Li, X., Xu, Y., Sun, Y., Huang, W., 2017. Gas flow behavior of nanoscale pores in shale gas reservoirs. Energies 10 (6), 751. Shen, W., Xu, Y., Li, X., Huang, W., Gu, J., 2016. Numerical simulation of gas and water flow mechanism in hydraulically fractured shale gas reservoirs. J. Nat. Gas Sci. Eng. 35, 726–735. Wang, H., 2014. Performance of multiple fractured horizontal wells in shale gas reservoirs with consideration of multiple mechanisms. J. Hydrol. 510, 299–312. Wang, X., Shi, F., Liu, H., Wu, H., 2016. Numerical simulation of hydraulic fracturing in orthotropic formation based on the extended finite element method. J. Nat. Gas Sci. Eng. 33, 56–69. Weng, X., Kresse, O., Chuprakov, D., Cohen, C.-E., Prioul, R., Ganguly, U., 2014. Applying complex fracture model and integrated workflow in unconventional reservoirs. J. Pet. Sci. Eng. 124, 468–483. Wu, Y., Li, J., Ding, D., Wang, C., Di, Y., 2014. A generalized framework model for the simulation of gas production in unconventional gas reservoirs. SPE J. 19 (05), 845–857. Xu, C., Li, P., Lu, D., 2017a. Production performance of horizontal wells with dendriticlike hydraulic fractures in tight gas reservoirs. J. Pet. Sci. Eng. 148, 64–72. Xu, D., Liu, Z., Zhuang, Z., Zeng, Q., Wang, T., 2017b. Study on interaction between induced and natural fractures by extended finite element method. Sci. China Phys. Mech. Astron. 60 (2), 024611. Yu, W., Luo, Z., Javadpour, F., Varavei, A., Sepehrnoori, K., 2014. Sensitivity analysis of hydraulic fracture geometry in shale gas reservoirs. J. Pet. Sci. Eng. 113, 1–7. Zeng, Q., Liu, W., Yao, J., 2018a. Hydro-mechanical modeling of hydraulic fracture propagation based on embedded discrete fracture model and extended finite element method. J. Pet. Sci. Eng. 167, 64–77. Zeng, Q., Yao, J., 2016. Numerical simulation of fracture network generation in naturally fractured reservoirs. J. Nat. Gas Sci. Eng. 30, 430–443. Zeng, Q., Yao, J., Shao, J., 2018b. Numerical study of hydraulic fracture propagation accounting for rock anisotropy. J. Pet. Sci. Eng. 160, 422–432. Zhang, F., Damjanac, B., Maxwell, S., 2019a. Investigating hydraulic fracturing complexity in naturally fractured rock masses using fully coupled multiscale numerical modeling. Rock Mech. Rock Eng. 1–24. Zhang, F., Dontsov, E., Mack, M., 2017. Fully coupled simulation of a hydraulic fracture interacting with natural fractures with a hybrid discrete‐continuum method. Int. J. Numer. Anal. Methods Geomech. 41 (13), 1430–1452. Zhang, R., Zhang, L., Tang, H., Chen, S., Zhao, Y., Wu, J., Wang, K., 2019b. A simulator for production prediction of multistage fractured horizontal well in shale gas reservoir considering complex fracture geometry. J. Nat. Gas Sci. Eng. 67, 14–29. Zhao, Y., Li, H., Zhang, L., Kang, B., 2017. Pressure transient analysis for off-centered fractured vertical wells in arbitrarily shaped gas reservoirs with the BEM. J. Pet. Sci. Eng. 156, 167–180. Zhao, Y., Zhang, L., Xiong, Y., Zhou, Y., Liu, Q., Chen, D., 2016. Pressure response and production performance for multi-fractured horizontal wells with complex seepage mechanism in box-shaped shale gas reservoir. J. Nat. Gas Sci. Eng. 32, 66–80. Zhuang, Z., Cheng, B.-B., 2011. Development of X-FEM methodology and study on mixedmode crack propagation. Acta Mech. Sin. 27 (3), 406–415. Zou, Y., Zhang, S., Ma, X., Zhou, T., Zeng, B., 2016. Numerical investigation of hydraulic fracture network propagation in naturally fractured shale formations. J. Struct. Geol. 84, 1–13.

(1) The equivalent method simplifies the flow in the formation containing the fracture network, and is feasible to assess the production performance of horizontal wells with hydraulic fracture network. (2) Fracture network can enhance the production capacity. And the enhancing effect will be better if the fracture network is more broadly distributed. (3) Crack aperture is one of the significant influencing factors for the stimulation effect. The apertures in crack tips are so small that the ends of cracks have little stimulation effect on the production capacity. Note that, our equivalent model adopts constant flux condition to assign the flow rate by the proportion of the length and crack segment aperture. That is, it is not applicable to the situation of constant wellpressure. Besides, as previously mentioned, the effect of fracturing process on the reservoir pore pressure is neglected. So, we would like to consider the case of constant well-pressure condition and the influence of fracturing process in future works. Acknowledgements This work is supported by the National Science and Technology Major Project of China (Grant No. 2017ZX05009005-002), the Natural Science Foundation of Shanghai, China (No. 19ZR1421400), and the China Postdoctoral Science Foundation (Grant No. 2018M643232). The authors are also indebted to the reviewers and the editors for their constructive comments on the manuscript. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.jngse.2019.102975. References AlTwaijri, M., Xia, Z., Yu, W., Qu, L., Hu, Y., Xu, Y., Sepehrnoori, K., 2018. Numerical study of complex fracture geometry effect on two-phase performance of shale-gas wells using the fast EDFM method. J. Pet. Sci. Eng. 164, 603–622. Blanton, T.L., 1986. Propagation of hydraulically and dynamically induced fractures in naturally fractured reservoirs. In: SPE 15261 Presented at Unconventional Gas Technology Symposium, Louisville, Kentucky, pp. 18–21 (May). Gale, J.F., Reed, R.M., Holder, J., 2007. Natural fractures in the Barnett Shale and their importance for hydraulic fracture treatments. AAPG (Am. Assoc. Pet. Geol.) Bull. 91 (4), 603–622. Jiao, C., Hu, Y., Xu, X., Lu, X., Shen, W., Hu, X., 2018. Study on the effects of fracture on permeability with pore-fracture network model. Energy Explor. Exploit. 36 (6), 1556–1565. Kresse, O., Weng, X., 2018. Numerical modeling of 3D hydraulic fractures interaction in complex naturally fractured formations. Rock Mech. Rock Eng. 1–19. Li, D., Zha, W., Liu, S., Wang, L., Lu, D., 2016. Pressure transient analysis of low

8