Analysis of Flow Behavior for Acid Fracturing Wells in Fractured-Vuggy Carbonate Reservoirs

This study develops a mathematical model for transient flow analysis of acid fracturing wells in fractured-vuggy carbonate reservoirs. This model considers a composite system with the inner region containing finite number of artificial fractures and wormholes and the outer region showing a triple-porosity medium. Both analytical and numerical solutions are derived in this work, and the comparison between two solutions verifies the model accurately. Flow behavior is analyzed thoroughly by examining the standard log-log type curves. Flow in this composite system can be divided into six or eight main flow regimes comprehensively. Three or two characteristic V-shaped segments can be observed on pressure derivative curves. Each V-shaped segment corresponds to a specific flow regime. One or two of the V-shaped segments may be absent in particular cases. Effects of interregional diffusivity ratio and interregional conductivity ratio on transient responses are strong in the early-flow period.The shape and position of type curves are also influenced by interporosity coefficients, storativity ratios, and reservoir radius significantly. Finally, we show the differences between our model and the similar model with single fracture or without acid fracturing and further investigate the pseudo-skin factor caused by acid fracturing.


Introduction
Acid fracturing is an important stimulation technique that has been widely used in developing fractured-vuggy carbonate reservoirs (FVCRs); thus, evaluating the flow behavior of acid fracturing wells in these reservoirs becomes necessary for further improving the wells' performance.
In general, fractured-vuggy carbonate rock can be considered as a triple-porosity medium, in which natural fractures, matrices, and vugs constitute the flow system [1][2][3].Recently, many scholars have made efforts to simulate the fluid transport behavior in such reservoirs through numerical and analytical methods.Wu et al. [2,3] proposed an analytical approach for transient flow analysis in the naturally fractured-vuggy reservoirs, focusing on single-phase flow and multiphase flow, respectively.Cobett et al. [4] developed a numerical model of fractured carbonate rocks and showed that numerical well testing had its limitations, especially when simulating disperse vugs.Gulbransen et al. [5] presented a multiscale mixed finite element (MsMFE) method for modeling of vuggy and naturally fractured reservoirs.Yao et al. [6] obtained the solution of the discrete fracture-vug network (DFVN) model through a standard Galerkin finite element method and provided a way of modeling realistic fluid flow in the fractured-vuggy porous media.Guo et al. [7] established a model for a horizontal well in a triple-porosity and dual-permeability system, and their results showed that type curves were dominated by external boundary as well as the permeability ratio of fracture system to the sum of fracture and matrix systems.These models provide a series of effective methods for treating the triple-porosity medium, and thus we can gain a comprehensive insight into the flow behavior in FVCRs.
However, dealing with acid fracturing wells in these reservoirs is not an easy work.Generally, carbonate minerals are usually composed of dolomite and calcite which are easy to be dissolved by hydrochloric acid [8,9].As one of the fundamental ways to stimulate carbonate reservoirs, acid fracturing alters the flow pattern in the reservoirs from a radial one to a linear one [10] and has dual characteristics of hydraulic fracturing and acidification.On the one hand, like hydraulic fracturing, multiple artificial fractures may occur near the wellbore when injecting fluids into a formation at a high pressure [11][12][13][14].Wang and Yi [15] investigated the flow behavior of a well with a finite-conductivity fracture in a carbonate reservoir.However, several field and laboratory observations have shown very complex fracture paths and presence of multiple fractures and multisegmented fractures near the wellbore after hydraulic fracturing [16][17][18][19][20]. On the other hand, the injected acid etches the face of the crack and these acid-etched patterns, namely, wormholes, remain to provide flow channels when hydraulic pressure is released [8,10,[21][22][23].Therefore, wormholes, together with artificial fractures, form the main channels of the fluid flowing into the wellbore.At present, acid fracturing has been widely implemented to enhance productivity of wells with damaged zones or produce from tight carbonate reservoirs [23][24][25].Unfortunately, complex seepage systems created by this stimulation pose challenges to study the transport processes in underground natural-resources recovery.
According to the concept of composite reservoir, when there is a discontinuity in the material properties in the radial direction from a well, the reservoir can be regarded as a composite reservoir [26][27][28].Considering that acid fracturing creates a high permeability zone near the wellbore, the case of a well with wormholes and artificial fractures in FVCRs can be represented by a composite reservoir.Therefore, the idea of composite reservoir can be used to deal with the case in this work.To the best of our knowledge, in the literature there is no similar composite model to study the flow behavior of vertical wells within wormholes and artificial fractures in FVCRs.
The purpose of this work is to conduct flow analysis of acid fracturing wells in FVCRs.We consider a composite reservoir with the inner region containing finite number of artificial fractures and wormholes, flow in both of which is linear flow, and the outer region showing a triple-porosity medium.Firstly, analytical solution of the proposed model is derived in the Laplace space, and then numerical solution of the same model is obtained with a finite difference method.The comparison results between the analytical and numerical solutions at different parameters verify the model accurately.Then, using the analytical solution, we plot the type curves and identify the main flow regimes.Effects of interregional diffusivity ratio, interregional conductivity ratio, fluid storativity ratio and interporosity coefficient of each medium, and reservoir radius on the transient response are further analyzed.Finally, we compare the differences between our model and the similar model with single fracture or without acid fracturing and also investigate the pseudo-skin factor caused by acid fracturing in detail.

Mathematical Model
FVCRs with multiple artificial fractures and wormholes are naturally structured by fracture system, vug system, and matrix system.Figure 1(a) shows a physical model of an acid fracturing well in a FVCR.The repetitive matrices and vugs are separated by a set of natural fractures.Using the concept of composite reservoir, the proposed model can be divided into two concentric regions: the inner region (Region I) and the outer region (Region II) (seen in Figure 1(b)).Figure 1(c) shows a sketch of the flow in the composite system.Similar to the dual-porosity model [29], a triple-porosity model is introduced into the outer region [3,4].Natural fractures are considered as the main flow pathways and connected with artificial fractures or wormholes directly and ultimately fluid flow from artificial fractures or wormholes into the wellbore.Wormholes are doped with artificial fractures near the wellbore, and it is difficult yet meaningless to determine their exact numbers.In practice, wormholes can be regarded as artificial fractures with larger aperture.By assuming Darcy linear flow in both kinds of flow channels, wormholes can be approximated to artificial fractures.
To simplify the mathematical problem, the basic assumptions are made: (1) In the inner region, there are a finite number, , of artificial fractures and wormholes with the same aperture, , permeability,  1 , and storage capacity, () 1 , which are conceptualized as vertical plates extending from the top to the bottom.This region's radius is   , equal to the length of artificial fractures or wormholes.
(2) In the outer region, the permeability and storage capacity for natural fractures are  2 and () 2 , for vugs are  2V and () 2V , and for matrices are  2 and () 2 .Natural fractures, vugs, and matrices are uniformly distributed in this region.
(3) All rock properties, such as permeability and porosity, are constant.The flow in the reservoir is considered to be a single-phase flow with slightly compressible fluid of constant viscosity  and constant formation volume factor .
(4) The reservoir is of uniform thickness with impermeable upper and lower boundaries and closed with a radius, equal to   .
(5) A vertical well is located in the center of the inner region and intercepted by these artificial fractures and wormholes.The well is produced at a constant rate  and no wellbore storage is considered.
(6) Isothermal and Darcy flow are assumed in two regions.

Establishment of Mathematical Model
2.1.1.Flow in Region I.In the inner region, flow in the artificial fractures and wormholes follows a linear flow pattern; thus, the governing equation for this region can be written as Initial condition is Well production condition for constant rate is

Flow in Region II.
Based on the triple-porosity model, any infinitesimal volume in the outer region contains a large proportion of three constitutive media.Each point in the model is assigned with three pressures.Each system interacts independently with the other two systems.The matrix and vug systems provide the main storage space but have no direct contributions to flow and transport.
The governing equations for natural fracture, vug, and matrix system are respectively.The subscripts 2, 2V, and 2 refer to natural fracture, vug, and matrix system, respectively; and  2 ,  2V , and  2V are the volumetric flux from matrix to natural Mathematical Problems in Engineering fracture, from vug to matrix, and from vug to natural fracture, respectively.Assuming that fluid flow between two systems follows the pseudo-steady state flow [29],  2 ,  2V , and  2V can be given as where   ,  V , and  V are the interporosity flow shape factors, which depend on the geometry of the interporosity flow and have dimensions of reciprocal area.By substituting ( 7)-( 9) into ( 4)-( 6), the previous governing equations can be simplified as respectively.
The initial and closed external boundary conditions in the outer region are Dimensionless pressure of natural fracture system is Dimensionless pressure of vug system is Dimensionless pressure of matrix system is Dimensionless time is Dimensionless radius is Dimensionless radius of external boundary is Dimensionless radius of wellbore is Interregional diffusivity ratio is Interregional conductivity ratio is Fluid storativity ratio of fracture system is Fluid storativity ratio of vug system is Fluid storativity ratio of matrix system is Interporosity coefficient of matrix system to fracture system is Interporosity coefficient of vug system to fracture system is Interporosity coefficient of vug system to matrix system is It is worth noting that there is a certain relationship among three storativity ratios:  2 + 2 + 2V = 1.In a FVCR, the interporosity coefficient of vug system to fracture system  V should be larger than that of matrix system to fracture system   .The differences in these interporosity coefficients lead to different response time of each system on type curves [1].
(1) Region I.The governing differential equation, initial condition, and inner boundary condition become respectively.
(2) Region II.The dimensionless governing differential equations of natural fracture, vug, and matrix system are respectively.
The initial and external boundary conditions become (3) Interface Conditions.The dimensionless equations of pressure and flux at the interface are given by

Laplace Transform of Dimensionless Mathematical Model.
The Laplace transform with respect to   is based on the following function: where  is the Laplace transform variable.
For the inner region, applying the Laplace transform to (31)-( 33), we have For the outer region, by employing the Laplace transform, ( 34)-( 37) can be simplified as where The external boundary condition in Laplace space is Similarly, taking Laplace transform for the dimensionless interface boundary conditions, we have =1 . (47)

Analytical Solution in Laplace
Space.The general solutions for ( 42) and ( 44) are given by where  0 and  0 are modified Bessel functions of zero order of the first and second kind, respectively.
Combining the general solutions with all the boundary conditions, we can get the analytical solution and then the wellbore pressure in Laplace domain can be obtained: where In real space, the dimensionless bottomhole pressure (BHP) can be obtained by using Stehfest numerical inversion [30] to convert p () to   (  ).

Finite Difference Form of Dimensionless Mathematical
Model.Though the analytical solution has been obtained in the previous section, we still seek to get the numerical solution of this model with a finite difference method [31], which can help us to verify the accuracy of the analytical solution.
As shown in Figure 2, we divide the reservoir into  segments in the radial direction, in which the inner region contains  segments of equal length and the outer region contains  −  segments of equal length.To ensure the computational efficiency and accuracy, the segment length of the inner region is shorter than that of the outer region.Meanwhile, in real space, time is discretized to obtain transient bottomhole pressure (BHP).

Finite Difference Model of Region I.
According to the uniform block-centered grid system shown in Figure 2, we can differentiate the governing equation (31) at segment  at a timestep  as which can be rearranged as where  1 = (Δ 2   /Δ  ) + 2 and    represents the dimensionless pressure at the th segment at a timestep  in the inner region.Equation (52) applies at all segments in the inner region,  = 2, 3, . . .,  − 1, except the first and last segments of the inner region.
Douglas and Peaceman [32] suggest treating no-flow boundaries by using a mirror image method.In the discretized system, taking the segment located at the wellbore ( = 1) as a mirror, a pseudo segment ( = 0) corresponding to the second segment ( = 2) can be obtained on the other side of the wellbore (shown in Figure 2).  0 ,   1 , and   2 represent the dimensionless pressure at the pseudo segment ( = 0), the first segment ( = 1), and the second segment ( = 2) at a timestep , respectively.Then, the inner boundary condition (see (32)) can be given as Rearranging (53) gives Applying (52) at the inner boundary of segment  = 1, we obtain Then, adding (54) and (55) together gives where  1 = 2Δ  .Equation ( 56) is the differentiated form of flow equation for the first segment ( = 1), which also corresponds to the inner boundary condition.

Finite Difference Model of Region II.
Similarly, the dimensionless governing equations of natural fracture system, vug system, and matrix system in the outer region can be differentiated as respectively.   , V   , and    represent the dimensionless pressure of fracture, vug, and matrix at the th segment at a timestep , respectively.
Similarly, using the mirror image method at the external boundary, a pseudo segment ( =  + 1) corresponding to the segment ( =  − 1) can be obtained at the outside of the outer region (shown in Figure 2).  +1 at a timestep , respectively.Then, the external boundary condition (see (38)) can be given as which can be rearranged as Equation ( 63) applies at the external boundary of segment  = , and the following equation is set up: Combining ( 66) and (67), the following equation can be obtained: Equation ( 68) is the differentiated form of flow equation for the last segment ( = ) in the outer region.

Finite Difference Treatment of Interface Condition.
Applying (52) for the segment ( = ) and (63) for the segment ( =  + 1) gives The differentiated form of the continuity of pressure at the interface (see (39)) can be easily got: Combining ( 71) and (69) will yield Equation ( 72) is the differentiated form of flow equation for the segment ( = ) at the inside of the interface, which just corresponds to the last segment of the inner region.
The equation of the continuity of flux at the interface (see (40)) can be given as which can be rearranged as Combining ( 70), (72), and (74) gives where Equation ( 75) is the differentiated form of flow equation for the segment ( =  + 1) at the outside of the interface, corresponding to the first segment of the outer region.

Numerical Solution of Dimensionless Mathematical
Model.Using the differentiated form for all the segments ( = 1, 2, 3, . . ., ,  + 1, . . ., ), a tridiagonal matrix with the rank of  ×  is formed, which is a linear system: ) . ( The initial conditions in the differentiated form for the inner and outer region are given as 0  = V 0  =  0  = 0  =  + 1,  + 2, . . ., . 1 represents the wellbore pressure (BHP) at a timestep .An iterative method is used to get the wellbore pressure at any timestep; thus, transient responses corresponding to different time can be obtained.In this paper, Thomas chasing method, one of the classical iterative methods, is used to solve the tridiagonal matrix.It is found that the solution does not change appreciably when the reservoir with a dimensionless radius of 100 is divided into 4000 segments.In all cases studied, 10 time intervals are taken in each log cycle of dimensionless time.

Model Verification
The model in this study can be verified by comparing the results of its analytical and numerical solution.Figure 3 shows the comparison results of the analytical and numerical solutions at different parameters.All the analytical solutions agree well with the corresponding numerical solutions within the most production life.At early time and some of intermediate time there is a slight discrepancy between the two kinds of solutions on pressure derivative curves, which may be caused by the long length of discrete segments, but the error is within an acceptable limit.These validations indicate that the analytical model in this work is reliable.In addition to acting as an effective way to verify the accuracy of the analytical solution, the numerical solution in this study has some advantages over the analytical solution in computing speed.However, the analytical solution also has its unique advantages in some respects: for example, wellbore storage and skin effect can be easily added to the analytical solution in Laplace domain [33][34][35], which may be extremely difficult for the numerical solution.

Results and Discussion
One objective of this research is to calculate type curves for various reservoirs and observe the flow behavior from these curves.By type curves matching, some reservoir parameters, such as permeability, skin factor, and reservoir drainage area, can be obtained [36,37].In this section, we discuss the characteristics of transient response of an acid fracturing well in a FVCR and also analyze the differences between our model and the similar model with single fracture or without acid fracturing and further investigate the pseudo-skin factor caused by acid fracturing.

Flow Regime Recognition.
In the carbonate reservoirs, either matrix system or vug system provides the major storage space; thus, these reservoirs can be subdivided into two categories: reserves dominated by matrix system and reserves dominated by vug system.Figures 4 and 5 show the complete transient flow processes in these two categories.The values of the parameters used in Figures 4 and 5 are listed in Tables 1  and 2, respectively.4, the flow process in this reservoir can be divided into eight flow regimes.Initially, there is a linear flow period characterized by a half-slope straight line on the pressure derivative curve.After a transition flow regime, the first interporosity flow develops and is followed with the second transition flow.Then, the second interporosity flow occurs and immediately is followed by the third transition flow.Eventually, a short radial flow regime may be observed before the system reaches a pseudo-steady flow.Interestingly, three characteristic V-shaped segments can be observed on the pressure derivative curve.
(1) Linear Flow Regime.In this regime, the production is dominated by the fluid stored in artificial fractures and wormholes.It occurs at very small values of dimensionless time.The pressure and pressure derivative curves are an upward straight line with a slope of 0.5.This regime is typical of fluid flow for acid fracturing wells.Before the interface is felt, the wellbore pressure solution in this study can be approximate to Corresponding to the square root of the dimensionless time, this equation implies that a log-log plot of function results a straight line with the slope of one-half, which is the typical characteristic of linear flow.This equation also indicates that the wellbore pressure in this regime depends on the values of Both the pressure and derivative curves show a unit slope straight line, which is the typical characteristic of closed reservoirs.By means of asymptotic analysis, the approximate solution in this regime can be yielded: This solution shows that wellbore pressure in pseudosteady regime only depends on natural fracture storativity ratio and reservoir radius (seen in Figures 6-14).

Case 2: Hydrocarbon Reserves Dominated by Vug System.
As is shown in Figure 5, flow in this reservoir can be divided into six flow regimes: linear flow, transition flow (a), interporosity flow, transition flow (b), radial flow, and pseudo-steady flow.The major difference of flow characteristics between two kinds of carbonate reservoirs lies in the interporosity flow period.The characteristics of the other flow regimes are the same with the corresponding regimes of the previous case, including the approximate solutions of linear flow and pseudo-steady flow.
Because the ability of vug system supplying fluid to natural fracture system is much stronger than that of matrix system supplying fluid to natural fracture system, the last two V-shaped segments are connected to each other and form one V-shaped segment.Consequently, only two V-shaped segments can be observed on the pressure derivative curve and the second V-shaped segment merges with the radial flow.

Parameter Sensitivity Analysis.
The shape rendered from real data may be distorted by noises in well testing, which makes it necessary to establish the stylized shapes under different parameter conditions [36].Figures 6-14 contain the type curves of pressure-transient response for an acid fracturing well in FVCRs.Later we will show how the shape of the type curves may be affected by changing the parameters.

Effect of Interregional Diffusivity Ratio.
According to the dimensionless definition, interregional diffusivity ratio  represents the ratio of fluid transmissivity ability between natural fractures in the outer region and artificial fractures or wormholes in the inner region.A larger  means a weaker diffusion ability of wormholes or artificial fractures.Figure 6 shows the effect of this parameter on pressure-transient response ( = 0.001, 0.005, 0.05, and 0.5).We can find that interregional diffusivity ratio mainly affects early-flow period.Before the first interporosity flow regime, a larger  leads to lower pressure depletion at the same production rate.Meanwhile, as this parameter increases, due to a weaker diffusion ability in the inner region, duration of the linear flow is prolonged, while duration of the first interporosity flow is shortened.When interregional diffusivity ratio gradually increases to 1, the first V-shaped segment becomes shallower and even disappears.All the curves for different interregional diffusivity ratios begin to normalize at the end of the first interporosity flow regime.

Effect of Interregional Conductivity Ratio.
According to the dimensionless definition, interregional conductivity ratio  represents a measure of the conductivity between wormholes or artificial fractures and natural fractures.In practice, a larger  means a lower conductivity of the inner region.Transient responses on the log-log plots for varied interregional conductivity ratio are shown in Figure 7 ( = 0.2, 0.4, 0.6, and 0.8).It can be concluded that this parameter mainly affects early-and middle-flow periods.A larger interregional conductivity ratio leads to more pressure depletion at the same production.Meanwhile, with the increase of this parameter, the hump corresponding to the first transition flow becomes more significant.Besides, the pressure curves for varied interregional conductivity ratio normalize in the late-flow period, while the pressure derivative curves begin to normalize from the first interporosity flow regime.

Effect of Storativity Ratios.
There are three storativity ratios in a FVCR.Each ratio represents the relative capacity of fluid stored in a corresponding system.Generally fluid stored in the natural fracture system is much less than that stored in the vug or matrix system.
Figure 8 presents the results of sensitivity analysis of fracture storativity ratio ( 2 = 0.001, 0.005, 0.025, and 0.125).It indicates that a larger  2 leads to more pressure depletion in the middle-and late-flow periods.Fracture storativity ratio begins to affect the transient response at the start of the first interporosity flow, which extends to the pseudostate flow period; thus, the curves do not normalize, but they are parallel to each other in the late-flow period.Due to a higher conductivity and diffusion capacity for natural fracture system provided by a larger fracture storativity ratio, when this ratio increases from 0.001 to 0.125, the last two Vshaped segments become shallower and even the second Vshaped segment may disappear.Moreover, the pseudo-steady flow begins earlier with the increase of this parameter.
When investigating the effect of another two storativity ratios, we keep  2 constant; thus, matrix storativity ratio has a reciprocal relationship with vug storativity ratio.Therefore, it is not necessary to discuss all the values in their own range.
Figure 9 shows the response of varying matrix storativity ratio,  2 = 0.6, 0.7, 0.8, and 0.9, which indicates that matrix system provides the majority of storage space.As is described in Figure 9, matrix storativity ratio mainly affects the middleflow period, and the pressure and derivative curves for different matrix storativity ratios overlap completely in the earlyand late-flow periods.But the shape of the last two V-shaped segments is affected by this storativity ratio significantly.With the increase of matrix storativity ratio, the second Vshaped segment becomes shallower, while the third V-shaped segment becomes deeper.In fact, matrix storativity ratio only affects the third V-shaped segment and vug storativity ratio only affects the second V-shaped segment.However, due to the reciprocal relationship between two storativity ratios, matrix storativity ratio affects the second V-shaped segment indirectly.
Figure 10 presents the results of sensitivity analysis of vug storativity ratio,  2V = 0.6, 0.7, 0.8, and 0.9, which indicates that the reserves are dominated by vug system.As shown in Figure 10, vug storativity ratio only has a slight effect on the middle-flow period and the related curves normalize in the late-flow period.With the increase of vug storativity ratio, the concave of the second V-shaped segment deepens slightly.

Effect of Interporosity
Coefficients.There are also three interporosity coefficients in the triple-porosity reservoir.A sensitivity analysis is also performed by varying these three coefficients.
Figure 11 demonstrates the effect of matrix-fracture interporosity coefficient on the transient behavior (  = 5 × 10 −5 , 5 × 10 −4 , 5 × 10 −3 , and 10 −2 ).It can be seen that this parameter mainly affects the middle-flow period and type curves for different interporosity coefficients normalize in the late-flow period.As this coefficient increases, the second Vshaped segment deepens, the third V-shaped segment moves to the left side, and the duration of radial flow regime is prolonged.This is mainly because a larger matrix-fracture interporosity coefficient means a stronger ability of matrix system supplying fluid for natural fracture system, which makes the arrival time of matrix-fracture interporosity flow regime earlier.
Figure 12 presents the results of sensitivity analysis of matrix-vug interporosity coefficient ( V = 0.001, 0.005, 0.025, and 0.05).This coefficient also mainly affects middleflow period, especially the matrix-fracture interporosity flow regime, corresponding to the third V-shaped segment.It can be interpreted by that this coefficient only affects the ability of vug system to supply fluid for matrix system.With the increase of this coefficient, the starting time of the second interporosity flow moves back and the corresponding Vshaped segment deepens greatly.Similarly, the pressure and derivative curves will normalize in the late-flow period.
Figure 13 shows the response of varying vug-fracture interporosity coefficient ( V = 0.01, 0.05, 0.25, and 0.5).This parameter mainly affects the middle-flow period, and the normalization can also be observed in the late-flow period.As this coefficient increases, the starting time of the first interporosity flow moves back and the third V-shaped segment becomes shallower.It can be interpreted by that the ability of vug system to supply fluid for fracture system becomes much stronger than that of matrix system supplying fluid for fracture system with the increase of this coefficient.

Effect of Dimensionless Reservoir
Radius. Figure 14 presents the results of sensitivity analysis of dimensionless reservoir radius (  = 30, 50, 70, and 90).There is no difference in the shape of type curves in the early-and middleflow periods under different reservoir radius; however, the pressure and pressure derivative curves show a unit slope straight line and go up rapidly in the late-flow period.In addition, with the increase of the dimensionless reservoir radius, the duration of radial flow regime is prolonged significantly.

6.3.
Comparison with Other Similar Models.Gringarten et al. [38] investigated unsteady-state pressure distributions created by a well with a single infinite-conductivity vertical fracture in a homogeneous reservoir.Further, we can easily obtain the transient-pressure solution for the same well type in a FVCR.There are a finite number of artificial fractures and wormholes in our model's inner region; thus, it is necessary to compare these two models' difference in flow behavior.Meanwhile, it is also significant to make a comparison of our model with the model not having acid fracturing to know the effect of acidification.In order to have a convenient comparison, we firstly make a comparison to the parameters used in three models and then compare their type curves by fixing a set of the same parameter values.Table 3 shows the differences between three models' parameters.It can be seen that the same parameters are  2 ,  2 ,  2V ,   ,  V ,  V ,   , and   , but the parameters  and  only belong to the model with acid fracturing.To make an effective comparison, we set the same values for the same parameters, which are shown in Table 3. Figure 15 shows the comparison result of type curves between three models, and we can observe the differences clearly: (i) in the early-flow period, the pressure and derivative curves for the model with acid fracturing and the model with single fracture coincide perfectly, respectively.It demonstrates that these two models' flow behavior is the same in this period.However, because drainage area of the well without acid fracturing is much smaller than that of the well with acid fracturing or with single fracture, the pressure depletion of the former case is faster than that of the latter two cases at the same group of production parameters; (ii) in the middle-flow period, the pressure depletion of the well without single fracture becomes faster than that of the well with acid fracturing, which also can be interpreted by the difference of drainage area; (iii) in the models with single fracture and without acid fracturing, only two V-shaped segments can be observed on the derivative curves, and the first V-shaped segment in the model with acid fracturing is replaced by a unit slope straight line or a horizontal line with the value of 0.5, which are the typical characteristics of linear flow and radial flow, respectively.

Parameter Influence on Pseudo-Skin Factor.
During the pseudo-steady flow period, for the same position, there is a difference of the dimensionless pressure between the models with acid fracturing and without acid fracturing.This difference remains constant at large value of dimensionless time and can be handled as a pseudo-skin factor [39,40] where Note that  1 and  2 have been given in (49).Because of higher connectivity between the wellbore and the reservoir caused by acid fracturing, the pseudoskin factor defined here is always negative and indicates an increase in well productivity.The larger the absolute value of the pseudo-skin factor, the better the stimulated effect of acid fracturing.
As is shown in Figure 16, under different interregional diffusivity ratios, the pseudo-skin factors at the same position near the wellbore are the same, which indicates that this parameter has no influence on the ability of enhancing the productivity in the late-flow period.In the pseudo-steady flow regime, the diffusion capacity of the inner region well satisfies the demand of fluid flow feeding the wellbore and the flow in the model reaches a dynamic balance; thus, the change of interregional diffusivity ratio does not affect the productivity.It also can be concluded that the closer the position to the wellbore, the better the stimulation effect.Figure 17 presents different pseudo-skin factors for the vertical well of varied interregional conductivity ratio.With the increase of interregional conductivity ratio, the absolute value of the pseudo-skin factor at the same position decreases, which implies that the ability of enhancing the well production decreases.As we know, a larger interregional conductivity ratio corresponds to a weaker conductivity of the inner region and also means the increase of the flow resistance, which explains the result above.Similarly, in this group, the closer the position to the wellbore, the better the stimulation effect.

Conclusions
This work presented a composite reservoir model to analyze the flow behavior of an acid fracturing well in a fracturedvuggy carbonate reservoir.Both the analytical and numerical solutions are obtained in this study.Based on this work, several important conclusions are obtained: (1) All the analytical solutions show well agreement with the corresponding numerical solutions within the most production life for different parameters.These verifications indicate that the analytical model in this work is reliable.The numerical solution has some advantages over the analytical solution in computing speed, but the analytical solution is convenient to handle the condition with wellbore storage or skin effect.
(2) Due to that reserves in the fractured-vuggy carbonate reservoir are dominated by either matrix system or vug system; for an acid fracturing well, flow in this composite model can be divided into six or eight flow regimes comprehensively.Linear flow regime is typical of fluid flow for acid fracturing wells in such reservoirs.For the two scenarios above, three and two characteristic V-shaped segments can be observed on their own pressure derivative curves, respectively.

Mathematical Problems in Engineering
(3) Effects of interregional diffusivity ratio and interregional conductivity ratio on transient responses are strong in the early-flow period.The existence and duration of linear flow and the first transition flow are mainly determined by these two parameters.As interregional diffusivity ratio increases, duration of linear flow is prolonged while that of the first transition flow is shortened.Besides, a lower interregional diffusivity ratio or a larger interregional conductivity ratio leads to more pressure depletion in the early-flow period at the same production.
(4) The shape and position of type curves are also dominated by storativity ratios and interporosity coefficients.In the late-flow period, for different interregional diffusivity ratio, interregional conductivity ratio, and interporosity coefficient, the pressure and pressure derivative curves will normalize.However, type curves for varied fracture storativity ratio and dimensionless reservoir radius are parallel to each other.
(5) Compared with the models with single fracture or without acid fracturing, the pressure depletion of the well with acid fracturing is slower than that of another two models at the same group of production parameters in the middleand late-flow period.Due to a higher connectivity caused by acid fracturing, the pseudo-skin factor defined in this work is always negative and indicates an increase in well productivity.Interregional diffusivity ratio has no influence on the ability of enhancing the productivity in the late-flow period.However, the ability of enhancing the well production decreases with the increase of interregional conductivity ratio.

Figure 1 :
Figure 1: Schematic diagram: (a) an acid fracturing well in a fractured-vuggy carbonate reservoir; (b) a composite reservoir model with two concentric regions; and (c) a sectional view of fluid transport in the composite system.

3 & 1 Figure 3 :
Figure 3: Comparison results of analytical and numerical solutions at different parameters.

Figure 15 :
Figure 15: The comparison of type curves between three different models.

1 r 4 Figure 16 :
Figure 16: Pseudo-skin factor versus dimensionless formation radius for an acid fracturing well at different interregional diffusivity ratios.

1 Figure 17 :
Figure 17: Pseudo-skin factor versus dimensionless reservoir radius for an acid fracturing well at different interregional conductivity ratios.

Table 1 :
The data used in Figure4.

Table 2 :
The data used in Figure5.

Table 3 :
The comparison of different parameters between different models.