Three-dimensional representation of discrete fracture matrix model for fractured reservoirs

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol

Three-dimensional representation of discrete fracture matrix model for fractured reservoirs

T

Yuyun Zhaoa,b, Hanqiao Jianga, Sheik Rahmanb, Yudong Yuanb, Lin Zhaoa,c, Jiawei Lib, Jiachao Geb, Junjian Lia,∗ a

State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing, China School of Minerals and Energy Resources Engineering, University of New South Wales, Sydney, Australia c Mathematics and Cybernetics, SINTEF Digital, Oslo, Norway b

ARTICLE INFO

ABSTRACT

Keywords: 3D fracture modeling Discrete fracture matrix model Field scale Production analysis

Accurate modeling of fractures is a matter of growing scientific interest. Several existing methods for 3D fractures modeling are available only for oversimplified geometries, e.g. rectangle. The identification of intersections of complex polygonal fractures with dip angles has not been clearly resolved, which results in inaccurate grid transmissibility calculation between fracture cells. This study aims at obtaining an efficient modeling strategy for 3D fractures with arbitrary spatial distributions and field-scale fracture modeling application. The dimensionality reduction method based on coordinate transformation is employed to identify the intersections of arbitrary shaped polygon fractures with any dip angles. Meanwhile, fine grid cells are generated in the vicinity of fractures to adapt high-pressure gradients, which is favorable to the simulation accuracy. For another, the local refinement allows transitions of cell sizes from the near-fracture region to far-fracture region, reducing grid number and further computational cost. The model validation is given in a tailored case where Cartesian grids are used to model vertical fractures with a local grid refinement by the industry-reference simulator ECLIPSE-E100 and a good match is achieved by our model. Three complex examples with different fracture densities are computed to demonstrate the robustness and stability of the simulator. The influence of 3D fracture parameters on the oil production, including the height, dip angle, area and shape, is analyzed.

1. Introduction Conventional oil and gas resources cannot satisfy the rising global demand, which stimulated the exploitation of unconventional reservoirs. Unconventional oil and gas, such as tight gas, shale gas, coalbed methane and tight oil, show a great potential under current economic and technological condition (Jiang et al., 2016; Zhao et al., 2018). Most unconventional oil and gas resources rely on the creation of induced fractures to be productive. The production of oil from fractured reservoirs has increased drastically worldwide recently, due to the development of advanced technologies such as horizontal drilling and multi-staged hydraulic fracturing (Abdelazim and Rahman, 2016). Therefore, the numerical simulation of fractured reservoirs has also attracted extensive attention and has been significantly developed. Among numerous models, the well-known dual-porosity (DP) model is the earliest model to deal with fractured reservoirs. This model was presented by Warren and Root (Warren et al., 1963) to model evenly distributed small-scale fractures with good connectivity. There are

∗

many developments based on DP: the dual-porosity dual-permeability (DPDP) model (Gilman et al., 1986), the multiple interacting continua (MINC) model (Pruess et al., 1985), the multi-porosity (MP) model (Yan et al., 2015), etc. However, these models are not suitable for sparse or poorly connected fracture networks. More precisely representing the spatial position and geometry of the fracture, discrete fracture matrix (DFM) model can provide better accuracy. A DFM model represents the fracture as a line segment in 2D space in the grid domain and treats it as a rectangle with the width equal to the aperture of the fracture in the computational domain (Kim and Deo, 2000; Juanes et al., 2002). DFM models usually use more flexible PEBI or tetrahedral (triangular in 2D) grids so that it can handle fractures of more complex patterns, such as arc fractures (Karimi-Fard et al., 2003). This approach has been presented in several studies (Li et al., 2006; Hoteit and Firoozabadi, 2008; Mallison et al., 2010; Jiang et al., 2016; Fang et al., 2018; Zhao et al., 2018). In the embedded discrete fracture matrix (EDFM) model, the fracture meshing is separate and that constitutes its main advantage, i.e. no need to tackle the

Corresponding author. E-mail address: [email protected] (J. Li).

https://doi.org/10.1016/j.petrol.2019.06.015 Received 21 November 2018; Received in revised form 3 June 2019; Accepted 4 June 2019 Available online 10 June 2019 0920-4105/ © 2019 Elsevier B.V. All rights reserved.

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

difficult problem of creating a high-quality, conforming unstructured grid (Lee et al., 2001; Li et al., 2006; Yan et al., 2016). Jiang et al. provided hybrid coupled DFM and multi-continuum models, in which multi-continuum model is used to simulate small scale fractures while DFM model is used to handle large scale fractures (Jiang et al., 2016). Discrete fracture network (DFN) model can also be used for heterogeneous reservoir simulation. It separates the grids into three types: fracture segment, fracture intersection (node) and rock matrix. However, the matrix grids are meshed roughly and approximated by rectangle with equivalent volume, and the flow in matrix is ignored (Boyle et al., 2012; McKoy and Sams, 1997). The enhanced discrete fracture network (EDFN) model takes the flow in matrix into account and presents a better accuracy (Mi et al., 2017). Due to the complexity of real fractured reservoirs, the 2D simplified models cannot accurately simulate 3D systems. In recent decades, some have studied the numerical simulation of three-dimensional fractured reservoirs. Wei et al. (1998) introduced a 3D numerical model to simulate pressure transient through fracture and matrix based on the DP model. Azim (2016) presented a 3D DFM model to evaluate the water coning in naturally fractured reservoirs. Azim and Sheik (Abdelazim and Rahman, 2016) obtained the 3D fracture map based on the pressure derivatives using this model. Moinfar (2013) developed a 3D EDFM model for simulating fluid flow in naturally fractured reservoirs. Fumagalli (Fumagalli et al., 2016) realized the automatic identification of the intersections of the bilinear surfaces (rectangular or parallelogram fractures), and provided an upscaling procedure for fractured reservoirs with the 3D EDFM model. Mallison and Hui made a lot of outstanding contributions in the field of construction and field application of the 3D DFM model with unstructured grids. They first introduced a realistic modeling method for fracture networks in a giant carbonate reservoir (Hui et al., 2007). Then they built a new grid generation algorithm for discrete fracture modeling workflows (Mallison et al., 2010), providing a simpler refinement scheme. This meshing method was used to construct a 2.5D DFM model for flow simulation in a giant carbonate reservoir (Hui et al., 2008), and for three-phase, fully compositional simulations of IOR/EOR processes (Hui et al., 2013b), and other flow simulation application (Hui et al., 2013a). Lim used this 2.5D DFM model to realistic field-scale models in practical simulation studies (Lim et al. (2009). They also provided a comprehensive methodology for gridding, discretizing, coarsening, and simulating 3D DFM models of naturally fractured reservoirs (Hui et al., 2018). However, Wei's model cannot simulate heterogeneous reservoirs. Azim and Moinfar's models cannot identify the intersections of fractures. Their treatments of the intersection of fractures are ambiguous, and the local grid refinement is not realized in their models. Fumagalli's model and Mallison and Hui's model limited to handle the rectangular fractures. Mallison and Hui's 2.5D DFM model did not consider the dip angle of the fracture in the vertical direction. However, true 3D fractures are composed of rectangular and elliptical fractures or other complex polygonal fractures. Even for simple rectangular fractures,

Fig. 2. Complex intersection of 3D fractures.

Fig. 3. Coordinate transformation.

complex polygonal shapes (e.g. trapezoid and triangle) may be caused by the influence of reservoir boundaries (Fig. 7). To our best knowledge, few papers have totally solved the problem of identifying the intersection of 3D fractures of arbitrary geometry shapes and any dip angles. And few papers have studied the effect of 3D parameters of fracture, such as dip angle, geometry and so on, on the oil productivity. In this paper, we develop a true 3D DFM model to simulate the elliptical and polygonal fractures in field-scale 3D fractured reservoirs. The dimensionality reduction method based on coordinate transformation is used to completely solve the automatic identification problem of the intersection of large number of arbitrarily shaped polygon fractures with any dip angles. Based on the mixed programming of MATLAB and C++, a faster simulator with a better compatibility is written to apply the 3D-DFM model to the realistic field scale. Our simulator is validated by a good match of the simulation result with the commercial numerical simulator ECLIPSE-E100 in the case of simple vertical fractures. Three complex examples with different fracture densities are calculated to demonstrate that our model is robust and practical for industry applications. Then the effect of 3D parameters of fracture, including the height, dip angle, area and geometry, on oil productivity is analyzed. The specific details of the meshing, which are foundation of our formulation, are shown in appendices.

Fig. 1. Geometrical description of 3D fractures. 887

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

Fig. 4. Identification of fracture intersections (black lines).

Fig. 5. Meshing of fracture constrained by intersecting line segments based on the dimensionality reduction method.

Fig. 6. Meshing of reservoir constrained by fractures.

Fig. 7. Complex polygonal fractures caused by the boundaries and meshing.

2. Theory

statistical regularity, which can be inferred from analysis of natural outcrops and core samples (McKoy and Sams, 1997). The vertical aspects of the fractures complicate the characterization process, which introduces the necessity of additional parameters including the dip angle and the height of the vertical direction. In this paper, we mainly consider the elliptical and rectangular fractures with central axes parallel to the horizontal plane. This feature enables the characterization by conventional parameters such as orientation and dip angle. Fig. 1 presents an illustration of two different fractures in three views. The length of the fracture is exactly that of central axis, denoted as L1 and L2 . Meanwhile, the angle between the central axis and the north direction in the clockwise direction is the

2.1. Geometrical description of 3D fracture From a geomechanical perspective, a fracture is any discontinuity of a formation where cohesion loss appears; in other words, a result of a rupture (van Golf-Racht, 1982). In terms of geometry, a fracture in the 2D space can be simplified to a line segment with a finite length. The main parameters of the fracture include the center-point position, orientation, length and aperture, first three of which can be described by the coordinates of the two endpoints. In the stratum, natural fractures do not exist in isolation, but cluster together, and exhibit certain 888

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

Table 1 Parameters for EF1, RF3 and RF5.

Fig. 8. Transmissibility between matrix grids.

Parameter

Value

Parameter

Value

Reservoir Size Reservoir Depth Initial Pressure

240m× 180m× 20 m 2000m 30 MPa

0.6 1 mm

Reference Pressure Matrix Porosity

23.4 MPa 0.1

Fracture Porosity Fracture Aperture Fracture Compressibility

Matrix Permeability Matrix Compressibility

5 × 10

0.1mD

3

MPa−1

Oil Density Oil Formation Volume Factor Oil Viscosity Oil Compressibility

5 × 10 3 MPa−1 800 kg/m3 1.065 5.0mPa s

5 × 10 3 MPa−1

(Persson and Strang, 2004) and TetGen (Hang, 2015). Fracture surfaces can be input as constraints to them. However, these pieces of software are not specifically written for fractured reservoir flow simulation, and can not automatically identify the intersection of fractures (Fig. 6-b). So it is necessary to determine the intersection before meshing. This is a very important process not clearly shown in some articles about the 3DDFM. 2.2.1. Identification of fracture intersection In 2D space, the intersection of fractures is simpler, which is a point where two line segments intersect. In 3D space, the intersection of fractures is significantly more complicated. The intersecting lines of the two whole planes can be simply calculated with formulas. It is more complicated to directly determine the intersection of two three-dimensional fractures, because fractures are polygonal or elliptical with boundary, rather than a whole plane. They may not intersect even if their whole planes have an intersecting line. In addition, 3D fractures do not intersect although their horizontal central axes intersect, which is shown in Fig. 2. Therefore, we introduce a dimensionality reduction method by coordinate transformation employing a matrix (M), as shown in Fig. 3. The first step is to determine the intersection of the planes of two fractures. Then coordinate transformation is performed to shift the fracture surface from 3D to 2D. The intersecting segments of a polygonal or elliptical fracture and the projection of intersecting line are identified in this 2D space (specified in Appendix A and Appendix B). Then the endpoints of this segment in the original 3D space can be

Fig. 9. Transmissibility between matrix and fracture grids.

orientation 1 and 2 . From the side view, the fracture resembles a segment, the angle between which and the horizontal line is the dip angle 1 and 2 . The height of this segment is the height of the fracture H1 and H2 . Some use the height of vertical projection, H 1, to represent the height of the fracture (McKoy and Sams, 1997). However this definition is only justified when the dip angle is 90°. If the dip angle is 90°, on the top view, the top edge V3 V4 and the bottom edge V1 V2 will coincide with the central axis of the fracture, which degenerates into the familiar 2D case. So the 3D characterization method described here is a natural extension of the traditional 2D method. 2.2. Meshing of 3D fractured reservoir There are some commercial and open software packages that can generate tetrahedral grids with some constraints, such as DistMesh

Fig. 10. Transmissibility between fracture grids and Star-Delta transformation.

Fig. 11. Three models with simple vertical fractures. 889

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

Fig. 12. Comparison of reservoir pressures (Eclipse vs. 3D-DFM model).

fracture. Then, the method of distmesh proposed in the literature (Persson and Strang, 2004; Persson, 2005) is used to obtain the grids in 2D space. The resulting nodes of gridblocks are transformed into the original 3D space by inverse coordinate transformation (Fig. 5).

M

{f , s } 3D

Fig. 13. Comparison of fracture pressures in EF1 (Eclipse vs. 3D-DFM model).

{f , s } 2D

distmesh

{P , t } 2D

M

1

{P , t } 3D

(2)

2.2.3. Meshing of reservoir constrained by fractures Once the grids of the fracture surface are obtained, their nodes and edges can be used as constraints for reservoir meshing. Local refinement is achieved by retaining the node with a probability , which is inversely proportional to the minimum distance (r) from the node to all fracture surfaces. If this probability is greater than the generated random number between 0 and 1, then this node will be conserved, otherwise it will be deleted as shown in equation (3). The proportion between and 1/ r 3 is taken as r03 , where r0 is the minimum distance from all nodes to fracture surfaces.

rand <

= r03/ ri3, keep ith node,

rand

= r03/ ri3, delete ith node,

(3)

where, for ith node ni , ri is the minimum distance of ni to all fracture surfaces, i.e. f j for all j, which can be represented as follows. Fig. 14. Comparison of the oil production rate (Eclipse vs. 3D-DFM model).

3D

M

{f , l } 2D

identifying intersections

{f

l} 2D

M

1

{f

r0 = min ri.

(5)

i

Then initial grids can be obtained by the Delaunay algorithm (Edelsbrunner, 2001). To modify the gridblocks closer to regular tetrahedrons, an optimization method is used to move the positions of the remaining nodes. The accuracy is improved as the segment between centroids of two adjacent grids is perpendicular to their contact surface. The edges of the grids (the connections between pairs of nodes) can be treated as bars, and the nodes are joints of the truss. Each bar has a force dependent on its length. The force balance of each truss is regarded as the optimization goal. The iteration in optimization is completed when all nodes are no longer moved, i.e. changes in position are all less than a preset threshold. The effectiveness of this method has been verified to be reliable (Persson and Strang, 2004; Persson, 2005). This meshing process can also be implemented by the open software distmesh or tetgen. However, distmesh is based entirely on MATLAB and is slower, not suitable for large-scale meshing at the field scale.

l} 3D

(4)

j

obtained by inverse coordinate transformation. Such intersecting line segments for both fractures are acquired, and consequently the common part of these two line segments is determined by judging the positional relationship of their four endpoints. If the common part is absent, these two fractures do not intersect. The identification results are shown in Fig. 4.

{f , l}

ri = min dist(ni , f j ),

(1)

2.2.2. Meshing of fractures constrained by intersecting line segments An intersection line between fractures is a line segment with finite length. These segments on a fracture are used as constraints in fracture meshing. First, with coordinate transformation, the fracture and the intersecting line segments are transformed onto the plane of the 890

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

Fig. 15. Three complex examples with different fracture densities to show the robustness.

Fig. 16. A realistic field-scale model with both natural (blackish green) and hydraulic (blue) fractures. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

2.3. Flow equations and discretization

Table 2 Parameters of natural fractures for the realistic field-scale model.

Set 1

Set 2

Parameter

distribution

Min/Mean

Max/Var

Center point Length Orientation Dip angle Center point Length Orientation Dip angle

Uniform Uniform Normal Uniform Uniform Lognormal Normal Uniform

(20 m, 20 m) 60 m N6°E 70° (20 m, 20 m) 80 m N135°E 70°

(1680 m, 980 m) 300 m 8° 90° (1680 m, 980 m) 20 m 8° 90°

Considering the single-phase fluid, by the law of conservation of mass, for any control volume , the governing equation is

(

) t

d +

v n dS =

sq

d ,

(6)

where is the boundary of control volume , n is the unit normal vector of dS , ρ is the fluid density subsurface and s is the fluid density on the surface (kg/m3), ϕ is the rock porosity, v is the Darcy velocity, q is the source/sink term (kg/s). For slightly compressible flow in a porous medium, both fluid density ρ and rock porosity ϕ are related to pressure p, and the relationships can be characterized by linear growth models as follows:

Tetgen is written in C++ and is faster and powerful. But it is not written specifically for DFM model. For penetrating fractures (i.e., fractures intersect with the top or bottom of the reservoir), the problem of self-intersection of the input facets needs to be addressed first. The partition results are shown in Fig. 6-a and Fig. 7. There are three types of cells in 2D-DFM mesh: 2D matrix cell, 1D fracture cell and 0D fracture intersection cell, respectively. Similarly, there are also three types of gridblocks in 3D-DFM grid, including 3D matrix cell, 2D fracture cell and 1D fracture intersection cell.

=

0 (1

+ cf (p

p0 )),

=

0 (1

+ cR (p

p0 )),

(7)

where p0 is the reference pressure, cf is the fluid compressibility, and cR is the rock compressibility, 0 , 0 are the rock porosity and fluid density at the reference pressure p0 , respectively. Assuming that the fluid in the reservoir follows Darcy's law, then 891

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

Fig. 17. Reservoir and fracture pressures of the realistic field-scale model shown in Fig. 16.

q=

PIi i (pi i

pw ),

(10)

w

where w is the set of well grids, PIi is the productivity index of the ith grid. So equation (6) can be discretized in the following matrix form,

C T U T (Cp

K ( p µ

g z ),

v n dS

ij Tij (pi j Ni

pj ) +

ij Tij ¯ij g (z i j Ni

zj ),

1 (PV t

PV0 0),

(11)

2.4. Transmissibility evaluation

(8)

where K is the permeability tensor, and z is the depth difference, which cannot be ignored in 3D simulation. The second term of equation (6) can be discretized as follows. i

pw ) =

where U is the upwind matrix mapping from a given face to its upwind grid. , T and PI are the diagonal matrix representing the mobility, the transmissibility and the productivity index, respectively. C is a sparse matrix used to calculate the gradient, whose i-row Ni1-column element value is 1 and i-row Ni2 -column element value is 1, where N is a matrix mapping a given face to its adjacent grids. The automatic differentiation framework of MRST is used to solve equation (11) in this paper (Lie et al., 2012; Krogstad et al., 2015; Lie, 2014).

Fig. 18. Effect of fracture height on oil production rate.

v=

¯ gCz ) + PI (p

2.4.1. Transmissibility between matrix cells The two-point flux-approximation (TPFA) scheme is used in this paper to calculate the transmissibilities. For incompressible and slightly compressible fluids, the transmissibility between matrix grids is the harmonic average of the half transmissibility values from the centroid of the grid to the interface, ignoring viscosity changes (Karimi-Fard et al., 2003; Chen, 2007; Lie, 2014), defined as below.

(9)

where i represents the grid cell considered, Ni denotes the set of indices of all adjacent grid cells of i , ij is the interface mobility, Tij is the transmissibility between grid cells i and j , obtained by the two-point flux approximation (TPFA) method, and ¯ij is the average fluid density of grid cells i and j . The right hand side of equation (6), the well flow, can be written as

Tij =

Ti =

Ti Tj Ti + Tj

Aij di

2

,

Ki d i ni ,

(12)

Tj =

Aij dj

2

Kj dj nj .

where, Aij is the area of interface between grids

Fig. 19. Simulation results with different dip angles for Case 1 (when reservoir thickness is low). 892

(13) i

and

j,

di is the

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

Fig. 20. Simulation results with different dip angles for Case 2 (when reservoir thickness is large).

vector from centroid of grid i to centroid of the interface, ni is the outward unit normal vector of the interface as a part of i (see Fig. 8).

3. Results

The first is a single vertical elliptical fracture (EF1) with a 50 m halflength and a 12 m height, and the other parameters are shown in Table 1. The second case has a primary rectangular fracture intersecting with two secondary rectangular fractures (RF3). The half-length and height of the primary fracture are 50 m, 20 m, respectively, and those of the secondary fractures are 24 m, 20 m. There are five parallel rectangular fractures (RF5) in the third case, with a half-length of 26 m, 34 m, 42 m, 34 m, 26 m and height of 20 m. The single horizontal well only perforates at the location of fractures. The cell size of the structured grid is 4 m (length of the cube edge), and the scale of local refined grid is 0.8 m. The minimum size of the unstructured grid in 3D-DFM model is 2 m (length of the tetrahedron edge), while the maximum size is 10 m. We simulate primary depletion on these two models and compare the reservoir pressures and oil production rates for these three cases. The comparisons between 3D-DFM model simulation results and ECLIPSE-E100 are shown in Figs. 12–14, respectively. The results show a good match between 3D-DFM model and ECLIPSE. The first case has the biggest maximum relative error (MRE) of 7.7%, and a maximum absolute error (MAE) of 0.65m3/day. The error mainly comes from the approximation to elliptical fracture by structured grids. The maximum relative errors of the second and third are both within 5%. The following three complex examples with different fracture densities demonstrate that our model is robust and practical for industry applications (Fig. 15).

3.1. Model validation

3.2. A realistic field-scale model

We validate the 3D-DFM model by comparing against an equivalent LGR-structured grid model, which is modeled by the commercial numerical simulator ECLIPSE-E100. Three models with simple vertical fractures are considered as shown in Fig. 11, involving different fracture shapes and intersection scenarios.

A realistic field-scale example with both natural and hydraulic fractures is given below (Fig. 16). The reservoir dimensions are 1700 m × 1000m × 40m . There are 1057 natural fractures and 20 hydraulic fractures. The parameters of hydraulic fractures are obtained from the interpretation of microseismic data. The geometrical parameters of

2.4.2. Transmissibility between matrix and fracture cells In grid domain, the fracture is a plane without thickness, while it is a block with a thickness of aperture wf in computational domain, as shown in Fig. 9. So the transmissibility between matrix and fracture cells can also be obtained by formula (12) and (13). 2.4.3. Transmissibility between fracture cells As mentioned in the description of the meshing process above, the intersection of multi-fractures is a line segment with finite length and is divided into small line segments of length h 0 . The well-known StarDelta transformation (Karimi-Fard et al., 2003) shown in Fig. 10 can be used to deal with the fracture intersection. Then the transmissibility between fracture grids can be calculated as

T12 =

T10 T20 T10 + T20 + T30

T1 T2 , T1 + T2 + T3

T13 =

T10 T30 T10 + T20 + T30

T1 T3 , T1 + T2 + T3

T23 =

T20 T30 T10 + T20 + T30

T2 T3 , T1 + T2 + T3

(14)

where Ti is given in equation (13).

Fig. 21. The extent of pressure depletion zone on the 11th day for Case 1 (when reservoir thickness is low).

893

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

Fig. 22. The extent of pressure depletion zone on the 11th day for Case 2 (when reservoir thickness is large).

4. Discussion This paper mainly focuses on the effect of 3D parameters including height, dip angle and geometry on the oil productivity. 4.1. The effect of fracture height on productivity The oil production rate in a primary depletion scenario is positively correlated with fracture height. Because a higher fracture can create a greater stimulated rock volume in the reservoir. An overestimation of fracture height in 2D DFM gridding results in more highly-permeable regions in the reservoir than what exists in reality. Fig. 18 compares the oil production rates with different heights based on RF5. 4.2. The effect of fracture dip angle on productivity Fig. 23. The volume of gridblocks of the extent of pressure depletion zone on the 11th day.

The oil production rates are calculated for different dip angles in two cases with various reservoir thicknesses in Figs. 21 and 22. The reservoir thickness is 20 m in Case 1 and 80 m in Case 2, and other parameters are kept the same in both cases. The results show that if the reservoir is not thick enough, then the oil production rate increases as the dip angle decreases (Fig. 19). Otherwise, the oil production rate is independent on the dip angle (Fig. 20). The extent of pressure depletion zone on the 11th day is calculated, which refers to all gridblocks with pressure less than the initial pressure pint . There is also a small amount of pressure change in the cells where the pressure has not spread due to the error of the numerical simulation. This is the so-called numerical dispersion (Hoteit et al., 2016). Therefore, we use the modification pintm = pint to replace pint , where is a small value and is set as 0.01 MPa in this paper. The total volume of these gridblocks is calculated with changing dip angles, shown in Fig. 23. The result exhibits the total volume changes slightly if the reservoir thickness is large enough (Case 2). Otherwise the volume decreases with an increasing of the dip angle (Case 1). Figs. 21 and 22 show the 3D view and side view of these cells for Case 1 and Case 2 respectively. It can be seen that the range of pressure propagation does not reach the boundary of the reservoir in Case 2 for all dip angles. So the flow contribution from matrix to fracture is same, then the oil production rate curves coincide. The range of pressure propagation reaches the boundary of the reservoir in Case 1, and the smaller the dip angle is, the larger the remaining part is, so the flow

Table 3 Parameters and results for 7 cases of different fracture geometries.

Case Case Case Case Case Case Case

1 2 3 4 5 6 7

Shape

a

Rectangle Ellipse Rectangle Rectangle Hexagon Circle Hexagon

160.00 50.93 80.00 40.00 24.82 22.57 20.00

a/m

b

b/m

10.00 10.00 20.00 40.00

Perimeter/m

SF

CP/m3

340.00 226.55 200.00 160.00 148.92 141.81 120.00

0.17 0.39 0.50 0.79 0.91 1.00 0.91

3646.50 3563.90 3525.30 3480.80 3478.10 3471.90 2291.10

a a is the length of the rectangle, the semi-major axis of the ellipse, the edge length of the hexagon or the radius of the circle. b b is the width of the rectangle or the semi-minor axis of the ellipse.

natural fractures obey some statistic laws shown in Table 2 which are obtained by natural outcrop data. The apertures of natural fractures are uniformly 50µm , while those of hydraulic fractures are 2 mm. Heights of fractures are all 40 m. The size of fracture cell is 5 m and the boundary cell is 25 m. There are a total of 3,675,933 cells. The reservoir and fracture pressures on the 511th day are shown in Fig. 17.

894

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

Fig. 24. Comparison of the intersection of fractures with different shapes.

contribution from the matrix to the fracture decreases as the dip angle increases. 4.3. The effect of fracture geometry on productivity In reservoir simulation, the fractures are most commonly represented as rectangles and ellipses. There are also articles that mention hexagonal fractures (Zihua et al., 2009). This paper studied the effect of fracture geometry on the oil productivity. We first studied the effect of a single fracture with different geometries and shape factors (SF) on the cumulative oil production (CP) under the same area (Case 1–6), and the effect of the fracture area on the CP under the same SF (Case 5 vs. Case 7). SF is defined by the following equation

SF =

4 A , P2

(15) Fig. 26. Comparison of the cumulative oil production for three cases shown in Fig. 25.

where A is the area, P is the perimeter. CP for 1000 days is simulated. The result in Table 3 shows that CP is positively correlated with the fracture area and negative correlated with the fracture SF. The effect of fracture area on CP is significant, while the effect of SF is very little. Differences in CP caused by square fracture, regular hexagonal fracture and circular fracture can be ignored in the case of disjoint. In fact, in the case where the fracture areas are equal, the smaller the SF of the fracture, the larger the perimeter and the more contact with the matrix, so more CP can be obtained. For square, regular hexagonal and circular fractures with the same area, their perimeters have

a small difference. Therefore, they cause a small difference in CP. If considering the intersection of fractures, in the case where the fractures have the same length and height, rectangular fractures are most likely to intersect others in 3D space, while hexagonal fractures and elliptical fractures are related to their shape factors. For instance, the rectangular fracture i in Fig. 24 intersects the three fractures 1, 2, and 3. The hexagonal fracture i intersects only fractures 1 and 3. And the elliptical fracture i intersects none of them. Both the rectangular

Fig. 25. Different representations of fracture geometry and the reservoir pressures on the 1111th day.

895

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

fracture ii and elliptical fracture ii intersect fractures 4 and 5. And the hexagonal fracture ii intersects none of them. A model of 274 natural fractures with unequal heights is considered. Rectangle, hexagon and ellipse are used to represent the geometry of fractures, respectively (Fig. 25). The number of intersecting segments is 1164, 762 and 864, and the total length is 14229 m, 7480 m and 8910 m, respectively. The result illustrated in Fig. 26 shows that the trend of the CP curve of the hexagonal fracture network is different from others because of the difference in fracture area and connectivity. Since parameters obtained from the field are always fracture length, height and so on, the case of equal areas is not discussed here. Obviously, in the case where the fracture area and height are equal, the elliptical fracture is longer, followed by the hexagonal fracture. This will affect the intersection of the fractures, and affect the oil production.

2. 3. 4.

5.

6.

5. Conclusion The main two contributions of this work are the development and the field-scale application of a complete 3D-DFM model for fractured reservoirs like shale oil and tight oil reservoir and the research of the effect of 3D parameters of fracture on oil productivity. The automatic identification problem of the intersection of large number of 3D arbitrarily shaped polygon fractures with any dip angles is completely solved by the dimensionality reduction method based on coordinate transformation. So that the conductivity between intersecting fracture cells can be accurately calculated by the ‘Y- ’ transformation method. The LGR is performed to capture near-fracture flow transients, as well as to reduce the number of global cells, both of which can provide a more accurate and efficient simulation. Based on the mixed programming of MATLAB and C++, a faster simulator with a better compatibility is written to apply the 3D-DFM model to the realistic field scale. The effect of 3D parameters of fracture on oil productivity for reservoirs with low matrix permeability can be concluded as follows.

contact area of the fracture with the matrix and the extent of pressure propagation are main factors affecting the oil production. The fracture area has a significant effect on the oil production. The larger the fracture area, the more the high permeable regions and the more contact with the matrix, the higher the oil production. The oil production rate is enhanced with a larger fracture height, as it creates larger fracture area. The effect of fracture dip angle on oil production rate is dependent on the fracture height and reservoir thickness. When the thickness of the reservoir is sufficiently larger than the fracture height, the dip angle has little effect on the oil production rate. Otherwise, the oil production rate is decreased with an increase in the dip angle. In the case of equal areas, the smaller the shape factor of the fracture, the larger the perimeter and the more contact with the matrix, so more oil production can be obtained. However, at the field scale, the effects of shape factors are largely negligible. The shape of the fracture affects their intersection, which affects the connectivity, then affects the oil production. Therefore, in the case of intersection, the influence of the shape of the 3D fracture on the oil production cannot be ignored.

We found that 3D fractures of any shape are difficult to characterize with the well-known parameters such as orientation, dip angle and so on. How to accurately represent fractures of arbitrary shape in 3D space is a complicated but valuable research topic. One of the future work will focus on the method using the centroid, topology, area and tensor related to the normal vector. Some problems in meshing, such as low grid quality caused by complex shapes like triangles, quadrilaterals, pentagons, and hexagons formed by several intersecting fractures, require further study. Acknowledgement This work is supported by the National Basic Research Program of China (973 Program, 2015CB250900). The authors gratefully acknowledge financial support from China Scholarship Council (No. 201706440045).

1. The greater stimulated rock volume provided by the fracture, the Appendix D. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.petrol.2019.06.015. Appendix A. Procedure of identifying fractures intersection Appendix A.1. Intersection of rectangular fractures In the specific meshing process, the coordinates of the four vertices, V1, V2 , V3 , V4 , are used as parameters for the rectangular fracture, and the position coordinates of the center point, C, the horizontal axis end point, A, and the incline axis end point, B, are used as parameters for the elliptical fracture. These parameters and the previous parameters can be converted to each other through a simple geometric relationship, and will not be described in detail herein. For a rectangular fracture, assume that its four vertices are Vi = (x i , yi , z i)T , k = 1,2,3,4 respectively, then its unit normal vector is

n = (n x , n y , n z )T =

V1 V2 × V1 V3 , V1 V2 × V1 V3

(A.1)

where, V1 Vk represents the vector from point V1 to point Vk , and

V1 Vk = (xk

x1, yk

y1 , z k

z1)T ,

(A.2)

k = 2,3.

So its plane equation is

n x (x

x1) + n y (y

y1) + nz (z

(A.3)

z1) = 0.

896

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

Fig. A.27. Dimensionality reduction method by coordinate transformation for rectangular fractures.

Let = n × , where = V1 V2/ V1 V2 , then a new 3D space x y z is formed based on vector groups { , , n} , which is shown in Figure. In x y z space, the fracture is a rectangle in a plane parallel to x y plane. Let matrix M denote the space-to-space coordinate transformation matrix from xyz space to x y z space, then

1 M = 0 , 0

0 M = 1 , 0

0 Mn = 0 , 1

(A.4)

So M [ , , n] = I , where I is the identity matrix. It's easy to know [ , , n] is an orthogonal matrix, so

M = [ , , n]

1

(A.5)

= [ , , n]T .

Any point V in x y z space can be back to xyz space by inverse transformation M

1

= [ , , n], that is, (A.6)

V = M 1V . At the same time, since M is an orthogonal matrix, the coordinate transformation is distance guaranteed, that is,

= M

TM T M

=

T

=

=

(A.7)

.

The coordinate transformation method used in this paper can transform the 3D case into a 2D case we are familiar with, which greatly reduces the mesh difficulty. Under the action of the matrix M, the positional relationship between the rectangle and the line in the 3D space can be transformed into the relationship between the rectangle and the line on the 2D plane. As shown in Figure A.27, assume that the rectangle Vi1 Vi2 Vi3 Vi 4 will be transformed to rectangle V i1 V i2 V i3 V i 4 by coordinate transformation, which means

V

ik

= MVik ,

(A.8)

k = 1,2,3,4.

Notice that

V

V

ik

i, k + 1

Vi, k + 1) = [ , , n]T (Vik

= M (Vik

Vi, k + 1),

(A.9)

k = 1,2,3,4.

Vi, k + 1) = 0 , that is, the z component of Vik Vi, k + 1 is 0, which means z components of where ni is the unit normal vector of fracture fi , so V i1 , V i2 , V i3 , V i 4 are equal, denote it as z i . The following of this paper uses V i1 , V i2 , V i3 , V i 4 represent the x y components of V i1 , V i2 , V i3 , V i 4 respectively without confusion. It is easy to know that the intersecting line lij of the planes where the fractures fi are located respectively is the solution of following equation niT (Vik

nix (x

x i1) + niy (y

yi1) + niz (z

z i1) = 0,

njx (x

xj1) + njy (y

yj1 ) + njz (z

z j1) = 0.

(A.10)

Denote the projection of this line in the plane of the fracture fi as l ij , and its equation as

Li (V ) = 0,

or

(A.11)

Li (x , y ) = 0.

It can be proved by derivation that Li has the following formula (see Appendix B)

Li (x , y ) = nTj i x + nTj i y + nTj ni niT Vi1 For vertices V

Li (V ik ) Li (V then V

ik

i, k + 1 ,

(A.12)

if (A.13)

i, k + 1) > 0,

and V

Li (V ik ) Li (V

and V

ik

nTj Vj1.

i, k + 1

i, k + 1)

must be on the same side of the line l ij , that's to say, line segment V

ik V i, k + 1

and line l ij do not intersect; if (A.14)

< 0,

then V ik and V i, k + 1 must be on both sides of the line l ij , that's to say, line segment V the line which V ik and V i, k + 1 are located and the line l ij ; if

ik V i, k + 1

and line l ij intersect, which is exactly the intersection of (A.15)

Li (V ik ) = 0, 897

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

then the vertex V ik is the intersection. There are three cases of the number of intersections between the straight line l ij and the rectangle V i1 , V i2 , V i3 , V i 4 , which are 0, 1, and 2. Only when the number of intersection points is 2, the intersecting line segment between the fractures can be formed. So the other two cases are not considered. Assuming that the two intersections are P i1 and P i2 in x y z space, the coordinates Pi1 and Pi2 in xyz space can be gotten by the inverse coordinate transformation, i.e.

Pik = Mi

1

P ik , zi

k = 1,2.

(A.16)

In the same way, the intersections Pj1 and Pj2 of the fracture f j and the straight line lij can be obtained. The positional relationship between Pi1, Pi2 and Pj1, Pj2 has the following five cases shown in Figure A.28 (if the order of i and j is not considered, otherwise, there are 9 cases, changing the order of i and j except for (5)). For the first case, where there are no common parts in the two intersecting line segments, the fractures fi and f j do not intersect. For the second case, the fractures fi and f j intersect at the line segment Pj1 Pj2 . For the latter three cases (3–5), they intersect at the line segment Pj1 Pi2 .

Fig. A.28. The positional relationship between Pi1, Pi2 and Pj1, Pj2 .2

Appendix A.2. Intersection of elliptical fractures As shown in Figure A.29, similarly to before, let = CA/ CA , = CB / CB , n = × , then a new 3D space x y z is formed based on vector groups { , , n} , which is shown in Figure. In x y z space, the fracture is an ellipse in a plane parallel to x y plane.

Fig. A.29. Dimensionality reduction method by coordinate transformation for elliptical fractures.3

After the coordinate transformation, we only need to consider the intersection of the ellipse and the straight line in the 2D plane, and it is easy to know that the intersection point is the solution of the following equation.

(x

C i1 ) 2 (y + a2

C i2 )2 = 1, Li (x , y ) = 0, b2

(A.17)

After obtaining the coordinates of the intersection point on the 2D plane, similar to the previous one, it is changed back to the space by inverse transformation, and the positional relationship of the four intersection points can be further discussed. Appendix A.3. Intersection of elliptical fractures and rectangular fractures In this case, the first step is to calculate the plane equation of the rectangular and elliptical fractures, then calculate the intersecting straight line of these two planes. The second step is to project this intersecting line onto the 2D fracture plane by the coordinate transformation, and calculate the common part of the ellipse and the projected line on the elliptical fracture plane. The same steps are repeated for the rectangular fracture. Then the two line segments are respectively transformed back to the original 3D space by the inverse coordinate transformations, respectively. At last, the common part of the two line segments is determined in the 3D space, resulting in the intersection of these two fractures.

898

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al.

Appendix B. Derivation of the straight line projection formula of intersecting of fractures In 3D xyz space, the equation of intersecting line lij between planes can be written in matrix form as follows.

x y z

niT nTj

niT 0T

0T Vj1 = 0. nTj

Vi1

(B.1)

It's easy to know that the equation of intersecting line l ij in x y z space has the following form. (B.2)

Ax + By + C = 0.

According to the coordinate transformation formula, for any point (x , y ) on the line l ij (recall that the z component is z i in x y z space), the coordinates of the corresponding point in the xyz space are

(x , y , z )T = Mi 1 (x , y , z i )T = [ i,

i,

(B.3)

ni](x , y , z i )T .

This point must be on the line lij in xyz space, so

niT nTj

[ i,

i,

x n i] y zi

niT 0T

0T Vj1 = 0. nTj

Vi1

(B.4)

Notice that

niT nTj

[ i,

i,

n i] =

0 nTj

i

0 nTj

i

1 . nTj ni

(B.5)

Substitute (B.5) into (B.4), the following can be obtained

nTj i x

+

zi nTj i y

niT Vi1 = 0, + nTj ni z i

nTj Vj1 = 0.

(B.6)

So for any point (x , y )T on line l ij , the following equation is true

nTj i x + nTj i y + nTj ni niT Vi1

nTj Vj1 = 0.

(B.7)

Then it has been proved that

Li (x , y ) = nTj i x + nTj i y + nTj ni niT Vi1

nTj Vj1.

(B.8)

Appendix C. Calculation of the distance from the node to a fracture The fracture of any shape will becomes a polygon after the triangulation. Therefor, only the distance from the node to a polygonal fracture needs to be considered. The first step is to calculate the foot Pf of the node P in the fracture plane, which is calculated as follows.

Pf = P

(nT P

(C.1)

nT V1) n,

where V1 is one of the vertices and n is the unit norm vector of the fracture. Then the coordinate transformation method is used to transform the foot and the fracture polygon into the fracture plane, and it is judged whether the foot falls within the fracture polygon (including the case on the boundary). If the foot is inside the polygon, the distance from the node P to the fracture is dist(P , Pf ) . Otherwise, the distance from the node to each edge of the fracture polygon is calculated, and the distance from the node P to the fracture is the minimum of these distances, ie mini dist(P , Edgei) . This can be written as follows.

dist(P , f ) =

dist(P , Pf ),

if Pf is in or on the polygon

min dist(P , Edgei),

otherwise

i

(C.2)

The distance from the node to each edge of the fracture can be calculated in a similar way.

framework for reservoirs with complex fractured media: theory, validation and case studies. J. Pet. Sci. Eng. 170, 945–957. Fumagalli, A., Pasquale, L., Zonca, S., Micheletti, S., 2016. An upscaling procedure for fractured reservoirs with embedded grids. Water Resour. Res. 52 (8), 6506–6525. Gilman, J., et al., 1986. An efficient finite-difference method for simulating phase segregation in the matrix blocks in double-porosity reservoirs. SPE Reservoir Eng. 1 (04), 403–413. Hang, S., 2015. Tetgen, a delaunay-based quality tetrahedral mesh generator. ACM Trans. Math Software 41 (2), 1–36. Hoteit, H., Chawathé, A., et al., 2016. Making field-scale chemical enhanced-oil-recovery simulations a practical reality with dynamic gridding. SPE J. 21 (06), 2–220. Hoteit, H., Firoozabadi, A., 2008. An efficient numerical model for incompressible twophase flow in fractured media. Adv. Water Resour. 31 (6), 891–905. Hui, M.H., Glass, J., Harris, D., Jia, C., 2013a. Discrete natural fracture uncertainty modelling for produced water mitigation: chuandongbei gas project, sichuan, China.

References Abdelazim, R., Rahman, S.S., 2016. Estimation of permeability of naturally fractured reservoirs by pressure transient analysis: an innovative reservoir characterization and flow simulation. J. Pet. Sci. Eng. 145, 404–422. Azim, R.A., 2016. Evaluation of water coning phenomenon in naturally fractured oil reservoirs. J. Petrol. Explor. Prod. Technol. 6 (2), 279–291. Boyle, E.J., Sams, W.N., et al., 2012. Nfflow: a reservoir simulator incorporating explicit fractures. In: SPE Western Regional Meeting. Society of Petroleum Engineers. Chen, Z., 2007. Reservoir Simulation: Mathematical Techniques in Oil Recovery, vol 77 Siam. Edelsbrunner, H., 2001. Geometry and Topology for Mesh Generation, vol 7 Cambridge University Press. Fang, W., Liu, C., Li, J., Jiang, H., Pu, J., Gu, H., Qin, X., 2018. A discrete modeling

899

Journal of Petroleum Science and Engineering 180 (2019) 886–900

Y. Zhao, et al. In: IPTC 2013: International Petroleum Technology Conference. Hui, M.-H., Kamath, J., Narr, W., Gong, B., Fitzmorris, R.E., et al., 2007. Realistic modeling of fracture networks in a giant carbonate reservoir. In: International Petroleum Technology Conference. International Petroleum Technology Conference. Hui, M.-H., Mallison, B.T., Fyrozjaee, M.H., Narr, W., et al., 2013b. The upscaling of discrete fracture models for faster, coarse-scale simulations of ior and eor processes for fractured reservoirs. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Hui, M.-H., Mallison, B.T., Lim, K.-T., et al., 2008. An innovative workflow to model fractures in a giant carbonate reservoir. In: International Petroleum Technology Conference. International Petroleum Technology Conference. Hui, M.-H.R., Karimi-Fard, M., Mallison, B., Durlofsky, L.J., et al., 2018. A general modeling framework for simulating complex recovery processes in fractured reservoirs at different resolutions. SPE J. 23 (02), 598–613. Jiang, J., Younis, R.M., et al., 2016. Hybrid coupled discrete-fracture/matrix and multicontinuum models for unconventional-reservoir simulation. SPE J. 21 (03), 1–009. Juanes, R., Samper, J., Molinero, J., 2002. A general and efficient formulation of fractures and boundary conditions in the finite element method. Int. J. Numer. Methods Eng. 54 (12), 1751–1774. Karimi-Fard, M., Durlofsky, L.J., Aziz, K., et al., 2003. An efficient discrete fracture model applicable for general purpose reservoir simulators. In: SPE Reservoir Simulation Symposium. Society of Petroleum Engineers. Kim, J.-G., Deo, M.D., 2000. Finite element, discrete-fracture model for multiphase flow in porous media. AIChE J. 46 (6), 1120–1130. Krogstad, S., Lie, K.-A., Møyner, O., Nilsen, H.M., Raynaud, X., Skaflestad, B., et al., 2015. Mrst-ad–an open-source framework for rapid prototyping and evaluation of reservoir simulation problems. In: SPE Reservoir Simulation Symposium. Society of Petroleum Engineers. Lee, S.H., Lough, M., Jensen, C., 2001. Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water Resour. Res. 37 (3), 443–455. Li, L., Lee, S.H., et al., 2006. Efficient field-scale simulation for black oil in a naturally fractured reservoir via discrete fracture networks and homogenized media. In: International Oil & Gas Conference and Exhibition in China. Society of Petroleum Engineers. Lie, K.-A., 2014. An Introduction to Reservoir Simulation Using Matlab: User Guide for the Matlab Reservoir Simulation Toolbox (Mrst). SINTEF ICT, Norway. Lie, K.-A., Krogstad, S., Ligaarden, I.S., Natvig, J.R., Nilsen, H.M., Skaflestad, B., 2012. Open-source matlab implementation of consistent discretisations on complex grids. Comput. Geosci. 16 (2), 297–322.

Lim, K.-T., Hui, M.-H., Mallison, B.T., et al., 2009. A next-generation reservoir simulator as an enabling technology for a complex discrete fracture modeling workflow. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Mallison, B.T., Hui, M.-H., Narr, W., 2010. Practical gridding algorithms for discrete fracture modeling workflows. In: ECMOR XII-12th European Conference on the Mathematics of Oil Recovery. McKoy, M.L., Sams, W.N., 1997. Tight Gas Reservoir Simulation: Modeling Discrete Irregular Strata-Bound Fracture Network Flow, Including Dynamic Recharge from the Matrix. Tech. Rep. EG and G Technical Services of West Virginia, Inc., Morgantown, WV (United States). Mi, L., Yan, B., Jiang, H., An, C., Wang, Y., Killough, J., 2017. An enhanced discrete fracture network model to simulate complex fracture distribution. J. Pet. Sci. Eng. 156, 484–496. Moinfar, A., 2013. Development of an Efficient Embedded Discrete Fracture Model for 3d Compositional Reservoir Simulation in Fractured Reservoirs. Persson, P.-O., 2005. Mesh Generation for Implicit Geometries. Ph.D. thesis. Massachusetts Institute of Technology. Persson, P.-O., Strang, G., 2004. A simple mesh generator in matlab. SIAM Rev. 46 (2), 329–345. Pruess, K., et al., 1985. A practical method for modeling fluid and heat flow in fractured porous media. Soc. Petrol. Eng. J. 25 (01), 14–26. van Golf-Racht, T.D., 1982. Fundamentals of Fractured Reservoir Engineering, vol 12 Elsevier. Warren, J., Root, P.J., et al., 1963. The Behavior of Naturally Fractured Reservoirs. Wei, L., Hadwin, J., Chaput, E., Rawnsley, K., Swaby, P., et al., 1998. Discriminating fracture patterns in fractured reservoirs by pressure transient tests. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Yan, B., Alfi, M., Cao, Y., An, C., Wang, Y., Killough, J.E., et al., 2015. Advanced multiple porosity model for fractured reservoirs. In: International Petroleum Technology Conference. International Petroleum Technology Conference. Yan, X., Huang, Z., Yao, J., Li, Y., Fan, D., 2016. An efficient embedded discrete fracture model based on mimetic finite difference method. J. Pet. Sci. Eng. 145, 11–21. Zhao, L., Jiang, H., Zhang, S., Yang, H., Sun, F., Qiao, Y., Zhao, L., Chen, W., Li, J., 2018. Modeling vertical well in field-scale discrete fracture-matrix model using a practical pseudo inner-boundary model. J. Pet. Sci. Eng. 166, 510–530. Zihua, Z., Ju, W., Rui, S., et al., 2009. Characteristics of fracture and 3d fracture network simulation of bs03. In: ISRM International Symposium on Rock MechanicsSINOROCK 2009. International Society for Rock Mechanics and Rock Engineering.

900