The Inversion Method of Shale Gas Effective Fracture Network Volume Based on Flow Back Data—A Case Study of Southern Sichuan Basin Shale

Abstract

Fracture network fracturing is pivotal for achieving the economical and efficient development of shale gas, with the connectivity among fracture networks playing a crucial role in reservoir stimulation effectiveness. However, flow back data that reflect fracture network connectivity information are often ignored, resulting in an inaccurate prediction of the effective fracture network volume (EFNV). The accurate calculation of the EFNV has become a key and difficult issue in the field of shale fracturing. For this reason, the accurate shale gas effective fracture network volume inversion method needs to be improved. Based on the flow back characteristics of fracturing fluids, a tree-shaped fractal fracture flow back mathematical model for inversion of EFNV was established and combined with fractal theory. A genetic algorithm workflow suitable for EFNV inversion of shale gas was constructed based on the flow back data after fracturing, and the fracture wells in southern Sichuan were used as an example to carry out the EFNV inversion. The reliability of the inversion model was verified by testing production, cumulative gas production, and microseismic results. The field application showed that the inversion method proposed in this paper can obtain tree-shaped fractal fracture network structure parameters, fracture system original pressure, matrix gas breakthrough pressure, fracture compressibility coefficient, reverse imbibition index, equivalent main fracture half length, and effective initial fracture volume (EIFV). The calculated results of the model belong to the same order of magnitude as those of the HD model and Alkouh model, and the model has stronger applicability. This research has important theoretical guiding significance and field application value for improving the accuracy of the EFNV calculation.

1. Introduction

The economical and effective development of shale gas plays an important role in the world’s energy structure. The shale gas revolution in the United States has achieved great success and thus the world energy pattern has been reshaped. In 2010, the first shale gas fracturing well was implemented in China. After 10 years of development, shale gas production has achieved a new leap [1]. The practice of shale gas fracturing has found that the information on the characteristics of complex fracture network hidden in the flow back data is very important. In the past, fracturing fluid flow back data were mostly ignored. In recent years, the shale gas effective fracture network volume and fracture characteristics information included in the fracturing fluid flow back data to evaluate hydraulic fracturing effects have attracted attention [2]. Due to the special two-phase seepage characteristics of shale gas, the fracturing evaluation method based on flow back is still in its infancy, and there are very few studies on the effective fracture network volume prediction of shale.

The fracture network that develops post-fracturing serves as the primary pathway for fluid flow back and gas production. It is of great significance to the evaluation of fracturing effect and productivity prediction. Effective fracture network volume information can be obtained by accurately interpreting flow back data [3]. It is of great significance to the evaluation of fracturing effect and productivity prediction. At present, microseismic monitoring is often used to analyze the effect of reservoir stimulation. The volume of reservoir stimulation determined by microseismic analysis rarely provides the actual recoverable reservoir volume of gas shale reservoirs. Reservoir stimulation volumes are overestimated or underestimated due to the occurrence of non-contributing microseismic events and the natural fractures [4]. A relatively new concept used to quantify the actual reservoir stimulated volume is the effective reservoir stimulated volume (ESRV) [5,6]. The ESRV is critical for estimating ultimate recovery, stimulated volume, optimal fracture length, and spacing. Some scholars have carefully identified high-frequency flow back data (flow rate and pressure) as a possible method to extract the properties of hydraulic fractures. A new mathematical model was used to approximate the effective volume of fractures and elucidate the mechanisms of early gas production [7] compared to previous flow back models [8,9,10,11]. The water and natural gas production data were analyzed in the Horn River shale, and the estimation of the effective fracture volume was based on the assumption of a two-phase tank model for the fracture system. In addition, due to the high uncertainty introduced by complex flow back models in the evaluation of fracture pore volume, a two-phase tank model was built for reducing parameter uncertainty in order to estimate the fracture pore volume accurately [12]. Comparing the estimated effective fracture pore volume with fracture design parameters (soaking time and proppant concentration) revealed that fracture closure drives single-phase flow back during further analysis of the effective fracture volume. Therefore, the effective fracture pore volume largely depends on the compressibility of the fracture [13]. The evaluation of fracture compressibility is essential for calculating the effective fracture volume, as well as evaluating the fracture volume change. It was found that fracture compressibility depends on how fracture porosity and pore size vary with effective stress. The fracture compressibility is mainly affected by porosity changes in propped fractures [14]. The harmonic decline (HD) model was used to estimate the water phase flow back and the initial effective fracture volume to study the change in the effective fracture volume. According to the findings, pressure depletion and fracture closure during early-time water flow back can lead to a loss of up to 30% of the effective fracture volume. When gas from the matrix enters the fracture network and offers adequate pressure support during late-time flow back, the rate of fracture volume loss reduces [15]. A gas–water two-phase flow model was established for a fractured horizontal well considering the influence of water saturation on gas–water permeability. Sensitive studies were performed to investigate the impact of various parameters on the productivity of a horizontal well [16,17]. The flow back period can be used to quantitatively identify the characteristics of the formed fracture network and calculate the effective fracture network volume; however, such data analysis is not commonly reported in the literature. Therefore, it is necessary to develop an appropriate model to match the data and acquire valuable information on the volume of the fracture network. In this paper, a tree-shaped fractal fracture network two-phase flow back model was first established based on the flow back data, and then the effective fracture network volume was inverted. And the reliability of the model was verified by cumulative production and microseismic data. The shale gas wells were taken as an example in the southern Sichuan Basin, to obtain the effective fracture network volume, which is useful for evaluating post-fracture effects. The flow back models are summarized in Table 1.

2. Model Development

2.1. The Model of Two-Phase Flow Back

In order to establish the flow back model, the assumptions were as follows:

• (1). The effective crack network of 1/2 single cluster is equivalent to the tree-shaped fractal fracture network shown in Figure 1.

• (2). Consider the reverse imbibition effect and the redistribution of original free gas in activated natural fractures.

• (3). Ignore the capillary pressure in the effective fracture network system and ignore the influence of gravity.

• (4). The matrix is only considered as the gas source. Material exchange occurs between the matrix system and the effective fracture network system through the cross-flow equation.

• (5). The effective fracture network system is an elastic porous medium, and assuming that its compressibility is much greater than that of the shale matrix, the tree-shaped fractal fracture permeability and effective fracture volume are pressure-dependent variables.

As shown in Figure 1, the flow rate of rectangular fractures at level k of 1/2 single-cluster tree-shaped fractal fracture network is described by Equation (1) according to the Hagen–Poiseuille equation:

(1)$Q_{\mathrm{k}}=\frac{W_{\mathrm{fk}}^3{h}_{\mathrm{fk}}}{12\mu }\frac{\Delta P_{\mathrm{k}}}{l_{\mathrm{k}}}$

where ΔP_k is the pressure difference of the k-th level fracture, Pa; μ is fluid viscosity, Pa s; l_k, W_fk, h_fk are the length, width, and height of the k-th level branch fracture, m.

The length, width, and height of cracks at level k are obtained by Equation (2).

(2)$\begin{array}{l}{l}_{\mathrm{k}}=l_0{R}_{\mathrm{L}}^k\\ W_{\mathrm{fk}}=W_{\mathrm{f}0}{R}_{\mathrm{W}}^k\\ h_{\mathrm{fk}}=h_{\mathrm{f}0}{R}_{\mathrm{h}}^k\end{array}$

where l₀, W_f0, and h_f0 are the length, width, and height of the initial tree-shaped fractal fractures, m; R_L, R_W, and R_h are ratio of length, width, and height of fractures, dimensionless.

From Equation (1), it can be known that the viscous resistance of fluid flow in a single fracture is

(3)$R_{\mathrm{k}}=\frac{12\mu l_{\mathrm{k}}}{\beta W_{\mathrm{fk}}^3{h}_{\mathrm{fk}}}=\frac{12\mu R_{\mathrm{L}}^k{l}_0}{\beta R_{\mathrm{W}}^{3k}{R}_{\mathrm{h}}^k{W}_{\mathrm{f}0}^3{h}_{\mathrm{f}0}}$

where β is correction coefficient, which can be corrected according to the conductivity experiment in practical applications.

According to the principle of fluid pressure drop in parallel and in series [18], the parallel model was used to convert each branch into a single channel, and then the series model was applied to convert all single channels into a channel. In this way, the total viscous resistance of the network can be calculated, and the total flow resistance of the tree-shaped fractal fracture network can be expressed as Equation (4):

(4)$R={\displaystyle \sum _{k=0}^m\frac{R_{\mathrm{k}}}{N_{\mathrm{k}}}}=R_0\frac{\left(1-{\left(\frac{R_{\mathrm{L}}}{n{R}_{\mathrm{W}}^3{R}_{\mathrm{h}}}\right)}^{m+1}\right)}{1-\frac{R_{\mathrm{L}}}{n{R}_{\mathrm{W}}^3{R}_{\mathrm{h}}}}$

where

(5)$R_0=\frac{12\mu l_0}{\beta W_{\mathrm{f}0}^3{h}_{\mathrm{f}0}}$

(6)$N_{\mathrm{k}}=n^k$

where n is the total number of branches fractal fracture, n = 2.

Then the flow rate is described by Equation (7):

(7)$Q=\frac{\Delta P}{R}=\frac{\beta W_{\mathrm{f}0}^3{h}_{\mathrm{f}0}\left(1-\frac{R_{\mathrm{L}}}{n{R}_{\mathrm{W}}^3{R}_{\mathrm{h}}}\right)\Delta P}{12\mu l_0\left(1-{\left(\frac{R_{\mathrm{L}}}{n{R}_{\mathrm{W}}^3{R}_{\mathrm{h}}}\right)}^{m+1}\right)}$

where ΔP is the total pressure difference of fractal fracture network, MPa.

(8)$\Delta P=P_{\mathrm{f}}-P_{\mathrm{wf}}$

As the fracturing fluid is produced, the fracture pressure decreases, and the fracture will be compressed under the closure stress. At this time, assuming that its height and length remain unchanged, and the width of W_fk at the k-th level has a width of W_fkc after compression, then the volume change is shown in Equation (9):

(9)$V_{\mathrm{fk}}-V_{\mathrm{fkc}}=C_{\mathrm{f}}{V}_{\mathrm{fk}}\Delta P_{\mathrm{f}}$

where C_f is fracture compression coefficient; V_fk is the k-level single fracture volume under original fracture pressure; V_fkc is the k-level single fracture volume under current fracture pressure; and ΔP_f is fracture system pressure drop, ΔP_f = P_fi − P_f.

(10)$V_{\mathrm{fk}}=W_{\mathrm{fk}}{h}_{\mathrm{fk}}{l}_{\mathrm{k}}$

(11)$V_{\mathrm{fkc}}=W_{\mathrm{fkc}}{h}_{\mathrm{fk}}{l}_{\mathrm{k}}$

where W_fkc is the width of the k-th level single fracture under the current fracture pressure.

Equations (10) and (11) are substituted into Equation (9), and then Equation (12) can be obtained.

(12)$W_{\mathrm{fkc}}=\left(1-C_{\mathrm{f}}\Delta P_{\mathrm{f}}\right)W_{\mathrm{fk}}$

The flow rate of 1/2 single-cluster tree-shaped fractal fracture network is expressed in Equation (13) by considering the fracture closure effect.

(13)$Q=\frac{\Delta P}{R}=\frac{\beta W_{\mathrm{f}0\mathrm{c}}^3{h}_{\mathrm{f}0}\left(1-\frac{R_{\mathrm{L}}}{n{R}_{\mathrm{W}}^3{R}_{\mathrm{h}}}\right)\Delta P}{12\mu l_0\left(1-{\left(\frac{R_{\mathrm{L}}}{n{R}_{\mathrm{W}}^3{R}_{\mathrm{h}}}\right)}^{m+1}\right)}$

where

(14)$W_{\mathrm{f}0\mathrm{c}}=\left(1-C_{\mathrm{f}}\Delta P_{\mathrm{f}}\right)W_{\mathrm{f}0}$

The two-phase flow depends on the relationship between gas–water relative permeability and water saturation in the fractal fracture system. In the flow model of a single-phase tree-shaped fracture network, the relative permeability is taken into account. Equation (15) illustrates the flow rate calculation model for a 1/2 single-cluster tree-shaped fracture network in gas/water two-phase flow.

(15)$Q_{\mathrm{i}}=\frac{\beta K_{\mathrm{ri}}\left(S_{\mathrm{w}}\right)W_{\mathrm{f}0}^3{h}_{\mathrm{f}0}\left(1-\frac{R_{\mathrm{L}}}{n{R}_{\mathrm{W}}^3{R}_{\mathrm{h}}}\right)\left(P_{\mathrm{f}}-P_{\mathrm{wf}}\right)}{12{\mu }_{\mathrm{i}}{B}_{\mathrm{i}}{l}_0\left(1-{\left(\frac{R_{\mathrm{L}}}{n{R}_{\mathrm{W}}^3{R}_{\mathrm{h}}}\right)}^{m+1}\right)}$

where P_f is the average pressure in the fracture system, which changes as fluid is produced, MPa; P_wf is the wellbore bottom flow pressure in a horizontal well, MPa; I is water or gas; B_i is volume coefficient, m³/m³; and K_ri(S_w) is the relative permeability of gas/water, dimensionless.

In this paper, the linear relative permeability model was adopted:

(16)$K_{\mathrm{rw}}=S_{\mathrm{w}}$

(17)$K_{\mathrm{rg}}=1-S_{\mathrm{w}}$

where S_w is the water saturation in the fracture.

Gas production superposition model is as follows:

(18)$Q_{\mathrm{g}}={\displaystyle \sum _{j=1}^{2N_{\mathrm{f}}}{Q}_{\mathrm{g}}^j}$

Water production superposition model is as follows:

(19)$Q_{\mathrm{w}}={\displaystyle \sum _{j=1}^{2N_{\mathrm{f}}}{Q}_{\mathrm{w}}^j}$

where N_f is the total clusters in horizontal well, which satisfies the following relationship:

(20)$N_{\mathrm{f}}=n_{\mathrm{f}}\cdot n_{\mathrm{CL}}$

where n_f is the number of operation stages, dimensionless; and n_CL is the clusters number, dimensionless.

2.2. Calculation of Fracture Network Volume

The volume of the fractal fracture network is

(21)$\begin{array}{l}{V}_{\mathrm{fi}}=2N_{\mathrm{f}}{\displaystyle \sum _{k=0}^m{N}_{\mathrm{k}}{V}_{\mathrm{k}}=}2N_{\mathrm{f}}{\displaystyle \sum _{k=0}^m{n}^k{W}_{\mathrm{fk}}{h}_{\mathrm{fk}}{l}_{\mathrm{k}}}\\ =2N_{\mathrm{f}}{V}_0\frac{\left(1-{\left(n{R}_{\mathrm{W}}{R}_{\mathrm{h}}{R}_{\mathrm{L}}\right)}^{m+1}\right)}{1-n{R}_{\mathrm{W}}{R}_{\mathrm{h}}{R}_{\mathrm{L}}}\end{array}$

Recently, the initial effective fracture volume V_fi has been used internationally as an important parameter for preliminary evaluation of shale gas fracture network fracturing [19,20]. In the above formula,

(22)$V_0=W_{\mathrm{f}0}{h}_{\mathrm{f}0}{l}_0$

where V_fi is initial volume of effective fracture, m³; and V₀ is initial fracture volume, m³.

It can be seen from Figure 1 that the longitudinal equivalent fracture half-length is

(23)$x_{\mathrm{f}}=l_0\left[1+\frac{R_{\mathrm{L}}\left(1-R_{\mathrm{L}}^m\right)\cos \theta }{1-R_{\mathrm{L}}}\right]$

where θ represents the bifurcation fractures, which varies from 0 to π/2.

As shown in Figure 1, the lateral extension width is given by

(24)$w_{\mathrm{f}}=l_0{R}_{\mathrm{L}}\sin \theta \frac{\left(1-R_{\mathrm{L}}^m\right)}{1-R_{\mathrm{L}}}$

The effective fracture network volume can be obtained by Equation (25):

(25)$\mathrm{EFNV}=N_{\mathrm{f}}{w}_{\mathrm{f}}{x}_{\mathrm{f}}{h}_{\mathrm{f}0}-V_{\mathrm{b}\text{\_}\mathrm{overlap}}$

where x_f denotes the longitudinal expansion degree and w_f denotes the lateral expansion degree. It is an important parameter commonly used in quantitative evaluation of shale gas fracture network fracturing [21,22], where

(26)$w_{\mathrm{f}}=l_0{R}_{\mathrm{L}}\sin \theta \frac{\left(1-R_{\mathrm{L}}^m\right)}{1-R_{\mathrm{L}}}$

(27)$V_{\mathrm{b}\text{\_}\mathrm{overlap}}=\left\{\begin{array}{cc}2S_{\mathrm{overlap}}{h}_{\mathrm{f}0}\left(N_{\mathrm{f}}-1\right)& w_{\mathrm{f}}>\frac{L_{\mathrm{w}}}{N_{\mathrm{f}}}\\ 0& w_{\mathrm{f}}\le \frac{L_{\mathrm{w}}}{N_{\mathrm{f}}}\end{array}\right.$

where V_{b_overlap} signifies the volume of overlap area within the EFNV.

In the matrix, the gas is entered the fracture network through the reverse imbibition of the fracturing fluid and the redistribution of the original free gas in the activated natural fractures during the process of fracturing and soaking. The reverse imbibition index (I_imb) is adopted to represent the combined effect. The primary origin of initial free gas in the effective fracture system is the reverse imbibition process [23]. The shale reverse imbibition index I_imb represents the proportion of free gas in the fracture network to the fracture volume, 0 ≤ I_imb ≤ 1. Obviously, there is the following relationship, the gas saturation of the original fracture is

(28)$S_{\mathrm{gi}}=I_{\mathrm{imb}}$

The volume of free gas within the subsurface fractures is

(29)$V_{\mathrm{gfi}}=I_{\mathrm{imb}}{V}_{\mathrm{fi}}$

Then the volume occupied by water within the fracture is

(30)$V_{\mathrm{wi}}=\left(1-I_{\mathrm{imb}}\right)V_{\mathrm{fi}}$

Water saturation of fractures is as follows:

(31)$S_{\mathrm{wi}}=\frac{V_{\mathrm{wi}}}{V_{\mathrm{fi}}}=1-I_{\mathrm{imb}}$

where S_gi represents the original gas saturation within fractures, dimensionless; V_gfi represents the volume of free gas within the subsurface fractures, m³; V_wi is water volume in fractures, m³; and S_wi is water saturation in fractures under initial conditions, dimensionless.

2.3. Material Balance Equation

The change of fracture volume is very complicated in the process of flow back. The pressure of the fracture system drops from the original fracture pressure P_fi to the current fracture pressure P_f, and the fracture pressure drop is ΔP_f = P_fi − P_f when a certain amount of fracturing fluid (W_p) is discharged from the fracture. The reduction of the fracture volume, the expansion of the free gas volume in the fracture, and the entry of matrix gas into the fracture will all reduce the volume of the fracturing fluid.

• (1). The reduction of fracture volume is calculated by

(32) $\Delta V_{\mathrm{f}}=V_{\mathrm{fi}}{C}_{\mathrm{f}}\Delta P_{\mathrm{f}}$

• • (a). In fracture network, the expansion of free gas can be expressed as

(33)$\Delta V_{\mathrm{gf}}=V_{\mathrm{gf}}-V_{\mathrm{gfi}}=I_{\mathrm{imb}}{V}_{\mathrm{fi}}\left(\frac{B_{\mathrm{gf}}}{B_{\mathrm{gfi}}}-1\right)$

where V_gf is the gas volume in the current fracture, m³; V_gfi is the gas volume in the original fracture, m³; B_gf is the gas volume factor in the fracture, m³/m³; B_gfi is the original gas volume factor in the fracture, m³/m³; P_fi is original fracture pressure, MPa; and P_f is current fracture pressure, MPa.

• • (b). Intrusion of matrix gas

The ground volume of intrusion of matrix gas channeling into the tree-shaped fractal fracture network is obtained using Equation (34) by considering the compressibility of matrix pores and the desorption effect of adsorbed gas in shale matrix.

(34)$G_{\mathrm{mf}}=\frac{V_{\mathrm{b}}{\varphi }_{\mathrm{m}}}{B_{\mathrm{gmi}}}-\frac{V_{\mathrm{b}}{\varphi }_{\mathrm{m}}\left(1-C_{\mathrm{m}}\Delta P_{\mathrm{m}}\right)}{B_{\mathrm{gm}}}+V_{\mathrm{b}}\left(1-{\varphi }_{\mathrm{m}}\right)\left(V_{\mathrm{Ei}}-V_{\mathrm{E}}\right)$

where B_gm and B_gmi are matrix gas volume coefficient under the condition of current matrix pressure and original matrix pressure, m³/m³; ϕ_m is matrix porosity of shale reservoir, dimensionless; C_m is rock compressibility coefficient of shale matrix, 1/MPa; V_Ei and V_E are unit shale adsorbed gas volume under original matrix pressure and current matrix pressure, m³/m³; and ΔP_m is matrix pressure drop, MPa.

V_b is the effective fracture network volume and the calculation method is shown in Equation (25). The average pressure drop of the matrix system and the desorption of matrix adsorbed gas are calculated as follows:

(35)$\Delta P_{\mathrm{m}}=P_{\mathrm{mi}}-P_{\mathrm{m}}$

(36)$V_{\mathrm{Ei}}=V_{\mathrm{L}}\frac{P_{\mathrm{mi}}}{P_{\mathrm{L}}+P_{\mathrm{mi}}}$

(37)$V_{\mathrm{E}}=V_{\mathrm{L}}\frac{P_{\mathrm{m}}}{P_{\mathrm{L}}+P_{\mathrm{m}}}$

The underground volume of matrix shale gas intrusion can be obtained as Equation (38) by combining it with the gas volume coefficient under the fracture pressure.

(38)$V_{\mathrm{mf}}=G_{\mathrm{mf}}{B}_{\mathrm{gf}}$

where P_mi is shale original matrix pressure, MPa; P_m is shale current matrix pressure, MPa; V_L is Langmuir volume, sm³/m³; P_L is Langmuir pressure, MPa; and G_mf is matrix shale surface volume of intrusion of gas channeling into the tree-shaped fractal fracture network, m³.

• • (c). Produced free gas (underground volume)

(39) $\Delta V_{\mathrm{gp}}=G_{\mathrm{p}}{B}_{\mathrm{gf}}$

The gas storage capacity of fractures is obtained by combining Equations (33), (38), and (39).

(40)$\Delta V_{\mathrm{g}}=\Delta V_{\mathrm{gf}}+V_{\mathrm{mf}}-\Delta V_{\mathrm{gp}}$

where G_p is cumulative gas production, m³.

• (2). The remaining volume of fracturing fluid in the fracture system

The volume of fracturing fluid in fractures will all be reduced due to the reduction of fracture pore volume, the expansion of fracturing fluid in fractures, and the gas storage capacity in fractures. Therefore, when the fracture original pressure P_fi decreases to P_f, the fracture fracturing fluid volume is

(41)$V_{\mathrm{w}}=V_{\mathrm{wi}}-\Delta V_{\mathrm{f}}-\Delta V_{\mathrm{g}}$

Then, Equation (42) can be obtained by combining Equations (32), (40), and (41):

(42)$V_{\mathrm{w}}=V_{\mathrm{wi}}-V_{\mathrm{fi}}{C}_{\mathrm{f}}\Delta P-I_{\mathrm{imb}}{V}_{\mathrm{fi}}\left(\frac{B_{\mathrm{gf}}}{B_{\mathrm{gfi}}}-1\right)-G_{\mathrm{mf}}{B}_{\mathrm{gf}}+G_{\mathrm{p}}{B}_{\mathrm{gf}}$

The remaining fracturing fluid volume is converted to surface conditions:

(43)$W_{\mathrm{res}}=\frac{V_{\mathrm{w}}}{B_{\mathrm{wf}}}=\left(\begin{array}{l}{V}_{\mathrm{wi}}-V_{\mathrm{fi}}{C}_{\mathrm{f}}\Delta P\\ -I_{\mathrm{imb}}{V}_{\mathrm{fi}}\left(\frac{B_{\mathrm{gf}}}{B_{\mathrm{gfi}}}-1\right)-G_{\mathrm{mf}}{B}_{\mathrm{gf}}+G_{\mathrm{p}}{B}_{\mathrm{gf}}\end{array}\right)/B_{\mathrm{wf}}$

The basic form of the fracturing fluid material balance equation in the fracture is as follows:

(44)$\frac{V_{\mathrm{wi}}}{B_{\mathrm{wi}}}=W_{\mathrm{p}}+W_{\mathrm{res}}$

where B_wi is the fracturing fluid volume coefficient under original pressure, dimensionless; W_p is current cumulative fracturing fluid volume of flow back m³; and W_res is the remaining fracturing fluid volume, m³.

The remaining fracturing fluid volume W_res is inputted into the above formula:

(45)$\frac{V_{\mathrm{wi}}}{B_{\mathrm{wi}}}=W_{\mathrm{p}}+\left(\begin{array}{l}{V}_{\mathrm{wi}}-V_{\mathrm{fi}}{C}_{\mathrm{f}}\Delta P\\ -I_{\mathrm{imb}}{V}_{\mathrm{fi}}\left(\frac{B_{\mathrm{gf}}}{B_{\mathrm{gfi}}}-1\right)-G_{\mathrm{mf}}{B}_{\mathrm{gf}}+G_{\mathrm{p}}{B}_{\mathrm{gf}}\end{array}\right)/B_{\mathrm{wf}}$

Then,

(46)$\begin{array}{l}{W}_{\mathrm{p}}{B}_{\mathrm{wf}}+G_{\mathrm{p}}{B}_{\mathrm{gf}}=V_{\mathrm{wi}}\left(\frac{B_{\mathrm{wf}}-B_{\mathrm{wi}}}{B_{\mathrm{wi}}}\right)\\ +V_{\mathrm{fi}}{C}_{\mathrm{f}}\Delta P_{\mathrm{f}}+I_{\mathrm{imb}}{V}_{\mathrm{fi}}\left(\frac{B_{\mathrm{gf}}}{B_{\mathrm{gfi}}}-1\right)+G_{\mathrm{mf}}{B}_{\mathrm{gf}}\end{array}$

where B_wf is the fracturing fluid volume coefficient under the current fracture pressure, m³/m³; and B_wfi is the fracturing fluid volume coefficient under the original fracture pressure, m³/m³.

The compressibility coefficients of gas and water are

(47)$\begin{array}{l}{C}_{\mathrm{g}}=\frac{B_{\mathrm{gf}}-B_{\mathrm{gfi}}}{B_{\mathrm{gfi}}\Delta P_{\mathrm{f}}}\\ C_{\mathrm{w}}=\frac{B_{\mathrm{wf}}-B_{\mathrm{wfi}}}{B_{\mathrm{wfi}}\Delta P_{\mathrm{f}}}\end{array}$

Then,

(48)$W_{\mathrm{p}}{B}_{\mathrm{wf}}+G_{\mathrm{p}}{B}_{\mathrm{gf}}=V_{\mathrm{fi}}\Delta P\left[\left(1-I_{\mathrm{imb}}\right)C_{\mathrm{w}\mathrm{f}}+C_{\mathrm{f}}+I_{\mathrm{imb}}{C}_{\mathrm{g}}\right]+G_{\mathrm{mf}}{B}_{\mathrm{gf}}$

The function h representing the fracture pressure and matrix pressure at the k + 1th time step is defined as follows:

(49)$\begin{array}{l}h\left(P_{\mathrm{f}}^{k+1},P_{\mathrm{m}}^{k+1}\right)=\\ V_{\mathrm{fi}}\left(P_{\mathrm{fi}}-P_{\mathrm{f}}^{k+1}\right)\left(\left(1-I_{\mathrm{imb}}\right)C_{\mathrm{w}\mathrm{f}}+C_{\mathrm{f}}+I_{\mathrm{imb}}{C}_{\mathrm{g}}\right)+G_{\mathrm{mf}}{B}_{\mathrm{gf}}\left(P_{\mathrm{f}}^{k+1}\right)\\ -W_{\mathrm{p}}{B}_{\mathrm{wf}}\left(P_{\mathrm{f}}^{k+1}\right)-G_{\mathrm{p}}{B}_{\mathrm{gf}}\left(P_{\mathrm{f}}^{k+1}\right)\end{array}$

where

(50)$B_{\mathrm{gf}}\left(P_{\mathrm{f}}^{k+1}\right)=\frac{p_{\mathrm{sc}}}{Z_{\mathrm{sc}}{T}_{\mathrm{sc}}}\frac{Z_{\mathrm{f}}{T}_{\mathrm{i}}}{P_{\mathrm{f}}^{k+1}}$

where C_wf is fracturing fluid compression coefficient in fractures, dimensionless; P_sc is ground pressure, MPa; Z_sc is gas deviation factor in ground conditions, dimensionless; T_sc is ground temperature, K; Z_f is gas deviation factor in fractures, dimensionless; and T_i is original formation temperature of reservoir, K.

Water saturation of fractures under current formation conditions is as follows:

(51)$S_{\mathrm{w}}=\frac{V_{\mathrm{w}}}{V_{\mathrm{w}}+V_{\mathrm{gfi}}+\Delta V_{\mathrm{g}}}$

The equation can be simplified to

(52)$S_{\mathrm{w}}=\frac{V_{\mathrm{w}}}{V_{\mathrm{fi}}-\Delta V_{\mathrm{f}}}$

The equations for fracture and matrix pressure at time k + 1 can be derived by setting Equation (49) equal to zero.

(53)$h\left(P_{\mathrm{f}}^{k+1},P_{\mathrm{m}}^{k+1}\right)=0$

The material balance equation needs to be established of the shale matrix system to solve this equation due to there being two unknowns, P_f^k+1 and P_m^k+1.

If the fracture pressure (P_f) falls below the matrix gas breakthrough pressure (P_BT), it will result in channeling between the matrix and fractures. The pressure P_BT is related to the porosity and permeability of the matrix. The channeling equation [24] is described by Equation (54).

(54)$q_{\mathrm{m}}=\frac{k_{\mathrm{m}}{\alpha }_{\mathrm{mf}}\left(P_{\mathrm{m}}-P_{\mathrm{f}}\right)}{{\mu }_{\mathrm{g}}}$

If P_BT ≤ P_f^k, then the channeling from the matrix to the fracture system from time k to k + 1 is 0:

(55)$\Delta G_{\mathrm{mf}}=0$

If P_BT > P_f^k, the channeling from the matrix to the fracture system from time k to k + 1 is

(56)$\Delta G_{\mathrm{mf}}=\Delta t{V}_{\mathrm{b}}{\varphi }_{\mathrm{m}}{q}_{\mathrm{m}}=\Delta t{V}_{\mathrm{b}}\frac{{\varphi }_{\mathrm{m}}{k}_{\mathrm{m}}{\alpha }_{\mathrm{mf}}\left(P_{\mathrm{m}}^k-P_{\mathrm{f}}^{k+1}\right)}{{\mu }_{\mathrm{g}}}$

where ΔG_mf is the channeling of gas from the shale matrix to the fracture system from time k to time k + 1, m³; Δt—the time from time k to time k + 1, s; V_b is shale reservoir stimulation volume, m³; ϕ_m is matrix porosity of reservoir, dimensionless; q_m is gas channeling flow from matrix to fracture system per unit time, s⁻¹; k_m is matrix permeability, um²; α_mf is matrix to fracture channeling factor [25], m⁻²; P_m^k is average pressure of shale matrix system at time k, MPa; P_f^k+1 is average pressure of fracture system at time k + 1, MPa; and u_g is shale gas viscosity, Pa·s.

The channeling flow diffusion from the matrix to the fracture system from time k to k + 1 can also be expressed as

(57)$\Delta G_{\mathrm{mf}}=G_{\mathrm{mf}}^{k+1}-G_{\mathrm{mf}}^k$

where

(58)$\begin{array}{l}{G}_{\mathrm{mf}}^k=\frac{V_{\mathrm{b}}{\varphi }_{\mathrm{m}}}{B_{\mathrm{gmi}}}-\frac{V_{\mathrm{b}}{\varphi }_{\mathrm{m}}\left(1-C_{\mathrm{m}}\left(P_{\mathrm{mi}}-P_{\mathrm{m}}^k\right)\right)}{B_{\mathrm{gm}}^k}\\ +V_{\mathrm{b}}\left(1-{\varphi }_{\mathrm{m}}\right)\left(V_{\mathrm{Ei}}-V_{\mathrm{E}}\left(P_{\mathrm{m}}^k\right)\right)\end{array}$

(59)$\begin{array}{l}{G}_{\mathrm{mf}}^{k+1}=\frac{V_{\mathrm{b}}{\varphi }_{\mathrm{m}}}{B_{\mathrm{gmi}}}-\frac{V_{\mathrm{b}}{\varphi }_{\mathrm{m}}\left(1-C_{\mathrm{m}}\left(P_{\mathrm{mi}}-P_{\mathrm{m}}^{k+1}\right)\right)}{B_{\mathrm{gm}}^{k+1}}\\ +V_{\mathrm{b}}\left(1-{\varphi }_{\mathrm{m}}\right)\left(V_{\mathrm{Ei}}-V_{\mathrm{E}}\left(P_{\mathrm{m}}^{k+1}\right)\right)\end{array}$

(60)$V_{\mathrm{E}}=\frac{V_{\mathrm{L}}P}{P+P_{\mathrm{L}}}$

where V_L is Langmuir volume, sm³/m³; P_L—Langmuir pressure, MPa. B_gm and B_gmi are matrix gas volume coefficient under current matrix pressure and original matrix pressure, m³/m³; C_m is rock compressibility coefficient of shale matrix, 1/MPa; and V_Ei and V_E are the unit shale adsorbed gas volume under original matrix pressure and current matrix pressure, m³/m³.

Equations (58) and (59) are substituted into Equation (57), then Equation (61) can be obtained.

(61)$\begin{array}{l}\Delta G_{\mathrm{mf}}=G_{\mathrm{mf}}^{k+1}-G_{\mathrm{mf}}^k\\ =\frac{V_{\mathrm{b}}{\varphi }_{\mathrm{m}}\left(1-C_{\mathrm{m}}\left(P_{\mathrm{m}}-P_{\mathrm{m}}^k\right)\right)}{B_{\mathrm{gm}}^k}\\ -\frac{V_{\mathrm{b}}{\varphi }_{\mathrm{m}}\left(1-C_{\mathrm{m}}\left(P_{\mathrm{m}}-P_{\mathrm{m}}^{k+1}\right)\right)}{B_{\mathrm{gm}}^{k+1}}\\ +V_{\mathrm{b}}\left(1-{\varphi }_{\mathrm{m}}\right)\left(V_{\mathrm{E}}\left(P_{\mathrm{m}}^k\right)-V_{\mathrm{E}}\left(P_{\mathrm{m}}^{k+1}\right)\right)\end{array}$

The function g representing the fracture and matrix pressure at the k + 1 time step is derived from Equation (62) by combining Equations (56) and (61):

(62)$\begin{array}{l}g\left(P_{\mathrm{m}}^{k+1},P_{\mathrm{f}}^{k+1}\right)=\\ \frac{{\varphi }_{\mathrm{m}}\left(1-C_{\mathrm{m}}\left(P_{\mathrm{mi}}-P_{\mathrm{m}}^k\right)\right)}{B_{\mathrm{gm}}^k}\\ -\frac{{\varphi }_{\mathrm{m}}\left(1-C_{\mathrm{m}}\left(P_{\mathrm{mi}}-P_{\mathrm{m}}^{k+1}\right)\right)}{B_{\mathrm{gm}}^{k+1}}\\ +\left(1-{\varphi }_{\mathrm{m}}\right)\left(V_{\mathrm{E}}\left(P_{\mathrm{m}}^k\right)-V_{\mathrm{E}}\left(P_{\mathrm{m}}^{k+1}\right)\right)\\ -\Delta t\frac{{\varphi }_{\mathrm{m}}{k}_{\mathrm{m}}{\alpha }_{\mathrm{mf}}\left(P_{\mathrm{m}}^k-P_{\mathrm{f}}^{k+1}\right)}{{\mu }_{\mathrm{g}}}\end{array}$

The equations for the fracture network and matrix pressure at time k + 1 can be derived by setting Equation (62) to zero:

(63)$g\left(P_{\mathrm{f}}^{k+1},P_{\mathrm{m}}^{k+1}\right)=0$

There are two unknowns, P_f^k+1 and P_m^k+1, in above equation, and the fracture pressure P_f^k+1 and matrix pressure P_m^k+1 at time k + 1 can be obtained by solving Equations (53) and (63), and the fracture pressure P_f^k and matrix pressure P_m^k at time k are known.

2.4. Model Solution

• (1). The tree-shaped fracture network structure parameters are preset, including l₀, W_f0, h_f0, R_L, R_w, R_h, θ, m, n, and C_f, and initial fracture network and matrix average pressure (P_f(k = 1) = P_fi, P_m(k = 1) = P_mi) and the bottom hole flow pressure P_wf, then Q_W(k) and Q_g(k) is calculated by Equation (15), k = 1, 2, …, Num;

• (2). The cumulative production of step k is solved, the cumulative flow back W_p = sum(Q_w(k)), and the cumulative gas production G_p = sum(Q_g(k));

• (3). If P_f(k) > =P_BT, then there is no matrix gas channeling into the effective fracture network, G_mf = 0, and then P_m(k + 1) = P_m(k). Based on the dichotomy method, the fracture pressure P_f(k + 1) is calculated by combining the fracture material balance Equation (49), and the network average pressure P_f(k + 1) ranges from [0, P_fi]. When P_f(k) < P_BT, the matrix gas channeling into fractures, the fracture pressure P_f(k + 1) range [0, P_fi], under the condition of known P_f(k + 1), the matrix pressure P_m(k + 1) range [0, P_mi]. The matrix pressure P_m(k + 1) is calculated in combination with the matrix material balance Equation (62), and then G_mf is calculated by using Equation (34), which is then brought into the Equation (49), and P_f(k + 1) is solved by the dichotomy method.

• (4). The P_f(k + 1) and P_m(k + 1) are assigned to P_f(k) and P_m(k), and steps (1), (2), and (3) are repeated until k = Num.

2.5. Flow Back Model Validation

The rate transient analysis (RTA) method was used to identify and calculate the flow back flow stage of early shale fracturing. Since a constant flow rate or bottom hole flow pressure was rarely used in the flow back, a small oil nozzle is usually used to control the discharge, step by step amplification, and adjust a stable drainage system in southern Sichuan. Therefore, the rate normalized pressure (RNP) and its derivative (RNP’) were used to observe the flow back characteristics. The equations are as follows [26]:

(64)$\mathrm{RNP}=\frac{P_{\mathrm{fi}}-P_{\mathrm{wf}}}{q_{\mathrm{w}}}$

(65)${\mathrm{RNP}}^{\text{'}}=\frac{\mathrm{dRNP}}{\mathrm{dln}{t}_{\mathrm{MB}}}$

(66)$t_{\mathrm{MB}}=\frac{W_{\mathrm{P}}}{q_{\mathrm{w}}}$

where P_fi is initial fracture network pressure, MPa; P_wf is bottom hole flow pressure, MPa; q_w is water production rate, m³/h; t_MB is equivalent time, h; and W_p is cumulative water production, m³.

The established shale gas fracturing two-phase flow back model was used to calculate the flow back data of 1000 h. The model parameter settings are shown in Table 2. As shown in Figure 2, in the initial state, the effective fracture network volume is completely filled with water (I_imb = 0) and the observed flow stages are consistent with field data [27] and analytical models [26] and Williams-Kovacs [28]. For the case of I_imb = 0, the two phases in RNP’ are characterized by a slope equal to 1. If the free gas is contained in the effective fracture network volume in the initial state, the more free gas content, the smaller the slope in stage 2 (for example, the slope of 0.70 at I_imb = 0.1 is greater than the slope of 0.65 at I_imb = 0.3), and the slope of stage 3 is always 1. Stage 2 and stage 3 last for a long time, and the characteristics of these two stages can be fully reflected by the flow back model proposed in this paper, so the established model is reasonable.

Some assumptions were made on the real seepage state of the fracturing fluid in order to be used in the inversion of the effective fracture network volume, so that the model will deviate from the real characteristics of the fluid flow back stage. But because the permeability of the effective fracture network is very large, the pressure wave propagates very fast, and this flow stage is usually impossible to observe in traditional low-frequency flow back data (transient flow in the fracture system may last only 10 s [29], and the absence of this stage has little error in the inversion). Stage 2 and stage 3 are the main stages in the production process of flow back, and the model proposed in this paper can fully reflect the characteristics of these two stages. It shows that the flow back model proposed in this paper is reasonable.

3. EFNV Inversion

In this section, a mathematical model was developed to invert the volume of an effective fracture network in shale gas using a tree-shaped fractal fracture flow back approach. Then, based on the flow back data of shale gas well, a set of genetic algorithm workflows suitable for the effective fracture network volume inversion of shale gas is developed and combined with a genetic algorithm. An example was taken from shale fracturing wells in southern Sichuan, where the inversion of effective fracture network volume was conducted for 162 wells. A basic data set of fracture network volume was obtained.

3.1. Statistical Evaluation Methods

The following four statistical descriptors are introduced to evaluate the reliability of the fitting between the flow back data (Q_wsc and Q_gsc) and in the predicted data of the model (Q_w and Q_g), y_i is the observed value, f_i is the predicted value, and N is the number of samples.

• (1). Coefficient of determination (R²)

The closer the R² value is to 1, the better the predictive performance of the model.

(67)$R^2=1-\frac{{\displaystyle \sum _{i=1}^N{\left(y_i-f_i\right)}^2}}{{\displaystyle \sum _{i=1}^N{\left(y_i-\overline{y}\right)}^2}}$

(68)$\overline{y}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{y}_i}$

• (2). Root mean square error (RMSE)

RMSE is the standard deviation of the error between the observed values (Q_wsc and Q_gsc) and the predicted values (Q_w and Q_g). The smaller the RMSE, the better the prediction performance of the model.

(69)$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(y_i-f_i\right)}^2}}$

• (3). Slope of regression line (K)

The K value is the slope of the fitted linear regression equation between the observed values (Q_wsc and Q_gsc) and predicted values (Q_w and Q_g).

(70)$K=\frac{{\displaystyle \sum _{i=1}^N{y}_i{f}_i}}{{\displaystyle \sum _{i=1}^N{\left(f_i\right)}^2}}$

• (4). Willmott’s index of agreement (IA)

The range of IA is [0,1]: IA = 1 indicates that the forecast is perfect, and IA = 0 indicates that the forecast is invalid.

(71)$IA=1-\frac{{\displaystyle \sum _{i=1}^N{\left(y_i-f_i\right)}^2}}{{\displaystyle \sum _{i=1}^N{\left(\left|f_i-\overline{y}\right|+\left|y_i-\overline{y}\right|\right)}^2}}$

According to statistical recommendations, a reliable result can be evaluated by R² > 0.64, 0.85 < K < 1.15, or IA > 0.80 [30].

3.2. Fitness Function

The goal was to find the most suitable fracture network structure through flow back data inversion, that is, the flow back water and gas production predicted by the model should be the same or very close to the actual monitored flow back data. Therefore, the fitness function was established based on the idea of maximizing the determination coefficient R² between the predicted and the observed value. The fitness function is as follows:

(72)$F\left(\overrightarrow{x}\right)=\left\{\begin{array}{cc}f\left(\overrightarrow{x}\right)& \overrightarrow{x}\in feasible\hspace{1em}\hspace{1em}region\\ f\left(\overrightarrow{x}\right)+penalty\left(\overrightarrow{x}\right)& \overrightarrow{x}\notin feasible\hspace{1em}\hspace{1em}region\end{array}\right.$

where

(73)$f\left(\overrightarrow{x}\right)=\frac{{\displaystyle \sum _{i=1}^N{\left(Q_{\mathrm{wsci}}-Q_{\mathrm{wpredi}}\right)}^2}}{{\displaystyle \sum _{i=1}^N{\left(Q_{\mathrm{wsci}}-{\overline{Q}}_{\mathrm{wsci}}\right)}^2}}+\frac{{\displaystyle \sum _{i=1}^N{\left(Q_{\mathrm{gsci}}-Q_{\mathrm{gpredi}}\right)}^2}}{{\displaystyle \sum _{i=1}^N{\left(Q_{\mathrm{gsci}}-{\overline{Q}}_{\mathrm{gsci}}\right)}^2}}$

(74)$penalty\left(\overrightarrow{x}\right)=\{\begin{array}{cc}{10}^8& p_{\mathrm{f}}\mathrm{or}{p}_{\mathrm{m}}\notin feasibleregion\\ 0& p_{\mathrm{f}}\mathrm{or}{p}_{\mathrm{m}}\in feasibleregion\end{array}$

The decision variables were as follows:

(75)$\overrightarrow{x}=\left(x_1,x_2,\cdot \cdot \cdot ,x_{12}\right)$

where

(76)$\begin{array}{cccc}\begin{array}{c}{x}_1=l_0\\ x_5=R_{\mathrm{h}}\\ x_9=p_{\mathrm{fi}}\end{array}& \begin{array}{c}{x}_2=W_{\mathrm{f}0}\\ x_6=R_{\mathrm{L}}\\ x_{10}={\alpha }_{\mathrm{mf}}\end{array}& \begin{array}{c}{x}_3=m\\ x_7=\theta \\ x_{11}=C_{\mathrm{f}}\end{array}& \begin{array}{c}{x}_4=R_{\mathrm{W}}\\ x_8=p_{\mathrm{BT}}\\ x_{12}=I_{\mathrm{imb}}\end{array}\end{array}$

The optimization goal was as follows:

(77)$\mathrm{minimize}F\left(\overrightarrow{x}\right)$

The variable upper and lower bounds were as follows:

(78)$\begin{array}{cc}L{B}_{\mathrm{i}}\le x_{\mathrm{i}}\le U{B}_{\mathrm{i}}& \left(i=1,2,\cdot \cdot \cdot ,12\right)\end{array}$

The constraints were as follows:

(79)$V_{\mathrm{fi}}\le TIV$

where

(80)$V_{\mathrm{fi}}=x_1{x}_2{h}_{\mathrm{f}0}\frac{\left(1-{\left(n{x}_4{x}_5{x}_6\right)}^{x_3+1}\right)}{1-n{x}_4{x}_5{x}_6}$

where Penalty(x) is the penalty function, which is taken as 10⁸ in this paper; Q_wsci is the flow back fluid volume actually measured at the i-th moment, m³/min; Q_wpredi is the flow back fluid volume calculated by the flow back model at the i-th moment, m³/min; Q_gsci is the gas production volume actually measured at the i-th moment, m³/min; Q_gpredi is the gas production volume calculated by the flow back model at the i-th time moment, m³/min; ${\overline{Q}}_{\mathrm{wsci}}$ is the average flow back fluid volume actually measured, m³/min; ${\overline{Q}}_{\mathrm{gsci}}$ is the average gas production volume actually measured by the shale gas well, m³/min; l₀ is the initial fracture length, m; W_f0 is the initial fracture width, m is the number of fractures; dimensionality; R_W is the fracture width ratio, dimensionless; R_h is the fracture height ratio, dimensionless; R_L is the fracture length ratio, dimensionless; θ is the fracture bifurcation angle, rad; p_BT is the matrix gas breakthrough pressure, MPa; p_fi is the initial fracture pressure, MPa; α_mf is the channeling factor, m⁻²; C_f is fracture compression coefficient, 10⁻⁹ Pa⁻¹; I_imb is the shale reverse imbibition index, dimensionless; LB_i is the lower limit of fitting parameters for flow back data inversion; UB_i is the upper limit of fitting parameters for flow back data inversion; V_fi is the initial effective fracture volume, m³; TIV is the total injected fracturing fluid volume, m³; h_f0 is the initial height of tree-shaped fractal fractures, m; and n is the bifurcation number of the tree-shaped fractal fracture, dimensionless.

3.3. Genetic Algorithm Workflow

The genetic algorithm (GA) was established by Holland [31] and inspired by natural selection in biological evolution. Three basic operators were used to generate high-quality solutions for optimization and search problems: selection, crossover, and mutation. The GA is used to optimize a given objective function (fitness function). The inversion work was conducted using shale gas well flow back data to minimize errors between the actual flow back data (Q_wsc and Q_gsc) and the predicted data from the fracture network model (Q_w and Q_g). Therefore, the fitness function included the optimization of the dual objectives of gas production and water production.

The genetic algorithm was used to predict the fracture network. The workflow for inversion of effective fracture network volume is shown in Figure 3. The inversion simulation process is as follows:

Step 1: The inversion simulation starts. Genetic algorithm parameters are shown in Table 3, and the upper and lower limits of the flow back model parameters are shown in Table 4;

Step 2: The basic parameters of shale gas wells are inputted, including flow back data, fracturing operation parameters, initial reservoir parameters, etc.;

Step 3: The bottom hole flowing pressure is calculated as the input value through the Beggs–Brill model [32];

Step 4: An initial total group of genetic algorithm fitting parameters is formed;

Step 5: The predicted values (Q_w and Q_g) are calculated through the flow back model;

Step 6: The fitness function value is calculated;

Step 7: Judging whether all the best fitting parameter values are found, or whether the termination condition is met;

Step 8: If the judgment result is “yes”, then the fitting parameters are stored; if the judgment result is “no”, then genetic algorithm selection operation, crossover operation, and mutation operation are performed to generate the next generation population of fitting parameters, and from step 5 one starts recounting.

3.4. Model Validation

3.4.1. Analysis of Production

The shale gas wells in southern Sichuan were evaluated by the established EFNV inversion method, and the EFNV data set of the block was obtained. The average EFNV was 1497.8 × 10⁴ m³, the standard deviation was 517.0 × 10⁴ m³. The relationship between EFNV and test production, 60-day cumulative gas production, 90-day cumulative gas production, 120-day cumulative gas production, and 20-year technically recoverable reserves were analyzed.

The relationships between the EFNV obtained from the inversion and the test production, 60-day cumulative gas production, 90-day cumulative gas production, 120-day cumulative gas production, and 20-year technically recoverable reserves are shown in Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8. It can be seen that the correlation decreases in the order (test production (R = 0.684) > 60-day cumulative gas production (R = 0.576) > 90-day cumulative gas production (R = 0.566) > 120-day cumulative gas production (R = 0.541) > 20-year technically recoverable reserves (R = 0.371)). The reason is that the EFNV obtained based on the flow back in this paper was mainly obtained from the inversion of the data within the first few hundred hours to two thousand hours of the flow back production, so its correlation with the test production is the highest.

The longer the production time is, the correlation with cumulative gas production gradually decreases. This is due to the limitations of the flow back data. However, the positive correlation between the effective fracture network volume and test production, 60-day cumulative gas production, 90-day cumulative gas production, 120-day cumulative gas production, and 20-year technically recoverable reserves (EUR) shows that if a larger effective fracture network volume is obtained, the probability of obtaining a larger gas production is greater. From the perspective of fracturing engineering, which can be used as an important parameter to evaluate fracturing effects and production prediction, this proves the reliability of the inversion results.

3.4.2. Microseismic Monitoring

Some wells were monitored by microseismic monitoring during the fracturing process in southern Sichuan. The relationship between the stimulated reservoir volume (SRV) of 20 microseismic monitoring wells and the EFNV obtained by flow back inversion was statistically analyzed. As shown in Figure 9, the correlation coefficient between them was R = 0.705, which belongs to a moderate positive correlation (0.5 < R < 0.8). The more fractured the reservoir detected by microseismic monitoring, the larger the SRV, and the greater the probability that the corresponding effective fracture network volume will be larger, which further indicates that the EFNV inversion method is reliable based on flow back, and the obtained EFNV can be used as an important parameter for evaluating fracturing effects and optimizing fracturing operation parameters.

4. Field Application

We gathered flow back data post-fracturing to input into the established tree-shaped fractal fracture flow back mathematical model for EFNV inversion. Then, we utilized the genetic algorithm workflow to process the data and predict fracture network parameters for optimization. We applied the inversion model to shale gas wells in the field, such as those in southern Sichuan Basin, to calculate the EFNV and key parameters for optimizing fracturing operations.

The fracturing fluid flow back analysis model and the EFNV inversion workflow were used for a field application. The flow back production data of a shale gas well H3 in southern Sichuan (Figure 10) was taken as an example. The statistics of the basic parameters are shown in Table 5.

The hydraulic fracturing fracture network characteristics inversion was performed on the flow back production data of Well H3 based on the proposed flow back model and the flow back data inversion genetic algorithm. The water transient history fitting of well H3 is shown in Figure 11. Among them, R²(Q_w) = 0.883, R²(Q_g) = 0.927, IA(Q_w) = 0.969, IA(Q_g) = 0.982, K(Q_w) = 1.03, and K(Q_g) = 0.95; according to the statistical recommendation standard, the inversion of fracture characteristics of Well H3 was reliable. The root mean square error RMSE of the predicted value and the measured value of Q_w and Q_g were 1.54 m³/h and 0.14 × 10⁴ m⁴/h, respectively. The changes in the bottom hole flow pressure, fracture network pressure, shale matrix pressure, and wellhead nozzle size during the flow back process are shown in Figure 12. It can be seen that the pressure drop in the fracture network is larger than that in the shale matrix after the well is opened. After the first well shut-in, the gas in the matrix channeled into the fractures due to the pressure difference between the fracture network system and the matrix system. In the network, the fracture network pressure gradually rises to be equal to the matrix pressure. After the second well opening, the fracture network pressure begins to decrease again until the second well shut-in, and the fracture pressure gradually rises again. The network pressure did not rise back to equal the matrix pressure due to the short shut-in time of the second well. It can be clearly seen from Figure 11 that the fitting effect of the flow back water data is average in the early part. It is speculated that the possible reason is that the flow back data records on site are not accurate due to frequent changes in the size of the choke nozzle, and the size of the choke nozzle changed four times before the first shut-in and changed nine times before the second shut-in (Figure 12).

The fracture characteristics of Well H3 obtained through the genetic algorithm workflow inversion are shown in Table 6, including tree-shaped fractal fracture network structure parameters, fracture system original pressure, matrix gas breakthrough pressure, fracture compressibility coefficient, reverse imbibition index, effective fracture network volume, etc. Two types of hydraulic fracture networks will be formed after fracture network fracturing. Some fractures contribute to flow conduction and others are isolated due to fracture closure and for other reasons. The fracture network with a conductivity contribution is a highly conductive fracture network that allows matrix shale gas to enter and flow into the wellbore along it; these are also called effective fractures. A law of harmonic decline (HD) from the relationship between water production and cumulative water production can be observed by Lee [33], which can estimate the effective fracture volume (Table 7 and Figure 13). It can be seen from Figure 13 that when the HD model was applied to Well H3, the previous water flow data were overestimated, so the effective initial fracture volume would be overestimated.

Alkouh [2] observed a linear relationship between flow normalized pressure (RNP) and mass balance time (t_MB) during the flow back process, and the effective initial fracture volume could be obtained (Table 7 and Figure 14). However, the fracture compressibility coefficient was ignored in the Alkouh model when the total compressibility coefficient was calculated, so the effective fracture volume would be overestimated.

As shown in Figure 15, the effective initial fracture volume calculated by the method proposed in this paper, in which the gas-water two-phase flow and fracture compressibility were taken into account, and the calculation results were of the same order of magnitude as those of the HD model and the Alkouh model, and were smaller than those of the two models. It shows that the inversion model of fracture network characteristics based on flow back is reasonable, and the applicability of the model is stronger. Shale gas fracturing generally has the problem of a low flow back rate [34], and the final estimated flow back rate of Well H3 was 20.8%. Most of the fracturing fluid is locked in the isolated ineffective fracture network with no conductivity, and it is difficult to flow back. This is the main reason for the low flow back rate.

The inversion effective fracture network volume fitting curves of some shale gas fracture network fracturing wells in southern Sichuan are shown in Figure 16. The model established can fit the actual field data very well in this paper, and it can be used to invert the EFNV efficiently after shale gas reservoir stimulation in southern Sichuan.

5. Conclusions

The inversion method of shale gas effective fracture network volume based on flow back data was established. Several effective fracture network volume inversion cases were simulated. The main conclusions are summarized below:

• (1). A two-phase flow back model of the tree-shaped fractal fracture network was established based on the flow back data, and the reliability of the model was verified by combining the RNP and RNP’. The model is reasonable, reliable, and has strong field adaptability.

• (2). Based on the principle of the genetic algorithm (GA), the GA workflow of EFNV inversion was established and combined with the two-phase flow back model, and then the EFNV inversion was carried out for shale gas wells in southern Sichuan. The results showed that the inversion of EFNV had a positive correlation with production and microseismic monitoring results, which confirmed the reliability of the inversion results.

• (3). The calculation showed that the results were in the same order of magnitude as those of the HD model and the Alkouh model, and were smaller than the calculation results of the two models, indicating that the fracture network inversion model proposed is reasonable and more applicable.

• (4). The fracture network feature inversion method can obtain tree-shaped fractal fracture network structure parameters and effective fracture network volume. It has important theoretical guiding significance and field application value for improving the calculation accuracy of EFNV and fracturing effect evaluation.

• (5). The model does not consider the effects of interference between wells and well channeling from adjacent wells, and accordingly the prediction of water production after well strings is relatively weak. It is recommended that the model can further consider the impact of simultaneous flow back of multiple wells and layers and the salinity of flow back water in the future.

Author Contributions

D.T.: conceptualization, methodology; J.W.: writing—review and editing; J.Z.: writing—original draft preparation; B.Z.: data collection; Y.S.: visualization, investigation; C.S.: formal analysis; L.R.: visualization; Y.H.: original draft preparation; Z.W.: software. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Acknowledgments

Thanks to all the authors for their contribution to the paper and their efforts.

Conflicts of Interest

Authors Dengji Tang, Jianfa Wu, Bo Zeng, Yi Song, Cheng Shen, and Yongzhi Huang are employed by the CNPC Southwest Oil and Gas Field Company; The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Tree-shaped fractal network of 1/2 single-cluster.

Figure 2. The relationship between RNP’, RNP, and tMB for flow back of shale gas wells.

Figure 3. Genetic algorithm workflow for inversion of effective fracture network volume.

Figure 4. Relationship between effective fracture network volume and test production.

Figure 5. Relationship between effective fracture network volume and 60-day cumulative gas production.

Figure 6. Relationship between effective fracture network volume and 90-day cumulative gas production.

Figure 7. Relationship between effective fracture network volume and 120-day cumulative gas production.

Figure 8. Relationship between effective fracture network volume and 20-year technically recoverable reserves.

Figure 9. Relationship between inversion EFNV and microseismic monitoring SRV.

Figure 10. Flow back production data of well H3.

Figure 11. The transient performance of water and gas fitted of well H3.

Figure 12. Transient performance of pressure and choke size of well H3.

Figure 13. HD model.

Figure 14. Alkouh model [2].

Figure 15. Comparison of initial effective fracture network volume estimated by different models.

Figure 16. Effective fracture network volume inversion fitting curve.

Figure 16. Effective fracture network volume inversion fitting curve.

Table 1

Summary of current flow back models.

Author	Year	Model	Work Summary	Inputs	Outputs
Williams-Kovacs and Clarkson [10]	2013	Two-phase flow back of multi-fractured horizontal wells stochastic mode.	Using high-frequency fluid production and flowing pressures (hourly or greater) to characterize hydraulic fracture or reservoir parameters.	Fractured stages, perforation clusters, et al.	Bulk permeability, effective fracture half-length, and production.
Alkouh et al. [2]	2014	Water flow back and long-term production simulated mode.	Procedures and examples are presented, including water flow back and water-production data, in the analysis of shale-gas wells using rate transient analysis.	Initial pressure, fracture porosity, etc.	Gas production and water production.
Fu et al. [13]	2015	A flowback model.	How flow back data can be interpreted to estimate effective fracture pore-volume and its relationship to fracture design parameters are illustrates.	Flow back rate, pressure date et al.	Effective fracture pore volume.
Ezulike et al. [12]	2016	A two-phase tank model for reducing parameter uncertainty.	The contributions of various drive mechanisms during flow back (fracture closure, gas expansion, and water depletion) are investigated.	Well length, total injected volume et al.	Fracture pore volume (PV), half-length, and permeability.
Xu et al. [15]	2017	An open-tank model	A comparative analysis to investigate the time variation of effective fracture volume during water flow back was conducted by this model.	Initial water saturation, rock compressibility et al.	Fracture–matrix interface area.
Hussain et al. [6]	2017	A fully 3-D simulation approach to estimate the SRV.	The real-time changes in the reservoir’s geomechanics as a function of fluid pressures are considered.	Fluid density, number of perforations et al.	SRV.
Qu et al. [5]	2019	Predict production model.	An integrated multiphysical model was developed to optimize the ESRV and predict the production.	Matrix permeability, bottom hole temperature, etc.	Production.
Xu et al. [14]	2020	Fracture compressibility model	The effects of rock and proppant parameters on fracture compressibility are investigated.	Fracture conductivity, etc.	The compressibility of fracture networks.
Umar et al. [4]	2021	A simple analytical model to compare the pressure/rate transient behavior of the three flow back cases.	The physics of flow back are understood by constructing basic diagnostic plots using two-phase flow back data from three multi-fractured horizontal wells.	Early-time flow back pressure and rate data, etc.	Single phase water rate and pressure.
Lan Ren et al. [3]	2023	Load recovery mode	A genetic expression programming is established to calculate flow back.	Geological engineering parameters.	Load recovery equation.

Table 2

Model calculation parameters.

Parameter	Unit	Value	Parameter	Unit	Value
P mi	106 Pa	50	T i	K	365.15
P fi	106 Pa	60	P BT	106 Pa	50
I imb	Dimensionless	0, 0.1, 0.3	C f	10−9 Pa−1	20
C m	10−9 Pa−1	0.1	C w	10−9 Pa−1	0.4
μ w	10−3 Pa·s	0.88	μ g	10−3 Pa·s	0.02
T pc	K	201	P pc	106 Pa	4.5
P L	106 Pa	2.95	V L	sm3/m3	0.3
L w	m	1200	h f	m	50
n f	Dimensionless	20	n CL	Dimensionless	3
ϕ m	Dimensionless	0.04	K m	10−9 μm2	200
B w	Dimensionless	1.05	α mf	m−2	200
W f0	m	0.009	l 0	m	30
R W	Dimensionless	0.7	R h	Dimensionless	0.7
R L	Dimensionless	0.7	θ	rad	π/6
m	Dimensionless	20	a	s	100,000

Table 3

Parameter settings of genetic algorithm for fitting flow back model parameters.

Genetic Algorithm Parameters		Value
Initial parameters	Population size	500
	Maximum algebra	1200
	Population initialization method	Random method
	Fitness function	Equation (72)
	Genetic selection	Roulette
	Termination condition	Maximum algebra
Genetic operator	Crossover probability	85%
	Mutation probability	15%

Table 4

Upper and lower limits of fitting parameters for flow back data inversion.

Parameter	Unit	Lower	Upper
l 0	m	1	300
W f0	m	0.000001	0.02
m	Dimensionless	1	500
R W	Dimensionless	0.00001	1
R h	Dimensionless	0.00001	1
R L	Dimensionless	0.00001	1
θ	rad	0	π/2
p BT	MPa	0	p mi
p fi	MPa	0	2pmi
α mf	m−2	0	1000
C f	10−9 Pa−1	0.1	1000
I imb	Dimensionless	0	1

Table 5

Basic parameters of shale gas well H3.

Parameter	Unit	Value	Parameter	Unit	Value
P mi	106 Pa	58.79	T i	K	367.8
TIV	m3	51,062.1	d casing	m	0.1143
L w	m	1500	h f0	m	43
C m	10−9 Pa−1	0.3	C w	10−9 Pa−1	0.46
μ w	10−3 Pa·s	0.28	μ g	10−3 Pa·s	0.042
T pc	K	190	P pc	106 Pa	4.61
P L	106 Pa	2.98	V L	sm3/m3	0.29
n f	Dimensionless	28	n CL	Dimensionless	3
ϕ m	Dimensionless	0.054	K m	10−9 μm2	380
B w	Dimensionless	1.05	n	Dimensionless	2

Table 6

Inversion results of genetic algorithm workflow.

Parameter	Unit	Value	Parameter	Unit	Value
l 0	m	50.9	W f0	m	0.0154
m	Dimensionless	24	R L	Dimensionless	0.618
R W	Dimensionless	0.614	R h	Dimensionless	0.616
θ	rad	1.07	P BT	106 Pa	58.79
P fi	106 Pa	58.81	α mf	m−2	575
C f	10−9 Pa−1	179	I imb	Dimensionless	0
EFNV	104 m3	1172	x f	m	90.85

Table 7

Results comparison of different models for estimating the initial EIFV.

Model	Unit	Value
HD model	m3	21,192
Alkouh model	m3	15,895
Established model	m3	10,643

References

1. Zhao, J.; Ren, L.; Jiang, T.; Hu, D.; Wu, L.; Wu, J.; Yin, C.; Li, Y.; Hu, Y.; Lin, R. et al. Ten years of gas shale fracturing in China: Review and prospect. Nat. Gas Ind. B; 2022; 9, pp. 158-175. [DOI: https://dx.doi.org/10.1016/j.ngib.2022.03.002]

2. Alkouh, A.; McKetta, S.; Wattenbarger, R.A. Estimation of effective-fracture volume using water-flowback and production data for shale-gas wells. J. Can. Pet. Technol.; 2014; 53, pp. 290-303. [DOI: https://dx.doi.org/10.2118/166279-PA]

3. Ren, L.; Wang, Z.; Zhao, J.; Lin, R.; Wu, J.; Song, Y.; Tang, D. Shale gas effective fracture network volume prediction and analysis based on flow back data: A case study of southern Sichuan Basin shale. Geoenergy Sci. Eng.; 2023; 228, 211963. [DOI: https://dx.doi.org/10.1016/j.geoen.2023.211963]

4. Umar, I.A.; Negash, B.M.; Quainoo, A.K.; Ayoub, M.A. An outlook into recent advances on estimation of effective stimulated reservoir volume. J. Nat. Gas Sci. Eng.; 2021; 88, 103822. [DOI: https://dx.doi.org/10.1016/j.jngse.2021.103822]

5. Qu, H.; Zhou, F.; Peng, Y.; Pan, Z. ESRV and production optimization for the naturally fractured keshen tight gas reservoir. Proceedings of the International Field Exploration and Development Conference 2017; Springer: Singapore, 2019; pp. 1576-1595.

6. Hussain, M.; Saad, B.; Negara, A.; Sun, S. Understanding the True Stimulated Reservoir Volume in Shale Reservoirs. Proceedings of the SPE Kingdom of Saudi Arabia Annual Technical Symposium and Exhibition; Dammam, Saudi Arabia, 24–27 April 2017; OnePetro: Richardson, TX, USA, 2017.

7. Adefidipe, O.A.; Xu, Y.; Dehghanpour, H.; Virues, C.J. Estimating effective fracture volume from early-time production data: A material balance approach. Proceedings of the SPE/CSUR Unconventional Resources Conference–Canada; Calgary, AB, Canada, 30 September–2 October 2014; OnePetro: Richardson, TX, USA, 2014.

8. Crafton, J.W. Oil and gas well evaluation using the reciprocal productivity index method. Proceedings of the Oklahoma City Oil and Gas Symposium/Production and Operations Symposium; Oklahoma City, Oklahoma, 9–11 March 1997; SPE: San Antonio, TX, USA, 1997; 37409.

9. Abbasi, M.A.; Ezulike, D.O.; Dehghanpour, H.; Hawkes, R.V. A comparative study of flowback rate and pressure transient behavior in multifractured horizontal wells completed in tight gas and oil reservoirs. J. Nat. Gas Sci. Eng.; 2014; 17, pp. 82-93. [DOI: https://dx.doi.org/10.1016/j.jngse.2013.12.007]

10. Williams-kovacs, J.; Clarkson, C.R. Stochastic modeling of two-phase flowback of multi-fractured horizontal wells to estimate hydraulic fracture properties and forecast production. Proceedings of the SPE Unconventional Resources Conference/Gas Technology Symposium; The Woodlands, TX, USA, 10–12 April 2013; SPE: San Antonio, TX, USA, 2013; 164550.

11. Clarkson, C.; Qanbari, F.; Williams-Kovacs, J. Innovative use of rate-transient analysis methods to obtain hydraulic-fracture properties for low-permeability reservoirs exhibiting multiphase flow. Lead. Edge; 2014; 33, pp. 1108-1122. [DOI: https://dx.doi.org/10.1190/tle33101108.1]

12. Ezulike, D.O.; Dehghanpour, H.; Virues, C.J.; Hawkes, R.V.; Jones, R.S., Jr. Flowback fracture closure: A key factor for estimating effective pore volume. SPE Reserv. Eval. Eng.; 2016; 19, pp. 567-582. [DOI: https://dx.doi.org/10.2118/175143-PA]

13. Fu, Y.; Ezulike, D.O.; Dehghanpour, H.; Steven Jones, R. Estimating effective fracture pore-volume from early single-phase flowback data and relating it to fracture design parameters. Proceedings of the SPE/CSUR Unconventional Resources Conference; Calgary, AB, Canada, 20–22 October 2015; OnePetro: Richardson, TX, USA, 2015.

14. Xu, Y.; Ezulike, O.; Dehghanpour, H. Estimating compressibility of complex fracture networks in unconventional reservoirs. Int. J. Rock Mech. Min. Sci.; 2020; 127, 104186. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2019.104186]

15. Xu, Y.; Dehghanpour, H.; Ezulike, O.; Virues, C. Effectiveness and time variation of induced fracture volume: Lessons from water flowback analysis. Fuel; 2017; 210, pp. 844-858. [DOI: https://dx.doi.org/10.1016/j.fuel.2017.08.027]

16. Jia, C.; Huang, T.; Yao, J.; Xing, H.; Zhang, H. Effect of isolated fracture on the carbonate acidizing process. Front. Earth Sci.; 2021; 9, 698086. [DOI: https://dx.doi.org/10.3389/feart.2021.698086]

17. Yao, J.; Ding, Y.; Sun, H.; Fan, D.; Wang, M.; Jia, C. Productivity Analysis of Fractured Horizontal Wells in Tight Gas Reservoirs Using a Gas–Water Two-Phase Flow Model with Consideration of a Threshold Pressure Gradient. Energy Fuels; 2023; 37, pp. 8190-8198. [DOI: https://dx.doi.org/10.1021/acs.energyfuels.3c00582]

18. Xu, P.; Yu, B.; Feng, Y.; Liu, Y. Analysis of permeability for the fractal-like tree network by parallel and series models. Phys. A Stat. Mech. Its Appl.; 2006; 369, pp. 884-894. [DOI: https://dx.doi.org/10.1016/j.physa.2006.03.023]

19. Fu, Y.; Dehghanpour, H.; Ezulike, D.O.; Jones, R.S., Jr. Estimating effective fracture pore volume from flowback data and evaluating its relationship to design parameters of multistage-fracture completion. SPE Prod. Oper.; 2017; 32, pp. 423-439. [DOI: https://dx.doi.org/10.2118/175892-PA]

20. Moussa, T.; Dehghanpour, H.; Fu, Y.; Ezulike, O. Dynamic fracture volume estimation using flowback data analysis and its correlation to completion-design parameters. Proceedings of the SPE Hydraulic Fracturing Technology Conference and Exhibition; The Woodlands, TX, USA, 5–7 February 2019; SPE: San Antonio, TX, USA, 2019; 194322.

21. Ren, L.; Lin, R.; Zhao, J.; Wu, L. Cluster spacing optimal design for staged fracturing in horizontal shale gas wells based on optimal SRV. Nat. Gas Ind.; 2017; 37, pp. 69-79. [DOI: https://dx.doi.org/10.1016/j.ngib.2017.10.001]

22. Lin, R. Dynamic Simulation of Stimulated Reservoir Volume (SRV) during Hydraulic Fracturing in Horizontal Shale Gas Well; Southwest Petroleum University: Chengdu, China, 2018.

23. Ezulike, O.D.; Dehghanpour, H. Modelling flowback as a transient two-phase depletion process. J. Nat. Gas Sci. Eng.; 2014; 19, pp. 258-278. [DOI: https://dx.doi.org/10.1016/j.jngse.2014.05.004]

24. Wang, F.; Li, B.; Zhang, Y.; Zhang, S. Coupled thermo-hydro-mechanical-chemical modeling of water leak-off process during hydraulic fracturing in shale gas reservoirs. Energies; 2017; 10, 1960. [DOI: https://dx.doi.org/10.3390/en10121960]

25. Bian, X.; Zhang, S.; Zhang, J.; Wang, F. A new method to optimize the fracture geometry of a frac-packed well in unconsolidated sandstone heavy oil reservoirs. Sci. China Technol. Sci.; 2012; 55, pp. 1725-1731. [DOI: https://dx.doi.org/10.1007/s11431-012-4775-z]

26. Song, B.; Ehlig-Economides, C.A. Rate-normalized pressure analysis for determination of shale gas well performance. Proceedings of the SPE Unconventional Resources Conference/Gas Technology Symposium; The Woodlands, TX, USA, 14–16 June 2011; SPE: San Antonio, TX, USA, 2011; 144031.

27. Abbasi, M.A. A Comparative Study of Flowback Rate and Pressure Transient Behaviour in Multifractured Horizontal Wells; University of Alberta: Edmonton, AB, Canada, 2013.

28. Williams-Kovacs, J.D. Quantitative Analysis of Multi-Phase Flowback from Multi-Fractured Horizontal Wells. Doctoral Thesis; University of Calgary: Calgary, AB, Canada, 2017.

29. Wylie, E.B.; Streeter, V.L. Fluid Transients; McGraw-Hill International Book Co.: New York, NY, USA, 1978.

30. Roy, P.P.; Roy, K. On some aspects of variable selection for partial least squares regression models. QSAR Comb. Sci.; 2008; 27, pp. 302-313. [DOI: https://dx.doi.org/10.1002/qsar.200710043]

31. Holland, J.H. Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence; MIT Press: Cambridge, MA, USA, 1992.

32. Beggs, D.H.; Brill, J.P. A study of two-phase flow in inclined pipes. J. Pet. Technol.; 1973; 25, pp. 607-617. [DOI: https://dx.doi.org/10.2118/4007-PA]

33. Lee, W.J.; Wattenbarger, R.A. Gas Reservoir Engineering; Henry, L. Doherty Memorial Fund of AIME Society of Petroleum Engineers: San Antonio, TX, USA, 1996.

34. Singh, H. A critical review of water uptake by shales. J. Nat. Gas Sci. Eng.; 2016; 34, pp. 751-766. [DOI: https://dx.doi.org/10.1016/j.jngse.2016.07.003]