Dynamic analysis of bottom hole pressure of a fractured horizontal well in fractal ow medium

: This paper addresses the problem of the bottom hole pressure (BHP) behavior of a fractured horizontal well in a fractal flow medium. A detailed step by step procedure for obtaining the semi-analytical solution is present. The proposed method is validated available by comparing previous analytical solution. Sensitivity analyses about the influential parameters including wellbore storage factor, fracture properties, skin coefficient, permeability modulus and horizontal wellbore length on BHP are investigated in detail. The fracture scale (fracture half-length, fracture width and fracture number) mainly affects the channeling flow from the non-fracture zone to the fracture zone. With the increase of the fracture scale, the intensity of the channeling flow decreases and its duration becomes shorter. The fracture fractal index does not affect channeling flow its appearing time and intensity, but its influence on typical curves is mainly reflected in the seepage stage of the fracture zone. The pressure and derivative curves gradually upward warp with the increase of time, and the upward warping amplitude increases with the increase of the fractal index. The advantage of this method is to incorporate well geometry (well storage factor) and reservoir structure (the fractal index). The findings of this study can help for better understanding of well test interpretation of fractal reservoir

are many types of wellbore storage systems, and the corresponding interpretation models are also very complicated. The traditional diffusion equation is applied to describe the BHP behavior, assuming that the homogenous reservoir has equal thickness and uniform original pressure distribution [2] . However, some recent research work has shown the homogeneity assumptions do not cover most situations [3] . Therefore, fractal geometry is used to describe the heterogeneous reservoir, which is called fractal reservoir [4] . Many scholars also use numerical method to solve the BHP in complex heterogeneous reservoirs [5][6] . However, the diffusion process is historically related in a fractal reservoir, and the abnormal diffusion characteristics can't be completely described through fractal models. In some of the literature, the fractal geometry and fractional calculus is introduced to describe fluid flow in fractal flow medium [7] .
In the past few years, how to improve the description of complex reservoirs, especially the fractal reservoirs, has always restricted the corresponding analysis accuracy of reservoir pressure transient behavior. Many scholars have established mathematical models according to different reservoir conditions and made research progress. In porous media, Razminia et al believe that moving particles have the characteristics of random distribution [8] . The diffusion equation to describe this medium is actually equivalent to the fracture distribution of the fractal reservoir [9]. In the fractal reservoir, fluid seepage conditions and reservoir heterogeneity are very important for well production. Metzler et al. first introduced time fractional derivatives to describe the fractal characteristics of the reservoir [10] . Then based on time fractional derivative model, Area et al. introduced a fractal-fractional diffusion equation to present the pressure transient behavior [11] . The fractional derivative technique is widely used to describe the fractal characteristics during the last two decades [12][13] .
Although there are many researches on cross-flow at present, the influence of fracture orientation is neglected. The application of fractal theory makes the characterization of heterogeneous reservoirs, especially fractures, more flexible and closer to reality. However, there are few studies on the application of fractal theory in fracture-vuggy reservoirs and it has not formed a system.
Besides, little research results have been reported at this point concerning interference between fractures in the fractional derivative technique. Well test analysis is actually based on the dynamic response of BHP. The diffusion process of conventional reservoirs is faster than that of fractal reservoirs, which means it is necessary to study more accurate mathematical models to characterize the complexity of fractal reservoirs. Therefore, for fractal reservoirs, it is necessary not only to improve the description accuracy of analysis characteristics, but also to establish a mathematical model to obtain BHP, so as to grasp the changing trend of BHP in real time and analyze the productivity of oil wells. This paper presents a new method to obtain the BHP in a fractured horizontal well by considering fracture orientation and the interference between fractures, which is the first application in fracture-vuggy fractal reservoirs. This paper proceeds as follows. Section 1 presents the assumptions and corresponding mathematical model. Section 2 presents the Galerkin method to solve the model and validate the model in this paper available with the published analytical solution, what's more, the flow period is divided. Section 3 presents a sensitivity research on the type curve. Section 4 presents the main conclusions.

Model establishing
In the development of fractal medium reservoirs, the formation pressure decreases and pore structure changes due to the production of underground oil. At the same time, due to the development of fractures in fractal reservoirs, artificial fractures are formed after hydraulic fracturing, which together with natural fractures constitute a complex fracture network system. In this paper, the fracture network system is called the fracture zone and the non-fracture network is called the non-fracture zone. Because the fractal reservoir is complex, we make the following assumptions before establishing the model: (1) The anisotropy of fracture permeability, skin coefficient and wellbore storage factor are taken into consideration.
(2) The oil well is in constant flow mode. The fractal reservoir is heterogeneous. Before production, the pressure is evenly distributed.
(3) The stress sensitivity of fracture is considered. What's more, we apply the fractal technology to describe complex fracture distribution in the reservoir. The interflow flow from the non-fracture zone to the fracture zone is nonlinear seepage. The fluid in the fractal reservoir is micro compressible with oil phase.
(4) The gravity of the fluid and capillary force can be ignored.
Assuming that the number of fractures is m, correspondingly, we divide the horizontal section of the horizontal well into m sections, as shown in figure 1. We set fracture i locates at (xi0， yi0， zi0). According to the equation of state and the principle of material balance, the mathematical model can be obtained as follow.
where ,, where kh / kv is the horizontal / vertical permeability.
( , , ) f x y z is to describe the interflow flow from the non-fractured zone to the fractured zone is nonlinear seepage, and its expression is: where δ is Dirac Delta function; ,,    is the angle between the fracture and the x, y, z axis, degree .
Since the entire reservoir is a fractal reservoir, when the reservoir pressure changes, it will cause changes in porosity. We assume that the change of porosity conforms to fractal characteristics. The fluid porosity is described by fractal dimension. The fluid porosity has an exponential relationship with the pressure [14][15] , and its expression with considering stress sensitivity is written as follow. Permeability not only depends on the way that the pore spaces are gathered, but also is inseparable from the connectivity between fractures (characterized by the abnormal diffusion coefficient). When the reservoir pressure changes, it will cause changes in fracture state, thereby affecting reservoir permeability. Predecessors have done a lot of well test research on sensitive reservoirs, and we find that there is also an exponential relationship between permeability and bedrock overlying pressure. On the basis of the definition of fractal porosity, the abnormal diffusion coefficient can be introduced to obtain the fractal permeability considering the stress sensitivity [16][17] . The expression is written as follows.
Substitute equations (4) and (5) into equation (1), we can get: ; k Ix , k Iy , k Iz are the permeability index in ,, x y z direction.
The initial condition: The inner boundary condition: The outer boundary condition: The equation (6) is the seepage control equation of a horizontal well in a fractal reservoir, including the usage of δ function to describe the channeling flow between the fracture zone and the non-fracture zone, and the application of fractal theory to describe the stress sensitivity and fractal characteristics of the fracture system.

Model solving and validation 2.1 Model solving
The non-dimensional transformation method is used to transform seepage governing equation into dimensionless equation [18] . Through dimensionless transformation, the equation variable is reduced, making seepage governing equation with symmetry and easy to use conventional methods to solve. The dimensionless transformation is introduced, as shown below.
where L is the horizontal wellbore section length. This paper aims at the fractal nonlinear seepage of a horizontal well in a fractured-cavity carbonate reservoir. We use the implicit method to solve seepage governing equation, and use the Galerkin method to obtain the finite element equation of the seepage system. After a series of assumptions, the dimensionless dimensional seepage control equation of a horizontal well in a fractal reservoir is: Initial condition: Inner boundary condition: Outer boundary condition: By integrating the seepage governing equation (16) in a single element, we can get: (22) where Ni is the component of the shape function.
1, 2,3...,  in . n is the cell nodes number. By applying Green's formula and integration by parts to Equation (22), we can get: On the boundary of the cell region  e , we can get: where ni is the cosine in the outer normal direction. For the internal elements and external boundary elements with closed conditions, the boundary integral term is 0, so the equation (23) can be rewritten as: We can write equation (25) in matrix form, that is: The matrix equation (26) can be further simplified as below: The previous work on horizontal well seepage law has indicated that it is very complicated to do research on the inner boundary condition. In most studies of analytical method to obtain pressure, attention has been given to simplify horizontal well as line source or plate source. Previous research has shown that Peaceman formula is generally applied to deal with the coupling relationship between well and grid in numerical well test of non-finite element method. However, the equivalent radius of Peaceman formula is related to the mesh size and is not applicable to the finite element method. Therefore, in this paper, the treatment method of the inner boundary in the finite element method is explored.
The inner boundary conditions of horizontal well test can be divided into two categories according to the calculation methods: uniform flow model and infinite conductivity model, as shown in figure 2. The uniform flow model means that the flow of a horizontal well is evenly distributed along the wellbore and the pressure of the well is unevenly distributed on the horizontal section. The infinite conductivity model refers to the fluid flow in the wellbore without pressure loss, and the pressure is evenly distributed on the horizontal section, but the flow is unevenly distributed on the horizontal section.
(a) uniform flow model (b)infinite diversion model Figure 2 Schematic diagram of inner boundary model Since the pressure loss of fluid flowing in the horizontal wellbore is very small, the assumption of infinite conductivity is closer to the actual situation, but its solution is relatively complicated. Therefore, the BHP is evaluated by the equivalent pressure points of uniform flow model and infinite conductivity model. Previous research [19] has shown that when the uniform flow model is applied to calculate the BHP, the equivalent pressure point is 0.65 ~ 0.71L (L is the length of horizontal section). In this paper, we take 0.68L as the equivalent pressure point. This paper first obtains the dimensionless pressure drop, then the pressure drop value is transformed into Laplace space, and the wellbore storage factor and skin coefficient are directly added into Laplace space. Finally, the dimensionless BHP is inversion into the real space to obtain the BHP of the horizontal well. The discrete numerical Laplace transform method is shown as below.
The solution of BHP in Laplace space is obtained through Duhamel principle, as shown below.
where S is the skin coefficient; CD is the wellbore storage factor; s is the Laplace variable.
The dimensionless BHP of horizontal wells is obtained through Stehfest numerical inversion [20] .
where I is an even number. Trial and error method can be used to determine I, usually I is 2~12.

Model validation
Bourdet, Ayoub and Pirard proposed an analytical solution for horizontal fractured wells with infinite conductivity [21] . This analytical solution is suitable for horizontal fractured wells with constant rate or pressure in anisotropic reservoirs. To verify this method, a special case in which the nonlinear coefficient and permeability modulus are zero, that is 0, 0 . Before the pressure wave reaches the reservoir boundary, the solution to horizontal fractured wells in this paper can be an infinite conductivity situation. From the literature of Bourdet, Ayoub and Pirard, we get some data (table 1). The results in the case of the traditional Darcy seepage model of infinite conductivity are shown in figure 3. Figure 3 shows that the working solution is well matched with Bourdet, Ayoub and Pirard solution, indicating that this method is scientific.

Figure3 Comparison chart of calculation results 2.3 Flow period division
In fractal media, due to the different physical properties of the fracture zone and the non-fracture zone, the flow rate or flow characteristics of the fluid flowing into the two zones are different. In the process of well test analysis, it is generally necessary to select the corresponding well test analysis model according to the flow characteristics of the actual data. Therefore, the division and description of flow characteristics are very important for well test interpretation. In this paper, we use equation (13) Stage A: Wellbore storage stage. The pressure and its derivative curve are coincident. The asymptotic analysis shows that the slope of the curve is 1. This stage is mainly affected by the wellbore storage effect.
Stage B: Skin effect transition stage. The pressure derivative curve presents "hump-shaped". The the rising speed of the pressure curve slows down. This stage is mainly affected by the comprehensive impact of skin coefficient and wellbore storage factor.
Stage C: Early radial flow stage. The pressure wave propagates to the radial flow of the horizontal axis before the external boundary. The pressure curve continues to rise, and the derivative curve is approximately straight line and its value is related to LD.
Stage D: Early linear flow stage. It reflects the linear flow characteristics around the fracture, and the pressure derivative curve shows an upward straight line with a slope of 0.5. The fractures along the horizontal well just like a sink. The fractures enforce the fluid in the reservoir to flow to the fractures with hyperbolic shape first (which lasts a very short time), Then the fluid flows in a linear flow to the fractures.
Stage E: Middle radial flow stage. This stage is observed after stage D, which lasts a very short time and is not easy to appear. The fluid in the fractures flows radial to the wellbore. During the stage E, the pressure curve and pressure derivative curve are straight lines, and the two are approximately parallel in a fractal medium.
Stage F: Quasi-steady state stage. When the flow between fracture zone and non-fracture zone is in equilibrium, the system reaches the quasi-steady state. The pressure derivative curve approximates a horizontal line, with a slope of 0.5. Please note that: The stage F of the fluid flow process in the reservoir would only be observed if the well test lasts a long period and the reservoir is infinite-acting.

Figure 4 The flow division schematic of type curves 3 Sensitivity analysis and application
Based on the semi-analytical solution for a fractured horizontal well in a fractal reservoir presented in the previous part, a sensitivity research on the parameters affecting the pressure and pressure derivative is carried out in detail. This section is to show the effect of these parameters including wellbore storage factor, fracture properties, skin coefficient, permeability modulus and horizontal wellbore length on the pressure-transient behavior of a fractured horizontal well in a fractal reservoir. The pressure and pressure derivative have different shapes for each combination of wellbore storage factor, fracture properties, skin coefficient, permeability modulus and horizontal wellbore length.

Wellbore storage factor
The effect of wellbore storage factor on type curves is present in figure 5. Figure  5 shows that the greater the wellbore storage factor is, the upward translation of the pressure curve goes, and the higher the "hump" of the pressure derivative curve is. When the wellbore storage factor decreases, the early radial flow stage lasts longer, and the linear flow appears later. When the well storage effect is serious, the early radial flow stage is covered up and the linear flow appears directly. The change of the well storage factor affects the duration of the wellbore storage stage, resulting in the corresponding translation of the type curves, however, when the tD>10 2 , the pressure/pressure derivative curve coincides.

Figure 5 The effect of wellbore storage factor on type curves 3.2 Fracture properties
The effect of fracture properties including fractal index, fracture number, fracture width and fracture length on type curves is present in figure 6-9.
Fractal index: Figure 6 shows it mainly affects the early radial flow stage and early linear flow stage. The pressure and its derivative gradually increase with the increase of time. The increasing amplitude increases with the increase of the fractal index. The fractal index has nothing to do with the time when the radial flow of the fracture system occurs, so the depth and duration of the upper and lower concave sections of the pressure derivative curve remain unchanged.
Fracture number: Figure 7 shows it mainly affects the early radial flow stage and early linear flow stage. When there are many fractures, they can form a fracture network, the early radial flow stage will be covered by early linear flow stage. The more developed the fracture network is, the seepage characteristics of the fractal reservoir are more similar to conventional reservoir. The non-fracture area in the reservoir directly supply fluid to the fractures and then flow to the wellbore. This flow mode is more conducive to pressure transmission and production increase.
Fracture length: Figure 8 shows it mainly affects the early radial flow stage, the early linear flow stage and middle radial flow stage. When the fracture length increases, the pressure drop becomes smaller. When the fracture length is very small (≤60m), the middle radial flow stage is more obvious. This is mainly because the fluid in non-fracture area cannot flow directly to the wellbore. It has to flow to the fracture system and then flows to the wellbore (Kambiz, et al., 2016). The drainage area is too limited for much fluid to flow. Therefore, in the design of hydraulic fracturing, the fracture length should be greater than 60m.
Fracture width: Figure 9 shows it mainly affects middle radial flow stage. With the values of the fracture width at the discontinuity increasing, the derivative pressure increases. The increase is mainly caused by the greater drainage area at the discontinuity. According to the results obtained, the effects of the fracture length at the discontinuity is much bigger than the effects of the fracture width.

The skin coefficient
The effect of the skin coefficient on type curves is present in figure 10. Figure 10 shows that with the skin coefficient increasing, the wellbore pollution is more serious, and the additional pressure drop around the wellbore is much bigger. What's more, with the values of the skin coefficient at the discontinuity increasing, the indicated hump in pressure derivative increases, which is mainly caused by much bigger pressure drop [22] . The skin coefficient affects the radial flow in the fracture zone.
When the sum of the total skin coefficient is constant (stotal=12), the influence of the skin coefficient distribution on type curves is shown in figure 11. Figure 11 shows that the skin coefficient causes different flow rates in each section, and the non-uniformly distributed skin coefficient produces different pressure responses at the bottom of the well. When the pollution at both ends of the horizontal well is relatively serious, the pressure drop loss caused is greater [23] , which is not conducive to improving the production of the horizontal well. Therefore, when using horizontal wells to implement reservoir stimulation measures, special attention should be paid to modifying the formations at both ends of the horizontal well to eliminate pollution near the wellbore [24] .

Figure 10
The effect of skin coefficient on type curves Figure 11 The effect of skin coefficient distribution on type curves

Permeability modulus
The effect of the permeability modulus on type curves is present in figure 12. Figure 12 shows that the permeability modulus mainly affects the early radial flow stage and the early linear flow stage. With the values of permeability modulus at the discontinuity decreasing, the indicated hump in derivative pressure increases first, and then the indicated concave also increase. The increase is mainly caused by the greater fluid fluidity enhancement at the discontinuity [24] .
According to the above, the effect of the permeability modulus is very similar to that of the fracture length. Therefore, it is very important to add the skin coefficient to the model description, which can help the model in this paper to avoid incorrect interpretation of well test parameters and obtain accurate well test parameter estimates.

Figure 12 The effect of permeability modulus on type curves 3.5 Horizontal wellbore length
The effect of the horizontal wellbore length on type curves is present in figure 13. Figure 13 shows that the horizontal wellbore length mainly affects the early radial flow stage, the early linear flow stage and the middle radial flow stage. When the horizontal wellbore length increases, the pressure and its derivative curve go much higher, and the early radial flow stage becomes more obvious, which cause the subsequent flow stage appears later. When the horizontal wellbore length is small to a certain extent, the early radial flow stage is easily concealed by the wellbore storage effect. Considering the characteristics of reservoir fractal and stress sensitivity, the pressure dynamic curves of horizontal wells do not converge together for a long period of time, but are distributed in parallel.

Figure 13
The effect of horizontal wellbore length on type curves 3.6 Application Genetic Algorithm (GA) is a computational model of biological evolution that simulates the natural selection and Genetic mechanism of Darwin's biological evolution. It is a method to search for the optimal solution by simulating the natural evolution process. Its essence is an efficient, parallel and global search method, which can automatically acquire and accumulate knowledge about search space in the process of search, and control the search process adaptively to get the best solution [25] . This paper adopts GA in curve fitting, and its process is shown in figure 14.

Figure 14 Curve fitting flow chart
Well A is located in block B03 of W oilfield. The middle depth of oil layer is 6950m, and the effective thickness of formation is 37 m. The porosity is 11.02%, and the comprehensive compression coefficient is 1.319 × 10 -3 MPa -1 . The crude oil viscosity is 1.07 MPa·s. The volume coefficient is 1.41. The well was shut in on April 7, 2021 and tested with pressure recovery. Before shutting in, the fluid production was 60.24 m 3 /d and the comprehensive water cut was 8.4%. The shut-in lasted 151 hours, and the shut-in pressure was 80.34MPa /6950m, and the pressure recovered to 83.08MPa /6950m. The pressure recovery curve showed a gradual upward trend and the shape was smooth and gentle, as shown in figure 15.
As for stage A, the fitting degree of the derivative curve between the actual