Simulation of a multistage fractured horizontal well in a tight oil reservoir using an embedded discrete fracture model

Abstract

As an unconventional resource, tight oil reservoirs hold abundant geological reserves and great development potential. However, tight oil reservoirs suffer from high heterogeneity and low permeability. To obtain commercial oil and gas flow, multistage hydraulically fractured horizontal well technology is often adopted. In this paper, an embedded discrete fracture model (EDFM) is adopted to study the fractured reservoir. First, the relationship between the EDFM connection list is studied and its coefficient matrix of the equation system is solved. EDFM connection list is a method to record the connections between the grids. Thereafter, the reasons for the sparsity of the coefficient matrix and the symmetry of the locations of nonzero elements in the matrix are explained, and a simple method to establish and assemble the coefficient matrix is proposed. Using this method, a simulator for tight oil reservoir with EDFM is developed. By comparing the simulation results of the self‐programming simulator with the commercial software SAPHIR, the usefulness of the simulator is verified. Finally, several simulations and sensitivity analysis are performed to determine the key factors affecting production. Modeling studies show the influence of fractures, the horizontal well, and stress sensitivity, providing very useful information for the design of hydraulic fracturing and production development programs.

INTRODUCTION

Tight oil is a petroleum resource accumulated in tight reservoirs that is interbedded with or adjacent to oil‐generating strata. The petroleum industry has not formed a unified standard for the definition of tight oil reservoirs.) Due to the reservoir and technical conditions of each country, British scholars define the permeability of tight oil reservoirs as <1 mD, whereas values of 0.6 mD and 0.1 mD are used in Germany and the United States, respectively. However, the consensus is that the permeability is extremely low because of the nanometer‐scale pore size and complex microstructure. It is also found that tight oil reservoirs have strong stress sensitivity, which means that, during development, the formation pressure will gradually decrease, and the already low permeability will drop even lower. To understand the complex fluid transport in unconventional reservoirs, many researchers have done a lot of significant work. Overall, tight oil reservoirs have poor physical properties. Generally, a single tight oil well has natural production capacities below the lower limit of commercial oil flow. However, because these reservoirs are interbedded with or adjacent to oil‐generating strata, in general, the oil properties are good, and amount of the original oil in place is large. In recent years, with the application of horizontal wells and hydraulic fracturing technology, the production of tight oil has been greatly improved.

Horizontal well multistage fracturing technology can stimulate the formation and increase the reservoir contact significantly. Tight reservoirs usually contain natural fractures, which lead to strong reservoir heterogeneity. Therefore, we must accurately model the complex flow of both hydraulic fractures and natural fractures with multiple orientations, lengths, and conductivities.

Barenblatt and Zheltov pioneered the concept of dual porosity. They superimposed the two media of the fracture and the matrix. The fluid flows relatively independently in each space, and exchange occurs between the two media. Based on this, Warren and Root developed the dual‐porosity model. They introduced the concept of interporosity flow to describe the fluid exchange between fractures and the matrix. Later, a “multiple interacting continua” (MINC) method was proposed to study fluid and heat flow in fractured porous media. The MINC method, which includes a series of “nested” matrix control volumes created in fracture control volumes, can successfully model the multiphase flow and, compared to explicit discretization, effectively reduce the extra cells introduced by subgridding. The traditional dual‐porosity modes ignored the flow between matrices. Blaskovich et al. and Hill and Thomas proposed an extended model, called the dual‐porosity dual‐permeability (DPDK) model. Compared to the traditional dual‐porosity mode, DPDK added the matrix‐matrix connection to account for the flow between matrix blocks. However, the dual‐porosity mode is only suitable for reservoirs with highly connected and small‐scale fractures. Moinfar et al. showed that the dual‐porosity model leads to unsatisfying solutions in the presence of highly localized anisotropy and preferential channeling.

Subsequently, Karimi‐Fard et al. developed the discrete fracture model (DFM) based on control‐volume formulations. They employed unstructured grids to characterize the geometry of the fracture networks. Compared to the dual‐porosity model, the DFM has the advantage that it can use realistic fracture geometry properties, so it can accurately describe the flow effects of each fracture. Zhang et al. studied fractured unconventional reservoirs using the DFM. However, the DFM also has limitations. The quality of the generated unstructured mesh is highly dependent on the mathematical algorithm applied, and the mesh can be very complicated for some situations. For reservoirs with numerous fractures, a large number of refined grids are generated around the fractures, which leads to a sharp increase in the calculation time.

Lee et al. proposed the embedded discrete fracture model (EDFM). The EDFM uses an orthogonal structured grid for the mesh matrix, and the fractures are divided into a series of segments based on the boundary of the matrix cell. These segments are the fracture grids in the EDFM. The EDFM is very flexible for handing various kinds of intersections, such as fracture‐fracture and fracture‐well intersections. The interporosity between two systems is calculated using the Peaceman‐like formulation, which appears as the source/sink term in the governing equations. Furthermore, the EDFM uses orthogonal structured grids, so the computational speed is greatly improved, especially when simulating flows of complex distributed fractures. Based on the above advantages, many scholars have performed numerous simulations using the EDFM, demonstrating its good adaptability for simulating fractured reservoirs.

The EDFM can describe the interporosity flow between different media, so each cell is connected with different numbers and types of cells. To understand and record these connections proposed the non‐neighboring connection (NNC) method, and Jiang et al. proposed the connection list, but they did not study their methods further. However, these connections are crucial in transforming the EDFM from a physical model into a mathematical and programming model, making clear connection theory important to the EDFM. Based on the above two methods, we studied the relationship between the connection of two cells and the coefficient matrix of the equation system to be solved, and we proposed a very useful method for establishing the coefficient matrix.

Using this method, we established an EDFM simulator to study tight reservoirs with multistage fractured horizontal wells. We verified the simulator results by using SAPHIR (a commercial numerical simulator). Our EDFM simulator can accurately simulate interporosity flow among the matrix, horizontal wells, hydraulic fractures, and natural fractures. To provide guidance for production optimization and to determine the influence of the reservoir, hydraulic fractures, natural fractures, and the horizontal well, we simulated the production process of multistage fracturing horizontal wells using two‐dimensional (2D) and three‐dimensional (3D) models with and without random natural fractures.

PHYSICAL MODELING of FLOW

Model assumptions

The basic assumptions for modeling multistage fractured horizontal wells in tight reservoirs are as follows: (a) The reservoir matrix is homogeneous and of equal thickness. (b) Oil flows isothermally, and we ignore the change in oil viscosity with pressure. (c) Gravity and the stress sensitivity of matrix permeability are considered. (d) The reservoir has closed boundaries. (e) Fractures are vertical.

Governing equations

There are two systems, the matrix and the fractures, that are fluid reservoirs and flow channels. Combined with Darcy's law and the law of conservation of mass, the fluid control equations in the two systems are [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] where ${\varphi }_m$ and ${\varphi }_f$ represent the porosities of the matrix and the fractures, respectively (both dimensionless), B_m and B_f represent the volume factors of oil in the matrix and in the fractures, respectively (in m³/m³), β_c (=86.4 × 10⁻⁶) represents the transmissibility factor, μ_m and μ_f represent the viscosities of oil in the matrix and in the fractures, respectively (in Pa·s), k_m and k_f represent the permeabilities of the matrix and the fractures, respectively (in D), p_m and p_f represent the pressures of oil in the matrix and in the fractures, respectively (in kPa), ${\rho }_m$ and ${\rho }_f$ represent the densities of oil in the matrix and in the fractures, respectively (in kg/m³), g (=9.8066 m/s²) represents the acceleration of gravity, Z_m and Z_f represent the relative heights in the matrix and in the fractures, respectively (in m), q_m and q_f represent the net flows (being positive when there is a net inflow) (in m³/day), q_mf and q_mw represent the matrix‐fracture and matrix‐well mass exchanges (being positive for net inflow to the matrix), respectively (in m³/day), and q_fm, q_ff, and q_fw represent the fracture‐matrix, fracture‐fracture, and fracture‐well mass exchanges (being positive for net inflow to the fractures), respectively (in m³/day). The matrix‐well mass exchange is calculated using the Peaceman model, and others are determined using the EDFM.

Stress sensitivity in tight reservoirs

Stress sensitivity is one of the important characteristics of tight reservoirs. To highlight the tight reservoir flow mechanism, we considered the stress‐sensitive effect of matrix permeability. Scholars have proposed different types of stress‐sensitive models, including binomial relations, logarithmic relations, power‐law relations, and exponential relationships. In our work, we adopted the index relationship model[Image Omitted. See PDF] where α represents the stress‐sensitive coefficient (in kPa⁻¹), p₀ represents the reference pressure, which is the initial reservoir pressure in our work (in kPa), and $k_m^0$ represents the matrix permeability at reference pressure (in D).

NUMERICAL MODEL AND SIMULATION APPROACH

Embedded discrete fracture model

The embedded discrete fracture model (EDFM) borrows the dual‐medium concept from conventional dual‐continuum models and also incorporates the effect of each fracture explicitly. Using the concept of the well index (WI) in the finite‐difference method, the matrix‐fracture (q_mf) and fracture‐fracture (q_ff) interporosity flows are treated as a part of the source/sink item.

In the EDFM, when a fracture passes through the boundary of the matrix grid, the fracture is divided into fracture segments of different shapes and sizes, based on the grid boundary. These new fracture segments are fracture grids, which have their own governing equations and grid parameters. As a consequence, the fracture grids have the same status as matrix grids when solving the fluid control equations.

Figure A shows a 2D EDFM example in which two vertical fractures communicate with each other and one of them is intersected by a horizontal wellbore. Figure B shows the connection list for this scenario. In addition to the traditional finite‐difference connection, we can see that each fracture grid is embedded in a matrix grid. Two fractures intersect at f2–f4, and the horizontal well is connected to a fracture at f1–w1. In Figure B, different line types show different types of connectivity. The idea of the EDFM connectivity list comes from Jiang et al. and is similar to the NNC concept. When solving the algebraic equations, we found that constructing such a connection list is extremely helpful for positioning the nonzero elements of the sparse matrix. We will illustrate how to understand and establish the sparse matrix in the next section.

Example of an EDFM connection list

After establishing the connection list, the flow rate between the new connections can be calculated using [Image Omitted. See PDF] [Image Omitted. See PDF] Where q_CON represents the flow rate (in m³/day), T_CON represents the transmissibility of the connection (in m³·s⁻¹·kPa⁻¹), Δp represents the pressure difference between the connections (in kPa), β_c (=86.4 × 10⁻⁶) represents the transmissibility factor, k_CON represents the harmonic mean permeability of the connection (in D), A_CON represents the contact area of the connection (in m²), d_CON represents the feature distances of the connection (in m), μ_CON and B_CON represent the viscosity (in Pa·s) and volume factor of the connection (in m³/m³), respectively, G_CON represents the transmissibility geometric factor (in D·m), and f_CON represents a weakly nonlinear transmissibility term (in Pa⁻¹·s⁻¹). For fractures not fully penetrating the grid, Li and Lee assumed that T_CON is linearly proportional to the facture length inside the grid.

In the connection list, as shown in Figure B, a new connection can be divided into four types. Their geometric factors can be obtained as follows.

CON‐I: the connection between the fracture grid and its embedded matrix grid, such as m1–f1 in Figure B: [Image Omitted. See PDF] where A_mf (in m²) represents the contact area between the fracture and its embedded matrix, which is 2A_f (because both sides of the fracture are in contact with the matrix), k_mf (in D) represents the harmonic mean permeability of the matrix and the fracture, which is approximately equal to k_m (because the fracture permeability is much greater than the matrix permeability), and d_mf represents the average normal distance (in m).

For d_mf, Lee et al. assumed that the fluid in the matrix perpendicularly flows to the fracture. Therefore, the pressure around a fracture is linearly distributed, and the average normal distance can be computed as [Image Omitted. See PDF] where dv and V represent the volumes of the microelement and grid, respectively (in m³), and ${\mathit{x}}_n$ is the normal distance between the microelement and the fracture (in m).

CON‐II: the connection resulting from the intersection of two different fractures, such as f2‐f4 in Figure B:[Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] Where L_int represents the length of the fracture intersection line (in m), ${\omega }_f$ and k_f represent the fracture aperture and permeability, respectively (in m and D, respectively), d_f represents the average of normal distances (in m) from the center of the fracture subsegments (located in each side of the intersection line) to the intersection line (with its microintegral method being similar to Equation ), and subscripts 1 and 2 refer to different fractures.

CON‐III: the connection between adjacent fracture grids inside the same fracture, such as f1–f2 in Figure B:[Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] where k_f represents the fracture permeability (in D), A_c represents the common area of adjacent fracture grids (in m²), d_seg represents the distance between the centroid of the fracture grid and the common surface of two fracture grids (in m), and subscripts 1 and 2 refer to different fracture grids inside the same fracture.

CON‐IV: the connection between the fracture grid and the well grid, such as w1–f1 in Figure B:[Image Omitted. See PDF] [Image Omitted. See PDF] where k_f represents the fracture permeability (in D), ${\omega }_f$ represents the fracture aperture (in m), L_f represents the length of the fracture grid (in m), h_f represents the height of the fracture grid (in m), r_e represents the equivalent well grid radius (in m), r_w represents the wellbore radius (in m), and Δθ represents the central angle of the well within the fracture (in rad). In general, the wellbore passes through one fracture grid, so $\mathrm{\Delta }\theta =2\pi$.

Constructing the coefficient matrix

Because the governing equations cannot be calculated directly, we must choose some methods to discretize and linearize their time term and space terms. In our work, we used the fully implicit method and single‐point upstream weighting method. We do not give the details of the derivation in this paper, but we will introduce a new method to understand the relationship between the coefficient matrix and the EDFM connection list, which can help beginners study EDFM effectively.

The discrete result for a matrix cell (i, j, k) can be expressed as[Image Omitted. See PDF] where $R_{m_{i,j,k}}$ represents a constant term, $\overline{\mathit{\delta }}{\mathit{p}}_{\mathit{m}{\mathit{m}}_{\mathit{i}\mathbf{,}\mathit{j}\mathbf{,}\mathit{k}}}$ and ${\mathit{C}}_{\mathit{m}{\mathit{m}}_{\mathit{i}\mathbf{,}\mathit{j}\mathbf{,}\mathit{k}}}$ represent the vectors of unknown pressure difference and its coefficient, respectively, which are derived from matrix‐matrix connections, $\overline{\mathit{\delta }}{\mathit{p}}_{\mathit{m}{\mathit{f}}_{\mathit{i}\mathbf{,}\mathit{j}\mathbf{,}\mathit{k}}}$and $C_{m{f}_{i,j,k}}$ represent the vectors of unknown pressure difference and its coefficient, respectively, which are derived from matrix‐fracture connections, and $\overline{\delta }{p}_{m{w}_{i,j,k}}$ and $W_{m{w}_{i,j,k}}$ represent the unknown pressure difference and its coefficient, respectively, which are derived from matrix‐well connection. If the well produces under constant pressure, $\overline{\delta }{p}_{w{f}_{i,j,k}}=0$. $\overline{\delta }{p}_{m_{i,j,k}}^{}$ and $C_{m_{i,j,k}}^{}$ represent the unknown pressure difference and its coefficient, respectively, which are derived from every connection associated with the matrix cell (i,j,k) and the accumulation of the matrix cell.

The discrete result for a fracture cell (u, v) can be expressed as [Image Omitted. See PDF] where $R_{f_{u,v}}$ represents a constant term, $\overline{\mathit{\delta }}{\mathit{p}}_{{\mathit{f}}_{\mathit{u}\mathbf{,}\mathit{v}}}^{\mathrm{ins}}$ and ${\mathit{C}}_{{\mathit{f}}_{\mathit{u}\mathbf{,}\mathit{v}}}^{\mathrm{ins}}$ represent the vectors of unknown pressure difference and its coefficient, respectively, which are derived from fracture‐fracture connections inside the same fracture, $\overline{\delta }{p}_{f{m}_{u,v}}$ and $C_{f{m}_{u,v}}$ represent the unknown pressure difference and its coefficient, respectively, which are derived from fracture‐matrix connection, $\overline{\mathit{\delta }}{\mathit{p}}_{\mathit{f}{\mathit{f}}_{\mathit{u}\mathbf{,}\mathit{v}}}$ and ${\mathit{C}}_{\mathit{f}{\mathit{f}}_{\mathit{u}\mathbf{,}\mathit{v}}}$ represent the vectors of unknown pressure difference and its coefficient, respectively, which are derived from fracture‐fracture connections in different fractures, and $\overline{\delta }{p}_{f{w}_{u,v}}$ and $W_{f{w}_{u,v}}$ represent the unknown pressure difference and its coefficient, respectively, which are derived from the fracture‐well connection. If the well produces under constant pressure, $\overline{\delta }{p}_{f{w}_{u,v}}=0$. $\overline{\delta }{p}_{f_{u,v}}^{}$ and $C_{f_{u,v}}^{}$ represent the unknown pressure difference and its coefficient, respectively, which are derived from every connection associated with the fracture cell (u,v) and the accumulation of the fracture cell.

After assembling discrete equations of all matrix cells and fracture cells, we can obtain a system of linear equations to calculate [Image Omitted. See PDF] where A represents the coefficient matrix (a large sparse matrix), $\overline{\mathit{\delta }}\mathit{p}$represents the vector of unknown pressure difference between two iteration steps, which is assembled from an unknown matrix and the fracture pressure difference ($\overline{\mathit{\delta }}{\mathit{p}}_m$ and $\overline{\mathit{\delta }}{\mathit{p}}_f$) in the EDFM, and b represents a constant vector.

For the model of Figure , we assumed constant‐pressure well production and used Equation 20 to get[Image Omitted. See PDF] (where only the nonzero elements are shown). From Equation 21, we can see that the coefficient matrix is sparse and the locations of nonzero elements in the matrix are symmetric. However, when the upstream method is employed for the flux term, the coefficient matrix is not symmetric. To be specific, the volume factor B at the interface of two cells is related to the larger pressure of the two cells and will contribute to the coefficient matrix.

Based on the physical meaning, this matrix can be divided into four parts. The upper left part is the traditional matrix part, called the M‐M part. The upper right part and lower left part are M‐F and F‐M parts, describing the connections between the matrix cells and the fracture cells. The lower right part is the F‐F part, which describes the connections of two fracture cells in a fracture or between two different fractures, which can be further subdivided based on different connections according to different fractures.

Algebraic equations for the ith cell are represented by the ith row in the coefficient matrix. The a(i, j) element, which is the element of the ith row and jth column in the coefficient matrix, represents the relationship between cell i and cell j, and the cell can be a matrix cell or a fracture cell. Therefore, if a(i, j) is not equal to zero, there is a connection between cell i and cell j. The location of a(i, j) in different parts of the matrix also means different connections. If cell i is connected to cell j, then cell j will also be connected to cell i, which makes both a(i, j) and a(j, i) nonzero. This makes the locations of nonzero elements symmetric in the matrix. Besides, a cell only connects to a few cells, which gives the matrix only a few nonzero elements per line, so the matrix is sparse.

Based on the above discussion of the connection list and matrix properties, we created a method to assemble this matrix. First, we sort the grid number of the matrix, fracture, and well cells in a natural order. Second, we construct the connection list, which means, for every single cell, find out which grids it is connected to and record the types of these connections. Then, we calculate the transmissibility and other factors that are required for the element in the matrix. Finally, according to the connection list, we assemble the coefficient matrix.

In this way, we built the relationships among the EDFM connectivity list, the physical meaning, and the location of the nonzero elements in the matrix. It then becomes very simple to understand and assemble the matrix. Similarly, this method can also be applied to a single‐porosity model and a dual‐porosity continuous medium model. For the single‐porosity model, there is only one kind of cell, the matrix cell. Therefore, we only need to construct the upper left part of the coefficient matrix in Equation 21. It is very simple. For the dual‐porosity continuous medium model (including DPDK and MINC), in addition to the connections between matrix cells, some matrix cells also connect to the fracture cell, we can also construct the coefficient matrix using the connection list, just like the EDFM.

Generation of natural fractures

Based on a comparison of random fracture generation methods, we used the most flexible method to randomly generate fractures. It is very easy to change the generation of fractures according to realistic distributions. Random fractures can be mathematically described as[Image Omitted. See PDF] where rand represents a random number in the interval [0,1], and its subscript is the number of times it generates the random number, P_x and P_y represent the coordinates of random point (x, y) in the reservoir, r_el and r_ew represent the length and width, respectively, of the rectangular reservoir (in m), θ_p and [θ] represent the random azimuth and the range of azimuth, respectively (in rad), and L_p and [L] represent the random length and the range of length, respectively (in m).

MODEL VERIFICATION

A horizontal well model with four hydraulic fractures was simulated using SAPHIR (a commercial numerical simulation software program) and our EDFM simulator. The simulation parameters are listed in Table . The EDFM uses equal time steps of 1 day. SAPHIR uses an adaptive step size, which is a combination of arithmetic and geometric progression. Figure A,B shows the EDFM and SAPHIR meshes. In this scenario, the total number of grid cells is about 1900 in Figure A, whereas the number is more than 5000 in Figure B. This is because, for the DFM, there are refined grids around wells and fractures, and there are only four fractures and one horizontal well in this simple scenario. However, the advantage of the EDFM is that it can handle a scenario with a large number of complex fractures without refined grids around each fractures, just like the scenarios as shown in Section 5.2.

Model simulation parameters

Parameter	Value	Unit
Reservoir dimensions (x, y, z)	1000, 500, 10	m
Grid block size (x, y, z)	10, 10, 10	m
Initial reservoir pressure	20	MPa
Horizontal wellbore length	660	m
Fracture half‐length	50	m
Fracture height	10	m
Fracture aperture	0.01	m
Fracture permeability	4 × 104	mD
Well radius	0.1	m
Initial oil volume factor	1.3	dimensionless
Oil density	0.8	g/cm3
Oil viscosity	20	mPa·s
Oil compressibility coefficient	1 × 10−3	MPa−1
Matrix porosity	0.1	Dimensionless
Matrix permeability	0.5	mD
Stress‐sensitive coefficient	10−3	MPa−1
Producer BHP	15	MPa

¹In model verification, the stress sensitivity is not considered.

A, Mesh for SAPHIR. B, Mesh for EDFM

For this scenario with only four fractures and one horizontal well, the productions simulated by the two simulators are almost the same, as shown in Figure . There is only a slight difference in the early stage, which is probably due to the simplification for interporosity flow for the EDFM. The EDFM assumes that the fluid in the matrix perpendicularly flows to the fracture. Based on the assumption, we can calculate the average normal distance d_mf in Equation  to represent the distance between a fracture cell and its embedded matrix cell. However, for the control‐volume finite‐difference formulation, the average pressures of the matrix cell and fracture cell are represented at the centers of the matrix cell and fracture cell. This simplification method for the pressure ignores the variation of pressure within the cells. Therefore, it is not very precise, and the average normal distance d_mf is not equal to the physical distance between matrix pressure point and fracture pressure point (the distance between the centers of the matrix cell and fracture cell), which causes a slight difference in production, as shown in Figure . However, overall, the difference is small, so we can conclude that the EDFM is consistent with the SAPHIR software.

Comparison of the EDFM and SAPHIR results

COMPREHENSIVE SENSITIVITY STUDY

Sensitivity analysis is a quantitative method to study the effect of parameter variation on the model results. It can identify the key reservoir and fracture parameters and provide important guidance for designing the development program. In our work, based on two tight oil reservoir scenarios, we conducted several simulation studies.

Model of tight oil reservoirs without natural fractures

We studied the sensitivity of the multistage fractured horizontal well model to hydraulic fracture number, hydraulic fracture half‐length, combination of fracture number and half‐length, fracture conductivity, stress sensitivity of the matrix permeability, and fracture height. A 3D model was adopted for the fracture height, and the other parameters used a 2D model. The formation and well parameters for the 2D model are listed in Table . In the 3D model, except for the formation thickness of 100 m, which is divided into 10 grids in the z direction, the simulated parameters were the same as in the 2D model.

Number of hydraulic fractures

We studied the effect of the number of hydraulic fractures (4, 6, 8, 12, and 14). The result was the same as our qualitative understanding. Production is positively correlated with the number of fractures. Figure A shows the relationship between daily and cumulative production and the number of fractures. Figure B shows that there is an “inflection point” in the relationship between fracture number and cumulative oil production (COP) between 6 and 8 fractures. Beyond 8, the effect of fracture number on productivity is reduced. Combined with the pressure profiles (Figure C), we can see that the pressure distributions at n = 12 and n = 14 are similar, and the sizes of their stimulated reservoir volumes (SRVs) (the oval pressure drop areas) are almost identical. Therefore, there is no need to design too many hydraulic fractures. Any number slightly greater than that at which the inflection point occurs is enough.

A, Hydraulic fracture number. B, Cumulative oil production. C, Pressure profiles at 400 days

Half‐length of hydraulic fractures

We varied the half‐length of hydraulic fractures (20 m, 50 m, 100 m, 200 m, and 250 m) to study the effect of half‐length on production. From the relationship of fracture half‐length and COP (Figure B), we found an inflection point between half‐lengths of 200 m and 250 m. Beyond 200 m, the effect of fracture half‐length on productivity is reduced. Moreover, when the fractures are completely through the reservoir, there is a theoretical upper limit for the influence of fracture half‐length on productivity, as shown by the black dotted line in Figure A. The pressure profiles (Figure C) show that increasing the length of hydraulic fractures can effectively expand the SRV and connect a larger reservoir area. Thus, longer hydraulic fractures should be designed as much as economically possible.

A, Hydraulic fracture half‐length. B, Cumulative oil production. C, Pressure profiles at 100 days

Half‐length and the number of hydraulic fractures

We analyzed the combined influence of fracture half‐length and the number of fractures. To control the variables, we set the total length of the fractures to 800 m, and the fracture combinations were 200 m × 4, 100 m × 8, 80 m × 10, and 50 m × 16. The simulation results (Figure A,B show that within 130 days, the difference in oil production is minor except for the combination of 200 m × 4. As the time increases, the difference gradually becomes smaller. By about 500 days, production from each combination is almost the same.

A, Hydraulic fracture half‐length. B, Cumulative oil production. C, Pressure profiles at 300 days

From the pressure profiles (Figure C), we can see that the matrix permeability is so low that if the factures are too far apart, the areas between them will not be effectively connected. It will take a long time to reach the quasi‐radial flow stage. Therefore, fracturing of the horizontal well should not be designed with a long half‐length and a wide fracture distance, which would lead to low initial production and is not conducive to cost recovery.

Conductivity of hydraulic fractures

A sensitivity study was performed on hydraulic fracture conductivities (fracture permeability × fracture aperture) of 100 mD·m, 200 mD·m, 400 mD·m, and 800 mD·m as well as infinite. According to the relationship between the oil production and the hydraulic fracture conductivity, as shown in Figure A,B, we can see that as long as the fracture conductivity exceeds 400 mD·m, the effect of fracture conductivity is minor. However, as for hydraulic fractures, there is always something (proppant, raised part of a rock face, etc.) to support fracture faces. So, in practice, the conductivity of the fracture will not be very low. From the pressure profiles (Figure C), we can see the fracturing conductivity mainly affects the pressure in the SRV during the early stage. However, the size of the SRV remains almost the same. Therefore, the difference in COP is not large, and hydraulic fracture conductivity does not have much impact on production.

A, Hydraulic fracture conductivity. B, Cumulative oil production. C, Pressure profiles at 200 days

Length of the horizontal well

We also studied the effect of the length of the horizontal well (400 m, 550 m, 660 m, and 750 m). From the simulation results (Figure A), we found that, in the early stages of production (about 0 days), well production is almost identical. After that, the longer the horizontal well is, the higher the production becomes. A length of 400 m had the biggest decline in daily production. According to the pressure profiles (Figure B), the pressure in the SRV for this case dropped to a very low value, because most of the oil has been produced. Therefore, we can conclude that in the later stage of production, the longer the horizontal well is, the longer the long axis of the ellipse quasi‐radial flow is, and the higher the production will be.

A, Length of horizontal well. B, Pressure profiles at 300 days

Stress sensitivity of matrix permeability

To investigate the effect of stress sensitivity of matrix permeability, we set sensitivity coefficients of 0 MPa⁻¹, 0.001 MPa⁻¹, 0.01 MPa⁻¹, 0.05 MPa⁻¹, 0.1 MPa⁻¹, and 0.2 MPa⁻¹. From the simulation results (Figure A,B, we can see that as the stress sensitivity coefficient increases, oil production decreases, which results in less pressure reduction in the SRV. This occurs because, during the development of a tight reservoir, the reservoir pressure decreases and the pore channels become smaller, resulting in an increase in flow resistance, a decrease in matrix permeability, and ultimately a decrease in production. Therefore, to effectively develop tight reservoirs, we should control the difference in production pressure and minimize the influence of stress sensitivity on production near the low‐pressure zone.

A, Stress sensitivity of matrix permeability. B, Pressure profiles at 400 days

Height of hydraulic fractures

We used a 3D model to analyze the sensitivity of hydraulic fracture height. The 3D model is basically the same as the 2D model, except that there are multiple layers in the z direction to account for the influence of gravity and the flow between longitudinal layers. For the 3D model, we set the total thickness of the oil layer to 100 m and fracture heights to 40 m, 60 m, 80 m, and 100 m, as shown in Figure A. To save calculation time, we set the time step to 2 days.

A, 3D mesh. B, Hydraulic fracture height. C, Pressure profiles at 100 days

From Figure B, we can see that the production rate is positively correlated with fracture height. Figure C shows the pressure profiles at 100 days. However, during the later production period, daily oil production is almost the same. We determined that due to the different sizes of the SRVs, there is a large difference in the initial production. During the later period, the oil inside the SRV is almost all recovered, and the oil outside the SRV pseudoradially flows into the SRV, mainly coming from the horizontal direction. For the stage of quasi‐radial flow, the main controlling factor is the extremely low matrix permeability, which makes the later production rate almost the same. Therefore, we can conclude that the height of the fractures significantly affects production during the early linear flow stage but has little effect on the quasi‐radial flow stage during the later stage of production.

Model of tight oil reservoirs with natural fractures

Many tight reservoirs contain several natural fractures that greatly enhance the flow capacity of the formation. Some natural factures are connected to each other, some are connected with the wellbore, and some are connected with the hydraulic factures. Therefore, the reservoir structure is very complicated, and there is strong heterogeneity. To accurately simulate the production and analyze the influence of natural fractures, a 2D EDFM was used to simulate the production of tight reservoirs with natural fractures. The parameters of the well, hydraulic fractures, and formation are the same as given in Table . The natural fracture parameters are listed in Table .

Natural fracture parameters

Parameter	Value	Unit
Natural fracture permeability	50	mD
Natural fracture aperture	0.01	m
Natural fracture height	10	m

Number of natural fractures

We considered 5, 10, 20, and 25 natural fractures in the reservoir. To control the variables, these generated natural fractures were not thought intersect with hydraulic fractures or the wellbore, as shown in Figure A. From the simulation results (Figure B), we can see that natural fractures can greatly increase oil production. According to the pressure profiles (Figure C), we can see that natural fractures can expand the pressure drop area of the reservoir, which means that they increase the flow capacity in the reservoir and expand the SRV, thereby increasing oil production.

A, Distribution of natural fractures. B, Number of natural fractures. C, Pressure profiles at 300 days

We also found an interesting phenomenon for which the natural fractures near and away from the hydraulic fracture zone have different effects on the shape of the isobaric surface. Nearby natural fractures cause the isostatic pressure to sag toward the well, whereas distant ones cause the isostatic pressure to bulge toward the outer boundary. Because the permeability of natural fractures is much greater than that of the matrix, the flow rate in the natural fractures is higher than that in the matrix, so as the fluid flows from the reservoir to the wellbore, the fluid in the matrix will preferentially flow to the nearby natural fractures. If the natural fracture is close to the boundary where the pressure difference is small, the fracture will not get enough energy (fluid) to replenish itself from the surrounding matrix. Therefore, it forms a phenomenon highlighted by the blue frame in Figure . In contrast, for the near‐well zone, the pressure difference is large and the fracture can get energy (fluid) and be replenished quickly, so the pressure remains high (as indicated by the red frame in Figure ).

Effects of different types of natural fractures

Number of natural fractures connected to hydraulic fractures

From the simulation described in Section 5.2.1, it was found that the natural fractures near the well can make the isobaric surface sag toward the well, which is beneficial to oil recovery. To further study the influence of natural fractures in the near‐well zone, we simulated production for different numbers of natural fractures connected to hydraulic fractures (5, 10, 15, and 20) under the total number of 25, as shown in Figure A.

A, Distribution of natural fractures. B, Different numbers of connections. C, Pressure profiles at 300 days

We can see from the simulation results (Figure B) that the connected natural fractures greatly improve oil production, compared with disconnected fractures (Figure B). From the pressure profiles (Figure C), we can see that as the number of connected fractures increases, the isobaric surface is more and more inwardly recessed and even “broken up,” which means there might occur multiple closed isobaric surfaces instead of one closed isobaric surfaces. This indicates that the connected natural fractures in the near‐well zone greatly improve the flow conditions of the tight reservoir, so that the production pressure drop can quickly be transmitted throughout the reservoir. Therefore, to improve the hydraulic fracturing effect, the horizontal well should be drilled in the zone where natural fractures are highly developed.

CONCLUSIONS

In our work, we studied the relationship between the EDFM connection list and the coefficient matrix. We explained the reasons for the sparsity of the coefficient matrix and the symmetry of the locations of nonzero elements in the matrix, and we proposed a simple method to establish the coefficient matrix. Based on this method, we developed an EDFM simulator for a multistage fractured horizontal well in a tight oil reservoir that includes consideration of the effects of generated natural fractures and the pressure sensitivity of the matrix. Through simulation verification, we confirmed that our EDFM simulator that uses structured grids can achieve accuracy similar to that of SAPHIR, which uses unstructured grids. Finally, we conducted comprehensive modeling studies to understand the key properties that affect production performance. From the simulated results, we reached the following conclusions:

1. In a certain range, increasing the number and half‐length of hydraulic fractures and the length of horizontal wells will greatly increase production. Beyond that range, however, the effect will be reduced. Therefore, in actual production, it is very important to find the inflection point of the effect.
2. When we set the total length of all hydraulic fractures to certain value, we find that, if the fracture combinations have a long half‐length with a wide fracture distance, the areas between fractures will not be effectively connected, which leads to very low initial production and is not conducive to cost recovery.
3. Hydraulic fracture conductivity does not significantly influence production unless the fracture conductivity is very low. Therefore, there is no need to design excessive hydraulic fracture conductivity.
4. The stress sensitivity of the matrix has a significant impact. We should control the difference in production pressure and minimize the influence of stress sensitivity in the low‐pressure zone near the well.
5. The height of the fractures significantly affects production during the early linear flow stage, but it has little effect on the quasi‐radial flow in the later stage.
6. Natural fractures have a great impact on production. Nearby natural fractures will cause isostatic pressure to sag toward the well, whereas distant natural fractures will cause the isostatic pressure to bulge toward the outer boundary. As the number of natural fractures connected to hydraulic fractures increases, the isobaric surface will be more and more inwardly recessed and even “broken up,” leading to a production pressure drop that can quickly be transmitted throughout the reservoir, which is very beneficial for production.

ACKNOWLEDGMENT

This work was supported by the National Natural Science Foundation of China (Grant No. 51704247 and 51874251), the National Natural Science Foundation of China (Key Program) (Grant No. 51534006), the PetroChina Innovation Foundation (No. 2018D‐5007‐0218), International S&T Cooperation Program of Sichuan Province (Grant No. 2019YFH0169), and the Deep Marine Shale Gas Efficient Development Overseas Expertise Introduction Center for Discipline Innovation (111 Center). The authors would also like to thank the reviewers and editors whose critical comments were very helpful in preparing this article.