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:
where Δ
Pk is the pressure difference of the k-th level fracture, Pa;
μ is fluid viscosity, Pa s;
lk,
Wfk,
hfk 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).
where
l0,
Wf0, and
hf0 are the length, width, and height of the initial tree-shaped fractal fractures, m;
RL,
RW, and
Rh 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
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):
where
where
n is the total number of branches fractal fracture,
n = 2.
Then the flow rate is described by Equation (7):
where Δ
P is the total pressure difference of fractal fracture network, MPa.
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
Wfk at the k-th level has a width of
Wfkc after compression, then the volume change is shown in Equation (9):
where
Cf is fracture compression coefficient;
Vfk is the k-level single fracture volume under original fracture pressure;
Vfkc is the k-level single fracture volume under current fracture pressure; and Δ
Pf is fracture system pressure drop, Δ
Pf =
Pfi −
Pf.
where
Wfkc 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.
The flow rate of 1/2 single-cluster tree-shaped fractal fracture network is expressed in Equation (13) by considering the fracture closure effect.
where
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.
where
Pf is the average pressure in the fracture system, which changes as fluid is produced, MPa;
Pwf is the wellbore bottom flow pressure in a horizontal well, MPa;
I is water or gas;
Bi is volume coefficient, m
3/m
3; and
Kri(
Sw) is the relative permeability of gas/water, dimensionless.
In this paper, the linear relative permeability model was adopted:
where
Sw is the water saturation in the fracture.
Gas production superposition model is as follows:
Water production superposition model is as follows:
where
Nf is the total clusters in horizontal well, which satisfies the following relationship:
where
nf is the number of operation stages, dimensionless; and
nCL is the clusters number, dimensionless.
2.2. Calculation of Fracture Network Volume
The volume of the fractal fracture network is
Recently, the initial effective fracture volume
Vfi has been used internationally as an important parameter for preliminary evaluation of shale gas fracture network fracturing [
19,
20]. In the above formula,
where
Vfi is initial volume of effective fracture, m
3; and
V0 is initial fracture volume, m
3.
It can be seen from
Figure 1 that the longitudinal equivalent fracture half-length is
where
θ represents the bifurcation fractures, which varies from 0 to π/2.
As shown in
Figure 1, the lateral extension width is given by
The effective fracture network volume can be obtained by Equation (25):
where
xf denotes the longitudinal expansion degree and
wf denotes the lateral expansion degree. It is an important parameter commonly used in quantitative evaluation of shale gas fracture network fracturing [
21,
22], where
where
Vb_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 (
Iimb) 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
Iimb represents the proportion of free gas in the fracture network to the fracture volume, 0 ≤
Iimb ≤ 1. Obviously, there is the following relationship, the gas saturation of the original fracture is
The volume of free gas within the subsurface fractures is
Then the volume occupied by water within the fracture is
Water saturation of fractures is as follows:
where
Sgi represents the original gas saturation within fractures, dimensionless;
Vgfi represents the volume of free gas within the subsurface fractures, m
3;
Vwi is water volume in fractures, m
3; and
Swi 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 Pfi to the current fracture pressure Pf, and the fracture pressure drop is ΔPf = Pfi − Pf when a certain amount of fracturing fluid (Wp) 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
- (a)
In fracture network, the expansion of free gas can be expressed as
where
Vgf is the gas volume in the current fracture, m
3;
Vgfi is the gas volume in the original fracture, m
3;
Bgf is the gas volume factor in the fracture, m
3/m
3;
Bgfi is the original gas volume factor in the fracture, m
3/m
3;
Pfi is original fracture pressure, MPa; and
Pf 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.
where
Bgm and
Bgmi are matrix gas volume coefficient under the condition of current matrix pressure and original matrix pressure, m
3/m
3;
ϕm is matrix porosity of shale reservoir, dimensionless;
Cm is rock compressibility coefficient of shale matrix, 1/MPa;
VEi and
VE are unit shale adsorbed gas volume under original matrix pressure and current matrix pressure, m
3/m
3; and Δ
Pm is matrix pressure drop, MPa.
Vb 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:
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.
where
Pmi is shale original matrix pressure, MPa;
Pm is shale current matrix pressure, MPa;
VL is Langmuir volume, sm
3/m
3;
PL is Langmuir pressure, MPa; and
Gmf is matrix shale surface volume of intrusion of gas channeling into the tree-shaped fractal fracture network, m
3.
- (c)
Produced free gas (underground volume)
The gas storage capacity of fractures is obtained by combining Equations (33), (38), and (39).
where
Gp is cumulative gas production, m
3.
- (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
Pfi decreases to
Pf, the fracture fracturing fluid volume is
Then, Equation (42) can be obtained by combining Equations (32), (40), and (41):
The remaining fracturing fluid volume is converted to surface conditions:
The basic form of the fracturing fluid material balance equation in the fracture is as follows:
where
Bwi is the fracturing fluid volume coefficient under original pressure, dimensionless;
Wp is current cumulative fracturing fluid volume of flow back m
3; and
Wres is the remaining fracturing fluid volume, m
3.
The remaining fracturing fluid volume W
res is inputted into the above formula:
Then,
where
Bwf is the fracturing fluid volume coefficient under the current fracture pressure, m
3/m
3; and
Bwfi is the fracturing fluid volume coefficient under the original fracture pressure, m
3/m
3.
The compressibility coefficients of gas and water are
The function
h representing the fracture pressure and matrix pressure at the
k + 1th time step is defined as follows:
where
where
Cwf is fracturing fluid compression coefficient in fractures, dimensionless;
Psc is ground pressure, MPa;
Zsc is gas deviation factor in ground conditions, dimensionless;
Tsc is ground temperature, K;
Zf is gas deviation factor in fractures, dimensionless; and
Ti is original formation temperature of reservoir, K.
Water saturation of fractures under current formation conditions is as follows:
The equation can be simplified to
The equations for fracture and matrix pressure at time k + 1 can be derived by setting Equation (49) equal to zero.
The material balance equation needs to be established of the shale matrix system to solve this equation due to there being two unknowns, Pfk+1 and Pmk+1.
If the fracture pressure (
Pf) falls below the matrix gas breakthrough pressure (
PBT), it will result in channeling between the matrix and fractures. The pressure
PBT is related to the porosity and permeability of the matrix. The channeling equation [
24] is described by Equation (54).
If
PBT ≤
Pfk, then the channeling from the matrix to the fracture system from time
k to
k + 1 is 0:
If
PBT >
Pfk, the channeling from the matrix to the fracture system from time
k to
k + 1 is
where Δ
Gmf is the channeling of gas from the shale matrix to the fracture system from time
k to time
k + 1, m
3; Δ
t—the time from time
k to time
k + 1, s;
Vb is shale reservoir stimulation volume, m
3;
ϕm is matrix porosity of reservoir, dimensionless;
qm is gas channeling flow from matrix to fracture system per unit time, s
−1;
km is matrix permeability, um
2;
αmf is matrix to fracture channeling factor [
25], m
−2;
Pmk is average pressure of shale matrix system at time
k, MPa;
Pfk+1 is average pressure of fracture system at time
k + 1, MPa; and
ug 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
where
where
VL is Langmuir volume, sm
3/m
3;
PL—Langmuir pressure, MPa.
Bgm and
Bgmi are matrix gas volume coefficient under current matrix pressure and original matrix pressure, m
3/m
3;
Cm is rock compressibility coefficient of shale matrix, 1/MPa; and
VEi and
VE are the unit shale adsorbed gas volume under original matrix pressure and current matrix pressure, m
3/m
3.
Equations (58) and (59) are substituted into Equation (57), then Equation (61) can be obtained.
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):
The equations for the fracture network and matrix pressure at time
k + 1 can be derived by setting Equation (62) to zero:
There are two unknowns, Pfk+1 and Pmk+1, in above equation, and the fracture pressure Pfk+1 and matrix pressure Pmk+1 at time k + 1 can be obtained by solving Equations (53) and (63), and the fracture pressure Pfk and matrix pressure Pmk at time k are known.
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]:
where
Pfi is initial fracture network pressure, MPa;
Pwf is bottom hole flow pressure, MPa;
qw is water production rate, m
3/h;
tMB is equivalent time, h; and
Wp is cumulative water production, m
3.
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 (
Iimb = 0) and the observed flow stages are consistent with field data [
27] and analytical models [
26] and Williams-Kovacs [
28]. For the case of
Iimb = 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
Iimb = 0.1 is greater than the slope of 0.65 at
Iimb = 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.