#### Abstract

For shale oil reservoirs, the horizontal well multistage fracturing technique is mostly used to reform the reservoir in order to achieve economic and effective development. The size of the reservoir reconstruction volume and the quantitative characterization of the fracture system are of great significance to accurately predict the productivity of shale oil wells. There are few flowback models for shale oil reservoirs. To solve this problem, first, a physical model of the simultaneous production of oil, gas, and water in the early flowback stage of shale oil development is established using the material balance equation for a fracture system. Second, the physical model of the underground fracture system is simplified, which is approximately regarded as a thin cylindrical body with a circular section. The flow of the fluid in the fracture system is approximately regarded as radial flow. In this model, the expansion of the fluid and the closure of the fracture are defined as integrated storage coefficients to characterize the storage capacity of the fracture system. Then, the curves illustrating the relationships between the oil-water ratio and the cumulative oil production and between the gas-water ratio and the cumulative gas production are drawn, and the curves are used to divide the flowback stage into an early stage and a late stage because the flowback process of shale oil wells exhibits obvious stage characteristics. Finally, the reservoir reconstruction volume and the related hydraulic fracture parameters are estimated based on the material balance method, and the rationality of the model is verified via numerical simulation. The interpretation results of this novel model are more accurate, making it an effective way to evaluate the hydraulic fracture parameters and transformation effect, and it has guiding significance for the evaluation of the hydraulic fracturing effect in the field.

#### 1. Introduction

In recent years, the exploration and development of unconventional oil and gas resources, such as shale oil, have increased, and unconventional oil and gas resources with rich reserves have been discovered all over the world. The broad prospects of the development of unconventional oil and gas resources have been demonstrated [1–5]. The properties of unconventional oil and gas reservoirs are significantly different from those of conventional oil and gas reservoirs. Unconventional reservoirs are generally characterized by low porosity and low permeability and the development of micro- and nanopores, and most of them do not have natural production capacity [6–8]. Currently, we mainly rely on multistage hydraulic fracturing of horizontal wells to form a complex fracture network in the reservoir and increase the effective drainage area in order to obtain industrial oil flow. Multistage fractured horizontal wells (MFHWs) have resulted in significant improvements in reservoir permeability and conductivity through horizontal well drilling and volumetric fracturing. During volumetric fracturing, hydraulically fractured fractures are formed perpendicular to the horizontal wellbore, and a complex fracture network is formed through the connection of the natural fractures inside the reservoir with the hydraulically fractured fractures [9–13]. The stimulated reservoir volume (SRV) is an important index used to evaluate the production capacity of MFHWs. Developing a method of effectively and correctly understanding the fracture network after volume fracturing has become a common hot topic of interest among scholars. Microseismic monitoring and postpressure evaluation techniques are currently used to interpret the parameters of MFHW fracturing modifications. However, the application of the above techniques is limited by cost factors and data interpretation techniques [14–18], and there has been a large gap between the actual application results and the awareness and needs of miners.

The drainage and production data in the initial stage of well opening after fracturing is the first-hand data that can best reflect the effect of fracturing. In recent years, these data have been deeply studied to evaluate the fracturing effect of MFHWs [19–25]. The fracture network formed by fracturing is initially filled with fracturing fluid. Therefore, analysis of the production data during early flowback should reveal the reservoir capacity and conductivity of the fractures. In the past, MFHWs were mostly used for the development of tight sandstone reservoirs. For tight sandstone reservoirs, there would be an obvious single-phase water flow stage during the initial flowback. It was found that by obtaining the flowback rate and flow pressure during the initial stage of flowback, combined with tracer data, the relevant information about the hydraulic fractures could be obtained from rate transient analysis of the flowback data. Abbasi et al. [26] proposed a single-phase rate transient analysis (RTA) model to analyze the water production data for hydraulic fractures and established a corresponding analytical mathematical model. Clarkson et al. [27, 28] divided the flowback stage of tight oil wells into preinvasion and postinvasion stages and established an analytical model to extract the fracture half-length and permeability from these two stages (permeability and fracture half-length).

Unlike tight sandstone reservoirs, shale gas wells are usually in the stage of rapid gas-water two-phase coproduction in the early stage of flowback. Postpressure shale gas wells are usually in a rapid gas and water production phase. Adefidipe et al. [29–32] divided the postpressure drainage of shale gas wells into an early gas production stage (EGP) and a late gas production stage (LGP) based on the turning point of the gas-to-liquid ratio curve. A volumetric analysis of the water and gas production data was performed by dividing the shale gas well production data phases, and a method of estimating the effective fracture volume by modeling the fracture system based on the two-phase material balance equation was proposed. The effective compression coefficient term and the linear dynamic relative permeability function were introduced to extend the existing single-phase analysis and two-phase simulation models to the analysis of shale gas reservoir flowback data. For shale oil reservoirs, Chen et al. [33–39] proposed a flowback model for complex fracture network shale oil reservoirs. During the numerical solution and inversion processes, it was necessary to use numerical solutions and semianalytical solutions to describe the flow of two-phase water and oil in a fracture and matrix system. In order to facilitate practical applications, in recent years, many researchers have used empirical decline curves to analyze multiflowback data. Jones and Blasingame [40, 41] applied hyperbolic and modified hyperbolic models to predict multiphase flow during reverse discharge. Fu et al. [42, 43] and Fu et al. [44] observed the harmonic decline behavior of the water return discharge and estimated the initial fracture volume versus final water production via decreasing curve analysis.

The goal of this study was to continue to develop an analytical model for the quantitative analysis of hydraulic fracturing flowback data. However, unlike in previous studies, the development of the mathematical model considers the transition from single-phase flow to multiphase flow when the formation fluid breaks through the fracture (three phases of oil, gas, and water). The expansion of the free gas in the fracture network, the expansion of the remaining fracturing fluid, and the closure effect of the fractures are considered in the model. These factors are considered as a comprehensive storage coefficient. Based on the single-phase diffusion equation, a complete mathematical model of shale oil flowback was established. Finally, the mathematical model was used to calculate the parameters related to the cracks in combination with the actual production data. The results obtained from the model were compared with the commercial software (Fracman) results, and the results were found to be relatively close. The calculated fracture data were also used for historical fitting of this field (using the commercial tNavigator software). The fitting results are good, which further demonstrates that the results of the model are reasonable. It is of great guiding significance for field production.

#### 2. Mathematical Model

##### 2.1. Mathematical Model of Material Balance

In order to relate the flow and pressure at the beginning of the flowback operation to the fracture and reservoir properties, in this study, a simple conceptual model was constructed for shale oil drainage return analysis (Figure 1). First, we described the conceptual model by constructing the material balance equation using the analytical equation. Next, we described the flow pattern of the fluid in the fractures at the beginning of the flowback by analyzing the production dynamic data. This led to obtaining the corresponding analytical solution by combining the continuity equation and the diffusion equation. Finally, the linear relationship between rate normalized pressure (RNP) and the matter equilibrium time (MBT) was established.

We made the following assumptions about this conceptual model. (1) The initial production during flowback occurs only in the fracture. (2) All of the fluid flows obey Darcy’s law. The effects of gravity and capillary pressure are negligible. (3) The contribution of the early production fluid from the matrix is negligible. (4) The fracture system contains an initial water (fracturing fluid) saturation of ; the initial free gas saturation is ; and the initial oil content saturation is . (5) The production from the horizontal wells is driven by fluid expansion within the fractures and closure of the fractures.

As can be seen from Figure 1, we approximate the entire crack system as a closed system. The equation describing the material balance in the crack obeys the law of conservation of mass. The fracture section is the SRV area formed by the hydraulic fracturing fractures, the secondary fractures generated by fracturing, and the natural fractures communicated after fracturing. The matrix section is a shale matrix connected to the fracture section. The wellbore section consists of a horizontal section and a vertical section. In this study, we defined a comprehensive compressibility factor, which takes into account gas expansion, oil expansion, fracturing fluid expansion, and fracture closure. The material balance equation is as follows:

, i.e., where is the volume factor, which is the ratio of the reservoir fluid volume under formation conditions to the fluid volume under surface conditions. is the volume of the gas phase in the fracture at any time when the fracture is filled with the oil, gas, and water phases. Similarly, is the volume of the oil phase in the fracture at any given time when the interior of the fracture is filled with the oil, gas, and water phases. This model assumes that the fluid flowing from the matrix into the fracture is negligible for a short period of time during the initial stage of the flowback. The variation in the fluid volume with time and the fluid density in the fractures and wellbore are described by the following equation:

For , solving Equations (1) and (2) simultaneously gives

For the convenience of calculation, Equation (3) is split into two parts:

According to the isothermal compressibility of the fluid,

Then, we can obtain the following formulas:

Substituting Equation (7) into Equations (4) and (5) and simplifying gives

Due to the fact that , , and are all functions of time , Equations (8) and (9) can be treated as follows:

Among them, , , and are considered to be functions of the cumulative oil production, cumulative gas production, and cumulative water production over time, respectively. , , and are defined as the original oil volume, original gas volume, and original water volume in the fracture, respectively. , , and are the oil, gas, and water outputs after blowout production, respectively.

Equations (12), (13), and (14) can be substituted into Equations (10) and (11) to obtain

Based on the equation for the compression coefficient of the fluid and using the chain derivation rule, we can rewrite Equations (15) and (16) as

By combining Equations (17), (18), (15), and (16), we obtain

The rate of change of the pressure in the wellbore with respect to time can be assumed to be approximately the same as the rate of change of the pressure in the fracture with respect to time during the early stage of flowback, and the volume of the fracture after fracturing is much larger than the volume of the wellbore. Therefore, can be neglected.

We can simplify Equation (19) to

The expression of the initial fracture volume is defined, which is a function of the initial oil saturation, initial gas saturation, and initial water saturation.

Both sides of the equation are normalized using Equation (21).

During the production process, the volume coefficient of the water hardly changes. By sorting out Equation (22), we can obtain

The total fluid produced in the initial stage of flowback is defined as . In addition, a comprehensive storage coefficient function is defined. Using these two newly defined functions, Equation (23) can be simplified as follows:

Each item in describes a different driving mechanism within the fracture system during the initial stage of flowback. It includes the effects caused by the fluid expansion and the crack closure with time.

##### 2.2. Model of Fracturing Fluid Flow in Fracture

Based on the material balance equation and diffusion equation, the mathematical model of the initial stage of flowback was established. There are many fracture forms after fracturing. When the crack height remains constant and the crack length increases, the longitudinal section of the crack is elliptical, and the crack shape is closer to the Perkins-Kern-Nordgren (PKN) model. In this study, we approximated the fracture surface sheet as a thin layer of cylinders, and the longitudinal section was approximately regarded as a circle. The flow of the underground fluid from the fracture space into the horizontal wellbore can be simplified as radial flow (Figure 2).

**(a)**

**(b)**

In this model, we considered the temperature and viscosity to be fixed values, ignored the influence of gravity, and simplified the fluid flow as single-phase radial flow of a slightly compressible fluid.

As is shown in the figure, the diffusion coefficient equation of the fluid flowing radially to the horizontal wells through the hydraulic fractures in the initial stage of flowback is

Substituting Equation (24) into Equation (27) gives

By integrating on both sides of Equation (26) at the same time, we obtain

When and , then is

Equation (29) can be written as

Integrating both the left and right ends of Equation (31) with respect to gives

When , , and we assume that , can be calculated at this time as follows:

Thus, Equation (32) can be written as

Here, we use Equation (35) to combine the average reservoir pressure with the pressure in the fracture and the bottom hole flow pressure: where the volume of the cylinder can be represented by

Substituting Equation (36) into Equation (35) gives

Substituting Equation (34) into Equation (37) gives

According to the partial integral, the solution of Equation (38) can be obtained as follows:

Equation (39) is sorted to obtain the following formula: where order and .

In the early stage of flowback, fluid expansion and fracture closure are the main factors affecting the fluid flow. According to the law of conservation of matter and the comprehensive storage coefficient,

By substituting the form of Equation (41) into Equation (40) and dividing both sides of the equation by , we obtain the rate normalization equation: where is the rate normalized pressure (RNP), and is regarded as the material balance time (MBT). Thus, we can rewrite Equation (42) as

can be calculated as follows:

In the early stage of flowback, assuming that the fracture stiffness is a constant value, the closed volume of the fracture can be expressed as

Substituting Equation (45) into Equation (44) and sorting it gives

Before the well is opened for production, the initial water volume in the fracture system and the total initial fracture volume obey the following relationship:

Based on Equation (21), the material balance equation for calculating the volume of the fracture system in the fracturing reconstruction can be obtained:

Due to the small compression coefficient of the water phase, the contribution of the water phase expansion in the fracture is negligible compared to the repulsion energy provided by the fracture closure. Therefore, Equation (48) can be further simplified to obtain the curve of the relationship between the total recovery volume and the pressure drop. We conducted linear regression of this curve to obtain the slope and intercept of the line. According to the obtained linear intercept, we can calculate the initial volume of the reconstructed fracture system:

The total surface area of the fracture system of the fracture modification can be obtained from the slope of the regression line, and the fracture stiffness of the shale can be obtained experimentally.

Equation (43) combines the material balance equation and diffusion equation, ignores the factors with relatively small influences, and finds a linear equation that is easy to apply to the actual production on site. By describing the linear relationship between the RNP and MBT, the calculations can be performed based on the production data measured in the early stages of rewiring. We regard Equation (43) as a univariate linear equation, where can be regarded as the slope of the curve. We can calculate the total storage coefficient of the fluid by calculating the slope of the curve. The intercept can be calculated to characterize the crack half-length of the radial flow under the PKN model.

#### 3. Results and Discussion

##### 3.1. Flowback Data Extraction

In order to establish a suitable mathematical model of early flowback in shale oil wells, in this study, the flowback data of fractured horizontal wells in the shale oil reservoir in block X in the Changqing Oilfield were analyzed, and the diagnostic curves of the changes in the oil, gas, and water productions with time and the changes in the oil-water ratio and gas-water ratio with cumulative gas production were drawn. The production data for this block were screened and sorted out, and several representative wells were selected for data analysis. By analyzing the flowback data, we obtained the oil production per hour, water production per hour, gas production per hour, total liquid production in different time periods, and bottom hole flow pressure in different time periods. Using the above data and formulas, we analyzed two typical wells and obtained the relevant parameters of the fracturing transformation.

In 2012, Abbasi et al. [26] drew related diagnostic charts by analyzing the flowback data for the tight oil and gas wells in Cardium. These diagnostic charts describe the relationships between the gas and water productions with time, as well as the relationship between the gas-water flow ratio and the cumulative gas produced. In this study, we used these diagnostic diagrams to describe the relationship between the oil-water ratio and the cumulative oil production and the relationship between the gas-water ratio and the cumulative gas production. According to the production data for the selected wells, a diagnostic map was drawn to determine the initial flowback time.

We selected the production data for two wells for analysis. Figures 3(a) and 4(a) show the oil production rate, water production rate, and oil-water ratio (OWR) of well A; and Figures 3(b), 4(b), and 5 show the oil production rate, gas production rate, water production rate, gas-water ratio (GWR), and OWR of well B. In the early stage of flowback, as well A was opened for flowback, oil was produced simultaneously with the fracturing fluid, and when well B was opened for flowback, oil, gas, and fracturing fluid were produced simultaneously. Therefore, it is necessary to consider multiphase simultaneous production in the flowback model. However, we can divide the flowback data for the shale oil into different regions according to the relationships shown in Figures 4 and 5, including early water production (EOP), early gas production (EGP), late oil production (LOP), and late gas production (LGP).

**(a)**

**(b)**

**(a)**

**(b)**

As the diagnostic rate diagram shows the instantaneous oil production, it can be inferred that the two phases (oil-water) or three phases (oil-gas-water) coexisted in the fracture network in the reconstruction area after hydraulic fracturing. When the fracture has just formed, it can be considered that the fracture is completely occupied by fracturing fluid, or there is a small amount of oil and gas only in the part connecting to the natural fracture. The fracturing construction period is long, and well blocking measures are usually taken after fracturing. During shut-in, under the action of gravitational differentiation and imbibition, the fluid in the matrix flows into the fracture network, leading to the coexistence of two phases (oil-water) or three phases (oil-gas-water). In the EOP and EGP regions, the trend line of the scattered points exhibits a negative slope, indicating that the initial oil saturation and gas saturation in the fracture network decreased as the production and development progressed. In the LOP and LGP areas, after part of the injected fracturing fluid flowed back, the water production gradually decreased. At this time, the oil or gas in the matrix flowed into the fractures, and the oil and gas production were supplemented. The slope of the trend line of the scatter plot of the oil-water ratio versus the gas-water ratio becomes positive and passes through the early flowback stage.

##### 3.2. Flowback Data Analysis

The early flowback data for well A and well B were analyzed, and the relationship between the total production volume and the pressure drop was drawn. It can be seen from Figures 6(a) and 6(b) that there is a strong linear correlation between the total recovery volume and the pressure drop. The initial crack system volume can be calculated using Equation (49). It is generally accepted that the volume of the hydraulic fractures formed during the hydraulic fracturing does not exceed the Total Injected Volume (TIV). The minimum value of the initial oil saturation in the hydraulic fractures can be inferred from Equation (49). The initial volume of the fracture system (including the hydraulic fractures and fracture modification zones that communicate with the natural fractures, which can be approximated as the SRV) is calculated by substituting the derived minimum initial oil saturation into Equation (48) and based on the total recovery volume versus pressure drop curve. Since the initial gas production from well B was small, the relevant parameters for well A can be used to estimate the initial volume of the oil phase in the SRV. Through the final calculation, it can be estimated that based on the data presented in Table 1, the SRV of well A was m^{3}, and the area of the reformed fracture system was m^{2}, while the SRV of well B was m^{3}, and the area of the reformed fracture system was m^{2}.

**(a)**

**(b)**

The linear relationship between the RNP and MBT was found by combining the EOP and EGP phases of the two completed MHFWs in Figure 7. The relationship between the fracture permeability, porosity, and fracture half-length was calculated using Equation (48), and the corresponding fracture half-length was obtained by estimating the fracture permeability [45]. It was calculated using the data given in Table 1. The hydraulic fracture half-length of well A was about 151 m, and the hydraulic fracture half-length of well B was 172 m.

**(a)**

**(b)**

##### 3.3. Numerical Simulation and Example Application

In order to verify the effectiveness of the calculation method developed in this study, hydraulic fracturing simulations and historical fitting of production performance data were carried out using the actual geological model, fracturing construction parameters, and production performance data for wells A and B. Figure 8 shows the half-length results of the fracturing simulation of the fractures in the two wells. Using the actual geological model and in situ stress parameters provided on site, an in situ stress model was established using the FracMan software, and the actual fracturing construction parameters were used for the fracturing simulation. Table 2 presents the fracturing simulation parameters of the two wells. The average half crack length of the simulated crack in well A is 162.3 m. The average half crack length of the simulated crack in well B is 189.6 m. These values are close to the calculation results.

Using the tNavigator reservoir numerical simulation software, the actual model of the well area was established. The relevant parameters of the fracture half-length and reconstruction area were calculated using this method, the dataset for the model, and the historical fitting. The production time of the two wells was from April 2013 to May 2021. As shown in Figure 9, the fitting results of the daily liquid production and daily water production of the two wells are good, which further demonstrates that the fracture parameters calculated using this method can reflect the actual transformation degree of the reservoir fracturing fractures better and have application value.

**(a)**

**(b)**

#### 4. Conclusions

(1)The material balance equation for the fracture system was used to establish a physical model for the simultaneous production of oil, gas, and water in the early flowback stage of a shale oil development. The physical model of the underground fracture was simplified and was approximately regarded as a thin-layer cylindrical body with a circular section, and the fluid flow was regarded as radial flow. Then, the novel mathematical model of multiphase flow flowback in the early stage of shale oil development was obtained by comprehensively considering the expansion of the fluid and the closure of the fractures in the model(2)The flowback stage of the shale oil exhibited obvious phase characteristics. According to the OWR and GWR, the flowback stage can be divided into an early stage and a late stage. Based on the early flowback stage, the SRV material balance model was established, and the size of the SRV was calculated based on the curve of the relationship between the total recovery volume and the pressure drop. Then, based on the mathematical model of the multiphase flowback in the early stage of shale oil development, the equation for the correlation between the RNP and MBT was obtained, and the relationship between the fracture permeability and the fracture half-length was calculated from it(3)The novel mathematical model was compared with commercial software results, and the resulting errors were found to be small. The obtained data were applied to a field example of historical fitting, and the fitting accuracy was high. The proposed model provides a novel and convenient calculation method for understanding the characteristics of hydraulic fracture parameters in production sites, and it has guiding significance for the development of unconventional reservoirs

#### Nomenclature

: | Rate (m^{3}/d) |

: | Density (kg·m^{−3}) |

: | Volume (m^{3}) |

: | Volume of fluid in wellbore (m^{3}) |

: | Volume of oil in wellbore (m^{3}) |

: | Volume of gas in wellbore (m^{3}) |

: | Time (d) |

: | Fluid viscosity (mPa·s) |

: | Volume compressibility (m^{3}/m^{3}) |

: | Compressibility (1/MPa) |

: | Total compressibility (1/MPa) |

: | Total storage coefficient |

: | Initial free oil volume (m^{3}) |

: | Initial free gas volume (m^{3}) |

: | Initial free water volume (m^{3}) |

: | Drainage radius (m) |

: | Wellbore radius (m) |

: | Flowing bottom hole pressure (MPa) |

: | Pressure of fluid in wellbore (MPa) |

: | Average reservoir pressure (MPa) |

: | Initial reservoir pressure (MPa) |

: | Production of fracturing fluid (m^{3}) |

: | Horizontal well length (m) |

: | Fracture length (m) |

: | Fracture bulk porosity (%) |

: | Total surface area of modified crack system (m^{2}) |

: | Fracture stiffness of shale (MPa/m) |

: | Crack width (μm). |

*Subscripts*

: | Water |

: | Gas |

: | Fracture |

: | Matrix. |

#### Data Availability

The data used to support the findings of this study are intersection within the article.

#### Conflicts of Interest

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This work was supported by the Chinese National Natural Science Foundation (No. 51774256 and No. 51974356).