The Laplace-transform embedded discrete fracture model for flow simulation of stimulated reservoir volume

Hydraulic fracturing is a commonly adopted approach to enhance the production of unconventional reservoirs; however, the resulted complex fracture network increases the difficulty in the prediction of the flow behaviors. In this work, for the first time, we propose the Laplace-transform embedded discrete fracture model focusing on the simulation of fluid transport in simulated reservoir volume (SRV). The main hydraulic fracture is represented explicitly by the fracture segments referring to the recently developed simulation scheme, namely the embedded discrete fracture model. The equivalent model and multiple continua model are employed to characterize the fluid flow in SRV. The Laplace-transform finite-different method is used to numerically model the flow among unstimulated reservoir volume, SRV, and main hydraulic fracture with sufficient flexibility to characterize the arbitrary fracture/SRV geometries. As the solution is performed in Laplace domain, the backward Euler difference for time discretization is not necessary. Thus, the stability and convergence problems caused by time discretization are avoided. A series of numerical test cases are discussed to examine the performance of the proposed model. The simulation workflow is validated through the comparison with discrete fracture model and embedded discrete fracture model in real space. It is demonstrated that the proposed Laplace-transform embedded discrete fracture model is accurate for single flow with complex SRV distribution. On the basis of the Laplace-transform embedded discrete fracture model, we provided three applications of the new method: (1) analysis of the SRV effect on fluid flow behavior, (2) pressure transient analysis considering reservoir heterogeneity, (3) production performance analysis with time-varying production rate.


Introduction
Fluid flow in fractured geological strata should be accurately predicted in many subsurface applications, including groundwater management, hydrocarbon exploitation, carbon dioxide sequestration, and geothermal energy extraction (Yang et al. 2019;Dong et al. 2019). For unconventional oil/ gas development, hydraulic fracturing is the most effective technique to achieve economic production (Thomas et al. 2017;Vidic et al. 2013;Zhang et al. 2017). The appropriate scale separation for the modeling of fluid transport is hardly possible due to the complexity and heterogeneity of the fracture network (Ren et al. 2018). Accurate fracture representation continues to be a big challenge for oil/gas production performance analysis. As such, accurate and efficient methods are essential for the modeling of fluid flow in hydraulically fractured media.
Seismic monitoring, a commonly used technology to detect fracture networks, can only provide limited information of fractures. Mayerhofer et al. (2010) put forward the concept of "stimulated reservoir volume (SRV)" based on microseismic-fracture-mapping, which is an effective method in hydraulic fracture diagnostic. It describes the distribution of complex fracture network around the horizontal wells formed after hydraulic fracturing treatment. The well production is closely determined by the SRV properties including the size, shape, and conductivity. Figure 1 shows 1 3 the microseismic-fracture-mapping of well XY with 24 fracture stages. The SRV has a very complex 3D structure, and the distribution is irregular both in horizontal and vertical planes. The complexity of SRV is strongly influenced by the preexisting natural fractures and in situ stresses in the formation (Weng et al. 2014). In formations that contain a large number of natural fractures and with mild stress anisotropy, the SRV is inclined to be large. While for the formation with no natural fracture development but strong stress anisotropy, the induced fractures are always in planar shape.
It is a great technical challenge to characterize the stimulated reservoir volume (SRV) and unstimulated reservoir volume (USRV). One of the major difficulties is that hydraulic fracturing causes significant level of fracture complexity in SRV region. Effective models capable of simulating the flow in different regions should be established carefully. Currently, the existing methods suitable for SRV characterization fall mainly into two classes, namely the equivalent continuum model and discrete fracture model. The equivalent continuum approach has been extensively used by the well performance analysis community for its straightforwardness, and various extended approaches have been proposed based on a similar philosophy. The dual continuum model is a traditional approach of modeling naturally fractured rocks in which the matrix is the main media for fluid storage, and the natural fracture is the main media contributing to the fluid conductivity (Warren and Root 1963;Russell and Truitt 1964;Ramey 1973, 1974). The dual continuum model requires the existence of representative elementary volume (REV), i.e., it assumed the formation has a large number of natural fractures and the fractured domain can achieve a good scale separation. A trilinear flow model coupling linear flows in three sequential flow regions was proposed to analyze the SRV effect on well performance. The inner formation between adjacent hydraulic fractures can be regarded as the SRV in which the dual continuum flow model is used Ozkan et al. 2011). The analytical trilinear flow solution to simulate the well behaviors of SRV is computationally convenient and capable of providing the dominant flow regimes. As an extension of the trilinear flow model, the five-region model is proposed to enable the subdivision of the reservoir into five regions. Accordingly, the model can simulate the hydraulic fractures that are surrounded by a stimulated region of a limited extent, while incorporating the other regions that remain unstimulated (Stalgorova and Mattar 2013). The SRV in this case is taken as a media with a single-porosity type. The equivalent hydraulic properties in SRV can be upscaled from the spatial information and conductivity of both fracture and matrix. The demand for accurate characterization of the transient flow between matrix and fracture call for the fine division of matrix. Therefore, the multiple interacting continua (MINC) model becomes increasingly popular Younis 2015, 2016;Wang et al. 2017). The concept of MINC model is to divide a matrix REV into a series of nested subcells to represent the multiple interconnected continua distributed all over the formation domain (Wu and Pruess 1988;Pruess 1985).
Because of the advantages in computational efficiency and simplicity of implementation, the equivalent continuum model is widely adopted for the analysis of subsurface fluid flow through fractured rock system (March et Kuhlman et al. 2015). However, the continuum model has shown some limitations because the model assumptions are not always appropriate (Karimi-Fard et al. 2006;Zhang et al. 2018). In certain cases, it is necessary to characterize the distribution of fracture network in detail, and the discrete fracture model (DFM) demonstrates great advantages in this respect. DFM is a typical discontinuum method in which the hydraulic fractures or natural fractures in SRV are represented explicitly. A series of DFMs are successfully developed and some showed good adaptability. A prerequisite for the DFM simulation is to determine the spatial distribution and conductivity of the fractures. One way to obtain the information of the hydraulic fracture networks is through direct hydraulic fracturing stimulation (Weng et al. 2015;Li et al. 2016). The major difficulty of this approach is to collect the information of preexisting natural fractures and in situ stresses in the formation. The uncertainty associated with the preexisting natural fractures and in situ stresses significantly reduces the accuracy of the hydraulic fracturing simulation. Another way is to extract the fracture information from microseismic observations, which, however, is still insufficient to achieve an accurate representation of fracture distribution. In order to simplify the simulation workflow, the SRV is assumed to be a group of fractures perpendicular to each other (Mayerhofer et al. 2006). When performing the simulation, the PEBI meshing technique is capable of capturing the local small-scale characteristics (Sun et al., 2016a(Sun et al., , 2016b. Through mesh refinement, Voronoi cells are constructed and optimized to accommodate the exact fracture geometry. However, the PEBI meshing technique always leads to a huge grid number and high computational cost. In recent years, the embedded discrete fracture model (EDFM) has drawn more and more attention. EDFM is firstly proposed as a hierarchical modeling method to treat small-, medium-, and large-scale fractures separately (Lee et al. 1999(Lee et al. , 2001Li and Lee 2008). EDFM does not employ the sugar cube approximation for matrix as in dual continuum model or need fine unstructured grid discretization as in conventional DFM model. Besides, test cases indicate that EDFM is computationally efficient. The classical embedded discrete fracture model needs to address different types of connections that may occur between different cells, such as connection between a fracture cell and its neighboring matrix grid, connection between two intersecting fracture cells, and connection between two grid domains of an individual fracture. The calculation rules for connection transmissibility are well-documented in the existing literature. A holistic review of existing studies on EDFM is given in Table 1. EDFM approach has been employed for compositional simulator development, fractured vuggy reservoir modeling, fracture inversion, and upscaling purposes. As conventional EDFM can only achieve good accuracy for the simulation of high conductivity fracture. An improved approach, known as the projection-based embedded discrete fracture model (pEDFM) (Ţene et al., 2017;Jiang and Younis 2017;Ren et al. 2018), accurately simulates the low conductivity fractures by adding extra non-neighbor connections (NNCs). As an effective tool, the improved pEDFM can produce good accuracy for both low conductivity fracture case and multiple fluid saturation case. Considering the merits and drawbacks of continuum model and DFM approaches, the combination of both methods proved to be more effective (Zhang et al. , 2020. DFM is used to treat the main hydraulic fracture and the continuum model is adopted to model the SRV region. In a recent study, the large-scale fractures were dealt with EDFM method, while micro-fractures were represented by DP model. Considering the low permeability of matrix, transient transfer model between matrix and natural fracture was adopted (Xu et al. 2019).
In this paper, we present a new workflow to predict the evolution of pressure field of unconventional oil/gas reservoirs with stimulated reservoir volume. The workflow is based on a newly developed method, i.e., the Laplace-transform embedded discrete fracture model method. The numerical method is firstly proposed as the Laplace-transform finite difference (LTFD) (Moridis and Reddell 1991;Moridis et al. 1994). When using this method, the time domain is solved semi-analytically so the accuracy and stability of the simulation will not be affected by the back finite difference of time in real space. The space derivatives in the governing equations are still approximated by numerical finite difference method. This method allows for the consideration of reservoir heterogeneities and flow simulation with a 3D SRV configuration. The method initially can only treat vertical well, but significant improvement was made by combining the LTFD concept with a discrete fracture model, i.e., the embedded discrete fracture model. Three family of methods to characterize the SRV are proposed. A general option is to treat all possible fractures with EDFM model instead of using the hybrid model, which, however, can substantially increase the calculation time. Thus, the MINC model and equivalent model are used in the modeling of fluid flow in SRV region, while the EDFM is used to deal with the main hydraulic fracture. We present some numerical test results compared with DFM. We compare the solutions and some integrated quantities by choosing different solution methods. Our results show that Laplace-transform EDFM can provide good accuracy. The paper is organized as follows. We present the detailed formulations of the described problem in Sect. 2, while Sect. 3 elaborates on the numerical solution strategy. The numerical results obtained from different methods are shown in Sect. 4 based on synthetic test cases, after which the analysis will be focused on application of LTEDFM.

Physical model description
A properly designed physical model is necessary for the accurate simulation of the fractured horizontal well with SRV. In view of the complexity of a real producing reservoir, a wide range of scales should be considered for the fractures. Hence, detailed modeling of fractures remains unattainable at present. The SRV concept was shown to be a simple and effective alternative method. As shown in Fig. 2, four cases (2) shows a fractured horizontal well with SRV around each hydraulic fracture where Ω SRV is the SRV domain; (3) shows a fractured horizontal well with SRV contains all hydraulic fractures; (4) shows the complex fracture system around the well explicitly 1 3 are generally considered: Case-A shows a hydraulically fractured horizontal well with four main fractures, around which no SRV was created; Case-B shows a hydraulically fractured horizontal well with SRV around each hydraulic fracture, the volume and conductivity may be different among the four SRVs; Case-C shows a hydraulically fractured horizontal well with all main hydraulic fractures contained in a single SRV, which means the SRVs are linked with each other after hydraulic fracturing treatment; Case-D shows the complex fracture system around the well, and this case explicitly presents all hydraulic fractures of different scales. We will analyze all the above cases using the newly developed numerical simulation model.

Mathematical model
The microseismic monitoring and fracture propagation modeling can be used to understand the hydraulic fracture distribution. Once the fracture information is determined (though in most cases only partially determined), the theoretical flow model should be established to take the complex fracture system into consideration. Figure 3 shows a simple reservoir example. The distribution of fractures around the horizontal well is complex, and the whole domain can be divided into three regions: the hydraulic fracture domain Ω F , the SRV domain Ω SRV , and the USRV domain Ω USRV . In the following analysis, case (2) in Sect. 2.1 will be chosen as the basic physical model with which the mathematical model is established. Then, we will show how to extend the model to the other three cases. Some assumptions should be followed to restrict the solution: (1) The flow domain is two-dimensional or three-dimensional, which means the main hydraulic fracture and SRV can have a 3D structure; (2) Three flow regions are considered from USRV to SRV to hydraulic fractures one by one; (3) Isothermal condition is considered; (4) Flow is single phase; (5) The fluid is slightly compressible.

Flow in hydraulic fracture
As shown in Fig. 3, the flow equation describing the fluid transport in hydraulic fractures can be expressed as: where F is the fracture porosity; is fluid density; q W F denotes the source and sink terms from fracture; v F is the flow velocity in fracture; q fs−F is the fluid flux between main fracture and secondary fracture in SRV.

Flow in SRV region
The SRV region contains two kinds of media: secondary fractures and matrix. The flow equation in secondary fractures can be expressed as: where fs is the secondary fracture porosity; v fs is the flow velocity in secondary fracture; q sf −m is the fluid flux between secondary fracture and matrix.
The flow equation of matrix in SRV can be expressed as: where m is matrix porosity; v m is the flow velocity in matrix. (1) Schematic representation of the physical model: four perforation points distribute along the well and complex hydraulic fractures form after fracturing fluid is injected. The red lines present the complex hydraulic fractures with high conductivity. The yellow regions are stimulated reservoir volume (SRV) containing the secondary fracture and matrix. The outer region is the unstimulated reservoir volume (USRV) which is modeled as single-porosity media. Fluids first flow into the wellbore from the SRV and then from USRV to SRV 1 3

Flow in USRV region
The flow equation of matrix in USRV can be expressed as: where USRV is matrix porosity; v USRV is the flow velocity in matrix. Darcy Law is applied to calculate the flow velocity in both SRV and USRV regions: where is the viscosity of fluid; k is the permeability; p is the fluid pressure.
The outer boundary condition is set to satisfy the Neumann boundary condition: where n Γ is the outward unit normal vector.
The continuity condition of the flow rate at the interface between the SRV domain Ω SRV and USRV domain Ω USRV is considered which can be described as:

Transformation of the flow equations
In order to achieve the solution in the whole domain, we transform the above mathematical model into Laplace space.
Generally, the liquid density can be defined as: where 0 R is the inverse of formation volume factor. Here, Subsequently, a general form of flow equation will be derived for main hydraulic fractures, SRV region, and USRV region.

R F is calculated as the following equation:
where c f is the fluid compressibility, p 0 is the reference pressure.
The porosity is a weak function of pressure which can be described as: Expanding the first term of the left-hand side of Eq. (1): For slightly compressible liquid, R ≈ 1 and C TF is reduced into Eq. (13) where where s is the Laplace space parameter; ΔR F (0) is the initial distribution of ΔR F at t = 0.

The flow equation in SRV region
R m and R fs are calculated as the following equation: Again, the porosity is a weak function of pressure: where 0,m and 0,fs are the porosity at the reference pressure for matrix and secondary fracture ( p 0,m and p 0,fs ), respectively. c fm,m is the matrix rock compressibility, and c fm,fs is the secondary fracture rock compressibility.
Expanding the first term of the left-hand side of Eq. (2): Expansion of the first term of the left-hand side of Eq.
(3) yields: where For slightly compressible liquid, R ≈ 1 and C T is reduced into Eq. (30) and (31): Expanding the second term of the left-hand side of Eq. 2 ∇ ⋅ v m : Expanding the second term of the left-hand side of Eq. 3 ∇ ⋅ v m : Considering the compressibility of rock and liquid. The general nonlinear equations, Eq. (2)  Equations (40) and (41) are the general flow equations in SRV.

The flow equation in USRV region
R USRV is calculated using the following equation: Similarly, the porosity-pressure relationship follows: Equation (4) can then be rewritten as:

Coupling model among different media
(1) Fluid transport in SRV The classical definition of shape factor that assumes a pseudo-steady mass transfer between fracture and matrix is not applicable for tight reservoir rock, so transient transfer should be considered instead. The MINC model is a reasonable method to describe the unsteady mass transfer between secondary fracture and matrix. Figure 4 shows the coupling strategy. The main hydraulic fracture is modeled using EDFM because of the high conductivity. The coupling of secondary fracture and matrix is modeled using MINC model.
In the MINC model, each matrix cell is subdivided into a series of nested subcells. We define the flux across each pair connection (subcell-subcell connection) as: Introducing Eqs. (8,9) transforms Eq. (52) as: Reformulating Eq. (53) in terms of ΔR: Laplace transform to Eq. (54) gives: MINC elements Fig. 4 Coupling model of different media in SRV: the red line represents the main hydraulic fractures, the blue line shows the secondary fractures. After upscaling, the coupling of secondary fracture and matrix can be described by the MINC model 1 3 Without considering 0 , the connection can be defined as: (2) Main fracture and secondary fracture coupling EDFM was proposed by Li and Lee (2008). Three connections should be defined in the implementation of EDFM: the connection between a fracture cell and its neighboring matrix grid, connection between two intersecting fracture cells in the same grid, and connection between two grids of an individual fracture. When considering the secondary fracture, the first connections should be replaced by a fracture cell and its neighboring secondary fracture grid. The average normal distance from secondary fracture to fracture should be defined and computed.
For EDFM, the flux at each pair connection is calculated by: For LTEDFM, we define the flux between each connected pair as: Introducing Eqs. (8) and (9): Taking the Laplace transform of Eq. (59): Without considering 0 , the NNC flux term becomes: As indicated above, LTEDFM introduced different connection lists compared with EDFM. Figure 5 shows a simple example for the connection list both with and without considering SRV. The whole domain is discretized into 25 matrix grids for the situation without SRV. The " + " shape fracture is discretized into six cells by the matrix grid boundaries. The total number of connections is 51, including five main fracture-main fracture cell connections, six matrix-fracture connections, and 40 matrix-matrix grid connections. The connection lists are the same as the classical EDFM. The connection calculation of classical EDFM can be used for the simulation of case 1 and case 4. However, the connection types are quite different when SRV is considered. As shown in Fig. 5(b), the whole space domain is discretized into 25 grids, including nine SRV grids and 16 matrix grids. The main hydraulic fracture cells are connected to the SRV grids directly. The connection list includes five main fracture-main fracture cell connections, 45 MINC subcell-subcell connections, six main fracture-SRV grid connections, 16 matrix-matrix grid connections, 12 SRV-SRV grid connections, and 12 matrix-SRV connections. It should be noted that the SRV grid must use the effective permeability and porosity when computing the connection between secondary fracture and matrix. Namely, in this example, grid 7, 8, 9, 12, 13, 14, 17, 18, and 19 should use the effective properties. The above connection calculation method can be used for the simulation of case 2 and case 3.

Solution strategy
EDFM are commonly employed to solve for pressure field directly in real space. We extended it to the Laplace space in Sect. 2. The LTEDFM system of simultaneous equations is: where the term M yields after dealing with all connection lists and the accumulation term. Ψ store the unknowns.D is a known right-hand side matrix.
The equation system can be solved to obtain the pressure field in Laplace space. Then, the Stehfest numerical inversion algorithm is used to calculate the pressure in real time space. Figure 6 shows the solution workflow. Some procedures to execute are summarized, as follows: (1)

Comparison with the existing solution methods
In this section, the developed new model is validated by comparing the pressure distribution of the discrete fracture method (DFM) and the EDFM method in real space. The discrete fracture method (DFM) proposed by Karimi-Fard and Firoozabadi (2001) is used to obtain the benchmark solution. This model is proposed for two-phase flow simulation initially and can be readily modified into single-phase flow simulation. The discretization is based on the standard Galerkin finite element method; thus, the unstructured triangular mesh should be used. To realize DFM, the software COMSOL Multiphysics is chosen as the mesh discretization tool. The dimension of fracture is reduced by one. Grid refinement is performed near the well and fractures to better capture the pressure profile. The weak form equation of fracture flow is derived from Eq. (1) and set in COMSOL directly. This strategy solves the problem in real space and is taken as the reference solution. The EDFM method in real space refers to our previous work (Xu et al. 2019). To achieve a robust verification and capture the characteristic of flow in complex fracture system, four models are considered corresponding to the cases in Sect. 2.1: case-A, case-B, case-C, and case-D.

) Case-A: fractured horizontal well with planar fractures
In this section, we validate the LTEDFM for fractured horizontal well with planar fractures. Figure 7(1) shows the fractured horizontal well benchmark case under the condition of constant production rate with no flow outer boundary. The size of 2D geometry is 600 m × 600 m, and the formation thickness is 10 m. Four planar fractures are located along the horizontal well. There is no SRV around the main hydraulic fractures. Thus, the mass transfer happens between the main fractures and matrix directly. The main fracture width is 0.0001 m, the absolute permeability is 7.5 × 10 (−10) m 2 calculated according to Eq. (63). Figure 7(2) shows the grid system for DFM model. The space domain is discretized with 13,340 triangular meshes. Figure 7(3) shows the grid system for EDFM model and LTEDFM model. The space domain is divided into 3600 square grids. Each main hydraulic fracture is discretized into ten cells. Totally, 40 fracture cells are divided. The simulation parameters of the problem are listed in Table 2.
To validate the LTEDFM model, we compare the fracture and matrix pressure with those obtained from DFM and EDFM in real space. Figure 8 shows the pressure The pressure distribution is symmetric with respect to the center of the reservoir. All three simulation methods are able to capture the typical fracture effect that the pressure spreads fast around the fracture. To allow for qualitative analysis between different methods, we compare all grid pressures between two different simulation methods in Fig. 9. The curves show an excellent agreement between the simulation results of the LTEDFM method and other methods. Figure 9 also shows that most grid pressures are above 5 MPa after 400 days of production. The calculation time for LTEDFM and EDFM is compared. A total of 110 time steps are simulated, with 253 Newton's iterations for EDFM in real space. The calculation time is 88.75 s. However, when the proposed LTEDFM is imposed to solve the same problem, the calculation time is 4.72 s. The new solution protocol becomes Newton's iteration free and the main computational expense is the Stehfest inversion process. Hence, the computational speed would be significantly enhanced.
(2) Case-B: SRV around each fracture Figure 10(1) shows a fractured horizontal well benchmark case considering SRV. In this case, we discuss the problem of hydraulic fracture-SRV-matrix coupling. There is a rectangular SRV around each hydraulic fracture. The size of SRV is 40 m × 160 m. Figure 10(2) shows the grid system for DFM model. The space domain is discretized  with 12,024 triangular meshes. All secondary fractures are explicitly displayed. Figure 10(3) shows the grid system for EDFM model and LTEDFM model. It should be noted that the secondary fracture permeability and porosity need to be converted to the corresponding effective continuum properties when using EDFM and LTEDFM. After conversion, the effective fracture porosity is 0.000008, and the effective fracture permeability is 8.33 × 10 (−16) m 2 . Figure 10(4) shows the MINC elements. Each matrix element is subdivided into five continua. The simulation parameters of the problem are listed in Table 3. Figure 11 shows the pressure distribution of case-B for different simulation models. The simulation results manifest the fast spread of pressure in SRV because of the high conductivity comparing to the USRV region. As a result, the pressure drop in SRV is much larger than USRV regions. The drawdown for the medium two SRVs is the largest to maintain at a fixed production rate. We compare all grid pressures between two different simulation methods, as shown in Fig. 12. The results show that the LTEDFM method produces good accuracy when considering SRV around main hydraulic fractures. The difference in pressure distribution is noticeable by comparing Case-A and Case-B. Comparing Fig. 8 and Fig. 11, the pressure drop for Case-B is much smaller at 400 days around the well. For Case-B, all grid pressures are above 10 MPa. The main reason is that the SRV significantly decreases the flow resistance near the wellbore region.
Another method to simulate the flow in fracture and matrix system is the equivalent model. Different from the MINC, the equivalent model resembles the single-porosity model by upscaling the matrix and fracture properties mainly including porosity and permeability. We compute the equivalent parameters for each grid within SRV region, which works out to be 0.1 for equivalent SRV porosity and 8.43 × 10 (−16) m 2 for equivalent SRV permeability. The permeability and porosity remain unchanged for USRV. Figure 13 shows the pressure distribution for DFM, EDFM, and LTEDFM using the equivalent model. The pressure distribution is consistent with that of EDFM + MINC model. Thus, the equivalent model is a reasonable choice for the SRV simulation problem, which is why some classical models, such as the five-region flow model (Stalgorova and Mattar, 2013), regarded SRV region as a single media. We also compare all grid pressures between two different simulation methods in Fig. 14. The results show that the LTEDFM method exhibits good calculation accuracy for the equivalent model.
(3) Case 3-SRV around all fractures Figure 15(1) shows a fractured horizontal well benchmark case considering a large size SRV. In this case, we set a rectangular SRV in which all main hydraulic fractures are located. The size of SRV is 400 m × 200 m. Figure 15(2) shows the grid system for DFM model. The space domain is discretized into 18,336 triangular meshes. All secondary fractures are explicitly displayed. Figure 15(3) shows the grid system for EDFM model and LTEDFM model. Figure 15(4) shows the MINC elements. Each matrix element is subdivided into five continua. To demonstrate the SRV effect on pressure distribution, we have simulated the flow of three different models. Figure 16 shows the pressure distribution for DFM, EDFM, and LTEDFM considering a large size SRV. We compare all grid pressures between two different simulation methods as shown in Fig. 17. The results show that the LTEDFM method exhibits good calculation accuracy when considering SRV considering a large size SRV. The simulation results indicate that all grid pressures are greater than 14 MPa after 400 days. The pressure drop for Case-C is the smallest among the first three cases in 400 days. The SRV greatly enhances the fluid mobility, and thus the fluid pressure can be maintained at a high level. We also run the simulations using the equivalent model. The upscaling properties is the same as case-B. Figure 18 shows the pressure distribution for DFM, EDFM, and LTEDFM using equivalent model. We also compare all grid pressures between two different simulation methods in Fig. 19, in (4) Complex fracture system Figure 20(1) shows the schematics of a fractured horizontal well with complex fracture system. In this case, all hydraulic fractures are explicitly characterized. Totally, there are eight fracture and two perforation points located at (215, 295) m and (345, 355) m. Figure 20(2) shows the grid system for DFM model. The space domain is discretized with 5688 triangular meshes. Figure 20(3) shows the grid system for EDFM model and LTEDFM model. All fractures are discretized into 110 cells. Figure 21 shows the pressure distribution for DFM, EDFM, and LTEDFM considering complex  Figure 22 demonstrates the comparison of the pressure distribution obtained from different simulation methods, and a good agreement is observed between LTEDFM and other methods. The simulation results indicate that the fluid flow is strongly dependent on the fracture geometry, as the pressure wave firstly diffuses along the fracture and then spread to far field.

SRV effect on flow performance
All simulation above are under the condition with a constant well production rate. A very interesting topic is to analyze the SRV effect on wellbore pressure performance. Figure 23 shows the pressure drop of horizontal well for case-A, B, C, and D. As can be seen, case-A without considering SRV yields a rapid pressure drop compared with the other three models. The pressure drop obtained from case-C and case-D is slow because the SRV/complex fracture reduces the flow resistance around the well. Comparisons between the simulation results of case-B and case-C indicate that a larger SRV size contributes to a higher conductivity for fluid flow.

Pressure analysis considering reservoir heterogeneity
In this part, we design a 2D case considering heterogeneity of both SRV and USRV. The object-based modeling is used to generate the permeability. The azimuth angle is set to 0, and the correlation length is set to 1000 m in both the major and minor correlation direction. The mean and deviation of the petrophysical parameter for each individual facies are given in Table 4. The vertical permeability field is set to equal to horizontal permeability. Figure 24 shows three classes of sampling result. The first class is for case-A and case-D without considering SRV. The second class is for case-B and the third class is for case-C. Figure 25 shows the pressure drawdown versus time for different simulation cases. The figures show that the pressure drop of different cases diverge obviously from each other when considering permeability heterogeneity. All curves show a rapid pressure drop initially, but the speed of pressure drawdown gradually decreases with the elapse of  time. For case-A, the pressure is much lower than other three cases because of the lowest conductivity around the well. It indicates that a larger-scale SRV is beneficial to maintain the pressure level. In Fig. 26, the effect of heterogeneity on pressure performance is assessed by different sampling with and without considering SRV. The pressure performance is notably different for different sampling models. After 1,000,0000 s production, the pressure has not yet spread to the reservoir boundary because of the low conductivity of USRV. When SRV exists, the pressure can diffuse to a relative farther distance from the well. For the same case, the pressure distribution varies from each other for different permeability heterogeneity. Therefore, neglecting the reservoir heterogeneity may lead to inaccurate forecast of pressure transient behavior.

Production rate variable versus time
The LTEDFM method can also be used to analyze the well production performance when the well rate varies with time.
Case-A is taken as the basic simulation model, and five time intervals are set according to Table 5. The pressure distribution is presented in Fig. 27. For case 1, the pressure spreads slower because of the low production rate at initial time, while for case 3, the pressure spreads faster because of the high production rate at initial time. At 500 days, the reservoir boundary has undergone notable pressure drawdown. Production rate always changes with time for field practice, and the LTEDFM can simulate this situation and facilitate the history matching process.

3D SRV
In this subsection, the problem with a 3D SRV structure is presented. We extend the 2D model into a five-layer 3D model. Figure 28 shows the physical model of the considered reservoir with a size of 600 × 600 × 50 m. As shown, no SRV is created in the top and bottom layers. The SRV size is 40 m × 160 m × 10 m for layer 2 and layer 4, while the SRV size is 400 m × 200 m × 10 m for layer 3. The background grid size is 10 × 10 × 10 m. We set four main hydraulic fractures which have the same profile as test case-A-D. The four fractures are all vertical and fully penetrates the whole five layers. The well perforation points locate at the third layer. The total background grid number is 18000, while the total main hydraulic fracture cell number is 200. The main hydraulic fractures are discretized by two-dimensional planar cells in space but expanded to equi-dimensional cells during computation. The width of the main hydraulic fracture is 0.0001 m, and each fracture cell volume is 0.001 m 3 . The connection number between fracture cell and fracture cell belonging to different matrix grids is 180 in the horizontal direction and 160 in the vertical direction. The connection number between the background grid and main hydraulic fracture cell is 200. The production rate for each fracture is 0.5 m 3 /day. Figure 29 shows the pressure distribution of five layers at t = 10 days. There is a contrast between the pressure distribution of different layers and the pressure diffusion in layer 3 in the fastest among all layers. The reason is that layer 3 has the biggest SRV volume; thus, the flow resistance is significantly reduced. As a comparison, the pressure spreads very slowly in layer 1 and layer 5 as there is no SRV created in both layers. Figure 30 shows the pressure distribution of five layers at t = 2000 days. The pressure distribution of five layers tends to be consistent, which indicates that the SRV in layer 3 influences the pressure diffusion not only in horizontal direction, but also in the vertical direction.

Conclusion
We developed a novel numerical simulation model for fluid production from a complex fracture network. The mathematical model for flow in the main hydraulic fracture, SRV, and USRV is derived. The Laplace-transform embedded discrete fracture model is developed to solve for the pressure evolution. We applied the model to four kinds of hydraulic fracturing situations, namely hydraulically fractured horizontal well without SRV, hydraulically fractured horizontal well with SRV around each hydraulic fracture, hydraulically fractured horizontal well with an SRV containing all hydraulic fractures, and complex fracture system around the well. Some important conclusions are drawn as follows: (1) The EDFM method is extended to Laplace space to capture the main fracture geometry. By coupling the EDFM and MINC, the flow characteristics in SRV is precisely captured. The present solution for a series of test cases is validated by comparing the simulation results with other simulation methods, including the EDFM and DFM.
(2) The detailed simulation of the four considered cases shows that SRV has a strong impact on pressure dis- tribution. Generally, a larger SRV size leads to a lower flow resistance that is beneficial for fluid production.
(3) Both equivalent model and multiple continuum model can be used for SRV description. The equivalent model concept has been used in some forward modeling stud-ies as it can group the fracture and matrix system into a single-porosity media. The equivalent model can obtain an accurate prediction of pressure distribution as the multiple continuum model does. (4) The proposed workflow can be extended to more complex problems, including pressure analysis considering reservoir heterogeneity, variable production rate, and 3D SRV simulation, which reveal the strong practicality of the Laplace-transform embedded discrete fracture model.
Funding This study was supported by the National Natural Science Foundation of China (51904323, U1762216), Natural Science Foundation of Shandong Province (ZR2019BEE030), and the Fundamental Research Funds for the Central University under 18CX02029A.

Conflict of interest
There is no conflict of interest.
Ethical approval I testify on behalf of all co-authors that our article submitted to Journal of Petroleum Exploration and Production Technology has not been published elsewhere. All authors have been personally and actively involved in substantive work leading to the manuscript, and will hold themselves jointly and individually responsible for its content.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.