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

^{1}State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, China^{2}Petroleum Systems Engineering, Faculty of Engineering and Applied Science, University of Regina, Regina, SK, Canada^{3}State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum (Beijing), Beijing, China^{4}Engineering Research Center of Gas Resource Development and Utilization of Ministry of Education, China University of Petroleum (Beijing), Beijing, China

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

With the increasing oil and gas consumption, more attentions have recently been directed to exploit the unconventional resources (Schmoker, 2002; Jia 2017; Yekeen et al., 2018; Gao et al., 2018; Chen et al., 2019). As an important and unconventional clean fossil fuel, tight oil is expected to play an increasingly important role in meeting the energy demand (Song et al., 2015; Yang et al., 2015a; Song and Yang, 2017; Wang M. et al., 2018; Amjed et al., 2019; Jiang et al., 2020a). Owing to the advancement of drilling horizontal wells with multistage hydraulic fractures (HFs), tight oil reservoirs have been exploited commercially though the flow behavior is quite different from those of the conventional reservoirs due to the extremely low permeability (Miller et al., 2008; Meyer and Bazan, 2011; Liu et al., 2018b; Liu et al., 2020a; Liu et al., 2020b, Liu et al., 2020c). To reduce operational costs and mitigate the negative environmental impact, multi-well pad schemes have proven to be more effective for exploiting tight oil/gas reservoirs with and without natural fractures (NFs) (Manchanda and Sharma, 2013; Meng et al., 2017; Xiao et al., 2018; Jiang et al., 2020b; Jiang et al., 2020c); however, multi-well pressure interference under multi-well pad schemes severely distorts the flow patterns. In addition to greatly affecting well production, such a distortion inevitably makes the transient pressure analysis more challenging while the stress-sensitivity effect and contribution of NFs are usually excluded (Liu et al., 2018a). Therefore, it is of fundamental and practical importance to accurately evaluate performance for a multi-well pad with multistage HFs and NFs in a tight formation by taking their stress-sensitive effects into account.

Considering its geological characteristics, fluid flow behavior in a tight formation with both NFs and HFs can be simulated with two main methods: the continuum model and the discrete fracture network (DFN) model. The most classical continuum model is the dual-porosity model proposed by Warren and Root (1963), assuming matrix and NFs as the main storage spaces and flow path. Then, the flow exchange between matrix subsystem and NF subsystem can be quantified (Kazemi 1969; Jalali and Ershaghi, 1987; Hassanzadeh et al., 2009; Wang, L et al., 2018). On the basis of the dual-porosity model, Wu et al. (2004) proposed a triple-porosity model to estimate the flow exchange among matrix, NFs, and HFs, which was subsequently improved with respect to algorithm convergence and computational expenses (Ozkan et al., 2009, 2011; Stalgorova and Mattar, 2012). Due to its simple assumptions and high computational efficiency, the continuum model has been widely used to describe fluid flow in tight oil/shale gas reservoirs.

Using the micro-seismic monitoring and surveillance techniques, it is found that HFs do not occupy the entire formation after being hydraulically fractured. In a DFN model, the HFs are usually described as high permeability channels by means of local grid refinement (Noorishad et al., 1982; Juanes et al., 2002; Rogers et al., 2010). As such, the DFN model is considered to be more rigorously to capture the details of flow behavior in complex fracture network, while it needs to be solved by either numerical difference or dimension reduction. Although the DFN method can be used to accurately describe the fluid flow, its computational expenses are too high (Cipolla et al., 2011).

Combining the advantages of the aforementioned two models, the latest development is to combine the dual-porosity model and the HF model (Zhao et al., 2014; Zhao et al., 2015; Chen et al., 2016; Jia et al., 2016, Jia et al., 2017; Morteza 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 secondary-fracture 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 Yang, 2018). Physically, more practical models together with computational algorithms are desired when the stress-sensitive effect needs to be taken into account (Liu et al., 2015; 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.

**FIGURE 1**. **(A)** 3D view and **(B)** Top view of multi-well pad with multiple fractures in a tight oil reservoir, and **(C)** Convergence flow effect around the wellbore in hydraulic fractures (HFs).

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, *γ*_{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 *γ*_{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 (Zhang et al., 2014),

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, *γ*_{hf} is the HF permeability modulus, and *k*_{hfmin} is the minimum HF permeability.

### Governing Equations

#### 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 *i*th 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 (Chen et al., 2016; 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:

#### Natural Fracture Subsystem

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:

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:

where

By integrating with the boundary condition, we can obtain

where

#### 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 2*N*_{1} × *M*_{1} + 2*N*_{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(2*N*_{1} × *M*_{1} + 2*N*_{2} × *M*_{2}) + 2]-order matrix from the pressure-continuity condition can be obtained.

Since pressure in NFs and HFs is the same (2*N*_{1} × *M*_{1} + 2*N*_{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, (2*N*_{1} × *M*_{1} + 2*N*_{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,

Therefore, a closed [2(2*N*_{1} × *M*_{1} + 2*N*_{2} × *M*_{2}) + 2]-order matrix is 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 transient-pressure 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 transient-production solution at a fixed pressure, which can be transformed through Eq. 24c in the Laplace domain.

where

### 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

**FIGURE 4**. Top view of commercial reservoir simulator simulator for the multi-well pad with multiple fractures.

**TABLE 1**. Physical properties of a multi-well pad with multiple fractures used in the commercial reservoir simulator.

**FIGURE 5**. Comparison between the results obtained from this work and numerical simulation from commercial reservoir simulator.

## 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, *γ*_{nfD} = 0.02, *γ*_{hfD} = 2, *ω* = 0.2, *λ* = 1 × 10^{−8}, *C*_{D} = 10, and *S* = 0.1. After considering the multi-well pressure interference and stress sensitivity, there are eight flow regimes, which can be detailed as follows:

(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 (Wang et al., 2017).

(8) The late-time pseudo-radial flow regime: Under the influence of matrix stress sensitivity, the DPD curve does not show a 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 *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, *γ*_{nfD} = 0.02, *γ*_{hfD} = 2, *ω* = 0.2, *λ* = 1 × 10^{−8}, *C*_{D} = 10, and *S* = 0.1. Comparison of different type curves illustrates that a smaller *ε* leads to a larger value of DPD and DPPD after the occurrence of multi-well pressure interference. 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 7**. Effect of well rate ratio on the transient pressure for multi-well pad with multiple fractures (MWPMF).

#### Effect of Hydraulic Fractures Permeability Modulus

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, *γ*_{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 *γ*_{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 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, *γ*_{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 *γ*_{nfD} is increased. For the case of *γ*_{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.

#### Effect of Hydraulic Fracture Spacings

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, *γ*_{nfD} = 0.02, *γ*_{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.

#### Effect of Hydraulic Fracture Length

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, *γ*_{nfD} = 0.02, *γ*_{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 (Chen et al., 2016; 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, *γ*_{nfD} = 0.02, *γ*_{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).

#### Effect of Minimum Fracture Conductivity

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, *γ*_{nfD} = 0.02, *γ*_{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.

#### Effect of Crossflow Coefficient

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, *γ*_{nfD} = 0.02, *γ*_{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 *λ* (Liu et al., 2015; Li et al., 2017).

#### Effect of Storage Ratio

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, *γ*_{nfD} = 0.02, *γ*_{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 well-rate 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.

## Data Availability Statement

The original contributions presented in the study are included in the article/supplementary materials, further inquiries can be directed to the corresponding author/s.

## Author Contributions

ZM: Conceptualization, Methodology, Formal analysis, Writing ‐ original draft; HL: Methodology, Validation; XT: Revision; GL: Resources, Supervision; LW: Validation, Formal analysis; DY: Conceptualization, Supervision, Funding acquisition, Resources, Writing ‐ review & editing.

## Funding

The authors are grateful to the State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation at the Southwest Petroleum University for their support of this work. The authors acknowledge a Natural Science Foundation of China (Grant No.: 51704246 and Grant No.: 51404282) for financial support. Also, the authors acknowledge a Discovery Grant and a Collaborative and Research Development (CRD) Grant awarded to D. Yang from the Natural Sciences and Engineering Research Council (NSERC) of Canada.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Glossary

*B* = volume formation factor, m^{3}/m^{3}*C*_{FD1} = dimensionless Well #1 fracture conductivity, fraction*C*_{FD2} = dimensionless Well #2 fracture conductivity, fraction*C*_{tnf} = total compressibility, MPa^{−1}*h* = reservoir thickness, m*h*_{hf} = hydraulic fracture height, m*k*_{hf} = hydraulic fracture permeability, D*k*_{nf} = natural fracture permeability, D*k*_{nfi} = initial natural fracture permeability, D*k*_{m} = matrix permeability, D*k*_{hfmin} = minimum hydraulic fracture permeability, D*L*_{H1} = horizontal length for Well #1, m*L*_{H2} = horizontal length for Well #2, m*L*_{f1} = Well #1 fracture half length, m*L*_{f2} = Well #2 fracture half length, m*M*_{1} = Well #1 fracture number, integer*M*_{2} = Well #2 fracture number, integer*N*_{1} = Well #1 discrete-fracture segments number, integer*N*_{2} = Well #2 discrete-fracture segments number, integer*p*_{nf} = natural fracture pressure, MPa*p*_{hf} = hydraulic fracture pressure, MPa*p*_{m} = matrix pressure, MPa*p*_{i} = reservoir initial pressure, MPa*q*_{Ni;1} = parameter defined in Eq. 8a*q*_{sc1} = Well #1 production rate, m^{3}/d*q*_{sc2} = Well #2 production rate, m^{3}/d*q*_{f} = flow-rate strength from reservoir into HFs, m^{3}/d*q*_{n} = flow-rate strength within HFs, m^{3}/d*r*_{w} = wellbore radius, m*s* = Laplace transformation variable*t* = time, hour*w*_{hf} = HF width, m*x*_{hf} = reference length, m*x*_{e} = reservoir length along *x* direction, m*y*_{e} = reservoir length along *y* direction, m*α* = shape factor, m^{−2}*ω* = storage ratio, fraction*λ* = interporosity coefficient, fraction*γ*_{nf} = NF permeability modulus, MPa^{−1}*γ*_{hf} = HF permeability modulus, MPa^{−1}*ε* = well-rate ratio, fraction*ρ* = fluid density, kg/m^{3}*µ* = viscosity, cp*ϕ*_{m} = matrix porosity, fraction

### Subscripts

*D*= dimensionless

*hf*= hydraulic fracture

*nf*= natural fracture

*m*= matrix

*i*= initial condition

### Superscript

- = Laplace transform## Appendix A: Derivation of Slab-Source Solution

A slab-source function can be obtained for the intersection of a slab-source in the *x*-direction and the *y*-direction (Zhang et al., 2014; 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,

Taking the Laplace transform, Eq. A1 can be written as follow,

where

where Eq. A2 can be rewritten as follow,

Since the slab source is assumed, the source-shape-dependent rate in terms of withdrawal rate is

The pressure-rate-source function in the Laplace domain can be obtained through Eqs A6 and A8,

According to the aforementioned definitions of dimensionless variables, Eqs A2–A7 can be rewritten as follows,

where

Then, the solution for multi-well pad schemes can be expressed as the follow:

## Appendix B: Extension to Multiple Wells

A closed [2(2*N*_{1} × *M*_{1} + 2*N*_{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(2*N*_{1} × *M*_{1} + 2*N*_{2} × *M*_{2})+2],

where *mn* means 2*N*_{1} × *M*_{1} + 2*N*_{2} × *M*_{2}, and

The submatrix *A*_{v,w} can be further written in Eq. B2 if *v* and *w* do not meet the previously discussed conditions:

The element of *a*_{i, j} in submatrix *A*_{v,w} is given as

where

The vectors of

The aforementioned equations show the coupled solution for the nonlinear systems of multi-well pad with two wells. It is worthwhile mentioning that one of the advantages of the newly developed model is that it can be extended to multiple wells.

Here, three horizontal wells are used as an example to demonstrate such an extension. A closed [2(2*N*_{1} × *M*_{1} + 2*N*_{2} × *M*_{2} + 2*N*_{3} × *M*_{3}) + 3]-order matrix from the pressure-continuity condition can then be obtained,

where matrix A is a coefficient matrix of dimension [2(2*N*_{1} × *M*_{1} + 2*N*_{2} × *M*_{2} + 2*N*_{3} × *M*_{3}) + 3],

where *xy* means 2*N*_{1} × *M*_{1} + 2*N*_{2} × *M*_{2} + 2*N*_{3} × *M*_{3}, and

The submatrix *A*_{v,w} can be further written in Eq. B2 if *v* and *w* do not meet the previously discussed conditions:

The element of *a*_{i, j} in submatrix *A*_{v,w} is given as

where

The vectors of

## References

Amjed, H., Mohamed, M., Abdulaziz, A., Mustafa, B., Salaheldin, E., Mohammed, B., et al. (2019). Gas condensate treatment: a critical review of materials, methods, field applications, and new solutions. *J. Petrol. Sci. Eng.* 177, 602–613. doi:10.1016/j.petrol.2019.02.089

Cao, L., Lu, L., Li, X., Wang, H., He, W., Xu, B., et al. (2018). Transient model analysis of gas flow behavior for a multi-fractured horizontal well incorporating stress-sensitive permeability. *J. Pet. Explor. Prod. Technol.* 9 (2), 855–867. doi:10.1007/s13202-018-0555-z

Chen, S., Li, H., and Zhang, Q. (2008). A New technique for production prediction in stress-sensitive reservoirs. *J. Can. Petrol. Technol.* 47 (3), 49–54. doi:10.2118/08-03-49

Chen, Z., Liao, X., Zhao, X., Lv, S., and Zhu, L. (2016). A semianalytical approach for obtaining type curves of multiple-fractured horizontal wells with secondary-fracture networks. *SPE J.* 21 (2), 538–549. doi:10.2118/178913-pa

Chen, Z., Liao, X., Zhao, X., Lyu, S., and Zhu, L. (2017a). A comprehensive productivity equation for multiple fractured vertical wells with non-linear effects under steady-state flow. *J. Petrol. Sci. Eng.* 149, 9–24. doi:10.1016/j.petrol.2016.09.050

Chen, Z., Liao, X., Zhao, X., Dou, X., Zhu, L., and Sanbo, L. (2017b). A finite-conductivity horizontal-well model for pressure-transient analysis in multiple-fractured horizontal wells. *SPE J.* 22 (4), 1112–1122. doi:10.2118/177230-pa

Chen, Y., Ma, G., Jin, Y., Wang, H., and Wang, Y. (2019). Productivity evaluation of unconventional reservoir development with three-dimensional fracture networks. *Fuel* 244, 304–313. doi:10.1016/j.fuel.2019.01.188

Cinco-Ley, H., and Meng, H. (1988). “Pressure transient analysis of wells with finite conductivity vertical fractures in double porosity reservoirs,” in SPE annual technical conference and exhibition, Houston, TX, October 2–5, 1988, Paper SPE-18172-MS.

Cinco-Ley, H., and Samaniego, -V., F. (1981). Transient pressure analysis for fractured wells. *J. Petrol. Technol.* 33 (9), 1749–1766. doi:10.2118/7490-pa

Cipolla, C. L., Lolon, E. P., and Erdle, J. C. (2010). Modeling well performance in shale-gas reservoirs. *SPE Reservoir Eval. Eng.* 13 (4), 638–653. doi:10.2118/125530-pa

Cipolla, C. L., Fitzpatrick, T., Williams, M. J., and Ganguly, U. K. (2011). “Seismic-to-simulation for unconventional reservoir development,” in SPE reservoir characterization and simulation conference and exhibition, Abu Dhabi, UAE, October 9–11, 2011, Paper SPE-146876-MS.

Feng, Q., Xia, T., Wang, S., and Singh, H. (2017). Pressure transient behavior of horizontal well with time-dependent fracture conductivity in tight oil reservoirs. *Geofluids* 2017, 5279792. doi:10.1155/2017/5279792

Gao, G., Zhang, W.-W., Ma, G.-F., Chen, G., Li, T., Hu, L.-Z., et al. (2018). Mineral composition and organic geochemistry of the lower cretaceous Xiagou formation source rock from the Qingxi Sag, Jiuquan Basin, Northwest China. *Petrol. Sci.* 15, 51–67. doi:10.1007/s12182-017-0213-y

Chen, H., Pooladi-Darvish, M., and Atabay, S. (2009). Shape factor in the drawdown solution for well testing of dual-porosity systems. *Adv. Water Resour.* 32 (11), 1652–1663. doi:10.1016/j.advwatres.2009.08.006

Jalali, Y., and Ershaghi, I. (1987). “Pressure transient analysis of heterogeneous naturally fractured reservoirs,” in SPE California regional meeting, Ventura, CA, April 8–10, 1987, Paper SPE-16341-MS.

Jia, P., Cheng, L., Huang, S., and Wu, Y. (2016). A semi-analytical model for the flow behavior of naturally fractured formations with multi-scale fracture networks. *J. Hydrol.* 537, 208–220. doi:10.1016/j.jhydrol.2016.03.022

Jia, P., Cheng, L., Huang, S., Xu, Z., Xue, Y., Cao, R., et al. (2017). A comprehensive model combining Laplace-transform finite-difference and boundary-element method for the flow behavior of a two-zone system with discrete fracture network. *J. Hydrol.* 551, 453–469. doi:10.1016/j.jhydrol.2017.06.022

Jia, C. (2017). Breakthrough and significance of unconventional oil and gas to classical petroleum geology theory. *Petrol. Explor. Dev.* 44 (1), 1–10. doi:10.1016/s1876-3804(17)30002-2

Jiang, L., Liu, T., and Yang, D. (2019a). Effect of stress-sensitive fracture conductivity on transient pressure behavior for a horizontal well with multistage fractures. *SPE J.* 24 (3), 1342–1363. doi:10.2118/194509-pa

Jiang, L., Liu, T., and Yang, D. (2019b). A semianalytical model for predicting transient pressure behavior of a hydraulically fractured horizontal well in a naturally fractured reservoir with non-darcy flow and stress-sensitive permeability effects. *SPE J.* 24 (3), 1322–1341. doi:10.2118/194501-pa

Jiang, L., Liu, J., Liu, T., and Yang, D. (2020a). A semianalytical model for transient pressure analysis of a horizontal well with non-uniform fracture geometry and shape-dependent conductivity in tight formations. *J. Petrol. Sci. Eng.* 195, 107860.

Jiang, L., Liu, J., Liu, T., and Yang, D. (2020b). Semi-analytical modeling of transient pressure behaviour for a fractured vertical well with hydraulic/natural fracture networks by considering stress-sensitive effect. *J. Nat. Gas Sci. Eng.* 82, 103477.

Jiang, L., Liu, J., Liu, T., and Yang, D. (2020c). Semi-analytical modeling of transient rate behaviour of a horizontal well with multistage fractures in tight formations considering stress-sensitivity effect. *J. Nat. Gas Sci. Eng.* 82, 103461.

Juanes, R., Samper, J., and Molinero, J. (2002). A general and efficient formulation of fractures and boundary conditions in the finite element method. *Int. J. Numer. Methods Eng.* 54, 1751–1774. doi:10.1002/nme.491

Kazemi, H. (1969). Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution. *SPE J.* 9 (4), 451–462. doi:10.2118/2156-a

Li, X., Cao, L., Luo, C., Zhang, B., Zhang, J., and Tan, X. (2017). Characteristics of transient production rate performance of horizontal well in fractured tight gas reservoirs with stress-sensitivity effect. *J. Petrol. Sci. Eng.* 158, 92–106. doi:10.1016/j.petrol.2017.08.041

Liu, M., Xiao, C., Wang, Y., Li, Z., Zhang, Y., Chen, S., et al. (2015). Sensitivity analysis of geometry for multi-stage fractured horizontal wells with consideration of finite-conductivity fractures in shale gas reservoirs. *J. Nat. Gas Sci. Eng.* 22, 182–195. doi:10.1016/j.jngse.2014.11.027

Liu, G., Meng, Z., Cui, Y., Wang, L., Liang, C., and Yang, S. (2018a). A semi-analytical methodology for multiwell productivity index of well-industry-production-scheme in tight oil reservoirs. *Energies* 11 (5), 1054. doi:10.3390/en11051054

Liu, G., Bai, Y., Gu, D., Lu, Y., and Yang, D. (2018b). Determination of static and dynamic characteristics of microscopic pore-throat structure in a tight oil-bearing sandstone formation. *AAPG (Am. Assoc. Pet. Geol.) Bull.* 102(9), 1867–1892. doi:10.1306/0108181613217061

Liu, G., Meng, Z., Luo, D., Wang, J., Gu, D., and Yang, D. (2020a). Experimental evaluation of interlayer interference during commingled production in a tight sandstone gas reservoir with multi-pressure systems. *Fuel* 262, 116557. doi:10.1016/j.fuel.2019.116557

Liu, G., Yin, H., Lan, Y., Fei, S., and Yang, D. (2020b). Experimental determination of dynamic pore-throat structure characteristics in a tight gas sandstone formation with consideration of effective stress. *Mar. Petrol. Geol.* 113, 104170. doi:10.1016/j.marpetgeo.2019.104170

Liu, J., Jiang, L., Liu, T., and Yang, D. (2020c). Modeling tracer flowback behaviour for a fractured vertical well in a tight formation by coupling fluid flow and geomechanical dynamics. *J. Nat. Gas Sci. Eng.* 84, 103656.

Manchanda, R., and Sharma, M. M. (2013). “Time-delayed fracturing: a new strategy in multi-stage, multi-well pad fracturing,” in the SPE annual technical conference and exhibition, New Orleans, LA, September 30-October 2, 2013, Paper SPE-166489-MS.

Meng, Z., Yang, S., Wang, L., Zou, J., Jiang, Y., Liang, C., et al. (2017). CO_{2} storage capacity for multi-well pads scheme in depleted shale gas reservoirs. *Energies* 10 (11), 1724. doi:10.3390/en10111724

Meyer, B. R., and Bazan, L. W. (2011). “A discrete fracture network model for hydraulically induced fractures-theory, parametric and case studies,” in SPE hydraulic fracturing technology conference, Woodlands, TX, January 24–26, 2011, Paper SPE-140514-MS.

Miller, B. A., Paneitz, J. M., Yakeley, S., and Evans, K. A. (2008). “Unlocking tight oil: selective multistage fracturing in the Bakken shale,” in SPE annual technical conference and exhibition, Denver, CA, September 21–24, 2008, Paper SPE-116105-MS.

Montgomery, C. T., and Steanson, R. E. (1985). Proppant selection: the key to successful fracture stimulation. *J. Petrol. Technol.* 37 (12), 2163–2172. doi:10.2118/12616-pa

Morteza, D., Hassan, H., and Chen, Z. (2018). Semi-analytical solution for pressure transient analysis of a hydraulically fractured vertical well in a bounded dual-porosity reservoir. *J. Hydrol.* 565, 289–301. doi:10.1016/j.jhydrol.2018.08.020

Noorishad, J., Ayatollahi, M. S., and Witherspoon, P. A. (1982). A finite-element method for coupled stress and fluid flow analysis in fractured rock masses. *Int. J. Rock Mech. Min. Sci. Geomech. Abstr.* 19, 185–193. doi:10.1016/0148-9062(82)90888-9

Ozkan, E., Brown, M., Raghavan, R., and Kazemi, H. (2009). “Comparison of fractured horizontal-well performance in conventional and unconventional reservoirs,” in SPE Western regional meeting, San Jose, CA, March 24–26, 2009, Paper SPE-121290-MS.

Ozkan, E., Brown, M. L., Raghavan, R., and Kazemi, H. (2011). Comparison of fractured-horizontal-well performance in tight sand and shale reservoirs. *SPE Reservoir Eval. Eng.* 14 (2), 248–259. doi:10.2118/121290-pa

Pedrosa, O. A. (1986). “Pressure transient response in stress-sensitive formations,” in SPE California regional meeting, Oakland, CA, April 2–4, 1986, Paper SPE-15115-MS.

Ren, Z., Wu, X., W., Liu, D., Rui, R., Wei, G., Chen, Z., et al. (2016). Semi-analytical model of the transient pressure behavior of complex fracture networks in tight oil reservoirs. *J. Nat. Gas Sci. Eng.* 35, 497–508. doi:10.1016/j.jngse.2016.09.006

Rogers, S., Elmo, D., Dunphy, R., and Bearinger, D. (2010). “Understanding hydraulic fracture geometry and interactions in the horn river basin through DFN and numerical modeling,” in SPE Canadian unconventional resources and international Petroleum conference, Calgary, AB, October 19–21, 2010, Paper SPE-137488-MS.

Schmoker, J. W. (2002). Resource-assessment perspectives for unconventional gas systems. *AAPG (Am. Assoc. Pet. Geol.) Bull.* 86 (11), 1993–1999. doi:10.1306/61EEDDDC-173E-11D7-8645000102C1865D

Song, C., and Yang, D. (2017). Experimental and numerical evaluation of CO_{2} huff-n-puff processes in Bakken formation. *Fuel* 190 (4), 145–162. doi:10.1016/j.fuel.2016.11.041

Song, Y., Li, Z., Jiang, L., and Hong, F. (2015). The concept and the accumulation characteristics of unconventional hydrocarbon resources. *Petrol. Sci.* 12 (4), 563–572. doi:10.1007/s12182-015-0060-7

Stalgorova, E., and Mattar, L. (2012). “Analytical model for history matching and forecasting production in multifrac composite systems,” in SPE Canadian unconventional resources conference, Calgary, AB, October 30-November 1, 2012, SPE 162516-MS.

Stehfest, H. (1970). Algorithm 368: numerical inversion of laplace transforms [D5]. *Commun. ACM.* 13 (1), 47–49. doi:10.1145/361953.361969

Van Everdingen, A. F., and Hurst, W. (1949). The application of the laplace transformation to flow problems in reservoirs. *J. Petrol. Technol.* 1 (12), 305–324. doi:10.2118/949305-g

Wang, H., Guo, J., and Zhang, L. (2017). A semi-analytical model for multilateral horizontal wells in low-permeability naturally fractured reservoirs. *J. Petrol. Sci. Eng.* 149, 564–578. doi:10.1016/j.petrol.2016.11.002

Wang, L., Yang, S., Meng, Z., Chen, Y., Qian, K., Han, W., et al. (2018). Time-dependent shape factors for fractured reservoir simulation: effect of stress sensitivity in matrix system. *J. Petrol. Sci. Eng.* 163, 556–569. doi:10.1016/j.petrol.2018.01.020

Wang, M., Chen, S., and Lin, M. (2018). Enhancing recovery and sensitivity studies in an unconventional tight gas condensate reservoir. *Petrol. Sci.* 15, 305318. doi:10.1007/s12182-018-0220-7

Warren, J. E., and Root, P. J. The behavior of naturally fractured reservoirs. *SPE J.* (1963). 3(3), 235–255. doi:10.2118/426-pa. CrossRef Full Text

Weaver, J. D., Rickman, R. D., and Luo, H. (2010). Fracture-conductivity loss caused by geochemical interactions between man-made proppants and formations. *SPE J.* 15 (1), 116–124. doi:10.2118/118174-pa

Wu, Y., Liu, H., and Bodvarsson, G. S. (2004). A triple-continuum approach for modeling flow and transport processes in fractured rock. *J. Contam. Hydrol.* 73, 145–179. doi:10.1016/j.jconhyd.2004.01.002

Xiao, C., Dai, Y., Tian, L., Lin, H., Zhang, Y., Yang, Y., et al. (2018). A semianalytical methodology for pressure-transient analysis of multiwell-pad-production scheme in shale gas reservoirs, Part 1: new insights into flow regimes and multiwell interference. *SPE J.* 23 (3), 0885–0905. doi:10.2118/187958-pa

Yang, D., Song, C., Zhang, J., Zhang, G., Ji, Y., and Gao, J. (2015a). Performance evaluation of injectivity for water-alternating-CO_{2} processes in tight oil formations. *Fuel* 139 (1), 292–300. doi:10.1016/j.fuel.2014.08.033

Yang, D., Zhang, F., Styles, J. A., and Gao, J. (2015b). Performance evaluation of a horizontal well with multiple fractures by use of a slab-source function. *SPE J.* 20 (3), 652–662. doi:10.2118/173184-pa

Yao, S., Zeng, F., and Liu, H. (2013). “A semi-analytical model for hydraulically fractured wells with stress-sensitive conductivities,” in SPE unconventional resources conference, Calgary, AB, November 57, 2013, Paper SPE-167230-MS.

Yekeen, N., Padmanabhan, E., and Idris, A. K. (2018). A review of recent advances in foam-based fracturing fluid application in unconventional reservoirs. *J. Ind. Eng. Chem.* 66, 45–71. doi:10.1016/j.jiec.2018.05.039

Zerzar, A., Tiab, D., and Bettam, Y. (2004). “Interpretation of multiple hydraulically fractured horizontal wells,” in SPE Abu Dhabi international conference and exhibition, Abu Dhabi, UAE, October 10–13, 2004, Paper SPE-88707-MS.

Zhang, F., and Yang, D. (2014). Determination of fracture conductivity in tight formations with non-Darcy flow behavior. *SPE J.* 19 (1), 34–44. doi:10.2118/162548-pa

Zhang, F., and Yang, D. (2018). Effects of non-Darcy flow and penetrating ratio on performance of horizontal wells with multiple fractures in a tight formation. *J. Energy Resour. Technol.* 140 (3), 032903. doi:10.1115/1.4037903

Zhang, Z., He, S., Liu, G., Guo, X., and Mo, S. (2014). Pressure buildup behavior of vertically fractured wells with stress-sensitive conductivity. *J. Petrol. Sci. Eng.* 122, 48–55. doi:10.1016/j.petrol.2014.05.006

Zhao, Y., Zhang, L., Luo, J., and Zhang, B. (2014). Performance of fractured horizontal well with stimulated reservoir volume in unconventional gas reservoir. *J. Hydrol.* 512, 447–456. doi:10.1016/j.jhydrol.2014.03.026

Keywords: transient pressure analysis, multi-well pad, stress-sensitivity, tight oil reservoir, slab source function

Citation: Meng Z, Lu H, Tan X, Liu G, Wang L and Yang D (2020) Effect of Stress-Sensitive Fracture Conductivity on Transient Pressure Behavior for a Multi-Well Pad With Multistage Fractures in a Naturally Fractured Tight Reservoir. *Front. Energy Res.* 8:600560. doi: 10.3389/fenrg.2020.600560

Received: 30 August 2020; Accepted: 16 October 2020;

Published: 19 November 2020.

Edited by:

Kaiqiang Zhang, Imperial College London, United KingdomReviewed by:

Yu Pang, University of Calgary, CanadaLu Wang, Chengdu University of Technology, China

Qian Wang, China University of Petroleum, China

Copyright © 2020 Meng, Lu, Tan, Liu, Wang and Yang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Guangfeng Liu, lgfmail@yeah.net; Daoyong Yang, Tony.Yang@uregina.ca

^{†}These authors have contributed equally to this work and share first authorship.