© L. Tian et al., published by IFP Energies nouvelles, 2018

Introduction

Developing unconventional resource has challenged conventional methodologies and engineers for the past few decades. Recently, advances in technology and new analysis methods have made Single Multi-Fractured Horizontal Wells (SMFHW) to become more economical. One of the major challenges is the determination of optimal stage and well spacing in a specific drilling area. Rate transient analysis is commonly used to assess the effectiveness of fracture stimulation of production wells, which is used to determine additional optimization of well/stage spacing. Beside rate transient analysis method, pressure-transient analysis of horizontal well with multi-stage fractures is of importance as well in reservoir engineering because it provides valuable information concerning well completion, well placement, well spacing, calculation of dynamic reserve, and most importantly, different flow regimes characterized by pressure-transient behavior provide estimations of the in-situ reservoir properties and fracture geometry and conductivity.

The presence of complex fracture topology has significant impacts on establishment of pressure-transient model for unconventional reservoirs (Khvoenkova and Delorme, 2011; Baroni et al., 2015; Noetinger, 2015). Originally, some analytical approaches have been used to model the transient flow behavior in such systems. With assuming uniform distribution of identical HF along horizontal well, Ozkan and Raghavan (1991) and Ozkan et al. (2011) utilized the concept of Lee and Brockenbrough (1986) of tri-linear model with inner reservoir of naturally fractured to represent the SMFHW performance in unconventional reservoirs. Brown et al. (2011) presented an analytical tri-linear flow model that incorporates transient fluid transfer from matrix to fracture to simulate the pressure transient and production behavior of fractured horizontal wells in unconventional reservoirs. However, those proposed analytical models are unable to explicitly represent the Hydraulic Fracture Networks (HFNW) system induced by interaction between Main Hydraulic Fractures (MHF) and Secondary Hydraulic Fractures (SHF). To fill this gap, Jones et al. (2013), Mirzaei and Cipolla, 2012, Cipolla and Wallace, 2014, Farah and Ding (2016) used an unstructured-grid simulator to analyze the type curves of HFNW system. This latest unstructured-grid technology can make the simulation of fracture complexity more accurate by refining the vicinity of the high conductivity fractures network with fine grids. However, it is inevitable to increase the complication and economical consumption of computation. Subsequently, several ingenious semi-analytical methods have been proposed to overcome the shortage of in-accurate analytical methods and time-consuming numerical methods. Zhou et al. (2014) and Yu et al. (2015) combined analytical reservoir solutions with numerical fracture network solution to characterize fracture network complexity. Similarly, Jia et al. (2015) and Chen et al. (2016) utilized star-transformation (Karimi-Fard et al., 2003) and boundary element integration methods to characterize fluid seepage within complex fracture network (induced hydraulic fractures and discrete natural fractures). These semi-analytical models provided quick insights into fracture-network performance and formed foundations to efficiently and accurately analyze transient pressure response.

Recently, Multi-Well-Pad-Production (MWPP) scheme has been in the center of attention as a promising technology to economically improve SG recovery (Awada et al., 2015; Guindon, 2015). Micro-seismic fracturing mapping shows hydraulic fractures extending between wells, gaining the existence of Inter-Well Pressure Interference (IWPI) (Farley and Hutchinson, 2014; Sardinha et al., 2014). Although recognizing IWPI within MWPP scheme can provide valuable insights to gain a better understanding of the fracture design (Soroush et al., 2013; Kaviani et al., 2010), enhancement of the possibilities of IWPI in MWPP scheme severely distorts flow regimes, which increases the burden of parameter estimations (Awada et al., 2015). Consequently, this technology strongly challenges the traditional pressure-transient analysis methods mentioned above, which focus on SMFHW without pressure interference. Therefore, a methodology to identify pressure-transient response of MWPP scheme without and with IWPI is of significance for reservoir engineers.

In this work, our objectives are to obtain better understanding of IWPI in MWPP scheme by addressing the following questions:

• What kinds of flow regimes for MWPP scheme can we obtain? And what is the difference between MWPP scheme and SMFHW in flow regimes?

• How to identify IWPI for MWPP scheme based on obtained flow regimes?

To answer questions outlined above, following three issues should be solved sequentially and systematically:

• Developing an efficient semi-analytical mathematical model for MWPP scheme;

• Based on proposed semi-analytical pressure-transient model, establishing methodology to identify different flow regimes in MWPP scheme;

• Based on proposed semi-analytical pressure-transient model, analyzing influences of key parameters on flow regimes.

Described below are the attributes of our methodology framework. Section 2 describes the conceptual model of HFNW system in MWPP scheme and development of semi-analytical pressure-transient model in detail. Section 3 systematically implemented model validation, methodology to identify the IWPI and propose some new insights into flow regimes in MWPP scheme. Finally, Section 4 summarizes our contribution and promising work in the future. More information on model derivation can be found in Appendix section.

1 Development of Semi-Analytical Pressure-Transient Model

Single phase gas is assumed to derive the semi-analytical model. We envisage two distinct flow regimes governed by different physics: SG reservoir flow and HFNW system flow. In the following sections, we first describe the conceptual model of SG formation and MWPP scheme used for presentation of our approach and then mathematically model these two flow processes, respectively. Finally, SG reservoir flow model and HFNW system flow model are coupled dynamically.

1.1 Conceptual Model

1.1.1 Multi-Well-Pad-Production (MWPP) Scheme

Figure 1a illustrates the placement of MWPP schemes in Shale Gas (SG) reservoir. We mainly focus on single MWPP scheme at the early-intermediate production period. Each MWPP scheme contains several SMFHW. Figures 1b and c illustrate the micro-seismic surveillance within one MWPP scheme. As Figure 1b, the hydraulic fractures system can be idealized to several regular fractures. As Figure 1c, the complex hydraulic fracture network system has been formed. Therefore, a general HFNW is idealized to characterize both of these situations (as shown in Fig. 1d), and the corresponding interference can be classified into following two types (Awada et al., 2015):

• Interference through HFNW. Interference directly through connected HFNW refers to communication in SG reservoir when HFNW connection is created between two wells, as illustrated in Figure 1d(2);

• Interference through reservoir. Interference through the SG reservoir would occur when HFNW are not directly connected between wells but are in close proximity. In Figure 1d, when the fracture conductivity of green zone tends to be zero, it shows the fracture configuration where communication through SG reservoir rock may be observed as illustrated in Figure 1d(1).

1.1.2 SG Reservoir Model

The conceptual model of a SG formation can be described as follows: The reservoir, within infinite boundary, is treated as a 2-D flat. SG reservoir is assumed to be isotropic, including Natural Fractures (NF) system and matrix system, and bounded by upper and lower impermeable strata. Three media exist in the SG reservoir: (1) the lowest permeable shale matrix, (2) the moderate permeable NF, and (3) the highest permeable HFNW connecting to wellbore.

To conveniently describe our methodology, we chose two wells to be our research objective. Two wells produced at constant gas rate q₁ and q₂. Fluid in SG reservoir flows into HFNW at varying flow-rate strength q_f along fractures. In addition, some other assumptions are made as follows:

• Two horizontal wells are intercepted by HFNW. The HFNW are assumed to fully penetrate SG reservoir;

• The HFNW of two wells has different fracture conductivity;

• SG reservoir has uniform thickness h. The initial pressure is P_i; and the initial temperature is T;

• Gas seepage within NF system meets Darcy's law. Gas unsteady-state diffusion within matrix is assumed to obey Second Fick's law. NF system is stress-dependent with initial permeability k_ri;

• Compressibility coefficient of the slightly compressible SG is constant;

• Impacts of gravity and capillary pressure are neglected;

• Gas absorption and adsorption meets Langmuir isotherm equation;

• Wellbore storage and skin factor are considered;

• No frictional pressure loss inside the wellbore is considered.

1.2 Mathematical Model

1.2.1 Hydraulic Fracture Networks Discrete

As Figure 2, we further classify the HFNW into MHF and SHF. The properties of MHF and SHF in HFNW for well1 include: permeability, k_mf1, k_sf1; fracture width, w_mf1, w_sf1; width of HFNW, W_f1, half-length of HFNW, L_f1. The properties of MHF and SHF in HFNW for well2 include: permeability, k_mf2, k_sf2; fracture width, w_mf2, w_sf2; width of HFNW, W_f2, half-length of HFNW, L_f2. The distance between two wells is L_w. The vertical distance of HFNW for two wells is L_vf12. The horizontal distance of HFNW for two wells is L_hf12. To establish a mathematical model, we first subdivide the HFNW systems. HFNW of well1 is divided into N₁ sub-fracture segments, NFNW of well2 is divided into N₂ sub-fracture segments. The hydraulic fracture number of well1 is M₁. The hydraulic fracture number of well2 is M₂. The length of sub-fracture segment of well1 and well2 can be presented as ΔL_f1, ΔL_f2, respectively. We can summarize fracture properties as follows:

• Well1: N₁, M₁, L_f1, ΔL_f1, W_f1, k_mf1, k_sf1, w_mf1, w_sf1

• Well2: N₂, M₂, L_f2, ΔL_f2, W_f2, k_mf2, k_sf2, w_mf2, w_sf2

1.2.2 Model of SG Flow in SG Reservoir System

To develop mathematical models in SG reservoir, the formula of NF system and matrix system can be established separately and then dynamically coupled (Tian et al., 2014, 2016; Wang, 2014; Liu et al., 2015). By applying the principle of integration, the pressure distribution of a random position (x_D, y_D) caused by one fracture segment (x_WD, y_WD) is given as follows (more information on dimensionless definition and model derivation can be found in Appendix A and Appendix B), (1) where

As one can seen from Figure 2, the HFNW system has been divided sub-fracture segments. According to Equation (1), we can obtain transient pressure response at the center of each segment of HFNW system in MWPP scheme by the superposition principle: (2)  where (3) (4)

is the pressure response at the oth fracture segment, caused by the flux of bth fracture segment in ath hydraulic fracture for well1; is the pressure response at the oth fracture segment, caused by the flux of jth fracture segment in ith hydraulic fracture for well2.

1.2.3 Model of SG Flow in Hydraulic Fracture Networks System

Currently, modeling of fluid flow in HFNW system mainly focus on two issues: flow state (compressible and incompressible flow) within independent fracture segment and fluid transfer at connecting point of two crossing fracture segments. For former issue about flow state, Jia et al. (2015), Zeng et al. (2012) analyzed the unsteady state flow with considering fluid compressibility, while Chen et al. (2016) and Zhou et al. (2014) analyzed the pseudo-steady state flow without considering fluid compressibility. In our paper, we establish a generic model by considering the fluid compressibility.

For the later issue about fluid transfer at interacting point of HFNW system. Zhou et al. (2014) and Chen et al. (2016) artificially determined the flow direction, while Jia et al. (2015) employed star-transformation to automatically determine flow direction. In this paper, the orthogonal hydraulic fractures system is idealized to characterize the geometry of HFNW, therefore, the star-transformation is utilized to automatically determine flow direction and solve the issue of fluid transfer at interacting point of HFNW system.

Independent Fracture Segment. Here, semi-analytical method (Zeng et al., 2012; Zhou et al., 2014; Chen et al., 2016) is applied to develop the model of independent fracture segment considering finite hydraulic fracture conductivity. Zeng et al. (2012) proposed a semi-analytical method which divided fracture system into several segments and each segment was still solved analytically. As one can see in Figure 3, we can chose ith fracture segment to analyze the flow equation. Fluid flow from position l_D1 to l_D2, because of the fluid supplement of segment flux q_fi, rate increases from q_ci1 to q_ci2 along l.

We take ith fracture segment as an example, the solutions of other fracture segment are similar. Deriving from the semi-analytical method proposed by Zeng et al. (2012), the governing equation of the ith independent fracture segment in Laplace space can be established as follows: (5)

Boundary conditions at ϵ_i1 and ϵ_i1 can be given by (6) (7)

Combining with Equations (5) and (6), we can obtain the pressure distribution within the ith.

Independent Fracture Segment. The solution in Laplace space for ith fracture segment can be mathematically characterized by flow rate q_ci1, q_ci2 at both sides of the ith fracture segment and the fluid influx q_fi from SG reservoir into this fracture segment. (8)

_{Therefore, pressure at the center of the ith fracture segment} (9)  _where (10) (11) (12)

Connecting Fracture Segment. Here, the Star-Delta transformation is adopted to solve the fluid transfer at the connecting fracture segment. Figure 4 depicts the transformation for four interconnected fracture segments. Karimi-Fard et al. (2003) used this transformation to eliminate intermediate control volume in discrete fracture network simulation. Taking the Star-Delta transformation shown in Figure 4 and eliminate intersection of cell 0, we make the four segments connecting directly. The transmissibility between two adjacent segments can be written as, (13) where, T_Di,0 is the dimensionless transmissibility between ith fracture segment and intersection 0. (14)

Fracture segment 4 in Figure 4 is taken as an example to illustrate the application of Star-Delta transformation. The flow equation of segment 4 can be given by (15)

Furthermore, Equation (15) can change as the following form after the transformation (16)

1.3 Solution of Transient-Pressure for MWPP Scheme

Considering the proposed equations, there are three unknowns, , for each fracture segment, additional two unknowns are bottom-hole pressure of two wells. Therefore, the total number of unknowns is equal to [3(M₁ × N₁ + M₂ × N₂) + 2]. A closed [3(M₁ × N₁ + M₂ × N₂) + 2]-order matrix from the following conditions can be obtained.

Combing Equations (2) and (9), pressure continuity between HFNW system and SG reservoir can be satisfied at the center of ith fracture segment, (17)

After rewriting Equation (17) at each fracture segment, we can obtain totally (M₁ × N₁ + M₂ × N₂)-equations.

For two connected fracture segments, ith segment and (i+1)th segment, the pressure and fluid rate at connecting point respectively are equal to each other, (18) (19)

After rewriting Equations (18) and (19) for each fracture segment, we can obtain another totally 2(M₁ × N₁ + M₂ × N₂)-equations.

Another two equations are required to form a closed matrix. Well1 and well2 are producing at constant rate q₁ and q₂, respectively. Here, we define a new variable ε, represents the ratio between q₁ and (q₁ + q₂), namely, ϵ = q₁/(q₁ + q₂). Then, (20) (21)

Finally, a closed [3(M₁ × N₁ + M₂ × N₂) + 2]-order matrix is formed. By applying Gauss elimination and Stehfest numerical algorithm (Stehfest, 1970), the bottom-hole pressure solution, and , can be solved. In Laplace domain, the wellbore storage effects can be easily added into the solution with Duhamel's theorem. For the definition of dimensionless pseudo-pressure (Eq. A1), the solution can be formulated as follows: (22) (23) where, CD1, CD2 are dimensionless wellbore storage coefficient for well1 and well2 which are defined as

By the Stehfest numerical invention algorithm (Stehfest, 1970), the solution in real space can be obtained. One can get the bottom-hole pressure for MWPP scheme in SG reservoir by taking the stress sensitivity of NF into consideration. (24) (25)

2 Results and Discussion

In this section, four aspects will be systematically analyzed: (1) influence of fracture discretization level on pressure response; (2) model validation by comparison between semi-analytical model and fully numerical simulation; (3) identification of flow regimes based on special pressure-transient characteristics; and (4) sensitivity analysis of pressure response to fracture parameters. The relevant parameters are shown in Table 1.

2.1 Discretization Level of HFNW System

Due to the inappropriate discretization of HFNW system, early transient pressure response can be distorted by artifacts of fracture subdivision (Chen et al., 2016). Our approach is partially dominated by the discretization level of HFNW system. The basic dimensionless parameters are as follows: C_1D = C_2D = 0, ζ_D = 0, λ = 0.002, γ = 0.15, ω = 0.0035, C_ηD = 10⁵, L_f1D = 250, L_f2D = 250, W_f1D = 500, W_f2D = 500, L_wD = 5000, L_vf12D = 4500, L_hf12D = 0, C_f1D = 125, C_f2D = 125. In Figure 5, ψ_wD represents the Dimensionless Pseudo-Pressure (DPP), dψ_wD/d(ln t_D) represents the Dimensionless Pseudo-Pressure Derivative (DPPD). As Figure 5, when the number of fracture segments, M, is larger than 12 for each hydraulic fracture, the calculating results do not change appreciably. Therefore, each fracture is divided into 12 segments in our work, including for the model validation and sensitivity analysis.

2.2 Model Validation

To our best knowledge, so far, no semi-analytical pressure-transient models for multi-wells have been developed. Therefore, a fully numerical simulation, modeled by a commercial simulator CMG-GEM module, is utilized to validate our proposed model for MWPP scheme. The top view is shown in Figure 6. Three representative HFNW system are modeled: Figure 6a models the non-connection between HFNW system of the two wells, Figure 6b models the direct connection between HFNW system of the two wells. Figure 6c models the regular HFNW with transverse MHF and without SHF. The input data are listed in Table 1. Figure 6d illustrates one representative segment to represent one part of the reservoir volume around a hydraulic fracture. Four MHF are orthogonal to the horizontal well at 250 m fracture spacing. Other SHF are orthogonal to main fractures. The model is a 2-D model with 100 grid cells in the x-direction, 50 grid cells in y-direction and only one grid cell in the z-direction. The multi-porosity model and Multiple INteracting Continua (MINC) method are applied to subdivide the matrix so that the transient diffusion in matrix can be simulated. Desorption phenomenon is characterized to be instant desorption model. DK-LS-LGR technology is employed to characterize HFNW system. Each hydraulic fracture is represented by a 3 × 3 locally refined grid with 0.025-m-wide. Then, the pressure solution is calculated under a constant production rate. We assume that ratio between q₁ and q₂ is 1:4, and the fracture properties of two wells are consistent. After that, the numerical solutions are compared with semi-analytical results calculated in this paper. As we can see in Figure 7, there is a good agreement between our results and CMG's, which indicates that our model is reliable.

2.3 Identification of Flow Regime

The main goal of our research is identifying flow regimes of MWPP scheme in SG reservoir. As the numerical validation section, we model three different kinds of HFNW system. For the previous two HFNW system, currently, it is still not very common to comprehensively describe the flow regimes for complex HFNW system although several contributions have been made. Chen et al. (2016) and Jia et al. (2015) separately developed mathematical models to identify the flow regimes in complex Hens system. Besides well-known classic flow regimes (linear-flow, bi-linear flow regimes, etc.), the authors also added some new flow regimes by themselves. For example, “fluid feed” flow regime is induced by fluid transfer between MHF and SHF. Pseudo Boundary Dominated Flow (PBDF) is induced the permeability contrast between ultra-low permeable SG reservoir and high-permeable Hens system. Figure 8 shows the pressure distribution within MWPI at different production time. As Figure 8, it is apparently to observe that the IWPI indeed occurs for these three kinds of Hens system. To obtain some practical and generic analysis, the flow regimes of Hens with transverse MHF and without SHF are identified in our work as Figure 6c.

We can set some common parameters for those two wells: M₁ = 4, M₂ = 4, ζ_D = 0.05, λ = 0.002, γ = 0.15, ω = 0.0035, C_ηD = 10⁵. Setting consistent MHF properties (conductivity and half-length) L_f1D = L_fD2 = 200, C_f1D = C_f2D = 10. Another two parameters, L_hf12D and L_vf12D, determine the different interference type illustrated as Figure 9, which can be mathematically described as follows:

• Figure 9(1): L_vf12D = L_f1D + L_fD2 + 100 = 500, and L_hf12D = 0;

• Figure 9(2): L_vf12D = L_f1D + L_fD2 + 100 = 500, and L_hf12D = 500.

DPP and the DPPD of WIPS scheme are shown in Figure 10. To clearly describe the type curves, we will compare the characteristics of pressure curves between WIPS and SMFHW. By comparing between SMFHWs and WIPS scheme, we can add some additional information and better explanation into these distorted flow regimes, which can be described in detail as follows:

Regime I: The pure wellbore storage period regime. DP curve and DPD curve align, and the slope of curves are equal to 1. This stage is mainly controlled by wellbore storage effect and difficult to be impacted by the MWPI. Thus, the type curve of WIPS and SMFHW overlap with each other.

Regime II: The transition flow regime. The early stage of this regime gradually derives from the straight line which has unit slope. This stage is mainly controlled by fluid properties and also difficult to be impacted by MWPI. Thus, the type curve of WIPS and SMFHW also overlap with each other.

Regime III: The linear flow regime within HFs. This stage is mainly dominated by fracture conductivity. At this linear flow regime, and also difficult to be impacted by MWPI. Thus, the type curve of WIPS and SMFHW also overlap with each other.

Regime IV: The bi-linear flow regime. This stage is mainly controlled by fracture length. At this bi-linear flow regime, we can start to detect the MWPI for type Figure 1(b1). The slope of the DPD curves is actually bigger than 0.25. The distortion degree of pressure curve for small gas rate is also more significant than that of big gas rate. But, it is still difficult to be impacted by the MWPI for type Figure 1(b2), Thus, the type curve of WIPS and SMFHWs still overlap with each other.

Regime V: The early pseudo-radial flow regime. This stage is mainly dominated by fracture spacing. At this pseudo-radial flow regime, we can start to detect the MWPI for type Figure 1(b2). The slope of the DPPD curves is actually bigger than 0. The distortion degree of pressure curve for small gas rate is more significant than that of big gas rate.

Regime VI: The intermediate-time linear flow regime. This stage is mainly dominated by wellbore length. At this intermediate-time linear flow regime, we also can detect the MWPI for type Figure 1(b1) and Figure 1(b2). The slope of the DPD curves is actually bigger than 0.5. The distortion degree of pressure curve for small gas rate is more significant than that of big gas rate.

Regime VII: The late-time pseudo-radial flow regime. The shape of DPD curve is a horizontal line. The value of this horizontal well is equal to 0.5.

2.4 Sensitivity Analysis

In this section, we conduct some sensitivity analysis on pressure-transient response for MWPP scheme in SG reservoir. The key factors that influence the transient pressure response for MWPP scheme include hydraulic fracture half-length L_fD, hydraulic fracture conductivity C_fD, hydraulic fracture spacing L_f12D, well spacing L_wD, ratio of well rate ε, stress sensitivity ζ_D. Here, we just take the Figure 2d(2) and traditional definition of DPP as an example. Our analysis is also on basis of assumption that fracture properties of two wells are consistent. Other situations can be analyzed similarly. Some dimensionless parameters can be: S = 0.2,C_D = 10, ζ_D = 0.05, λ = 0.002, γ = 0.15, ω = 0.0035, C_ηD = 10⁵, L_f1D = 2000, L_f2D = 2000, L_wD = 3000, L_f12D = 1000. C_f1D = 50,C_f1D = 50. The results are discussed in detail as follows:

Ratio of well rate, ε. We set ε to be 0.1, 0.3, 0.5 respectively. Figure 11 shows effects of ratio of well rate on pressure performance for MWPP scheme. We can judge the occurrence of MWPI by whether pressure curves of two wells overlap together. MWPI starts form first radial-flow regimes. Subsequently, first radial-flow regime and second linear-flow regime are distorted severely. When the shape of pressure curves are distorted by the MWPI, the smaller the well rate, the more severely the pressure curves are distorted. Moreover, the smaller the well rate, the bigger the DPP and DPPD. Therefore, we can judge the well rate of two wells based on the relative position of the DPP curves and DPPD curves. We also can summarize that the ratio of well rate ε approximately has no any influence on the time when the MWPI occurs. This phenomenon can be strictly explained as follows: the time when the MWPI occurs is mainly relied on pressure wave propagation within reservoir and fracture system, which is dominated by fracture and rock properties, such as rock compressibility, fracture conductivity, half-length and well spacing. Thus, if we want to adjust the occurrence of MWPI, changing the production rate is useless.

Well spacing, L_wD. We set L_wD to be 450, 650, 900 respectively, and we also set ε to be 1:4. Figure 12 illustrates the impacts of well spacing L_wD on pressure performance for MWPP scheme. Similarly, we also can judge the occurrence of MWPI by whether pressure curves of two wells overlap together. On condition of different well spacing L_wD, MWPI basically starts form first radial-flow regimes. As the well spacing L_wD increases, the occurrence of MWPI becomes later. Besides, pressure curves will be split at part of flow regimes (such as first radial-flow regime and second linear-flow regime), the pressure curves will overlap again subsequently. We also can clearly observe another phenomenon that well spacing almost does not distort the shape of pressure curves, the pressure curves just move upward or downward (the slope of pressure curves keeps constant). When the gas rate is big enough, the impacts of well spacing on multi-well interference can hardly be identified (the dot line). Therefore, if we can find that the pressure curves are not distorted apparently, we can make an original judgment that the well rate is relative big.

Hydraulic Fracture spacing, L_f12D. We set L_f12D to be 100, 500, 1000 respectively, and we also set ε to be 1:4. Figure 13 illustrates the impacts of hydraulic fracture spacing L_f12D on pressure performance for MWPP scheme. Similar to the effects of well spacing on pressure curves, hydraulic fracture spacing mainly impacts the occurrence of MWPI and has no any influence on the distortion of flow regimes, the pressure curves just move upward or downward (the slope of pressure curves keeps constant). Different from the effects of well spacing on pressure curves, the impact of fracture spacing is more significant than that of well spacing. When L_f12D = 100, MWPI basically starts from first linear-flow regime when L_f12D = 1000, MWPI basically starts from first radial-flow regimes. As the fracture spacing L_f12D increases, the occurrence of MWPI becomes later. Therefore, if we can find that the pressure curves are distorted at the early flow regime, we can make an original judgment that the stimulated HF may be closed to each other. Similarly, fracture spacing also just splits part of flow regimes (such as first linear-flow regime, bi-linear flow regime, first radial-flow regime and second linear-flow regime), the pressure curves will overlap again subsequently. Similarly, when the gas rate is big enough, the impacts of fracture spacing on MWPI also cannot be observed from pressure curves.

Hydraulic fracture length, L_f1D, L_f2D. We set L_f1D = L_f2D = L_fD to be 1000, 1500, 2000 respectively and ε to be 1:4. Figure 14 illustrates the impacts of hydraulic fracture length L_fD on pressure performance for MWPP scheme. We can systematically analyze the impacts of L_fD on pressure performance from three aspects: (1) before the occurrence of MWPI, for a certain fracture length, the pressure curves of two wells will overlap together. However, when fracture length is varying, the pressure curves will paralleled move upward or downward. As the fracture length increases, the pressure curves will paralleled move downward; (2) when the MWPI occurs, as the fracture length L_fD increases, the occurrence of MWPI becomes earlier. For example, when L_fD = 2000, MWPI basically starts form first linear-flow regime. when L_fD = 1000, MWPI basically starts form first bi-linear-flow regime; (3) when MWPI reaches certain degree, the pressure curves will overlap again subsequently. The bigger the gas rate, the more lately the pressure curves overlap. For example, for well1, the pressure curves will overlap at second linear-flow regime, for well2, the pressure curves will overlap at pseudo-steady diffusion regime. In conclusion, fracture length impacts the whole flow regimes for MWPP scheme.

Stress sensitivity coefficient, ζ. We set ζ_D to be 0, 0.03, 0.05 respectively and ε to be 1:4. Figure 15 illustrates the impacts of stress sensitivity coefficient ζ_D on pressure performance for MWPP scheme. Different from the previous factors, including well spacing, ratio of gas rate, fracture spacing and fracture length, it is almost impossible to find that stress sensitivity coefficient can induce influences on the occurrence of MWPI. However, stress sensitivity coefficient can distort flow regimes at an inverse direction. For example, when ζ = 0.03, radial-flow regimes and pseudo-steady diffusion regimes are distorted. When ζ = 0.05, radial-flow regimes, pseudo-steady diffusion regimes and second linear-flow regime are distorted. As the ζ increases, the distortion of pressure curves becomes severe, and more flow regimes will be distorted. It is also found that the pressure curves do not overlap again due to the existence of stress sensitivity. This meaningful finding can assist us to make a preliminary judgment that the gas rate of two wells is different or the SG reservoir is stress sensitivity, although these two factors distort the flow regimes due to different mechanisms.

Conclusion

To gain better understanding about well performance of MWPP scheme, in this paper, we develop a new semi-analytical pressure transient model in Laplace domain to identify flow regimes without and with IWPI. Model validation is implemented using CMG numerical simulator, and sensitivity analysis is also conducted. Some meaningful conclusions are summarized as following:

• There is good agreement between our model and numerical simulation, moreover, and our approach also gives a much faster calculation speed compared to numerical simulation, both of which demonstrate the accuracy and efficiency of our method;

• Some expected flow regimes are apparently distorted by IWPI. The slope of type curves which characterizes the linear or bi-linear flow regime is no longer equal to 0.5 or 0.25. The horizontal line which characterize radial flow regime is no longer equal to 0.5. For different interference type, IWPI can distort different flow regimes. Interference directly through HF is more rapid than interference through reservoir;

• Well rate and stress sensitivity coefficient mainly determine the distortion of pressure curves. As the well rate decreases or stress sensitivity coefficient increases, the distortion of pressure curves will become severe. Well rate will distort pressure curves when IWPI occurs, on the contrary, stress sensitivity coefficient can distort pressure curves at an inverse direction which is from radial flow to diffusion flow regime, the bigger the stress sensitivity coefficient, the more flow regimes will be distorted;

• Fracture length, well spacing, fracture spacing mainly determine when the IWPI occurs. As the well spacing increases, fracture length decreases, fracture spacing decreases, the occurrence of IWPI becomes later. For well spacing, fracture spacing, when IWPI occurs, pressure curves split, and then overlap again. For fracture length, pressure curves will always split until IWPI reaches certain degree, which is the ending of IWPI flow regime.

Nomenclature

MWPP: Multi-Well-Pad-Production

SG: Shale Gas

IWPI: Inter-Well Pressure Interference

SMFHW: Single Multi-Fractured Horizontal Wells

NFs: Natural Fractures

Hens: Hydraulic Fracture Networks; DPP Dimensionless Pseudo Pressure

DPPD: Dimensionless Pseudo Pressure Derivation

MHFs: Main Hydraulic Fractures

SHFs: Secondary Hydraulic Fractures

T: formation temperature, K

P_i: initial formation pressure, MPa

ψ_i: initial pseudo-pressure, MPa²/(mPa · s)

T_sc: temperature under standard condition, K

P_sc: pressure at standard condition, MPa

C_t: total compressibility, MPa⁻¹

ψ: pseudo-pressure, MPa²/(mPa · s)

ψ_f: fracture pseudo-pressure, MPa²/(mPa · s)

µ: viscosity, mPa · s

h: formation thickness, m

L_ref: reference length, m

Φ: porosity, fraction

L_f1: hydraulic fracture half length of well1, m

L_f: hydraulic fracture half length of well2, m

L_hf12: horizontal distance of the Hens between two wells, m

L_vf: vertical distance of the Hens between two wells, m

t: time, h

x, y: coordination, m

r: radial coordinate in NF system, m

r: radial coordinate in matrix system, m

l: coordination of hydraulic fracture, m

v: integration variable

q₁: well production rate of well1, m³/d

q₂: well production rate of well2, m³/d

q_sc: well production rate under standard condition, m³/d

k_ri: initial permeability of NF system, D

k_f: permeability of hydraulic fractures for well1, D

k_f2: permeability of hydraulic fractures for well2, D

R_m: matrix radius, m

V: gas concentration, sm³/m³

w_f: width of hydraulic fractures for well1, m

w_f2: width of hydraulic fractures for well2, m

W_f1: width of Hens system for well1, m

W_f2: width of Hens system for well2, m

ζ: stress sensitivity coefficient, (mPa · s)/MPa²

ρ: density,g/cm³

C_g: gas compressibility, MPa⁻¹

M₁: total number of hydraulic fracture for well1, integer

M₂: total number of hydraulic fracture for well1, integer

M_w: total number of wells, integer

t_D: dimensionless time

q_fD: dimensionless flux rate

q_cD: dimensionless fracture rate

x_D, y_D: dimensionless space

r_D: dimensionless radial coordinate in NF system, m

r_mD: dimensionless radial coordinate in matrix system, m

L_f1D: dimensionless hydraulic fracture half length of well1, m

L_f2D: dimensionless hydraulic fracture half length of well2, m

L_f12D: dimensionless fracture distance between two wells, m

L_wD: dimensionless distance between two wells, m

C_f1D: dimensionless hydraulic fracture conductivity for well1

C_f2D: dimensionless hydraulic fracture conductivity for well2

C_1D: dimensionless wellbore storage coefficient for well1

C_2D: dimensionless wellbore storage coefficient for well2

u: Laplace variable

Subscript

D: dimensionless

f 1: well1

f 2: well2

Superscript

—: Laplace transform

Acknowledgments

The authors are grateful to the anonymous reviewers for their insightful and constructive comments for improving the manuscript. The authors acknowledge a fund from the National Union Foundation (NUF) of China (No. U1562102) and supports from the MOE Key Laboratory of Petroleum Engineering.

Appendix A Dimensionless definitions

In our research, we introduce a new definition of dimensionless pseudo-pressure based on the total rate of all wells. Therefore, dimensionless pseudo-pressure for SG system and HFNW system can be given as follows, (A1)

For the dimensionless time (A2) where:

For the dimensionless spacing and fracture length (A3) (A4)

For the dimensionless gas rate influx, gas concentration and gas flow rate can be defined, respectively: (A5)

For the fracture-flow model, the hydraulic fracture conductivity can be assumed to be uniform for the same well, respectively, at the same time, the fracture conductivity can be varying for every well. Therefore, dimensionless fracture conductivity can be defined as (A6)

Dimensionless transmissibility coefficient of HF can be presented as: (A7)

Dimensionless storage ratio of SG formation can be presented as: (A8)

Dimensionless adsorption index which denotes the SG desorption ability can be defined as, (A9)

Dimensionless diffusion coefficient which denotes the SG transferring from matrix into natural fracture can be defined as, (A10)

Dimensionless stress sensitivity coefficient can be defined as, (A11)

Appendix B Derivation of line source solution

To develop the mathematical models in SG reservoir, the mathematical formula for NF system and matrix system can be established respectively and then be dynamically coupled. The model development can be given as follows:

Seepage model for SG in NF system. To begin with, flow in natural fractures is assumed to be single phase by obeying Darcy's law. Combining with gas state equation and motion equation, the SG flow in natural fracture system can be described as following diffusion equation with consideration of gas adsorption behavior: (B1)

The pseudo-pressure function was used to account for the pressure dependent gas properties, the governing Equation (B1) with the formula of pseudo-pressure is as follows: (B2)

The stress-dependent permeability can be described by introducing a permeability modular ζ, the relationship between permeability and pseudo pressure can be (Pedrosa, 1986): (B3) where, k_fi is initial permeability of natural fracture under initial pressure condition, Darcy unit; ζ is permeability modular, (mPa · s)/MPa²; ψ_i is the initial pseudo-pressure, MPa²/(mPa · s).

Submitting Equation (B3) into Equation (B2), the final formation of governing formula for natural fracture system can be transformed as follows: (B4)

The initial condition: (B5)

Based on the theory of line sink, the inner condition is presented as follows: (B6)

The reservoir is assumed to be infinite, and the outer boundary is as follows: (B7) where, k_r is the permeability of SG formation, D; ψ is the pseudo-pressure in fracture system, MPa²/(mPa · s); µ is the viscosity of shale gas, mPa · s; Φ is porosity of fracture system, fraction; C_t is the total compressibility coefficient, 1/MPa; t is the time, h; r is radial distance, m. P_sc is the pressure at standard condition, MPa; T_sc is the temperature at standard condition, K; q is production rate of line sink, sm³/d; h is the thickness of the reservoir, m.

For the convenience of solution, some dimensionless variables are defined previously. With the definition of these dimensionless variables, Equations (B1)–(B7) with the dimensionless formation can be presented as follows: (B8) (B9) (B10) (B11)

Equation (B8) shows that the seepage model is strongly nonlinear and some addition method needs to be applied to obtain analytical solution. Here, the perturbation technology and the Presoda transformation are applied to linearize the equations (Pedrosa, 1986): (B12)

According to the theory implemented by Pedrosa (1986), performing a parameter perturbation in ζ_D by defining the following series: (B13a) (B13b) (B13c)

Considering the facts that the ζ_D, dimensionless stress sensitivity coefficient, is always small, thus the zero-order perturbation solution can greatly meets the requirements. The final formation of Equations (B8)–(B11) are as follow: (B14) (B15) (B16)

The Laplace transformation with respect to t_D is then used to deal with Equations (B14)–(B16), the Laplace transform is based on t_D and functions as follows: (B17)

And we can obtain the formation of governing equation of SG reservoir system in the Laplace domain, (B18) (B19) (B20)

Seepage model for SG in matrix system. Due to the ultra-low permeability of matrix, the fluid flowing in the matrix system is treated as unsteady-steady state flow (de Swaan, 1990; Noetinger et al., 2001; Landereau et al., 2001). Thus, we assume that SG flow within shale matrix obeys Second Fick’s law. The shale matrix can be treated as spherical geometry, the SG flowing in the matrix system can be described as follows (Wang, 2014): (B21)

The diffusive flow in the spherical SG matrix blocks is symmetric, thus the center of matrix blocks can be treated as a no-flow boundary, which gives the following inner boundary condition, (B22)

Adsorptive SG concentration on the external surface of the matrix blocks can be evaluated at the gas pressure in the natural fracture system, so the outer boundary condition can be described as follows (B23)

Furthermore, adsorption behavior of SG can be described by Langmuir isotherm equation (Langmuir, 1918). Thus, (B24)

And the initial condition (B25) where, D is SG diffusion coefficient, m²/s; R_m is the radius of spherical matrix, m; V_i is the initial concentration in matrix, sm³/m³; V_E is the SG concentration at the surface of matrix, sm³/m³.

With the dimensionless definitions given in Appendix A, Equations (B21)–(B22) can be written in the following dimensionless formula in Laplace space, (B26) (B27) (B28)

The solution of Equations (B26)–(B28) can be given by (Wang, 2014) (B29)

Coupling seepage model for SG system. The SG flow in NFs system is represented by Equation (B14), the internal source item representing the effects of desorption of SG is given by the following equation based on the spherical matrix blocks assumption, (B30)

Combination of Equation (B29) and Equation (B30) in Laplace domain gives (B31)

Finally, substituting Equation (B31) into Equation (B18), one can get the following governing equation coupling NFs system and matrix system, (B32) where, 

The general solution of Equation (B32) for the infinite-acting SG reservoir can be given by (Ozkan and Raghavan 1991, Xiao et al., 2016), (B33) By the requirement that P_D vanish at infinity, we must have A=0 in equation. From the condition given by Equation (B20), we can obtain (B34)By applying the principle of integration, the pressure distribution of a random position (x_D, y_D) caused by one fracture segment (x_WD, y_WD) is given by the following equation.

(B35)

References

All Tables

All Figures