Characteristics of transient pressure performance of horizontal wells in fractured-vuggy tight fractal reservoirs considering nonlinear seepage

Since the classical seepage theory has limitations in characterizing the heterogeneity of fractured-vuggy tight reservoirs, well test interpretation results are not consistent with actual production by far. Based on the nonlinear percolation theory, a new nonlinear seepage equation considering the boundary layer and yield stress was derived to describe the seepage characteristics of dense matrix blocks and the stress sensitivity and fractal features of fracture systems were characterized by applying the fractal theory. Thus, the nonlinear model of a horizontal well in a fractured-vuggy tight fractal reservoir was established naturally. Then the finite element method was applied to solve the bottom hole pressure based on the processing of internal boundary conditions. After solving the model, the seepage characteristics of different models were summarized by analyzing the bottom hole pressure dynamic curves and the sensitivity analysis of multiple parameters such the nonlinear parameter and fractal index were conducted. Finally, the practicality of the model was proved through a field application. The results show that the pressure dynamic curves can be divided into nine flow stages and the increase of the nonlinear parameter will cause the intensity of the cross flow from matrix blocks to the fracture system to decrease. The fractal index is irrelevant to the intensity of the cross flow while it decides the upwarping degree of the curve at the middle and late flow stages. On the basis of the results of the field application, it can be concluded that the model fits well with actual production and the application of this model can improve the accuracy of well test interpretation.

Abstract. Since the classical seepage theory has limitations in characterizing the heterogeneity of fracturedvuggy tight reservoirs, well test interpretation results are not consistent with actual production by far. Based on the nonlinear percolation theory, a new nonlinear seepage equation considering the boundary layer and yield stress was derived to describe the seepage characteristics of dense matrix blocks and the stress sensitivity and fractal features of fracture systems were characterized by applying the fractal theory. Thus, the nonlinear model of a horizontal well in a fractured-vuggy tight fractal reservoir was established naturally. Then the finite element method was applied to solve the bottom hole pressure based on the processing of internal boundary conditions. After solving the model, the seepage characteristics of different models were summarized by analyzing the bottom hole pressure dynamic curves and the sensitivity analysis of multiple parameters such the nonlinear parameter and fractal index were conducted. Finally, the practicality of the model was proved through a field application. The results show that the pressure dynamic curves can be divided into nine flow stages and the increase of the nonlinear parameter will cause the intensity of the cross flow from matrix blocks to the fracture system to decrease. The fractal index is irrelevant to the intensity of the cross flow while it decides the upwarping degree of the curve at the middle and late flow stages. On the basis of the results of the field application, it can be concluded that the model fits well with actual production and the application of this model can improve the accuracy of well test interpretation.

Nomenclature D f , h
Fractal dimension, anomalous diffusion coefficient C mt , C ft , C vt Comprehensive compression coefficient of matrix, fracture and cave system, MPa À1 r w , r e Radius of the wellbore and reservoir, m K fh , K fv Horizontal and vertical permeability of fracture system, mD K m , K v Permeability of matrix and cave system, mD p m , p f , p v Pressure of matrix, fracture and cave system, MPa À1 p D Dimensionless bottom hole pressure without considering wellbore reservoir coefficient and skin factor p wD , p wD Dimensionless bottom hole pressure considering wellbore storage coefficient and skin factor in Laplace space and real space / fw , / w Porosity of fracture system and entire reservoir at the point of the production well K fw , K w Permeability of fracture system and entire reservoir at the point of production well, mD / m , / f , / v Porosity of matrix, fracture and cave system l Fluid viscosity, mPa s C f , a K Porosity compressibility coefficient and permeability modulus of fracture system, MPa À1 r 0 , L 0 Capillary radius and length, m n 0 Number of capillaries per unit crosssectional area d Boundary layer thickness, lm

Introduction
The fractured-vuggy carbonate reservoir has a large number of caves that occur in the interior of the bedrock and its fracture distribution is quite complex. Therefore, the concept of triple-medium is applied to characterize the complex pores of such reservoirs (Ge and Wu, 1982;Wang Y. et al., 2017). Studying the seepage model of fractured-vuggy tight carbonate reservoirs requires full consideration of the nonlinear flow of matrix blocks and the stress sensitivity and fractal characteristics of a fracture system Li et al., 2017;Lin and Myers, 2018;Wang W.D. et al., 2017;Zeng et al., 2011). The classical seepage theory has many limitations in characterizing the seepage mode of such reservoirs. A group of international and domestic academics have done lots of research on the seepage model of fracturedvuggy reservoirs. In order to use simple form and definite physical meaning of mathematical language to solve practical production problems, Yao et al. (2004) conducted the well test interpretation of triple-medium reservoirs considering the variable well storage, innovating the well test theory of triple-medium reservoirs. Camacho-Velazquez et al. (2002) proposed a triple-medium model characterizing a fracture-vuggy carbonate reservoir and analyzed the model curves in detail by considering the solutions of the situation that the cave system connects with the wellbore. Popov et al. (2009) proposed a unified approach to the problem, which is relevant in flow simulations of fractured-vuggy reservoirs, where caves are embedded in a porous rock and are connected via fracture networks on multiple scales by using the Stokes-Brinkman equations on the fine scale. Nie et al. (2011) established a new flow model in a triple media carbonate reservoir which considers the unsteady inter-porosity flow based on the hypothesis that the shape of the cave is spherical. Gomez et al. (2006) applied the global optimization solution to research the well test method of fracture-vuggy reservoirs. Wang and Yi (2018) presented a coupling mathematical model for the fractured horizontal well in a triple media carbonate reservoir by conceptualizing caves as spherical shapes in which the infinite conductivity of the acid fractures was taken into account.
The classical seepage theory is based on Euclidean space and has a certain characteristic scale, so Darcy's law cannot describe the complexity and heterogeneity of the fracturevuggy carbonate reservoir well. It is worth mentioning that applying the fractal theory to petroleum industry is a great innovation, which can effectively characterize the permeability and porosity of heterogeneous reservoirs. Chang and Yortsos (1990) presented a formulation for a fractal fracture network embedded into Euclidean matrix blocks. Single-phase flow in the fractal object was described by an appropriate modification of the diffusivity equation. In addition, the viability of the model in explaining parameters such the transmissibility, storativity and fractal dimension of a complex naturally fractured geothermal reservoir was analyzed through the pressure transient test data (Aprilian et al., 1993). What is more, the fractal theory was applied to real well tests in various fractured reservoirs and the physical meaning of the fractal parameters was presented in the context of well testing in 1995 (Acuna et al., 1995). From then on, seepage models of multiplemedium reservoirs developed rapidly since the application of the fractal theory (Kong et al., 2008;Lu et al., 2017;Razminia et al., 2015;Wang et al., 2015;Xu et al., 2008). In addition to that, Xu et al. (2017) developed a fractal dual-porosity model to study the single-phase fluid flow through fractured porous media, which can capture the statistical characteristics of fractures and shed light on the transport mechanism of fractured porous media compared with empirical formulas for effective permeability.
Except applying the fractal theory to deal with the fracture system, studying the seepage model of fracturedvuggy tight fractal reservoirs should also take the nonlinear flow of matrix blocks into consideration (Chen and Yao, 2017;Xu et al., 2013Xu et al., , 2019. In recent years, models of low-speed and non-Darcy seepage flow developed rapidly (Ezulike and Dehghanpour, 2013;Vijayalakshmi et al., 2009;Wu et al., 2019). When it comes to the mathematical descriptions of a low-speed nonlinear seepage phenomenon, the commonly adopted models are mainly divided into three categories: quasi start-up pressure gradient model, section model and continuous model. Quasi start-up pressure gradient model ignores the bending section of the seepage curve, so the model narrows the scope of fluid flow in the low-permeability reservoirs. Section model does not reflect the nonlinear section of the actual seepage condition, and it is difficult to judge the critical point of linear and nonlinear seepage. Continuous models are all based on the seepage experiment and use different mathematical functions to fit the experimental data, but their physical meanings are not rich enough.
Although the predecessors have carried out a lot of research and summarized the seepage characteristics of fractured-vuggy tight fractal reservoirs, the established seepage models ignore the nonlinear and fractal characteristics, which are too simplified and cannot reflect the heterogeneity characteristics of such reservoirs. This paper derives a new nonlinear seepage equation considering the boundary layer and yield stress to describe the seepage characteristics of dense matrix blocks and considers fractal characteristics of the fracture system, so that a seepage model of the horizontal well in a triple-medium reservoir is established. In the conceptual mathematical model, fractured vuggy rock is considered as a triple-continuum medium, consisting of fractures, rock matrix, and vugs. There is a strong demand for processing the inner boundary in order to obtain the bottom hole pressure solution based on the finite element principle. After comparing the seepage laws of different models and conducting sensitivity analysis of multiple parameters, the model is applied to conduct a well test interpretation of an actual reservoir.
2 Seepage mechanism of fractured-vuggy tight fractal reservoirs 2.1 Nonlinear seepage mechanism of a matrix system 2.1.1 Microscopic mechanism of the nonlinear seepage The research shows that the main reasons of nonlinear seepage flow in low permeability reservoirs are described as below (Jiang et al., 2012;Xu et al., 2019): 1. Low permeability reservoirs have low porosity and permeability, small pores and complex microscopic pore structure. These characteristics cause the fluid to have a micro-scale flow effect on the low-permeability reservoir. 2. The effect of the fluid yield stress and the adsorption of the boundary layer cannot be ignored when the size of pore throat is tiny, belonging to the micron level. The thickness of boundary layer becomes thinner with the displacement pressure rising, the percentage of flowable fluid increasing, as well as the viscosity of the fluid increasing, so the fluid in the pore needs to overcome greater fluid yield stress. 3. Any kind of the fluid has a certain yield stress, and only when the displacement pressure can overcome the yield stress will the fluid begin to flow. In lowpermeability reservoirs, small pore throats result in a greater yield stress, and non-Newtonian enhancements make the starting pressure gradient more significant.
Yield stress, adsorption of the boundary layer and microscale effects are not alone but mutually affected and coupled. The presence of microscale effects makes the influence of surface force is much larger than the volume force and needs to be regarded as a dominant force. The larger the surface force is, the stronger the reservoir fluid is adsorbed by the boundary constraints. While the presence of the adsorbed boundary layer makes the seepage channel smaller and exacerbates the micro-scale flow effect, resulting in nonlinear seepage strengthening (Cao et al., 2016;Cossio et al., 2013).
This paper combines the micro-seepage mechanism of low-permeability reservoirs, the boundary layer flow and capillary seepage model to derive a new low-speed nonlinear seepage model. It is worth noting that the values of the boundary layer and the starting pressure gradient in lowpermeability reservoirs are based on experimental data. The new model of this paper can deal with this problem flexibly, which the paper will introduce in detail later.

Derivation of a nonlinear seepage new model
The radius of pore throat in dense matrix blocks is tiny, which can be measured by micrometer, so the effect of the fluid yield stress and the adsorption of the boundary layer cannot be ignored (Wu et al., 2017). Considering the adsorbed boundary layer and the fluid yield stress, the modified form of the Hagen-Poiseuille law can be expressed as: where d is the thickness of the boundary layer, reflecting the influence of the boundary layer on seepage flow. s 0 is the fluid yield stress, reflecting the influence of fluid non-Newtonian characteristics on seepage flow. The reservoir core is equivalent to a set of parallel capillary tubes with radius r 0 embedded in solids. If there are n 0 such capillaries per unit of cross-sectional area, the specific discharge through the porous media block is: In equation (2): In the same capillary tube, the larger the pressure gradient is, the thinner the boundary layer is. Therefore, it is assumed that d/r 0 = a 1 /rp. If the yield stress value of the same fluid remains unchanged, then 8s 0 /3r 0 can be regarded as a constant value, namely 8s 0 /3r 0 = a 2 (the values of a 1 and a 2 are obtained by experimental fitting) (Celata et al., 2006;Xu and Yue, 2007). Equation (2) can be converted into: Since higher order term can be ignored and the formula can be simplified as: c 1 , c 2 , and c 3 are applied to denote the different polynomials of a 1 , a 2 ðc 1 ¼ 4a 1 þ a 2 , c 2 = a 1 , c 3 ¼ 6a 1 2 þ 3a 1 a 2 ). In other words, the value of c 1 , c 2 , and c 3 can only be obtained by experimental fitting. In order to simplify the formula, c 3 can be represented by c 1 and c 2 through exploring the relationship between c 1 , c 2 , and c 3 . Therefore, the new model can be simplified as: Substituting the condition "rp = 0, v = 0" into equation (5), then c 3 can be represented by c 1 and c 2 (c 3 = c 1 c 2 ). Therefore, the new model can be simplified as: In equation (5), the parameter c 1 reflects the yield stress of the fluid and the influence of the boundary layer while parameter c 2 mainly reflects the influence of boundary layer on the seepage flow. It means c 1 /(rpc 2 ) is equivalent to the starting pressure gradient term. When c 1 = 0, the above model is simplified as the Darcy model. But if c 2 = 0, the above model is simplified as the quasi-starting pressure gradient model. When v ! 0, then rp ! c 1 + c 2 , so the real minimum starting pressure gradient is rp min = c 1 + c 2 .
When c 1 + c 2 = 0, the rp min = 0, if c 1 + c 2 > 0, then rp min > 0. When there is a true starting pressure gradient in the seepage, the seepage curve does not pass through the origin and the model is expressed as: Under this circumstance, the seepage curve does not pass through the origin point. As the pressure gradient increases, the curve tends to be straight because when the pressure gradient increases to a certain extent, all pore throats participate in the flow. The thickness of the boundary layer no longer changes, and the quasi-linear flow occurs.

Fractal characteristics of the fracture system
Fractal theory thinks that if fractures have fractal characteristics in fractured-vuggy reservoirs, the permeability or porosity is not only a function of pressure, but also relates to the fractal dimension and anomalous diffusion coefficient. D f denotes fractal dimension which is a symbol of the complexity of a fractal system and it can reflect the geometric characteristics of the reservoirs. h is anomalous diffusion coefficient which depicts the connectivity and geometric characteristics of a fractal fracture network. It is reasonable to characterize such reservoirs by applying the fractal theory, which can accurately reflect the heterogeneity of the fracture morphology so that the fractal structure can be more in line with the actual reservoir. The porosity of the reservoir is exponentially related to the effective overburden pressure, so the expression of the porosity with stress sensitivity can be obtained: The fracture porosity is mainly affected by the aggregation mode of the pore space (characterized by the fractal dimension) and exponentially related to the effective overburden pressure. As for fracture permeability, excepting the aggregation mode of the pore space, the connectivity between fractures (characterized by anomalous diffusion coefficients) also plays an important role (Karimpouli and Tahmasebi, 2016). Pedrosa (1986) found that the permeability and the overburden pressure of the bedrock are exponentially related in the study of unstable test wells in stress-sensitive formations. Fractal permeability considering stress sensitivity can be obtained by introducing the anomalous diffusion coefficient based on the definition of fractal porosity: If we consider the fractal of the fracture in time, it will pose a greater challenge to the calculation of the nonlinear model of horizontal wells in fractured-vuggy tight fractal reservoirs (Albinali et al., 2016;de Swaan, 2016;Noetinger et al., 2001Noetinger et al., , 2016. The reason why we apply fractal theory to deal with fracture system is to study the initial distribution of fractures in space, so we considered the fractal characteristics of the fracture in space, but not in time.
3. The model considers the skin factor, wellbore storage coefficient and anisotropy of the permeability. The horizontal and vertical permeability of the fracture system are K fh and K fv , respectively. The permeability of the matrix and cave system are K m and K v respectively. 4. The horizontal well production is Q, located at the center of the reservoir plane. 5. The influence of gravity and the capillary force was ignored. 6. The stress sensitivity, fractal characteristics of the fracture system and the nonlinear seepage mechanism of the matrix were considered.
What we study in this paper is a carbonated salt reservoir which is well-developed fractures and vugs, the vug in full contact with the fracture is more in line with the actual situation. Lots of current researches on triple medium are based on this model which vugs are in full contact with fractures, ignoring the case that vugs are with limited contact with fracture or within blocks (Kang et al., 2006;Wu et al, 2007). In this study, we focus on singlephase transient flow and fractured vuggy rock, conceptualized also as a triple-or multiple-continuum medium, consisting of highly permeable fractures, low-permeability rock matrix, and vugs.

Mathematical model
In a Cartesian coordinate system, the mass conservation equation of single-phase fluid flow in the model can be combined with the nonlinear motion equation and state equation to obtain a nonlinear seepage model of the horizontal well in a fractured-vuggy carbonate fractal reservoir: where, k vf is the cross flow factor of cave system to fracture system; k mv is the cross flow factor of matrix system to cave system; k mf is the cross flow factor of matrix system to fracture system; r vf , r mv and r mf are the interporosity flow shape factors. The shape factor for m-v or m-f is defined by Warren and Root (1963): For v-f interaction, the shape factor for vugs is defined as: where A vf is the total fracture and vug connection area per unit volume of rock (m 3 /m 3 ) and l vf is characteristic length. If (p AD À p mD ) < (c 1D + c 2D )/2, d LVDA = 0, else d LVDA = 1 À c 1D /[2(p AD À p mD ) À c 2D ] (A = f, v) The initial condition: The closed outer boundary condition: The constant pressure boundary condition: The inner boundary condition: Each dimensionless quantity is defined as follows:

Model solving
Analytic solutions are highly demanding on mathematics and rely heavily on idealized assumptions which are actually very different from the actual engineering problems, and once the engineering problem is slightly more complicated, we cannot directly get the analytical solution, or the analytical solution error is too large.
The finite element method is used to solve the problem of Partial Differential Equations (PDE) with specific boundary conditions in mechanics and mathematics. It requires a large scale of equations, resulting in a large amount of input data and preparation for calculation. The float of the accuracy based on the method is relatively large. The finer the unit division, the more accurate the approximation result, accompanied by a large increase in the amount of calculation. However, the finite element method can simulate infinite complexes with finite, interrelated elements, no matter how complex the geometry can be simplified with the corresponding elements. In that way, results could be calculated and complex engineering problems could be simplified through model analysis.
The process of finite element solution can be described as: total structure discretizationunit mechanics analysis unit assemblytotal structure analysisapplication of boundary conditionstotal structure responsereaction analysis of a unit inside the structure. The fractal nonlinear seepage model established in the paper is complex and suitable for numerical solution rather than analytical solution.
Considering that the infinite conductivity model is more difficult to solve than the uniform traffic model, this paper selects 0.7 L as the equivalent pressure point under the two models, so that the uniform traffic model can be used to evaluate the bottom hole pressure.
The horizontal well is considered as the source-sink term of the unit and being integrated: The horizontal well is divided into n nodes and the finite element formula of the element source-sink term can be obtained by applying the Delta function: The dimensionless form of the formula (17) is: For the fracture system, the finite element equation of the fracture system can be obtained by applying the Galerkin method: On the basis of Green's theorem, the finite element equation of the inner element and the outer boundary element under the closed condition can be obtained through integration by parts: Thus, the matrix form of the finite element equation is: ZZZ Assuming E is an identity matrix, the expressions of other matrices are as follows: Further simplification of the matrix equation can be obtained: where, the expression forms of K e and F e are followed as: Similarly, the finite element equation of the matrix system can be expressed as: where, the expression forms of K e and F e are followed as: The finite element equation of the cave system can be expressed as: where, the expression forms of K e and F e are as follows: Equations (22), (25), and (28), respectively, denote the finite element balance equation of the fracture system, matrix system, and cave system. Where K e is the nonlinear coefficient matrix, F e is the unit load vector, the superscript n represents the time variable. The solution process acquires giving the initial value firstly, then conducting linear treatment to obtain the pressure solution.
The wellbore storage coefficient and skin factor are taken into consideration through the Duhamel's principle.
Then, the final form of the bottom hole pressure solution is obtained after conducting the Stehfest numerical inversion: In equation (32):

Model comparison and law summary
The The value of the cross flow factor determines the time of occurrence of the concave section while the value of the elastic storativity ratio determines the depth and width of the concave. If vugs are not considered in the model, then there is only one concave on the pressure derivative curve. In general, the cross flow factor of cave system is larger, so the cross flow stage of cave system to fracture system appears earlier.
The research shows that there are still nine type flow stages in the pressure and its derivative curves even under different models, the specific flow stage division and the comparison of different models are shown in Table 1. The fractal nonlinear model contains all the characteristics of the fractal model and nonlinear model, so the model is more complicated and more factors need to be discussed (Kazemi, 1969).

Sensitivity analysis
Firstly, the sensitivity analysis of the nonlinear parameter and fractal index are conducted. The values of the known parameters are as follows: C D = 100, S = 0, x f = 0.1, x v = 0.5, k m = 0.001, k v = 0.1, L D = 5, h D = 200, a KD = 0.02. Figure 3 shows that the nonlinear parameter mainly affects the cross flow stage of the matrix system to the fracture system. The parameter decides the ability of cross flow from the matrix system to the fracture system and it can cause a later appearance of the concave of the derivative curve. In other words, the bigger the nonlinear parameter is, the shorter time the cross flow lasts, the narrower and shallower the V-shape concave of the pressure derivative curve. Even so, the nonlinear parameter does not affect the end time of the cross flow. It is obviously shown that the variation of the depth of the concave gets smaller gradually with the increasing nonlinear parameter. Figure 4 indicates that the influence of the fractal index on the pressure and derivative curves are mainly reflected in the stages after the cross flow stage from cave system to fracture system. Starting from this stage, both the pressure and derivative curves turn upward gradually with the increasing  Curves of pseudo-pressure and derivative coincide and exhibit an upward line with a slope of "1".
The models keep consistent.
2 Skin effect transition stage The pressure derivative curve exhibits like a "hump".
-Nonlinear model is consistent with Darcy model. -The pressure derivative curves of the fractal model and fractal nonlinear model are lower than the other two models. 3 Early radial flow The pressure curve continues to rise and the pressure derivative curve is shown as a level straight line, whose value is related to L D .
The results of model comparison are the same as the second flow stage. 4 Early linear flow The pressure derivative curve manifests as a straight line with a slope of "0.5", reflecting linear flow characteristics around the reaction fracture.
-Nonlinear model is consistent with Darcy model. -The slope of the pressure derivative curves of the fractal model and fractal nonlinear model are slightly larger than the other two models.

Mid-term radial flow
The pressure derivative curve behaves as a horizontal line with the value of "0.5", indicating the oil stored in the fracture system flows radially to the horizontal well.
-Nonlinear model is consistent with Darcy model. -The pressure derivative curves of the fractal model and nonlinear model still rise slowly, and the height of the derivative curves begin to exceed the other models. 6 Cross flow stage of cave system to fracture system The pressure derivative curve is characterized by a concave, which is the typical response of the oil transfers from cave system to fracture system. 7 Pseudo radial flow of cave system and fracture system The pressure derivative curve shows a "0.5" level line without considering the stress sensitivity.
-Nonlinear model is consistent with Darcy model. -The pressure and derivative curves of fractal model turn upward. The pressure curve of the fractal nonlinear model is similar to fractal model except that its derivative curve shows a more obvious upwarping. 8 Cross flow stage of matrix system to fracture system The flow stage is similar to the sixth stage, except that the cross flow coefficient of matrix system to fracture system is smaller and the concave section of derivative curve appears later.
The results of model comparison are similar to the sixth flow stage except that the concave of nonlinear model becomes narrower and shallower. 9 Late radial flow The pressure derivative curve shows as a horizontal line with the value of "0.5" without considering the stress sensitivity. It means that the production has gone up to a unified state of dynamic balance.
-Nonlinear model is consistent with Darcy model.  Figure 5 reflects that if the permeability modulus increases, the upwarping amplitude of the pressure and derivative curves of the stages after the linear flow increase because the stronger the stress sensitivity of the reservoir is, the more difficult the fluid flows, and the larger drawdown pressure is needed. Both the pressure and derivative curves turn upward gradually with the increasing of time. When the value of the permeability modulus is greater than 0.1, the upwarping amplitude increases sharply. Figure 6 shows the influence of the horizontal well position on the well test curves. When other parameters remain unchanged, if the value of z wD is close to 0.5, that is to say, the horizontal well is in the middle of the reservoir, the radial flow stage occurs earlier and the time of the early linear flow stage becomes longer. However, the location of the horizontal well does not influence other flow stages. Figure 7 shows that the flow stages after the early linear flow would be delayed because the longer the horizontal well is, the lower the pressure and derivative curves are, the longer time the early radial flow and lasts. Consequently the early linear flow regime will last longer meanwhile the appearance of the next flow periods will be postponed. It is obviously shown that the early radial flow will be covered by the wellbore storage effect if the horizontal well length is short enough. Considering the fractal characteristics and stress sensitivity of the fracture system, the pressure or derivative curves of different horizontal well lengths manifest as several parallel lines instead of converging into a single line in the late radial flow stage.

Example verification
The standard Genetic Algorithm (GA) is used to fit the measured pressure curve and output the fitting value of the parameter. First, the evolutionary algebra counter, the maximum evolution algebra and the range of parameter values are set to prepare an initial population of several initial solutions for genetic manipulation. Then, the evaluation function value (adaptation) is used to evaluate the merits of the solution or the individual and it is served as the basis for subsequent genetic operations. In addition, selection, replication and mutation operations are used to simulate the natural selection and genetic mechanism of the biological evolution (Yu et al., 2015). Finally, the optimal solution is searched by simulating the natural evolution process, the theoretical curve and the measured curve are drawn and the fitted values are output. The specific processes are as shown in Figure 8. S63 is a production well of X oilfield which was tested in September 2005. The fractal nonlinear model was selected to conduct the well test interpretation and the known parameters of the model are as shown in Table 2.  The analysis of the pressure derivative curve of the well S63 (Fig. 9) shows that the wellbore storage coefficient changes from small to large during well testing, resulting in a poor fit between the actual and theoretical pressure derivative curves in the initial flow stage. The pressure derivative curve is separated from the pressure curve and its slope is less than 1 in wellbore storage stage. In the middle and late flow stages, the pressure derivative curve of S63 well generates two "concave" which clearly reflects the characteristics of the triple-medium reservoir. The actual value fits well with the calculated value of the fractal nonlinear model accept the wellbore storage stage although the well test data is incomplete. The seepage parameters explained by the model agree with the actual situation of the well which are shown in Table 3, indicating that the fractal nonlinear model can be used for well test interpretation and seepage law analysis without considering the variation of wellbore storage coefficient.

Conclusion
1. A nonlinear seepage model of a horizontal well in a triple-medium fractal reservoir was derived by orderly adopting a series of treatments, including proposing a new nonlinear seepage model to characterize the seepage situation of matrix blocks and applying the fractal theory to describe the stress sensitivity and fractal characteristics of the fracture system. 2. To divide the flow stage and summarize the characteristics of different seepage flow models based on the bottom hole pressure curves, the processing of the internal boundary condition was conducted and the finite element principle was applied to solve the models. 3. The nine stages were summarized in the seepage process namely pure wellbore storage stage, skin effect transition stage, early radial flow, early linear flow, mid-term radial flow, cross flow stage of cave system to fracture system, pseudo radial flow of matrix system and cave system, cross flow stage of matrix system to fracture system and late radial flow. 4. Typical factors influencing the transient pressure performance were analyzed to better understand the flow characteristics and the systematic summary was made based on the sensitivity analysis of the factors.
5. The triple-medium fractal nonlinear model was applied to conduct the well test interpretation of an actual reservoir. The field application indicates that the measured value and the theoretical calculation value are better fitted if the variation of wellbore storage coefficient ignored. Since the explained seepage parameters are in line with the actual situation, it's proved that the triple-medium fractal nonlinear model can be used to conduct well test interpretation and seepage law analysis.