Numerical Modeling of Hydraulic Fracturing Interference in Multi-Layer Shale Oil Wells

Abstract

Multi-layer horizontal well development and hydraulic fracturing are key techniques for enhancing production from shale oil reservoirs. During well development, the fracturing performance and well-pad production are affected by depletion-induced stress changes. Previous studies generally focused on the stress and fracturing interference within the horizontal layers, and the infilled multi-layer development was not thoroughly investigated. This study introduces a modeling workflow based on finite element and displacement discontinuity methods that accounts for dynamic porous media flow, geomechanics, and hydraulic fracturing modeling. It quantitatively characterizes the in situ stress alteration in various layers caused by the historical production of parent wells and quantifies the hydraulic fracturing interference in infill wells. In situ stress changes and reorientation and the non-planar propagation of hydraulic fractures were simulated. Thus, the workflow characterizes infill-well fracturing interferences in shale oil reservoirs developed by multi-layer horizontal wells. Non-planar fracturing in infill wells is affected by the parent-well history production, infilling layers, and cluster number. They also affect principal stress reorientations and reversal of the fracturing paths. Interwell interference can be decreased by optimizing the infilling layer, infill-well fracturing timing, and cluster numbers. This study extends the numerical investigation of interwell fracturing interference to multi-layer development.

1. Introduction

The stimulation of horizontal wells by large-scale and multi-stage hydraulic fractures is a key technique in the commercial development of shale oil worldwide [1,2,3,4]. Many major shale oil reservoirs have great payzone thicknesses [5,6]. As shown in

Table 1. Comparison of the area, lithology, facies, and cumulative thickness among several shale oil plays (data collected from [5]).

	Eagle Ford Shale	Bakken Shale	Jimsar Shale
Lithology and facies	Marine shale	Marine marlite	Lacustrine shale
Estimated area	2 × 104 km2	7 × 104 km2	3 × 104 km2
Cumulative payzone thickness	30–90 m	2–20 m	80–200 m

, US unconventional plays like Eagle Ford and Bakken are typically in marine facies, while other shale plays such as the Jimsar shale in the Junggar Basin in China are typically in lacustrine facies. These shale plays have great cumulative payzone thicknesses. As a result, effective development of such reservoirs can be obtained by multi-layer development, which involves the placement of horizontal wells with tight vertical spacings, and infill-well drilling and completion can increase the contact area with the matrix and improve production. However, the fracture quality of such vertically infilled wells is governed by the vertical interwell interference and fracture interference [6,7], and accurate quantification of these interferences is critical for the efficient development of shale oil reservoirs associated with relatively thick payzones. When multiple hydraulic fractures propagate in layered shale reservoirs, the interaction between them can negatively affect fracturing qualities and production efficiencies [8]. The geomechanical heterogeneity is a key factor in this process.

It is recognized that the interwell hydraulic fracturing interference is caused by stress and pressure disturbances in shale reservoirs [7,9,10]. Major shale plays such as Eagle Ford, Haynesville, and Marcellus have reported impacts of interwell interference on fracture quality, casing damage, and well-pad production, and the fracture connection between neighboring wells has been identified as a major cause. It is also noticed that the fracture connection can either improve or be damaged, and a case-by-case study is required for each scenario [6,11,12]. Based on more than 3000 datasets from five major US unconventional basins, Miller et al. [13] revealed that hydraulic fracturing interference is more likely to improve well-pad performance in Haynesville and Bakken, while it is more likely to damage it that in Woodford and Niobrara. Lindsay et al. [14] then indicated that discrepancies in parent and infill-well fracture geometries can cause different production decline patterns in these wells. Interwell fracture connections result in the fact that parent and infill wells compete for production, which leads to sudden changes in bottom-hole pressure and production. In addition, dense well placement and large-scale fractures increase reservoir depletion and further the interwell interference. In order to understand the mechanism of interwell fracturing interference and to optimize the overall fracture network quality, Gupta et al. [15] proposed an integrated workflow and indicated that accurate quantification of the temporal and spatial evolution of stress fields is fundamental in the optimization of well spacing and timing for infill-well completion and refracturing. Due to the effect of interference, infill-well fracture geometries can be greatly differ from parent-well fracture geometries, and governing factors include the stress contrast, fracturing fluid properties, parent-well production constraints, rock mechanical properties, infill-well locations, and infill-well fracture designs [16,17,18]. In addition, fluid loss, phase changes, and complex fluid physical changes in such processes should also be taken into account to improve the authenticity of the model strategy [19,20,21].

Depletion-induced in situ stress alterations and hydraulic fracturing modeling in complex stress fields are key components in the evaluation of infill-well fracturing interference. Biot’s consolidation theories and related poromechanical theories are the foundation for the coupled flow and geomechanics modeling of in situ stress evolution caused by hydrocarbon production [22,23]. In the coupled flow and geomechanics modeling, porous media flow is typically coupled with solid mechanics through three techniques: full coupling, explicit coupling, and sequential coupling. Thus, the pore pressure and displacement of rock can be computed, and the depletion-induced stresses can be quantified. Among the coupling techniques, full coupling is deemed as the most stable method that guarantees accuracy and convergence [24,25]. In order to quantify depletion-induced rock deformation, momentum balance is employed to establish the stress equation, and stress–strain constitutive equations are employed to obtain the rock displacement, strain, and stress solutions, which can be used to characterize the temporal evolution of the 3D stress fields in the reservoirs [26,27]. Apart from the accurate modeling of pressure and stress fields, it is also critical to characterize fracture propagation in complex and heterogenous stress conditions [28]. Natural fractures, bedding layers, anisotropy, and fracturing fluid properties all affect fracture propagation patterns and the final fracture geometries [29]. During the propagation of hydraulic fractures, plasticity can be induced near fracture tips, which affects the fracture height growth [30,31,32]. Lee et al. [33] employed the phase-field method and revealed that the fluid phase in the fracture has an impact on the propagation process. The order of the fracturing stages have been identified as another relevant parameter affecting hydraulic fracture qualities [34].

The existing studies of interwell fracturing interference have been insightful. However, the mechanism of hydraulic fracturing interference for infill wells in the multi-layer development of shale oil reservoirs has not been thoroughly discussed. Therefore, a numerical investigation into the relevant parameters’ effects on interwell interference in the multi-layer development of shale oil reservoirs is needed as a reference for the efficient development of such reservoirs.

This study introduces an in-house simulation workflow that includes finite element methods and displacement discontinuity methods and quantifies the porous flow–geomechanics–fracturing physics. It characterizes the production-induced stress evolution and the corresponding non-planar hydraulic fracturing in infill wells in various candidate layers. Based on a realistic shale oil dataset, the interwell fracturing interference between layers of various depths is quantified, and the corresponding interwell fracture networks are simulated. Effects of relevant parameters, including infilling layers, fracturing timing in infill wells, and cluster numbers, are also evaluated. This study provides insight into the understanding of the mechanism of interwell hydraulic fracturing interference in the development of shale oil reservoirs through multi-layer horizontal wells.

2. Coupled Flow and Geomechanics Model

Depletion-induced stress evolution can be calculated by the coupling of fluid flow in porous media with solid mechanics. Mass conservation is employed in the flow equations, and the momentum balance with the Cauchy stress tensor is employed in solid mechanics.

Biot’s consolidation theories express the correlation between pore pressure, total stress, and effective stress [22] as follows:

$${{\sigma }_{\mathrm{ij}}}^{\prime }={\sigma }_{\mathrm{ij}}-b{\delta }_{\mathrm{ij}}p$$

(1)

$$b=1-{\displaystyle \frac{K_d}{K_s}}$$

(2)

where ${{\sigma }_{\mathrm{ij}}}^{\prime }$ is the effective stress; ${\sigma }_{\mathrm{ij}}$ is the total stress; $b$ is the Biot coefficient; ${\delta }_{\mathrm{ij}}$ is the Dirac delta function; $p$ is the pore pressure; $K_d$ is the drained modulus; and $K_s$ is the rock-skeleton modulus.

Based on mass conservation, the water-phase and oil-phase flow equations involve the mass accumulation, flux, and sink/source as follows:

$${\displaystyle \frac{d}{dt}}{\int }_{\Omega }{m}_{i}d\mathsf{\Omega }+{\int }_{\Gamma }{f}_i\cdot \mathrm{nd}\mathsf{\Omega }={\int }_{\Omega }{q}_{i}d\mathsf{\Omega }$$

(3)

and the flux term $f_i$ can be extended as

$$f_i=-{\rho }_i{\displaystyle \frac{k}{{\mu }_i}}(\nabla p-{\rho }_{i}g)$$

(4)

where $m$ is the mass of a phase; $q_i$ is the flow rate; ${\rho }_i$ is the density; $k$ is the permeability; ${\mu }_i$ is the viscosity; $g$ is the gravity.

For phase i, the compressibility is used to correlate pore pressure and density as

$$c_i={\displaystyle \frac{1}{{\rho }_i}}{\displaystyle \frac{\partial {\rho }_i}{\partial p}}$$

(5)

where $c_i$ represents the compressibility.

Based on the quasi-static assumption and Cauchy stress, the balance between traction boundaries and normal stresses and shear stresses in a three-dimensional space is presented as

$$\nabla \cdot \sigma +{\rho }_{b}g=0$$

(6)

Considering the poromechanical effects in saturated porous media, the relationship between stress, pore pressure, and strain is written as

$$\sigma ={\sigma }_0-b\left(p-p_0\right)\mathbf{1}+C_{dr}:\epsilon$$

(7)

where $C_{\mathrm{dr}}$ is the stiffness matrix, and ε is the displacement.

Based on infinitesimal transformation, the displacement of the rock skeleton and the strain are expressed as

$$\epsilon ={\displaystyle \frac{1}{2}}({\nabla }^{T}u+\nabla u)$$

(8)

where u is the displacement.

Full coupling of flow and geomechanics is employed in this study. In the flow equation, the fluid velocity of a certain phase is regarded as an interstitial velocity, which consists of the relative velocity between fluid and solid rock, and the rock phase velocity representing deformation as

$$v_i=v_{is}+v_s$$

(9)

where the rock phase velocity is very small as it represents the solid rock deformation rate. Porosity is time-dependent, and Biot’s coefficient is used to couple pore pressure effects and the volumetric strain as

$${\displaystyle \frac{\partial \varphi }{\partial t}}={\displaystyle \frac{b-\varphi }{K_s}}{\displaystyle \frac{\partial p}{\partial t}}+\left(b-\varphi \right){\displaystyle \frac{\partial {\epsilon }_v}{\partial t}}$$

(10)

where $\varphi$ is the porosity; $t$ is the time; and ${\epsilon }_v$ is the volumetric strain.

In this full coupling strategy, the dependent variables (displacement, velocity, and pressure) are solved simultaneously in the same numerical system at each individual time step. This guarantees the stability and convergence rate. In contrast, in other coupling strategies such as explicit coupling and sequential coupling, dependent variables are solved in different steps and different numerical systems, and the simulation may not reach convergence or stable solutions [35,36].

Mixed finite element methods, including the standard Galerkin method and the discrete Galerkin method, are used for the space discretization of the governing equations. The finite difference method is used for time-stepping with the backward Euler method to guarantee the numerical stability and convergence. The Newton–Raphson method is used to address nonlinearity in the solution. Specifically, the standard Galerkin method is used for the solution of pressure and saturation in the water and oil phases in the porous media; the discrete Galerkin method is used for the solution of deformation. The Newton–Raphson method iteratively linearizes the solution for dependent variables, where partial derivatives of residuals in Jacobian matrices are solved.

3. Non-Planar Hydraulic Fracture Model

In the multi-layer development of shale oil reservoirs, production in tightly spaced horizontal wells alters the magnitude and orientation of in situ stress, which increases the complexity of the stress field and leads to complex interwell fracture networks. In this work, a pseudo-3D non-planar fracture model based on the displacement discontinuity method is employed to simulate the simultaneous fracture propagation, fluid flow in fractures and in the wellbore, shear failure and tensile failure of fractures, and stress shadow effects.

The shear failure and tensile failure are solved by the tangential and normal displacement discontinuities as

$${\sigma }_s^i=\sum _{j=1}^n{A}_{s,s}^{ij}{D}_s^i+\sum _{j=1}^n{A}_{s,n}^{ij}{D}_n^i$$

(11)

$${\sigma }_n^i=\sum _{j=1}^n{A}_{n,s}^{ij}{D}_s^i+\sum _{j=1}^n{A}_{n,n}^{ij}{D}_n^i$$

(12)

where ${\sigma }_s$ is the traction boundary for shear stress; ${\sigma }_n$ is the traction boundary for normal stress; A is the elastic-influence matrix; and D is the displacement discontinuity.

Traction boundaries can be further expressed by considering the fluid pressure in fractures and maximum and minimum principal stresses as

$${\sigma }_n^i=p_{elem}^i-{\sigma }_n^{r,i}$$

(13)

$${\sigma }_s^i=-{\sigma }_s^{r,i}$$

(14)

where $p_{elem}$ is the pressure in a fracture element.

The fluid flow in fractures and the wellbore is considered in the fracture model, and the first and the second Kirchhoff’s laws are used to describe the fluid flow distributions in simultaneously propagated fractures. The pressure drop is calculated by lubrication theories. In the model, the total pumping rate is a constant, and the dynamic distribution of fracturing fluids in different fracture clusters is calculated by considering the fractions in the wellbore, perforations, and fractures [37,38,39,40] as

$$p_h=\mathsf{\Delta }{p}_{wf,m}+\mathsf{\Delta }{p}_{pf,m}+p_{ff,m}$$

(15)

where $p_h$ is the pressure at the heel of the horizontal wellbore; $\mathsf{\Delta }{p}_{wf,m}$ is the friction in the wellbore; $\mathsf{\Delta }{p}_{pf,m}$ is the friction in the perforations; and $\mathsf{\Delta }{p}_{ff,m}$ is the friction in the fractures.

The limitation of the derivation above is that the fracture closure effect and the associated proppant transport effect cannot be quantified, and a constant conductivity is assumed in the study.

In the implementation of DDM, mechanical displacements along the fracture are described. This includes the opening and sliding displacements and the fluid flow dynamics within the fractures and wellbore, adhering to Kirchhoff’s laws for fluid conservation and using the lubrication theory for pressure calculations. The implementation of the solution algorithm involves using the backward Euler method for time-stepping to enhance stability and employing the Newton–Raphson method to handle nonlinearities in the system.

4. Numerical Study and Discussion

First, the coupled flow and geomechanics model is validated against the widely recognized Terzaghi’s problem for consolidation [41,42]. This focuses on the 1D consolidation in homogeneous, isotropic, and saturated porous media. In this validation, the porous media are loaded vertically at the top. The top boundary also has a drainage condition for pressure dissipation. The bottom boundary is impermeable and fixed. The initial excess pore pressure is 100 kPa, the coefficient of consolidation is 1 × 10⁻⁶ m²/s, and the total length of the 1D domain is 1 m. The initial excess pore pressure here represents the additional pressure due to the application of the load at the top. The validation results are shown in Figure 1. The deformation results and pressure distribution results at various locations and time steps are compared between the lab results and the numerical modeling results. The total settlement is the greatest deformation accumulated at the top boundary, which keeps increasing as time progresses. The pressure values are greater at and near the top boundary as the excess pressure is induced by the compressive top boundary. However, pressure tends to drop with time, as the top boundary is also a drainage boundary, where pressure can dissipate freely. Based on the results, the model is validated against the Terzaghi problem.

The numerical study is carried out based on a simplified single-stage model with reservoir and rock mechanical parameters from the Jimsar shale oil reservoirs. It is located in the Junggar Basin in northwestern China, and it is one of China’s first commercially developed shale oil plays. Payzones are identified in the Upper and Lower Permian Lucaogou Formations with great thicknesses and different depths. In order to effectively exploit hydrocarbons in these formations, multi-layer horizontal wells with hydraulic fracturing are necessary. The placement of infill wells in layers next to layers already under horizontal well development is regarded as a strategy for further improving well-pad production, and an interwell interference study focused on the vertical scale for the upper and lower formations is beneficial for the strategy design.

The numerical study consists of modeling the stress evolutions induced by parent-well depletion and the hydraulic fracturing simulation for infill wells placed at several candidate locations in layers at different depths. In the study, planar fractures are assumed in the parent well. Local grid refinements are used to denote the complex fracture networks.

A dataset from this field is used to construct the numerical model. Microseismic monitoring for a hydraulically fractured horizontal well is carried out, which shows the fracture network geometry for the parent well as in Figure 2. In this study, the production of the last stage is considered in the simulation of the parent-well production. Its fracture dimension is shown in

Table 2. Fracture dimensions of the last stage in the parent well.

Fracture Half-Length	Fracture Height	Fracture Network Width
100 m	15 m	20 m
Fracture direction	Number of events	Main fracture number in a stage
NE 26°	62	3

. Based on the dimensions of the last fracture stage, a single-stage hydraulic fracture model for the parent well is constructed as in Figure 3, with the respective reservoir and fracture parameters recorded in

Table 3. Reservoir parameters and parameters for the fracture stage in the parent well.

Parameter	Value
Matrix permeability	0.01 mD
Matrix porosity	12%
Young’s modulus	25 GPa
Poisson’s ratio	0.22
Oil density	0.9 g/cm3
Initial reservoir pressure	27 MPa
Minimum principal stress	30 MPa
Maximum principal stress	35 MPa
Oil saturation	75%
Fracture half-length	100 m
Fracture height	15 m
Fracture number	3
Fracture spacing	10 m
Cumulative thickness	65 m
Location of the horizontal wellbore	X = 0 m, −100 m < y < 100 m, z = 0 m

.

4.1. Depletion-Induced Stress Evolutions in Different Layers

The production of the parent well alters the in situ stress magnitude and orientation, which induces interlayer stress interference and affects the fracture quality of the infill well in an adjacent layer [43,44] Therefore, the temporal and spatial evolutions of in situ stresses in different layers are first characterized.

Figure 4 shows the simulated pore pressure and maximum principal stress orientation distributions after 5 years of parent-well production in the layer where the parent well is located. Figure 4a–c show the results after 1 year, 2 years, and 5 years, respectively. The results indicate that pressure depletion is mainly distributed in the stimulated reservoir volume (SRV), and it is not significant outside the SRV. In Figure 4c, 5-year production drops the pressure to 12 MPa in the fractures, and stress reorientation mainly occurs at and around the fracture tips, indicating stress interference within the same layer. This interference is not strong after the 1-year production, as shown in Figure 4a, while it becomes more noticeable as the production time increases to 5 years. The stress reorientation can be expressed by

$$\mathit{tan}2\theta ={\displaystyle \frac{2\tau }{{\sigma }_x-{\sigma }_y}}$$

(16)

Figure 5 shows the minimum principal stress during the 5-year parent-well production in the vertical plane in the x–z dimension, which is located 50 m away from the parent-well fracture tips in the y direction. In the specific result, the vertical scale is focused between −140 m and +140 m. In Figure 5a, results after 1 year of production are shown; in Figure 5b,c, results after 2 years and 5 years of production are shown.

The results in Figure 5c show that the depletion-induced stress disturbance is within 0.2 MPa, and the disturbance is restricted within the same width as the parent-well SRV. The disturbance does not change with depth. This indicates that the width of the parent-well SRV is a key factor affecting the depletion-induced stress alterations outside the SRV. In addition, the stress disturbance increases with time. This indicates that although pore pressure is decreased during the depletion caused by the production, effective stress is increased as the compressive state is strengthened in the solid rock skeleton. The spectrum-like area indicates that this type of elevation is highly correlated with the fracture geometry.

Figure 6 shows the depletion-induced minimum principal stress changes in the x–z plane crossing the SRV after 5 years of parent-well production. In the presented time steps, the decrease in the maximum principal stress is significant in the SRV. Increases in the maximum principal stress are observed outside the top and the bottom of the SRV. The stress changes in the x direction outside the SRV are not as significant as those in the z direction. In Figure 6a, the lowest stress is around 28 MPa after 1 year of production, while in Figure 6b, the lowest stress drops to 25.5 MPa; in Figure 6c, the lowest stress drops to 25 MPa, and the extent of the pressure drop is the greatest among the three different time steps. Also, it is noted that from the 1-year to the 5-year production, the areas of stress increase near the top and bottom boundaries also expand.

The results in Figure 5 and Figure 6 indicate that the minimum principal stress in vertical planes crossing the parent-well SRV is affected by both the width and the height of the SRV, while it is mainly affected by the SRV width for x–z planes outside the SRV. Furthermore, increases in the minimum principal stress are observed outside the SRV, which are explained by the effective stress: in these areas, the contribution of the decrease in total stress caused by pressure depletion is smaller than the contribution of the increase in effective stress induced by hydrocarbon production, which results in the increase in total stress.

Figure 7 further presents the 1D distributions of the minimum principal stress at different layers in vertical planes at different time steps.

In the vertical plane outside the parent-well SRV, the minimum principal stress increases from 30 MPa to 30.2 MPa during the production of the parent well, which can also be explained by the effective stress theory. Also, the minimum principal stress profiles at various depths are generally the same at the same time step, and the insignificant discrepancy increases with production time, indicating that the depth effect is not strong outside the SRV.

The observations are different in the vertical plane crossing the parent-well SRV, where the minimum principal stress is affected by both production time and the depths of the layers. The most significant decrease in ${\sigma }_{hmin}$ is in the layer at z = 0 m (parent-well layer), which decreases from 26.6 MPa in the first year to 25.5 MPa in the fifth year. The ${\sigma }_{hmin}$ decrease profile is also greatly affected by the fracture geometry. For the layer at z = 15 m, ${\sigma }_{hmin}$ increases with the production time, and the profile is governed by discrete fractures. The ${\sigma }_{hmin}$ finally increases to 31.3 MPa after a 5-year depletion. For the layer at z = 30 m, ${\sigma }_{hmin}$ still slightly increases with the production time, while the profile is not greatly affected by the parent-well fracture geometry, as the effect of SRV depletion on the layers away from the parent well becomes limited.

4.2. Non-Planar Hydraulic Fracturing in Vertically Infilled Wells

After the quantitative characterization of depletion-induced in situ stress alterations, hydraulic fracturing modeling is carried out for the vertically infilled horizontal well at several candidate depths to numerically investigate the interwell fracturing interference in the multi-layer development of shale oil reservoirs with great cumulative thicknesses. As shown in Figure 8, three depths in the vertical plane 50m away from the parent-well fracture tips in the y direction are used as candidate infill-well locations, and the hydraulic fracturing is simulated correspondingly. Then, the effects of fracturing timing, infill-well layers, and cluster numbers on the interwell fracturing interference are quantitatively investigated.

Figure 9 demonstrates the infill-well fracturing results at three candidate well depths after the 1-year parent-well production. Note that planar fractures in the parent well are at depth z = 0 m, and the infill-well non-planar fractures are simulated at various depths of 0 m, 15 m, and 30 m. The fractures presented in the results are projected in the z direction to exhibit their interwell extent in the x–y planes. Two cluster densities are simulated: 3 clusters and 4 clusters per stage. The results show that the non-planar feature of infill-well fractures is strongest when the infill well is in the same layer as the parent well. When the infill well moves 30m away from the parent-well layer, the non-planar feature of the infill-well fractures becomes the weakest. Stress reorientation (as presented in Figure 4) near the parent-well fracture tips make the outer infill-well fractures divert when they approach the parent-well fractures, which avoids further interwell fracture interference between parent and infill-well fractures. The inner infill-well fracture(s) show non-planar propagations, while the stress reorientation does not divert them away from parent-well fractures. They lead to the crossing of parent and infill-well projections on the x–y plane, which are observed for inner fracture(s) in both the 3-cluster and 4-cluster cases in Figure 9. In the 4-cluster cases, the layer at z = 0 m does not exhibit fracture hits in the z-direction projection, while layers at z = 15 m and 30 m exhibit hits of fracture projections in the z direction. Considering the fact that the layer at z = 0 m endures the strongest depletion-induced stress alteration and layers at z = 15 m and 30 m show weaker depletion-induced stress alterations, it is concluded that when the parent-well depletion is strong enough, the sufficiently induced stress reorientation can divert infill-well fractures and protect parent-well fractures while fracture hits or hits of fracture projections in the z direction occur in layers at other depths enduring intermediate or low depletion-induced stress alterations. In addition, it is noted that as the cluster number increases from three to four, the non-planar feature of the infill-well fractures becomes stronger, indicating stronger interwell interference.

The results in Figure 9 show that the infill location and cluster number jointly affect the interwell and interlayer interference. It is also concluded that sufficient parent-well-depletion-induced stress reorientation can, in fact, protect parent-well fractures, while insufficient depletion-induced stress reorientation cannot divert infill-well fractures or avoid fracture hits.

Table 4. Infill−well fracture length comparison at candidate locations at various depths (z = 0 m, 15 m, and 30 m) after a 1-year parent−well production (3 clusters and 4 clusters per stage).

Infill-well Fracture	Fracture Length (from Left to Right)
z = 0 m, 3-cluster	202.24 m, 179.44 m, 201.93 m
z = 15m, 3-cluster	199.07 m, 179.47 m, 199.93 m
z = 30 m, 3-cluster	199.95 m, 178.85 m, 201.18 m
z = 0 m, 4-cluster	210.67 m, 182.50 m, 181.04 m, 211.69 m
z = 15 m, 4-cluster	209.40 m, 181.39 m, 183.05 m, 211.92 m
z = 30 m, 4-cluster	209.44 m, 180.21 m, 181.43 m, 209.78 m

documents the lengths of individual fractures in all the scenarios.

Then, Figure 10 presents the infill-well fracturing results at candidate well depths after the 3-year parent-well depletion. As the parent-well production time increases, depletion-induced stress alterations become stronger, and the non-planarity of infill-well fractures becomes stronger as well. This indicates that the time of legacy production in parent wells is another key factor affecting the interwell and interlayer interference. Compared with the results in Figure 9, the increase in the magnitude of interference is weaker for 3-cluster cases and stronger for 4-cluster cases; hits of parent-well and infill-well fractures or fracture projections become stronger or more noticeable in 4-cluster cases.

Subsequently, Figure 11 presents the results after the 5-year parent-well production. Comparatively, the infill-well fractures’ non-planarity becomes even stronger for both 3-cluster cases and 4-cluster cases. For 4-cluster cases, the stronger stress reorientation around parent-well fracture tips pushes infill-well fractures away and effectively avoids interwell fracture connections. However, interwell and interlayer interferences in these cases are so strong that inner fractures hit outer fractures at the same stage in infill-well fracturing, indicating intense inter-cluster fracturing interference.

Table 5. Infill−well fracture non-planarity data after 1 year, 3 years, and 5 years of production. The 3-cluster scenario is used for comparison.

	1-yr, z = 0 m	1-yr, z = 15 m	1-yr, z = 30 m	3-yr, z = 0 m	3-yr, z = 15 m	3-yr, z = 30 m	5-yr, z = 0 m	5-yr, z = 15 m	5-yr, z = 30 m
Cluster 1	2.16	2.19	2.15	2.17	2.17	2.17	2.21	2.2	2.18
Cluster 2	1.78	1.79	1.78	1.79	1.78	1.79	1.81	1.79	1.79
Cluster 3	2.09	2.13	2.13	2.12	2.12	2.14	2.14	2.14	2.13

summarizes the non-planarity of each individual fracture in the 3-cluster cases after 1 year, 3 years, and 5 years of production. The non-planarity here is quantified by the ratio of a curve (fracture)‘s total length to the straight-line distance between its endpoints. Greater values indicate stronger non-planarity magnitudes.

The results presented above show that infill-well fracturing timing, cluster design, and infill-well depths all play important roles in affecting the interwell and interlayer fracturing interference when developing a shale oil reservoir with multi-layer horizontal wells. The interference is the strongest when the parent and infill wells are in the same layer. As the infill-well candidate layer moves away from the parent-well layer, the interference decreases. Also, the interference increases with the increase in cluster number. However, it is noted that the stronger interwell and interlayer interference does not always indicate a higher likelihood for the occurrence of interwell and interlayer fracture hits or connections. When the interwell and interlayer interference is sufficient, the induced stress reorientation can divert infill-well fractures away from parent-well fractures and, in fact, protect parent-well fractures.

In addition, the diverted fracture paths also have an impact on infill-well productivity. When infill-well fracture lengths are similar, the pressure depletion magnitude at and around the fractures becomes the biggest key factor affecting infill-well productivity. When infill-well fractures are diverted away from parent-well fractures, the infill-well fractures tend to control more undepleted reservoir areas where the oil saturation and pressure gradient (once drained) are higher. Also, the effect of this type of interference on the productivity of the parent well cannot be neglected. In reservoir areas simultaneously drained by infill-well fractures (new fractures) and parent-well fractures (legacy fractures), productivity competition between the new and legacy fractures is generated, and the production efficiency in both infill and parent wells can be impacted.

5. Conclusions

This work presents a modeling workflow that includes a coupled flow and geomechanics simulation for depletion-induced in situ stress alterations and a non-planar hydraulic fracture simulation for infill-well fractures and the corresponding interwell and interlayer fracturing interference. The contribution of the work is the extension from interwell interference studies in the same horizontal layers to the interwell and interlayer interference caused by infill-well fracturing in various layers at different depths. Therefore, the numerical analysis can be a reference for the development of shale oil reservoirs with relatively thick payzones where multi-layer horizontal wells and hydraulic fracturing techniques should be used.

In conclusion:

(1)

Pressure depletion induced by parent-well production is significant within the SRV, while the correspondingly induced in situ stress alterations are observed both inside and outside the SRV.

(2)

The interwell and interlayer interference is jointly affected by infill-well fracturing timing, cluster design, and the candidate infill well’s depth. The interference is the strongest when the parent and infill wells are in the same layer. In addition, the interference is augmented by the prolonged legacy production in the parent well and by an increased cluster number per fracture stage in the infill well.

(3)

Stronger interwell and interlayer interference does not always indicate a greater tendency for fracture hits or connections between parent- and infill-well fractures. When the parent-well-depletion-induced in situ stress alteration is sufficient, infill-well fractures propagated at various depths can be pushed away by the stress reorientation, and parent-well fracture protection can be achieved.

(4)

Based on the observations in this specific numerical study, it is recommended that a combination of coupled flow and geomechanical simulation and infill-well fracturing prediction is carried out before the placement and completion of infill wells. It helps to address the non-monotonic relationship between interwell interference and the fracture hit potential.