Analysis of bottom hole pressure of unsteady flow in dual-medium fracturing vertical well

In order to improve the validity of bottom hole pressure model, and simplify its calculation process, a mathematical model of instantaneous pressure for unsteady flow was established by considering the crossflow between the fractures and matrix. Different conditions, including the reservoir top has constant pressure, were considered. The basis for obtaining bottom hole pressure is to solve diffusivity equation with the integration of axisymmetric transformation and similar methods, which is presented for the first time. Different from the traditional method of using the Green’s function and source solution, this paper uses Laplace transformation, axisymmetric transformation and similar methods, separation of variables to obtain the analytical solution of Laplace domain. Then, the Stephenson Numerical method was used to obtain the numerical solution in a real domain. The results of this method agree with the numerical simulations and actual test data, suggesting the validity and accuracy of this method. Finally, the sensitivity analysis revealed that the pressure curve can be divided into eight stages, namely, early linear flow, continuous flow transition section, fracture linear flow, formation linear flow, crossflow, transitional flow, pseudo-radial flow and boundary control flow. The advantage of the analytical solution utilized in this paper is to incorporate exchange coefficient and skin factor efficiently, providing a theoretical basis for optimizing production pressure difference and determining the reasonable productivity.


Introduction
The concept of tight oil can be divided into broad and narrow senses. Tight oil in a broad sense refers to the petroleum resources stored in tight reservoirs, and its development requires fracturing technology. Tight oil in a narrow sense refers to the petroleum resources stored in tight reservoirs other than shale, and its exploitation also requires fracturing technology. The broad concept of tight oil includes oil resources in shale reservoirs and emphasizes the tightness of the reservoirs. This paper mainly focuses on tight oil in a narrow sense.
Currently, the development of tight oil reservoirs largely depends on hydraulic fracturing technology. Through hydraulic fracturing, fractures can be built in vertical wells, which not only change the way of fluid seepage near the wellbore, but also expand the swept area in the process of oil reservoir exploitation and enhance oil recovery. Therefore, hydraulic fracturing technology has broad application prospects in tight oil reservoir development, and how the bottom hole pressure changes along with the fracturing process is of great importance for the tight oil exploitation.
Formation pressure refers to the pressure acting on the fluid in the pore space of the rock in the reservoir formation, which reflects the energy in oil reservoir, and it is the driving force of the crude oil flow. In the oil reservoir development, bottom hole pressure is the most direct response to the formation pressure. Generally, there are two methods for obtaining pressure, one is direct measurement, and the other is indirect acquisition through well test, as shown in Fig. 1.
For the indirect method, many researchers focus on the analytical model. They have done a lot of researches on the unsteady flow pressure of fracturing vertical wells and have made great progress (Cai et al. 2008;Wang et al. 2014;Zhao et al. 2015;Liu 2016;Zhu et al. 2017;Shen et al. 2018;Jiang et al. 2019). However, previous classical methods still have limitations in the application and guidance of tight oil reservoir development. It is mainly showed as follows: (1) The seepage model is ideal and does not fully consider the crossflow between the fracture and matrix. (2) In order to obtain the solution to the pressure, some approximate treatments were made, which affected the accuracy of the results. (3) The derivation process is complicated, and the calculation is inconvenient. (4) Green's function is suitable for closed boundary reservoirs with considering the upper and lower boundaries, so the method has certain limitations (Ramey and Gringarten 1973;Ramey and Gringarten 1974).
We attempt to give some new insights in understanding the bottom hole pressure in a circular reservoir. This paper presents an analytical solution that describes pressure transient behavior for unsteady flow by considering the crossflow between the fracture and matrix, and it is successfully applied to examine effects of the fracture orientation, fracture half-length, elastic reservoir capacity ratio, exchange coefficient, skin factor and permeability ratio on performance of fracturing vertical well in a tight oil reservoir with closed boundaries based on pressure and pressure derivative concepts. More specifically, this paper uses Laplace transformation, axisymmetric transformation and similar methods, separation of variables to obtain the analytical solution of Laplace domain. Then, the Stephenson Numerical method (Stehfest 1970) was used to obtain the numerical solution in a real domain. The results agree with the numerical simulations and actual test data, suggesting this method is scientific and available. The advantage of the solution devotes to its efficiency and accuracy.

Mathematical model
Basing on production testing characteristics and well test analysis, it is recognized that the dual media reservoir has a linear flow after fracturing (Bourdet 1985;Liu 1993), that is, the fracture can stably supply fluid to the wellbore. To simplify the description of the subsurface fluid by the mathematical model, the following assumptions are made as below: (1) The fluid is a micro-compressible fluid with two media, namely matrix and fracture. (2) The crossflow between the matrix and the fracture is considered. (3) The gravity and capillary force are not taken into consideration. In the center of the reservoir, there is a fractured vertical well, as shown in Fig. 2. The diffusivity equation of the unsteady flow is given as follows:  (1) h m , h f are the effective thickness in the matrix and the fracture, m.
t is the production time, s. q it the fluid production, m 3 /d. c tm , c tf are the compressibility of the matrix and the fracture, MPa −1 . m , f are the porosity of matrix and fracture. S is the skin factor. The subscripts f and m represent fracture and matrix, respectively. Equations (1)-(4) constitute the instantaneous pressure mathematical model of the unsteady flow in a dual-medium fracturing vertical well.

Model solving
At present, there are mainly five methods to solve the unsteady flow seepage control equation. (1) The Green function: Without considering the propagation time of the seepage pressure, it is used to obtain the analytical solution of the diffusivity equation in the Laplace space domain (Feng and Ge 1985;Alpheus and Tiab 2008). (2) The steady-state successive replacement method (SSR): The pressure propagation radius is subdivided into multiple segments, and the pressure in each segment can be expressed as a function of time. SSR is used to obtain the numerical solution of the diffusion equation (Li and Liu 1997;Buhidmal and Raghavan 1980). (3) The numerical approximation algorithm: Basing on the series theory, it is used to obtain the numerical solution of the diffusion equation (Feng and Fu 1997;Liu et al. 2015). (4) The perturbation method: It is used to obtain the approximate solution of the diffusion equation in the case of constant production (Kale and Mattar 1980). (5) The Boltzmann transformation method: It is used to obtain the solution of the diffusion equation under the condition of constant production in an infinite reservoir (Peres et al. 1987).
The above five methods have a good effect on solving the self-similarity problem. However, there are still some limitations in solving the nonlinear seepage diffusion equation. The basis for obtaining bottom hole pressure is to solve diffusivity equation with the integration of axisymmetric transformation and similar methods, which is presented for the first time. Different from the traditional method of using the Green's function and source solution, this paper uses Laplace transformation, axisymmetric transformation and similar methods, separation of variables to obtain the analytical solution of Laplace domain. Then, the Stephenson Numerical method was used to obtain the numerical solution in a real domain.
The Laplace transform of Eqs.
(1)-(4) are According to the solution form of the axisymmetric two-dimensional seepage equation, suppose the solution of Eq. (5) is given as follows: w h e r e : (5), and we can get: Equation (10) has a solution, which must satisfy the following conditions:

Which is
The solution of Eq. (12) is given as follows: Substitute Eq. (13) into Eq. (9), and we can get: where: According to Eq. (14), at the interface between fracture and matrix, the pressure is equal, and the analytical expression of bottom hole pressure can be solved as follows: Equation (15) is the analytical solution of the Laplace domain, and the pressure value in the real domain is obtained by Stehfest numerical inversion method (Stehfest 1970).

E300 simulation
The latest E300 module in the Eclipse2014 software is well applied in dual-medium heterogeneous reservoirs. This article uses this module to simulate bottom hole pressure. First, the oil reservoir is set as a circular reservoir, and there is a vertical well producing in the center of the reservoir. After hydraulic fracturing, 5 fractures are built, and they are all symmetrically distributed about the wellbore. The reservoir boundary has a constant pressure. In order to facilitate the description of the interaction between the fractures and the matrix around the wellbore, corner grids are used for grid division, and each fracture has at least three grids. The reservoir radius is 1000 m. The radial length step is 25 m, and the direction step is 10°. In the vertical direction, 5 simulation layers are planned, so the total number of grids is 40 × 36 × 5 = 7200 grids. The grid division diagram is shown in Fig. 3. The other required parameters are shown in Table 1. The solution of the E300 simulation and this paper is shown in Fig. 4. We can see a very good agreement between the solution in the E300 simulation and this paper, and the relative error is very small, indicating that the method used in this paper produces reliable bottom hole pressure.

The actual data
A new sweet spot was discovered in a tight oil reservoir in 2012 (Liu et al. 2012). After the well in the sweet spot center was hydraulically fractured, 24 fractures were built. Well tests showed that the sweet spot was approximately a circular reservoir. The reservoir radius is about 1932 m, and its effective thickness is 16.5 m. The bottom hole pressure was monitored during oil production. Based on the actual reservoir parameters, the model in this paper is used to calculate the bottom hole pressure, as shown in Fig. 5. Figure 5 shows that the two methods are in good agreement, indicating the validity and accuracy of this method. In summary, both show that the method used in this paper produces reliable bottom hole pressure.

Flow period division
By using Stehfest numerical inversion method to solve Eq. (15), the pressure values in the real domain were obtained, as shown in Fig. 6. Figure 6 shows that pressure curves can be divided into eight stages.
Stage A: Early linear flow. The fluid in the wellbore flows to the wellhead. This stage is mainly affected by the wellbore storage effect, and the curve shows that the pressure and the derivative of the pressure coincide with each other. The slope of this curve section is 1.
Stage B: It is called a continuous flow transition section, and the fluid in the formation begins to flow into the wellbore.   Stage E: Crossflow. The fluid in the matrix begins to flow to the fracture, and the pressure derivative appears sag. This stage is an unsteady flow stage. Due to the interfacial flow between the fracture and the matrix, the seepage velocity of the fluid in the two media is obviously different, and the pressure derivative curve shows the trend of concave downward. The progressive analysis shows that the greater the permeability difference between fracture and matrix, the more obvious the tendency of pressure derivative curve is.
Stage F: Transitional flow. In this stage, the pressure derivative has a curved curve. If the fluidity in the fracturing unaffected area becomes better, the pressure derivative will drop. If the fluidity in the fracturing unaffected area becomes worse, the pressure derivative will rise.
Stage G: Pseudo-radial flow. The flow between the fracturing affected areas and unaffected areas reaches a pseudosteady state. The pressure curve is 1 when the system is in quasi-steady state.
Stage H: The boundary control flow. This stage is related to the boundary type of the reservoir itself. When the boundary type of the reservoir is a constant pressure boundary, from the perspective of energy, the reservoir has sufficient external energy supply, and the pressure can be continuously conducted until the whole system is in quasi-steady state.

Sensitivity analysis
Based on the numerical solution for fracturing vertical well in a circular reservoir presented in the previous part, a sensitivity study for the parameters affecting the pressure and pressure derivative in the model is carried out. The intention of this study is to show the effect of each of these parameters on the dynamic behavior of fracturing vertical well in a circular reservoir. We evaluate the pressure transient solution by varying the values of six parameters including the fracture orientation, fracture half-length, elastic reservoir capacity ratio, exchange coefficient, skin factor and permeability ratio. The other parameters required for sensitivity analysis are shown in Table 1.

Fracture orientation
Fractures formed by hydraulic fracturing are random. Fracture orientation can be simply divided into two categories, that is, the fracture distribution is symmetrical about the wellbore (Fig. 7a) or not symmetrical about the wellbore (Fig. 7b). Figure 8 shows the fracture orientation mainly affects fracture linear flow. When the fracture distribution is symmetrical about the wellbore, the fluid in the reservoir can flow faster and more quickly to the wellbore, and the fracture acts as a highway connecting the reservoir and the wellbore rather than expanding the drainage area, on par with the other fracture orientation type.

Fracture half-length
As the scale of hydraulic fracturing is different, the fracture half-length formed in the reservoir is also different. The  The schematic diagram of fracture orientation Fig. 8 The effect of fracture orientation on pressure curves effect of the fracture half-length on pressure is shown in Fig. 9. Figure 9 shows that the fracture half-length mainly affects the stage A and B. As the fracture half-length decreases, the pressure propagation slows down, and the bottom hole pressure decreases. During the stage A, as the fracture half-length increases, the time of the stage A increases relatively, that is, the larger the fracture half-length is, the larger the seepage cross-sectional area of the fluid is, and the larger the liquid supply amount to the wellbore is, and the longer the stage A lasts under constant production condition. In addition, the fracture half-length only changes the amount of liquid supply, and it does not affect the speed of the liquid supply. Therefore, the pressure curves appear to be nearly parallel in stage A and B. The slope of the pressure curves is 1.

Elastic storage ratio
The essence of the elastic storage ratio is the ratio of the fluid volume in the matrix to the fracture. That is, the smaller the elastic storage ratio is, the smaller the percentage of oil and gas equivalent in the matrix is. The effect of the elastic storage ratio on pressure curves is shown in Fig. 10. Figure 10 shows that the elastic storage capacity ratio mainly affects the stage C, D and E, showing the width and depth of the "concave" curve. The larger the elastic storage ratio is, the less obvious the "concave" of the pressure curves is. When the elastic storage ratio is 0.05, the "step" segment in the pressure curves is almost invisible, which indicates almost no crossflow. The pressure curves approximate the characteristics of a single medium. It is mainly because when the elastic storage ratio tends to 1, the crossflow is not obvious, and the fluid flows to the wellbore through fracture.

Exchange coefficient
The exchange coefficient indicates the strength of fluid seepage between the matrix and the fracture. The larger the exchange coefficient is, the more easily the fluid flow from the matrix to the fracture. The effect of the exchange coefficient on pressure curves is shown in Fig. 11. Figure 11 shows that the exchange coefficient mainly affects the appearance time in stage E. The larger the exchange coefficient is, the earlier the "concave" section of pressure curves appears. Due to the wellbore storage effect, the exchange coefficient has no effect on pressure curves. With oil well producing, the fracture and matrix begin to supply fluid to the wellbore. The larger the exchange coefficient is, the fluid flows from the matrix to the fracture becomes much easier, and the crossflow occurs earlier. When the pressure of the fracture equals to the pressure of the matrix, the crossflow reaches a quasi-steady state, and there is no fluid exchanging between the matrix and fracture. The whole system is in a quasi-steady state. The pressure curves appear to coincide with each other in the later stages.

Skin factor
The skin factor indicates the degree of pollution or modification near the wellbore. When it is more than zero, it Fig. 9 The effect of fracture half-length on pressure curves Fig. 10 The effect of elastic storage ratio on pressure curves Fig. 11 The effect of exchange coefficient on pressure curves indicates the reservoir near the wellbore has been damaged. On the contrary, the reservoir near the wellbore has been improved. The effect of the skin factor on pressure curves is shown in Fig. 12. Figure 12 shows that the skin factor mainly affects the stage C. The larger the skin factor is, the higher the "hump" peak of pressure derivative curve is. It is mainly because that as the skin factor increases, the pollution of the reservoir near the wellbore becomes more serious, and the additional pressure difference is greater. Only under the condition that the bottom hole pressure must be much bigger enough, can the fluid overcome the resistance of the additional pressure difference and successfully flow to the wellbore. However, the skin factor does not affect the speed of the radial flow, indicating that the pressure curves appear parallel to each other.

Permeability ratio
The permeability ratio of fracture to matrix is generally different, and it can be up to 10,000. In this paper, the permeability of the fracture is set to be constant, and only the permeability of the matrix is changed. The effect of permeability ratio on pressure curves is shown in Fig. 13. Figure 13 shows that with the permeability ratio decreasing, the peak of the "hump" of the pressure derivative curve becomes smaller. It is mainly because that the permeability of the matrix increases, the fluidity of the fluid in the reservoir becomes stronger, and the reservoir fluid can flow under a smaller displacement pressure difference. Therefore, on the pressure curve, the bottom hole pressure gradually decreases as the permeability ratio increases. It should be noted that when the permeability ratio is 1, there is no difference in fluidity between fracture and matrix, and the crossflow disappears.

Conclusion
A mathematical model of instantaneous pressure for unsteady flow was established here by considering the crossflow between the fractures and matrix. Different conditions, including the constant pressure condition, were considered. The basis for obtaining bottom hole pressure is to solve diffusivity equation with the integration of axisymmetric transformation and similar methods, which is presented for the first time. Different from the traditional method of using the Green's function and source solution, this paper uses Laplace transformation, axisymmetric transformation and similar methods, separation of variables to obtain the analytical solution of Laplace domain. Then, the Stephenson Numerical method was used to obtain the numerical solution in a real domain. The results of this method agree with the numerical simulations and actual test data, suggesting this method is scientific and available. We also found that the pressure curves can be divided into eight stages, namely, early linear flow, continuous flow transition section, fracture linear flow, formation linear flow, crossflow, transitional flow, pseudo-radial flow and boundary control flow. The fracture orientation and skin factor mainly affect fracture linear flow. The fracture half-length mainly affects the early linear flow. The elastic storage capacity, exchange coefficient and permeability ratio mainly affect the crossflow. The validity and efficiency of this method suggest that it has potential to be used in practical tight oil exploitation.
Funding information This paper is supported by the project "Shunbei fault-solution reservoir description and reserve classification and evaluation".

Fig. 12
The effect of skin factor on pressure curves Fig. 13 The effect of permeability ratio on pressure curves