Research Article  Open Access
Mingxian Wang, Zifei Fan, Xuyang Dong, Heng Song, Wenqi Zhao, Guoqing Xu, "Analysis of Flow Behavior for Acid Fracturing Wells in FracturedVuggy Carbonate Reservoirs", Mathematical Problems in Engineering, vol. 2018, Article ID 6431910, 20 pages, 2018. https://doi.org/10.1155/2018/6431910
Analysis of Flow Behavior for Acid Fracturing Wells in FracturedVuggy Carbonate Reservoirs
Abstract
This study develops a mathematical model for transient flow analysis of acid fracturing wells in fracturedvuggy 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 tripleporosity 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 loglog type curves. Flow in this composite system can be divided into six or eight main flow regimes comprehensively. Three or two characteristic Vshaped segments can be observed on pressure derivative curves. Each Vshaped segment corresponds to a specific flow regime. One or two of the Vshaped segments may be absent in particular cases. Effects of interregional diffusivity ratio and interregional conductivity ratio on transient responses are strong in the earlyflow 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 pseudoskin factor caused by acid fracturing.
1. Introduction
Acid fracturing is an important stimulation technique that has been widely used in developing fracturedvuggy 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, fracturedvuggy carbonate rock can be considered as a tripleporosity medium, in which natural fractures, matrices, and vugs constitute the flow system [1–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 fracturedvuggy reservoirs, focusing on singlephase 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 fracturevug network (DFVN) model through a standard Galerkin finite element method and provided a way of modeling realistic fluid flow in the fracturedvuggy porous media. Guo et al. [7] established a model for a horizontal well in a tripleporosity and dualpermeability 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 tripleporosity 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–14]. Wang and Yi [15] investigated the flow behavior of a well with a finiteconductivity 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–20]. On the other hand, the injected acid etches the face of the crack and these acidetched patterns, namely, wormholes, remain to provide flow channels when hydraulic pressure is released [8, 10, 21–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–25]. Unfortunately, complex seepage systems created by this stimulation pose challenges to study the transport processes in underground naturalresources 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–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 tripleporosity 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 pseudoskin factor caused by acid fracturing in detail.
2. 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 dualporosity model [29], a tripleporosity 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.
(a)
(b)
(c)
To simplify the mathematical problem, the basic assumptions are made:
In the inner region, there are a finite number, , of artificial fractures and wormholes with the same aperture, , permeability, , and storage capacity, , 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.
In the outer region, the permeability and storage capacity for natural fractures are and , for vugs are and , and for matrices are and . Natural fractures, vugs, and matrices are uniformly distributed in this region.
All rock properties, such as permeability and porosity, are constant. The flow in the reservoir is considered to be a singlephase flow with slightly compressible fluid of constant viscosity and constant formation volume factor .
The reservoir is of uniform thickness with impermeable upper and lower boundaries and closed with a radius, equal to .
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.
Isothermal and Darcy flow are assumed in two regions.
2.1. 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
2.1.2. Flow in Region II
Based on the tripleporosity 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 arerespectively. The subscripts , , and refer to natural fracture, vug, and matrix system, respectively; and , , and are the volumetric flux from matrix to natural fracture, from vug to matrix, and from vug to natural fracture, respectively. Assuming that fluid flow between two systems follows the pseudosteady state flow [29], , , and can be given aswhere , , and 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 asrespectively.
The initial and closed external boundary conditions in the outer region are
2.1.3. Interface Conditions
Natural fractures in the outer region are connected with artificial fractures and wormholes in the inner region directly. Because of the continuity of pressure and flux at the interface, the following conditions must hold along the surface:
2.2. Dimensionless Mathematical Model
2.2.1. Dimensionless Definitions
Dimensionless pressure of artificial fracture and wormhole system in the inner region is
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: . In a FVCR, the interporosity coefficient of vug system to fracture system 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].
2.2.2. Dimensionless Mathematical Model in Real Space
Introducing the dimensionless variables into (1)–(3) and (10)–(14), we obtain the following dimensionless equations.
(1) Region I. The governing differential equation, initial condition, and inner boundary condition becomerespectively.
(2) Region II. The dimensionless governing differential equations of natural fracture, vug, and matrix system arerespectively.
The initial and external boundary conditions become
(3) Interface Conditions. The dimensionless equations of pressure and flux at the interface are given by
3. Analytical Solution
3.1. 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 aswhere
The external boundary condition in Laplace space is
Similarly, taking Laplace transform for the dimensionless interface boundary conditions, we have
3.2. Analytical Solution in Laplace Space
The general solutions for (42) and (44) are given bywhere and 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 to .
4. Numerical Solution
4.1. 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).
4.1.1. Finite Difference Model of Region I
According to the uniform blockcentered grid system shown in Figure 2, we can differentiate the governing equation (31) at segment at a timestep aswhich can be rearranged aswhere 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, , except the first and last segments of the inner region.
Douglas and Peaceman [32] suggest treating noflow boundaries by using a mirror image method. In the discretized system, taking the segment located at the wellbore () as a mirror, a pseudo segment () corresponding to the second segment () can be obtained on the other side of the wellbore (shown in Figure 2). , , and represent the dimensionless pressure at the pseudo segment (), the first segment (), and the second segment () 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 , we obtain
Then, adding (54) and (55) together giveswhere . Equation (56) is the differentiated form of flow equation for the first segment (), which also corresponds to the inner boundary condition.
4.1.2. 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. , , and represent the dimensionless pressure of fracture, vug, and matrix at the th segment at a timestep , respectively.
Combining (58) and (59), the following equations can be derived:where
Then, substituting (60) and (61) into (57) giveswhere
Equation (63) applies at all segments in the outer region, , except the first and last segments of the outer region.
Similarly, using the mirror image method at the external boundary, a pseudo segment () corresponding to the segment () can be obtained at the outside of the outer region (shown in Figure 2). , , and represent the dimensionless fracture pressure at segment , , and at a timestep , respectively. Then, the external boundary condition (see (38)) can be given aswhich 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.
4.1.3. Finite Difference Treatment of Interface Condition
Applying (52) for the segment () and (63) for the segment () 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 aswhich can be rearranged as
Combining (70), (72), and (74) giveswhere
Equation (75) is the differentiated form of flow equation for the segment () at the outside of the interface, corresponding to the first segment of the outer region.
4.2. Numerical Solution of Dimensionless Mathematical Model
Using the differentiated form for all the segments (), 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
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.
5. 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–35], which may be extremely difficult for the numerical solution.
(a) , 0.5
(b) , 0.5
(c) , 0.5
(d) , 0.9
(e) , 0.7
(f) , 10−2
(g) , 0.05
(h) , 0.1
6. 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 pseudoskin factor caused by acid fracturing.
6.1. 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.


6.1.1. Case 1: Hydrocarbon Reserves Dominated by Matrix System
As is described in Figure 4, the flow process in this reservoir can be divided into eight flow regimes. Initially, there is a linear flow period characterized by a halfslope 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 pseudosteady flow. Interestingly, three characteristic Vshaped 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 toCorresponding to the square root of the dimensionless time, this equation implies that a loglog plot of function results a straight line with the slope of onehalf, which is the typical characteristic of linear flow. This equation also indicates that the wellbore pressure in this regime depends on the values of interregional diffusivity ratio and interregional conductivity ratio (seen in Figures 6 and 7).
(2) Transition Flow Regime (a). The first transition flow takes place when the interface takes effect. The derivative curve in this regime shows the first Vshaped segment, which represents fluid flow from natural fractures into artificial fractures or wormholes. It can be interpreted by the fact that as the pressure of the inner region depletes gradually, natural fracture system responds first owing to the direct connection to the inner region. The shape of this segment is also controlled by interregional diffusivity ratio and interregional conductivity ratio, and this Vshaped segment may not be observed at particular parameters (seen in Figures 6 and 7).
(3) Interporosity Flow Regime (a). The first interporosity flow regime corresponds to the second Vshaped segment on the derivative curve, which represents the crossflow between vug system and natural fracture system. This interporosity flow takes place first because vug permeability is larger than matrix permeability. This segment is predominantly determined by storativity ratios of natural fracture system and matrix system, and interporosity coefficient between vugs and natural fractures, and also may be absent from the curves at particular parameters (seen in Figures 8–10, and 13).
(4) Transition Flow Regime (b). With the depletion of vug system, the matrix system begins to supply fluid to natural fracture system gradually and a short transition flow occurs. The duration of this regime depends on the difference of the ability supplying fluid between matrix system and vug system; thus, this regime may be so short that it can hardly be identified from the derivative curve.
(5) Interporosity Flow Regime (b). In this interporosity flow regime, the production is mainly dominated by the fluid flow from matrix system to fracture system. This regime is characterized by the third Vshaped segment on the derivative curve, which is governed by two physical processes: one is interporosity flow from matrix system to fracture system and the other is interporosity flow from vug system to matrix system. The former has greater influence on the shape of this regime than the latter (seen in Figures 11 and 12).
(6) Transition Flow Regime (c). With the depletion of matrix system, fluids far away from the wellbore begin to flow. Like the second transition flow, another short transition flow regime may develop before the dynamic balance reaches the whole reservoir. It also may not be observed on the pressure derivative curve.
(7) Radial Flow Regime. When the matrix system depletes fast and the reservoir radius is large enough, fluid flow in the whole reservoir is likely to reach a dynamic balance; thus, a radial flow occurs after the third transition flow. In this regime, the slope of pressure derivative curve is zero, and the pressure derivative converges to “0.5 lines.” As shown in Figure 14, reservoir radius has a great influence on the duration of radial flow.
(8) PseudoSteady Flow Regime. For closed boundary, the longtime behavior corresponds to a pseudosteady flow. 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).
6.1.2. 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 pseudosteady 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 pseudosteady 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 Vshaped segments are connected to each other and form one Vshaped segment. Consequently, only two Vshaped segments can be observed on the pressure derivative curve and the second Vshaped segment merges with the radial flow.
6.2. 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 pressuretransient 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.
6.2.1. 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 pressuretransient response (, 0.005, 0.05, and 0.5). We can find that interregional diffusivity ratio mainly affects earlyflow 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 Vshaped 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.
6.2.2. 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 loglog plots for varied interregional conductivity ratio are shown in Figure 7 (, 0.4, 0.6, and 0.8). It can be concluded that this parameter mainly affects early and middleflow 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 lateflow period, while the pressure derivative curves begin to normalize from the first interporosity flow regime.
6.2.3. 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 ( = 0.001, 0.005, 0.025, and 0.125). It indicates that a larger leads to more pressure depletion in the middle and lateflow 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 lateflow 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 pseudosteady flow begins earlier with the increase of this parameter.
When investigating the effect of another two storativity ratios, we keep 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, = 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 early and lateflow periods. But the shape of the last two Vshaped segments is affected by this storativity ratio significantly. With the increase of matrix storativity ratio, the second Vshaped segment becomes shallower, while the third Vshaped segment becomes deeper. In fact, matrix storativity ratio only affects the third Vshaped segment and vug storativity ratio only affects the second Vshaped segment. However, due to the reciprocal relationship between two storativity ratios, matrix storativity ratio affects the second Vshaped segment indirectly.
Figure 10 presents the results of sensitivity analysis of vug storativity ratio, = 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 middleflow period and the related curves normalize in the lateflow period. With the increase of vug storativity ratio, the concave of the second Vshaped segment deepens slightly.
6.2.4. Effect of Interporosity Coefficients
There are also three interporosity coefficients in the tripleporosity reservoir. A sensitivity analysis is also performed by varying these three coefficients.
Figure 11 demonstrates the effect of matrixfracture 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 middleflow period and type curves for different interporosity coefficients normalize in the lateflow period. As this coefficient increases, the second Vshaped segment deepens, the third Vshaped segment moves to the left side, and the duration of radial flow regime is prolonged. This is mainly because a larger matrixfracture interporosity coefficient means a stronger ability of matrix system supplying fluid for natural fracture system, which makes the arrival time of matrixfracture interporosity flow regime earlier.
Figure 12 presents the results of sensitivity analysis of matrixvug interporosity coefficient ( = 0.001, 0.005, 0.025, and 0.05). This coefficient also mainly affects middleflow period, especially the matrixfracture interporosity flow regime, corresponding to the third Vshaped 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 lateflow period.
Figure 13 shows the response of varying vugfracture interporosity coefficient ( = 0.01, 0.05, 0.25, and 0.5). This parameter mainly affects the middleflow period, and the normalization can also be observed in the lateflow period. As this coefficient increases, the starting time of the first interporosity flow moves back and the third Vshaped 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.
6.2.5. 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 lateflow 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 unsteadystate pressure distributions created by a well with a single infiniteconductivity vertical fracture in a homogeneous reservoir. Further, we can easily obtain the transientpressure 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 , , , , , , , 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 earlyflow 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 middleflow 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 Vshaped segments can be observed on the derivative curves, and the first Vshaped 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.
6.4. Parameter Influence on PseudoSkin Factor
During the pseudosteady 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 pseudoskin factor [39, 40] that can be defined as where and represent the dimensionless pressure near the wellbore for the model with acid fracturing and without the stimulation, respectively.
Applying the Laplace transform to (82), we havewhereNote that and 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 pseudoskin factor, the better the stimulated effect of acid fracturing.
As is shown in Figure 16, under different interregional diffusivity ratios, the pseudoskin 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 lateflow period. In the pseudosteady 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 pseudoskin factors for the vertical well of varied interregional conductivity ratio. With the increase of interregional conductivity ratio, the absolute value of the pseudoskin 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.
7. 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:
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.
Due to that reserves in the fracturedvuggy 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 Vshaped segments can be observed on their own pressure derivative curves, respectively.
Effects of interregional diffusivity ratio and interregional conductivity ratio on transient responses are strong in the earlyflow 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 earlyflow period at the same production.
The shape and position of type curves are also dominated by storativity ratios and interporosity coefficients. In the lateflow 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.
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 middle and lateflow period. Due to a higher connectivity caused by acid fracturing, the pseudoskin 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 lateflow period. However, the ability of enhancing the well production decreases with the increase of interregional conductivity ratio.
Nomenclature
Field Variables:  Permeability, D 
:  Formation thickness, cm 
:  Fluid viscosity, cp 
:  Pressure, atm 
:  Radial distance, cm 
:  Halflength of artificial fracture or wormhole, cm 
:  Aperture of artificial fracture or wormhole, cm 
:  Time variable, s 
:  Porosity, fraction 
:  Compressibility factor, 1/atm 
:  Total rate at the wellbore, cm^{3}/s 
:  Volume factor, dimensionless 
:  Shape factor, dimensionless 
:  Fluid transmissivity factor, dimensionless 
:  Number of artificial fractures and wormholes, dimensionless 
:  3.1415926 
:  Dimensionless radius 
:  Dimensionless drainage radius 
:  Dimensionless wellbore radius 
:  Dimensionless grid step in the inner zone 
:  Dimensionless grid step in the outer zone 
:  Dimensionless time step in the inner zone 
:  Dimensionless time 
:  Dimensionless pressure 
:  Dimensionless pressure in the inner zone 
:  Dimensionless wellbore pressure 
:  Dimensionless pressure of natural fracture system 
:  Dimensionless pressure of vug system 
:  Dimensionless pressure of matrix system 
:  Dimensionless pressure derivative. 
:  Time variable in Laplace domain 
:  Relation used in Laplace space to distinguish triporosity properties 
:  Corresponding coefficient in Laplace domain 
:  Corresponding coefficient in Laplace domain 
:  Dimensionless pressure in Laplace domain 
:  Dimensionless pressure in Laplace domain 
:  Dimensionless wellbore pressure in Laplace domain 
:  Dimensionless pressure in Laplace domain 
:  Dimensionless pressure in Laplace domain 
:  Dimensionless pressure in Laplace domain. 
:  Interregional diffusivity ratio 
:  Interregional conductivity ratio 
:  Fluid storativity ratio 
:  Interporosity coefficient. 
:  Corresponding coefficient in the tridiagonal matrix 
:  Corresponding coefficient in the tridiagonal matrix 
:  Total number of discrete segments in the reservoir 
:  Number of discrete segments in the inner region. 
:  Modified Bessel function (1st kind, zero order) 
:  Modified Bessel function (1st kind, first order) 
:  Modified Bessel function (2nd kind, zero order) 
:  Modified Bessel function (2nd kind, first order) 
:  Exponential integral function 
:  Hyperbolic sine function 
:  Hyperbolic cosine function. 
:  Dimensionless 
:  Wellbore 
:  Reservoir external boundary 
:  The inner region 
:  Natural fracture system 
:  Vug system 
:  Matrix system 
:  Initial or ordinal. 
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge that this study was funded by the National Science and Technology Major Project of China (no. 2017ZX05030).
References
 J. Liu, G. S. Bodvarsson, and Y.S. Wu, “Analysis of flow behavior in fractured lithophysal reservoirs,” Journal of Contaminant Hydrology, vol. 6263, pp. 189–211, 2003. View at: Publisher Site  Google Scholar
 Y. S. Wu, G. Qin, R. E. Ewing, Y. Efendiev, Z. Kang, and Y. Ren, “A multiplecontinuum approach for modeling multiphase flow in naturally fractured vuggy petroleum reservoirs,” in Proceedings of the International Oil & Gas Conference and Exhibition in China, SPE 104173MS, p. 5, Beijing, China, 2006. View at: Google Scholar
 Y. Wu, C. A. EhligEconomides, G. Qin et al., “A triplecontinuum pressuretransient model for a naturally fractured vuggy reservoir,” in Proceedings of the SPE Annual Technical Conference and Exhibition, SPE 110044MS, Anaheim, Calif, USA, 2007. View at: Publisher Site  Google Scholar
 P. W. M. Cobett, S. Geiger, L. Borges, M. Garayev, and J. G. Gonzalez, “Limitations in the numerical well test modeling of fractured carbonate rocks,” in Proceedings of the SPE EUROPEC/EAGE Annual Conference and Exhibition, SPE 130252MS, Barcelona, Spain, 2010. View at: Publisher Site  Google Scholar
 A. F. Gulbransen, V. L. Hauge, and K.A. Lie, “A multiscale mixed finiteelement method for vuggy and naturally fractured reservoirs,” SPE Journal, vol. 15, no. 2, pp. 395–403, 2010. View at: Publisher Site  Google Scholar
 J. Yao, Z. Huang, Y. Li, C. Wang, and X. Lv, “Discrete fracturevug network model for modeling fluid flow in fractured vuggy porous media,” in Proceedings of the International Oil and Gas Conference and Exhibition in China, SPE130287MS, Beijing, China, 2010. View at: Publisher Site  Google Scholar
 J.C. Guo, R.S. Nie, and Y.L. Jia, “Dual permeability flow behavior for modeling horizontal well production in fracturedvuggy carbonate reservoirs,” Journal of Hydrology, vol. 464465, pp. 281–293, 2012. View at: Publisher Site  Google Scholar
 C. N. Fredd, “Dynamic model of wormhole formation demonstrates conditions for effective skin reduction during carbonate matrix acidizing,” in Proceedings of the SPE Permian Basin Oil and Gas Recovery Conference, PE59537MS, Midland, Tex, USA, 2000. View at: Publisher Site  Google Scholar
 K. Bybee, “Acid Fracturing a Carbonate Reservoir,” Journal of Petroleum Technology, vol. 56, no. 07, pp. 49–52, 2015. View at: Publisher Site  Google Scholar
 D. E. Nierode, B. R. Williams, and C. C. Bombardieri, “Prediction of stimulation from acid fracturing treatments,” Journal of Canadian Petroleum Technology, vol. 11, no. 04, pp. 30–41, 1972. View at: Publisher Site  Google Scholar
 K. D. Mahrer, W. W. Aud, and J. T. Hansen, “Farfield hydraulic fracture geometry: a changing paradigm,” in Proceedings of the SPE Annual Technical Conference and Exhibition, SPE36441MS, Denver, Colo, USA, 1996. View at: Publisher Site  Google Scholar
 Y.L. Zhao, B.C. Shan, L.H. Zhang, and Q.G. Liu, “Seepage flow behaviors of multistage fractured horizontal wells in arbitrary shaped shale gas reservoirs,” Journal of Geophysics and Engineering, vol. 13, no. 5, pp. 674–689, 2016. View at: Publisher Site  Google Scholar
 Y. Zhao, L. Zhang, Y. Xiong, Y. Zhou, Q. Liu, and D. Chen, “Pressure response and production performance for multifractured horizontal wells with complex seepage mechanism in boxshaped shale gas reservoir,” Journal of Natural Gas Science and Engineering, vol. 32, pp. 66–80, 2016. View at: Publisher Site  Google Scholar
 A. A. Daneshy, “Offbalance growth: a new concept in hydraulic fracturing,” Journal of Petroleum Technology, vol. 55, no. 4, pp. 78–85, 2003. View at: Publisher Site  Google Scholar
 Y. Wang and X. Yi, “Transient pressure behavior of a fractured vertical well with a finiteconductivity fracture in triple media carbonate reservoir,” Journal of Porous Media, vol. 20, no. 8, pp. 707–722, 2017. View at: Publisher Site  Google Scholar
 T. L. Blanton, “Propagation of hydraulically and dynamically induced fractures in naturally fractured reservoirs,” in Proceedings of the SPE Unconventional Gas Technology Symposium, SPE15261MS, Louisville, Ky, USA, 1986. View at: Publisher Site  Google Scholar
 H. D. Murphy and M. C. Fehler, “Hydraulic fracturing of jointed formations,” in Proceedings of the International Meeting on Petroleum Engineering, SPE14088MS, Beijing, China, 1986. View at: Publisher Site  Google Scholar
 X. Zhu, J. Gibson, N. Ravindran, R. Zinno, and D. Sixta, “Seismic imaging of hydraulic fractures in carthage tight sands: a pilot study,” Leading Edge, vol. 15, no. 3, pp. 218–224, 1996. View at: Publisher Site  Google Scholar
 M. K. Fisher, C. A. Wright, B. M. Davidson et al., “Integrating fracture mapping technologies to optimize stimulations in the barnett shale,” in Proceedings of the SPE Annual Technical Conference and Exhibition, SPE77441MS, San Antonio, Tex, USA, 2002. View at: Publisher Site  Google Scholar
 M. J. Mayerhofer, E. P. Lolon, J. E. Youngblood, and J. R. Heinze, “Integration of MicroseismicFractureMapping Results With Numerical Fracture Network Production Modeling in the Barnett Shale,” in Proceedings of the SPE Annual Technical Conference and Exhibition, SPE102103MS, San Antonio, Tex, USA, 2006. View at: Publisher Site  Google Scholar
 K. M. Hung, A. D. Hill, and K. Sepehrnoori, “Mechanistic model of wormhole growth in carbonate matrix acidizing and acid fracturing,” Journal of Petroleum Technology, vol. 41, no. 1, pp. 59–66, 1989. View at: Publisher Site  Google Scholar
 Y. K. AlDuailej, H. T. Kwak, S. Caliskan, and I. S. AlYami, “Wormhole Characterisation Using NMR,” in Proceedings of the International Petroleum Technology Conference, Beijing, China, 2013. View at: Publisher Site  Google Scholar
 M. Liu, S.C. Zhang, J.Y. Mou, and F.J. Zhou, “Wormhole propagation behavior under reservoir condition in carbonate acidizing,” Transport in Porous Media, vol. 96, no. 1, pp. 203–220, 2013. View at: Publisher Site  Google Scholar
 Z. M. Ahmed, M. Mofti, E. Fidan et al., “Complex hydraulic acid fracturing field application in appraisal and development of tight, sour, overpressured, and tectonically active, HPHT carbonates in North Kuwait,” in Proceedings of the Abu Dhabi International Petroleum Exhibition & Conference, SPE182930MS, Abu Dhabi, UAE, 2016. View at: Google Scholar
 P. Guizada, “Competent utilization of multistage acid fracturing technology to increment well productivity during stimulation of low permeability gas carbonate reservoirs – Case study,” in Proceedings of the SPE Russian Petroleum Technology Conference and Exhibition, SPE182116MS, Moscow, Russia, 2016. View at: Publisher Site  Google Scholar
 D. C. C. Poon, Pressure Transient Analysis of a Composite Reservoir with Uniform Fracture Distribution, SPE13384MS, 1984, https://www.onepetro.org/general/SPE13384MS.
 J. S. Olarewaju and W. J. Lee, “A comprehensive application of a composite reservoir model to pressure transient analysis,” in Proceedings of the SPE California Regional Meeting, SPE16345MS, Ventura, Calif, USA, 1987. View at: Publisher Site  Google Scholar
 J. F. Stanislav, G. V. Easwaran, and S. L. Kokal, “Analytical solutions for vertical fractures in a composite system,” Journal of Canadian Petroleum Technology, vol. 26, no. 5, pp. 50–56, 2013. View at: Publisher Site  Google Scholar
 J. E. Warren and P. J. Root, “The behavior of naturally fractured reservoirs,” SPE Journal, vol. 3, no. 3, pp. 245–255, 2013. View at: Publisher Site  Google Scholar
 H. Stehfest, “Numerical inversion of Laplace transforms,” Communications of the ACM, vol. 13, no. 1, pp. 47–49, 1970. View at: Publisher Site  Google Scholar
 C. P. Chiang and W. A. Kennedy, “Numerical simulation of pressure behavior in a fractured reservoir,” in Proceedings of the Fall Meeting of the Society of Petroleum Engineers of AIME, SPE3067MS, Houston, Tex, USA, 1970. View at: Publisher Site  Google Scholar
 J. Douglas and D. W. Peaceman, “Numerical solution of two‐dimensional heat‐flow problems,” AIChE Journal, vol. 1, no. 4, pp. 505–512, 1955. View at: Publisher Site  Google Scholar
 A. F. Van Everdingen and W. Hurst, “The application of the Laplace transformation to flow problems in reservoirs,” Journal of Petroleum Technology, vol. 1, no. 12, pp. 305–324, 1949. View at: Publisher Site  Google Scholar
 E. Ozkan and R. Raghavan, “New solutions for welltestanalysis problems: Part 1: Analytical considerations,” SPE Formation Evaluation, vol. 6, no. 3, pp. 359–368, 1991. View at: Publisher Site  Google Scholar
 E. Ozkan and R. Raghavan, “New solutions for welltestanalysis problems: Part 2: Computational considerations and applications,” SPE Formation Evaluation, vol. 6, no. 3, pp. 369–378, 1991. View at: Publisher Site  Google Scholar
 R.S. Nie, Y.F. Meng, J.C. Guo, and Y.L. Jia, “Modeling transient flow behavior of a horizontal well in a coal seam,” International Journal of Coal Geology, vol. 92, pp. 54–68, 2012. View at: Publisher Site  Google Scholar
 L. Wang, X. Wang, J. Li, and J. Wang, “Simulation of pressure transient behavior for asymmetrically finiteconductivity fractured wells in coal reservoirs,” Transport in Porous Media, vol. 97, no. 3, pp. 353–372, 2013. View at: Publisher Site  Google Scholar
 A. C. Gringarten, H. J. Ramey, and Raghavan. R., “Unsteadystate pressure distributions created by a well with a single infiniteconductivity vertical fracture,” SPE Journal, vol. 14, no. 4, pp. 347–360, 1974. View at: Publisher Site  Google Scholar
 H. CincoLey, Unsteadystate pressure distributions created by a slanted well, or a well with an inclined fracture, [PhD Thesis], Stanford University, 1974.
 A. V. Dinh and D. Tiab, “Transientpressure analysis of a well with an inclined hydraulic fracture using type curve matching,” SPE Reservoir Evaluation & Engineering, vol. 13, no. 06, pp. 845–860, 2013. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Mingxian Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.