Hydraulic Fracturing in Formations with Permeable Natural Fractures

The recently developed Unconventional Fracture Model (UFM*) simulates complex hydraulic fracture network propagation in a formation with pre-existing closed natural fractures, and explicitly models hydraulic injection into a fracture network with multiple propagating branches [1]. The model predicts whether a hydraulic fracture front crosses or is arrested by a natural fracture it encounters, which defines the complexity of the generated complex hydraulic fracture network.

While taking into account the leakoff of the fracturing fluid into the formation, the leakoff into the natural fractures should also be considered, especially in low-matrix permeability conditions. The transmissibility of natural fractures can become significant, and the fracturing fluid can penetrate into natural fractures. Different regions can coexist along the invaded natural fracture: hydraulically opened region filled with fracturing fluid, region of still closed natural fracture invaded by fracturing fluid due to natural fracture permeability, and the region of natural fracture filled with original reservoir fluid.

Explicit modelling of hydraulic fractures interacting with permeable natural fractures becomes extremely complicated with the necessity to account for conservation of fluid mass, pressure drop along natural fractures, leak-off into the formation from natural fracture walls, pressure sensitive natural fracture permeability, properties of natural fractures, fluid rheology, while tracking the interface of each region along invaded natural fracture. A main challenge is integrating this hydraulic fracture/natural fracture interaction modelling into the overall hydraulic fracture network propagating scheme without losing model effectiveness and CPU performance.

The updated UFM model with enhancement to account for leakoff into the natural fractures will be presented.

1. Introduction

It is believed that complexity of the resulting fracture network during hydraulic fracturing treatments in formations with pre-existing natural fractures is caused mostly by the interaction between hydraulic and natural fractures. Natural fractures can be important for hydrocarbon production in the majority of low-permeability reservoirs, particularly where the permeability of the rock matrix is negligible. Understanding and proper modelling of the mechanism of hydraulic-natural fractures interaction is a key to explain fracture complexity and the microseismic events observed during HF treatments, and therefore to properly predict production.

When hydraulic fracture (HF) intercepts natural fracture (NF) it can cross the NF, open (dilate) the NF, or be arrested at NF. If hydraulic fracture crosses natural fracture, it remains planar, with a possibility to open the intersected NF if the fluid pressure at the intersection exceeds the effective stress acting on the NF. If the HF does not cross the NF, it can dilate and eventually propagate into the NF, which leads to more complex fracture network.

The interaction between HF and NF depends on in-situ rock stresses, mechanical properties of the rock, properties of natural fractures, and hydraulic fracture treatment parameters including fracturing fluid properties and injection rate. During the last decades, extensive theoretical, numerical, and experimental work has been done to investigate, explain, and use the rules controlling HF/NF interaction [3-17]. A new crossing model [2] recently implemented in UFM is able to predict the crossing behaviour of HF at the NF accounting for the effects of fluid properties and NF permeability [18].

One of the important effects of natural fractures is enhanced leakoff, which can lead to a premature screenout during proppant injection. In a formation with low-matrix-permeability, the transmissibility of natural fissures can be significantly higher than that of the reservoir matrix. The fracturing fluid can readily penetrate into natural fissures during the fracturing process and maintain a pressure nearly equal to the pressure in the primary fracture [19].

The concept that natural fractures (fissures) could alter leakoff has been a subject of numerous studies [20-24] with considering the fissure opening conditions or pressure-sensitive leakoff conditions. It was often reported that permeability of natural fractures is pressure dependent [23, 25, 26].

The ways that elevated pressure could affect natural fractures have been described in [27]. Fissures with rough surfaces and minimal mineralization are most likely highly sensitive to the net stress pushing on them. Under virgin reservoir conditions (when the pressure p within the fissure equals the initial reservoir pressure p _ini ), the effective stress is fairly high and the open channels formed by mismatched fracture faces are most likely deformed and nearly closed. As the pressure in the fissure increases because of leakoff of the high-pressure fracturing fluid (p > p _ini ), the net closure stress is reduced and the fissure porosity opens. In this regime, the leakoff coefficient is highly pressure dependent. As the pressure exceeds the closure stress on the fissure (p > p _fo ), the entire fissure opens, yielding an accelerated leakoff condition. The estimation of the critical pressure in PKN-type HF (exceeding closing pressure p > p _fo ) to open a vertical fissure has been given through the function of the principal horizontal stresses and Poisson ratio [20].

A more detailed description of the effects from natural fissures in reservoirs where natural fissures are the primary source of permeability is provided in [23]. The enhanced rate of fluid loss throughout the treatment is predicted, with leakoff accelerating as the fracturing pressure increases. The increase in fluid pressure in the fissures reduces the effective normal stress acting to close the fissures and hence increases their permeability. For hydraulic fracturing purposes, the effect of the magnified permeability is reflected as an increase in fluid-leakoff coefficient. The fluid-leakoff in the presence of natural fissures could be as high as 2 to 3 times that for normally occurring pressure – dependent leakoff behaviour, even under the net pressure conditions.

For slightly elevated pressures NF porosity begins to open as the pore pressure increases because the elevated pressure relieves some of the net stress on the asperity contacts. Several models of this process have been developed. For example, [26] predicts the change in NF permeability resulting from changes in stress and pressure. This model have been validated and used in numerous studies [23].

Among existing HF models accounting for the permeability of intercepted natural fractures mention [10] which couples fluid flow, elastic deformation, and frictional sliding to obtain a solution which depends on the competition between fractures for the permeability enhancements. The effect of initially closed but conductive fracture is specifically addressed. The possible scenarios for evolution of fracture opening and fluid transport in closed NFs implemented in [10] are shown in Figure 1. The initial aperture w ₀ along a closed pre-existing NF corresponds to its residual conductivity. It is equal to the effective aperture for the parallel plate model. The initial conductivity of a closed natural fracture arises from the fact that its surfaces are rough and mismatched at fine scale, i.e. the aperture w ₀ is related to the fracture porosity. With increasing the fluid pressure, the hydraulic aperture will slightly change due to micro structural change in the natural fracture, although the fracture still remain closed and carries some contact stresses. In the end, fracture will be opened mechanically as the fluid pressure exceeds the normal stress acting on the fracture. In this case, the effective hydraulic conductivity is equal to the sum of both hydraulic aperture and mechanical opening since the fracture opening augments the initial hydraulic aperture, as shown in Figure 1. Zhang’s model also considers the possibility of frictional sliding through the Coulomb frictional law, and accounts for three types of contact behaviour at fracture surface: fracture is opened, fracture is closed but surface is in sticking mode, and fracture is closed but in sliding mode.

Figure 1.

Evolution of natural fracture opening [10]

The HF models [28-32] do not account for permeability of natural fractures explicitly. The 2D model in [33] uses approach from [13] to simulate interaction between induced propagating fracture and natural fracture. A modified leak-off model for an intersecting fracture based on poro-elasticity was introduced to account for the increased leakoff at the intersections. A poro-elastic solution for the stresses in the HF/NF interaction zone has been used as a basis for hydraulic/natural fracture interaction criteria. A fully coupled finite element based approach was used to simulate HF propagation in a poroelastic formation with existing natural fractures.

The approach given in [10] is based on boundary element method and rigorously models HF interaction with permeable NF. It is computationally expensive, and is applicable for analysis of limited (small) number of HF/NF interactions. For a more general complex fracture network model like UFM which deals with a large number (order of thousands) of natural fractures, the CPU time is important and model should be computationally efficient while still being physically correct.

The important aspect of HF/NF interaction is shear slippage of NF faces. The possibilities of shear slippage in natural fractures due to change of stress field (in isolated natural fractures) or during HF/NF interactions, and the influence of shear slippage on fracture aperture change and dilation have been a subject of experimental and numerical studies [10, 15, 17, 30, 34-37]. The conditions for shear slippage and the corresponding shear displacement (apertures) have been investigated [36], and the estimation of permeability of NF with changing effective normal stress is done by [11,38]. The shear slippage effect during HF/NF interaction is also included in current approach.

This paper describes how leakoff into the natural fractures during HF/NF interaction (crossing or arresting before NF opens) is integrated into the complex hydraulic fracture model UFM.

2. UFM model specifics

A complex fracture network model, referred to as Unconventional Fracture Model (UFM), had recently been developed [1, 39, 40]. The model simulates the fracture propagation, rock deformation, and fluid flow in the complex fracture network created during a treatment. The model solves the fully coupled problem of fluid flow in the fracture network and the elastic deformation of the fractures, which has similar assumptions and governing equations as conventional pseudo-3D fracture models. Transport equations are solved for each component of the fluids and proppants pumped. A key difference between UFM and the conventional planar fracture model is being able to simulate the interaction of hydraulic fractures with pre-existing natural fractures, i.e., determine whether a hydraulic fracture propagates through or is arrested by a natural fracture when they intersect and subsequently propagates along the natural fracture.

To properly simulate the propagation of multiple or complex fractures, the fracture model takes into account the interaction among adjacent hydraulic fracture branches, often referred to as “stress shadow” effect. It is well known that when a single planar hydraulic fracture is opened under a finite fluid net pressure, it exerts a stress field on the surrounding rock that is proportional to the net pressure. The details of stress shadow effect implemented in UFM are given in [40].

The branching of hydraulic fracture when intersecting natural fracture gives rise to the development of a complex fracture network. A crossing model that is extended from the Renshaw-Pollard [12] interface crossing criterion, applicable to any intersection angle, has been developed, validated against the experimental data [16, 17], and was integrated at first in the UFM. The crossing model, showing good comparison with existing experimental data, did not account for the effect of fluid viscosity and flow rate on the crossing pattern. More recently a new advanced OpenT crossing model, taking into the account the impact of fluid and NF properties, have been developed [2] and integrated in UFM [18].

The modelling approach used in UFM to predict the leakoff into the NFs is presented below.

3. Modeling leakoff into permeable NF in UFM

The main assumptions for the current modelling of leakoff from HF into the intercepted NF in UFM are given below.

• The rock formation contains vertical discrete deformable fractures (NFs), which are initially closed but conductive because of their pre-existing apertures (due to surface roughness, etc).

• The propagation direction of the hydraulic fractures is not affected (unless intercepted) by closed and not invaded natural fractures. The intercepted NF could affect HFN propagation even from closed parts (when shear slippage takes place)

• The natural fractures are assumed to contain pore space and are permeable.

• The original fluid inside the natural fractures and in the reservoir (oil, gas, water) is compressible and Newtonian.

• Fracturing fluid can be incompressible or compressible and its rheology can be Newtonian or power law.

• The rock material is assumed to be permeable and elastic.

• When intercepted by the main hydraulic fracture, a natural fracture may remain closed, while still being able to accept fracturing fluid, or may be mechanically opened by fracturing fluid pressure depending on the magnitude of the fracturing fluid pressure, confining stresses applied on natural fracture, and frictional properties of natural fracture.

• The flow inside NFs is assumed to be 1D.

• The natural fractures can be opened by fluid pressure that exceeds the normal stress acting on them and/or experience Coulomb type frictional slip.

• The original natural fracture has width w ₀ and are filled with reservoir fluid with pressure equal to pore pressure p ₀ =p _res .

• Fluid flow invaded into the natural fracture develops along NFs. Invaded fracturing fluid into the NF may reach the end of NF, break the rock, and start to propagate into the rock accordingly to previously implemented propagation rules (only if the NF is opened). Fracture re-initiation from other points along the NF other than its ends (offsets) is not modeled at this time.

Hydraulic fractures propagating in the rock are modeled in accordance with existing approach in UFM model. The Schematic of the complex HF interaction with permeable NF is shown in Figure 2 and Figure 4.

Figure 2.

Hydraulic fracture intercepting natural fracture and possible situation to model

Fracturing fluid invasion into the two wings of the NF needs to be considered separately. Four possible regions can co-exist in each wing of the NF encountered by HF (Figure 4):

1. Opened part filled with invaded fracturing fluid (fluid pressure exceeds the normal effective stress on NF), with length of opened part L _opened >0

2. Invaded closed part of NF (filtration zone) filled with fracturing fluid (fluid pressure above pore pressure but below the closure stress) with length L _filtration >0

3. Closed pressurized part filled with pressurized original reservoir fluid (fluid pressure above the pore pressure) with length L _pressurized >0

4. Closed undisturbed part of NF filled with reservoir fluid under original pore pressure conditions.

When a natural fracture is intercepted by the hydraulic fracture, the fluid pressure in the hydraulic fracture transmits into the natural fracture. If the fluid pressure is less than the normal effective stress on the natural fracture, the natural fracture remains closed. Even closed natural fractures may have hydraulic conductivities much larger than the surrounding rock matrix, and in this case fracturing fluid will invade the natural fractures more than leakoff into the surrounding matrix. If the portion of injected fluid is lost into closed natural fractures from the main HF, the HF growth could be affected.

For a closed fracture, the equivalent fluid conductivity is expected to change with the fluid pressure since contact deformation is a function of effective normal stress. This pressure-induced dilatancy and the associated increase in conductivity are important in increasing leakoff. Also, any reduction in effective contact stress may result in fracture sliding, which can lead to local stress variations and slip induced fracture dilation, which can in turn change the overall conductivity of fracture networks.

The governing processes in first three regions listed above should be modeled, and the modeling approaches in different regions (also referred to as zones in the following context) are different due to different flow behaviors and rock/fluid properties.

4. Basic governing equations

4.1. Continuity of fluid volume (mass)

The equation for the continuity of incompressible fluid volume has the form

$$\frac{\partial q_{NF}}{\partial s}+\frac{\partial A}{\partial t}+q_L=0,q_L=\frac{2h{C}_{tot}^{rock}}{\sqrt{t-\tau (s)}},A=\varpi h$$ E1

where

q _NF (t) - volumetric flow rate through a cross section of area A of natural fracture [m³/s]

A - cross sectional area of the natural fracture

q _L - the volume rate of leakoff per unit length

$$\varpi$$ - average hydraulic fracture width (different from w, fracture opening by fluid pressure exceeding normal stress)

h - fracture height

C _tot - total leakoff coefficient from the wall of natural fracture

More generally in the case of compressible fluid the equation (1) should account for fluid density ρ _f and mass flux q _m . Considering the rate of change of fluid mass per unit length in a fracture $$\dot{m}$$ , continuity of fluid mass in the fracture is governed by equation

$$\frac{\partial q_m}{\partial s}+\dot{m}+{\rho }_f{q}_L=0$$ E2

or along the fracture of constant length

$$\frac{\partial q_m}{\partial s}+\frac{\partial ({\rho }_f\varpi h)}{\partial t}+{\rho }_f{q}_L=0,q_L=\frac{2h{C}_{tot}^{rock}}{\sqrt{t-\tau (s)}}$$ E3

Here s is coordinate along NF, and total leakoff coefficient from the walls of the natural fracture $$C_{tot}^{rock}$$ is equal to combined leakoff coefficient [41]

$$C_{tot}^{rock}=C_{vc}^{rock}=\frac{2C_v^{rock}{C}_c^{rock}}{C_v^{rock}+\sqrt{C_v^{rock2}+4C_c^{rock2}}}$$ E4

Where leakoff coefficient for the filtration zone in the rock and leakoff coefficient for the reservoir zone as shown in equations (5a - 5b)

$$\begin{array}{l}{C}_v^{rock}=\sqrt{\frac{k_r{\phi }_r\Delta p}{2{\mu }_f}}=\sqrt{\frac{k_r{\phi }_r\Delta p}{2{\mu }_f}},\Delta p=p_f-p_s\text{a}\\ C_c^{rock}=\sqrt{\frac{k_r{\phi }_r{c}_T}{\pi {\mu }_r}}\Delta p,\Delta p=p_f-p_r\text{b}\\ {\overline{C}}_v^{}=\frac{2{\left(C_v^{}\right)}^2\sqrt{t}}{V_L}\text{c}\end{array}$$ E5

with

ϕ _r - reservoir porosity

c _T - total compressibility of reservoir

k _r - permeability of rock matrix

µ _r - reservoir fluid viscosity in the porous media

µ _f - filtrate fluid viscosity

ρ _f - filtrate fluid density

p _r - reservoir pressure

In the case of multiple fluids, the invaded zone can be described by replacing $$C_{\nu }$$ with equivalent term (see (5c)) [42] Where $${\overline{C}}_v^{}$$ is calculated using the average viscosity and relative permeability of all the filtrate fluids leaked off up to the current time, and V _L is the fluid volume per unit area that previously leaked off into the reservoir.

4.2. Pressure drop along closed NF

The pressure drop along closed NF can be expressed from Darcy’s law

$$q_{NF}=-\frac{k_{NF}A}{{\mu }_f}\frac{\Delta p}{L(t)},A=\varpi h$$ E6

Or for mass flux

$$q_m=-{\rho }_f\frac{k_{NF}}{{\mu }_f}A\frac{\partial p}{\partial s}$$ E7

Then the pressure can be calculated from

$$\frac{\partial p}{\partial s}=-\frac{{\mu }_f}{{\rho }_f{k}_{NF}A}{q}_m=-\frac{{\mu }_f}{{\rho }_f{k}_{NF}\varpi h}{q}_m,\text{at the inlet : }p=p_{in}(t)$$ E8

Here

k _NF - permeability of natural fracture

µ _f - filtrate fluid viscosity

ρ _f - filtrate fluid density

A - cross sectional area of closed NF

p _in - fluid pressure at the inlet

4.3. Change of NF permeability due to stress and pressure changes

The leakoff into the natural fracture, or permeability of the natural fracture, is highly pressure dependent when pressure of invading fluid exceeds reservoir pressure but still is below the closure pressure. In general, permeability of natural fracture is a function of normal stress on NF, shear stress (or shear displacement due to shear slippage), and fluid pressure, and can be represented as a combination of permeability due to normal stress and permeability due to shear slippage [26]:

$$\begin{array}{l}{k}_{NF}=f(k_{NF}^n,k_{NF}^s)\\ k_{NF}^n=f_1(k_o,{\sigma }_n,p_{}),k_{NF}^s=f_2(u_s,{\varphi }_{dil})\\ k_{NF}^n=k_o{\left\{C\ln \left[\frac{\sigma *}{{\sigma }_n-p_{}}\right]\right\}}^3\end{array}$$ E9

Where constants C and σ ^* (reference stress state) are determined from field data, k _o is the initial NF permeability (reservoir permeability under in-situ conditions), σ _n is the normal stress on the NF, p is the pressure in the NF, and u _s is shear-induced displacement (slippage).

4.4. The width of closed invaded NF

The width $$\varpi$$ of closed NF invaded by the treatment fluid (hydraulic aperture) is related to the pressure-dependent permeability as [10]

$$\varpi =\sqrt{12k_{NF}}$$ E10

The hydraulic width $$\varpi$$ can be evaluated from Barton-Bandis model following approach [11,36,38], i.e. directly from Eq. (11) for given effective normal stress σ _eff , reference effective stress $${\sigma }_n^{ref}$$ , initial hydraulic fracture aperture $${\varpi }_o$$ (related to the roughness of fracture surface), shear displacement $$u_s$$ and dilation angle $${\varphi }_{dil}$$

$$\begin{array}{l}\varpi =\frac{{\varpi }_0}{1+9\frac{{\sigma }_{eff}}{{\sigma }_n^{ref}}}+{\varpi }_s+{\varpi }_{res},{\varpi }_s=\left|u_s\right|\tan ({\varphi }_{dil}^{eff}),{\varphi }_{dil}^{eff}=\frac{{\varphi }_{dil}^{}}{1+9\frac{{\sigma }_{eff}}{{\sigma }_n^{ref}}}\\ {\sigma }_{eff}={\sigma }_n-p_f(s)\end{array}$$ E11

4.5. Shear failure

The second term in Eq. (11) represents the shear-induced dilation which contributes to the NF permeability (Eq.9). Shear induced dilation is related to the frictional slip which occurs when the shear stress reaches the frictional shear strength of the natural fractures $${\tau }_s=\lambda ({\sigma }_n-p)$$ . In the present study the NF propagation due to the shear induced slip is not considered, but the contribution of shear slippage zone to the NF enhanced permeability is evaluated based on the enhanced 2D DDM approach [40,43,44].

$$\begin{array}{l}{\sigma }_n^i-p_{}^i={\displaystyle \sum _j{A}^{ij}{C}_{nn}^{ij}{D}_n^j}+{\displaystyle \sum _j{A}^{ij}{C}_{ns}^{ij}{D}_s^j}\\ {\tau }_{}^i={\displaystyle \sum _j{A}^{ij}{C}_{sn}^{ij}{D}_n^j}+{\displaystyle \sum _j{A}^{ij}{C}_{ss}^{ij}{D}_s^j,}{A}^{ij}=1-\frac{d_{ij}^{2.3}}{{\left(d_{ij}^2+h^2\right)}^{2.3/2}}\end{array}$$ E12

The fracture surface slip (shear displacement) $$u_s$$ can be found as shear displacement discontinuity $$D_s^{}$$ calculated for the closed sliding elements following Coulomb frictional law $$\tau \ge {\tau }_s=\lambda (\sigma -p_{})$$ from elasticity equations in the stress shadow calculation approach with accounting for the mechanical opening in HFN from Eq. (12).

Figure 3.

Stress Shadow Effect from opened (blue) HFN and closed parts (grey) of intercepted NFs

For any element j in the opened part of HFN (including opened part of intercepted NFs) the input normal displacement discontinuity $$D_n^j$$ in Equation (12) is given by known fracture aperture (width) $$D_n^j=w^j\ne 0$$ , and the shear stress is zero $${\tau }_{}^j=0$$ . Along the closed part of intercepted NFs the mechanical the opening is zero, $$D_n^j=0$$ , and the pressure and normal stress should be tracked to detect (find) elements sliding in shear. If $${\tau }^i\ge {\tau }_s^i$$ then element i is sliding and dilating in shear, and the shear stress $${\tau }^i={\tau }_s^i$$ for this element is used to find fracture surface slip $$u_s^i=D_s^i$$ from the second equation (12). So, equation (12) could be solved for shear displacements $$D_s^j$$ along opened and sliding in shear parts of total HFN and intercepted NF as schematically shown in Fig.3. In Eq.(12) C ^ij are 2D, plane strain elastic influence coefficients [43] defining interactions between the elements i and j, and A ^ij are 3D correction factors [44] accounting for the 3D effect due to fracture height h depending on the distance between elements d _ij .

4.6. Fluid density as function of pressure and temperature

Density of gas as a function of pressure and temperature has form

$${\rho }_{gas}=\frac{p{m}_m}{ZRT}$$ E13

Where p is pressure, T is temperature, Z is compressibility, R is gas constant, and m _m is molar mass. The changes in pressure and temperature with time produce changes in gas density. Density of a compressible fluid as function of pressure

$$\frac{d\rho }{dt}=\frac{\rho }{B}\frac{dp}{dt}=\rho c_f\frac{dp}{dt}$$ E14

where B is bulk modulus (fluid elasticity) in Pa and c _f =1/B is fluid compressibility in Pa^-1. More generally, change in fluid density due to changes in pressure and temperature for the low compressibility fluids through the known values of density ρ ₀ and temperature T₀ at pressure p ₀ has form (β here is volumetric expansion coefficient)

$$\rho =\frac{{\rho }_0}{\text{1-}\beta \text{(T}-T_0)}\times \frac{1}{1-\frac{p-p_0}{B}}$$ E15

4.7. Fluid flow in opened NF

The fluid flow in the opened NF (p _f >σ _n ) will be handled as fluid flow in HFN and have been described before [1] depending on the flow regime:

Laminar fluid flow: Poiseuille Law [45]

$$\frac{\partial p_{}}{\partial s}=-{\alpha }_0\frac{1}{{\overline{w}}^{2n\text{'}+1}}\frac{q}{H_{fl}}{\left|\frac{q}{H_{fl}}\right|}^{n\text{'}-1}$$ E16

$${\alpha }_0=\frac{2K\text{'}}{\phi {\left(n\text{'}\right)}^{n\text{'}}}{\left(\frac{4n\text{'}+2}{n\text{'}}\right)}^{n\text{'}};\phi \left(n\text{'}\right)=\frac{1}{H_{fl}}{\displaystyle \underset{H_{fl}}{\overset{}{\int }}{\left(\frac{w(z)}{\overline{w}}\right)}^{\frac{2n\text{'}+1}{n\text{'}}}}dz$$

Where $$\overline{w}$$ is average fracture opening, and n’ and K’ are fluid power law exponent and consistency index.

Turbulent fluid flow (N_Re>4000):

$$\frac{\partial p_{}}{\partial s}=-\frac{f{\rho }_f}{{\overline{w}}^3}\frac{q}{H_{fl}}{\left|\frac{q}{H_{fl}}\right|}^{},q=H_{fl}{\left(-\frac{{\overline{w}}^3}{f{\rho }_f}\frac{dp}{ds}\right)}^{1/2}$$ E17

With Reynolds number (N _Re) for the power law fluid between parallel plates and Fanning friction factor (f) defined as

$$N_{\mathrm{Re}}=\frac{3^{1-n\text{'}}{2}^{2-n\text{'}}\rho V^{2-n\text{'}}{\overline{w}}^{n\text{'}}}{K\text{'}{\left(\frac{2n\text{'}+1}{3n\text{'}}\right)}^{n\text{'}}},f\approx \frac{1}{16{\left[{\log }_{10}{\left(\frac{\epsilon }{7.4\overline{w}}\right)}^{}\right]}^2};$$ E18

Where V is fluid velocity and ε is surface roughness height.

Darcy fluid flow through proppant pack of height h

$$\frac{\partial p_{}}{\partial s}=-\frac{q{\mu }_{fl}}{kh{\overline{w}}^{}}$$ E19

will take place if the height of the fluid in fracture element become smaller than minimum fluid height. The minimum fluid height is calculated from the condition that pressure drop is equal to pressure drop due to the Darcy flow. The minimum height for turbulent and laminar flow is defined as

$$\begin{array}{l}\text{Laminar flow: } h_{fl}^{\min }={\left(\frac{k{\alpha }_0(n\text{'})H_{fl}}{{\overline{w}}^{2n\text{'}}{\mu }_{fl}}\right)}^{1/n\text{'}}\frac{1}{q^{\frac{1}{n\text{'}}-1}}\\ \text{Turbulent flow: } h_{fl}^{\min }={\left(\frac{k{H}_{fl}f\rho q}{{\overline{w}}^2{\mu }_{fl}}\right)}^{1/2}\end{array}$$ E20

Where $${\mu }_{fl}$$ is fluid dynamic viscosity, and k is proppant pack permeability. The boundary conditions at the inlet and tip of opened fracture

$$p=p_f(t),p_{tip}={\sigma }_n^{tip}(t)$$ E21

where p _f is known fluid pressure at the intersection with natural fracture.

5. Combined fluid flow into the opened and closed parts of invaded NF

As I mentioned before, natural fracture can be closed, closed but invaded with fracturing fluid, closed and filled with pressurized reservoir fluid, or opened (Figure 4). The partially opened NF can contain opened, invaded, pressurized and closed parts which are dynamically changing with time. When fronts (positions of the boundaries between co-existing parts in invaded NF) change or propagate, the velocity of each propagating front can be considered as velocity of the corresponding fluid front in the relevant part of the invaded NF.

To properly model invasion of fracturing fluid into the NF, the propagation of each front should be modelled, and in different parts of the invaded NF different governing equations should be satisfied.

Figure 4.

Details on different possible zones in intercepted permeable NF

5.1. Opened part of NF

If the NF is opened at intersection element i with HF, then fluid pressure at intersection exceed the local normal stress on NF:

$$p_f(i)>{\sigma }_n^{NF}(i)$$ E22

The tip of opened part of NF or the intersection element with HF (if NF is closed) becomes the inlet (injection point) into the closed invaded part of NF (filtration zone) $$p_{in}^{filtr}=p_{tip}^{open}={\sigma }_n^{NF}$$ . Mention that if NF is completely opened, then it is a part of total hydraulic fracture network (HFN) and handled based on the approach described previously in [1], see also equations (16)-(21).

If NF is closed at intersection element i with HF, then fluid pressure is below the local normal stress, but can still be higher than reservoir (pore) pressure

$$p_o<p_f(i)\le {\sigma }_n^{NF}(i)$$ E23

5.2. Filtration zone (closed part of NF invaded by fracturing fluid)

The fluid pressure, width and flow rate along the filtration zone can be calculated from given pressure $$p_{in}^{filtr}=p_f(i)$$ which satisfies condition (23). The flow rate along filtration zone can be iteratively solved from the system of equations (24) with flow rate from the inlet to filtration zone $$q_{in}^{filtr}$$ which is a part of total solution

$$\{\begin{array}{l}\frac{\partial p}{\partial s}=-\frac{{\mu }_f}{k_{NF}{\varpi }_{filtr}h}q,\text{at the inlet : }p=p_{in}^{filtr}(t)\\ \frac{\partial q}{\partial s}+\frac{\partial ({\varpi }_{filtr}h)}{\partial t}+q_L=0,q_L=\frac{2h{C}_{tot}^{rock}}{\sqrt{t-{\tau }_o(s)}}\\ k_{NF}=\frac{{\varpi }_{filtr}{}^2(p)}{12}\\ {\varpi }_{filtr}(p)=\frac{{\varpi }_0}{1+9\frac{{\sigma }_{eff}}{{\sigma }_n^{ref}}}+{\varpi }_s+{\varpi }_{res},{\varpi }_s=\left|u_s\right|\tan ({\varphi }_{dil}^{eff}),{\sigma }_{eff}={\sigma }_n-p(s),\end{array}$$ E24

The length of filtration zone can be calculated by the tracking the volume of fracturing fluid leaked into the NF by marching from the inlet along the NF. At the end of filtration zone mass balance should be satisfied and fluid pressure will be higher or equal to the reservoir pressure

$$\begin{array}{l}{q}_{in}^{filtr}dt=dVo{l}_{frac}^{filtr}(dt)+dVo{l}_{leak}^{filtr}(dt)\\ p_r<p_{}^{endfiltrzone}<{\sigma }_n^{NF}\end{array}$$ E25

The flowrate (filtration/pressurized front velocity) at the last element of filtration zone is used for calculations in the pressurized zone. The position of the front between the filtration zone and pressurized zone should be tracked, giving velocity of filtration front and pressure $$p_{in}^{pres}$$ as input for solution in pressurized zone.

If NF is partially opened, then in the solution scheme for the filtration zone the intersection (inlet) element is replaced by the tip element i of the opened part of NF with $$p(i)={\sigma }_n^{NF}(i)$$ .

5.3. Closed pressurized part of NF (filled with reservoir fluid).

Mention that if the leakoff coefficient for filtration zone $$C_{\nu }^{rock}=0$$ , then $$p_r=p_f^{endfiltrzone}<{\sigma }_n^{NF},$$ and there will not be pressurized zone in NF ( $$p_r=p_f^{endfiltrzone}<{\sigma }_n^{NF},L_{pressurized}=0$$ ).

For the general case of compressible reservoir fluid and non-zero total leakoff from the walls of pressurizes NF into the rock, the leakoff to the rock from the NF part filled with pressurized reservoir fluid is defined by the compressibility controlled leakoff coefficient

$$C_c^{rock}=\sqrt{\frac{k_r{\phi }_r{c}_T}{\pi {\mu }_r}}\Delta p,\Delta p=p_{NF}-p_r$$ E26

The governing equations to calculate fluid pressure, width, and flow rate along the pressurized zone from known influx q _filtr/pres and pressure p _filtr/pres are similar to Equations (24) for filtration zone with replacing fracturing fluid with reservoir fluid and using $${\varpi }_{pres}$$ as hydraulic width of pressurized zone of invaded NF, and c _T as reservoir fluid compressibility.

Notice that the pressure at the front between opened and filtration zones along invaded NF (or at HF/closed NF intersection point) and the time step are the inputs for the new pressure, width, and flow rate calculation in closed invaded part of intercepted NF. Initial influx q _in can be prescribed based on the pressure in NF from the previous time step, and then the solution scheme will be applied from the intersection towards the end of NF with tracking the incremental mass balance (fluid injected to NF at current time step) to define the end of filtration front, and tracking pressure in the rest of NF (not invaded part, filled with reservoir fluid) to track the end of pressurized zone (fluid pressure equal to reservoir pressure). The flow rate and pressure are tracked and corresponding front positions are to be updated iteratively until in the pressurized zone of NF the following condition is satisfied (which indicates the position of the end of the pressurized zone)

$$\begin{array}{l}q(L_{open}+L_{filtr}+L_{pres})=0\\ p(L_{open}+L_{filtr}+L_{pres})=p_r\end{array}$$ E27

At the end of pressurized zone pressure is equal to the reservoir pressure and flow rate is zero. If the end of NF is reached and pressure is above the reservoir pressure, then Equation (27) is replaced with condition q(L _NF )=0. During pumping the pressurized zone extends until the end of NF, and then gradually shrinks, while the lengths of invaded zone and opened zones increase. Eventually whole NF will be filled with fracturing fluid, so the NF will contain only filtration and/or opened zones. When NF is completely opened, it becomes a part of total HFN. The elements in closed part of invaded NF can slip under some conditions (for example, due to the stress field change from stress shadow), influencing the calculations of total NF permeability k _NF .

Each element in the closed part of intercepted NF is checked for shear slip possibility. The fracture surface slip (the shear displacement $$u_s$$ ) can be found as shear displacement discontinuity calculated for the closed sliding elements satisfying Coulomb frictional law $$\tau \ge {\tau }_s=\lambda (\sigma -p_f)$$ , from elasticity equations (12) accounting for the mechanical opening in HFN.

If some elements in closed not disturbed part of NF are sliding then the corresponding shear stress is involved in the stress shadow calculations and thus influences simulations results.

6. Numerical approach description

The treatment of permeable natural fractures is a part of UFM model. It is possible to model leakoff from HF into the NF form UFM in different ways, depending on the importance of required accuracy, numerical stability and CPU time.

The most computationally expensive but at the same time most accurate is a fully coupled numerical approach. The fully coupled approach means to discretise different parts of the invaded NFs and numerically solve the pressure as a part of total system of equations for the fracture network using iterative solution scheme. Another, more CPU efficient approach, is a decoupled numerical approach, where pressure and width along the invaded but closed part of NF are calculated separately based on the results from previous time step calculations along the HFN and corresponding pressure at HF/NF intersections (inlets).

Two approaches have been considered to model interaction of hydraulic fracture with permeable natural fractures.

6.1. Decoupled numerical approach

• When NF is intercepted by HF, create elements along whole NF to be used at next time step

• Evaluate initial guess of flow rate into the NF based on the pressure at the intersection and old pressure profile along closed NF

• Check for a possibility of frictional sliding along the closed parts of NFs to evaluate pressure-dependent permeability and conductivity along NF

• Iteratively calculate pressure, hydraulic width, flow rate and length of each zone along NF with using corresponding equations (for filtration and pressurized zones) by marching from intersection till the end of NF until condition (27) is satisfied indicating final results and final positions of zones’ fronts

• Track the end of filtration zone for each pressure and flow rate iterations by checking the volume of fluid injected into NF at given time step

• Save invaded volume for volume balance, influx for mass balance, pressure at intersection, and time step to be used at the next time step

• Track the pressure at intersections and/or tips of opened HF part in NF to capture opened zones.

• If intersection is opened, use the pressure at the tip of opened HFN part along invaded NF for calculations along corresponding closed NF part

• Apply rules to treat special situations (intersecting NFs, etc)

• Save elements information (volume, pressure, flow rates) for the next time step

6.2. Fully coupled numerical approach

Figure 5.

The proposed algorithm to account for leakoff from HF into the permeable NFs (i=N indicates the end of NF)

• Discretize the NF when HF intercepts it

• Make the NF elements a part of total network (HFN) and include pressure calculations along different zones of NF into the whole iterative scheme to calculate pressure, time step, and front positions.

• Calculate the flow rate and volume of fluid injected to NFs at each pressure iteration, and update/calculate pressure along HFN using existing scheme for opened elements and new equations (described above) for elements in closed invaded and pressurized NF parts

• From the calculated pressure in NF elements iteratively update the positions of the propagating fronts for opened, filtration, and pressurized zones in NF during pressure/time step iterations

• Track the pressure at intersections to capture when NF start to open

• Include the elements in the closed parts of invaded NFs into the stress shadow calculation scheme

The fully coupled NF modeling approach is heavier and more CPU expensive than decoupled approach. The Decoupled Numerical Approach has been selected as basic approach and it is described schematically on Figure 5 as a part of total solution

7. Conclusions

The new approach developed to account for the complex processes due to NF permeability accompanying HF/NF interaction have been presented in detail. This approach accounts for the important physical processes taking place during HF and permeable NF interaction, and will be implemented in UFM.

The next step is to evaluate the influence of leakoff into the natural fractures during HFN simulations on the total HFN footprint and production forecast. The approach will be also validated against existing numerical, experimental and field data.