Abstract
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.
1. 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–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–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–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–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 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–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.
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 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.
(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 single-phase 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 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 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 pseudo-steady 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 block-centered 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 no-flow 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 pseudo-skin 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 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 toCorresponding 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 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 V-shaped 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 V-shaped 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 V-shaped 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 V-shaped 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) Pseudo-Steady Flow Regime. For closed boundary, the long-time behavior corresponds to a pseudo-steady 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 pseudo-steady 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 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.
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 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.
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 pressure-transient response (, 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.
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 log-log 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 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.
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 late-flow periods. Fracture storativity ratio begins to affect the transient response at the start of the first interporosity flow, which extends to the pseudo-state 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 V-shaped segments become shallower and even the second V-shaped 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 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 middle-flow period, and the pressure and derivative curves for different matrix storativity ratios overlap completely in the early- and 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 V-shaped 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, = 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.
6.2.4. 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 V-shaped 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 ( = 0.001, 0.005, 0.025, and 0.05). This coefficient also mainly affects middle-flow 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 V-shaped 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 ( = 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.
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 middle-flow 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 , , , , , , , 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.
6.4. 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] 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 pseudo-skin 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.
7. Conclusions
This work presented a composite reservoir model to analyze the flow behavior of an acid fracturing well in a fractured-vuggy 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 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.
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.
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.
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 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.
Nomenclature
Field Variables| : | Permeability, D |
| : | Formation thickness, cm |
| : | Fluid viscosity, cp |
| : | Pressure, atm |
| : | Radial distance, cm |
| : | Half-length 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, cm3/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).