Investigation of Flowback Behavior for Multi-Fractured Horizontal Wells in Gulong Shale Oil Reservoir Based on Numerical Simulation

1. Introduction

As human demand for energy has increased, conventional oil and gas resources have been unable to meet people's demand. Shale oil represents a type of unconventional energy resources with huge reserves, and the research and development of shale oil and gas have attracted much attention. In a broad sense, shale oil refers to oil and gas resources in organic-rich shale and adjacent non-source rock interlayers. Zou et al. [1] believe that shale oil constitutes mature oil contained within organic-rich shale formations characterized by nanoscale pore sizes. Organic-rich shale develops laminar structure, and micron-nano-sized pore pipes and micro-fractures are the main storage space. Jiang et al. [2] believe that shale oil refers to the liquid hydrocarbons in free, dissolved or adsorbed state in the effective source rock shale, which is the residual retention and accumulation after hydrocarbon generation and shale expulsion. Shale oil flows little or not far from the source rock. Shale oil is characterized by low porosity and low permeability, tight reservoirs and strong heterogeneity, which often result in difficult transformation of shale oil during the production process, low oil rate per well, rapid decline of oil rate and low or no natural production The production of shale oil depends on the combination of horizontal drilling technology and hydraulic fracturing technology. Horizontal drilling in shale oil reservoirs increases the contact area between the reservoir and the wellbore, thereby increasing productivity and recovery ratio. Hydraulic fracturing is a method of fracturing rock by injecting fluids at high pressure to increase the flow of reservoir fluid.

During fracturing, shale oil develops complex fracture networks. The hydraulic fractures formed during fracturing interlock with the original fractures in the shale, and the original fractures are extended to different degrees. Evaluating the conductivity of fractures is very important for the shale oil flowback process. Palmer and Mansoor [3] incorporated pore compressibility and permeability as functions of effective stress and matrix shrinkage into a single equation, resulting in a widely adopted theoretical permeability model. Schutjens, P.M.T.M. et al. [4] proposed the relationship between the degree of compaction, the decrease of porosity based on the linear pore elasticity theory of small strain, and proposed that decrease in permeability of sandstone near the elastic range is primarily dependent on the increase in average effective stress. Legrand, N. et al. [5] studied the production mechanism of tight oil in Yemen, and the crude oil flow mainly depended on matrix microfractures and fracture channels. The crude oil flows through the matrix microfractures into the fracture channel and then is produced. Lu [6] established an elastic-plastic deformation model of self-supporting fractures after hydraulic fracturing in deep shale, and studied the flow capacity and characteristics of self-supporting fractures.

A mathematical model is a method used to study the flowback rule after fracturing. Fang et al. [7] developed a shale gas productivity prediction model for multi-stage fractured horizontal wells, which considers the flow of fracture systems, shale matrix compressibility, gas slip, Knudsen diffusion and surface diffusion. The semi-analytic model is solved using Laplace integral transformation and the point source method. The research results of Du [8] provide a new perspective for understanding the flow law of shale gas reservoirs with dual media. His non-Darcy flow model takes into account a multitude of complex factors, including adsorption, desorption and slippage effects, and has a higher degree of accuracy and reliability than the traditional Darcy flow model. Furthermore, his model can be solved by Laplace spatial solutions under different boundary conditions, and has been verified by drawing well test curves. Andersen [9] developed a shale gas productivity prediction model that considers the compressibility of the shale matrix and solves it using a semi-analytical method. This enables the spatial distribution of pressure, adsorption, porosity and apparent permeability to be determined at a given time. Gao et al. [10] developed a mathematical model for the nonlinear flow of shale gas, taking into account the random distribution of non-uniform equivalent porosity and permeability. They proposed a semi-analytical method and an explicit iterative scheme based on the Boltzmann transform and local homogenization approximation. Ren et al. [11] constructed an unstable flow mathematical model of stagewise fractured horizontal wells in low permeability reservoirs. They solved this model by combining perturbation transform, Laplace transform, and superposition principle with a semi-analytical method. However, when the mathematical model is used to describe the fluid flow, the flow mechanism is limited. For the mathematical model of a simple mechanism, it can be solved quickly and has great application value. However, for a complex mechanism, it is often impossible to solve. When a discrete crack is used, the calculation is large and the time is long. When the equivalent fracture zoning model is employed, the fundamental characteristics of the fracture network cannot be reflected. Furthermore, the fracture distribution is challenging to characterize, the difference in liquid production between the fractures cannot be considered, the yield is evenly divided, and the irregular reconstruction area is difficult to address.

The fracture network structure of fractured reservoir is complex and the reservoir contains multi-phase fluids such as oil, water and gas. The mechanism of flowback stage after fracturing is very complicated. The numerical simulation method can accurately simulate the fracture distribution and multiphase fluid flow of stage-fractured horizontal Wells, which provides an effective way to study the flowback mechanism and production prediction after fracturing. Clarkson C. R. [12] established a shale gas production model considering adsorption and slippage effects, which is highly applicable to shale reservoirs without fractures. Takuto Sakai et al. [13] developed a new type of reservoir simulator to predict the flow behavior of shale gas/oil, which can simulate three-dimensional and three-phase flow behavior. The simulator is able to reproduce the flow of fluid between induced hydraulic fractures and natural fractures. Fan et al. [14] established a multistage fractured horizontal well composite model considering SRV based on the characteristics of low porosity and low permeability in tight oil gas reservoirs, and used the finite element method to solve the model. The fluid flow around the hydraulic fracture is divided into 5 flow states, and the larger the SRV, the shorter the time to reach quasi-steady flow. Zhao et al. [15] used CMG software to build a single well productivity model of shale oil fracturing. The fluid model adopted a component model to study and explore the effects of compaction, gas phase non-Darcy effect, and nanoporous fluid phase state on productivity. Through the inversion of reservoir parameters through historical fitting, the final recovery ratio is predicted. The model does not take into account the stress sensitivity of opened beddings and the threshold pressure gradient in the matrix.

A suitable flowback strategy can enhance the production of shale oil. The application of measures such as soaking and controlled pressure flowback has been demonstrated to be an effective means of increasing shale oil production in practice. In a study by LaFollette [16], multiple linear regression and an enhanced tree model were employed to analyze a range of reservoir quality indicators, well structure, completion, stimulation measures and production parameters in the Barnett shale gas and Middle Bakken oil areas. The study concluded that the quality of the reservoir is a principal determinant of tight oil production. Wang et al. [17] developed a dual-medium two-phase flow model of shale reservoirs, which considers the flow mechanism of shale reservoir fluids. The flow behavior throughout the entire cycle of drilling and production is simulated, and the impact of flowback performance, including soaking time and pressure drop rate, on shale oil production is studied. Zhu [18] conducted a comprehensive investigation into the nine factors that influence shale oil production capacity, utilizing layer system analysis and the entropy weight method. Furthermore, from the perspective of enhancing the energy efficiency of fracturing fluids, he put forth an efficacious system for shale oil production, comprising post-fractured soaking and the regulation of the decline in pressure during the flowback phase. Nevertheless, a transparent system for controlling the rate of pressure decrease and an optimization method for shale oil reservoirs have yet to be established. Karantinos et al. [19] proposed a coupling model between the wellbore and the reservoir. This model employs the black oil model in the reservoir and multiphase flow analysis in the wellbore. The model permits the choke size to be altered at different times, and it is possible to study the influence of a choke size management strategy on shale production. This approach is anticipated to assist operators in selecting the optimal choke size, enhancing well performance and efficiency, and reducing the risk of wellbore or completion failure. Bagci et al. [20] developed a methodology for selecting the optimal choke size following fracturing. The impact of optimized choke size and flowback time on proppant production was demonstrated. The production of proppant was found to decrease with an increase in choke size and flowback period. The results demonstrate that the optimal choke size is contingent upon pressure drop, fluid volume, wellhead pressure, and fracture geometry parameters. Liu [21] employed the EDFM to construct a discrete fracture model, considered the fracture network formed by real reservoir fracturing and the migration of proppant, proposed three choke systems, studied the changes in fracture conductivity and production under different choke systems, and concluded that radical choke systems lead to a rapid reduction in bottomhole flow pressure and proppant production, thereby reducing fracture conductivity and resulting in a lower recovery ratio. Nevertheless, this model does not sufficiently consider the law of fluid flow within the shale matrix. In general, a realistic choke management analysis is lacking in unconventional reservoirs.

In this paper, based on a case of shale oil reservoir in Gulong, Daqing, extensive research has been conducted on the flowback performance attributes and factors influencing shale oil reservoirs. The production system for shale oil is reasonably optimised. Firstly, the analysis focuses on the reservoir and fluid properties of the Gulong shale reservoir, and the effect of opened bedding by hydraulic fracturing on flowback is considered. A numerical model of multiphase flowback of the Gulong shale reservoir was established based on the flow mechanism, including PTPG and stress sensitivity effects. This model was utilized to study the influence of geological parameters, flowback modes, and other factors on the shale oil flowback performance. Secondly, the effects of different factors (PTPG, stress sensitivity, opened bedding permeability, matrix permeability) on the flowback performance of the Gulong shale reservoir were studied using the numerical model established earlier. Subsequently, a summary of production data from Daqing Gulong shale oil wells is presented, with nine representative wells selected for analysis. Furthermore, three modes of shale oil flowback are distinguished: the fast-slow mode, the slow-fast-slow mode, and the fast-slow-fast mode. A numerical model is employed to examine the impact of various flowback modes on the flowback performance of shale reservoirs. Finally, the mode was optimized from the perspective of total oil and recovery ratio. The study is of great significance for understanding the flowback performance and for the beneficial development of Gulong shale.

2. Characteristics of Gulong Block

The Qingshankou Formation within the Songliao Basin boasts abundant shale oil resources, with its geological resources amounting to 15.1 billion tons, which is the main direction for developing Daqing Oilfield [22,23]. Daqing Gulong shale oil distributes in the northern part of Songliao Basin, Qijia-Gulong depression, which is a typical continental deposit shale oil, covering an area of 2,778km2. The primary petroliferous strata of Daqing Gulong Shale oil is the Qingshankou strata, which developed an extensive depositional environment of semi-deep lakeand deep lake phase. Many dark mud shales rich in organics gathered in this area. The thickness of the sweet spot strata is 40-120 m [24], and the maturity of Gulong shale oil is about 1.20%-1.67%. Daqing Gulong shale oil characterizes by mature crude oil, highly developed beddings, and a high gas-oil ratio. In the continuous development of Gulong shale oil, it found that the complexity of flowback performance after Gulong shale oil soaking has become a significant cause restricting the development of Gulong shale oil.

2.1. Reservoir Characteristics of the Gulong Block

The microporous structure analysis of the Gulong shale reservoir core have shown that that Gulong shale oil mainly stored in pores of clay minerals and shale beddings. The total porosity of Gulong shale matrix is about 5%-12%. And the pore throat diameter ranged from 10-100 nm [25]. The matrix permeability of Gulong shale oil reservoir is about 5×10-4mD. And there are about 500-3000 beddings per meter [24]. Extensively developed shale beddings result in significantly higher horizontal direction permeability than vertical permeability.

The indoor physical simulation study of Gulong shale reservoir shows that, due to the development of beddings in shale, the fracturing fluid will flow into the shale beddings cemented weakly during the fracturing process, opening numerous horizontal opened beddings along the extension direction of the shale. Part of the beddings cemented weakly will widen. And at the same time, the vertical hydraulic fractures will connect with the opened beddings. The fractures show the structure of the joint network as shown in Figure 1. The existence of the fractures network increases the mobilizability of shale oil, but also leads to a complex flow state, including flow in porous media and pipe in the volume containing the fractures. The mechanism of the post-fractured flowback stage after soaking is understood poorly.

2.2. Fluid Characteristics of the Gulong Block

Measured Crude oil from well A was after degassing: the viscosity of the crude oil was 2.77-15.6 mPa∙s, the density of the crude oil was 0.79-0.82 g/cm3, and the production gas-oil ratio was 370 m3/m3. Pang et al [24] measured the PVT parameters of Gulong shale oil, and obtained that the bubble-point pressure of the crude oil was 13.31 MPa, the coefficient of compressibility was 1.3625 GPa-1. The viscosity, density, and gas-oil ratio of this crude oil under reservoir conditions were 0.80 mPa∙s, 0.724 g/cm3, and 76.46 m3/m3, respectively. Overall, Gulong shale oil characterized by low viscosity, low density, medium-high gas-oil ratio, and high reservoir pressure [24].

The escape of too much gas will cause the viscosity of oil to become large [25], which is not conducive to the oil flow. And the existence of the free gas in shale also increases the flow resistance of the oil phase. Therefore, the impacts of gas escape and gas flow must be considered, in the study of Gulong shale oil flowback performance.

3. Model Descriptions

3.1. Model Assumption

After fracturing, the hydraulic fracture permeability is approximately 10D. The permeability of the opened bedding is approximately 10mD. The permeability of the matrix is approximately 5.8×10-4mD [26]. The order of permeability magnitudes leads to the conclusion that the flow capacity of the hydraulic fractures and opened beddings is significantly greater than that of the matrix. As the well begins to produce, the fluids inside the hydraulic fractures and beddings are rapidly recovered, resulting in a rapid decrease in pressure within the hydraulic fractures and bedding. The fluids in the reservoir rock matrix will initially flow into the opened beddings and subsequently into the hydraulic fractures.

A physical model of matrix-bedding-hydraulic fractures is established, as shown in Figure 2. The matrix region between two adjacent opened beddings is taken as the study area. Due to the proximity of the two opened beddings and the weak influence of gravity, the fluid in the matrix in this area flows linearly to both sides of the opened beddings, which can be considered as symmetric. In order to ensure the calculation efficiency for productivity evaluation and interpretation, only half of this region can be studied. The middle of the matrix between two opened beddings is considered as a closed boundary, which means that the difference between multi-opened beddings is ignored. Similarly, the opened beddings are divided into two parts, with the assumption that the fluid flows from the matrix into the opened beddings and flows on the side near the matrix. This allows for the study of half of the opened beddings. Under this physical model, the yield of a single cluster is the model capacity multiplied by the number of beddings multiplied by two. The internal flow pattern of the shale oil can be described as follows: the oil initially flows into the matrix, then into the opened bedding, and subsequently into the hydraulic fractures, before finally entering the horizontal wellbore.

In constructing the numerical model, the following assumptions are made:

(1) Given the dominant role of the main fractures in the fracture system, the complex fracture network formed by hydraulic fractures and natural fractures is simplified. Fluid production is assumed to be uniform along the wellbore and vertically.

(2) The hydraulic fracture is considered to have infinite conductivity, with the flow process in the artificial fractures to the wellbore being ignored. Only the flow process in the matrix and opened beddings is considered.

(3) Opened bedding penetrates the reservoir in the fracture length direction.

(4) It is assumed that the original structure of the shale reservoir is not damaged.

(5) Following fracturing, a water-encroaching layer is observed in the matrix in the vicinity of the opened bedding, with a water saturation of 100%.

(6) In the absence of a soak or a brief soak period following fracturing, the pressure in the seam is comparable to the fracture pressure, with a water saturation of 100% in the opened bedding.

3.2. Reservoir Model

In the model established, the artificial hydraulic fractures can be regarded as having an infinite conductivity. The matrix and the opened bedding are the main factors affecting the fluid recovery. Therefore, according to the fluid flow mode assumed in [27], the matrix and the opened bedding are established to simulate the fluid flow from the matrix to the opened bedding and the fluid flow from the opened bedding to the hydraulic fractures. In the matrix, the flow is approximated to linear flow. A rectangular coordinates system is used for the purposes of meshing.

Niu et al. [28] determined that the degree of opened bedding development was 0.2-0.36 per meter in the Gulong block following fracturing by analyzing the flowback of fractured horizontal wells. As shown in Figure 3, the model comprises 51 grids, with one 100m×0.01m×30m opened bedding grid situated at the top and fifty 100m×0.1m×30m matrix grids at the bottom. Within the opened bedding, a larger stress sensitivity coefficient was assigned to simulate the decrease in permeability resulting from the closure of the opened bedding with decreasing pressure. A smaller stress sensitivity coefficient was assigned to the matrix in order to simulate the reduction of microfractures in matrix permeability. Furthermore, a specific PTPG was set for the shale matrix.

3.3. Fluid Parameters

Gulong shale oil represents a type of unconventional oil and gas resource that is distinguished by its high maturity, elevated content of light hydrocarbons, and a notably high gas-to-oil production ratio. The Figure 4 illustrates the relationship between pressure and the oil viscosity, oil volume coefficient, and dissolved gas-oil ratio of Gulong shale oil. The saturation pressure of Gulong shale oil is 31.3MPa. When the pressure exceeds the saturation pressure, there is no gas escaping as the pressure decreases, resulting to a gradual increase in oil volume coefficient and a decrease in the oil viscosity. Conversely, when the pressure is below the saturation pressure, a significant increase in oil viscosity occurs as the dissolved gas escapes with decreasing pressure, while the oil volume coefficient decreasing.

Hydraulic fracturing is a common practice employed to enhance reservoir permeability, particularly in instances where natural productivity is exceedingly low during production. Following the fracturing operations, it is customary to conduct well soaking procedures with the objective of enhancing the effectiveness of the fracturing process. This often results in a significant increase in the water content of the reservoir. In order to accurately reflect the complex flow dynamics of oil, gas, and water in such production conditions, the model must consider the three-phase flow of these fluids. Figure 5 depicts the oil-gas-water three-phase permeability curve of Gulong shale oil, which illustrates the correlation between the relative permeability and saturation of the fluid. In the relative permeability curve shown in Figure 5, the irreducible water saturation is about 10% and the residual oil saturation is about 40%. As the water saturation increases, the relative permeability of the oil phase decreases and the relative permeability of the water phase increases. As the gas saturation increases, the relative permeability of the oil phase decreases.

3.4. Flow Mechanism

3.4.1. Pseudo-Threshold Pressure Gradient

Shale characterized by tight reservoir, low porosity, small pore throat size and low permeability. In Gulong shale, the diameter of the matrix throat is between 10 and 100 nm. The permeability is ultra-low. Previous studies have found that fluids in tight reservoirs conform to low-velocity non-Darcy flow [29]. In the flow regime curve representing low-velocity non-Darcy flow, there is a nonlinear relationship between velocity and pressure when velocity falls below the threshold pressure. As shown in Figure 6, the linear segment intersects the x-axis at a point defined as the pseudo-threshold pressure gradient (PTPG).

The use of a PTPG reduces the difficulty of numerical simulation and has good applicability. Therefore, the PTPG is generally used to study low-velocity non-Darcy flows. In practical scenarios, the PTPG often replaces the threshold pressure gradient under the assumption that fluid flow stops when the pressure drops below the PTPG.

3.4.2. Stress Sensitivity

In shale reservoir development, the primary controlling factor of shale reservoir permeability is the permeability of reservoir fractures, including hydraulic fractures, opened beddings, microcracks, and so on. As fluid production progresses, the effective stress exerted on the reservoir rock rises, which leads to the closure of microcracks, large pores, or fractures, the rise of reservoir flow resistance, and the decline of reservoir permeability. Regarding stress sensitivity, numerous fitting formulas have been proposed by previous researchers, including the exponential model, power-law model, pore-shell model, logarithmic model, and plate fracture model. The exponential model is the most commonly utilized in shale reservoir simulations. The coefficient of the exponential model, referred to as the stress sensitivity coefficient, exerts a significant influence on the stress sensitivity effect of the fractures.

In setting up the model in this paper, the study area consists of open bedding and matrix. The opened beddings refer to the beddings that are widened under the action of high-pressure fluid during fracturing. The permeability of open bedding varies significantly with the reduction of internal pressure, and they have a large stress sensitivity coefficient. Within the matrix region, microcracks are present, and the permeability of these microcracks gradually decreases with decreasing pressure. The matrix region has a smaller stress sensitivity coefficient. The experiments on the Gulong shale cores were conducted by previous researchers. They measured the permeability under different effective stress, concluded that the stress sensitivity coefficient of Gulong shale open bedding was 0.33-0.47 [26]. The stress sensitivity coefficient of Gulong shale matrix is 0.010 [15], as shown in Figure 7.

3.5. Establishment of Numerical Model

A numerical model was constructed based on the assumption of the flow process of Gulong shale oil, utilizing the tNavigator software. The influence of the PTPG and stress sensitivity effect was considered in the modeling process. The results of the mesh generation are presented in Section 3.2, while the fluid parameters are shown in Section 3.3.

The PTPG of the grid in tNavigator can be defined using the keyword PTHRESHI. The ARITHMATIC keywords in tNavigator enable the real-time modification of matrix and bedding mesh permeability. The real-time permeability of the bedding mesh and matrix mesh is calculated using real-time pressure data. In the SCHEDULE part, the permeability of the opened bedding mesh and matrix mesh is updated after each step to simulate the stress sensitivity during the development process. The simulation of the subsequent time step will employ the updated permeability, and this process can be repeated to simulate stress-sensitive phenomena during development.

3.6. Verification of Model

A comparison of the flowback performance curve of the Gulong shale oil well A1 with the simulated flowback performance curve of the established multiphase flowback numerical model, as shown in Figure 8, reveals that while there are some deviations in the early period, the two curves match well in the later period. The matching degree is inadequate due to the effect of wellbore and hydraulic fractures, where has been ignored in our model, storage in the early flowback stage. Fluids within the wellbore and hydraulic fractures are discharged elastically immediately upon opening the well. However, over time, the impact of the wellbore and hydraulic fractures on the flowback performance gradually diminished. The degree of matching gradually becomes more accurate. This paper posits that the aforementioned model is applicable to shale oil in Gulong. The established flowback model allows for the realization of production dynamic adjustment, fracture volume inversion, and productivity prediction.

Table 1. Critical parameters for numerical simulation.

Parameter	numerical value	Parameter	numerical value
Thickness of reservoir (m)	40.0	Length of shale beddings (m)	100.0
Number of fracturing stages	27	Degree of opened bedding development (1/m)	0.2
Length of horizontal section(m)	2700.0	Grid permeability of opened beddings (mD)	20.0
Permeability of hydraulic fracture (D)	10.0	Width of opened beddings (m)	0.01
Porosity of matrix	0.10	Permeability of matrix (mD)	$5.8\times$10-4
Porosity of opened beddings	0.90	The initial pressure of the reservoir (MPa)	34.9
Temperature of reservoir(K)	395.0	Viscosity of crude oil (mPa·s)	2.19
Density of crude oil (kg/m3)	762.3	Water saturation of matrix	0.10
Water saturation of opened beddings (%)	100.0	Pseudo-Threshold Pressure Gradient (MPa/m)	10.0
Stress sensitivity coefficient of opened beddings	0.47	Stress sensitivity coefficient of matrix	0.10

Table 1. Critical parameters for numerical simulation.

Parameter	numerical value	Parameter	numerical value
Thickness of reservoir (m)	40.0	Length of shale beddings (m)	100.0
Number of fracturing stages	27	Degree of opened bedding development (1/m)	0.2
Length of horizontal section(m)	2700.0	Grid permeability of opened beddings (mD)	20.0
Permeability of hydraulic fracture (D)	10.0	Width of opened beddings (m)	0.01
Porosity of matrix	0.10	Permeability of matrix (mD)	$5.8\times$10-4
Porosity of opened beddings	0.90	The initial pressure of the reservoir (MPa)	34.9
Temperature of reservoir(K)	395.0	Viscosity of crude oil (mPa·s)	2.19
Density of crude oil (kg/m3)	762.3	Water saturation of matrix	0.10
Water saturation of opened beddings (%)	100.0	Pseudo-Threshold Pressure Gradient (MPa/m)	10.0
Stress sensitivity coefficient of opened beddings	0.47	Stress sensitivity coefficient of matrix	0.10

The discrepancy between the model outcomes and the operational performance of well A1 was subjected to a quantitative analysis, as shown in Figure 9. Prior to the 100th day, the actual production exceeded the simulated result due to the influence of wellbore accumulation. The model is more focused on the study of oil production in the later period, and thus the fitting results after the 100th day are used for quantitative analysis. From the 100th day to the 300th day, the average relative error between the oil rate and the gas rate was 2.51% and 11.69%, respectively, while the relative error of the water rate was -13.98%. The degree of correspondence between the 100th-day and 300th-day simulation results was calculated. The goodness of matching for the oil phase was found to be 0.93, while that for the gas phase was 0.72. It can be observed that the model demonstrates greater accuracy in forecasting oil and gas production in the latter stages.

4. Sensitivity Study of Flowback Performance

A numerical model was constructed to consider multiple mechanisms, and a sensitive factor analysis was conducted to address the issue of difficult prediction of oil breakthrough time and production in the Gulong shale oil reservoir. This chapter addresses the effects of PTPG, stress sensitivity, and permeability on flowback and production performance, as well as the oil breakthrough time. The influence of various factors on the yield at the model scale is shown in figures 10, 12, 14, and 16, respectively. The actual horizontal well production is calculated as the model production multiplied by the number of fracturing stages multiplied by the number of opened beddings, and then further multiplied by two.

4.1. Pseudo-Threshold Pressure Gradient Analysis

The time required for oil to breakthrough varies with different PTPG. As the PTPG increases, the oil breakthrough time exhibits a delayed trend. As shown in Figure 10.a, during the increasing stage, the rate of increase of the oil rate curves significantly slows down with the increase of the PTPG. Concomitantly, the peak oil rate also declines significantly, and the time corresponding to the peak oil rate exhibits a delayed trend. In the decreasing stage, the rate of decline of the oil rate curve significantly slows down with the increase of the PTPG. The total oil content at the 600th day exhibited a gradual decline as the PTPG increased. This is due to the fact that an increase in the PTPG results in a reduction in the capacity for oil to flow through the matrix, an increase in the pressure gradient within the matrix, and a narrowing of the oil supply range. As shown in Figure 10.b, PTPG exerts a profound influence on water rate. The smaller the PTPG, the larger the peak water rate, the earlier the peak water rate time, and the lower the water rate in the later period. With the increase of PTPG, the water rate curve tends to shift towards the lower right side. The larger the PTPG, the lower the peak water rate yield, the later the peak water rate time, and the smaller the fluctuation of water rate. As shown in Figure 10.c, the gas rate decreases with the increase of PTPG. As the initial reservoir pressure surpasses the saturation pressure, there is an absence of free gas within the matrix, with gas production primarily stemming from the degassing of crude oil. Consequently, the higher the PTPG, the lower the crude oil production, the lower the corresponding oil degassing, and therefore the lower the daily gas production.

The analysis of this feature indicates that the fluid primarily flows through the opened bedding during the initial stage. Due to the substantial quantity of fracturing fluid present in the opened bedding, the oil saturation in the opened bedding is relatively low, resulting in a period of low production. Additionally, the permeability of the opened bedding is considerable, rendering it unaffected by the PTPG. Consequently, the production curve remains largely consistent during the low production stage. In the increasing stage, oil and gas begin to flow and produce near the opened bedding. The remarkably low permeability of the matrix significantly influences fluid flow within it, predominantly governed by the PTPG. A substantial PTPG amplifies the fluid's flow resistance within the matrix, leading to a pronounced pressure gradient along the flow pathway (Figure 11). The small propagation range of pressure waves leads to the slow growth of the oil supply area. Subsequently, the rate of production growth slows and the peak production level declines with the increase in PTPG during the initial stages of production.

The peak production time is also advanced as the PTPG decreases, because the production is sensitive to the change in the PTPG, which leads to the difference of several times the oil production in the early stage. The peak production is the maximum value of the production and the node where the production increases and decreases. Due to the limited spread range of the reservoir and the limited recovery rate, the larger the spread range in the early stage, the larger the corresponding production, and the faster the fluid is extracted from the formation, the earlier the peak production time. The oil production in the early stage of small PTPG is much larger than that in the early stage of large PTPG, so the peak production time also advances as the PTPG decreases.In the decline stage, as the PTPG increases, the rate of fluid production in the matrix slows down. Consequently, the influence of fluid production on the production decline is minimal, resulting in a long, stable production period and a slow decline rate.

4.2. Opened Bedding Stress Sensitivity Analysis

As shown in Figure 12.a, as the stress sensitivity coefficient (S in Figure 11) increases, the increasing rate of the oil rate curve slows down, and the oil rate peak shows a decreasing trend. When the stress sensitivity coefficient goes up, the time to oil rate peak is constantly advanced. During the declining stage, there is a rising tendency in the rate of oil decline as the stress sensitivity coefficient increases. With an increase in the stress sensitivity coefficient, the total oil at the 600th day decreases. Influenced by the exponential relationship between effective stress and permeability of beddings, the declining of total oil production is accelerating, when the coefficient increases from 0.47 to 0.6. As recovery proceeds, the BHP decreases, and the opened beddings with larger stress sensitivity coefficients close more quickly. Their permeability decreases rapidly, and the pressure gradient in the shale beddings is more significant, resulting in a reduced supply range and lower production.

As shown in the Figure 12.b, during the initial stages, there was negligible disparity in the water rate curve. However, the peak water rate diminishes with the escalation of the stress sensitivity coefficient. The difference increased after the peak water rate. Large stress sensitivity coefficients lead to the opened bedding permeability decrease faster and the water rate decrease faster after water rate peak. As shown in Figure 12.c, the increase of stress sensitivity coefficient leads to a significant decrease in gas rate.

The analysis of this feature shows that in the increasing stage, the opened bedding pressure decreases with the production of the fluid in bedding and matrix. However, the pressure drop in the early stage is small, so the stress sensitive effect is weak, and the pressure drop propagation range is small. The permeability of the opened bedding decreases exponentially. The higher the stress sensitivity coefficient, the faster the permeability decline, resulting in slower production increasing and lower peak production with the increase of stress sensitivity coefficient. In the decline stage, due to the low opened bedding pressure, the opened bedding permeability is also very low. There is a large pressure gradient in the opened bedding (Figure 13), which greatly limits the subsequent pressure diffusion in the matrix. As the stress sensitivity coefficient increases, the oil supply area diminishes, leading to a faster rate of oil production decline.

4.3. Opened Bedding Permeability Analysis

As shown in Figure 14.a, the differences of oil rate begin to appear in the later incremental stage. As the permeability of the opened beddings (Kf in Figure14) increases, the oil rate peak shows an increasing trend; the difference in the time of the peak is not apparent. In the decreasing stage, the decreasing rate of the oil rate curve decreases with the increased opened bedding permeability. The higher opened bedding permeability means the greater the opening degree of the opened beddings. Therefore, its closing speed is slower under the stress sensitivity effect, and the corresponding permeability after closing is larger, which leads to a higher total oil.

As shown in the Figure 14.b, the water rate curves flow under different opened bedding permeability is similar to that under different stress sensitivity coefficients. There is no significant difference before reaching the water rate peak, and the difference of water rate peak is also small. The water rate is slightly larger with the higher the permeability of opened bedding in the decline stage. As shown in Figure 14.c, the decrease of opened bedding permeability leads to the decrease of gas production. However, the decline rate of gas rate curve under different opened bedding permeability is basically the same in the decline stage.

According to the analysis of this feature, in the increasing stage, with the production of fluid in the opened bedding and matrix, the pressure of the opened bedding decreases rapidly, and its permeability decreases exponentially. When the opened bedding pressure is equal, the larger the original permeability of opened bedding is, the larger the permeability after decreasing, so the greater the production peak can be reached. In the decline stage, although the original permeability of opened beddings varies, their stress sensitivity coefficients remain consistent, resulting in a similar trend in fracture permeability changes. When the permeability of opened beddings is high (>50×10^-3${\mu m}^2$), the yield difference during the decline stage is minimal initially but gradually increases with production progress. The variation in yield across different opened bedding permeability is not significant. However, differences mainly arise from challenges in fluid production due to fracture permeability within the supply area.

Figure 15. Formation pressure distribution on the 600th day at different opened bedding permeability (opened bedding permeability is respectively 10, 50, 100, 150, 200mD).

Figure 15. Formation pressure distribution on the 600th day at different opened bedding permeability (opened bedding permeability is respectively 10, 50, 100, 150, 200mD).

4.4. Matrix Permeability Analysis

As shown in Figure 16.a, the four curves are almost identical in the low yield stage. Differences appeared in the incremental stage. As the matrix permeability (Km in Figure16) increased, the increasing rate of the oil rate curve accelerated. The oil rate peak also showed a rising trend, and the time to oil rate peak was earlier. In the decreasing phase, the rate of decrease in the oil rate curve accelerated with increasing matrix permeability. Matrix permeability reflects the resistance to fluid mobilization within the matrix, resulting in a smaller supply range and lower oil rate, which is consistent with the impact of the PTPG.

As shown in the Figure 16.b, there are obvious differences in water production under different matrix permeability. With the increase of matrix permeability, the shape of water rate curve becomes narrower and higher. Higher matrix permeability can result in a longer increasing stage and a larger water rate peak. These reservoirs also tend to experience shorter high production stages and faster decline rates after reaching the water rate peak. As shown in Figure 16.c, the gas rate increases with the increase of matrix permeability.

The characteristics of liquid production are analyzed. In the incremental stage, with the production of fluid in the opened bedding and adjacent matrix, more fluid in the matrix is affected by the pressure drop and flow to the opened bedding. Greater matrix permeability correlates with reduced fluid mobilization and flow resistance within the matrix, consequently resulting in a decreased pressure gradient in the matrix. This results in a wider range of pressure drops that can reach, leading to faster production increasing rates and the ability to reach larger peak production.

Similar to the PTPG, the early production of shale oil is also very sensitive to the matrix permeability, and the pressure spread range in the early under different matrix permeability is very different, resulting in a large difference in the early oil production and the change time of production from increasing to decreasing. Therefore, the higher the matrix permeability, the faster the fluid is extracted from the matrix, and the corresponding advance of the peak production time.

In the decline stage, when the matrix permeability is large, the pressure wave can be diffused to a longer range. This means that in the decline phase, the fluid is mainly supplied by the matrix that is farther away from the opened bedding. Due to the large permeability of the matrix, fluid flow in the matrix is relatively easy, so the production decline is slow. These changes indicate that matrix permeability has a significant impact on both the increase and decline stages of shale oil production. The increase of permeability will enhance the fluid flow capacity and increase the oil supply area, thus promoting the increase of production. At the same time, the increase of matrix permeability will slow down the production decline rate and increase the cumulative production.

Figure 17. Formation pressure distribution on the 600th day at different matrix permeability (matrix permeability is respectively 1.8, 3.8, 5.8, 7.8, 9.8×10^-4mD).

Figure 17. Formation pressure distribution on the 600th day at different matrix permeability (matrix permeability is respectively 1.8, 3.8, 5.8, 7.8, 9.8×10^-4mD).

5. Optimization of Flowback Modes

A reasonable flowback strategy is a crucial factor in the efficient production of shale oil. In the context of oilfield operations, the flowback strategies of shale oil wells are frequently modified through the replacement of different-sized chokes. The study of choke size on shale oil flowback performance is limited to the regulation of the shale oil rate. The varying pipe resistance under different choke sizes leads to different pressure fields, which in turn affects the shale oil rate. It is commonly accepted that casing pressure and flowing bottom hole pressure exert a significant influence on oil rate. The bottom hole flow pressure exerts a direct influence on the inflow performance relationship; however, it is typically not feasible to measure this parameter directly in the field. In this chapter, we establish the relationship between the choke size and the casing pressure and between the casing pressure and the bottom hole flow pressure (BHP). This allows us to indirectly obtain the link between the choke size and the BHP, which is necessary to achieve a thorough study of the shale oil flowback strategies.

5.1. Flowback Modes for Shale Oil Wells

A selection of oil wells within the Gulong block were subjected to analysis. The study on the parameters of choke size, casing pressure, and BHP revealed that the decreasing rates of casing pressure and BHP vary with choke sizes. Consequently, three shale oil flowback modes are distinguished by the rate of decline in BHP [21], as shown in Figure 18.

As shown in Figure 18a, the fast-slow mode is distinguished by a rapid increase in the choke size during the initial phase, which is then maintained at a constant level. The rate of decline in casing and bottom hole flow pressure is initially rapid, followed by a slower rate of decline. As shown in Figure 18.b, the fast-slow-fast mode is defined by a rapid increase in choke size during the initial phase, followed by a medium choke size in the intermediate stage and a large choke size in the late stage. The rate of decrease in casing and bottom hole flow pressure is initially rapid, then gradual, and finally rapid again. As shown in Figure 18.c, the slow-fast-slow mode is characterized by a gradual increase in the size of the choke throughout the stage. The rate of decline in casing pressure and BHP is initially slow, followed by a period of accelerated decline, and then a subsequent slowdown.

The BHP exerts a direct influence on the inflow performance relationship (IPR). If the BHP is excessive, the production pressure differential will be insufficient, resulting in a reduction in the oil rate. Conversely, if the BHP is insufficient, the production pressure differential will increase. Consequently, the oil rate will be larger, accompanied by a series of problems such as crude oil degassing and reservoir stress sensitivity [31]. Consequently, the method of controlling BHP is employed to categorize the three flowback modes.

5.2. Optimization Method of Flowback Modes

In this paper, we employed a constructed numerical model to optimize three distinct flowback modes. Due to the inability to directly simulate the oil choke as a variable, we employed the BHP as a control variable and conducted simulations with variable BHP. This implies that control was achieved by modifying the lower limit of the BHP on a daily basis. During the simulation, the reservoir and fluid properties parameters were maintained at their original values, and a 1,000-day production simulation was conducted. The objective of this study was to evaluate the three different flowback modes based on the final oil production. By analyzing and comparing the simulation results, it is possible to determine the ultimate oil production and use this as a criterion for selecting the preferred flowback mode. The optimization process facilitates a more comprehensive understanding of the most suitable flowback strategy, which in turn leads to enhanced oil production and improved production efficiency.

5.3. The Impact of Different Flowback Modes

5.3.1. The Impact on Oil Break-through Time

The oil rate curves of the three flowback modes exhibit distinct characteristics, with the primary difference being the oil breakthrough time. The order of oil breakthrough time from shortest to longest is as follows: fast-slow-fast mode, fast-slow mode, slow-fast-slow mode. The oil breakthrough times for the aforementioned flowback modes are 29, 63, and 75 days, respectively. The analysis indicates that the oil breakthrough time is primarily influenced by the BHP at the outset of the flowback period. The presence of a PTPG within the shale matrix results in a differential pressure drop rate for the BHP during the early stages of flowback, contingent upon the specific flowback mode. A more rapid decline in pressure during the initial stages of the flowback period results in a greater pressure gradient within the matrix, which in turn leads to a shorter oil breakthrough time.

5.3.2. The Impact on Flowback Performance

The decline rate of bottom hole pressure (BHP) during the initial stages of production directly affects the rate of oil increase during that period. Consequently, the decline rate of BHP in the late stage determines the rate of oil decrease during that phase.

In the slow-fast-slow mode, the decline rate of BHP is relatively modest during the initial stages. Furthermore, there is a prolonged period of time during which no oil is produced. Upon transitioning to the fast stage, the BHP undergoes a pronounced increase. Consequently, the rate of oil increase is rapid. The time required to reach the production peak is the shortest, and the oil rate peak is the largest. The period of stabilized production is shorter, and the decrease in oil rate is faster after the period of stabilized production.

In the fast-slow mode, the rate of decrease in the bottom hole flow pressure is considerable during the initial stages. Concomitantly, the oil rate increases at a rapid pace. Subsequently, as the process enters the slow stage, the oil rate increases at a gradual pace, eventually reaching a point of decline. The slower rate of oil decrease is primarily attributable to the elevated BHP at the onset of the slow stage. Furthermore, the bottom hole flow pressure can continue to decrease in the subsequent slow stage, and the oil supply area can be further expanded. The period of stabilized production is lengthy, yet the stable oil rate is the lowest.

The fast-slow-fast mode is characterized by a significant decrease in BHP in the early stages, which facilitates the breakthrough of oil. Subsequently, the oil rate increases at a gradual pace as the well transitions into the slow stage. Upon re-entry into the fast stage, a notable increase in oil rate is observed. The period of stabilized production is relatively brief, yet the stabilized oil rate is high.

Figure 19. Oil rate curves of three flowback modes.

Figure 19. Oil rate curves of three flowback modes.

A comparison of the water-cut rate curves presented in Figure 20 reveals notable discrepancies in the water-cut rate during the medium stage. At the outset of the flowback process, the water-cut rate of the slow-fast-slow mode exhibited the most rapid decline, while the water-cut rate of the fast-slow-fast mode exhibited the slowest decline. As the development time progresses, the disparity in water-cut rate between the three modes gradually diminishes. As production progresses, the water-cut rate of the fast-slow mode gradually becomes the largest, although the difference between the rates remains minimal.

The data presented in Figure 21 indicates that various flowback modes exert a considerable influence on the gas-oil ratio. Given the sensitivity of the dissolved gas-oil ratio to pressure changes, particularly below the bubble point pressure, fluctuations in bottom-hole flow pressure have a notable impact on the gas-oil ratio. During the initial production phase, the gas-oil ratio of the fast-slow mode exhibits an initial increase, and the time required to maintain a low gas-oil ratio is longer. In the slow phase (200th day), the gas-oil ratio exhibits a rapid increase to a higher level. In the fast-slow-fast mode, the production process commences with a low gas-oil ratio until the fast phase (200th day), at which point the gas-oil ratio rises rapidly to a higher level. In the slow-fast-slow mode, the gas-oil ratio rises rapidly in the fast phase (150th day) and remains at a high level. However, after the 300th day, the gas-oil ratio decreases at a faster rate than in the other two modes. In conclusion, the choice of production system model will have a significant impact on the gas-oil ratio. According to the characteristics of these models, a suitable system model can be selected to control the gas-oil ratio in order to achieve more efficient production.

5.4. Optimization Results and Discussion

As shown in Figure 19, the initial production discrepancy between the three flowback modes is considerable. As shown in Figure 22, the total oil differences between the three flowback modes are the greatest in the initial stages. Over time, the discrepancy in total oil between the three modes diminishes gradually. The fast-slow mode exhibited the greatest total oil on the 1000th day.

As shown in Figure 23, the three flowback modes were optimized using the total oil and recovery ratio on the 1000th day. The total oil produced by the three flowback modes, namely fast-slow, slow-fast-slow, and fast-slow-fast, was 9.6211 m³, 9.4783 m³, and 9.4331 m³, respectively. The corresponding recovery ratio of oil was 1.2231%, 1.2050%, and 1.1992%, respectively. From the perspective of the most significant total oil production and a highest recovery ratio, the results of the three flowback modes are as follows: the fast-slow mode is the most optimal, the slow-fast-slow mode is the second most optimal, and the fast-slow-fast mode is the least optimal.

A correlation has been identified between choke size, casing pressure, and BHP, which allows for the derivation of the optimal choke mode among the three modes. The most optimal choke mode is characterised by a rapid increase in choke size in the initial stage, followed by a maintenance of a medium choke size.

In the case of shale, due to its extremely low permeability and significant PTPG in the matrix, the movable oil in shale reservoirs is primarily distributed within the volume containing fractures. Furthermore, within the volume containing fractures, only the oil in close proximity to the fractures exhibits a superior mobility effect. As oil is recovered near the fractures, the rate of shale oil production begins to decline. The primary focus of production is on the expansion of the oil supply area and the further recovery of crude oil in the vicinity of the fractures, which are in a decline phase.

There is a relationship between oil rate and pressure spread range, the larger the pressure spread range, the higher the oil production. As shown in Figure 24, Figure 25 and Figure 26, the production performance of the three models was evaluated based on the 300th day and 1000th day pressure distribution maps. Upon analysis of the pressure distribution maps at the 300th and 1000th days, it was observed that the pressure spread range exhibited considerable variation across the three production modes, ranging from large to small in the following order: slow-fast-slow, fast-slow, fast-slow-fast. Conversely, at the 1000th day, the pressure spread range exhibited minimal variation, with the three modes displaying a similar pattern of large to small in the following order: fast-slow, slow-fast-slow, fast-slow-fast. The outcomes of the pressure spread range are in accordance with those of the oil rate.

In slow-fast-slow mode, the initial pressure spread range is the largest, the crude oil production range is the largest, and the initial daily oil rate and total oil are the largest. However, the rapid reduction of bottom hole pressure leads to the rapid closure of the opened beddings, which in turn increases the flow resistance of the opened bedding. A large amount of crude oil degasses, and the oil rate fluctuation is large, which is not conducive to the later crude oil recovery. In the latter stages of production, the pressure spread is to some extent affected. The late pressure spread range is relatively narrow, the oil supply range is similarly limited, and the oil rate is similarly low.

In fast-slow-fast mode, the initial pressure spread range is the smallest, the crude oil supply range is similarly limited, and the oil rate and total oil are similarly low in the early stages. The initial production is similarly low, which makes the oil rate more stable. In the latter stages of production, a precipitous decline in bottomhole pressure results in the abrupt closure of previously opened beddings. The pressure spread range is narrow, the oil production range is limited, and the oil rate in the latter stages is minimal.

In the fast-slow mode, the initial pressure spread range is situated between the fast-slow mode and the slow-fast-slow mode. The initial oil rate is lower than that observed in the slow-fast mode, and the degasser phenomenon is less pronounced in the reservoir. The period of stable oil rate is the longest. In the subsequent production phase, the pressure drop spread range is the greatest, during which the oil supply area within the matrix expands effectively and the oil recovery area is also high. In the fast-slow mode, the BHP drops more rapidly in the initial stages. Furthermore, the fracturing fluid can be recovered with minimal delay following the fracturing process. The well proceeded rapidly to the production stage. The oil recovery ratio of the matrix in the vicinity of the fracture is high, and the oil recovery period is prolonged. The stable choke size is conducive to stable oil rate and the longest stable oil rate period. Consequently, the fast-slow mode yields the highest total oil yield.

6. Conclusions

• The multiphase flowback numerical model of the Gulong shale reservoir, which consides PTPG and stress sensitivity, has been established with high applicability in the later flowback stage. From the 100th day to the 300th day, the average relative error between oil rate and gas rate is 2.51% and 11.69%. The goodness of matching of oil is 0.93, and that of gas is 0.72.
• PTPG and matrix permeability are indicative of the difficulty of crude oil flow in the shale matrix during the flowback period, which results in significant fluctuations in oil rate. The PTPG is responsible for the alteration of oil breakthrough time, with the larger the PTPG, the longer the oil breakthrough time. The increase of PTPG and the decrease of matrix permeability result in a larger pressure gradient in the matrix and a smaller pressure spread range. This, in turn, leads to a reduction in the peak oil rate, an advance in the peak time, and an acceleration in the oil decline rate.
• The permeability and stress sensitivity of the opened bedding are indicative of the difficulty of crude oil flow in opened bedding. The stress sensitivity of the shale is indicative of the speed at which the opened bedding is closing, which in turn affects the conductivity and pressure gradient within the opened bedding. An increase in the stress sensitivity coefficient of opened bedding results in a reduction in conductivity and an increase in the pressure gradient within the opened bedding. Consequently, the peak oil rate declines, accompanied by an increase in the rate of oil depletion. The permeability of opened bedding is also similar.
• The production data from Daqing Gulong shale oil wells indicates that three flowback modes can be identified and classified. These are the fast-slow mode, the slow-fast-slow mode, and the fast-slow-fast mode, which correspond to different choke modes. The total oil and recovery ratio are used to evaluate the production of different modes. The results of a 1000-day numerical simulation indicate that the flowback characteristics of the fast-slow mode are as follows: oil is observed earlier in the early stage, production is larger in the early stage, and the stable oil rate period is the longest. The flowback characteristics of the fast-slow-fast mode are as follows: the earliest oil breakthrough time, the lowest production in the early stage, and the highest production in the stable production stage. The slow-fast-slow mode flowback characteristics are as follows: the latest oil breakthrough time, the largest production in the early stage, and the production fluctuation.
• The total oil recovered from the three modes of flowback, namely fast-slow, slow-fast-slow, and fast-slow-fast, was 9.6211 m³, 9.4783 m³, and 9.4331 m³, respectively. The corresponding recovery ratio of oil was 1.2231%, 1.2050%, and 1.1992%, respectively. In order to achieve the greatest total oil yield, the optimal bottom hole pressure mode is identified as the fast-slow mode. The optimal choke flowback mode is characterised by a rapid increase in choke size during the initial stage, followed by a maintenance of a medium size. The optimal fast-slow mode is characterised by an early oil breakthrough time, a high oil production rate in the initial stages, a long stable production period, a low gas and oil production ratio, and a large oil supply area. This study has certain significance for determining reasonable flowback strategies in the Gulong block.