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

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.

oil-generating strata, in general, the oil properties are good, and amount of the original oil in place is large. 8 In recent years, with the application of horizontal wells and hydraulic fracturing technology, the production of tight oil has been greatly improved. 9 Horizontal well multistage fracturing technology can stimulate the formation and increase the reservoir contact significantly. [10][11][12][13][14] 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 15 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 16 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 [17][18][19] 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, 20 effectively reduce the extra cells introduced by subgridding. The traditional dual-porosity modes ignored the flow between matrices. Blaskovich et al. 21 and Hill and Thomas 22 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. 23 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. 24 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. 25 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. 26,27 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 28 . 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, 23,[29][30][31][32][33][34][35][36][37][38][39][40] many scholars have performed numerous simulations using the EDFM, demonstrating its good adaptability for simulating fractured reservoirs. 23,[29][30][31][32][33][34][35][36][37][38][39][40] 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 30 proposed the non-neighboring connection (NNC) method, and Jiang et al. 31 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.

| 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 where m and 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 3 /m 3 ), β c (=86.4 × 10 −6 ) 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), m and f represent the densities of oil in the matrix and in the fractures, respectively (in kg/m 3 ), g (=9.8066 m/s 2 ) 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 3 /day), q mf and q mw represent the matrixfracture and matrix-well mass exchanges (being positive for net inflow to the matrix), respectively (in m 3 /day), and q fm , q ff , and q fw represent the fracture-matrix, fracturefracture, and fracture-well mass exchanges (being positive for net inflow to the fractures), respectively (in m 3 / day). The matrix-well mass exchange is calculated using the Peaceman model, 41 and others are determined using the EDFM.

| Stress sensitivity in tight reservoirs
Stress sensitivity is one of the important characteristics of tight reservoirs. 42 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. 3,[42][43][44] In our work, we adopted the index relationship model where α represents the stress-sensitive coefficient (in kPa −1 ), p 0 represents the reference pressure, which is the initial reservoir pressure in our work (in kPa), and k 0 m represents the matrix permeability at reference pressure (in D).

| 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. 30 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 1A 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 1B 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 1B, different line types show different types of connectivity. The idea of the EDFM connectivity list comes from Jiang et al. 31 and is similar to the NNC concept. 30,45 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.
After establishing the connection list, the flow rate between the new connections can be calculated using Where q CON represents the flow rate (in m 3 /day), T CON represents the transmissibility of the connection (in m 3 ·s −1 ·kPa −1 ), Δp represents the pressure difference between the connections (in kPa), β c (=86.4 × 10 −6 ) 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 2 ), 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 3 /m 3 ), respectively, G CON represents the transmissibility geometric factor (in D·m), and f CON represents a weakly nonlinear transmissibility term (in Pa −1 ·s −1 ). For fractures not fully penetrating the grid, Li and Lee 28 assumed that T CON is linearly proportional to the facture length inside the grid.
In the connection list, as shown in Figure 1B, a new connection can be divided into four types. Their geometric factors can be obtained as follows. 30,45 CON-I: the connection between the fracture grid and its embedded matrix grid, such as m1-f1 in Figure 1B: where A mf (in m 2 ) 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. 26,27 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 where dv and V represent the volumes of the microelement and grid, respectively (in m 3 ), and 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 1B: Where L int represents the length of the fracture intersection line (in m), 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 1B: where k f represents the fracture permeability (in D), A c represents the common area of adjacent fracture grids (in m 2 ), 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 1B: where k f represents the fracture permeability (in D), 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 Δ = 2 .

| 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 where R m i,j,k represents a constant term, p mm i,j,k and C mm i,j,k represent the vectors of unknown pressure difference and its coefficient, respectively, which are derived from matrixmatrix connections, p mf i,j,k and C mf i,j,k represent the vectors of unknown pressure difference and its coefficient, respectively, which are derived from matrix-fracture connections, and p mw i,j,k and W mw 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, p wf i,j,k = 0. 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 represent the vectors of unknown pressure difference and its coefficient, respectively, which are derived from fracture-fracture connections inside the same fracture, 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, p ff u,v and C ff u,v represent the vectors of unknown pressure difference and its coefficient, respectively, which are derived from fracture-fracture connections in different fractures, and 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, 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 where A represents the coefficient matrix (a large sparse matrix), ̄ prepresents the vector of unknown pressure difference between two iteration steps, which is assembled from an unknown matrix and the fracture pressure difference (̄ p m and ̄ p f ) in the EDFM, and b represents a constant vector.
For the model of Figure 1, we assumed constant-pressure well production and used Equation 20 to get (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 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 1. The EDFM uses equal time steps of 1 day. SAPHIR uses an adaptive step MPa size, which is a combination of arithmetic and geometric progression. Figure 2A,B shows the EDFM and SAPHIR meshes. In this scenario, the total number of grid cells is about 1900 in Figure 2A, whereas the number is more than 5000 in Figure 2B. 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. 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 3. 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 3. However, overall, the difference is small, so we can conclude that the EDFM is consistent with the SAPHIR software.

SENSITIVIT Y 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 1. 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 4A shows the relationship between daily and cumulative production and the number of fractures. Figure 4B 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 4C), 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.

| 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 5B), 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 halflength on productivity, as shown by the black dotted line in Figure 5A. The pressure profiles ( Figure 5C) 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.

| 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 6A,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. From the pressure profiles ( Figure 6C), 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 7A,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 7C), 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.

| 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 8A), 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 8B), 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 quasiradial flow is, and the higher the production will be.

| Stress sensitivity of matrix permeability
To investigate the effect of stress sensitivity of matrix permeability, we set sensitivity coefficients of 0 MPa −1 , 0.001 MPa −1 , 0.01 MPa −1 , 0.05 MPa −1 , 0.1 MPa −1 , and 0.2 MPa −1 . From the simulation results ( Figure 9A,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.

| 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  Figure 10A. To save calculation time, we set the time step to 2 days.
From Figure 10B, we can see that the production rate is positively correlated with fracture height. Figure 10C 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 1. The natural fracture parameters are listed in Table 2.

| 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 11A. From the simulation results ( Figure 11B), we can see that natural fractures can greatly increase oil production. According to the pressure profiles ( Figure 11C), 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.
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 12. 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 12).

| 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 13A.
We can see from the simulation results ( Figure 13B) that the connected natural fractures greatly improve oil production, compared with disconnected fractures ( Figure 11B). From the pressure profiles ( Figure 13C), 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.