Fracture Spacing Optimization Method for Multi-Stage Fractured Horizontal Wells in Shale Oil Reservoir Based on Dynamic Production Data Analysis

Abstract

In order to improve the shale oil production rate and save fracturing costs, based on dynamic production data, a production-oriented optimization method for fracture spacing of multi-stage fractured horizontal wells is proposed in this study. First, M. Brown et al.’s trilinear seepage flow models and their pressure and flow rate solutions are applied. Second, deconvolution theory is introduced to normalize the production data. The data of variable pressure and variable flow rate are, respectively, transformed into the pressure data under unit flow rate and the flow rate data under unit production pressure drop; and the influence of data error is eliminated. Two kinds of typical curve of the normalized data are analyzed using the pressure and flow rate solutions of M. Brown et al.’s models. The two fitting methods constrain each other. Thus, reservoir and fracture parameters are interpretated. A practical model has been established to more accurately describe the seepage flow behavior in shale oil reservoirs. Third, using Duhamel’s principle and the rate solution, the daily and cumulative production rate under any variable production pressure can be obtained. The productivity can be more accurately predicted. Finally, the analysis method is applied to analyze the actual dynamic production data. The fracture spacing of a shale oil producing well in an actual block is optimized from the aspects of production life, cumulative production, economic benefits and other influencing factors, and some significant conclusions are obtained. The research results show that with the goal of maximum cumulative production, the optimal fracture spacing is 5.5 m for 5 years and 11.4 m for 10 years. All in all, the fracture spacing optimization and design theory of multi-stage fractured horizontal wells is enriched.

1. Introduction

The global shale reservoir resources are rich and valuable to exploit [1]. However, shale reservoirs have a low permeability, which makes crude oil flow and reservoir production difficult [2]. Therefore, shale reservoirs are often exploited by the multi-stage fracturing of horizontal wells [3]. Through hydraulic fracturing, the reservoir can be stimulated, several main fractures are formed perpendicular to the horizontal wells and the fracture network is formed around the main fractures. Therefore, the stimulated reservoir volume area is formed. Reservoir permeability and porosity are increased, oil flows more easily and oil recovery efficiency is improved. The higher the fracture numbers, the higher the production rate [4], but the higher the production cost [5,6]. The fracture spacing not only needs to meet the requirements of a high production rate, but also needs to make the production cost not too high. Therefore, it is very necessary to optimize the fracture spacing [7]. The schematic diagram of fracture spacing optimization is shown in Figure 1. By optimizing the fracture spacing, efficient reservoir development can be achieved.

In recent years, scholars have carried out a series of studies for optimizing fracture spacing [8,9,10,11,12,13,14,15]. From 2017 to 2022, some scholars took fracture propagation in the process of hydraulic fracturing as the research object. They presented the numerical model considering elastic fluid mechanics and stress disturbances and different fracture flow distributions [8], the computational model of embedded discrete fractures [9], the mathematical models of the coupling effect of rock and fluid dynamics [10], the computational optimization model based on intelligent variable-fidelity radial basis function [11], the mathematical model of fully coupled deformation and seepage flow in porous media [12], the 3D solid mechanics model of hydraulic fracturing [13], the numerical model considering fracture geometry and proppant flow dynamics [14] and the numerical model considering stress variation with depth [15] in order to optimize the fracture spacing in shale reservoirs from different research angles. However, these models involve fluid viscosity, reservoir permeability, reservoir porosity, reservoir thickness, fracture width, fracture length, fracture number, total stress, effective stress, rock storage coefficient and other parameters, which are difficult to obtain in the actual production process.

With the development of field measurement technology, a large amount of on-site pressure and production rate data can be obtained. The reservoir and fracture parameters can be interpreted by inversion method using the dynamic production data. Then the production rate can be calculated through the forward computation of the seepage model by using the interpreted values of reservoir and fracture parameters. Then the fracture spacing can be optimized. In 2006, through modeling the flow in each streamline independently in real time, the Wang–Kovscek [16] streamline method for production data inversion had been improved by Vegard R. Stenerud and Nut-Andreas Lie [17]. The results showed that this method had better matching and faster convergence rate. In 2016, aiming at the problem that nonlinearity and variable production rate should be considered when interpreting production data of shale gas reservoirs, a classical trilinear flow model was modified and a method for comprehensively analyzing variable production rate data was proposed by Wu et al. [18]. In this method, considering the desorption and the nonlinearity of compressibility, modified material balance equation and material balance time were used to process the production data. It was proved through a field case that this method could more accurately interpret the production data. In 2017, aiming at the problem of fracture inversion, a fracture inversion method was proposed by using production data based on Griffith failure criterion and ground stress correlation by Zhang et al. [19]. Theoretical examples showed that this method was effective for the accurate inversion of fractures, but as the fracture numbers are more, the inversion results become worse. In 2018, aiming at the problem of the significant discontinuities in production data caused by frequent shut-ins, a new production data analysis method for solving the discontinuous problem based on pseudo time was proposed by Li et al. [20]. Duhamel’s principle, Laplace transform and inversion and the Newman method were used to solve the model used for production data analysis, and the analytical and numerical solutions were verified. The results showed that this method had great potential in estimating formation parameters and predicting the well production dynamics more effectively. In 2021, an improved spatial inversion method of data was proposed by Liu et al. [21]. The reservoir state fields can be quickly predicted by observing the production data. The method was also tested in the field. The results showed that this method had high computational efficiency and accuracy. The above studies [16,17,18,19,20,21] proposed some new inversion methods of dynamic production data or improved the existing methods for interpretating reservoir parameters or fracture parameters after fracturing, and a certain theoretical basis for the dynamic production data inversion technology of multi-stage fractured horizontal wells was provided. However, due to the dramatic changes in flow pressure and the production rate in unconventional oil and gas production data with large errors [22], the normalized typical data points in the above dynamic production data inversion method were scattered, smooth typical curves were difficult to obtain and the data fitting effect was also poor, which resulted in great uncertainty in the fitting results. In addition, the interpreted post-hydraulic fracturing models of seepage flow during production in the aforementioned studies was rarely further applied to the optimization of productivity enhancement in the oil field.

In 2020, Mohammed and Joseph combined data analysis with theoretical models to establish a hybrid hydraulic fracturing model that combines data and theory. The results showed that the hybrid model has higher accuracy. It is feasible to combine data analysis with theoretical models [23]. Therefore, based on the dynamic production data inversion, a new production-oriented optimization method for the fracture spacing of multi-stage fractured horizontal wells in shale oil is proposed. In particular, deconvolution algorithm [24,25,26,27,28,29,30,31] is introduced to normalize the pressure data. Not only can the data of variable pressure and variable flow rate be directly converted into pressure data under unit flow, but also the regularization of deconvolution calculation can be performed, which can eliminate the influence of data error and expand the investigation distance of production data analysis. As a result, more information for production data analysis can be provided, and thus the fitting effect is improved and the uncertainty of parameter interpretation is reduced. The main research contents of this study include: first, a three-linear seepage mathematical model for multi-stage fractured horizontal wells in shale reservoir, the pressure solution under constant flow rate and the flow rate solution under constant production pressure in Laplace space are introduced [32]. Second, the abundant on-site production data are fully utilized for dynamic production data inversion, and the deconvolution theory to normalize the production data is also introduced. By referring to the specific algorithm of pressure deconvolution for data normalization, the data of variable pressure and variable flow rate are converted into the pressure data under unit flow rate, and the influence of data errors is also eliminated. According to the pressure under unit flow rate, the typical curve analysis of the pressure data under unit flow rate is carried out. The reservoir parameters and fracture parameters after hydraulic fracturing are interpreted, so that the interpreted seepage model is more in line with the reality and the seepage flow behavior can be represented more accurately. Then, the Duhamel’s principle and the analytical solution of the interpretation model are used to calculate the flow rate per unit production pressure drop. The daily and cumulative production rate of horizontal wells under any production pressure system can be obtained, which can predict the productivity more accurately and efficiently. According to the productivity obtained, the fracture spacing can be optimized, and an optimization method for the fracture spacing of multi-stage fractured horizontal wells is proposed. Finally, the proposed fracture spacing optimization method was used to analyze the dynamic production data of a shale oil production well in the actual block. The fracture spacing was optimized from the aspects of production life, cumulative production, total economic benefit [33], balance of payments, fracturing cost, oil price and other influential factors [34]. The optimization method of fracture spacing proposed in this paper has its own advantages compared with the optimization method of fracture spacing based on fracture propagation in solid mechanics [35], and they can complement each other. By comprehensively utilizing these two methods, better fracture spacing can be obtained. Significant reference for the design of adjacent well fracture spacing in the same block in the future is provided. Some technical guidance is provided for later production and secondary fracturing of reservoirs.

2. Evaluation of Optimal Fracture Spacing for Multi-Stage Fractured Horizontal Wells in Shale Oil

2.1. Mathematical Model of Shale Oil Seepage Flow in Multi-Stage Fractured Horizontal Wells

The physical model of trilinear seepage flow in multi-stage fractured horizontal wells in shale reservoirs [32] is shown in Figure 2. The wellbore direction of the horizontal well is parallel to the Y-axis and the radius of the wellbore is r_w. A number of equally spaced primary fractures have been hydraulically fractured perpendicular to the wellbore (i.e., along the x direction). All fractures penetrate the reservoir completely. The top and bottom of the reservoir are closed and the ambient temperature is constant. Fluid flow in the reservoir is divided into reservoir flow area, inter-fracture flow area and main fracture flow area. The fluid first flows from the reservoir flow area into the inter-fracture flow area, and then flows from the inter-fracture flow area into the main fracture and finally flows through the main fracture into the wellbore of the horizontal well. The horizontal well is in production with constant flow rate and variable pressure. Due to geometric symmetry, the inter-fracture interference is not considered.

In 2011, a mathematical model of trilinear seepage flow in multi-stage fractured horizontal wells in shale reservoir based on the above physical model was established by M. Brown et al. [32]. The dimensionless variables in M. Brown et al.’s model is defined as follows:

$$\begin{array}{l}{C}_{\mathrm{FD}}=\frac{k_{\mathrm{F}}{w}_{\mathrm{F}}}{k_{\mathrm{I}}{x}_{\mathrm{F}}}; y_{\mathrm{e}}=\frac{d_{\mathrm{F}}}{2}; {\eta }_{\mathrm{FD}}=\frac{{\eta }_{\mathrm{F}}}{{\eta }_{\mathrm{I}}}; {\eta }_{\mathrm{OD}}=\frac{{\eta }_{\mathrm{O}}}{{\eta }_{\mathrm{I}}}; x_{\mathrm{D}}=\frac{x}{x_{\mathrm{F}}}; y_{\mathrm{D}}=\frac{y}{x_{\mathrm{F}}}; x_{\mathrm{eD}}=\frac{x_{\mathrm{e}}}{x_{\mathrm{F}}}; \\ w_{\mathrm{D}}=\frac{w_{\mathrm{F}}}{x_{\mathrm{F}}}; q_{\mathrm{F}}=\frac{q}{n_{\mathrm{F}}}; {\eta }_{\mathrm{F}}=\frac{k_{\mathrm{F}}}{{\varphi }_{\mathrm{F}}{c}_{\mathrm{t}}{}_{\mathrm{F}}\mu }; {\eta }_{\mathrm{I}}=\frac{k_{\mathrm{I}}}{{\varphi }_{\mathrm{I}}{c}_{\mathrm{t}}{}_{\mathrm{I}}\mu }; {\eta }_{\mathrm{O}}=\frac{k_{\mathrm{O}}}{{\varphi }_{\mathrm{O}}{c}_{\mathrm{t}}{}_{\mathrm{O}}\mu }; y_{\mathrm{eD}}=\frac{y_{\mathrm{e}}}{x_{\mathrm{F}}}; \\ t_{\mathrm{D}}=\frac{{\eta }_{\mathrm{I}}}{x_{\mathrm{F}}{}^2}t; p_{\mathrm{OD}}=\frac{2\mathsf{\pi }{k}_{\mathrm{l}}h}{q_{\mathrm{F}}B\mu }\left(p_{\mathrm{ini}}-p_{\mathrm{o}}\right); p_{\mathrm{ID}}=\frac{2\mathsf{\pi }{k}_{\mathrm{l}}h}{q_{\mathrm{F}}B\mu }\left(p_{\mathrm{ini}}-p_{\mathrm{I}}\right); \\ p_{\mathrm{FD}}=\frac{2\mathsf{\pi }{k}_{\mathrm{l}}h}{q_{\mathrm{F}}B\mu }\left(p_{\mathrm{ini}}-p_{\mathrm{F}}\right); C_{\mathrm{D}}=\frac{C}{2\mathsf{\pi }{\left(\varphi c_{\mathrm{t}}\right)}_{\mathrm{I}}h{x}_{\mathrm{F}}{}^2}.\end{array}$$

where C_tO is the comprehensive compression coefficient of seepage flow in the reservoir flow area, atm⁻¹; C_tI is the comprehensive compression coefficient of seepage flow in the inter-fracture flow area, atm⁻¹; C_tF is the comprehensive compression coefficient of seepage flow in the main fracture flow area, atm⁻¹; ϕ_O is the porosity in the reservoir flow area; ϕ_I is the porosity in the inter-fracture flow area; ϕ_F is the porosity in the main fracture flow area; k_O is the permeability in the reservoir flow area, D; k_I is the permeability in the inter-fracture flow area, D; k_F is the permeability in the main fracture flow area, D; n_F is the number of primary fractures; q is the horizontal well production, cm³/s; B is the volume coefficient; p_w is the bottom hole pressure, atm; w_F is the fracture width, cm; t_D is dimensionless time; x_D is the dimensionless distance in the x direction; x_eD is the dimensionless outer boundary distance; η_OD is a defined dimensionless variable; y_D is the dimensionless distance in the y direction; w_D is the dimensionless fracture width; y_eD is the dimensionless distance at 1/2 of the fracture spacing; η_FD is a defined dimensionless parameter; C_FD is a defined dimensionless quantity; and C_D is the dimensionless wellbore storage factor. The unit system of the formulas in the manuscript is the Darcy unit system.

The Laplace transformation solution of dimensionless pressure in reservoir flow area is as follows [32]:

$${\tilde{p}}_{\mathrm{OD}}={{\tilde{p}}_{\mathrm{ID}}|}_{x_{\mathrm{D}}=1}\cdot \frac{\cosh \left[\sqrt{s/{\eta }_{\mathrm{OD}}}\left(x_{\mathrm{eD}}-x_{\mathrm{D}}\right)\right]}{\cosh \left[\sqrt{s/{\eta }_{\mathrm{OD}}}\left(x_{\mathrm{eD}}-1\right)\right]},$$

(1)

where ${\tilde{p}}_{\mathrm{OD}}$ is Laplace transformation of dimensionless pressure p_OD in the reservoir flow area; ${\tilde{p}}_{\mathrm{ID}}$ is the Laplace transformation of dimensionless pressure p_ID in the inter-fracture flow area; s is the complex variable of Laplace transformation.

The Laplace transform solution of dimensionless pressure in the inter-fracture flow area is [32]:

$${\tilde{p}}_{\mathrm{ID}}={{\tilde{p}}_{\mathrm{FD}}|}_{y_{\mathrm{D}}=w_{\mathrm{D}}/2}\cdot \frac{\cosh \left[\sqrt{{\alpha }_{\mathrm{O}}}\left(y_{\mathrm{eD}}-y_{\mathrm{D}}\right)\right]}{\cosh \left[\sqrt{{\alpha }_{\mathrm{O}}}\left(y_{\mathrm{eD}}-\frac{w_{\mathrm{D}}}{2}\right)\right]},$$

(2)

where ${\tilde{p}}_{\mathrm{FD}}$ is Laplace transformation of dimensionless pressure p_FD in main fracture flow area; α_O is a defined dimensionless parameter.

The Laplace transform solution of dimensionless pressure in main fracture flow area is as follows [32]:

$${\tilde{p}}_{\mathrm{FD}}=\frac{\mathsf{\pi }}{C_{\mathrm{FD}}s\sqrt{{\alpha }_{\mathrm{F}}}}\cdot \frac{\cosh \left[\sqrt{{\alpha }_{\mathrm{F}}}\left(1-x_{\mathrm{D}}\right)\right]}{\sinh \left(\sqrt{{\alpha }_{\mathrm{F}}}\right)},$$

(3)

where α_F is a defined dimensionless parameter.

According to Equation (3), the dimensionless bottom hole pressure ${\tilde{p}}_{\mathrm{wD}}$ at a constant flow rate can be obtained as follows [32]:

$${\tilde{p}}_{\mathrm{wD}}={\tilde{p}}_{\mathrm{FD}}\left(x_{\mathrm{D}}=0\right)=\frac{\mathsf{\pi }}{C_{\mathrm{FD}}s\sqrt{{\alpha }_{\mathrm{F}}}\tanh \left(\sqrt{{\alpha }_{\mathrm{F}}}\right)},$$

(4)

where ${\tilde{p}}_{\mathrm{wD}}$ is the dimensionless bottom hole pressure.

The flow in the fracture of the trilinear seepage flow model is one-dimensional linear flow. However, the fluid flow along the fracture surface into the wellbore of the horizontal well will produce the radial flow. For solving the contradiction, a choking skin factor for approximating the choking resistance generated by the radial flow was proposed by Mukherjee and Economides [36]. The choking skin factor s_c can be calculated as follows [36]:

$$s_{\mathrm{c}}=\frac{k_{\mathrm{I}}h}{k_{\mathrm{F}}{w}_{\mathrm{F}}}\left[\ln \left(\frac{h}{2r_{\mathrm{w}}}\right)-\frac{\mathsf{\pi }}{2}\right],$$

(5)

According to Equations (4) and (5), the dimensionless bottom hole pressure with constant flow rate considering the radial flow at the fracture surface can be obtained as follows [32]:

$${\tilde{p}}_{\mathrm{wD}}=\frac{\mathsf{\pi }}{C_{\mathrm{FD}}s\sqrt{{\alpha }_{\mathrm{F}}}\tanh \left(\sqrt{{\alpha }_{\mathrm{F}}}\right)}+\frac{s_{\mathrm{c}}}{s},$$

(6)

In order to make the solution more practical, the wellbore storage effect should be considered. ${\tilde{p}}_{\mathrm{wD}}$ in Equation (6) can be substituted into the following convolution expression in the Laplacian domain [32]:

$${\tilde{p}}_{\mathrm{wD},\mathrm{storage}}=\frac{{\tilde{p}}_{\mathrm{wD}}}{1+C_{\mathrm{D}}{s}^2{\tilde{p}}_{\mathrm{wD}}},$$

(7)

Applying Duhamel’s principle, Equation (7) can be used to obtain the dimensionless flow rate solution per unit production pressure drop as follows [32]:

$$\widetilde{q_{\mathrm{D}}}=\frac{1}{{\tilde{p}}_{\mathrm{wD}}}\cdot \frac{1}{s^2}\cdot p_{\mathrm{wD},\mathrm{const}},$$

(8)

where $\widetilde{q_{\mathrm{D}}}$ is the Laplace transformation of dimensionless flow q_D corresponding to fixed bottom hole pressure; P_wD,const is a fixed dimensionless bottom hole pressure.

The dimensionless variables in the above model solution are defined as follows:

$$p_{\mathrm{wD},\mathrm{const}}=\frac{2\mathsf{\pi }{k}_{\mathrm{l}}h}{q_{\mathrm{F}}B\mu }; q_{\mathrm{D}}=\frac{q_{\mathrm{u}}}{q_{\mathrm{F}}}$$

where q_u is the flow rate corresponding to the fixed bottom hole pressure, cm³/s; q_D is the dimensionless flow rate corresponding to the fixed bottom hole pressure.

Equations (7) and (8) are, respectively, the pressure under constant flow and the flow rate under constant pressure in Laplace space, but not the solution in real space. Therefore, the Stehfest algorithm [37,38] is introduced to invert the solution of Laplace space to obtain the solution of the real space [39]. The calculation formula of Stehfest algorithm is as follows [37,38]:

$$f\left(T\right)=\frac{\ln 2}{T}{\displaystyle \sum _{i=1}^N{V}_i}\overline{f}\left(\frac{\ln 2}{T}i\right),$$

(9)

where:

$$V_i={\left(-1\right)}^{N/2+i}{\displaystyle \sum _{k=\left[\frac{i+1}{2}\right]}^{Min\left(i,N/2\right)}\frac{k^{N/2+1}\left(2k\right)!}{\left(N/2-k\right)!k!\left(k-1\right)!\left(i-k\right)!\left(2k-i\right)!}},$$

(10)

where the larger value of N is, the more accurate the calculation will be; generally, N is an even integer between 6 and 18 [39].

2.2. Deconvolution

In the actual production process of shale oil, due to the low permeability of shale reservoir and the change in flow dynamics in the production process, it is difficult to keep the pressure and flow rate of dynamic production data constant in reality. However, the mathematical model used in this article is of constant flow rate or constant pressure. The deconvolution algorithm is introduced to normalize the bottom hole flow pressure data, and then the data of variable pressure and variable flow rate are transformed into the pressure data under unit flow rate; the influence of data errors can also be eliminated by the data normalization process, so more information for production data analysis can be provided and ultimately the fitting effect can be improved. The deconvolution principle for well test interpretation or production dynamics data analysis is as follows:

According to Duhamel’ s principle, the pressure-flow rate convolution relation is obtained as follows [25]:

$$p_{\mathrm{ini}}-p_{\mathrm{wD}}\left(t\right)={\displaystyle {\int }_0^{t}q\left(t-\tau \right)}\frac{\partial p_{\mathrm{u}}\left(\tau \right)}{\partial t}\mathrm{d}\tau ,$$

(11)

where p is the bottom hole pressure under variable flow rate, atm; t is the production time, s; p_u is the flow response per unit flow rate, atm.

According to Duhamel’s principle, the flow rate function under variable bottom hole pressure is obtained as follows [40]:

$$q\left(t\right)={\displaystyle {\int }_0^t\Delta p_{\mathrm{wD}}\left(t-\tau \right)\cdot q^{\prime }{}_{\mathrm{u}}\left(\tau \right)\mathrm{d}\tau }$$

(12)

where Δp_wD is the production pressure drop, atm.

In the case of known variable flow rate q and bottom hole pressure p under the variable flow q rate, Equation (11) can be used to obtain the transient pressure response p_u in the oil reservoir for the whole production time [24]). It is worth noting that when using deconvolution calculations in practical applications, it is necessary to exclude the influence of stimulation measures during production, and the changes in reservoir physical property, fluid property and variable wellbore storage effect [25,26]. A series of studies on the deconvolution algorithm for inversion of reservoir production data [25,26,27,28,29,30,31] were carried out by many scholars. In 2004, based on Tikhonov regularization objective function [27], a new deconvolution algorithm was proposed by Schroeter et al. [28]. In 2006, the practical application of a B-spline-based deconvolution algorithm in well test analysis was investigated by Ilk et al. [30,31]. From 2017 to 2018, an improved B-spline deconvolution algorithm was proposed by Liu et al. [25,26]. Adding a nonlinear regularization based on curvature minimization, the stability of B-spline deconvolution algorithm was improved by Liu et al. [25,26]. Some theoretical bases for the practical application of deconvolution algorithm were provided by the above researchers. In this study, the improved B-spline deconvolution algorithm [25,26] is used, which has fast calculation speed and good stability. Through this deconvolution algorithm, the variable flow rate and variable pressure data can be transformed into the pressure data under unit flow rate and flow rate data under unit production pressure drop, and the error can be eliminated [26].

The specific pressure deconvolution algorithm [25,26] is as follows:

(1)

The pressure derivative per unit flow rate is reconstructed by employing the IlK second-order B-spline function weight sum.

(2)

The convolution integral property is adopted so that the sensitivity matrix for the deconvolution calculation can be solved quickly and analytically.

(3)

The idea of curvature minimization is introduced to increase nonlinear constraints, which reduces errors and improves stability.

The specific algorithm can be referred to in the literature [25,26].

The deconvolution algorithm applied to production decline [26] is as follows:

(1)

The flow rate derivative per unit production pressure drop is reconstructed by employing the IlK second-order B-spline function weight sum.

(2)

The convolution integral property is adopted so that the sensitivity matrix for the deconvolution calculation can be solved quickly and analytically by piecewise integration according to the pressure drop section.

The specific algorithm can be referred to in the literature [26].

2.3. Optimization Method of Fracture Spacing in Multi-Stage Fractured Horizontal Well of Shale Oil

In the actual production process of the reservoir, the pressure and flow rate in the formation are not constant. However, the flow model introduced in this study has a constant flow, so a method based on deconvolution is proposed to optimize the fracture spacing in shale oil. The operational and constant parameters are the reservoir parameters and fracture parameters explained through production data, fracture pacing and number of fractures, fracturing cost and oil price. The optimal fracture spacing can be determined with production and economic benefits as the objective functions. The operational and constant parameters can be obtained through seismic data, production data analysis and experimental testing. Firstly, by using the actual production data in the field and Duhamel’s principle, the deconvolution algorithm is used to normalize these actual production data by Equations (11) and (12), so that the actual production data with variable flow rate and variable pressure can be transformed into pressure data under unit flow rate and flow rate data under unit production pressure drop, and the influence of data error can also be eliminated. Then, based on the theoretical model of seepage flow and the pressure solution under unit flow rate introduced above (i.e., Equation (7)), the pressure data (pressure drop and pressure derivative) under unit flow rate calculated by deconvolution is analyzed by double logarithmic typical curve fitting method, so as to interpret the reservoir parameters and fracture parameters. At the same time, based on the theoretical model of seepage flow and the flow rate data under unit production pressure drop introduced above (i.e., Equation (8)), the flow rate data under unit production pressure drop calculated by deconvolution are analyzed by Blasingame production decline typical curve fitting method, so as to interpret the reservoir parameters and fracture parameters. The two fitting methods can constrain each other and significantly reduce the uncertainty of model interpretation results. In addition, during the fitting process of the feature curve, the B-spline cardinality and smoothing factor are used as constraints, so that the normalization parameter tuning and the theoretical model calculation parameter tuning are mutually constrained. At the same time, seismic data and on-site data are used as conditional constraints, and through the combined action of multiple constraints, the multiplicity of interpretation results is greatly reduced, and the fitting degree of the double logarithmic typical curve is improved. Therefore, the interpretation results have high accuracy. Finally, using Duhamel’s principle (i.e., Equation (12)) and the model analytical solution (i.e., Equation (8)) to calculate the flow rate under the unit production pressure drop, the daily production rate q and cumulative production rate of horizontal well under any production pressure control can be calculated. Then productivity can be predicted more accurately and efficiently. Based on the model productivity calculation, under different production life, fracturing cost and oil price, the fracture spacing is optimized with the goal of maximum cumulative production and break-even, respectively. The flow chart of fracture spacing optimization method for multi-stage fractured horizontal wells in shale oil is shown in Figure 3.

There are some precautions to apply this calculation method to practical engineering. The optimization method of fracture spacing is only applicable to single-phase flow and cannot occur in oil gas two-phase flow. The pressure of the reservoir cannot be lower than the bubble point pressure. The fracture spacing has little effect on the fracture propagation ability. In addition, obtaining a certain amount of data is necessary.

3. Practical Application

In this section, the fracture spacing of a multi-stage fractured horizontal well in a shale oil block at a China oilfield is optimized. Due to the low water production rate during the long-time production period of this well, fluid flow is considered as single-phase oil flow. And due to the absence of other wells around the well, interference between wells is not considered. Therefore, this well is suitable for the mathematical model in this article. It is known that the initial pressure of reservoir is 12.5 MPa, the length of the horizontal wellbore is 1740 m, the wellbore radius is 0.076 m, the number of main fractures is 60, the width of the main fracture is 0.001 m, the effective reservoir thickness is 9.9 m, the porosity of the shale matrix is 10%, the fluid viscosity of oil is 0.5 mPa∙s, the average water cut of the production well is 0.35, the shale oil density is 850 kg/m³, the fracture cost of hydraulic fracturing is 160,000 Yuan per cluster and the current shale oil price is 3800 Yuan/ton. The matrix permeability is less than 1.0 mD and the volume coefficient is 1.3. The information can be used as constraint for dynamic production data analysis.

3.1. Dynamic Production Data Analysis

The multi-stage fractured horizontal well in Ordos Basin was tested for long-time flowing pressure without well shut-in. The dynamic data of bottom hole pressure and daily production rate are shown in Figure 4.

No other stimulation measures are implemented during the well’s production; furthermore, the seepage flow process can be approximated as a single-phase oil flow. Thus, the production data can meet the requirements for deconvolution application. The production data were analyzed according to the aforementioned analysis method based on deconvolution.

First of all, the dynamic data of variable pressure and variable production rate in Figure 4 can be normalized by the application of deconvolution algorithm and Equation (11). As a result, the deconvolved bottom hole pressure data per unit flow rate are obtained, which is shown in Figure 5. It can be seen from Figure 5 that the application of deconvolution eliminates the impact of data errors and a smooth pressure drop curve is obtained.

The typical curve analysis (pressure drop and pressure drop derivative) can be performed using the unit-rate bottom hole pressure solution (i.e., Equation (7)) of the model obtained. The analysis result is shown in Figure 6. It can be seen from the Figure 5 and Figure 6 that the model obtained fits the normalized production data very well; the application of deconvolution eliminates the impact of data error, and data divergence is effectively prevented. Smooth typical curves are obtained, and the bottom hole pressure drop behavior in the reservoir development can be clearly reflected.

Then, normalizing the dynamic data of variable pressure and variable production rate in Figure 4 by the application of deconvolution algorithm and Equation (12), the deconvolved production rate data per unit production pressure drop is obtained, which is shown in Figure 7. It can be seen from Figure 7 that the impact of data errors is eliminated and a smooth production decline curve is obtained.

The typical curve analysis of production decline under unit production pressure drop can be performed using the flow rate solution (i.e., Equation (8)) under unit production pressure drop of the model obtained. The analysis result is shown in Figure 8. It can be seen from Figure 7 and Figure 8 that the obtained model fits well with the normalized production data; the application of deconvolution eliminates the impact of data errors and effectively prevents data divergence. A smooth typical curve has been obtained, which can clearly reflect the production rate behavior during reservoir development.

The fitting method in Figure 6 is performed from the perspective of pressure drop and pressure drop derivative. The fitting method in Figure 8 is performed from the perspective of production rate. The two fitting methods can constrain each other and significantly reduce the uncertainty of model interpretation results.

Based on Duhamel’s principle and Equation (12), the historical fitting data of production rate are obtained, which are shown in Figure 9. It can be seen from Figure 9 that the fitting effect of productivity history data is very good. The dynamic geological reserve is evaluated as 3.17 × 10⁶ m³.

The seepage flow behavior can be characterized more accurately. The reservoir parameters and fracture parameters are shown in

Table 1. Interpretation results of reservoir parameters and fracturing parameters for a shale oil reservoir.

Parameters	Values Obtained by Typical Curve Analysis	Nuit
Initial reservoir pressure	12.5	MPa
Reservoir thickness	10	M
Half-length of fracture	50	M
Fracture width	0.001	M
Fracture conductivity	30	mD∙cm
Shale matrix permeability	0.75	mD
Fluid viscosity	0.5	mPa∙s
Dynamic geological reserves	3.17 × 106	m3

.

3.2. Optimization of Fracture Spacing for Multi-Stage Fractured Horizontal Wells in Shale Oil

The production pressure drop is set as 7.5 Mpa in the future for the well. Then according to Duhamel’s principle (i.e., Equation (12)) and the flow rate solution under unit production pressure drop calculated by the interpretation model (i.e., Equation (8)), the daily production rate of horizontal well under any production pressure control can be obtained quickly. The cumulative production rate of horizontal well can be obtained through an integral calculation. Then the fracture spacing can be optimized according to the well productivity.

In the following, the fracture spacing are optimized from the aspects of production life, cumulative production, total economic benefits (it is equal to the production multiplied by oil price), balance of payments (the total economic benefit is equal to the fracturing cost.), fracturing cost, oil price and other influencing factors. Research on the variation in daily production rate and cumulative production rate with different fracture spacing (or different fractures number under the same horizontal well length) with production time is conducted. The effect of fracture spacing (or fracture numbers) on cumulative production under different production life is studied. The optimal fracture spacing (or fractures number) is determined with the goal of the maximum cumulative production and the balance of payments. By changing the oil price, the effect of oil price on the optimal fracture spacing is studied, with the goal of maximum total cumulative production and the balance of payments. By changing the fracturing cost, the effect of fracturing cost on the optimal fracture spacing is studied, with the goal of maximum total cumulative production and the balance of payments. Finally, the results of fracture spacing optimization are compared and analyzed.

3.2.1. Production Rate Change with Production Time under Different Fracture Spacing

The variation in daily oil production rate with production time under different fractures number is shown in Figure 10. It can be seen from Figure 10 that daily shale oil production rate declines continuously as production time goes on, with a rapid decline rate in the first year and then a relatively slow decline rate. The reason is that with the continuous production of the reservoir, the reservoir pressure gradient gradually decreases; then, the flow rate of shale oil and the recovery rate of the reservoir decreases. In the early stage of reservoir exploitation, the more the fracture number, the smaller the fracture spacing, the greater the daily shale oil production rate; however, in the later stage of exploitation, the more the fracture number and the smaller daily shale oil production rate. The reason is that in the early stage of production, the more the fracture numbers, the greater the stimulated reservoir volume area that can introduce a higher rate of oil recovery. Therefore, the daily production rate is higher. However, with the increase in production time, the more the hydraulic fractures, the smaller the fracture spacing, the easier the reservoir is to be mined out and the faster the production decline.

The variation in cumulative shale oil production with production time under different fracture numbers is shown in Figure 11. It can be seen from Figure 11 that as production time goes by, the cumulative production of shale oil gradually increases, but the growth rate gradually decreases and the cumulative production growth rate is faster in the early stage of exploitation. This is because the daily production rate of shale oil is large in the early stage of exploitation, and the daily production rate of shale oil gradually decreases with the growth of time. At the same time, the more the fracture number, the higher the cumulative production of shale oil; however, with the increase in the fracture number, the increase in the cumulative production becomes less. This indicates that with the increase in the fracture number (that is, the fracture spacing decreases), the effect of the increase in the fracture number on the increase in cumulative shale oil production will become lower and lower, and the economic benefits brought by the increase in the fracture number will become less and less.

3.2.2. Relationship between Cumulative Production and Fracture Spacing under Different Production Life

The curves of cumulative production and fracture number under two different production lives are shown in Figure 12. The reservoir settings, fracture conductivities and economic parameters under different production life are same. It can be seen from Figure 12 that the cumulative production gradually increases with the increase in the fracture number (the fracture spacing decreases), and the cumulative production basically remains unchanged when the fracture number increases to a certain extent. In the case of 5 years of production, when the fracture number is more than 320, the cumulative shale oil production remains basically unchanged, which can be regarded to be very close to the well-controlled reserve, and the increase in the fracture number cannot bring more economic benefits. In the case of 10 years of production, when the fracture number is more than 160, the cumulative shale oil production remains basically the same, and the increase in the fracture number does not bring any more economic benefits. And the fracture number corresponding to the maximum cumulative production for 5 years of production is greater than that for 10 years of production.

The cumulative production difference between 5 years and 10 years is equal to the vertical distance between the two curves in Figure 12. Therefore, it can be concluded that as the fracture number increases (the fracture spacing decreases), the cumulative production difference decreases gradually. And as the fracture number reaches 320, almost no difference exists. This indicates that when the fracture number is greater than or equal to 320, almost no oil can be produced in the second 5 years of production, which means that the reservoir has been fully exploited in the first 5 years. As the fracture number increases, the 5-year cumulative production increases faster than the 10-year cumulative production. If more oil is needed in the short term, the fracture number can be increased.

3.2.3. Fracture Spacing Optimization under Different Production Life

(1)

Cumulative Production as the Optimization Objective

The curve of the relationship between cumulative production and fracture spacing under different production life is shown in Figure 13. It can be seen that with the increase in fracture spacing, the cumulative production gradually decreases. When the fracture spacing is less than 5.5 m, the 5-year cumulative production stays at the maximum, which means the reservoir is almost fully developed based on 5-year production life. When the fracture spacing is less than 11.4 m, the reservoir is almost fully developed based on 10-year production life. Since the cost of fracturing increases with the fracture number, the fracture number should not be too high, that is, the fracture spacing should not be too small. Therefore, the minimum spacing of fractures should be 5.5 m, which corresponds to a 5-year cumulative production of 7394 tons; and the minimum spacing of fractures should be 11.4 m, which corresponds to a 10-year cumulative production of 7430 tons. The minimum fracture spacing for 5 years is smaller than that for 10 years. If the production period is shorter, the maximum fracture spacing should be smaller.

(2)

Balance of Payments as the Optimization Objective

The curves of relationship between total economic benefit under different production life and fracturing cost and fracture spacing is shown in Figure 14. The total economic benefit is obtained by multiplying the cumulative production with the shale oil price of 3800 Yuan/ton. The higher the cumulative production, the higher the total economic benefit. With the increase in fracture spacing, the total economic benefit gradually decreases. This is because as the fracture spacing increases, the fracture number decreases, the mining rate decreases, the cumulative production decreases and the total economic benefits also decrease. The larger the fracture spacing, the smaller the fracture number and the smaller the total fracturing cost. The fracture cost is selected as 160,000 Yuan per cluster. When the fracture spacing is too small, the fracture number is too high, which will lead to higher production cost than the total economic benefit, that is, a deficit emerges. When the fracture spacing is 10.3 m, the 5-year payments are exactly balanced, and the total economic benefit for 5 years and fracturing cost are both 27.08 million Yuan. When the fracture spacing is 9.8 m, the 10-year payments are exactly balanced, and the total economic benefit for 10 years and fracturing cost are both 28.16 million Yuan. The longer the production life, the smaller the minimum fracture spacing and the higher the total economic benefit and fracturing cost.

The curves of relationship between total economic benefit under different production lives and fracturing cost and fracture spacing at different oil price are shown in Figure 15; the oil price is set as 3300 Yuan/ton, 3800 Yuan/ton and 4300 Yuan/ton. It can be seen that as the oil price rises, economic benefit will also increase. In the equilibrium state of payments, the higher the oil price, the smaller the fracture spacing and the more the produced oil. The fracture spacing under the zero profit constraint for 10 years is smaller than that for 5 years.

The curves of relationship between total economic benefit under different production life and total fracturing cost and fracture spacing at different fracturing cost per cluster are shown in Figure 16. The fracturing cost is set as 140,000 Yuan per cluster, 160,000 Yuan per cluster and 180,000 Yuan per cluster, respectively. The oil price keeps constant. It can be seen that as the fracturing cost increases, the fracture spacing corresponding to the equilibrium state of payments is larger and the produced cumulative oil is less. With the advancement of technology, the fracturing cost will reduce largely and then the fracture number corresponding to the equilibrium state of payments will increase, the fracture spacing will reduce and, thus, more produced cumulative oil will be obtained. The fracture spacing under the zero profit constraint for 10 years is smaller than that for 5 years, but their difference is not big.

(3)

The Optimization Result Analysis

The statistical results of fracture spacing optimization are shown in

Table 2. Statistical results of fracture spacing optimization.

Production Lifetime	Optimal Fracture Spacing Aiming at Maximum Cumulative Production	Maximum Cumulative Production	Fracture Spacing under the Zero Profit Constraint	Cumulative Production under the Zero Profit Constraint	Effect of Oil Price on Fracture Spacing under Zero Profit Constraint	Effect of Fracturing Cost on Fracture Spacing under the Zero Profit Constraint
5 years	5.5 m	7394 tons	10.3 m	7126 tons	For a 500 Yuan increase in oil price per ton, the fracture spacing should be reduced by 1.6 m, and the total oil production will be increased by 140 tons.	For a 20,000 Yuan increase in average fracturing cost per cluster, the fracture spacing should be reduced by 1.5 m, and the total oil production will be increased by 121 tons.
10 years	11.4 m	7430 tons	9.8 m	7410 tons	For a 500 Yuan increase in oil price per ton, the fracture spacing should be reduced by 1.4 m, and the total oil production will be increased by 3 tons.	For a 20,000 Yuan increase in average fracturing cost per cluster, the fracture spacing should be reduced by 1.3 m and the total oil production will be increased by 1.3 tons.

. According to Table 2, in comparison with 5-year production life, the optimal fracture spacing aiming at maximum cumulative production corresponding to 10-year production life is larger. There is little difference between the maximum cumulative production for 5 years and 10 years. The optimal fracture spacing aiming at the highest profit is equal for 5 years and 10 years. The fracture spacing under the zero profit constraint for 10 years is smaller than that for 5 years, but their difference is not big. The cumulative production under the zero profit constraint for 10 years is larger than that for 5 years. When the oil price increases by 500 Yuan per ton, in comparison with the case of 5 years of production, fracturing spacing under the zero profit constraint for the case of 10 years of production needs to be reduced by a smaller value and the total production increase is also smaller. And when the average fracturing cost of each cluster is reduced by 20,000 Yuan, the reduction in 10-year fracture spacing under the zero profit constraint is small. The smaller the production life, the greater the impact of increasing the same oil price or fracturing cost on the optimal fracture spacing, and the higher the sensitivity.

4. Conclusions

In this study, an optimization method for fracture spacing of multi-stage fractured horizontal well based on dynamic production data inversion is proposed by making full use of the abundant production data in shale oil field. A deconvolution algorithm is applied to normalize the production data in order to transform the data of variable pressure and variable flow into the pressure data under the unit flow rate and flow rate data under unit production pressure drop. And the influence of data error can be eliminated.

First, a trilinear seepage flow mathematical model of multi-stage fractured horizontal well in shale oil reservoirs and its pressure solution under the unit rate in Laplace domain are introduced. Then, the pressure solution under the unit flow rate is used to fit the normalized pressure data under the unit flow rate by double logarithmic typical curve fitting method. The flow rate solution under unit production pressure drop is used to fit the normalized data under the unit production pressure drop by a Blasingame production decline typical curve fitting method. And some parameters of reservoir and fracture are interpreted. The two fitting methods can constrain each other and significantly reduce the uncertainty of model interpretation results. The interpreted mathematical model is more in line with the reality and the seepage flow behavior can be depicted more accurately. Furthermore, using the Duhamel’s principle and the rate solution under unit production pressure drop of the interpreted model, the daily production rate and cumulative production of horizontal well under any production pressure regime can be obtained quickly, and, thus, the productivity can be predicted more accurately and efficiently. Based on the productivity calculation of the model, the fracture spacing is optimized from the perspective of productivity.

Finally, based on the above optimization method for fracture spacing, the fracture spacing of a production well in an actual shale oil block is optimized from the aspects of production life, cumulative oil production, total economic benefit, net profit, fracturing cost, oil price and other influencing factors. The research results show that if the maximum cumulative production is taken as the goal, the optimal fracture spacing is 5.5 m for 5 years and 11.4 m for 10 years. The maximum cumulative production for 5 years is 7394 tons. The maximum cumulative production for 10 years is 7430 tons. The optimal fracture spacing for 5 years is less than that for 10 years. If the production period is shorter, the maximum fracture spacing should be smaller. And if the balance of income and expenditure is taken as the constraint, the optimal fracture spacing is 10.3 m for 5 years and 9.8 m for 10 years. The cumulative production under the zero profit constraint for 5 years is 7126 tons. The cumulative production under the zero profit constraint for 10 years is 7410 tons. The shorter the production life, the more sensitive the effect of oil price and fracturing cost towards the selection of the optimal fracture spacing. All in all, a significant reference value for hydraulic fracturing and the optimization of fracture spacing of adjacent wells in the same shale oil block is provided in this study, and some technical guidance for the later production and secondary fracturing of reservoir is also provided.