**1. Introduction**

Multistage fracturing of horizontal wells is widely used in the exploration and development of shale gas. It typically takes tens of thousands of square meters of fracturing fluids and thousands of tons of sand to fracture a shale gas horizontal well, and the flowback of the fracturing fluids directly affects shale gas production. Hence, the flowback curve of fracturing fluids is a key issue in shale gas exploration [1]. The pores and fractures in shale gas reservoirs have different scales [2], and new multiscale pore-split systems are created by multistage fracturing. Thus, the shale gas flow mechanism is very complex, including pressure sensitivity, adsorption, diffusion, slip, imbibition, and seepage [3–6]. Currently, no consensus has been reached on this flow mechanism. Moreover, there is gas–water two-phase flow in the formation after large-scale fracturing. At present, there is no commercial software to realize the numerical simulation of flowback curve. It may take

**Citation:** Guo, W.; Zhang, X.; Kang, L.; Gao, J.; Liu, Y. Investigation of Flowback Behaviours in Hydraulically Fractured Shale Gas Well Based on Physical Driven Method. *Energies* **2022**, *15*, 325. https://doi.org/ 10.3390/en15010325

Academic Editor: Maxim Tyulenev

Received: 10 November 2021 Accepted: 23 December 2021 Published: 4 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

several years to develop such software because, when developing such calculation programs, the solutions of the flow equations, unstructured meshing, equation discretization, and large-scale irregular sparse matrix are very complex. The gas flow equation for shale gas involves many flow mechanisms. The solution of this equations requires many known parameters which are difficult to obtain directly or measure experimentally. For example, shale gas is a multicomponent mixture containing hydrocarbons and nonhydrocarbon gases (CO2 and N2) [7], and its molecular free path at a high temperature and pressure cannot be obtained for unknown component proportions. Even if there is a perfect numerical simulation program, inaccurate parameters will affect the reliability of the results.

In recent years, with the development of big-data-related technology, researchers began to introduce machine learning into the research of oil and gas development and have mainly used it as a prediction tool. Kohli et al. [8] took well-log data as input parameters and trained the multilayer forward neural network by using the least square Levenberg–Marquardt optimization algorithm to predict the formation permeability; the predicted permeability is consistent with field data. Jia et al. [9] studied the water channel problem resulting from long-term water injection. A density peak clustering algorithm based on streamline clustering was used to quantify the flow area for water flooding in the oil reservoir and thereby effectively identify the invalid water injection circulation channels between the injection and production wells, as well as areas with development potential. Adibifard et al. [10] carried out a Chebyshev polynomial interpolation of pressure derivative data and input the results to an artificial neural network to estimate reservoir parameters. Ghaffarian et al. [11] processed the pseudo pressure of gas wells and used the pseudo pressure derivative data as the input of single and coupled multilayer perceptron network to identify the condensate gas reservoir model. Tian et al. [12] used recurrent neural network learning to train the data collected by permanent downhole pressure gauge (PDG) for the inversion of reservoir permeability and other parameters and production prediction. Hung et al. [13] introduced the application of Gaussian process regression (GPR), support vector machine (SVM), and random forest (RF) to predict CO2 trapping efficiency in saline formations.

As the fracturing fluids flowback is influenced by many factors, big data technology has attracted increasing attention. Bai et al. [14] used exponential and harmonic functions to fit the water production in the flowback stage and produced the water stages of 32 wells and developed a prediction tool for accumulated water production based on the fitting results. Zhou et al. [15] used various binary and multivariate methods to analyse the data from 187 wells in the Marcellus block in the USA. The flowback rate during the first three weeks of data collection was found to increase with the thermal maturity and decrease with increasing strata thickness. Lin et al. [16] used a back propagation neural network to model the static data and flowback rate for 74 wells in Sichuan Province, China. The flowback rate and gas production in the first month were estimated using six parameters with relatively high weights in the model. Liu et al. [17] developed two neural networks with different structures to predict the flowback rate at specific time points and compared the predictions. In summary, current big data-based methods used to determine the flowback pattern focused on static or dynamic data of many wells for a given block, and a combined analysis of static and dynamic data has not been conducted. Some studies have been performed by using a fixed data point or a periodic average of dynamic data, which is a static data methodology in essence. All these lead to incomplete data utilization and limited application scope of the method.

In order to simplify the time-dependent flowback curve into simple parameters, as the target parameters of big data analysis and flowback influencing factor analysis, this paper combined the flowback dynamic data and static data of shale gas fracturing horizontal wells to study the flowback curve. The approaches and processes used in this study are shown in Figure 1. The first part of this paper briefly introduces the geological background of the research area and the relevant data collected. In the second part, the multi-stage fracturing physical model of shale gas horizontal wells was established. The flow equation in Laplace

space was obtained from the convolution formula, and the time-dependent flowback rate was obtained. In the third part, the flowback data of 214 horizontal wells in the Weiyuan block were fitted to obtain the characteristic parameters. The calculation formula of characteristic coefficient including characteristic parameters was established. Then, the correlation analysis was carried out by using a flowback characteristic coefficient, and the influencing factors of the flowback characteristic coefficient were comprehensively studied. Through the method established in this paper, the shale gas flowback curve containing complex flow mechanism can be abstracted into simple characteristic parameters and characteristic coefficients, and the relationship between static data and dynamic data is established, which can help to establish a machine learning method for predicting the flowback curve of shale gas horizontal wells.

**Figure 1.** Schematic diagram of the approach used in the study 2. Research area background and analysis data.

Weiyuan shale gas block is located in the southwest of Sichuan Basin of China. The whole Weiyuan block has a monoclinic structure inclined to the southeast. The Longmaxi Formation (LF) is the main exploration target layer. The LF can be divided into two layers: L1 and L2. The first layer of L1 can be subdivided into two sublayers, namely, L11 and L12. The L11 sublayer can be further subdivided into four layers, namely, L111, L112, L11<sup>3</sup> and L114, from bottom to top. The middle and lower parts of L11<sup>1</sup> have the best reservoir quality.

Table 1 shows the data from the study area. A total of 20 different types of data were collected from 282 wells in Weiyuan shale gas block, including the fracturing lateral length, the number of fracturing stages, first-year average daily production rate, vertical depth, TOC, porosity, high-quality reservoir thickness, gas saturation, pressure coefficient, brittle mineral content, average fracturing stage interval, fracturing fluid intensity, proppant intensity, average hydraulic fracturing fluid displacement, drilling length in high-quality reservoir, EUR (Estimated Ultimate Recovery), 30 day flowback rate, 90 day flowback rate, 180 day flowback rate, 360 day flowback rate, and peak gas production flowback rate. The data volume ranges from 214 to 282. The minimum data volume is 214, and the corresponding data of these 214 wells is used in subsequent analysis. The 30 day flowback rate, 90 day flowback rate, 180 day flowback rate, 360 day flowback rate, peak gas production flowback rate, and first-year average daily production rate are taken from the daily production reports of the horizontal wells. The fracturing lateral length, number of fracturing stages, average fracturing stage interval, fracturing fluid intensity, proppant intensity, and average hydraulic fracturing fluid displacement are taken from the drilling and completion reports of horizontal wells. The vertical depth, TOC, porosity, high-quality reservoir thickness, gas saturation, pressure coefficient, brittle mineral content and drilling length in high-quality reservoirs are taken from well logging interpretation.


**Table 1.** Statistical tables of data for the study area. Copyright permission: The copyright of Table 1 belongs to the publisher and authors.

#### **2. Basic Theory**

The time-dependent daily flowback volume curve for the Weiyuan shale gas block in southern Sichuan Province, China, shows that the daily flowback volume is large during the early stages of continuous production and decreases to a stable value during the late stage of production. The daily flowback volume during early-stage production is two to three orders of magnitude higher than that during the later stage, and the relationship between daily flowback and time is generally exponential. The flowback of a fracturing fluid after large-scale fracturing is essentially the seepage of fluid from the stimulated reservoir volume (SRV) area to the wellbore. The corresponding law can also be obtained from the seepage equation of fracturing fluid.

#### *2.1. Bottomhole Flow Equation*

The formula derivation in this paper considers single-phase liquid seepage, assumes constant flowing bottom hole pressure (FBHP) and continuous production, and is based on the multi-stage fracturing seepage equation of horizontal well [18]. The following dimensionless parameters are defined:

Dimensionless pressure:

$$P\_D = \frac{2\pi k h [p\_i - p(r, t)]}{qB\mu} \tag{1}$$

Dimensionless time:

$$t\_D = \frac{kt}{\phi \mu \mathbf{C}\_t L^2} \tag{2}$$

Dimensionless distance:

$$\mathbf{x}\_D = \frac{\mathbf{x}}{L} \tag{3}$$

$$y\_D = \frac{y}{L} \tag{4}$$

where *K* denotes the permeability, m2; *φ* denotes the porosity; *μ* denotes the viscosity, Pa.s; *Ct* denotes the comprehensive compressibility coefficient, 1/Pa; *xfi* denotes the half-length oftheithfracture,m;and*L*= *n* ∑*xfi*denotes thesumoverthehalf-lengthsofnfractures,m.

*i*=1 For multistage fractured horizontal wells, the total flowback volume is the sum of the flowback volumes of all the individual fractures. The flowback volume can be related to the FBHP in Laplace space as follows:

$$\sum\_{j=1}^{n} \overline{q}\_{Dj} = \frac{1}{s} \tag{5}$$

$$\overline{\mathcal{P}}\_{\text{WDi}}(f(\mu)) = \sum\_{j=1}^{n} s \overline{\eta}\_{Dj} \overline{\mathcal{P}}\_{Di^j}(f(\mu)) \tag{6}$$

where:

$$\overline{P}\_{Dij}(f(u)) = \pi \int\_0^\infty P\_D(t\_D) e^{-f(u)t\_D} dt\_D \tag{7}$$

$$P\_{\rm D}(t\_{\rm D}) = \int\_{0}^{t\_{\rm D}} \mathbf{G}\_{\rm xD}(\tau\_{\rm D}) \mathbf{G}\_{\rm yD}(\tau\_{\rm D}) d\tau\_{\rm D} \tag{8}$$

The expressions for *GxD* and *GyD* vary with the external boundary conditions. For enclosed strata with rectangular boundaries, *GxD* and *GyD* are defined as follows [19]:

$$G\_{\rm xD}(\tau\_{\rm D}) = \frac{1}{\chi\_{\rm eD}} \left[ 1 + 2 \sum\_{n=1}^{\infty} \exp\left( \frac{-n^2 \pi^2 (\tau\_{\rm D})}{2 \chi\_{\rm eD}^2} \right) \cos \frac{n \pi \chi\_{\rm wD}}{\chi\_{\rm eD}} \cos \frac{n \pi \chi\_{\rm D}}{\chi\_{\rm eD}} \right] \tag{9}$$

$$\mathcal{G}\_{yD}(\tau\_{\rm D}) = \frac{2}{y\_{\rm eD}} \left[ 1 + \frac{2y\_{\rm eD}}{\pi} \sum\_{n=1}^{\infty} \frac{1}{n} \exp\left(\frac{-n^2 \pi^2 (\tau\_{\rm D})}{2y\_{\rm eD}^2}\right) \sin\frac{n\pi}{y\_{\rm eD}} \cos\frac{n\pi y\_{\rm wD}}{y\_{\rm eD}} \cos\frac{n\pi y\_{\rm D}}{y\_{\rm eD}} \right] \tag{10}$$

Equations (5) and (6) can be used to obtain the following linear equation:

$$\begin{bmatrix} f(u)\overline{\mathcal{P}}\_{D1,1}f(u)\overline{\mathcal{P}}\_{D1,2}f(s)\overline{\mathcal{P}}\_{D1,3}\cdots f(u)\overline{\mathcal{P}}\_{D1,n} \\ f(u)\overline{\mathcal{P}}\_{D2,1}f(u)\overline{\mathcal{P}}\_{D2,2}f(s)\overline{\mathcal{P}}\_{D2,3}\cdots f(u)\overline{\mathcal{P}}\_{D2,n} \\ \vdots \\ f(u)\overline{\mathcal{P}}\_{Dk,1}f(u)\overline{\mathcal{P}}\_{Dk,2}f(u)\overline{\mathcal{P}}\_{Dk,3}\cdots f(u)\overline{\mathcal{P}}\_{Dk,n} \\ \vdots \\ f(u)\overline{\mathcal{P}}\_{Dn,1}f(u)\overline{\mathcal{P}}\_{Dn,2}f(u)\overline{\mathcal{P}}\_{Dn,3}\cdots f(u)\overline{\mathcal{P}}\_{D,n} \\ \end{bmatrix} \begin{bmatrix} \overline{q}\_{1D} \\ \overline{q}\_{2D} \\ \vdots \\ \overline{q}\_{kD} \\ \overline{q}\_{nD} \\ \overline{\mathcal{P}}\_{\overline{\mathcal{P}}\_{\overline{\mathcal{P}}\_{\overline{\mathcal{P}}}}} \\ \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \dots \\ 0 \\ \dots \\ 0 \\ 0 \\ 1 \end{bmatrix} \tag{11}$$

The solution of Equation (11) without considering the wellbore storage and the skin factor yields the FBHP as the dimensionless pressure *PWD*. The dimensionless pressure *PCWD* is obtained considering the skin factor and wellbore storage for multistage fracturing of horizontal wells in low permeability reservoirs as:

$$\overline{P}\_{\rm CWD} = \frac{f(u)\overline{P}\_{\rm WD}(f(u)) + S}{\mu \{1 + \mathcal{C}\_{\rm D} f(u) S \overline{P}\_{\rm WD}(f(u)) + S\}} \tag{12}$$

where *C* is the wellbore storage constant, Pa/m3; *CD* = *C* 2*πφCthL*<sup>2</sup> denotes the dimensionless wellbore storage constant; *S* is the total skin factor of the fracture system; *u* is the Laplace transform variable; and *f*(*u*) is a function that characterizes the properties of the strata. For a homogeneous stratum, *f*(*u*) *= u*.

For production wells with a variable flowback volume, the dimensionless pressure can be expressed in the following convolution form:

$$P\_D(\mathbf{x}\_D, y\_D, z\_D, t\_D) = \int\_0^{t\_D} Q\_D(t\_D - \tau) \frac{\partial P\_{\rm u}(\mathbf{x}\_D, y\_D, z\_D, \tau)}{\partial \tau} d\tau \tag{13}$$

where *Pu*(*xD*, *yD*, *zD*, *t*) denotes the dimensionless pressure of multistage fractured horizontal wells per unit flowback volume; *QD*(*tD*) = *q*(*t*) *qI* denotes the dimensionless flowback volume; and *qI* represents the unit flowback volume.

Using the properties of Laplace transform for convolutions yields the relationship between the pressure and flowback volume in Laplace space:

$$
\overline{\mathcal{P}}\_{\mathcal{D}}(\mathbf{x}\_{\mathcal{D}}, y\_{\mathcal{D}}, z\_{\mathcal{D}}, u) = \overline{\mathcal{Q}}\_{\mathcal{D}}(u) u \overline{\mathcal{P}}\_{u}(\mathbf{x}\_{\mathcal{D}}, y\_{\mathcal{D}}, z\_{\mathcal{D}}, u) \tag{14}
$$

Equation (14) can be used to obtain the flowback volume in Laplace space:

$$\overline{Q}\_{D}(u) = \frac{\overline{P}\_{D}(x\_{D}, y\_{D}, z\_{D}, u)}{\overline{\mu}\_{\text{\textquotedblleft}}(x\_{D}, y\_{D}, z\_{D}, u)}\tag{15}$$

Equation (15) can also be expressed in terms of the FBHP as follows:

$$\overline{\mathcal{Q}}\_D(u) = \frac{\overline{P}\_{wfD}(u)}{\overline{\mu}\_{wfu}(u)}\tag{16}$$

where *Pwfu*(*u*) is the *PWD* obtained by solving Equation (11) for a unit flowback using *GxD* and *GyD*, which are both exponential functions of time. Therefore, the flowback can be simplified to *Ql*(*t*) = *<sup>α</sup>exp*(−*βt*). Thus, the time-dependent flowback curve is defined by two characteristic parameters, *α* and *β*. Figure 2 shows the flowback volume versus time for horizontal segments of a 1000-m horizontal well with 20 stages, an effective fracture half-length of 40 m, and an effective fracture permeability of 0.25 mD. Fitting the curve yields the following equation for the flowback volume:

$$q(t) = 10\theta.23 \exp(-0.007t) \tag{17}$$

where *t* denotes the flowback time, (d); and *q*(*t*) is the daily flowback volume, (m3/d).

**Figure 2.** Theoretical flowback curve for multistage fracturing of horizontal wells.

Equation (9) accurately describes the characteristics of the flowback curve of horizontal wells, which indicates that α and β can be effectively used to parametrize the dynamic flowback curves. In addition, the flowback pattern is consistent with field measurements (Figure 3). The field data in Figure 3 are from one of the 282 Wells in Table 1. The abscissa is production days, the ordinate is daily fracturing fluid flowback volume, and the blue dot represents the daily flowback volume of a specific day.

**Figure 3.** Field flowback curve for multistage fracturing of horizontal wells.
