Semi-analytical solution for bottomhole pressure transient analysis of a hydraulically fractured horizontal well in a fracture-cavity reservoir

This paper study the role of hydraulic fracture properties on the transient bottomhole pressure (BHP) behavior of a horizontal well producing from a tight fracture-cavity reservoir. A combination of point source function, Laplace transformation and Perturbation transformation are used to obtain BHP step by step. Through literature comparison and numerical simulation, the results of BHP have a good consistency, which indicates the proposed method is scientific and reasonable. We divide the fluid flow into five stages, namely the wellbore storage stage, the karst cave fluid flows to the fracture stage, the radial flow stage of karst cave and fracture system, the matrix fluid flows to the fracture stage and the quasi-steady state stage. We come to the conclusion that the number of fractures and fracture direction mainly affect radial flow stage. In contrast, the length of horizontal subsection and skin factor mainly affect the karst cave fluid flows to the fracture stage. The matrix fluid flows to the fracture stage is more obvious when the fracture half-length and the horizontal segment spacing of the horizontal well are small. The study believes special attention should be paid to reforming the formations at both ends of the horizontal well. The advantage of this method is to incorporate well geometry (skin factor) and hydraulic fracturing design (fracture parameters), which is useful for well test interpretation through generating a new set of type curves. What’s more, this new method has the characteristics of easy calculation. The findings of this study can help for better understanding of well test analysis in fracture-cavity reservoir. However, the limitation of this study is that it is only suitable for this situation the horizontal well does not encounter karst caves and the karst caves in the reservoir are connected to the wellbore through fractures.

www.nature.com/scientificreports/ same, and the orientation of fractures is generally treated as simple, such as fractures parallel to each other or perpendicular to the horizontal section of the horizontal well. Fractures are of various orientations and sizes, and skin factor varies from fracture to fracture. Therefore, it is necessary to consider the actual situation of fracture and establish a pressure-solving model. The point source function is an effective tool to solve the fracture heterogeneity. So far, the fracture network system no longer has the characteristics of embedded matrix and fractures at two different scales. This paper presents a new way to obtain the point source function and a new method of "three steps" to obtain BHP of a hydraulically fractured horizontal well under the constant pressure boundary, as shown below: Step 1: Through using the triple media model, the Pedrosa permeability calculation formula is applied to establish the seepage model of the triple media reservoir considering the formation stress sensitivity.
Step 2: By perturbation transform and Laplace transform, the point source function considering stress sensitivity in fracture-cavity is obtained in Laplace space. The point source function in the infinite plate reservoir is obtained by the principle of mirror image and superposition.
Step 3: The method of giving BHP under the constant pressure boundary is established through using the new point source function. This paper proceeds as follows. First, Section "Model establishing" presents the physical model, and the assumptions and corresponding mathematical model is also present in section "Model establishing". Then, section "Model solving" presents the process to solve the model including solution research, dimensionless transformation, equation solving, and discussion part. Section "Model Validation" validate the model in this paper available with the published analytical solution and numerical simulation, and the flow period is divided. Section "Sensitivity analysis" presents sensitivity research on the type curves. At last, section "Summary and Conclusions" presents the main conclusions.

Model establishing
Physical model. The fracture-cavity reservoir develops fractures, and the distribution of fractures is complex, as shown in Fig. 2. The fracture-cavity reservoir has a triple medium property consisting of natural fractures, karst caves, and matrix. In Fig. 2a, at least three vertical wells are required to develop the reserves of five karst caves. In Fig. 2 (b), by hydraulic fracturing, only one horizontal well can develop the reserves of these five caves. Compared with vertical well, horizontal wells have lower cost and control larger drainage area in the exploitation of fracture-cavity reservoir. www.nature.com/scientificreports/ In Fig. 2b, the well is drilled into the matrix of the reservoir, hydraulic fracturing must be used to build the artificial fractures, which can establish the communication channel between karst caves and wellbore. This paper aims to solve the dynamic response of BHP under the condition that the horizontal section of the well is located in the matrix of the tight fracture-cavity reservoir.
For the Shunbei Oilfield, the fracture-cavity body of the reservoir is relatively developed. Its bedrock porosity is minor, and the permeability is low. It is a somewhat typical tight limestone with a very high carbon-oxygen ratio. The fluid flow in the reservoir mainly depends on the communication between fractures and karst caves 29 . Due to the random distribution of fractures and karst caves, it is generally impossible to form a good fracturecavity system, so the heterogeneity of fracture-cavity reservoirs is particularly strong 30 . The Shunbei Oilfield is a particular fracture-cavity reservoir and has some unique features 31 . On one hand, many unexposed karst fracture-cavity reservoirs are developed along the fault zone, on the other hand, the reservoir storage space, storage type, fluid properties and distribution all show diversity and complexity. The horizontal heterogeneity is substantial, but the vertical connectivity is good, and the matrix has no storage and permeability capabilities. Due to the development of faults in the Shunbei oilfield, it is common to set wells in the reservoir matrix and communicate karst caves and natural fractures by hydraulic fracturing in the drilling design.
Fractured-cavity reservoirs have diverse reservoir spaces and complex fluid flow. The actual fluid flow diagram in the reservoir is shown in Fig. 3a. To simplify calculations and facilitate mathematical modeling, the basic flow of the reservoir is simplified as the flow shown in Fig. 3b. Meanwhile, the following assumptions are made: 1. The oil reservoir is composed of the karst cave system, matrix system, and fracture system. The karst cave system is the main storage space, and the fracture system is the main flow channel. All the well production comes from the influx of the fracture system; 2. Considering the permeability sensitivity of natural fracture system, it is assumed that the permeability of matrix system is constant; 3. There is a point source in the oil reservoir. The initial pressure is equal everywhere in the oil reservoir, and the reservoir temperature is constant during the production process of the oil well. 4. The reservoir fluid is slightly compressible, and the influence of gravity and capillary force is not considered.
The porosity of the reservoir and fluid viscosity is constant.
Since the orientation of fractures is not uniform, the number of fractures in the reservoir is assumed to be m, and the horizontal section of the horizontal well is divided into m sections accordingly. The center of each horizontal well segment is parallel to the fracture center along the z-axis. After division, the schematic diagram is shown in Fig. 4. Without considering the channeling between fractures, each fracture can be regarded as the source phase of the corresponding horizontal well segment. Then, considering the fluid flow of a fracture and the related horizontal section, the pressure model under a single fracture can be established. Finally, the pressure distribution of the entire horizontal well can be obtained according to the pressure superposition.
Mathematical model. Based on the above assumptions, the center of a single fracture is located at (x 0 , y 0 , z 0 ), and the governing equation of the fracture is shown as follows 18 :  www.nature.com/scientificreports/ f(x,y,z) is the source (sink) term, and its expression is shown as below.
If p = p i − p , then the change trend of p and p is consistent. The pressure term p in the governing Eq. (1) can be directly replaced with pressure difference term p.  The governing equation of the karst cave system is obtained as the Eq. (6) When considering the stress sensitivity of the natural fracture system, the natural fracture permeability decreases with the decrease of formation pressure. The natural fracture permeability can be expressed as:

Model solving
Solution research. Literature survey found that there are four main methods to solve the Eq. (4). First: assuming that the reservoir pressure reaches the bottom of the well instantaneously, based on the Green's function, the analytical solution of the pressure in the Laplace space domain is obtained 5,8,11 ; Second: a numerical solution to the pressure in the real domain is obtained based on the separation of variables and the numerical difference discrete method 3,12,13 ; Third: Assume that the bottom hole flow reaches a pseudo-steady state, ignore the pressure and time changes, only consider the relationship between pressure and space, and use the successive steady-state method to obtain the pressure distribution 2,9 ; Fourth: get the pressure in the Laplace space through dimensionless change, Laplace transform, and Fourier transform and then obtain the real domain pressure numerical solution through numerical inversion 14,17,18 . This paper considers the triple-media and stress-sensitive characteristics of fracture-cavity reservoirs, applying perturbation transformation, Laplace transforms, mirror image principle and Poisson summation formula, a new point source function for fracture-cavity reservoirs is established. Simultaneously, considering the problem that the conventional source function cannot calculate the formation stress sensitivity, the numerical solution of the BHP of the horizontal well is solved. www.nature.com/scientificreports/ Equation solving. To transform a practical physical problem into a pure mathematical one, we must first transform dimensional quantities (physical quantities) into dimensionless quantity (mathematical quantity), so that all analysis, calculation is applicable to all the scale of measurement. This process is called dimensionless processing. The dimensionless treatment of this paper is shown in Appendix 1. Through dimensionless transformation of Eq. (4), we can get: Laplace transformation and Perturbation transformation are used to solve Eq. (8) and obtain the point source function for the problem under study. The detail of the derivation is shown in Appendix 2.
The pressure distribution calculation formula of instantaneous point source function per unit strength is obtained as follow: Please note that the point source is located at the origin of coordinates.
Discussion. Point source function is not at the origin of coordinates. If the point source function is not at the origin of coordinates, but at (x wD , y wD , z wD ) , first define the following expression: Then the calculation formula of pressure distribution generated by the unit source function is as follows: If: Then, the time-space solution of the second part in Eq. (11) is: By applying the superposition principle, the solution of continuous point source function can be obtained as follows: Based on the property of convolution, Laplace transform of Eq. (13), and we can get: So, the Eq. (11) can be rewritten as: Upper and lower closed boundary. The pressure of the fracture-cavity reservoir is solved based on the mirror image principle. It is assumed that the point source is in a plate fracture-cavity reservoir with closed upper and lower boundaries. The reservoir height is h . The position of the point source is z w . Schematic diagram of closed boundary mirroring principle, as shown in Fig. 5.
After mirror inversion, the position of the point source is: Based on the superposition principle, the point source function of upper and lower closed plate fracturecavity reservoir can be expressed as: When z = 0 is the closed boundary, and z = h is the constant pressure boundary, the solution of the point source function is: The boundary types of fracture-cavity reservoirs are diverse, so we will not discuss them one by one here.
Solution of obtaining BHP. Based on the above three types of point source functions with different upper and lower boundary properties, seepage solution equations can be constructed under different conditions. For the physical model in this article, the length of the horizontal section of the horizontal well is L, and the BHP solving steps are as follows: www.nature.com/scientificreports/ ① Firstly, it is assumed that the length and height of a certain fracture are 2L fi and h wi . Then, according to the boundary type of fracture-cavity reservoir (the above pressure boundary is taken as an example to illustrate), Formula (20) integrates x in the interval of (x w -L f , x w + L f ) , and then integrates z in the interval of(z w -h w /2, z w + h w /2) . The expression of the pressure of a single fracture in the horizontal section of the horizontal well can be obtained. ② According to the superposition principle, the pressure generated by each fracture in the corresponding horizontal section is superimposed to obtain the Laplace space solution of the total pressure. ③ According to the solution of Laplace space, the Stehfest method is used to carry out numerical inversion, and the numerical solution of real space is obtained.

Model validation
Literature comparison and validation. Riley 20 . For the convenience of comparison, it is assumed that the fracture in the physical model of horizontal well in this paper is one. The angle between the fracture and the horizontal section is 90°. The fracture locates in the center of the horizontal wellbore. To be consistent with the conditions in the literature, the influence of karst caves is not considered. Other basic reservoir parameters are shown in Table 1.
The calculation results of Riley et al. and those in this paper are shown in Fig. 6. Figure 6 shows that the results of the two have a good consistency. Riley et al. is a special case of the model in this paper, so the model in this paper is correct.
The numerical simulation validation. For comparison purposes, the number of fractures is assumed to be 5. The relevant parameters of fractures are shown in Table 2. The reservoir is a plate reservoir with upper and lower boundary with constant pressure. By simulating the productivity of horizontal wells with different fracture numbers and fracture locations, the BHP of the corresponding horizontal wells (the simulation solution) was obtained and compared with the pressure obtained from numerical inversion in this paper (this paper solution). The E300 module in Eclipse 2017 is designed for triple-media fracture-cavity reservoirs. To meet the assumptions of the model derived in this paper, the numerical model is set as follows: The width and length of the www.nature.com/scientificreports/ heterogeneous fracture-cavity reservoir are 1000 m, and there is a horizontal well in the center of the reservoir with a horizontal section length of 200 m, as shown in Fig. 7. The triangulation mesh method is adopted to divide meshing and ensure that each fracture has at least three meshing to describe the heterogeneity of formation fluid. If the fracture is short, local mesh encryption should be performed at the fracture. The plane mesh step size is 20 m, and the vertical mesh step size is 5 m. Then the total number of grids is 50*50*16 = 40,000. Other parameters required for numerical simulation are shown in Table 3. The numerical simulation and calculation results in this paper are shown in Table 4. Table 4 shows that the relative error is controlled within 5% under the basic error, which is consistent with the allowable error range 19 , indicating that the method we provide is reliable.   www.nature.com/scientificreports/ Flow characteristics analysis. A fracture-cavity reservoir with a rectangular outer boundary has a horizontal well in the center. The number of fractures is 4. The length of the horizontal section of the horizontal well is 400 m. k x = k y , k z /k x = 100. The wellbore storage factor is 0.001, and the skin factor is 1. Based on the basic parameters in Table 1, the dimensionless pressure and pressure derivative were calculated, as shown in Fig. 8. Figure 8 shows that the fluid flow process of the triple media fracture-cavity reservoir can be divided into five stages, considering fracture orientation and stress sensitivity.
Stage A: The wellbore storage stage. The pressure and pressure derivative curves coincide, and the asymptotic analysis shows that the slope of the curve is about 1. This stage is mainly affected by the wellbore storage effect. Stage B: The karst cave fluid flows to the fracture stage. Generally, the permeability of the karst cave system is greater than that of the fracture system. As the fracture system supplies fluid to the wellbore, the pressure  www.nature.com/scientificreports/ of the fracture system drops. Then the karst cave system first supplies fluid to the fracture to supplement the pressure of the fracture system. A "dent" appears in the curve. Stage C: Radial flow stage of karst cave and fracture system. The slope of the pressure derivative curve is about 0.5 when the fracture fluid supplied by the karst cave system is in equilibrium with the wellbore fluid supplied by the fracture system. This stage is short or infrequent due to the presence of a matrix, which also supplies fracture fluid. Stage D: Matrix fluid flows to the fracture stage. When the fluid in the karst cave and fracture system is extracted, the pressure of the two systems decreases continuously until it is less than the pressure of the matrix system, and the fluid in the matrix begins to flow to the fracture and karst cave under the action of pressure difference. Stage E: Quasi-steady state stage. When the flow between matrix-fracture, karst cave-fracture, and fracture-wellbore is in equilibrium, the system reaches the quasi-steady state. The pressure derivative curve approximates a horizontal line, and the asymptotic analysis shows that the horizontal line is about 0.5.

Sensitivity analysis
According to the method and process of solving BHP in this paper, the sensitivity analysis of the influencing factors of pressure and pressure derivative is carried out based on the principle of the single factor variable. The influences of fracture number, fracture angle, fracture half-length, skin factor, and horizontal well segment length, and horizontal well segment spacing on pressure and pressure derivative are analyzed in detail.
Horizontal segment spacing of horizontal well. When the total length of the reservoir shot by the horizontal well is fixed, the entire drainage area communicating with the fractures outside the wellbore is certain. The influence of horizontal segment spacing on the type curves is shown in Fig. 9. Figure 9 shows that horizontal segment spacing mainly affects stages B and C. The larger the interval between the horizontal segments of the horizontal well, the longer it takes for fractures to supply the wellbore. The more difficult the channeling flow between the karst cave and fracture will be. When the horizontal segment spacing is large, a new "platform" appears on the pressure derivative curve in stage C, with a horizontal line of about 0.25. The asymptotic analysis shows that the new "platform" height depends on the number of segments in the horizontal well. In general, when the number of segments in the horizontal well is n, the corresponding horizontal line of the new "platform" is about 1/n. www.nature.com/scientificreports/ Infield well testing, it isn't easy to have stage E if the test time is not long enough. Therefore, during well test interpretation, the actual permeability can be obtained by multiplying the permeability k explained in stage C by the number of segments of horizontal wells.
Horizontal well segment length. The influence of different horizontal segment lengths on the type curves is shown in Fig. 10. Figure 10 shows that segment length mainly affects stage B. The longer the segment length is, the lower the "hump" is, and the earlier stage B is. However, the shorter the duration of stage B is, and the lower the value of the pressure derivative is and more quickly recovered to 0.5. The larger the reservoir interval opened by horizontal wells, the more favorable the flow of reservoir fluid to the wellbore, and the earlier the well reservoir effect will disappear.
When the total length of the horizontal segment is fixed, the influence of each segment's length and distribution mode on the type curves is shown in Fig. 11. Figure 11 shows that when the section length at both ends of the horizontal well increases, the pressure drop generated becomes smaller, conducive to horizontal well production. Therefore, in the actual perforation process, it is suggested to increase the production section at the toe and foot ends of the horizontal well to increase the supply of fractures to the horizontal wellbore and improve the horizontal wellbore productivity. Skin factor. The influence of the skin factor on the type curves is shown in Fig. 12. Figure 12 shows that the larger the skin factor is, the more serious the wellbore pollution is. What's more, the longer the duration of stage A is, and the later the appearance of stage B is, and the higher the "hump" is. Stage C is covered in varying degrees, and the skin factor affects the radial flow of the fracture.
When the sum of the total skin factor of each segment is constant, the influence of the size and distribution of the skin factor of each segment on the type curves is shown in Fig. 13. Figure 13 shows that the skin factor causes different flow rates in each section, and the non-uniformly distributed skin factor 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, 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.
Fracture number. The influence of the fracture number on the type curves is shown in Fig. 14. Figure 14 shows that the fracture number mainly affects stage C. When there are many fractures, and the fractured-vuggy www.nature.com/scientificreports/ reservoir has a fracture network, stage D will be covered by stage C. To a certain extent, the matrix does not directly supply the karst caves. Instead, the matrix directly supplies the fractures, mainly caused by the development of a fracture network in the reservoir. The more developed the fracture network, the seepage characteristics of the triple-medium fracture-cavity reservoir are more similar to the dual-media. The matrix/cavities in the reservoir directly supply the fractures and then flow to the wellbore. This flow mode is more conducive to pressure transmission and production increase.
Fracture direction. When the angle between the fracture and the horizontal wellbore changes, different fracture angles have different effects on the type curves, as shown in Fig. 15. Figure 15 shows that the fracture direction mainly affects stage C. When the fracture is perpendicular to the wellbore, the pressure drop is smaller, and the pressure transmission is faster, which is more conducive to improving the productivity of horizontal wells. Because the drainage area is larger, and more fluid in the reservoir flows through the fractures to the wellbore at the same time. The perforating gun should be parallel to the horizontal wellbore during fracturing as much as possible, and then fractures perpendicular to the wellbore should be generated.
Fracture half-length. The effect of the fracture half-length on the type curves is shown in Fig. 16. Figure 16 shows that the fracture half-length mainly affects stage B and stage C. When the fracture half-length is longer, the pressure drop is smaller, which is most conducive to production. When the fracture half-length is very small (≤ 50 m), stage D is more obvious because the fracture half-length is small, the fluid in the matrix cannot flow directly to the fracture. It must flow to the karst cave first and then flows to the fracture through the karst cave. The drainage area in the reservoir is limited, and the pressure spreads more slowly. Therefore, in the design of hydraulic fracturing, the fracture half-length should be greater than 50 m as far as possible.
Permeability modulus. The effect of the permeability modulus (α) on the type curves is shown in Fig. 17. Figure 17 shows that the permeability modulus mainly affects stage B and stage C. The larger the permeability modulus is, the greater the minimum value of the curve in stage B is. The permeability modulus also affects the maximum value of stage C. The smaller the permeability modulus is, the greater the maximum value of the curve in stage C is. Permeability modulus has similar properties to rock compressibility. Since the compressibility is of very small order of magnitude (generally in the order of 10 -4 ), it can be regarded as a constant in engineering applications. However, the order of magnitude of permeability modulus is much larger than that of compression coefficient, and its variation should be considered according to the actual situation in application, which is not simply regarded as a constant in fracture-cavity reservoirs. www.nature.com/scientificreports/ Well test interpretation. Through conventional linear analysis of well test, wellbore storage coefficient can be obtained from the curve of early pure wellbore storage stage, and formation coefficient, formation effective permeability, skin coefficient and reservoir pressure can be obtained from the curve of late radial flow stage. Some initial values obtained by conventional well test interpretation are substituted into the well test model, and the initial parameters are constantly modified through using the genetic simulated annealing algorithm, until the theoretical curve completely fits the measured curve. The whole well test interpretation process is a semiautomatic analysis process. For the triple media reservoir model with fracture and wellbore connection, the parameters to be explained in well test are mainly wellbore storage coefficient, skin coefficient, fracture permeability, channeling-flow coefficient between matrix and fracture, channeling-flow coefficient between karst cave and fracture, elastic storage ratio of matrix and elastic storage ratio of fracture. Genetic Algorithm (GA) is a computational model of biological evolution that simulates the natural selection and Genetic mechanism of Darwin's biological evolution 32 . 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 33 . This paper adopts GA in curve fitting, and its process is shown in Fig. 18.
Well S1 is located in the southern part of the Shunbei field, as shown in Fig. 19. The well was drilled on November 12, 2016, and completed on January 6, 2017, with a depth of 5703 m. The length of horizontal section is 639.2 m, with a closed azimuth of 301.4°.
After the well test, well S1 was shut in and pressure recovered was measured on January 17, 2017. The bottom hole pressure was 61.83 MPa/5703 m when shut in and 73.15 MPa after shut in 110 h. The well test model fitting is shown in Fig. 20. Figure 20 shows that the wellbore storage coefficient of well S1 has changed from small to large, resulting in poor fitting of the pressure derivative curve at the beginning. At the same time, the pressure curve does not coincide with the pressure derivative curve, and the slope of the curve is less than 1. In the middle and late stages of percolation, a "concave" appeared in the pressure derivative curve of well S1, reflecting the supply of liquid from karst caves to fractures. Although the well test data is incomplete, the actual value is in good agreement with the calculated value of the model in this paper.
This case confirms that the measured value and the theoretical calculation value are in good agreement. If we ignore the change of the wellbore storage coefficient, the fitting in the early stage will also be better. After considering the fracture orientation of the reservoir, the model in this paper is closer to the actual situation of the reservoir. www.nature.com/scientificreports/ Please note that: The biggest drawback of the model in this paper is that it cannot completely solve the bottom hole pressure problem of the multi-cavity system. When the horizontal well does not encounter karst caves, the karst caves in the reservoir are connected to the wellbore through fractures. The model in this paper can solve this situation. However, when a horizontal well is drilled into multiple karst caves, because the karst caves directly supply fluid to the horizontal wellbore, instead of first flowing to the fractures and then from the fractures to the wellbore, it will interfere with the bottom hole pressure response characteristics. In this case, the model in this paper needs to be further improved, and a triple-media seepage model needs to be established to solve this problem.

Summary and conclusions
A model is derived mathematically to investigate the role of hydraulic fracture properties on the transient bottomhole pressure (BHP) behavior of a horizontal well producing from a tight fracture-cavity reservoir. First, a new point source function for triple-medium fractured-vuggy reservoirs is established. Then, being based on the new point source function, the steps for solving the transient bottomhole pressure (BHP) of horizontal wells are given. Through literature comparison and numerical simulation methods, the rationality of the method in this paper is verified. By drawing typical pressure and pressure derivative curves, the seepage laws of horizontal wells in a tight fracture-cavity reservoir are summarized, and the characteristics of each flow stage are analyzed in detail.
The results identify five flow regimes, namely the wellbore storage stage (p wD ∼ t D and dp wD /dlnt D ∼ t D ), the karst cave fluid flows to the fracture stage, the radial flow stage of karst cave and fracture system, the matrix fluid flows to the fracture stage and the quasi-steady state stage (p wD ∼ t D and dp wD /dlnt D ∼ t D ).
The presented results reveal that the influences of fracture number, fracture angle, fracture half-length, skin factor, horizontal well segment length, permeability modulus and horizontal well segment spacing play significant roles on transient bottomhole pressure (BHP) behavior of a horizontal well producing from a tight fracturecavity reservoir.
The number of fractures and fracture direction mainly affect radial flow stage. In contrast, the length of horizontal subsection and skin factor mainly affect the karst cave fluid flows to the fracture stage. The matrix fluid flows to the fracture stage is more obvious when the fracture half-length and the horizontal segment spacing of the horizontal well are small.
The study indicates that special attention should be paid to reforming the formations at both ends of the horizontal well to eliminate pollution near the wellbore When using horizontal wells to implement reservoir production stimulation measures. When designing a fracturing operation, the perforating gun should be made parallel to the horizontal wellbore as much as possible, and the fracture half-length formed should be greater www.nature.com/scientificreports/ than 50 m. When designing the fracturing of horizontal wells, it is necessary to increase the heel and toe production sections of horizontal wells, expand the drainage area, and accelerate the supply of fractures to horizontal wellbore, to improve the productivity of horizontal wells. The developed solution in this paper is useful for well test interpretation through generating a new set of type curves and the accuracy of well test interpretation is improved to a certain extent. It can also be used to analyze the seepage law and explain the reservoir parameters, including fracture number, fracture angle, fracture halflength, skin factor, horizontal well segment length, permeability modulus and horizontal well segment spacing. In other words, the results of this study have practical significance, especially for well testing of hydraulic fractured horizontal wells in fractured-vuggy tight reservoirs.

Appendix 1 Dimensionless transformation
The dimensionless transformation is defined as follows:  www.nature.com/scientificreports/ Based on the dimensionless transformation above, Eq. (4) can be summarized as: where f(r D ) is Matrix system dimensionless equation: (32) f (r D ) = 1 8w xD w yD w zD v D cos α cos β cos χ · δ[(r D − r 0D )] 3 dv D (33) (x 0D − w xD ) y 0D − w yD (z 0D − w zD ) ≤ v D ≤ (x 0D + w xD ) y 0D + w yD (z 0D + w zD )  www.nature.com/scientificreports/ Equation (43) is a strongly nonlinear partial differential equation, and a perturbation transformation is introduced: The derivative of Eq. (45) with respect to r D can be obtained as follows: The derivative of Eq. (46) with respect to t D can be obtained as follows: Substitute Eqs. (46) and (47) into Eq. (43), we can get: In the form of power series, theψ and 1/(1/αψ) can be rewritten as Because of the small permeability modulus ψ , scholars believe that the 0-order perturbation solution can meet the needs of engineering calculation. Take the 0-order perturbation solution in Eqs. (49) and (50), and substitute it into Eq. (48), we can get: If: Then Eq. (51) can be rewritten as: The Equation (53) is a non-homogeneous second-order differential equation. The idea of solving it is: Firstly, find the general solution of the corresponding homogeneous second-order differential equation; secondly, find a special solution of the non-homogeneous second-order differential equation; finally the general solution + special solution is the non-homogeneous solutions of second-order differential equation  www.nature.com/scientificreports/ Then, the solution of instantaneous point source in the Laplace space of fracture-cavity triple medium reservoir is: Received: 21 October 2021; Accepted: 4 May 2022