Effect of Stress-Sensitive Fracture Conductivity on Transient Pressure Behavior for a Multi-Well Pad With Multistage Fractures in a Naturally Fractured Tight Reservoir

This paper presents a semianalytical solution for evaluating transient pressure behavior of a multi-well pad with multistage fractures in a naturally fractured tight reservoir by considering the stress-sensitive effect imposed by both natural and hydraulic fractures. More specifically, the model pertaining to matrix/natural fractures is considered as a dual-porosity continuum, while its analytical flow model can be obtained by use of a slab-source function in the Laplace domain. The hydraulic fracture model is solved by discretizing each fracture into small segments to describe the flow behavior, while stress sensitivity in both the natural fracture (NF) subsystem and hydraulic fracture (HF) subsystem has been taken into account. To validate the newly developed semianalytical model, its solution has been obtained and compared with those of a commercial numerical simulator. By generating the type curves, there may exhibit eight flow regimes: pure wellbore storage, skin effect transition flow, linear flow regime within HFs, early radial flow, biradial flow, transition flow, pseudo-steady diffusion, and the late-time pseudo-radial flow. Furthermore, late-time flow regimes are found to be significantly distorted by the multi-well pressure interference. The smaller the well-rate ratio is, the more distorted the pressure and pressure derivative curves will be. In addition, well spacing and fracture length are found to dominate the flow behavior when multi-well pressure interference occurs. As the well spacing is decreased, the fracture length is increased, and thus occurrence of multi-well pressure interference is initiated earlier. Permeability moduli of NFs and HFs impose no impact on the multi-well pressure interference; however, it can distort flow regimes, leading to a severe distortion of pressure and pressure derivative curves. Similarly, the effect of HF permeability modulus on the flow in a hydraulic fracture, the minimum fracture conductivity is another key factor affecting the “hump” on the pressure curve. As the crossflow coefficient is increased, flow exchange between matrix and NFs is increased. With an increase in the storage ratio, flow exchange lasts longer and the second “dip” on the pressure curve becomes deeper.


INTRODUCTION
Combining the advantages of the aforementioned two models, the latest development is to combine the dualporosity model and the HF model (Zhao et al., 2014;Zhao et al., 2015;Chen et al., 2016;Jia et al., 2016Morteza et al., 2018). Such combined models can not only accurately describe fluid flow within the HFs at different angles, but also consider the contribution of the secondaryfracture networks (Ren et al., 2016;Chen et al., 2017b). In practice, either a line-source function or a slab-source function can be used to obtain the semianalytical solution for such complex flow systems under different conditions, though convergence difficulties are encountered and computational expenses are high (Zhang and Yang, 2014;Yang et al., 2015b;Wang et al., 2017;Liu et al., 2018a, Zhang and. Physically, more practical models together with computational algorithms are desired when the stress-sensitive effect needs to be taken into account Chen et al., 2016;Wang et al., 2017). Usually, the stress sensitivities of both NFs and HFs are assumed to be the same while this is found to be not reasonable in reality (Jiang et al., 2019b). In fact, it has been found from laboratory experiments that proppant has a great influence on stress sensitivity of HFs, i.e., the HF permeability modulus is larger than that of NFs (Montgomery and Steanson, 1985;Feng et al., 2017). So far, no attempts have been made to examine the stress sensitivity for NFs and HFs with different contributions in a multi-well pad with consideration of the multi-well pressure interference.
In this paper, a new semianalytical model has been developed and validated to examine the stress-sensitive fracture conductivity on transient pressure response for a multi-well pad in a naturally fractured tight reservoir. More specifically, matrix/NFs is considered as a dual-porosity continuum, which can be solved analytically by use of a slab-source function in the Laplace domain, while HFs can be solved by discretizing each fracture into small segments. Then, the perturbation technique together with superposition principle and the Stehfest algorithm has been employed to solve such a highly nonlinear equation matrix with consideration of different stress-sensitivity for NFs and HFs. In order to validate the semianalytical model, its solution has been compared with that of a commercial simulator. With the assistance of this new approach, not only can flow regimes of multi-well pad with multiple fractures be identified, but also sensitivity analysis has been made to illustrate how key parameters affect the transient pressure responses in such a complex flow system.

Theoretical Formulations
In this study, two hydraulically fractured horizontal wells are considered as a multi-well pad, producing at two different surface flow rate q sc1 and q sc2 see (Figure 1). It is worthwhile mentioning that this two-well model can be extended to multiple wells, though computational expenses will be correspondingly increased (Appendix B). For convenience, the tight reservoir including matrix subsystem, NF subsystem, and HF subsystem is assumed to be isotropic, while other main assumptions are listed as follows: • The reservoir upper and lower boundaries are impermeable and the lateral boundary is infinite. • Two hydraulic fracturing horizontal wells are parallel to the upper and lower boundaries and penetrate completely, assuming the fractures have finite conductivity. • Stress-sensitive effects in the NFs and HFs are considered, while their stress sensitivity coefficients are different. • A pseudo-steady crossflow occurs between the matrix and the NFs, while the HFs are the only pathway for the fluid flow from the NFs to the wellbore. • Oil flow in the HFs and NFs obeys the Darcy law, and incompressible flow is assumed in the HFs because the fracture volume is negligible compared with that of the NFs (Cinco-Ley and Meng, 1988;Jiang et al., 2019a). • Single phase and isothermal flow is considered but gravity and capillary effects are neglected.
With the aforementioned assumptions, we first discretize the HFs (Figure 2) each of which for Well #1 and Well #2 is respectively divided into 2*N 1 and 2*N 2 segments, assuming the number of HFs for Well #1 and Well #2 is M 1 and M 2 , respectively. The superposition principle can then be used to calculate the pressure responses for these sub-fracture segments. Finally, the dual-porosity flow model and HF flow model are dynamically coupled and solved. In this study, the stress-sensitive effects in the NFs and HFs are considered, while their stress-sensitivity coefficients are different.
The NF permeability modulus, c nf , is defined as the exponential relationship between NF permeability and its pressure (Chen et al., 2008), i.e.,   where k nf is the NF permeability, p nf is the NF pressure, p i is the initial pressure, k nfi is the initial NF permeability, and c nf is the NF permeability modulus. As for HFs, the HF permeability tends to have a minimum value (i.e., k hfmin ) due to the influence of proppants when the effective stress increases continuously (Weaver et al., 2010). The stress sensitivity of hydraulic fractures can be adjusted by changing the ratio of k hfmin to k hfi , which is usually obtained by experiments. The larger the ratio of k hfmin to k hfi is, the weaker the impact of effective stress on HF permeability will be. As such, the relationship between HF permeability and its pressure can be described as follow , where k hf is the HF permeability, p i is the initial reservoir pressure, p hf is the HF pressure, k hfi is the initial HF permeability, c hf is the HF permeability modulus, and k hfmin is the minimum HF permeability.

Matrix Subsystem
Since the matrix permeability of a tight oil reservoir is extremely low, the fluid flow in the matrix can be neglected, and the matrix is considered as evenly distributed sources when the oil in the matrix flows into natural fractures (Li et al., 2017). In this way, the governing equation of the matrix subsystem can be described by: where ρ is the fluid density, ϕ m is the matrix porosity, k m is the matrix permeability, p m is the matrix pressure, μ is the fluid viscosity, α is the shape factor, t is time, and q ex is the flow rate between the matrix and NFs. Because of no flow in the matrix, the initial condition can be described as the following:

Natural Fracture Subsystem
Considering the stress sensitivity of NFs, the governing equation of NFs can be expressed as follows: where c tnf is the total NF compressibility. Then, the initial condition and boundary conditions can be expressed as follow:

Hydraulic Fracture Subsystem
To date, as for fluid flow in the HFs, both compressible and incompressible flows have been extensively studied. Since the reservoir volume is much larger than the hydraulic fracturing volume in a tight oil reservoir (Cinco-Ley and Samaniego, 1981), fluid compressibility is not considered in this study. Considering the stress sensitivity of HFs, the governing equation of HFs can be written as follows: where w hf and h hf are the HF width and height, respectively, µ is fluid viscosity and q f is the flow rate from reservoir into HFs. The initial condition and boundary conditions are expressed as follow: where q Ni;1 is the flow rate from the N node of the ith HF to the wellbore and M is the number of hydraulic fractures.

Non-Dimensionalization
For the convenience of derivation, analysis, and improving the applicability of the model proposed in this study, dimensionless variables are defined as follows Wang et al., 2017;Jiang et al., 2019a;Jiang et al., 2019b):

Matrix Subsystem
According to the aforementioned definitions of dimensionless variables, Eqs 3a-4b can be rewritten as follows: With the aforementioned definitions, Eqs 5-6b can be rewritten as follows:

Hydraulic Fracture Subsystem
Similarly, Eqs 7-8c can be rewritten as follows:

Solutions for a Multi-Well Pad
Matrix/Natural Fractures.
After considering the NF stress sensitivity, it is difficult to obtain analytical solutions for the mathematical models of matrix and NFs. In this paper, we use the perturbation theory (Pedrosa, 1986) to obtain semianalytical solutions of matrix and NFs, i.e., p nfD can be rephrased as follows: According to perturbation theory, η nfD can be rewritten as follow: Numerous studies have found that the zero-order approximate solution completely meets the solution requirements (Jiang et al., 2019b). Therefore, Eqs 5-b can be rewritten as follows: (14c) And the initial condition can be rewritten as follows: Then, the equation can be solved in the Laplace domain. In this paper, a slab-source function is derived and presented in Appendix A, while the solution for multi-well pad schemes can be expressed as the follow:

Hydraulic Fractures
Considering the HF stress sensitivity, the corresponding flow equation becomes highly nonlinear and is solved by discretizing each fractures into small segments. According to the dimensionless definitions, Eqs 12e-13a can be rewritten as follows: Therefore, the governing equation for segments can be rewritten: with the boundary condition and initial condition: By integrating with the boundary condition, we can obtain

Coupling Matrix/Natural Fractures and Hydraulic Fractures
The solution of transient flow responses is obtained by coupling the Matrix/NFs and HFs dynamically. According to the previous model description, HFs are discretized into 2N 1 × M 1 + 2N 2 × M 2 segments, each of which has three variables: q fDi,j , q NDi,j , and p hfDi,j . In addition, two unknown variables are the flowing bottomhole pressure p wD1 and p wD2 . A closed [2(2N 1 × M 1 + 2N 2 × M 2 ) + 2]-order matrix from the pressure-continuity condition can be obtained. Since pressure in NFs and HFs is the same (2N 1 × M 1 + 2N 2 × M 2 ) equations can then be obtained, Because it is assumed that fluid in HFs incompressible, inflow and outflow should be conserved according to the conservation of mass, (2N 1 × M 1 + 2N 2 × M 2 ) equations can be obtained, Either Well #1 or Well #2 is produced at a constant rate, and ratio of the oil production for these two wells is denoted by ε. Then, two equations can be obtained, In addition, skin effect and wellbore storage can be introduced into transient-solution directly in the Laplace domain (Cao et al., 2018),

Coupling Solutions
An iterative method proposed by Van Everdingen and Hurst (1949) is applied here to solve the nonlinear matrix, while the detailed flowchart is shown as Figure 3. First, the transientpressure solution (Eq. 24a) can be obtained with iterations in the Laplace domain, and then stable solution (Eq. 24b) obtained in the real-time domain by using the Stehfest inversion algorithm (Stehfest, 1970). It is worthwhile mentioning that the new semianalytical model can be used to obtain the transientproduction solution at a fixed pressure, which can be transformed through Eq. 24c in the Laplace domain.

Model Validation
In this section, we used a commercial reservoir simulator (CMG, version 2015) to simulate the performance of a multi-well pad with multiple fractures. In order to simulate the performance of an infinite reservoir, we intentionally enlarged the reservoir scale. The whole grid system is 300 × 400 × 5, and the dimension of each grid is 20 m × 20 m × 2 m. A "double Klinkenberg permeability, logarithmic spacing and local refined grid (DK-LS-LRG) method" (Cipolla et al., 2010) is used to characterize the HFs. The top view of the multi-well pad with multiple fractures is shown in Figure 4, and the physical properties used for model validation are tabulated in Table 1. As can be seen in Figure 5, there is a good agreement between the numerical simulation and the results obtained from the newly developed model in this study, confirming that our approach is reliable.

RESULTS AND DISCUSSION
In this section, type curves of a multi-well pad with multiple fractures are obtained, and the flow regimes are identified. In addition, sensitivity analysis is conducted to examine the effects of some main parameters, including ratio of well rate, hydraulic fracture spacing, hydraulic fracture length, well spacing, permeability modulus of NFs and HFs, minimum fracture conductivity, storage ratio, and crossflow coefficient.

Identification of Flow Regime
Dimensionless pressure (DP) and the dimensionless pressure derivative (DPD) of a multi-well pad with multiple fractures (MWPMF) are shown in Figure 6. The following dimensionless parameters are used to generate the type curves: L f1D L f2D 75, D f1 D f2 500, M 1 M 2 4, N 1 N 2 5, C fDi 10, C fDmin /C fDi 0.4, c nfD 0.02, c hfD 2, ω 0.2, λ 1 × 10 −8 , C D 10, and S 0.1. After considering the multi-well pressure interference and stress (1) The pure wellbore storage period regime: It can be identified by a slope of unity on the log-log plots, when DP curve and DPD curve align with each other. This period is mainly dominated by wellbore storage effect (Liu et al., 2018a). (2) The skin effect transition flow regime: The most obvious feature of this stage is a "hump" in the pressure derivative curve. This period is mainly dominated by fluid properties (Cao et al., 2018). (3) The linear flow regime within HFs: Usually, this period is characterized by a slope of 1/2 on the pressure derivative curve; however, the new model is affected by dynamic conductivity, and the slope of the pressure derivative curve is larger than 1/2. Moreover, the variation rate of the conductivity of two wells is different, which results in that the DP curve and DPD curve of two wells do not overlap at this stage (Yao et al., 2013). (4) The early radial flow regime: It is characterized by a slope of 0 on the pressure derivative curve. This period is mainly dominated by hydraulic fracture spacing. This regime will occur only when the hydraulic fracture spacing is appropriate (Chen et al., 2016. (5) The biradial flow regime: It is characterized by a slope of 1/3 on the pressure derivative curve. This period is mainly dominated by wellbore length (Zerzar et al., 2004). (6) The transition flow regime: This is the transition regime between the biradial flow and pseudo-steady diffusion flow (Cinco-Ley and Meng, 1988). (7) The pseudo-steady diffusion regime: The most obvious feature of this regime is a "dip" in the pressure derivative curve. In this period, mass transfer occurs between matrix and NFs .  horizontal line, and its value is greater than 0.5. Because of the existence of multi-well pressure interference, the smaller the oil rate is, the greater the distorted value will be (Liu et al., 2018a).

Sensitivity Analysis
Effect of Well-Rate Ratio Figure 7 demonstrates the effect of well-rate ratio on type curves of MWPMF. The basic data used to generate the type curves are FIGURE 5 | Comparison between the results obtained from this work and numerical simulation from commercial reservoir simulator.
FIGURE 6 | Flow regimes of a multi-well pad with multiple fractures identified by the type curve.
In addition, it can be concluded that the well-rate ratio ε does not affect the onset time of multi-well pressure interference (Liu et al., 2018a). This is mainly because the onset time of multi-well pressure interference depends on the velocity of pressure wave propagation, which is only related to reservoir properties and well  parameters, mainly including matrix, fracture permeability, fracture half-length, wellbore length, and well spacing, but has nothing to do with the well-rate ratio. Figure 8 illustrates the effect of HF permeability modulus on type curves of MWPMF. Basic data used to generate the type curves are L f1D L f2D 75, D f1D D f2D 500, M 1 M 2 4, N 1 N 2 5, C fDi 10, C fDmin /C fDi 0.4, c nfD 0.02, ω 0.2, λ 1 × 10 −8 , ε 0.2, C D 10, and S 0.1. As illustrated in Figure 8, the HF permeability modulus is found to be one of the key factors affecting the "hump," which only affects the flow in a hydraulic fracture. The slope of the hump becomes larger as the HF permeability modulus is increased (Yao et al., 2013). Certainly, the greater the HF permeability modulus is, the easier it approaches C fDmin , and it is easier to move to the next regime (i.e., early radial flow regime). For an extreme case c hfD 0, the early radial flow regime will be masked by the linear flow regime within HFs. This is because C fD maintains a constant value and the linear flow lasts longer, and then it is more difficult for the early radial flow to occur.

Effect of Hydraulic Fractures Permeability Modulus
Effect of Natural Fracture Permeability Modulus Figure 9 displays the effect of NF permeability modulus on type curves of MWPMF. Basic data used to generate the type curves are L f1D L f2D 75, D f1D D f2D 500, M 1 M 2 4, N 1 N 2 5, C fDi 10, C fDmin /C fDi 0.4, c hfD 2, ω 0.2, λ 1 × 10 −8 , ε 0.2, C D 10, and S 0.1. As can be seen in Fig. 9, a larger NF permeability modulus leads to a larger value of DPD and DPPD after the occurrence of multi-well pressure interference (Liu et al., 2018a;Xiao et al., 2018). Besides, the NF permeability modulus can distort flow regimes, and the distortion of pressure curves become severer as the c nfD is increased. For the case of c nfD 0, the DPPD value in the late-time pseudo-radial flow regime is equal to 0.5. In the early regime of flow, the flow regime is independent of the NF permeability modulus before the multi-well pressure interference occurs. Figure 10 presents the effect of HF spacing on type curves of MWPMF. Basic data used to generate the type curves are L f1D L f2D 75, M 1 M 2 4, N 1 N 2 5, C fDi 10, C fDmin /C fDi 0.4, c nfD 0.02, c hfD 2, ω 0.2, λ 1 × 10 −8 , ε 0.2, C D 10, and S 0.1. Comparison of different type curves illustrates that the HF spacing affects the occurrence time and duration of the biradial flow regime, but has no effect on the distortion of the convection pattern. DPD and DPPD curves only shift from left to right or vice versa (the curve slope remains unchanged) (Liu et al., 2018a). A smaller D fD leads to the earlier appearance of the biradial flow regime. After the disappearance of this flow regime, the pressure drop and pressure derivative curves overlap again under different well spacings. Figure 11 plots the effect of HF length on type curves of MWPMF. Basic data used to generate the type curves are D f1D D f2D 500, M 1 M 2 4, N 1 N 2 5, C fDi 10, C fDmin /C fDi 0.4, c nfD 0.02, c hfD 2, ω 0.2, λ 1 × 10 −8 , ε 0.2, C D 10, and S 0.1. It is shown that HF length only affects the early flow regimes, including the skin effect transition flow regime, the linear flow regime within HFs, and the early radial flow regime. As the HF length is decreased, pressure drop and pressure derivative curves are shifted upwards 2Chen et al., 2017a;Cao et al., 2018). With the decrease of HF length, the duration of linear flow regime within HFs becomes short. Similar to the effect of HF spacing, after the completion of the early radial flow regime, the pressure drop and pressure derivative curves overlap again under different HF lengths.  Effect of Well Spacing Figure 12 depicts the effect of well spacing on type curves of MWPMF. Basic data used to generate the type curves are L f1D L f2D 75, D f1D D f2D 500, M 1 M 2 4, N 1 N 2 5, C fDi 10, C fDmin /C fDi 0.4, c nfD 0.02, c hfD 2, ω 0.2, λ 1 × 10 −8 , ε 0.2, C D 10, and S 0.1. By observing whether the pressure derivative curves overlap, we can judge when the multi-well pressure interference occurs. As shown in Figure 12, well spacing affects flow regimes after multi-well pressure interference starts. Comparison of different type curves illustrates that a smaller D wD leads to the earlier occurrence of multi-well pressure interference (Liu et al., 2018a;Xiao et al., 2018). It is worthwhile mentioning that well spacing also does not distort the shape of pressure drop and pressure derivative curves, and the  curve is only shifted from left to right or vice versa (i.e., the curve slope remains unchanged). Figure 13 illustrates the effect of minimum fracture conductivity on type curves of MWPMF. Basic data used to generate the type curves are L f1D L f2D 75, D f1D D f2D 500, M 1 M 2 4, N 1 N 2 5, C fDi 10, c nfD 0.02, c hfD 2, ω 0.2, λ 1 × 10 −8 , ε 0.2, C D 10, and S 0.1. The higher the value of C fDmin /C fDi is, the stronger the effect of proppant in HFs will be, and the higher the fracture conductivity can be maintained when the effective pressure is increased. Similar to the effect of HF permeability modulus on the flow in a hydraulic fracture, C fDmin /C fDi is another key factor affecting the "hump," whose slope becomes smaller as the minimum fracture conductivity is decreased. As can be seen in Figure 13, when the value of C fDmin /C fDi reaches 0.6, the hump basically disappears. This is because a higher C fDmin /C fDi means that the hydraulic fracture loses a small part of conductivity, and it is difficult to distinguish the difference among type curves. Figure 14 reveals the effect of crossflow coefficient on type curves of MWPMF. Basic data used to generate the type curves are L f1D L f2D 75, D f1D D f2D 500, M 1 M 2 4, N 1 N 2 5, C fDi 10, C fDmin /C fDi 0.4, c nfD 0.02, c hfD 2, ω 0.2, λ 1 × 10 −8 , ε 0.2, C D 10, and S 0.1. As described by Eq. 14a, with a larger λ, the magnitude of the difference between the NF permeability and the matrix permeability becomes smaller. This means that the ability of NFs to maintain production is weakened. Comparison of different type curves illustrates that the larger the λ is, the flow transfers from the matrix to the NFs easier and the "dip" comes earlier. Especially for an oil well with low production, it is more unlikely to have crossflow at a smaller λ Li et al., 2017). Figure 15 shows the effect of storage ratio on type curves of MWPMF. Basic data used to generate the type curves are L f1D L f2D 75, D f1D D f2D 500, M 1 M 2 4, N 1 N 2 5, C fDi 10, C fDmin /C fDi 0.4, c nfD 0.02, c hfD 2, λ 1 × 10 −8 , ε 0.2, C D 10, and S 0.1. As shown in Figure 15, storage ratio ω mainly determinates the duration and the depth of the second "dip" on the pressure derivative curve, while it also has a significant effect on the occurrence of multi-well pressure interference (Li et al., 2017;Jiang et al., 2019a). Comparison of different type curves illustrates that a smaller ω results in the earlier occurrence of the multi-well pressure interference with an increase in the duration of flow exchange. The slope of the second dip becomes larger with a decrease in the storage ratio.

CONCLUSIONS
In this paper, a new semianalytical model has been developed and validated to describe the flow behavior of a multi-well pad with multistage fractures in a naturally fractured reservoir, while various stress-sensitive effects in the NFs and HFs are examined. The main conclusions are summarized as follows: (1) Type curves of a multi-well pad with multiple fractures are identified, and there may exhibit eight flow regimes: pure wellbore storage, skin effect transition flow, linear flow regime within HFs, early radial flow, biradial flow, transition flow, pseudo-steady diffusion, and the late-time pseudo-radial flow. In addition, changing the location of two wells in a multi-well pad may mask some flow regimes.
(2) Under the multi-well pad schemes, multi-well pressure interference is inevitable and will significantly distort the late-time flow regimes. It is found that the smaller the wellrate ratio is, the more distorted the curve will be. As the well spacing is decreased, the fracture length is increased, and thus the multi-well pressure interference occurs earlier.
(3) Stress-sensitive effects in the NFs indicate the permeability damage which occurs at intermediate times, resulting in an increasing pressure drop. Stress-sensitive effects in the HFs and minimum fracture conductivity only affect the flow regime in hydraulic fractures. The sensitivity analysis shows that the "hump" becomes stronger as the minimum fracture conductivity is decreased. A slab-source function can be obtained for the intersection of a slab-source in the x-direction and the y-direction Jiang et al., 2019a). Then, the instantaneous source function for a multi-well pad with stress-sensitive effect in both NFs and HFs can be expressed as follows, S x, y, x w , y w , x f , y f , t VII x, x w , x f , t VII y, y w , y f , t

APPENDIX B: EXTENSION TO MULTIPLE WELLS
A closed [2(2N 1 × M 1 + 2N 2 × M 2 ) + 2]-order matrix from the pressure-continuity condition can be obtained as, where matrix A is a coefficient matrix of dimension [2(2N 1 × M 1 + 2N 2 × M 2 )+2], X → is the unknown vector, and B → is the known vector. The matrix and vectors are further expanded as the submatrix, Frontiers in Energy Research | www.frontiersin.org November 2020 | Volume 8 | Article 600560