Semi-Analytical Model for Two-Phase Flowback in Complex Fracture Networks in Shale Oil Reservoirs

Abstract

Flowback data is the earliest available data for estimating fracture geometries and the assessment of different fracturing techniques. Considerable attention has been paid recently to analyze flowback data quantitively in order to obtain fracture properties such as effective half-length and effective conductivity by simply assuming fractures having bi-wing planar geometries and constant fracture compressibility. However, this simplifying assumption ignores the complexity of fracture networks. To overcome this limitation, we proposed a semi-analytical method, which can be used as a direct model for fast inverse analysis to characterize complex fracture networks generated during hydraulic fracturing. A two-phase oil–water flowback model with a matrix oil influx for wells with bi-wing planar fractures is also presented to identify limitations of the former solution. Since most available flowback studies use constant fracture properties and the assumption of planar fractures, considering variable fracture properties and complex fracture geometries gives this model more robustness for modeling fracture flow during flowback, more realistically. The proposed models have been validated by numerical simulations. The presented procedure provides a simple way for modeling early flowback in complex fracture networks and it can be used for inverse analysis.

1. Introduction

Reservoir simulation is a powerful tool to predict reservoir performance under different operating conditions, however, when it comes to fractured wells, it requires detailed knowledge of induced fractures. Routinely, well testing and long-term rate-transient analysis (RTA) are used for post-fracture analysis to determine fracture properties [1,2]. Recent advances in field data acquisition and development of powerful data analytic techniques promise new venues for unlocking fracture properties such as fracture conductivity and complexity of induced fracture networks by analysis of flowback data if a good physical model for flowback is available. Flowback data is the earliest part of production data that operators may acquire to determine fracture properties for building reservoir models and assess the stimulation design.

The first attempt to quantitively analyze flowback data uses a single-phase reciprocal productivity Index (RPI) procedure [3]. This method provides an estimation of the effective permeability-thickness, apparent fracture half-length or effective wellbore radius for vertical wells. The model is derived based on the bi-wing planar fracture geometries and assumes that fracture volume does not change during flowback. Flowback production time in shale wells is divided into three periods: (1) the pre-breakthrough water dominated region, (2) the transition region, and (3) the post-breakthrough hydrocarbon dominated region [4]. The authors developed a single-phase flowback model for the first period only and assumed a linear correlation between rate-normalized pressure and material balance flowrate to calculate total storage coefficients and effective fracture permeability. A semi-analytical model for flowback in tight oil reservoirs uses the concept of a dynamic drainage area (DDA) to capture the transient behaviors [5]. A bi-wing planar fracture is assumed in the model. However, due to the presence of natural fractures in the formation and their interaction with advancing hydraulic fractures, the geometry of induced fractures can be very complicated [6]. The outcrop of cross-cutting joints found in the Marcellus shale exposed in Oatka Creek also shows well-developed natural fracture systems [7]. To accommodate the reactivation of pre-existing natural fractures, they defined an enhanced fracture region (EFR) along the primary fracture, inside which the permeability is enhanced. However, this approach indicates that overall permeability is enhanced in the direction perpendicular to the primary fracture, which misleads the true directionality of fracture permeability in complex fracture networks [8]. A complex fracture network is desirable from the design perspective as it increases the effective surface area of the wellbore and enable more-efficient linear flow to predominate over the radial flow. Production from reservoirs with complex fracture networks is often modeled by a numerical simulator, which can be computationally expensive for use in inverse analyses, unless some meshless methods are utilized for flow through fractures [9]. Most available flowback models incorporate fracture closure by using constant fracture compressibility for wells with bi-wing planar fractures [4,5,8,10]. However, the outcome can be misleading for two main reasons: (1) compressibility of propped fracture is expected to evolve with effective pressure based on the analytical solution for compliance of proppant beds. Especially during early flowback, rapid pressure drawdown may lead to significant decrease in fracture compressibility. Flowback analysis based on constant compressibility will overestimate the fracture storage capacity at the end of flowback. (2) Fracture width is decreasing from the wellbore to the tip, knowing that fracture compressibility is a function of proppant density, we do not expect that fracture compressibility remains constant. This assumption can introduce misleading impacts on determining fracture properties. To incorporate compressibility changes over time, we initially proposed analytical solutions for single and two-phase water–oil flowback using eigenfunction expansion methods. Then, we develop a semi-analytical procedure to incorporate variations of fracture width and fracture compressibility along the fracture, which can be utilized for complex fracture networks and long-planar fractures. We also proposed a two-phase oil-water flowback solution with matrix oil influx for wells with bi-wing planar fractures, which can be applied to formations without significant natural fractures. The proposed model is verified against numerical simulations implemented in COMSOL Multiphysics (COMSOL, Burlington, Massachusetts, USA) Multiphysics and CMG (Computer Modelling Group LTD., Calgary, Alberta, Canada). The goals of this paper are set to (1) provide a quick and accurate direct flowback model for fluid flow through the fracture network, (2) investigate the influence of non-uniform fracture closure and variable fracture compressibility on flowback characteristics.

2. Governing Equations for Flowback Models

Flowback from shale oil wells often initially show only a typical single-phase region. Assuming there is no matrix flow, the dimensionless form of the governing equation for 1D single-phase flowback is given by

$$\frac{\partial P_D}{\partial t_{aD}}=\frac{{\partial }^2{P}_D}{\partial x_D^2},$$

(1)

where dimensionless variables are defined as

$$P_D=\frac{P}{P_i},$$

(2)

$$x_D=\frac{x}{L},$$

(3)

$$t_{aD}=\frac{1}{{\mu }_w{L}^2}{{\displaystyle \int }}_0^t\frac{k_f\left(\overline{p}\right)}{w_{f_D}\left(\overline{p}\right){\varphi }_f\left(\overline{p}\right)c_f\left(\overline{p}\right)}dt,$$

(4)

where $P$ is fluid pressure along the fracture, $P_i$ is the initial fracture pressure, $\overline{p}$ is the average pressure in the fracture, $L$ is the length of the fracture segment, $\varphi$ is porosity, ${\mu }_w$ is water viscosity, $k_f$ is fracture permeability, $w_{f_D}$ is dimensionless fracture width which is defined as $\frac{w_f}{w_{f_i}}$, x is the coordinate along the fracture direction.$c_t$ is the total compressibility which is the summation of fracture compressibility $c_f$ and water compressibility $c_w$, which are defined as follows

$$c_f=\frac{1}{V_f}\frac{d{V}_f}{dP},$$

(5)

$$c_w=-\frac{1}{V_w}\frac{d{V}_w}{dP}.$$

(6)

In the case of wells with an extended shut-in period, a small amount of oil may enter the fracture network. The following assumptions should be considered to develop the two-phase oil–water model: (1) capillary pressure is neglected inside the fractures, (2) oil and water are both incompressible, and (3) no matrix flow is considered. The average saturation of the fracture domain does not change due to the assumption of incompressible fluids as the coefficients of $\frac{d{S}_w}{dt}$ terms collapse to zero in the two-phase governing equation. The mass balance equation for the fracture system is the same as Equation (1), except the pseudo-time is defined differently as

$$t_{aD}=\frac{{\lambda }_t}{L^2}{{\displaystyle \int }}_0^t\frac{k_f\left(\overline{p}\right)}{w_{f_D}\left(\overline{p}\right){\varphi }_f\left(\overline{p}\right)c_f\left(\overline{p}\right)}dt,$$

(7)

where ${\lambda }_t$ is the summation of ${\lambda }_o$ and ${\lambda }_w$, which are oil mobility and water mobility defined as $\frac{k_{ro}}{{\mu }_o}$ and $\frac{k_{rw}}{{\mu }_w}$. The straight-line relative permeability model is used for fractures in this study i.e.,

$$k_{rw}=S_w,$$

(8)

$$k_{ro}=1-S_w,$$

(9)

where $S_w$ is water saturation. Detailed derivation of governing equations for a two-phase flow without matrix influx is provided in Appendix B.

During the flowback period, initially, water production dominates (region 1) which is gradually replaced by oil production (region 3) after oil breakthrough the fracture as shown schematically in Figure 1. Since the permeability of shale reservoirs is extremely low, no matrix flow from the reservoir is considered for early-time flowback. As a result, the proposed flowback model will be restricted to region 1. Due to the extremely small pore diameter of the shale formation, the large capillary effect is expected. Therefore, breakthrough pressure can be much lower than initial reservoir pressure, which makes the closed-chamber solution works for an extended period of time.

Even though the solution developed by assuming no matrix influx is accurate in early flowback periods, it cannot predict flowrate accurately for the late flowback period after a significant amount of oil enters the fracture network. To test the applicability of the above-mentioned model without considering matrix flow, we proposed another solution that can consider matrix flow but only limited to wells with bi-wing planar fractures. In the fracture governing equation, the matrix influx is considered as a source term, which is given by

$$\frac{{\partial }^2{P}_p}{\partial x_D^2}=\frac{\partial P_p}{\partial t_{aD}}-q_m,$$

(10)

where

$$P_p={{\displaystyle \int }}_P^{P_i}\frac{{\lambda }_t{k}_f\left(\overline{p}\right)}{{\varphi }_f{c}_f\left(\overline{p}\right)}dp,$$

(11)

$$t_{aD}=\frac{1}{L^2}{{\displaystyle \int }}_0^t\frac{{\lambda }_t}{w_{fD}{\varphi }_f{c}_f\left(\overline{p}\right)}dt,$$

(12)

$$q_m=\frac{L^2{v}_y\left(t\right)}{w_f{P}_i}.$$

(13)

The matrix oil velocity $v_y\left(t\right)$ is calculated from the matrix governing equation

$$\frac{{\partial }^2{P}_D}{\partial y_D^2}=\frac{\partial P_D}{\partial t_D},$$

(14)

where

$$P_D=\frac{P_m}{P_i},$$

(15)

$$y_D=\frac{y}{D},$$

(16)

$$t_D=\frac{k_m}{{\mu }_o{D}^2{\varphi }_m{c}_f}t.$$

(17)

where $P_m$ represents matrix pressure, D represents the length of stimulated reservoir volume (SRV) or half of the fracture spacing for multi-stage hydraulic fractures in horizontal wells. Although $\frac{\partial S_w}{\partial t}$ the term is canceled out, water saturation inside the fracture should be updated from the volume balance equation to update the relative permeability values.

3. Analytical Solution of Governing Equations

The eigenfunction expansion methods are used to solve the diffusivity equation with time-dependent boundary conditions. The analytical solutions of Equation (1) can be classified into two cases based on the boundary conditions. One case is variable pressure boundary conditions at both fractures’ ends and the other case is variable pressure at one end and no flow at the other end. Detailed derivation of the solution is provided in Appendix A. Initial condition of Equation (1) for is given by

$$P_D\left(x_D,0\right)=1.$$

(18)

First, we present the solution for variable pressure boundary conditions at both ends of the fracture without matrix influx, time dependent pressure boundary conditions are prescribed at both ends as

$$P_D\left(0,t_{aD}\right)=P_1\left(t_{aD}\right),$$

(19)

$$P_D\left(1,t_{aD}\right)=P_2\left(t_{aD}\right).$$

(20)

Using an eigenfunction expansion for the solution $p$ as shown in Appendix A, analytical solution of Equation (1) with two Dirichlet boundary conditions are given by

$$P_D\left(x_D,{\mathrm{t}}_{\mathrm{aD}}\right)=\left(P_2\left({\mathrm{t}}_{\mathrm{aD}}\right)-P_1\left({\mathrm{t}}_{\mathrm{aD}}\right)\right)x_D+P_1\left({\mathrm{t}}_{\mathrm{aD}}\right)+{{\displaystyle \sum }}_{n=1}^{\infty }\left({{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(\tau \right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right)sin\left({\lambda }_n{x}_D\right),$$

(21)

where $c_n$ is the Fourier series obtained from the initial conditions, ${\widehat{S}}_n\left(t_{aD}\right)$ is the term that accounts for the time derivative of the two Dirichlet boundary conditions, ${\lambda }_n$ is the eigenvalue:

$$c_n={{\displaystyle \int }}_0^1\left[1-P_2\left(0\right)+P_1\left(0\right)\right]\sin \left({\lambda }_n{x}_D\right)d{x}_D,$$

(22)

$${\widehat{S}}_n\left(t_{aD}\right)=-2{{\displaystyle \int }}_0^1\sin \left({\lambda }_n{x}_D\right)\left(\frac{d{P}_1}{d{t}_{aD}}-\frac{d{P}_2}{d{t}_{aD}}\right)d{x}_D,$$

(23)

$${\lambda }_n=\frac{n\pi }{L}.$$

(24)

In the single-phase solution, water rate can be evaluated by taking derivative of Equation (21) with respect to x as

$${q_w|}_{{\mathrm{x}}_{\mathrm{D}}=0}=\frac{w_f{k}_{f}h{P}_i}{{\mu }_{w}L}\left[P_2\left(t_{aD}\right)-P_1\left(t_{aD}\right)+{{\displaystyle \sum }}_{n=1}^{\infty }\left({{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(\tau \right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right)\right],$$

(25)

$${q_w|}_{{\mathrm{x}}_{\mathrm{D}}=1}=\frac{w_f{k}_{f}h{P}_i}{{\mu }_{w}L}\left[P_2\left(t_{aD}\right)-P_1\left(t_{aD}\right)+{{\displaystyle \sum }}_{n=1}^{\infty }\left({{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(\tau \right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right)\right]\cos \left({\lambda }_n\right).$$

(26)

Similarly, in the two-phase solution, the total flow rate at the inlet and outlet are

$${q_t|}_{{\mathrm{x}}_{\mathrm{D}}=0}=\frac{{\lambda }_t{w}_f{k}_{f}h{P}_i}{L}\left[P_2\left(t_{aD}\right)-P_1\left(t_{aD}\right)+{{\displaystyle \sum }}_{n=1}^{\infty }\left({{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(\tau \right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right)\right],$$

(27)

$${q_w|}_{{\mathrm{x}}_{\mathrm{D}}=1}=\frac{w_f{k}_{f}h{P}_i}{{\mu }_{w}L}\left[P_2\left(t_{aD}\right)-P_1\left(t_{aD}\right)+{{\displaystyle \sum }}_{n=1}^{\infty }\left({{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(\tau \right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right)\right]\cos \left({\lambda }_n\right).$$

(28)

For solution with variable pressure boundary condition at one end and no-flow boundary at the other end without matrix influx, time-dependent pressure boundary conditions are prescribed at both ends as

$$P_D\left(0,t_{aD}\right)=P_1\left(t_{aD}\right),$$

(29)

$$\frac{\partial P_D\left(1,t_{aD}\right)}{\partial x_D}=0.$$

(30)

Using an eigenfunction expansion for the solution $p$ as shown in Appendix A, the analytical solution of Equation (1) with one Dirichlet boundary condition and one Neumann boundary condition is

$$P_D\left(x_D,{\mathrm{t}}_{\mathrm{aD}}\right)=P_1\left(t_{aD}\right)+{\displaystyle \sum }_{n=1}^{\infty }\left[{{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(t_{aD}\right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right]sin\left({\lambda }_n{x}_D\right),$$

(31)

where

$$c_n={{\displaystyle \int }}_0^1\left(1-P_1\left(0\right)\right)\sin \left({\lambda }_n{x}_D\right)d{x}_D,$$

(32)

$${\widehat{S}}_n\left(t_{aD}\right)=-2{{\displaystyle \int }}_0^1\sin \left({\lambda }_n{x}_D\right)\frac{d{P}_1\left(t_{aD}\right)}{d{t}_{aD}}d{x}_D=-\frac{2}{{\lambda }_n}\frac{\partial P_{\mathrm{D}1}}{\partial t_{aD}},$$

(33)

$${\lambda }_n=\frac{\left(2n-1\right)\pi }{2}.$$

(34)

In the single-phase solution, water rate can be evaluated by taking the derivative of pressure, i.e.,

$${q_w\left(t_{aD}\right)|}_{{\mathrm{x}}_{\mathrm{D}}=0}=\frac{w_f{k}_{f}h{P}_i}{{\mu }_{w}L}\left[{{\displaystyle \sum }}_{n=1}^{\infty }\left({{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(\tau \right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right)\right].$$

(35)

In the two-phase solution, total flow rate can be evaluated as

$${q_t\left(t_{aD}\right)|}_{{\mathrm{x}}_{\mathrm{D}}=0}=\frac{{\lambda }_t{w}_f{k}_{f}h{P}_i}{L}\left[{{\displaystyle \sum }}_{n=1}^{\infty }\left({{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(\tau \right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right)\right].$$

(36)

For two-phase flowback with matrix influx, the concept of pseudo-pressure is used to capture the saturation change inside the fracture. Initial condition of Equation (10) for fracture flow is given by

$$P_p\left(x_D,0\right)=0.$$

(37)

The time-dependent pressure boundary conditions are prescribed at both ends as

$$P_p\left(0,t_{aD}\right)=P_{p1}\left(t_{aD}\right),$$

(38)

$$\frac{\partial P_p\left(1,t_{aD}\right)}{\partial x_D}=0.$$

(39)

Using the eigenfunction expansion for homogenized pseudo-pressure, the analytical solution of Equation (10) will be

$$P_p\left(x_D,{\mathrm{t}}_{\mathrm{aD}}\right)=\left(t_{aD}\right)+{\displaystyle \sum }_{n=1}^{\infty }\left[{{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(t_{aD}\right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right]sin\left({\lambda }_n{x}_D\right),$$

(40)

where

$$c_n={{\displaystyle \int }}_0^1\left(1-P_1\left(0\right)\right)\sin \left({\lambda }_n{x}_D\right)d{x}_D,$$

(41)

$${\widehat{S}}_n\left(t_{aD}\right)=-2{{\displaystyle \int }}_0^1\sin \left({\lambda }_n{x}_D\right)\left(\frac{d{P}_{p1}}{d{t}_{aD}}+q_m\right)d{x}_D=-\frac{2}{{\lambda }_n}\left(\frac{d{P}_{p1}}{d{t}_{aD}}+q_m\right),$$

(42)

$${\lambda }_n=\frac{\left(2n-1\right)\pi }{2}.$$

(43)

Taking the derivative of the pseudo-pressure, the total flow rate at the wellbore can be evaluated as

$${q_t\left(t_{aD}\right)|}_{{\mathrm{x}}_{\mathrm{D}}}=\frac{w_f{k}_{f}h{P}_i}{L}{\displaystyle \sum }_{n=1}^{\infty }\left[{{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(t_{aD}\right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right]cos\left({\lambda }_n{x}_D\right).$$

(44)

Water rate and oil rate can be evaluated as,

$${q_w\left(t_{aD}\right)|}_{x_D}=\frac{{\lambda }_w}{{\mathsf{\lambda }}_t}{q}_t\left(t_{aD}\right),$$

(45)

$${q_o\left(t_{aD}\right)|}_{x_D}=\frac{{\lambda }_o}{{\mathsf{\lambda }}_t}{q}_t\left(t_{aD}\right).$$

(46)

Since the source term in Equation (10) is obtained by solving the matrix flow equation, matrix oil velocity also needs to be updated. The computational cost of the solution is O(M), where M is the number of total segments. To balance the accuracy and computational costs, we divided each half of the fracture into four segments Figure 2 shows hydraulic fractures in a given stimulated reservoir volume (SRV), which is assumed to have a rectangular shape. Matrix oil flux is calculated separately for the four fracture segments. Initial and boundary conditions of the matrix oil flow equation is given by

$$P_D\left(0,t_D\right)={\overline{P}}_{fD,k}k=1,2,3,4$$

(47)

$$\frac{\partial P_D\left(1,t_D\right)}{\partial y_D}=0,$$

(48)

$$P_D\left(y_D,0\right)=1,$$

(49)

where ${\overline{P}}_{fD,k}$ represents the average pressure at segment k. The analytical solution of Equation (14) is given by

$$P_{D, k}\left(x_D,{\mathrm{t}}_{\mathrm{aD}}\right)={\overline{P}}_{fD,k}+{{\displaystyle \sum }}_{n=1}^{\infty }\left[{{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(t_{aD}\right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right]sin\left({\lambda }_n{x}_D\right),$$

(50)

$${v_y\left(t_{aD}\right)|}_{{\mathrm{x}}_{\mathrm{D}}=0}=\frac{k_f{P}_i}{D}{\displaystyle \sum }_{n=1}^{\infty }\left[{{\displaystyle \int }}_0^{t_{aD}}{e}^{-{\lambda }_n^2\left(t_{aD}-\tau \right)}{\widehat{S}}_n\left(t_{aD}\right)d\tau +e^{-{\lambda }_n^2{t}_{aD}}{c}_n\right].$$

(51)

Average water saturation in the fracture segment k is determined using volume balance as

$$S_{w,k}=S_{wi}+\frac{\left(q_{w,in,k}-q_{w,out,k}\right)dt}{V_{p,k}},$$

(52)

where $S_{wi}$ is the initial water saturation, dt is time step and $V_{p,k}$ is the volume of the fracture segment at the current time. $q_{w,in,k}$ and $q_{w,out,k}$ is the inflow and outflow water rates, respectively. The continuous water saturation profile inside the fracture is approximated by a second-order polynomial interpolation between nodal saturations $S_{w,k}.$ The two-phase flowback solution with matrix oil influx can be described as following:

• For the first time step, the average saturation in each segment is equal to initial water saturation ($S_{w,k}=S_{wi}$). Solve Equation (10) by Equations (40) and (44) assuming no oil influx at the initial time ($v_o=0$) and calculate average fracture pressure ${\overline{P}}_{f,k}$ for each fracture segment.
• Based on ${\overline{P}}_{f,k}$, calculate average oil influx $v_{y,k}$ using Equation (51).
• For the second time step, update oil rate $v_{y,k}$ calculated from step 2 in Equation (10) and solve it again. Average water saturation for each segment is calculated from the material balance Equation (52). Update pressure and repeat steps 2–3 till reaching the end of the simulation time.

We also show the workflow of the solution in Figure 3.

The analytical solutions presented in this section is only applicable to a single fracture segment. however, fractures may exist in the form of complex networks. In the following sections, we will show how we can discretize the complex networks and solve the systems of equations simultaneously for the whole system.

4. Fractal Fracture Network and Fracture Discretization

Two different fracture network geometries are considered in this study: tree-shape fracture networks and orthogonal fracture networks. Fracture re-initiation at the intersection with natural fractures is a possible branching mechanism during hydraulic fracturing [11]. Tree-shape fracture networks are commonly found as shown in Figure 4b from simulation results. Researchers used a neural network model to process microseismic data to estimate fracture geometry and also believe that, the fracture network of multiple fractured horizontal wells in highly brittle shale formation is usually tree-shaped [12]. An orthogonal fracture network is another common situation in complex fracture networks. In naturally fractured formations, if the joint sets are orthogonal to each other, like the situation in Marcellus Shale [11], induced hydraulic fractures may likely form an orthogonal mosaic network by activating natural fractures as shown in Figure 4b.

Fracture discretization is the first step to simplify calculations. Fracture networks can be divided into multiple segments and nodes as illustrated in Figure 4a and Figure 5a for tree-shape and orthogonal fractures, respectively. The tree-shape network is divided into 18 segments and 8 nodes. The orthogonal fracture network is divided into 24 segments and eight nodes. It is worth noting that in the case of symmetric tree-shape or symmetric orthogonal fracture networks, one may take advantage of the symmetry and reduce the number of nodes and segments to reduce the computational costs.

To generalize fracture networks, tree-shape fracture networks and orthogonal fracture networks can be synthesized using fractal theory as shown in Figure 6. In a tree-shape fracture network (Figure 4a), each fracture segment can be modeled with segment length $L_i$, width $w_{f_i},$ and compressibility $c_{f_i}$, where index i denotes the fracture-segment generation level. The oldest or initial segment is obviously near the well (i = 0).

The structure and properties of the network can be characterized by fractional changes after each generation based on fractal theory i.e.,

$$L_i=r_L{L}_0,$$

(53)

$$w_{f_i}=r_{w_f}{w}_{f0}.$$

(54)

In orthogonal fracture networks, the width of the fracture segments decreases with the distance from the wellbore, which can be easily verified by stress analysis [13]. The structure and properties of the orthogonal network can also be characterized by fractional reduction after each intersection away from the wellbore as shown in Equations (53) and (54). The compressibility of the activated natural fractures which are labeled with level number larger than 0 in Figure 6 because width reduction also indicates less layers of proppants. The permeability of the activated natural fractures is assumed to be the same.

For the case of long planar fractures, non-uniform fracture width and non-uniform fracture closure can be simulated by discretizing the fracture into several segments as shown in Figure 7. The width of fractures decreases away from the wellbore. For each segment, the width of the segment can be set with an average value based on experience or hydraulic fracturing simulation results [13,14,15].

To model fracture closure more accurately, understanding of fracture properties is of great importance. Current models for flowback analysis assume that fracture properties, such as fracture permeability and fracture compressibility, do not change over time. In this paper, we will estimate the evolution of these properties with time-based on the mechanical interaction between the formation rock and proppants. For propped fractures, compressibility is contributed by two components: fracture aperture changes and fracture porosity change [16], i.e.,

$$c_f=c_{\varphi }+c_{width}.$$

(55)

To incorporate fracture aperture changes during flowback, we adopt the analytical solution for normal compliance of the propped fracture [17]. Assuming there is no slippage between the particles and deformation between the proppant particles is elastic, the compressibility of the fracture due to its width change $c_{width}$ can be obtained by

$$c_{width}=\frac{1}{w_f}\left[\frac{4}{\sqrt[3]{3}}{\left(\frac{1-{\nu }_f^2}{E_f}+\frac{1-{\nu }_p^2}{E_p}\right)}^{\frac{2}{3}}+\left(n-1\right){\left(\frac{1-{\nu }_p^2}{E_p}\right)}^{\frac{2}{3}}\right]R{\sigma }_n^{-\frac{1}{3}},$$

(56)

where n is the number of proppant layers, ${\sigma }_n$ is the effective normal stress on proppants, $E_f$ and $E_p$ are the Young’s Moduli of the formation and proppants, respectively. ${\nu }_f$ and ${\nu }_p$ are Poisson’s ratios of the formation rock and proppant particles, respectively. We adopt a model to estimate the porosity of the proppants pack [18], that can be simplified as

$$\varphi ={\varphi }_i{\left[1-C_o{\left(\frac{{\sigma }_n}{E_p}\right)}^{\frac{2}{3}}\right]}^3\mathrm{for}{\sigma }_n<<E_p,$$

(57)

where ${\varphi }_i$ is the initial porosity of the proppant pack, $C_o$ is a constant depending on the packing structure. Hence, porosity compressibility is estimated from the model by assuming that fracture porosity is equal to the proppant pack as,

$$c_{\varphi }=\frac{2C_o}{{\sigma }_n}\frac{\frac{{\sigma }_n}{E_p}}{1-C_o{\left(\frac{{\sigma }_n}{E_p}\right)}^{\frac{2}{3}}}.$$

(58)

Permeability can also be modeled as in [18]

$$k=k_i{\left[1-C_o{\left(\frac{{\sigma }_n}{P_o}\right)}^{\frac{2}{3}}\right]}^4.$$

(59)

Thus, total fracture compressibility can be expressed as a function of effective normal stress, proppant properties and mechanical properties of the formation,

$$c_f=\frac{2C_o}{{\sigma }_n}\frac{\frac{{\sigma }_n}{E_p}}{1-C_o{\left(\frac{{\sigma }_n}{E_p}\right)}^{\frac{2}{3}}}+\frac{1}{w_f}\left[\frac{4}{\sqrt[3]{3}}{\left(\frac{1-{\nu }_f^2}{E_f}+\frac{1-{\nu }_p^2}{E_p}\right)}^{\frac{2}{3}}+\left(n-1\right){\left(\frac{1-{\nu }_p^2}{E_p}\right)}^{\frac{2}{3}}\right]R{\sigma }_n^{-\frac{1}{3}}.$$

(60)

In this study, we compare the flowback response for fracture propped with two kinds of proppants, sand, and ceramics whose properties are shown in

Table 1. Properties of two types of proppants used in this study.

Properties	Ceramics	Sand
${\nu }_p$	0.2	0.2
$E_p$(Gpa)	70	20
R (mm)	0.3	0.3
Co	2	2

. In comparison to sand, ceramics have a higher strength and is less compressible. Young’s modulus and Poisson’s ratio of the formation used in this study is assumed to be 25 Gpa and 0.27, respectively. The size of the proppants used in this study is close to 40–70 mesh, we consider an average radius of the proppant is 0.3 mm in the analytical model on fracture compressibility. However, our proposed model is not limited to a specific mesh size and the same procedure can be applied for different proppant mesh sizes.

It is worth noticing that in the case with partial proppants coverage, a considerably high fracture compressibility can be used to model the rapid closure of the unpropped part of the fracture.

Upon discretization of the complex fracture networks, the analytical solution for each discretized fracture segment will be solved simultaneously. The next section is dedicated to introducing the numerical technique we utilized for the semi-analytical approach.

5. Workflow of Semi-Analytical Solution for Fracture Networks

Since flowrate is calculated by using the pressure gradient, the pressure boundary condition at each node is solved using the mass balance equation

$${\left(q_{in}\right)}_i={\left(q_{out}\right)}_{i}i=1\dots M.$$

(61)

In Equation (35), for the single-phase region, water rate is used while total saturation is used for the two-phase region. The corresponding residual function is defined as

$$R_i={\left(q_{in}\right)}_i-{\left(q_{out}\right)}_{i}i=1\dots M.$$

(62)

The Newton-Raphson method is utilized to solve the above equations. It should be noted that using this method with a tolerance value of 1 × 10⁻¹⁶, convergence is achieved by only two to four iterations at each time step. The system of equations is expressed as

$$JdP=-R,$$

(63)

where

$$J=\left[\begin{array}{cccc}\frac{\partial R_1}{\partial P_1}& \frac{\partial R_1}{\partial P_2}& \dots & \frac{\partial R_M}{\partial P_M}\\ \frac{\partial R_2}{\partial P_1}& \frac{\partial R_2}{\partial P_2}& \dots & \frac{\partial R_2}{\partial P_2}\\ \dots & \dots & \dots & \dots \\ \frac{\partial R_M}{\partial P_1}& \frac{\partial R_M}{\partial P_2}& \dots & \frac{\partial R_M}{\partial P_M}\end{array}\right],$$

(64)

$$R=\left[\begin{array}{c}{R}_1\\ R_2\\ \dots \\ R_M\end{array}\right].$$

(65)

Pressure is updated after an arbitrary iteration k as

$$P_{k+1}=P_k+d{P}_k.$$

(66)

For the two-phase water-oil solution, pseudo-pressure is updated as

$$P_p{}_{k+1}=P_p{}_k+d{P}_p{}_k.$$

(67)

A schematic diagram of the semi-analytical workflow is shown in Figure 8.

6. Results and Discussions

To test the accuracy of the proposed semi-analytical model for flowback analysis of the wells completed in shale oil reservoirs, several examples are generated and compared with results from full-fledged numerical simulations with COMSOL Multiphysics and CMG for the two different topologies for fracture networks used in this study. In COMSOL Multiphysics, the diffusivity equation is solved by the built-in two-phase Darcy’s flow module using the finite element methods, which allows the generation of fracture networks with non-orthogonal topology. About 3000 quadratic elements is used to solve the two-dimensional, two-phase oil–water flow in both fractal tree-shape and orthogonal fracture networks without matrix influx. Only half of the fracture is modeled due to the symmetry of the problem. A pressure boundary condition is applied at the mouth of the fracture. In CMG, only a quarter of the reservoir is modeled due to the symmetry of the problem with finite difference methods. About 10,000 elements including fracture elements and matrix elements are used to model the two-phase flow with matrix influx. A layer of the elements represents one-half of the fracture, which has higher permeability in comparison to the other elements. The element at the end of the fracture elements is representing the well. The H-adaptive mesh is adopted i.e., elements are coarsening as moving from the fracture to the distant reservoir.

A summary of the reservoir and fracture network properties used in our examples is provided in

Table 2. Input parameters for validation of the semi-analytical solutions.

Properties	Tree-Shape	Orthogonal
Initial fracture pressure (Pa)	3 × 107	3 × 107
Initial Fracture water saturation (fraction)	0.95	0.95
Fracture porosity (fraction)	0.3	0.3
Fracture height (m)	30	30
Fracture width of first level (m)	0.01	0.01
Fracture length of first level (m)	100	50
Width reduction ratio	0.6	0.75
Length reduction ratio	0.5	1
Proppant type	sand	sand
Fracture maximum permeability (m2)	1 × 10−12	1 × 10−12
Water viscosity (Pa-s)	0.001	0.001
Simulation Time (days)	3	3
Number of fracture stages	20	20

. The fracture networks shown in Figure 4 are used for validation purposes. For the orthogonal fracture network, starting from the wellbore, fractures’ width and permeability reduce 25% after each intersection away from the wellbore. For the tree-shape fracture network, fracture segments’ width reduces 50% after each generation as indicated by earlier geomechanical simulations [13]. In this study, we presumed a smooth pressure drawdown characteristic at the wellbore during early flowback using the below form function,

$$P_{wf}\left(t\right)=\frac{1}{1.5}({\mathrm{e}}^{-\frac{t}{a}}+0.5)P_i,$$

(68)

where a is the drawdown coefficient with the value of 100,000 for our validation example, t is time (in seconds).

A good match for total production between the semi-analytical model and commercial numerical simulations has been obtained to verify the proposed analytical solution as shown in Figure 9. The computational cost of the semi-analytical model in MATLAB (MathWorks, Natick, MA, USA) and the numerical model in COMSOL Multiphysics is compared in Figure 10. In order to examine the influence of the fracture network topology on the flowback behavior, we compare flowback rate from a tree-shape and an orthogonal network fractures sets with the same total length and total fracture volume as shown in Table 2. The reader may notice that the flowrate in the orthogonal fracture network reaches a higher peak in flowrate in comparison to tree-shape fracture but declines faster. The probable reason for this phenomenon is higher conductivity in the orthogonal fracture set in comparison to the tree-shape network with the same volume and length. Different flowback behaviors in the two fracture networks indicate that flowback analysis based on the material balance could be misleading since it cannot capture the effect of the network topology.

For the purpose of doing a sensitivity analysis, we considered the cases described earlier as the base case, and only change a key parameter each time with respect to the base case.

The total fracture volume of the fracture network can be calculated using the length of fracture segments and their corresponding widths. Since it is not realistic to assign a single width to all fractures and on the other hand, it would be computationally very expensive to consider a variable width along all each fracture segment, we use a simplifying hypothesis consistent with geomechanical analyses. We are presuming a value for the width of the fracture segment connected to the well and then assuming a fixed ratio for the fracture width reduction at each intersection away from the wellbore. Therefore, we can investigate the effect of the fracture width distribution by changing the maximum width at the wellbore and the ratio of width changes. In general, increasing fracture width or its length increases the volume inside the fracture network. The influence of the fracture volume on flowback looks similar in both fracture networks.

6.1. Flowback Responses with Different Fracture Volume

Figure 11 compares the flowback rate for different maximum fracture width (different $w_{fi}$) in tree-shape networks and orthogonal networks. We can see that in both cases, the flowback rate is nearly proportional to the maximum fracture width (otherwise known as average fracture width). It is also notable that the peak time for the flowrate does not shift with changes in the maximum fracture width as the topology of the fracture network remains the same. In shale gas wells due to nonlaminar effects, this observation may not be necessarily true. Figure 12 compares the flowback rate for different ratios of the fracture width reduction (different $r_{w_f}$). The result shows that the flowback rate looks similar initially during the transient flow period. Though after a while, the flowrates slightly increase with higher $r_{w_f}$. We also noticed that changing $r_{w_f}$ shifts the peak time for the flowback rate. The occurrence of peak flowrate delays slightly with higher value of $r_{w_f}$. A similar phenomenon may also happen with the increase of the total fracture length. Since several parameters affect the flowback behaviors together, the inverse analysis using the proposed model can be non-unique. Hence, constraining one of these values would be helpful to reach a unique fracture network geometry during inverse analyses. For instance, knowledge of micro-seismic data and injection volume can be helpful for estimating the total fracture volume. Fracture width reduction can also be estimated from a coupled analysis of the hydraulic fracturing treatment. Proppant concentration, fracking fluid viscosity, viscosity degradation of the fracking fluids, natural fractures direction with respect to minimum horizontal stress and injection rates are among different factors that may change the fracture width reduction ratio.

6.2. Flowback Responses with Different Proppants

Since fracture closure due to in-situ tectonic stresses could be one of the main mechanisms for the flowback of the fracking fluid from induced fracture networks, we investigate the possible effect of the fracture compressibility on the flowback behavior. Figure 13 compares the flowback response in different types of proppants i.e., sand and ceramics. We found that the impact of the proppant type on the flowback is almost the same for the two different fracture networks. Ceramic proppants are less compliant than the regular sand proppants. In both tree-shape and orthogonal networks, flowback rate and recovery can be highly impacted by the proppant type and its corresponding compressibility. As expected, since sand has a lower Young’s modulus and a higher compressibility, the water recovery rate increases in comparison to the case with the ceramic proppants. We also noticed that changing proppant type shifts the peak time for the flowback rate. The occurrence of peak flowrate delays slightly with less compressible material. The results suggest that using more compressible proppants may benefit the fracturing fluid recovery, however, it may not provide a high-enough hydraulic conductivity for future production due to the permeability loss driven by compaction.

6.3. Flowback Responses with Different Drawdown Strategies

Botthomhole pressure is used as the inner boundary condition. Three different operational schemes are assumed in this study: 1) normal-pressure drawdown, 2) slow pressure drawdown, and 3. Fast pressure drawdown. The drawdown coefficient a in Equation (68) is ranging from 20,0000, 10,0000, to 50,000 to represent slow, normal and fast drawdown cases, respectively. Simulation results for the two fracture networks are shown in Figure 14. The influence of botthomhole pressure drawdown on flowback is similar in the two fracture networks. The initial flowback rate is high, but it also declines faster when the bottomhole pressure drawdown is implemented more aggressively. A faster pressure drawdown results in a higher cleanup rate. However, even though an initial aggressive flowback rate may accelerate the onset of hydrocarbon production, excessively high initial flowback rates may cause proppant flowback and damage to the fracture conductivity as a result of large drag forces driven by the high pressure-gradient forming along the fracture. It can be seen in Figure 14 that using the normal drawdown scheme, we can still achieve a good fracking fluid recovery in the late flowback period in comparison to the aggressive drawdown case, however, a concrete conclusion requires coupling with oil flow from the matrix.

6.4. Flowback Responses with Isotropic Fracture and Anisotropic Fracture Sets

In general, it is expected that the fracture’s width decreases with distance from the wellbore. In orthogonal fracture network, while the difference between maximum horizontal stress ($S_{Hmax}$) and minimum horizontal stress ($S_{hmin}$) is not significant, fractures’ width is expected to remain almost the same in different directions. However, when $S_{Hmax}$ is considerably larger than $S_{H,min}$, it is reasonable to expect that opening for fractures perpendicular to the maximum horizontal stress to be smaller than the one for fractures perpendicular to the minimum horizontal stress. Thus, it makes sense to examine the effect of the difference between maximum and minimum horizontal stresses on the flowback behavior in orthogonal fracture networks (like Figure 5b). The total fracture volume is the same as the volume in the base case, however, the fractures’ width is assumed to be different depending on the fracture direction. We assumed that the fracture width in one direction to be 10% higher than the isotropic case while the fracture width in the other direction is 10% lower than the isotropic case. We observe that tectonic stress anisotropy would decrease the flowrate slightly as shown in Figure 15, which is due to the fact that increasing anisotropy can reduce the overall connectivity. However, a concrete conclusion required coupled geomechanical analyses. The difference between the flowrate in the two cases gradually increases till reaching the peak rates and then starts to decrease. The peak of flowback rate decreased about 10% and average flowback rate decreased about 5% in comparison to the case in which the difference between $S_{Hmax}$ and $S_{hmin}$ is not significant. This indicates that the water recovery is slightly higher in the formation which has a smaller difference between $S_{Hmax}$ and $S_{Hmin.}$

6.5. Flowback Responses with Constant Compressibility and Variable Compressibility Models

In current works, flowback analysis is based on the assumption that fracture compressibility remains constant during flowback and production, which indeed according to Equation (60) should be a function of the fracture width and the effective proppant stress. Here, we examine the error that may arise from this assumption. Figure 16 shows that in both fracture networks, the flowback by assuming constant fracture compressibility leads to underestimating flowrate peak and overestimating later production. The result suggests that the flowback analysis using constant compressibility can overestimate the fracture compressibility at the end of the flowback production. By using the proposed semi-analytical model, we can model flowback rate more accurately by dropping the invalid assumption of constant compressibility that can benefit subsequent inverse analyses.

In addition, all the above case studies indicate that with the same fracture volume, an orthogonal fracture network will provide higher flowback production at an early time while the depletion is faster as well. This is because fracture connectivity is higher in orthogonal fracture network.

6.6. Applicability/Limitation of the Proposed Closed-Chamber Solution

The proposed semi-analytical solution presented in this paper only considers fluid flow inside the fracture network. Basically, the matrix flow is neglected in this solution as it takes a while for the reservoir fluid to enter the fracture. In order to gain an insight into the limitation of the proposed closed-chamber solution, we derived a two-phase flowback solution with a matrix influx for bi-wing fractures. The proposed solution with matrix inflow has been validated by CMG as shown in Figure 17 with the input data listed in

Table 3. Input parameters of validation case for a two-phase oil–water solution with the matrix influx.

Properties	Value
Initial fracture pressure (Pa)	3 × 107
Initial Fracture water saturation (fraction)	1
Fracture porosity (fraction)	0.6
Fracture height (m)	30
Fracture width (m)	0.02
Fracture length (m)	100
Fracture permeability (${\mathrm{m}}^2$)	1 × 10−12
Matrix permeability $({\mathrm{m}}^2)$	1 × 10−18
Matrix porosity(fraction)	0.05
Total matrix compressibility (1/Pa)	1 × 10−9
Drawdown coefficient	1 × 106
Proppant Type	Sand
Oil viscosity (Pa·s)	0.001
Water viscosity (Pa·s)	0.001
Number of fracture stages	20

. To test the accuracy of the proposed model in the extreme case when oil flux is significant, the oil saturation of the matrix is set to be one. Drawdown coefficient is 50,000 to represent an aggressive choke strategy. To capture the worst-case scenario to limit our flowback solution, the reservoir breakthrough pressure is set to be equal to the initial fracture pressure in the model with a formation influx. However, the true reservoir pressure is much lower than the fracture pressure at the beginning of flowback, additionally, there would be a huge capillary effect which blocks the formation fluid from entering the fracture. Therefore, lower breakthrough pressure should be used in practical analyses. The result indicates that the two-phase flowback solution with matrix oil influx is accurate while its computational cost is much lower than numerical simulation. However, it is worth noting that the model has more error when the oil velocity in matrix is high. This error is rooted in the definition of our two-phase pseudo-pressure, in which the saturation gradient inside the fracture cannot be captured accurately as a continuous function in the case of the large pressure gradient. However, this situation rarely happens during flowback in shale oil reservoirs due to the extremely low permeability of the matrix. This solution is only developed for simple fracture geometries and despite the first solution, this solution cannot be applied to complex fracture networks. But this solution can still serve as a tool to gain insight on the time limitation of the semi-analytical model, i.e., determine when matrix inflow starts to dominate the flowback flow.

We did three case studies to examine the allowable time interval to use the proposed close-chamber solution. All other inputs are similar to the validation case shown in Table 3. The criteria that we used to define transition time from the fracture linear flow to the matrix linear flow is the time that the flowback rate from the proposed solution deviates more than 20% from the flowback calculated from the solution with matrix oil influx, the model is considered inapplicable hereafter. Figure 18 shows the comparison of the flowrates with and without matrix flow. The results show that the time when the proposed solution without oil matrix influx becomes invalid is 25, 27, and 35 hours for ases 1 to 3, respectively. Case 1 has the highest oil production among all the cases which makes the applicable time relatively short in comparison to the other two cases. Breakthrough pressure is 2.3 × 10⁷ Pa. In case 2, we consider the formation has a smaller pore diameter (i.e., less porosity, less permeability, and higher capillary pressure). Porosity and permeability used is 0.03 and 1 × 10⁻¹⁹ ${\mathrm{m}}^2$, respectively. The breakthrough pressure is decreased to 1.9 × 10⁷ Pa. The significant capillary effect postpones the oil breakthrough time. The result indicates the model can be applied for an extended time in shale reservoirs with extremely small pore sizes. In case 3, we consider the formation of oil is more viscous, so oil viscosity is increased to 0.002 Pa·s. When the reservoir fluid is more viscous, matrix oil influx will be slower and the flowback model can be used for a longer time. From the applicability analysis, we show that the proposed semi-analytical solution is working for at least one day or longer times depending to reservoir properties. The allowable time window for this analysis is long enough to make data collection feasible. However, we suggest increasing data collection frequency for more accurate inverse analysis.

6.7. Field Example

To test the practical applicability of the proposed model, we use the flowback data from a multi-fractured horizontal well completed in a tight oil reservoir in Western Canada [5]. First, using the proposed solution for planar fractures with matrix influx, we showed that a similar production match can be achieved as shown in Figure 19. As we did not include resolved gas in oil in our model, the history matching is based on water production data only.

As indicated in their study, the authors pointed out that a larger than typical value of fracture width is used in the history matching to account tor activated secondary fractures during stimulation.

7. Conclusions

We proposed a semi-analytical solution for two-phase oil-water flowback within the complex fracture network in shale oil wells. The model can capture non-uniform fracture openings within the fracture network and incorporates situations that some oil exists inside the fractures before initiating flow in the well due to imbibition. The solution collapses to the single-phase solution when initial water saturation is one. The results of the semi-analytical solution show a good match with numerical simulations. However, due to the assumption of no matrix influx, the model is only applicable for early-time flowback. The error may arise when a significant amount of oil enters into fractures. As a result, the operators may only apply early data for the purpose of inverse analysis to avoid significant error. We have also proposed a two-phase flowback solution with a matrix oil influx to test the limitation of the closed-chamber solution. Analysis shows that the proposed semi-analytical closed-chamber model is applicable for at least one day for a well with quick oil breakthrough which ensures enough time for the data collection. For wells with bi-wing fractures for formations without a significant amount of natural fractures, the later model with matrix influx can be used to access reservoir properties (e.g., formation permeability, porosity) as well. The semi-analytical model with complex fracture networks can also be further improved by including matrix flow in future works. The proposed models provide a quick tool for inverse analysis which will be addressed in future studies.