Three-Dimensional Geomechanical Modeling and Analysis of Refracturing and “Frac-Hits” in Unconventional Reservoirs

Abstract

Field experience has demonstrated that infill well fractures tend to propagate towards the primary well, resulting in well-to-well interference, or so-called “frac-hits”. Frac-hits are a major concern in horizontal well refracturing because they adversely affect the productivity of both wells. This paper provides a 3D geomechanical study of the problem for the first time in order to better understand frac-hits in horizontal well refracturing and its mitigating solutions. To our knowledge, this is the only refracturing study focused on fracture mechanics and within the context of coupled proroelasticity using a single model. The modeling is based on the fully coupled 3D model, GeoFrac-3D, which is capable of simulating multistage fracturing of multiple horizontal wells. The model couples pore pressure to stresses, and makes it possible to create dynamic models of fracture propagation. The modeling results show that production from production well fractures leads to a nonuniform reduction of the reservoir pore pressure around the production well and in between fractures, leading to an anisotropic decrease of the total stress, potentially resulting in stress reorientation and/or reversal. The decrease in the total stress components in the vicinity of the production fractures creates an attraction zone for infill well hydraulic fractures. The infill well fractures tend to grow asymmetrically towards the lower stress zone. The risk of frac-hits and the impact on the “parent” and “infill” well production vary according to the reservoir stress regime, in situ stress anisotropy, and production time. By optimizing well and fracture spacing, fracturing fluid viscosity, and the timing of refracturing job, frac-hit problems can be minimized. The simulation results demonstrate that the risks of frac-hits can be potentially mitigated by repressurization of the production well fractures before fracturing the infill well.

1. Introduction

In unconventional reservoir development, the refracturing of horizontal wells is a cost-effective technique for production enhancement. The refracturing technique can restore or increase well productivity and may provide additional reserves by improving hydrocarbon recovery. Refracturing the original production well or drilling an infill (or “child”) well near the primary production (or “parent”) well are commonly used refracturing techniques in the oil and gas industry (e.g., [1,2,3]). The reservoir depletion near the production well creates a localized “pressure sink” zone with reduced total stresses. If an infill well is drilled in the vicinity of this reduced zone and is stimulated, its hydraulic fractures tend to grow towards the depleted zones (e.g., [4,5,6]), which may result in communication between the production and infill wells, which is usually called “frac-hits”. In addition, fracture fluid and proppant communication may occur between the neighboring wells, which also contributes to an ineffective stimulation of the infill well. Several causes of frac-hits have been reported, including reservoir stress reversal, the tendency of infill well fractures to grow towards existing well fracture networks, and improper well spacing. Longer duration of production from the parent well also leads to frac-hits. It has been reported that the longer the initial well has been in production, the greater the likelihood of an induced stress change, resulting in a higher probability of problems with frac-hits. Frac-hits usually reduce the productivity of the production wells, and in certain conditions, they may affect the efficacy of the parent well. Frac-hits also impact the production potential of subsequent infill wells, since both wells drain from the same reservoir area. Repressurization of the production well before infill well stimulation has been suggested as a mitigation method for frac-hits (e.g., [1,2]).

In the last three decades, several numerical and experimental studies have been presented to demonstrate and quantify the effects of fluid injection or reservoir depletion on reservoir pore pressure and stress state change [7,8,9,10,11,12,13,14]. Perkins and Gonzalez [7] studied thermoelastic effects on stress reorientation/reversal and fracture propagation. Since then, several numerical and experimental investigations have demonstrated the effects of injection/production on reservoir pore pressure and stress [9,10,11,12,13,15]. Several field tests conducted by Siebrits et al. [16] showed a significant increase in production due to refracture reorientation in Barnett shale. Recently, several modeling and field studies have been carried out to study the impact of reservoir depletion on pore pressure and stress variation in the infill well zones [1,2,3]. Rezaei et al. [17,18] presented a refracturing analysis using a two-dimensional (2D) poroelastic displacement discontinuity model built using the work of Ghassemi and Zhang [14], and Chun [19]. However, a 2D model cannot present a complete picture of horizontal well refracturing because of the 3D nature of the stress perturbations and fracture propagation. Roussel et al. [20] and Safari et al. [20] addressed aspects of infill well fracturing in the case of closely spaced multiple horizontal wells using 3D models. Through extensive parametric analysis, Safari et al. [21] demonstrated that protecting parent wells from infill well fractures has a greater chance of success in a reservoir with low stress anisotropy, low hydrocarbon viscosity, low bottomhole production pressure, high reservoir pore pressure gradient, high reservoir permeability, and high Biot’s constant. Sangnimnuan et al. [22] developed a sequentially coupled reservoir-geomechanics model for reservoir stress and pore pressure analyses. An embedded discrete fracture model (EDFM) was used to model the fractures, and a 2D DD-based model was used to model the infill well fracture propagation. Most of these earlier refracturing simulations used two different models for the reservoir depletion analysis and subsequent fracturing of the infill wells, potentially leading to relatively high computational cost and data interpolation errors. In contrast, the GeoFrac-3D model used in this study can simulate both aspects of the refracturing simulation using a single model.

Understanding the 3D stress state resulting from depletion-induced pore pressure change is vital for a successful refracturing treatment. This paper presents a 3D geomechanical analysis of production-induced pressure sink and local stress redistribution to address frac-hits. A key parameter identified is the in situ stress difference. Possible fracture-hit mitigation is proposed based on numerical simulations. In this paper, we use a fully coupled 3D poroelastic model GeoFrac-3D for the numerical simulations. The model is based on the 3D poroelastic displacement discontinuity (DD) method for the solid rock deformation and pore fluid diffusion in the rock matrix. The fluid flow inside the fractures is simulated using the Galerkin’s finite element method (FEM), and fracture propagation is treated in the framework of the linear elastic fracture mechanics approach (LEFM). The GeoFrac-3D model is capable of simulating multistage fracturing of horizontal wells. Both aspects of refracturing simulation (i.e., reservoir depletion analysis and subsequent infill well fracture creation) are included in GeoFrac-3D. Also, to investigate the impact of reservoir layering heterogeneity on depletion-induced reservoir pore pressure and stress change, we also used a state-of-the-art coupled Galerkin’s FEM model and compared its results to those from GeoFrac-3D. A brief description of the mathematical model is presented first, followed by a detailed analysis of reservoir depletion and infill well stimulations

2. Theory and Governing Equations

2.1. Poroelastic Deformation of Rock Matrix and Pore Fluid Diffusion

The porous rock deformation and pore fluid diffusion in the rock matrix are governed by linear poroelasticity theory [23,24,25,26]. Linear poroelasticity theory describes the relationship between (i) the rate of change of bulk volume due to pore pressure change, and thus, the rate of loading, and (ii) pore pressure changes due to mean stress change [26]. The constitutive relationships, transport laws, and conservation laws are usually used to explain the behavior of porous rock. A brief description of each component is provided in the following sections.

2.1.1. Constitutive Equations

The constitutive equations for a poroelastic rock mass under isothermal condition are givens as [25,26,27]:

$${\sigma }_{ij}=2G{\epsilon }_{ij}+\frac{2G\nu }{\left(1-2\nu \right)}{\delta }_{ij}{\epsilon }_{kk}-\alpha {\delta }_{ij}p;\hspace{1em}\hspace{1em}(i,j=1,2,3)$$

(1)

$$\zeta =\alpha {\epsilon }_{kk}+\frac{1}{M}p$$

(2)

where${\epsilon }_{ij}$ is the strain tensor, ${\sigma }_{ij}$ is the total stress component, ${\sigma }_{kk}$ and ${\epsilon }_{kk}$ represents the summation of the stresses and strains, respectively, p is the pore pressure,$\zeta$ is the fluid content variation per unit volume of the porous rock, G is the shear modulus, $\nu$ $is$ Poisson’s ratio, ${\delta }_{ij}$ is the Kronecker delta function and $\alpha$ is the Biot’s effective stress coefficient. In above equation, sign convention is tension positive. M represents the Biot’s modulus, defined as [24]:

$$\frac{1}{M}=\frac{\varphi }{K_f}+\frac{\alpha -\varphi }{K_s}; \frac{1}{K_s}=\frac{1-\alpha }{K}$$

(3)

where $K_f$ is the fluid bulk modulus, $K_s$ is the bulk modulus of the solid grains, $\varphi$ is the formation porosity and K is the formation bulk modulus.

2.1.2. Transport Equation

The pore fluid diffusion in a porous rock matrix is governed by Darcy’s law as [25]:

$$q_i=-\left(\frac{k}{\mu }\right)\left(p_{,i}-f_i\right)$$

(4)

where $q_i$ is the fluid flux defined as the rate of fluid volume flow across a unit area, k is the intrinsic reservoir permeability, μ is the fluid viscosity, $p_{,i}$ represents the pressure gradient in the direction “i”, and $f_i={\rho }_f{g}_i$ is the body force per unit volume of fluid, where ${\rho }_f$ is the fluid density, and $g_i$ is the gravity component in the direction “i”.

2.1.3. Conservation Laws

The conservation laws for poroelastic rocks are given using two equilibrium conditions. The static equilibrium condition leads to the local stress balance equation or equilibrium, as:

$${\sigma }_{ij,j}=-F_i$$

(5)

where $F_i$ is the body force per unit volume in the direction “i”. The mass conservation condition for a compressible fluid gives the local continuity equation as:

$$\frac{\partial \zeta }{\partial t}+q_{i,i}=\gamma$$

(6)

where γ is the intensity of the fluid source or sink.

2.1.4. Field Equations

The field equations for poroelastic rock deformation and pore fluid diffusion in the rock matrix can be given based on the constitutive relationships, fluid diffusion equation, and the conservation laws. If, there is no external fluid source or body force acting on the rock mass, the substitution of the stress constitutive equation (Equation (1)) into static equilibrium conditions (Equation (5)) results in the Navier’s equations with coupling terms for the poroelastic rock matrix deformation as [25,27]:

$$G{\nabla }^2{u}_i+\frac{G}{\left(1-2\nu \right)}{u}_{k,ki}-\alpha p_{,i}=0$$

(7)

where $u_i$ is the solid rock displacement in direction “i”. The diffusion equation of the pore fluid is obtained by substituting Darcy’s law (Equation (4)) and the constitutive pore pressure relationship (Equation (2)) into the continuity equation (Equation (6)) as [25,27]:

$$\frac{1}{M}\frac{\partial p}{\partial t}=\left(\frac{k}{\mu }\right)\nabla {}^{2}p-\alpha \frac{\partial \epsilon }{\partial t}$$

(8)

The above equations present the governing equations for the linear poroelastic analysis of a porous rock mass. Depending on the nature of the given problem, these equations can be solved for the mechanical deformation of the rock mass and pore fluid diffusion in the rock matrix using the appropriate initial and boundary conditions.

2.2. Fluid Flow Inside a Fracture

For the fluid flow simulation, a hydraulic fracture can be assumed as a large channel separated by two smooth parallel plates with a very narrow and smoothly varying aperture. The “parallel plate” or lubrication theory is a fundamental assumption that is made by nearly all hydraulic fracture models [28]. The “parallel plate” assumption is valid for hydraulic fractures, since a hydraulic fracture usually has a height and length that are significantly greater than its aperture. Based on these assumptions, the lubrication theory or parallel plate model, or “cubic law”, is applicable for the fluid flow simulation inside a hydraulic fracture [29]. Recent field evidence from a hydraulic fracturing test site (HFTS) [30] showed that the hydraulic fractures were, in fact, nearly smooth or not varying extensively. The fluid flow inside a hydraulic fracture is governed by the Navier-Stokes’ equation, which gives the relationship among the fracture aperture, fluid flux, and the pressure gradient using the momentum balance condition, as [31]:

$$\mathit{q}\left(s,t\right)=-\frac{w^3\left(s,t\right)}{12\mu }\nabla p_f\left(s,t\right);\hspace{1em}s\in A$$

(9)

where $\mathit{q}\left(s,t\right)$ is the fluid flux, $w\left(s,t\right)$ represents the fracture hydraulic aperture,$\nabla$ is a 2D gradient operator, $p_f\left(s,t\right)$ is the fluid pressure inside the fracture, µ is the fluid viscosity, $s=\left(\mathbf{x}{}^{\prime },y{}^{\prime },0\right)$ represents the local coordinates of fracture plane (A), and t is the fluid injection time. In addition to the momentum balance condition, the fluid continuity condition is required to describe the behavior of fluid flow inside the fractures. The continuity condition for an incompressible fluid flow inside a fracture is given as [32,33,34]:

$$\nabla ·q\left(s,t\right)=-\frac{\partial D_n\left(s,t\right)}{\partial t}-2q_L\left(s,t\right)+\delta \left(s-s_i\right)·Q{}_i-\delta \left(s-s_e\right)·Q{}_e$$

(10)

where $\partial D_n/\partial t$ is the fracture volume change rate per unit area, $D_n$ represents fracture mechanical aperture, $Q_i$ and $Q_e$ are the fluid injection and extraction rates, respectively, $\delta$ is the Dirac delta function, $s_i$ and $s_{e }$ are the local coordinates of the injection and extraction wells, respectively, and $q_L$ is the fluid leak-off rate from one side of the fracture surface into the reservoir matrix, which is governed by Darcy’s law, as in Equation (4). The multiplication factor 2 in the leak-off term $q_L$ corresponds to two surfaces of a fracture, which is only valid for symmetric fractures in homogeneous rock, as in this study. However, in cases where the reservoir rock is not homogeneous, two components for the fluid leak-off rate should be used for each side of the fracture. The DD model presents a fully coupled poroelastic model, so that fluid diffusion or leak-off into the rock matrix through fracture surfaces is governed by Darcy’s law, as in Equation (4). The fundamental DD solution is poroelastic, so that pressure variations in the reservoir and their impact on the leak-off rate are implicitly included. In this way, we also do away with the need for an empirical leak-off coefficient, which is another relative novelty of our model. By combining Darcy’s law for the fluid diffusion in Equation (4), the momentum balance condition (Equation (9)), and the fluid continuity condition in Equation (10), the governing equation for an incompressible fluid inside a hydraulic fracture can be written as:

$$\begin{array}{ll}\nabla ·\{\frac{w^3\left(s,t\right)}{12\mu }\nabla p_f& \left(s,t\right)\} + 2\left(\frac{k}{\mu }\right)·\frac{\partial p_f\left(s,t\right)}{\partial n}\\ & =\frac{\partial D_n\left(s,t\right)}{\partial t}+\delta \left(s-s_i\right)·Q{}_i-\delta \left(s-s_e\right)·Q{}_e\end{array}$$

(11)

where n represents the unit normal vector perpendicular to fracture (s). For the solution of the above equation, the injection rate $Q_i\left(t\right)$ or injection pressure $p_i\left(t\right)$, and production rate $Q_e\left(t\right)$ or production pressure $p_e\left(t\right)$ can be prescribed as the boundary conditions, depending on the nature of the problem. In this study, we used a constant production pressure $p_e\left(t\right)$ at the wellbore as a boundary condition to simulate parent well fracture deformation due to production, whereas a constant injection rate $Q_i\left(t\right)$ was used to simulate the propagation of the infill well fractures that gave rise to frac-hits as they propagated. So, we have two different problems; each was solved according to its own boundary conditions. The boundary condition at the wellbore (i.e., the fluid flux equal to the injection rate) is already included in the above equation as the fluid injection rate. For the unique solution of the above equation, an additional boundary condition was required at the fracture fronts, which was usually assumed to be a no-flow boundary condition or zero fluid flux condition [34]. Reservoirs with high in-situ stress, which are prevalent in hydrocarbon reservoirs, have almost zero lag between the fluid front and the crack front; hence, the two fronts effectively coalesce into a single front denoted by the crack front [35]. Hence, the assumption of a zero-fluid flux boundary condition at the fracture front is justified for unconventional reservoirs. Equation (11) can be used to simulate the fluid flow inside the “mechanically open” fracture and propagating fractures (i.e., infill well fractures in this case) and “mechanically closed” previously propagated fractures (i.e., production well fractures in this case).

3. Numerical Methodology

In this study, we use GeoFrac-3D for the numerical simulation of the porous rock mass deformation and pore fluid diffusion into the rock matrix and the fracture propagation. The GeoFrac-3D model comprises two main components: a 3D boundary element method-based model to simulate solid rock deformation in an infinite domain and fluid flow in the porous rock matrix, and a finite element model for the fluid flow simulation inside the fractures. Galerkin’s finite element method is used for the numerical solution of incompressible fluid flow inside the fractures [36,37]. We used the linear elastic fracture mechanics approach (LEFM) for fracture propagation. A brief description of the numerical implementation procedure for each component is presented in the following sections.

3.1. Coupled Solution of Solid Rock Deformation, Pore Fluid Diffusion, and Fluid Flow in a Fracture

The rock matrix and fracture deformation, and fluid diffusion into the rock matrix are solved together using the 3D poroelastic displacement discontinuity (DD) method, which is an indirect boundary element method (BEM). A fracture in a poroelastic medium can be simulated by distributing displacements and discontinuous fluid flux, and by applying the principle of superposition to sum their effects. Using the corresponding singular solution, boundary integral equations for the resultant displacements, fluid fluxes, stresses, and reservoir pore pressure can be derived. The boundary integral representations for the stresses and pore pressure at any point in the reservoir or on the fracture surface can be expressed as (Curran and Carvalho [38]; Cheng [25]; Ghassemi and Zhou [33]; Safari and Ghassemi [39]; Kumar and Ghassemi [40]):

$${\sigma }_{ij}(\mathbf{x},t)-{\sigma }_{ij}^0(\mathbf{x})={\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{A}{\int }\left\{{\sigma }_{ijkn}^{id}(\mathbf{x}-\mathbf{x}{}^{\prime },t-\tau )D_{kn}(\mathbf{x}{}^{\prime },\tau )+{\sigma }_{ij}^{is}(\mathbf{x}-\mathbf{x}{}^{\prime },t-\tau )D_f(\mathbf{x}{}^{\prime },\tau )\right\}}}dA(\mathbf{x}{}^{\prime })d\tau$$

(12)

$$p(\mathbf{x},t)-p_0(\mathbf{x})={\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{A}{\int }\left\{p_{ij}^{id}(\mathbf{x}-\mathbf{x}{}^{\prime },t-\tau )D_{ij}(\mathbf{x}{}^{\prime },\tau )+p^{is}(\mathbf{x}-\mathbf{x}{}^{\prime },t-\tau )D_f(\mathbf{x}{}^{\prime },\tau )\right\}dA(\mathbf{x}{}^{\prime })d\tau }}$$

(13)

where x = (x,y,z) and $\chi$ = (x’,y’,z’) represent the local coordinates of the source and field points, respectively, t is the current time, $\tau$ is the time when the location $\chi$ receives the fluid first time, $D_{kn}$ represents the DD vector, $D_f$ is the fluid source intensity vector,${ \mathsf{\sigma }}_{ijkn}^{id}$,${ \mathsf{\sigma }}_{ij}^{is}$,$p_{ij}^{id}$, and $p^{is}$ are the instantaneous fundamental solutions for the stresses and pore pressure due to a unit impulse of the displacement discontinuity “id ” in the “kn” direction and a unit impulse of the fluid source intensity “is”, ${\sigma }_{ij}^0\left(\mathbf{x}\right)$ are the in situ stress components, and $p_0\left(\mathbf{x}\right)$ is the in situ reservoir pore pressure.

The Galerkin’s finite element method approach is used for discretization of the fluid flow equation. The FEM discretization of Equation (11) can be written as (Ghassemi et al. [34]):

$${\displaystyle \sum _{e=1}^E{\displaystyle \underset{A_e}{\int }{\nabla }^T{\mathsf{\phi }}^{(e)}\nabla {\mathsf{\phi }}^{(e)}dA\cdot \tilde{\mathbf{p}}}}+{\displaystyle \sum _{e=1}^E{\displaystyle \underset{Ae}{\int }{\mathsf{\phi }}^{T(e)}\nabla {\mathsf{\phi }}^{(e)}dA\cdot {\tilde{D}}_f}}=-{\displaystyle \sum _{e=1}^E{\displaystyle \underset{Ae}{\int }\left(\frac{D_n^{(k)}-D_n^{(k-1)}}{\Delta t}\right)}{\mathsf{\phi }}^{T(e)}dA}+{\displaystyle \underset{\partial A_p}{\int }{\mathsf{\phi }}^{(i)}{Q}_i(t)ds}-{\displaystyle \underset{\partial A_p}{\int }{\mathsf{\phi }}^{(p)}{Q}_e(t)ds}$$

(14)

where E is the number of the elements on the fracture plane and $A_e$ represents the area of an element “e”, ${\mathsf{\phi }}^{\left(e\right)}$ are the shape functions, $\tilde{p}$ and ${\tilde{D}}_f$ are the vectors of nodal pressure and fluid source intensity, respectively,$D_n^{\left(k\right)}$ and $D_n^{\left(k-1\right)}$ are the fracture mechanical opening at the current time step (k) and previous time step (k − 1), respectively, and $\Delta t$ is the time increment. The detailed mathematical and numerical procedure for the implementation of the coupled 3D poroelastic DD model and fluid flow inside the fracture can be found elsewhere [21,34,39,40,41]. The coupled solution of fracture aperture, fluid flux, and fluid pressure inside the fracture presents a highly nonlinear problem, which needs to be solved in a coupled manner. We used a fully coupled solution scheme as suggested by Ghassemi et al. [34] and modified by Kumar and Ghassemi [40,42].

Both hydraulic and natural fracture problems can be simulated in the GeoFrac-3D model. The numerical procedure to simulate production fracture deformation using the joint element model is similar to that used by Zhou and Ghassemi [41], Safari and Ghassemi [39], and Kamali and Ghassemi [43,44]. The infill well fracture propagation is simulated using the elastic module [42,45]. For hydraulic fractures, the effective stress on the fracture surface becomes zero as the two surfaces of the crack separate from each other. That is, the total normal stress on the fracture surface is equal to the fluid pressure. Additionally, the hydraulic fracture surfaces will have zero shear stresses, since there is no contact between the surfaces. Hence, for a hydraulic fracture, the total stresses on the fracture surfaces are given as (Safari and Ghassemi [39]):

$${\sigma }_n=-p_f; {\sigma }_{s1}=0 ; {\sigma }_{s2}=0$$

(15)

where ${\sigma }_n$ is the total normal stresses on the fracture surface, and ${\sigma }_{s1}$ and ${\sigma }_{s2}$ represent the total in-plane and out-of-plane shear stresses, respectively. Whereas, a different approach is used for the closed fractures or joint deformation simulation, as the joint normal and shear stiffness are nonzero and the effective normal and shear stresses on a joint surface change with its normal and shear displacements. For a linear elastic joint, the total normal and shear stresses on the fracture surface can be given as (Crouch and Starfield [46]; Safari and Ghassemi [39]):

$${\sigma }_n=-k_n{D}_n-p_f; {\sigma }_{s1}=-k_s{D}_{s1} ; {\sigma }_{s2}=-k_s{D}_{s2}$$

(16)

where $k_n$ and $k_s$ are the normal and shear stiffness of the joint filling material, respectively, $D_n$ represents the normal displacement discontinuity component, and $D_{s1}$ and $D_{s2}$ represents the in-plane (mode-II) and out-of-plane (mode-III) shear displacement discontinuity components, respectively. In this study, the positive values of the normal displacement discontinuity represent fracture opening, whereas the negatives values correspond to the fracture closure. Using the normal displacement discontinuity component or fracture mechanical opening, the fracture hydraulic aperture for the fluid flow Equation (11) is updated as:

$$w^k=w^{k-1}+\Delta D_n^k$$

(17)

where $w^{\left(k\right)}$ and $w^{\left(k-1\right)}$ represents the fracture hydraulic aperture at current time step (k) and previous time step (k − 1), respectively, and $\Delta D_n^{\left(k\right)}$ represents the change in the normal displacement discontinuity. By updating fracture aperture using Equation (17), Equation (11) can be used to simulate fluid flow either in the hydraulic (i.e., mechanically open fracture) or the natural fracture (i.e., mechanically closed fracture or joint).

3.2. Fracture Propagation

Dynamic fracture propagation is implemented in the framework of the linear elastic fracture mechanics approach (LEFM), based on the Griffith’s energy balance theory for the failure of brittle materials (Griffith [47]). Irwin [48] introduced a modification of Griffith theory and suggested that the stress field in the vicinity a crack front can be quantified in terms of three global parameters, i.e., stress intensity factors (SIF): opening mode or Mode-I, in-plane shear displacement or Mode-II, and out-of-plane shear displacement or Mode-III. Based on the Irwin’s criterion, the LEFM postulate that a crack will advance when the stress intensity factor at the crack front reaches a value equal to fracture toughness, which is a material property. The modified maximum circumferential stress criterion for the 3D fracture propagation, as suggested by Schöllmann et al. [49], is used. Based on the near crack tip stresses, the maximum principal stress is defined as ([49]):

$${\sigma }_1=\frac{{\sigma }_{\theta }+{\sigma }_z}{2}+\frac{1}{2}\sqrt{{\left({\sigma }_{\theta }-{\sigma }_z\right)}^2+4{\tau }_{\theta z}{}^2}$$

(18)

where ${\sigma }_{\theta }$, ${\sigma }_z$, and ${\tau }_{\theta z}$ are components of the stress tensor in a cylindrical coordinate system defined at the crack front. The expressions for these stresses and the maximum principal stress (${\sigma }_1$) in the cylindrical coordinates system $\left(r, \theta ,\mathsf{\phi }\right)$ in terms of the stress intensity factors can be found in Schöllmann et al. [49] and Richard et al. [50]. Schöllmann’s criterion also assumes that ${\sigma }_z$ has no contribution to the kink angle (i.e., ${\sigma }_z$ is assumed equal to zero). Based on the assumption that a crack grows perpendicular to the maximum principal stress, the fracture propagation angle (${\theta }_0$) is estimated by maximizing the maximum principal stress (${\sigma }_1$) as:

$${\frac{\partial {\sigma }_1}{\partial \theta }|}_{\theta ={\theta }_0}=0 ; {\frac{{\partial }^2{\sigma }_1}{{\partial }^2\theta }|}_{\theta ={\theta }_0}=0$$

(19)

The expressions for the first and second order derivatives of the maximum principal stress (${\sigma }_1$) in Equation (18) with respect to fracture propagation angle can be found in Schöllmann et al. [49]. Since there is no closed-form solution for the above equation, the fracture propagation angle (${\theta }_0$) needs to be estimated numerically using a root finding technique. In this study, the Newton-Raphson method is used to calculate the fracture propagation angle. Using the maximum principal stress (${\sigma }_1$) as defined in Equation (18) with the assumption of stress component ${\sigma }_z=0$, Schöllmann et al. [49] defined an equivalent stress intensity factor as:

$$K_{eq}=\frac{1}{2}cos\frac{{\theta }_0}{2}.\left\{K_{I}co{s}^2\frac{{\theta }_0}{2}-\frac{3}{2}{K}_{II}sin\frac{{\theta }_0}{2}+\sqrt{{\left[K_{I}co{s}^2\frac{{\theta }_0}{2}-\frac{3}{2}{K}_{II}sin\frac{{\theta }_0}{2}\right]}^2+4K_{III}{}^2}\right\}$$

(20)

The fluid injection time increment is adjusted in such a manner that at least one node on the fracture front should satisfy the local equilibrium condition, that is the equivalent mode-I stress intensity factor ($K_{eq}$) is equal to the mode-I fracture toughness ($K_{IC})$ within a prescribed tolerance limit (i.e., $K_{eq}\approx K_{IC}$). Once this local equilibrium condition has been achieved, the crack tip nodes which satisfy the equilibrium condition are propagated using scaling, as suggested by Gupta and Durate [51].

One of the important aspects of a hydraulic fracture simulation is to incorporate the time dependent variation of the fracture geometry. The fracture geometry is an essential parameter for simulation affecting the stability, accuracy, robustness, and the computational time. The adaptive remeshing of the fracture geometries is a critical component of the hydraulic fracture simulation process, since the final fracture geometry is much larger than that of the initial geometry. The configuration or the arrangement of the elements will be significantly changed as the fracture propagates. The three-dimensional fracture problem with moving fracture fronts requires an effective remeshing scheme for problem stability and cost-effective simulations. A remeshing scheme for the unstructured quadrilateral element mesh was development using the GMSH open source code (Geuzaine and Remacle [52]).

4. Stress Reorientation and “Frac-Hits”

In this section, the concepts of stress reorientation/reversal due to reservoir depletion and “frac-hits” in horizontal well refracturing are discussed using numerical examples.

4.1. The Concept of Stress Reorientation/Reversal

Production from a hydraulic fracture redistributes the pore pressure in an ellipsoidal region around the fracture surfaces. Figure 1a shows a plan-view of nonisotropic induced reservoir pore pressure changes around a production fracture after 2 years. The reservoir rock properties are as follows: Young’s modulus = 53.57 GPa, Poisson’s ratio = 0.29, reservoir permeability = 3.0 micro-Darcy, reservoir porosity = 4.5%, and Biot’s effective stress coefficient = 0.65. The in situ values are: minimum horizontal stress ${\sigma }_h=56.0 \mathrm{MPa}$, maximum horizontal stress${\sigma }_H=57.5 \mathrm{MPa}$, vertical stress${\sigma }_V=73.12 \mathrm{MPa}$, and reservoir pore pressure p = 50.9 MPa. The pore pressure reduction accompanies stress state variations in the vicinity of the fracture. In fact, the horizontal stress component parallel to the production fracture (in this case,${\sigma }_H$) decreases more than the horizontal stress component perpendicular to the fractures (i.e.,${\sigma }_h$ in this case), which is associated with a relatively higher reservoir pore pressure gradient in this direction. This induced-stress anisotropy modifies principal stress magnitudes and orientations. If the induced stress changes are large enough to overcome the initial horizontal stress difference (i.e.,${\sigma }_H-{\sigma }_h)$, complete stress reorientation/reversal occurs. The coupling between thermo-poroelastic processes is also important [15,33,34,39], because the thermally induced stresses and pore pressures also impact the reservoir stress state. If the production well is refractured after the stress reversal has occurred, the fracture will tend to propagate perpendicular to the initial propagation direction until the boundary of the stress reversal zone (see, Figure 1d) is reached. Beyond this zone, the new hydraulic fracture would re-align with the original maximum horizontal stress direction. The stress reversal zone evolves as a function of pore pressure depletion and varies with the reservoir’s poroelastic properties and the in situ stress differential.

4.2. The “Frac-Hits” in Horizontal Well Refracturing

The coupled poroelastic processes control the frac-hits phenomenon. Depending on the reservoir initial stress state, a frac-hit can occur mainly in two ways: (a) fracture-to-well interference, where the infill well fractures interact with the production well, as shown in Figure 2a; and (b) fracture-to-fracture interference, where the infill well fractures interact with the production well fractures, as shown in Figure 2b. For both examples, the reservoir rock properties are as follows: Young’s modulus = 53.57 GPa, Poisson’s ratio = 0.29, reservoir permeability = 3.0 micro-Darcy, reservoir porosity = 4.5%, and Biot’s effective stress coefficient = 0.65. The initial production fracture has a uniform aperture equal to 2.0 (mm and fracture height and half-length are equal to 60.0 (m) and 90.0 (m), respectively. The bottomhole pressure is equal to 12.0 MPa. For the first case, consider production from the production well in a normal faulting stress regime (i.e.,${\sigma }_V$ > ${\sigma }_H$ >${\sigma }_h$) with an in situ horizontal stress contrast equal to 3.42 (MPa). The in situ values are taken as the minimum horizontal stress ${\sigma }_h=54.3 \mathrm{MPa}$, the maximum horizontal stress${\sigma }_H=57.72 \mathrm{MPa}$, vertical stress ${\sigma }_V=73.12 \mathrm{MPa}$, and reservoir pore pressure p = 50.6 MPa. The distribution of reservoir pore pressure and the local maximum principal stress trajectories around production well fractures after 2 years of production are shown in Figure 2a. Note that because of significant reorientation of the local maximum principal stress near the tips of the production well fractures, the infill well fractures will curve to avoid the production well fracture tips. Beyond the production well fracture tips region, the infill well fractures will turn towards the production well to the reduced stress zones (Path-“A” in Figure 2a) and potentially communicate with the production well.

For the second case, we consider production from the production well fractures in a strike-slip faulting stress regime (i.e., ${\sigma }_H$ > ${\sigma }_V$ > ${\sigma }_h$) with a relatively high in situ horizontal stress difference of 32.9 (MPa). The in situ values are taken as the minimum horizontal stress ${\sigma }_h=41.80 \mathrm{MPa}$, the maximum horizontal stress ${\sigma }_H=74.70 \mathrm{MPa}$, vertical stress ${\sigma }_V=65.80 \mathrm{MPa}$, and reservoir pore pressure p = 32.5 MPa. Figure 2b shows the reservoir pore pressure distribution and the local maximum principal stress trajectories around production well fractures after 2 years of production. In this case, because of the large in situ horizontal stress differential, the local maximum principal stresses are not significantly reorientated (see, Figure 2b), and the infill well fractures will tend to propagate along the direction of maximum horizontal stress. However, the infill well fracture will experience low resistance and grow asymmetrically towards the production well fractures (due to the reduced stress zones around the parent well fractures). As a result, the infill well fractures may communicate with the production well fractures, as indicated by Path-“B” in Figure 2b. In what follows, we investigate these scenarios by dynamic simulations.

5. Model Verification

The GeoFrac-3D model was verified for penny-shaped fractures in poroelastic media subject to either uniform normal or shear traction by Zhou and Ghassemi [53]. Ghassemi et al. [34] verified the model by comparing pressurized penny-shaped fracture widths and induced normal stress with analytical solutions. Safari and Ghassemi [54] validated the model by analyzing the pressure history data in an injection/extraction (huff and puff) experiment in a geothermal field in Germany, and found good agreement between their numerical results and field measurements. Kumar and Ghassemi [45] verified a propagating penny-shaped fracture in elastic media. The same authors [40] also compared the results for the propagation of a penny-shaped fracture (both in elastic and poroelastic media) against the asymptotic solution in the viscosity-toughness transition regime as given by Savitski and Detournay [55], which characterizes the fracture propagation regime based on a parameter termed “dimensionless fracture toughness” (κ) (i.e., a parameter defined using rock properties and fluid injection parameters). For completeness, two additional verifications of the GeoFrac-3D model are presented in the following sections. The following reservoir rock properties were used for both verification examples: Young’s modulus = 53.57 GPa, Poisson’s ratio = 0.29, reservoir permeability = 3.0 micro-Darcy, reservoir porosity = 4.5%, and Biot’s effective stress coefficient = 0.65.

5.1. Penny-Shaped Fracture Propagation in Toughness Dominated Regime

Savitski and Detournay [55] defined a dimensionless fracture toughness ($\kappa )$ parameter to characterize the propagation regime, which is given as:

$$\kappa =K^{{}^{\prime }}{\left(\frac{t^2}{{{\mu }^{{}^{\prime }}}^5{Q}_{inj}^3{E^{{}^{\prime }}}^{13}}\right)}^{1/18}$$

(21)

where the parameters in the above equation are defined as:

$${\mu }^{{}^{\prime }}=12\mu ; E^{{}^{\prime }}=\frac{E}{\left(1-\nu {}^2\right)}; K^{{}^{\prime }}=4{\left(\frac{2}{\pi }\right)}^{1/2}K{}_{IC}$$

(22)

where μ is the fluid viscosity, E is Young’s modulus, ν is the Poisson’s ratio, $K_{IC}$ is the mode-I fracture toughness, t is the fluid injection time, and $Q_{inj}$ represents the fluid injection rate. The cases in which $\kappa <<1$ represent the viscosity-dominated regime, whereas those where $\kappa >>1$ are the toughness-dominated regime. Detournay [56] demonstrated the dependence of the propagation regime on dimensionless toughness ($\kappa )$, and categorized a viscosity-dominated regime if $\kappa <1$, and a toughness-dominated regime if $\kappa >3.5$. In Equation (21), the dimensionless fracture toughness is a time-dependent parameter; hence, the propagation behavior moves from a viscosity-dominated regime towards a toughness-dominated regime over time. Consider a penny-shaped fracture propagating in an impermeable reservoir rock. The fluid injection rate is 0.01 (m³/s) and fluid viscosity is equal to 0.001 (Pa.s). In this case, a high fracture toughness value (i.e., K_Ic = 10.0 MPa.m^0.5) is assumed to ensure that the fracture is propagating in the toughness-dominated regime. The distributions of the fracture aperture and fluid pressure at the end of injection (i.e., after 1300 (s) of injection) are shown in Figure 3a,b, respectively. In this case, since the fracture is propagating in a strong toughness dominated regime, the resultant fluid pressure is almost constant (i.e., variation of approximately 2 kPa from the injection well to the fracture tip), which is consistent with the asymptotic solution of the fracture fluid pressure put forward by Savitski and Detournay [55]. Comparisons of time variation of the dimensionless fracture toughness, fracture aperture at wellbore, fracture radius, and the net injection pressure from the GeoFrac-3D model and analytical solution are shown in Figure 4a–d, respectively. The temporal variation of dimensionless fracture toughness in Figure 4a demonstrates a toughness dominated propagation regime (i.e., $\kappa >3.5)$. The numerical results of fracture aperture at the wellbore (Figure 4b), and fracture radius (Figure 4c), and the net injection fluid pressure (Figure 4d) from the GeoFrac-3D model show very close agreement with the analytical solutions. The maximum relative errors in the fracture width and radius are in the range of 2.5–3%, whereas the fluid injection pressure variation shows an error of less than 0.5%.

5.2. Injection Induced-Pore Pressure Change around a Fracture

Kumar and Ghassemi [40] derived the numerical solution for the induced stress change due to the short-term pressurization (i.e., for 30 (s) of pressurization) and long-term pressurization (i.e., for 1 year of pressurization) of a circular fracture, and verified the accuracy of GeoFrac-3D results against the analytical solution presented by Sneddon [57]. To verify the reservoir pore pressure change due to fluid diffusion, an additional case of a pressurized 3D flat elliptical fracture is presented here. The change of reservoir pore pressure in the direction perpendicular to the fracture surface is approximately given as [58,59]:

$$P\left(\lambda ,\epsilon \right)=P_0+\left(P_f-P_0\right)e^{-\xi \lambda \left(1+\sqrt{\epsilon +b^2}\right)}$$

(23)

where ɛ and λ represent the ellipsoidal coordinates, P_f is the average pressurization pressure, and P₀ is the reservoir in situ pore pressure. The leak-off factor ξ is defined as:

$$\begin{gathered} \xi ={\left(\frac{g\pi }{2}\right)}^2 \\ g=\sqrt{\frac{\varphi \mu c}{kt}} \end{gathered}$$

(24)

where ϕ represents the reservoir porosity, µ is the fluid viscosity, c is the fluid compressibility, k is the reservoir permeability, and t is the pressurization time. Consider an elliptical fracture aligned in the YZ-plane with a half-length equal to 90.0 (m) and a height equal to 60.0 (m). The reservoir rock properties and in situ conditions are listed in

Table 1. Physical properties of the reservoir rocks in different layers.

Parameter	Layer-1	Layer-2	Layer-3	Layer-4	Layer-5
Layer thickness, d (m)	3.35	8.5	26.2	33.8	28.0
True vertical depth, h (m)	2696.22	2699.57	2708.07	2734.27	2773.07
Young’s modulus, E (GPa)	57.28	56.88	53.57	53.64	60.1
Poisson’s ratio (ν)	0.29	0.29	0.29	0.29	0.29
Permeability, k (μd)	0.07	1.84	3.04	2.96	1.84
Porosity, ϕ (%)	2.1	3.6	4.5	4.0	3.3
Biot’s coefficient (α)	0.65	0.65	0.65	0.65	0.65

and

Table 2. Rock mechanical properties and in-situ stresses.

Parameter	Value
Fluid viscosity, μ (Pa.s)	0.00025
Fluid compressibility, $c_f$ (MPa−1)	5.0 × 10−4
Vertical stress (${\sigma }_V$) gradient (MPa/m)	0.0266
Max. horizontal stress (${\sigma }_H$ ) gradient (MPa/m)	0.021
Min. horizontal stress (${\sigma }_h)$ gradient (MPa/m)	0.0198
Reservoir pore pressure (p) gradient (MPa/m)	0.0184
Fracture toughness, $K_{IC}$ (MPa.m0.5)	2.0
Fracture normal stiffness,$k_n$ (GPa/m)	300.0
Fracture shear stiffness,$k_s$ (GPa/m)	3000.0

. A constant uniform fluid pressure of 55.0 (MPa) is applied to the fracture surface. The change of reservoir pore pressure in the direction perpendicular to the fracture surface and its comparison with the analytical solution are shown in Figure 5a–d for pressurization times of 1 month, 6 months, 1 year, and 2 years, respectively. It can be observed that the numerical results from the GeoFrac-3D model show close agreement with the analytical solution, with a maximum error of 3.5%, and that the fluid diffusion front is moving further with pressurization time, as shown by the “double-sided arrow” mark.

6. Reservoir Rock Properties and In-Situ Stresses

The numerical examples in the following sections use the reservoir rock properties and in situ stresses for a shale formation having five layers, as shown in Figure 6. The rock mechanical and reservoir properties of the layered reservoir are summarized in Table 1 [5]. The reservoir rocks in each layer are assumed to be homogeneous and isotropic with uniform permeability. For numerical simulations using the DD model, we used the layer thickness-weighted average reservoir properties of all the layers [60], as the DD model is based on homogeneous rock mass assumption which cannot account for the properties of different layers, whereas the FEM model accounts for distinct layer properties, as given in Table 1. It is important to mention here that as the objective of this study is to investigate the impact of reservoir depletion on the pore pressure and stress changes using the field scale production fractures (i.e., fracture half-lengths equal to 90.0 m and height equal to 60.0 m), the localized impact of the drilling-induced reservoir pore pressure and in situ stresses change near the wellbore are not considered at this time. Therefore, the initial conditions for all the examples are the same as the in situ reservoir conditions listed in Table 2. The production well has a true vertical depth (TVD) of 2743.07 m.

7. Impact of Reservoir Layers Properties

The GeoFrac-3D model, being a boundary element method-based model, works only for a homogeneous infinite reservoir with constant rock properties. However, unconventional reservoirs usually have a finite thickness with different over- and under-burden formations. To study the impact of the layered structure on the stresses and pore pressure changes due to depletion, we carried out a finite element simulation using our state-of-the-art coupled Galerkin’s finite element model (FEM), which accounts for distinct layer properties. A 3D poroelastic FEM model was developed and used to simulate the time-dependent distributions of pore pressure and stresses. The model is based on the poroelastic governing equations, i.e., Equations (7) and (8). By applying the Green-Gauss theorem and Galerkin’s weighted residual method to equilibrium and mass balance equations, the weak form of these equations is obtained [61,62,63]. The implicit finite difference method was used for time discretization, which is first-order accurate and unconditionally stable. The system of algebraic equations thus obtained are solved for values of field variables at each node of the finite element mesh. The numerical mesh is highly refined near the fractures to resolve the strong pore pressure gradients in these locations. A stimulated reservoir volume (SRV) region is assumed around the production fractures. This is a common approach (e.g., [64]) when using a finite element method without a zero-thickness element to describe the fracture zone.

In the FEM simulations, each fracture is explicitly modeled with an unstructured fine grid so that the model is capable of capturing the fracture flow and the large pressure gradients involved. Linear 8-node hexahedral elements are used for domain discretization. An equivalent fracture conductivity equal to 2.0 (md.ft) is taken for the production fractures which, with a fracture aperture of 2.0 (mm), translates to an equivalent fracture permeability of 304.0 (md). Isothermal, single-phase, slightly compressible fluid is assumed. Due to symmetry in the YZ- and XZ-planes, only 1/4th of the geometry is modeled, and the plane of symmetry boundary conditions are applied on the XZ- and YZ-planes, respectively. The initial conditions are set assuming a linear pore pressure gradient of 0.0184 (MPa/m). For the mechanical boundary conditions, a minimum horizontal stress gradient ${\sigma }_h$ = 0.0198 (MPa/m) and a maximum horizontal stress gradient ${\sigma }_H$ = 0.021 (MPa/m) are applied on far-field boundaries, as shown in Figure 7c,d. In the vertical direction, the model extends from 2443.07 (m) to 3043.07 (m) in depth. Constant vertical traction is applied on the top face of the domain to account for the overburden load, which is equal to lithostatic pressure at a depth of 2443.07 (m). The zero normal displacements are applied to the two planes of symmetry and the bottom boundary. For the fluid flow, a pressure boundary condition is applied at the wellbore with a bottomhole pressure of 10.5 (MPa) prescribed at each well node intersecting the fractures. A no-flow boundary condition is applied at all the remaining boundaries. A similar approach for the boundary conditions was followed in Gray et al. [65], Feng et al. [66], and Rutqvist et al. [67]. To eliminate any potential boundary effects on the reservoir pore pressure and stresses in the study area, the boundaries of the FEM model are kept sufficiently far from the production fractures; the modeling domain dimensions are 1200 m × 1200 m × 600 m, as shown in Figure 7c,d while the production fracture half-length and height are 90 m × 60 m. A stimulated reservoir volume (SRV) is modeled with dimensions of 175 m × 230 m × 60 m, and a permeability 15 times higher than the reservoir layer permeabilities is assigned to the SRV region to capture the enhanced permeability in the region next to propped fractures. In this simulation example, an optimum configuration for the SRV region is considered where all the areas between hydraulic fractures implicitly contain activated fracture networks.

Comparisons of the reservoir pore pressure and stresses distribution from the FEM model (with layered reservoir properties) and from the DD model (with equivalent reservoir properties) are shown in Figure 8 and Figure 9, respectively. The distributions from both methods are in close agreement. Comparisons of the reservoir pore pressure and the horizontal stress components from the DD and FEM models after 2 years of production along a line in the central XY-plane (dotted “red line” in Figure 7b) are shown in Figure 10a–c, respectively. The maximum difference in pore pressure is about 6.50%, and stress variations show close agreement. Hence, using equivalent reservoir properties in the DD model is an acceptable alternative with a significant reduction in computational effort. The FEM simulations were performed on the supercomputer using 1 node and a total of 24 cores (2.3 GHz) using 64 Gb of RAM, as reservoir volume discretization is required. The discretization of the FEM model resulted in 852,640 8-node hexahedral elements, and it took around 4 h to run the code. In contrast, the DD model simulation was performed on a Dell Precision T7600 (Quad-core processor, 2.30 GHz) desktop machine and took around 45 min. In the DD model, only the fracture surface discretization is required to get unknown variables such as displacement discontinuities and fluid source intensities. Once these unknown variables are obtained for the fracture surfaces, the change in the reservoir pore pressure and stresses are calculated in the postprocessing. In the DD model, for fracture surface discretization, we used a total of 900 4-noded quadrilateral elements and 16,000 field points for the reservoir pore pressure and stress change calculations.

8. Numerical Examples

In this section, we present a simulation of the undepleted zones among closely spaced horizontal wells and subsequent fracture propagation from the infill well in the altered stress state. The numerical simulations are carried out in two stages. For the first stage, reservoir depletion analysis for a certain period of production from the primary wells is considered, and in the second stage, the infill well fracture propagations are simulated. The fracture propagation from the infill well is considered in two scenarios before and after the repressurization of the primary well fractures and in two different stress regimes (i.e., the normal faulting and strike-slip faulting stress regimes).

8.1. Reservoir Depletion Due to Production from Parent Wells

A 3D schematic of two horizontal production wells, each having five stages with one cluster per stage, along with the proposed location of the infill well, is shown in Figure 11. Each cluster is assumed to have one dominating fracture. The dash-dotted “blue line” shows the potential infill well hydraulic fracture propagation path. Consider the in situ stress state of the reservoir in the normal faulting regime. The minimum horizontal (${\sigma }_h)$ and maximum horizontal (${\sigma }_H$) in situ stresses are acting along the x- and y-axes, respectively, and the vertical in situ stress (${\sigma }_V$) is acting along the z-axis. The production and infill wells are drilled parallel to the minimum horizontal stress (${\sigma }_h)$ direction (i.e., along the x-axis in this case) and stacked in the horizontal plane. The production fractures are assumed to be planar and orthogonal to the wellbore axis (i.e., parallel to the YZ-plane). The production fractures are assumed to behave like joints or natural fractures while closing, to account for the presence of proppant. The fracture half-length and height are assumed to be equal to 90.0 (m) and 60.0 (m), respectively. The spacing between the horizontal wells is equal to 400.0 (m), and spacing among stages is assumed to be equal to 90.0 (m). The production fractures are assumed to be fully propped, with the initial uniform fracture aperture equal to 2.0 (mm), which is a justified value in the subsurface, as its value typically ranges from 0.1 (mm) to 3.0 (mm) [68]. In a recent numerical study, Sesetty and Ghassemi [69] showed that the fracture aperture value for hydraulic fractures in unconventional reservoirs has a range of 1.9–3.0 (mm). Both production wells are producing at the constant bottomhole pressure equal to 10.5 (MPa). The in situ stresses and mechanical properties of the reservoir rocks are listed in Table 2. The numerical simulations are carried out in two stages: (a) reservoir depletion analysis; and (b) subsequent infill well fracturing. These two stages of the simulation are presented in the next sections. For the first stage of modeling, the production from the production well fractures is carried out at a constant bottomhole pressure (BHP) of 10.5 (MPa) for 2 years. In this study, the positive sign of stresses represents compression and the negative sign represents tension. The depletion-induced reservoir pore pressure change and mechanical closure of the production fractures will affect both the total and effective reservoir stresses; these parameters are studied next.

8.1.1. Impact of Reservoir Depletion on the Total Stresses

The reservoir pore pressure and the total stress change after 2 years of production are shown in Figure 12a–d, respectively. Due to symmetry, the reservoir pore pressure and the total stress change are shown only for one-half of the reservoir volume in the vicinity of the parent well fractures. The reservoir pore pressure distribution is shown in Figure 12a, which demonstrates an ellipsoidal depletion zone formed around the production fractures. Note that the reservoir pore pressure is decreased to 10.5 (MPa) in the depletion zone. The impact of reservoir depletion on the total horizontal stress component ${\sigma }_{xx}$ is shown in Figure 12b, indicating a decrease of 11.3 (MPa) (i.e., from the in situ value of 54.30 to 43.0 MPa) in the vicinity of parent well fractures. The impact of reservoir depletion on the total horizontal stress component ${\sigma }_{yy}$ is shown in Figure 12c, which shows that its value has decreased to 46.8 (MPa) in the depletion zone. Figure 12d shows that the vertical stress ${\sigma }_V$ decreased to 65.60 (MPa) from 73.12 (MPa) in the depletion zone. In general, the vertical stress is a function of density of the overburden layers, and is independent of the reservoir pressure. However, the production from a hydraulic fracture redistributes the local pore pressure in an ellipsoidal region around the fracture surfaces, and hence, impacts reservoir stress components in all directions. Therefore, we can observe a reduction in the vertical stress also in this case. Again, in boundary element methods, the boundary is at infinity, so there would be no effects. A boundary, as in earth’s surface, may have some effect; however, the depth of the simulation is relatively large, so it is, in effect, infinite. Depending on the reservoir initial stress state, depletion-induced stress changes might alter the stress regime. A check for the status of the stress regime after 2 years of production time is carried-out. The differences between the vertical stress and the minimum horizontal stress (${\sigma }_V-{\sigma }_h$), and the maximum horizontal stress (${\sigma }_V-{\sigma }_H$) are shown in Figure 13a,b, respectively. Positive values of the difference between the total vertical and horizontal stresses show that the reservoir stress state is still within a normal faulting stress regime (i.e., ${\sigma }_V>{\sigma }_H>{\sigma }_h).$

Another perspective of the reservoir pore pressure and the total stresses change is presented using a horizontal cross-section in the central XY-plane, as shown in Figure 14. The reservoir pore pressure distribution and the local maximum principal stress trajectories after 2 years of production are shown in Figure 14a. The pore pressure distribution plot shows an elliptical depletion zone around the production fractures. A significant reorientation of the stress state in the vicinity of production fractures occurred. The total horizontal stress components (${\sigma }_{xx}$) and $({\sigma }_{yy}$) are decreasing from the in-situ values around the production fractures (see, Figure 14b,c). The deformation of the parent well fractures is simulated in a fully coupled solution system; hence, the stress shadowing impact from the outer fractures on the inner fractures is already accounted for in this study. This stress shadowing impact can be observed in the stress component ${\sigma }_{xx}$ plot in Figure 14b, where the inner fractures show relatively higher stress changes compared to the outer fractures. The plot of horizontal stress contrast ($\Delta \sigma ={\sigma }_{yy}-{\sigma }_{xx}$) in Figure 14d highlights the stress reversal. Negative stress contrast values (dark blue regions in Figure 14d) show stress reversal zones, where the magnitude of the maximum horizontal stress became less than the minimum horizontal stress.

8.1.2. Impact of Reservoir Depletion on the Effective Stresses

The impact of reservoir depletion on the effective stresses are presented in this section. A 3D visualization of the effective horizontal stress components ${\sigma }_{xx}^{{}^{\prime }}$ and ${\sigma }_{yy}^{{}^{\prime }}$, and the effective vertical stress ${\sigma }_{zz}^{{}^{\prime }}$ are shown in Figure 15a–c, respectively. Due to reservoir depletion, both the pore pressure and the total stress components decrease; however, the rate of pore pressure decrease is higher than the total stress reduction. Hence, the resultant effective stress components (${\sigma }_{ij}^{{}^{\prime }}={\sigma }_{ij}-{\delta }_{ij}\alpha p)$ are increasing. For example, in this case, after 2 years of production, the maximum value of the effective horizontal stress component${\sigma }_{xx}^{{}^{\prime }}$ has increased to 37.0 (MPa) from its in-situ value, i.e., 21.41 (MPa), the maximum value of the effective horizontal stress component ${\sigma }_{yy}^{{}^{\prime }}$ has increased to 37.6 MPa from its in-situ value, i.e., 24.83 (MPa), and the maximum value of effective vertical stress ${\sigma }_{zz}^{{}^{\prime }}$ has increased to 56.6 (MPa) from its in-situ value, i.e., 40.23 (MPa). The increase of the effective stresses results in an increase in reservoir rock displacements around the parent well fractures. The reservoir rock displacements in the x- and y-directions are shown in Figure 16a,b, respectively. In Figure 16c, the dotted “black” lines show the resultant displacement directions. From Figure 16c, it can be observed that the reservoir rock is compressed in both the x- and y-directions.

8.1.3. Impact of Production Time on the Reservoir Pore pressure and Total Stresses

For this example, consider the impact of production time on the reservoir pore pressure and stresses change. The reservoir pore pressure distributions and trajectories of the local maximum principal stress (“dotted black lines”) after 1 month, 1 year, 2 years, and 5 years of production times are shown in Figure 17a–d, respectively. It can be observed that with production time, the fluid diffusion front advances; after approximately 2 years of production time, we can notice horizontal stress reversal near the infill well zone.

Figure 18 shows the temporal variation of the total horizontal stress components ${\sigma }_{xx}$ and ${\sigma }_{yy}$ along the infill well axis in the central XY-plane. With production time increase, the reservoir pore pressure decreases, resulting in a decrease of both horizontal stress components. However, the stress component parallel to the production fractures (${\sigma }_{yy})$ decreases at a higher rate compared to the horizontal stress component perpendicular to the production fractures (${\sigma }_{xx})$. This contrast in the reduction of horizontal stresses will result in the reorientation of the principal stresses, which may lead to the development of stress reversal zones, as indicated in Figure 18d. It should be noted that in this case, the reservoir pore pressure and stress changes occur due to the combined effect of the mechanical closure of the production fractures and the depletion-induced poroelastic effect. From Figure 18a, it can be observed that in the early production time (i.e., 1 month), the stress component perpendicular to the parent well fractures ${\sigma }_{xx}$ has increased from its in situ value near the production fracture tips, and later, is showing peak and valleys behavior. Since the mechanical closure of the production hydraulic fractures is considered in this study using a linear joint model, the combined effect of stress shadowing among the production fractures and the mechanical closure of the fracture tips causes an increase in compressive stress near the tips. As a result, the horizontal stress component ${\sigma }_{xx}$ has increased near the production fractures tips. In Figure 18c, it can be seen that after 2 years of production time, the horizontal stresses ${\sigma }_{xx}$ and ${\sigma }_{yy}$ have approximately similar magnitudes near the infill well location, which suggests that beyond 2 years of production time, stress reversal will occur in this zone.

To investigate the impact of depletion-induced stress on subsequent infill well fracture propagation paths, consider Figure 19, where the time variations of the reservoir pore pressure, the total horizontal stress components ${\mathsf{\sigma }}_{xx}$ and ${\mathsf{\sigma }}_{yy}$, and the horizontal stress contrast in the central XY-plane and parallel to the parent well fractures are presented (i.e., along the dotted “blue line” in Figure 11b). It can be observed that with production time, an anisotropic reduction in the reservoir pore pressure is occurring (see Figure 19a), which causes a heterogeneous reduction of the total horizontal stresses (see Figure 19b,c). However, the stress component parallel to the production fractures (${\sigma }_{yy})$ decreases at a higher rate compared to the stress component perpendicular to the production fractures (${\sigma }_{xx})$. The horizontal stress contrast value (i.e., Figure 19d) decreases with production time, which will promote stress reorientation and may lead to a reversal. As a result of this reorientation, the infill well fractures will turn away from the parent well fracture tips.

8.2. Fracture Propagation from the Infill Well

The reservoir pore pressure and stress analysis in the previous section has demonstrated that production from the parent wells creates a “pressure sink” or reduced stress zone near the production fractures. These reduced stress zones will create an attraction zone for the subsequent infill well fractures, and an asymmetric propagation of the hydraulic fractures from the infill well will occur. Two main reasons for the preferred growth of the infill well fractures toward the parent well are: firstly, the reduced reservoir pore pressure of the depleted zone causes most of the fracturing fluid to flow towards the reduced pressure zone; and secondly, the reservoir pore pressure depletion causes a decrease in the total stresses, which facilitates paths of least resistance for fracture propagation. However, in this case, the reduction of total reservoir stresses is the dominating factor which leads to preferred growth of the infill well fracture towards the parent wells. A single fracturing stage with one fracture cluster with a single propagating fracture from the infill well has been simulated in this case. It is important to mention that to simulate the fracture propagation from the infill well, we used the elastic module of the model, and the poroelastic module was turned off. For all infill well fracturing simulations, a constant fluid injection rate equal to 0.05 (m³/s) is used and the fracturing fluid is slick water with a viscosity equal to 0.001 (Pa.s). The infill well fractures are also fully contained in the pay zone, similar to production well fractures. To study the impact of in situ stress regime, the infill well fracture propagations before and after the repressurization of the parent well fractures are simulated both in the normal faulting and the strike-slip stress regime conditions.

8.2.1. Infill Well Fracture Propagation in the Normal Faulting Stress Regime

For the first example, consider infill well fracturing in a normal faulting stress regime after 1.5 years of production. The initial stress state and reservoir pore pressure for this case are listed in Table 2. The created fracture geometry and its location in the reservoir after 1.5 years of production are shown in Figure 20a,b, respectively. In this case, the infill well fracture shows more growth towards the production well due to reduced stress zones and turning away from the parent well fractures (see Figure 20b). In this case, due to the reorientation of the local maximum principal stress near the parent well’s fracture tips, the infill well fractures will curve and avoid the parent well’s fracture tips. Beyond the parent well fracture tips region, the infill well fractures will again turn towards the parent well due to the presence of reduced stress zones, as indicated by “Path-A” in Figure 2a, which potentially leads to communication with the production wells.

For the second case, consider infill well fracture creation in the normal faulting stress regime after 3 years of production. The created fracture geometry and its location in the reservoir after 3 years of production are shown in Figure 21a,b, respectively. As the in-situ horizontal stress contrast, in this case, has a relatively low value of 3.42 (MPa) at the wellbore level, after 3 years of production, we can notice a complete stress reversal at the infill well area (see, Figure 21b). In this case, due to the modified stress state, longitudinal fractures will be created. As can be seen from Figure 21a, a nonplanar fracture propagation initiates from the infill well which tends to align with the infill well axis (i.e., along the x-axis in this case, which is the modified orientation of maximum horizontal stress). It should be noted that the infill well fracturing in this situation will not lead to frac-hits. However, the hydrocarbon production from these infill well fractures will be significantly less than from the corresponding transverse fractures, since these fractures will have limited contact surface areas in the reservoir. Hence, for the horizontal well refracturing, there is a time window in which to optimize production from the infill fractures.

8.2.2. Infill Well Fracture Propagation in the Strike-Slip Stress Regime

For the second example, consider the infill well fracturing in a strike-slip stress regime. In this case, the maximum horizontal stress (${\sigma }_H$) gradient is equal to 0.0271 (MPa/m), the minimum horizontal stress (${\sigma }_h$) gradient is equal to 0.0151 (MPa/m), the vertical stress (${\sigma }_V$) gradient is equal to 0.0238 (MPa/m), and the reservoir pore pressure gradient is equal to 0.0122 (MPa/m). The created fracture geometry and its location in the reservoir are shown Figure 22a,b, respectively. The in situ horizontal stress contrast, in this case, has a relatively high value equal to 32.9 (MPa) at the wellbore level; hence, after 3 years of production, there is not any significant stress reorientation in the infill well area (see, Figure 22b). In this case, the infill well fracture shows asymmetric growth towards the parent wells due to the creation of reduced stress zones. Because of depletion, stress decreases around the production wells and the infill well fracture experiences reduced stress zones and shows more growth towards the production wells. It should be noted that infill fractures also show asymmetric growth in the vertical direction due to the presence of in situ stress gradients in the reservoir. This asymmetric growth will potentially lead to communication with the parent wells, and most probably, will result in interference of the infill well fracture with the parent well fractures, as demonstrated in Figure 2b.

8.2.3. Infill Well Fracture Propagation in the Strike-Slip Stress Regime

For the second stage of the simulation analysis, consider the repressurization of the production well fractures before the fracturing of the infill well. In this case, the production well fractures are pressurized for 72 h at a constant bottomhole pressure of 45.0 (MPa). To avoid any potential propagation of the parent well fractures, this repressurization pressure is kept below the stress component normal to the fracture surface (i.e., less than the minimum horizontal stress, equal to 54.3 MPa in this case). Due to repressurization, the production fractures expand and the fluid from the fractures diffuses into the reservoir, causing an increase in the reservoir pore pressure and stresses. The reservoir pore pressure, total horizontal stress components ${\sigma }_{xx}$ and ${\sigma }_{yy}$, and the vertical stress ${\sigma }_{zz}$ are shown in Figure 23a–c and Figure 24d, respectively. It can be observed that the reservoir pore pressure and stress state have approximately returned to their original in situ values. In this updated stress state, the infill well fractures will not experience any reduced stress zones and will not show preferred growth towards the parent wells. Hence, there will be a lower risk for a frac-hit to occur. The created fracture geometry and its location in the reservoir for the normal faulting stress regime are shown in Figure 24a,b, respectively. The infill well fracture shows symmetric growth towards production wells and asymmetric growth in the vertical direction (i.e., z-axis) due to the in situ stress gradient (see, Figure 24a). Figure 24a shows that the maximum fracture width is not at the well where the pressure is maximum, owing to the presence of an in situ stress gradient which increases with the formation depth; therefore, the upper part of the fracture experiencies relatively lower in situ stresses, and hence, reports higher fracture apertures in this region. The results demonstrate that by the repressurization of the production well fractures before infill well fracturing, well communication and frac-hits problem can be mitigated.

9. Conclusions

A geomechanical analysis of well-to-well interference or frac-hits in horizontal well refracturing is presented using a fully coupled 3D poroelastic simulator model GeoFrac-3D. The model can simulate both the hydraulic and natural fracture deformation in the elastic and poroelastic reservoirs. The poroelastic DD model used in this study is a computationally efficient numerical technique for porous rock deformation and matrix fluid flow simulations. The DD model reduces the problem dimensionality by one, and only the fracture surface (s) discretization is required. A comparison of the depletion induced reservoir pore pressure and stresses change from the GeoFrac-3D model with average reservoir properties and the FEM model with layered reservoir properties shows close agreement. Hence, the DD model, which is computationally less expensive, can be used to effectively simulate reservoir depletion and refracturing analyses, even for the layered reservoirs. The results show that depletion-induced pore pressure and stress state changes have strong effects on subsequent infill well fracturing. The production from the production well fractures creates a local “pressure sink” or stress reduction zones, which leads to a high risk of frac-hits problem or well-to-well communication. It is observed that if we repressurize the production well fractures before infill well fracturing, the risk of frac-hits can be mitigated. The infill well fracture simulation examples show that this aspect could be improved if carried out between the start of production and the onset time of stress reversal. The simulation results demonstrate that the in situ reservoir stress regime and horizontal stress anisotropy have a strong impact on the infill well fracturing. The model presented in this study can simulate both aspects of refracturing simulation, i.e., reservoir depletion analysis and subsequent infill well fracture propagation modeling, using a single model. In contrast, most of existing refracturing simulation models use two separate models, i.e., one for the reservoir depletion analysis and the second for the subsequent infill well fracturing.