Abstract
Triple-porosity model is usually adopted to describe reservoirs with multiscaled pore spaces, including matrix pores, natural fractures, and vugs. Multiple fractures created by hydraulic fracturing can effectively improve the connectivity between existing natural fractures and thus increase well deliverability. However, little work has been done on pressure transient behavior of multistage fractured horizontal wells in triple-porosity reservoirs. Based on source/sink function method, this paper presents a triple-porosity model to investigate the transient pressure dynamics and flux distribution for multistage fractured horizontal wells in fractured-vuggy reservoirs with consideration of stress-dependent natural fracture permeability. The model is semianalytically solved by discretizing hydraulic fractures and Pedrosa’s transformation, perturbation theory, and integration transformation method. Type curves of transient pressure dynamics are generated, and flux distribution among hydraulic fractures for a fractured horizontal well with constant production rate is also discussed. Parametric study shows that major influential parameters on transient pressure responses are parameters pertinent to reservoir properties, interporosity mass transfer, and hydraulic fractures. Analysis of flux distribution indicates that flux density gradually increases from the horizontal wellbore to fracture tips, and the flux contribution of outermost fractures is higher than that of inner fractures. The model can also be extended to optimize hydraulic fracture parameters.
1. Introduction
The concept of dual-porosity media was first proposed by Barenblatt et al. [1] to describe naturally fractured media. Warren and Root [2] then extended this model to analyze pressure transient dynamics for vertical wells in naturally fractured oil reservoirs. Recently, Cai et al. [3, 4] analyzed the imbibition mechanism in fractured-porous media based on fractal geometry. Although the dual-porosity model can describe most naturally fractured reservoirs, there are still some drawbacks when applying the dual-porosity model to describe reservoirs with multiple pore media, such as reservoirs with natural fractures and vugs (different from matrix pores).
In order to better describe reservoir heterogeneity in this type of naturally fractured-vuggy reservoirs, the concept of triple-porosity model was proposed. The first triple-porosity model in the literature was proposed by Abadassah and Ershaghi [5], in which the reservoir was represented by a combination of two matrix systems with different properties and uniformly distributed natural fractures. The model was then further extended by Al-Ghamdi and Ershaghi [6] to represent reservoirs with two fracture systems (microfractures and macrofractures) and one matrix system. After that, over the past several decades, numerous efforts have been made to investigate transient pressure behavior in so-called triple-porosity reservoirs. However, almost all the research focuses on vertical wells [7–19] or horizontal wells [20–22] in triple-porosity reservoirs, and few studies have been done on complex well-reservoir configuration, such as multistage fractured horizontal well in triple-porosity reservoirs.
Field practices have proven that horizontal wellbore combining multistage hydraulic fracturing can effectively enhance the permeability in the vicinity of horizontal wellbore. Almost all the triple-porosity models [23, 24] for fractured horizontal wells assume linear flow, which cannot completely characterize the pressure responses during all flowing periods (such as pseudoradial flow period). On the other hand, in most existing models for fractured horizontal wells in triple-porosity reservoirs three porous media denote hydraulic fractures, natural fractures, and matrix pores, which means the reservoir itself is actually considered as a dual-porosity medium.
This study develops a semianalytical model based on source/sink function theory for analyzing transient pressure responses and flux distribution of naturally fractured-vuggy reservoirs, in which the reservoir itself is conceptualized as triple-porosity media and the interaction between hydraulic fractures and reservoir is considered. The model presented here can completely reflect transient pressure characteristics during all the possible flowing period as well as flux distribution among multiple fractures. In addition, this model also takes into account the stress-sensitivity of permeability caused by closure of natural fractures during production.
2. Physical Model of a Multistage Fractured Horizontal Well in a Triple-Porosity Reservoir
As shown in Figure 1, a horizontal well with multiple hydraulic fractures, producing at a constant rate () in a reservoir with multiscale storage spaces, is considered here. Both the cap rock and underlying rock formation of the reservoir are usually shales with extremely low permeability, which can be conceived as no-flow boundaries. The reservoir thickness and oil viscosity are assumed to be and , separately. Other basic assumptions involved in the derivation of the triple-porosity model are as follows:(1)The reservoir is assumed to be composed of three contiguous porous media, which are matrix, natural fractures, and vugs (representing larger voids in the reservoir which are not part of natural fractures).(2)Natural fractures are assumed to be directly connected to hydraulic fractures or the horizontal wellbore, while the fluid stored in matrix pores or vugs does not have direct access to hydraulic fractures or the horizontal wellbore. Figure 2 presents a simple illustration of the fluid flowing path in triple-porosity reservoirs.(3)The permeability of natural fractures is considered stress-sensitive to incorporate the effect of partial or complete closure of natural fractures during production.(4)Both oil and rock are considered slightly compressible, and the reservoir has uniform pressure distribution at time zero.(5)Consider single-phase isothermal flow and negligible gravity and capillary effects.
To obtain the pressure responses caused by the well-reservoir configuration shown in Figure 1, a mathematical model together with corresponding boundary and initial conditions is required. In the problem addressed here, the multiple fractures inner boundary condition is so complex that it is rather difficult to directly describe it with a mathematical equation. In order to solve this problem, source function theory, a powerful tool to solve complex transient flow problems in reservoirs, is adopted in this work. We first studied the transient pressure responses caused by a continuous line-sink in triple-porosity reservoirs and then adopted the line-sink solution with superposition principle to obtain the pressure transient dynamics and flux distribution of multistage fractured horizontal wells in triple-porosity reservoirs.
3. Line-Sink Model in Stress-Sensitive Triple-Porosity Reservoirs
3.1. Physical Model for a Line-Sink in Triple-Porosity Reservoirs
As shown in Figure 3, a continuous line-sink in a triple-porosity reservoir with stress-dependent fracture permeability is considered in this section. The strength of the line-sink is represented by .
3.2. Mathematical Model for a Line-Sink in Triple-Porosity Reservoirs
(1) Governing Equations. Imposing mass-balance equation over a control volume in the triple-porosity reservoir and substituting equation of state and equation of motion, we can get the governing equations for matrix, natural fracture, and vug systems.
For natural fracture system with stress-dependent permeability, the governing equation is as follows (i.e., (A.26b) in Appendix A):
For matrix system, the governing equation is as follows (i.e., (A.30) in Appendix A):
For vug system, the governing equation is as follows (i.e., (A.34) in Appendix A):(2) Boundary Conditions. According to the derivation in Appendix B, the inner boundary condition of the line-sink model is (i.e., (B.5) in Appendix B)
Assuming infinitely large reservoirs, the outer boundary condition can be expressed as(3) Initial Condition. The initial pressure distribution in the reservoir is assumed to be uniform, which is
3.3. Dimensionless Form of the Line-Sink Model
For the purpose of simplifying derivation and comparing between different reservoirs, the mathematical model is derived and solved in dimensionless form in this work. The pressure drops involved in the mathematical model are , , and , and relevant dimensionless variables used in the mathematical model are listed in Table 1.
Obviously, , , and in Table 1 satisfy the following equation:
With the definitions in Table 1, the dimensionless forms of (1)~(6) can be obtained as follows:
3.4. Solution of the Dimensionless Line-Sink Model
3.4.1. Pedrosa’s Linearization
The stress-sensitivity of natural fracture permeability makes the above seepage model strongly nonlinear which cannot be solved analytically. Here, following Pedrosa Jr. [25], we introduced the following equation to alleviate the nonlinearity of (8) and (11):where is an intermediate dimensionless parameter which is also a function of radial distance and time.
3.4.2. Perturbation Technique
According to the perturbation theory, the terms , , and in (15) through (20) can be expanded as power series in dimensionless permeability modulus, and we can get the following equations:
Given that the dimensionless permeability modulus is usually small (much smaller than 1), the zero-order approximate solution can satisfy the requirements of engineering precision, so (15)~(20) become
Equation (22) is the linearized line-sink model in triple-porosity reservoirs.
3.4.3. Laplace Transformation
Take the following Laplace transformation:where is the Laplace transformation variable.
Then (22) becomes
3.4.4. Line-Sink Solution
Equations (25) and (26) can be rewritten as
Substituting (29) into (24) yields
If we definethen (30) can be expressed as follows:
The general solution of (32) can be easily obtained as follows:where and are Bessel functions and and are unknown coefficients which can be determined by corresponding boundary conditions.
If the line-sink is located at the origin of coordinates (center of the triple-porosity reservoir), then the dimensionless radial distance in (31) can be calculated by the following equation:
If the line-sink is located at instead, then is calculated by the following equation:where
Substitution of (33) into the inner boundary condition (27) yields
According to the properties of Bessel’s function, which are and , (39) becomes
Similarly, with the outer boundary condition given in (28) and the properties of Bessel’s function, and , the coefficient in (33), , should satisfy the following equation:
With (40) and (41), we can get the final form of (33) as follows:
Equation (42) is the basic zero-order perturbation solution for a line-sink in a stress-sensitive triple-porosity reservoir. As mentioned above, zero-order perturbation solution is enough to approximate the exact solution of (15) to (20); that is, . With the definition of in Table 1, (42) can be expressed as
Equation (43) can be rewritten aswhere
Equation (44) is the basic line-sink solution for stress-sensitive triple-porosity reservoirs. With the basic line-sink solution and the superposition principle, we can obtain the pressure responses caused by multistage fractured horizontal wells in stress-sensitive triple-porosity reservoirs.
4. Pressure Responses for Multistage Fractured Horizontal Wells in Stress-Sensitive Triple-Porosity Reservoirs
Figure 4 shows the schematic of a horizontal well with multiple (a total of ) hydraulic fractures in a stress-sensitive triple-porosity reservoir. The -axis is set to be aligned with the horizontal wellbore. The intersection of the th hydraulic fracture and the -axis is represented by , and the distance between every two adjacent fractures is .
To successfully obtain the pressure responses caused by the multiple fractures with the line-sink solution derived in the above section, it is necessary to discretize the multiple hydraulic fractures. As shown in Figure 4, each hydraulic fracture is divided into segments along its length. The coordinate of the midpoint (i.e., discrete node) of the th segment on the th fracture is denoted by , and the coordinates of the two corresponding endpoints are denoted by and .
With the discretization scheme shown in Figure 4, the coordinates of endpoints of discrete segments can be calculated by the following equations:
The coordinates of discrete nodes (midpoint of each segment) can be calculated by the following equations:where is the length of the left wing of th fracture, m; is the length of the right wing of th fracture, m. In the model presented in this paper, and can be different and can vary from fracture to fracture.
The distance between and is
It is obvious that the flux strength is different along the fracture length. However, in the same discrete segment (assuming the discrete segment is small enough), the flux strength can be considered as constant. If we use to represent the flux density per unit length, then, with the superposition principle, the pressure response at in the reservoir caused by the discrete segment can be obtained by integrating the basic line-sink along the discrete segment, which is
With the definitions of and , that is, (35) and (37), (49) can be changed into
If we define dimensionless flux density per unit length as and take Laplace transformation of (51),
Equation (50) can be rewritten as
The pressure response at caused by segments can be obtained by applying the principle of superposition over all discrete fracture segments:
Thus, the pressure response at the discrete segment can be obtained as follows:
Since the permeability of hydraulic fractures is always much higher than the original reservoir permeability, the pressure drop in hydraulic fractures is much smaller than the pressure drop caused by oil flow in the reservoir. Thus, the hydraulic fractures can be considered as infinitely conductive, meaning the pressure in hydraulic fractures is equal to the bottom-hole pressure:
Combining (55) and (56), we can get
The subscripts and in (57) denote that this equation is written for discrete segment . By writing (57) for all discrete segments, that is, letting , we can get equations. It should be noted that there are unknowns in these equations, which are and . To compose a closed system of equations, the assumption of constant total flow rate can give the following equation: where is the length of the discrete segment () and can be calculated by
Taking Laplace transformation of (58), one can get
The dimensionless form of (60) iswhere
Equations (57) and (61) represent a system of () linear algebraic equations relating () unknowns, which can be solved by direct methods, such as Gaussian elimination method, Gauss-Jordan reduction method.
calculated by (57) and (61) does not take into account the effects of wellbore storage and skin. According to van Everdingen [26] and Kucuk and Ayestaran [27], the effects of wellbore storage and skin can be incorporated in the calculated bottom-hole pressure response bywhere is the dimensionless wellbore storage coefficient and is the skin factor defined by the following, respectively: where is the wellbore storage coefficient, m3/Pa; is the sandface flow rate of the multistage fractured horizontal well in the triple-porosity reservoir, m3/s; is an extra pressure drop near the wellbore, Pa.
in (63) is the bottom-hole pressure response obtained in the Laplace domain, and by Stehfest [28] numerical inversion algorithm we can calculate the pressure responses in real time domain. Then, with (66), we can obtain the bottom-hole pressure response for multistage fractured horizontal wells in triple-porosity reservoirs incorporating the stress-sensitivity of natural fracture system:
5. Results and Discussion
In this section, a horizontal well with three hydraulic fractures (i.e., ) in a triple-porosity reservoir is considered, and the model introduced above is applied to investigate the dimensionless pressure responses of this well-reservoir configuration.
5.1. Flow Periods and Effect of Stress-Sensitivity
Figure 5 presents the dimensionless pressure responses and corresponding pressure derivatives of a multistage fractured horizontal well in an infinitely large triple-porosity reservoir with consideration of stress-sensitivity of natural fractures. It can be observed from the figure that as many as ten flow periods can be identified at different time scales. Different flow periods and their corresponding characteristics are explained as follows.
Period 1. It is the early-time wellbore storage period. In this period, both pressure curves and their derivative curves exhibit unit slope on the log-log plots.
Period 2. It is the transition flow period after wellbore storage period. The pressure derivative curves exhibit a “hump,” representing the effect of skin.
Period 3. It is the first linear flow period. This period corresponds to the linear flow in the direction perpendicular to hydraulic fracture length as illustrated in Figure 6(a). The pressure derivative curves exhibit a straight line with slope equal to “1/2.”
(a)
(b)
(c)
(d)
Period 4. It is the first pseudoradial flow period. This period can occur when the distances between adjacent fractures are relatively large. Oil flows in the reservoir in a way as shown in Figure 6(b), and pseudoradial flow occurs around each fracture. This period corresponds to a nearly horizontal line with a value of “” in the derivative curve, where is the number of hydraulic fractures.
Period 5. It is the second linear flow period. This period corresponds to linear flow perpendicular to the horizontal wellbore as shown in Figure 6(c). Another straight line with a slope of about “1/2” can be observed in the derivative curve. When the permeability of natural fracture system remains constant during the whole production, corresponding to the no stress-sensitivity case, the slope of the derivative curves during this period is “1/2.” When taking into account the stress-sensitivity of natural fracture system, the corresponding pressure derivative curve slightly deviates from the one-half-slope line and exhibits a slight upward tendency.
Period 6. It is the second pseudoradial flow period, corresponding to radial flow in natural fracture system. During this period, multiple hydraulic fractures and the horizontal wellbore behave like an enlarged wellbore. The derivative curves exhibit a horizontal line with an upward tendency depending on the value of dimensionless permeability modulus. In the case that , the derivative curves are horizontal with a value of 0.5 on the -axis. With the increase of the value of , the derivative curves turn upward gradually.
Period 7. It is the first interporosity flow period. This period reflects mass transfer between vugs and natural fracture system. Due to the oil flowing from vugs to natural fractures, the bottom-hole pressure drops slower and an obvious “dip” can be observed in the derivative curves.
Period 8. It is the third pseudoradial flow period, corresponding to compound radial flow in natural fracture system and vugs. Due to the existence of stress-sensitivity of natural fracture system, the pressure derivative curves are no longer horizontal. The position of the derivative curves becomes higher with larger value of .
Period 9. It is the second interporosity flow period. This period reflects oil flowing from matrix to natural fracture system. Another characteristic “dip” can be observed in the derivative curves.
Period 10. It is the late-time compound pseudoradial flow period. This period reflects compound radial flow in natural fracture system, matrix, and vugs as shown in Figure 6(d). Similarly, the derivative curves exhibit upward tendencies due to the existence of stress-sensitivity of natural fracture system.
5.2. Effect of Skin Factor
Figure 7 shows the effect of skin factor on dimensionless pressure curves and corresponding derivative curves. As stated above, the skin factor mainly affects the shape of type curves during the transition period (Period 2) after wellbore storage period. With the increase of the value of skin factor, the position of the “hump” in the derivative curves becomes higher, representing larger pressure drop in the formation.
5.3. Effect of Storativity Ratio of Natural Fracture System
Figure 8 shows the effect of storativity ratio of natural fracture system on dimensionless pressure and pressure derivative curves. It is obvious in Figure 8 that when the storativity ratio of matrix (i.e., ) and other parameters are kept constant, the storativity ratio of natural fracture system mainly affects the transient pressure dynamics during the first linear flow period, the first pseudoradial flow period, the second linear flow period, second pseudoradial flow period, and the first interporosity flow period. On the one hand, the position of the dimensionless pressure curves becomes higher during the abovementioned flowing periods with the decrease of . On the other hand, with the decrease of , the position of the pressure derivative curves during the first linear flow period, the first pseudoradial flow, and the second linear flow period becomes higher, and the first “dip” in pressure derivative curves becomes wider and deeper.
5.4. Effect of Storativity Ratio of Matrix
Figure 9 shows the effect of storativity ratio of matrix on dimensionless pressure and derivative curves. It can be observed that, with fixed , the storativity ratio of matrix has primary effect on pressure dynamics during two interporosity flowing periods. With the decrease of the value of , the fist “dip” in the pressure derivative curves becomes wider and deeper, while the second “dip” in the pressure derivative curves exhibits the opposite trend.
5.5. Effect of Interporosity Flow Coefficient between Vugs and Natural Fractures
Figure 10 shows the effect of interporosity flow coefficient between vugs and natural fractures () on dimensionless pressure and pressure derivative curves. As expected, mainly affects the occurrence time of the first interporosity flow period, that is, the appearance time of the first “dip” in pressure derivative curves. The smaller the value of is, the later the first “dip” can be observed in the pressure derivative curves. When decreases to a certain extent, the third pseudoradial flow period (Period 8 in Figure 5) may not be observed in the type curves. In addition, when increases to a certain extent, the second pseudoradial flow period (Period 6 in Figure 5) and even the second linear flow period (Period 5 in Figure 5) may be masked by the following first interporosity flow period.
5.6. Effect of Interporosity Flow Coefficient between Matrix and Natural Fractures
Figure 11 shows the effect of interporosity flow coefficient between matrix and natural fractures () on dimensionless pressure and pressure derivative curves. It can be seen that mainly affects the occurrence time of the second interporosity flow period, that is, the appearance time of the second “dip” in pressure derivative curves. The smaller the value of is, the later the second “dip” can be observed in the pressure derivative curves. With the increase of , the reflection of the third pseudoradial flow period on the derivative curves (Period 8 in Figure 5) may not be observed. When increases to a certain degree, the second “dip” may merge with the first “dip,” meaning that simultaneous interporosity mass transfers happen between natural fractures and matrix as well as natural fractures and vugs.
5.7. Effect of Length of Hydraulic Fractures
Figure 12 presents the effect of half-length of hydraulic fractures on type curves. For the convenience of discussion, is assumed here. It can be observed that the length of hydraulic fracture has primary effect on the first linear flow period and the subsequent first pseudoradial flow period. With the increase of fracture length, the first linear flow period lasts longer with lower position of pressure derivative curves, and the first pseudoradial flow period occurs later with shorter duration. When the half-length of hydraulic fractures continues to increase, the horizontal line in the pressure derivative curves reflecting the first pseudoradial flow period around each fracture will gradually disappear.
5.8. Effect of Hydraulic Fracture Spacing
Figure 13 shows the effect of the distance between adjacent hydraulic fractures. The case discussed in Figure 13 assumes equal fracture spacing for the convenience of discussion. It can be observed that the distance between adjacent hydraulic fractures mainly affects the first pseudoradial flow period and the subsequent second linear flow period. The larger the distance between adjacent fractures, the longer the duration of the first pseudoradial flow period, and the later the occurrence of the second linear flow period. Another point should be addressed here is that when increases to a certain degree, the second pseudoradial flow period may not be observed in the type curves; that is, the first interporosity flow period appears after the second linear flow period.
5.9. Flux Distribution
Figure 14 presents the dimensionless flux distribution along hydraulic fractures at a given time based on the calculation results obtained with the model proposed in this paper. In the case presented in Figure 14, three hydraulic fractures with equal half-length and fracture spacing are considered, and each hydraulic fracture is discretized into ten segments. It can be observed that, for the same hydraulic fracture, the flux distribution is symmetrical with respect to the horizontal wellbore, and the flux density at fracture tips is larger than that in the middle of the fracture. In addition, for different hydraulic fractures, it can be found that the flux contribution from hydraulic fracture in the middle (Fracture 3 in Figure 14) is always smaller than that from outer fractures (Fractures 1 and 2 in Figure 14).
6. Conclusion
Based on the mathematical model and discussion in this paper, the following conclusions can be warranted:(1)A mathematical model is presented for investigating transient pressure dynamics as well as flux distribution of multistage fractured horizontal wells in stress-sensitive triple-porosity reservoirs. The model can better describe fluid flow in naturally fractured reservoirs with multiple types of storage spaces and complex well type.(2)Semianalytical solution is developed in the Laplace domain by the integral transformation method and discretization of hydraulic fractures. With the semianalytical solution and the numerical inversion algorithm proposed by Stehfest, dimensionless pressure responses and corresponding pressure derivatives as well as flux distribution among all the fractures can be obtained in real time domain. Computation results prove that the model presented in this paper is stable and accurate.(3)Analysis of transient pressure and derivative responses indicates that as many as ten flow periods can be identified for a multistage fractured horizontal well in triple-porosity reservoirs. Major factors affecting the transient pressure responses include skin factor, storativity ratios of different porous media type, interporosity flow coefficients, and parameters related to the geometry and placement of hydraulic fractures. With different combinations of these parameters, the reflection of some flow periods may be masked in the pressure derivative curves.(4)Two characteristic “dips” can be observed in the pressure derivative curves for triple-porosity reservoirs, reflecting the mass transfer between natural fracture-vugs and natural fracture-matrix, separately.(5)The stress-sensitivity of natural fracture system results in larger pressure drop during intermediate and late flowing periods which is reflected by upward tendencies in both dimensionless pressure curves and corresponding derivative curves.(6)Analysis of flux distribution among multiple hydraulic fractures indicates that the flux contribution from outer fractures is higher, and for a given fracture the flux contribution from fracture tips is higher.(7)The presented model can be used to interpret pressure signal for fractured horizontal wells in triple-porosity reservoirs and provide some important dynamic parameters for reservoir development.
Appendices
A. Derivation of Governing Equations for Natural Fracture System, Matrix System, and Vug System in Triple-Porosity Reservoirs
Governing equations for natural fracture system, matrix system, and vug system in triple-porosity reservoirs can be obtained by combining equation of motion, equation of state, and mass conservation equation.
(1) Equation of Motion. Assuming radial flow in triple-porosity reservoirs, the oil flow velocity in natural fracture system is given by where is the radial oil velocity in natural fracture system, m/s; is the permeability of natural fracture system under pressure , m2; is the oil viscosity, Pa·s; is the pressure of natural fracture system, Pa; is the radial distance, m.
During the production process, the reservoir pressure gradually decreases and the effective stress increases. Consequently, the aperture of natural fractures gradually decreases, resulting in the reduction in the permeability of natural fracture system, which is so-called stress-sensitivity of permeability of natural fracture system. To account for this, a stress-dependent permeability is adopted in the natural fracture flow model.
Following Kikani and Pedrosa [29], the permeability modulus , which is used to describe the stress-sensitivity of natural fracture system, is defined bywhere is the permeability modulus, Pa−1.
Equation (A.2) can be further written aswhere is the permeability of natural fracture system under initial condition, m2; is the reservoir pressure under initial condition, Pa.
Solving (A.3), one can get
Equation (A.4) is one of the most common relationships which are used to describe stress-dependent permeability. This exponential relationship corresponds to Type I rocks discussed by Raghavan and Chin [30].
Substituting (A.4) into (A.1) yields the equation of motion in natural fracture system with consideration of stress-sensitivity effect as follows:(2) Equation of State
(a) Equation of State for Oil. Oil in natural fracture system can be described by the following equation of state:where is the oil density in natural fracture system, kg/m3; is the oil density at the reference pressure , kg/m3; is the oil compressibility, Pa−1.
The term in (A.6a) can be expanded into Maclaurin series. Considering that , the second and higher orders of the expansion can be neglected; thus we can get
Similarly, oil in matrix system can be described by the following equation of state:orwhere is the pressure of matrix system, Pa.
Oil in vug system can be described by the following equation of state:orwhere is the pressure of vug system, Pa.
(b) Equation of State for Porous Media. The equation of state for natural fracture system can be written aswhere is the porosity of natural fracture system, fraction; is the porosity of natural fracture system at the reference pressure , fraction; is the compressibility of natural fracture system, Pa−1.
Maclaurin series can also be adopted to approximate (A.9a). Considering that , the second and higher order terms can be neglected, and one can get
By following the same procedure, one can get the equation of state for matrix system as follows:orwhere is the porosity of matrix system at the reference pressure , fraction; is the compressibility of matrix system, Pa−1.
Similarly, the equation of state for vug system can be described as follows:orwhere is the porosity of vug system at the reference pressure , fraction; is the compressibility of matrix system, Pa−1.
(3) Mass Conservation Equation. By applying the mass conservation law on a representative elemental volume in a triple-porosity reservoir, we can get the following equation for natural fracture system:
Equation (A.12a) can also be written in the following equivalent form:where is the mass flow rate between matrix system and fracture system per unit-volume reservoir, kg/(m3·s), and it can be expressed by (A.13); is the mass flow rate between vug system and fracture system per unit-volume reservoir, kg/(m3·s), and it can be expressed by (A.14)where is the shape factor of matrix blocks, m−2; is the shape factor of vug system, m−2; is the permeability of matrix system, m2; is the permeability of vug system, m2.
By using the same method, we can get the following mass conservation equations for matrix system and vug system:
Here, following the assumption adopted in the dual-porosity model proposed by Warren and Root [2], natural fractures are assumed to be oil flowing channels, and matrix and vugs mainly serve as storage space; thus the flowing terms in (A.15) and (A.16) can be reduced to
With (A.17), (A.15) and (A.16) can be rewritten as
Equations (A.12a), (A.12b), (A.18), and (A.19) are mass conservation equations for natural fracture system, matrix system, and vug system in triple-porosity reservoirs, separately.
(4) Governing Equation for Natural Fracture System. Substitution of (A.15) into the mass conservation equation for natural fracture system, that is, (A.12b), yields
With (A.6a), (A.6b), and the transformation formula , the first term in the left-hand side of (A.20) can be written as
Following the same procedure, the second term in the left-hand side of (A.20) changes into
With (A.6b) and (A.9b), the product term within the bracket in the right-hand side of (A.20) can be expanded as follows:
Given that both and are very small terms, their product will yield a much smaller term which is negligible (), and (A.23a) can be simplified as
Thus the first term in the right-hand side of (A.20) can be written aswhere is the total compressibility of natural fracture system, including the compressibility of oil and natural fractures, and
Substituting (A.13), (A.14), (A.21), (A.22), and (A.24) into (A.20), one can get
After deriving governing equations for matrix system and vug system, (A.26a) can be further written as
Equation (A.26b) is the final form of governing equation for natural fracture system incorporating the pressure-dependent permeability, that is, (1) in the main body of this paper.
(5) Governing Equation for Matrix System. With (A.7b) and (A.10b), the product of oil density and porosity of matrix system can be written as
Analogously, because both and are very small terms, their product can be neglected, and (A.27a) can be simplified to obtain
Therefore,where is the total compressibility of matrix system, including the compressibility of oil and matrix, and
Substitution of (A.13) and (A.28) into the mass conservation equation for matrix system, that is, (A.18), yields
Equation (A.30) is the final form of governing equation for matrix system in triple-porosity reservoirs, that is, (2) in the main body of this paper.
(6) Governing Equation for Vug System. The derivation of governing equation for vug system is similar to that applied to matrix system. With (A.8b) and (A.11b), the product of oil density and porosity of vug system can be written as
Given that both and are much smaller than 1, their product is negligible, and (A.31a) reduces to
Therefore,where is the total compressibility of vug system, including the compressibility of oil and vug, and
Substituting (A.14) and (A.33) into mass conservation equation for vug system, that is, (A.19), one can get
Equation (A.34) is the final form of governing equation for vug system in triple-porosity reservoirs, that is, (3) in the main body of this paper.
B. Derivation of Inner Boundary Condition for the Line-Sink
A vertical line-sink can be treated as the limiting case of a vertical cylinder-sink with radius . The flow rate across the cylinder wall with radius can be expressed aswhere is the lateral area of the cylinder-sink, m2, which can be calculated by (B.2); is oil velocity at the cylinder wall, m/s, which can be calculated by (B.3); is flow rate of the continuous line-sink under standard condition, m3/s; is the oil formation volume factor, dimensionless:
Substitution of (A.4) into (B.3) yields
Substituting (B.2) and (B.4) into (B.1), we can get
Equation (B.5) is the inner boundary condition for a line-sink in triple-porosity reservoirs, that is, (4) in the main body of this paper.
Nomenclature
| Oil formation volume factor, m3/sm3 | |
| : | Wellbore storage coefficient, m3/Pa |
| : | Dimensionless wellbore storage coefficient, dimensionless |
| : | Oil compressibility, Pa−1 |
| : | Compressibility of natural fracture system, Pa−1 |
| : | Compressibility of matrix system, Pa−1 |
| : | Compressibility of vug system, Pa−1 |
| : | Total compressibility of natural fracture system, Pa−1 |
| : | Total compressibility of matrix system, Pa−1 |
| : | Total compressibility of vug system, Pa−1 |
| : | Reservoir thickness, m |
| : | Permeability of natural fracture system, m2 |
| : | Permeability of natural fracture system under initial condition, m2 |
| : | Permeability of matrix system, m2 |
| : | Permeability of vug system, m2 |
| : | Modified Bessel function of the second kind, zero order |
| : | Modified Bessel function of the second kind, first order |
| : | Modified Bessel function of the first kind, zero order |
| : | Reference length, m |
| : | Length of horizontal wellbore, m |
| : | Total number of hydraulic fractures |
| : | Number of discretized segments for each fracture |
| : | Pressure of natural fracture system, Pa |
| : | Pressure of matrix system, Pa |
| : | Pressure of vug system, Pa |
| : | Dimensionless pressure of natural fracture system, dimensionless |
| : | Dimensionless pressure of matrix system, dimensionless |
| : | Dimensionless pressure of vug system, dimensionless |
| : | Initial reservoir pressure, Pa |
| : | Reference pressure, Pa |
| : | Extra pressure drop caused by skin effect near the wellbore, Pa |
| : | Production rate of the line-sink under standard condition, m3/s |
| : | Dimensionless production rate of the line-sink under standard condition, dimensionless |
| : | Constant surface production rate of the multistage fractured horizontal well, m3/s |
| : | Flux density of the th segment in the th hydraulic fracture, m3/(s·m) |
| : | Laplace transformation of |
| : | Dimensionless flux density of the th segment in the th fracture, dimensionless |
| : | Laplace transformation of |
| : | The sandface flow rate of the multistage fractured horizontal well, m3/s |
| : | Mass flow rate between matrix and natural fractures per unit-volume reservoir, kg/(m3·s) |
| : | Mass flow rate between vug system and fracture system per unit-volume reservoir, kg/(m3·s) |
| : | Radius of horizontal wellbore, m |
| : | Radial distance, , m |
| : | Dimensionless radial distance, dimensionless |
| : | Skin factor, dimensionless |
| : | Time, s |
| : | Dimensionless time, dimensionless |
| : | Variable of Laplace transformation, dimensionless |
| : | Radial velocity for oil flow in natural fracture system, m/s |
| , : | - and -coordinates, m |
| , : | Dimensionless - and -coordinates, dimensionless |
| , : | - and -coordinates of the line-sink, m |
| Dimensionless - and -coordinates of the line-sink, dimensionless | |
| , : | - and -coordinates of the th discrete segment in the th fracture, m |
| : | Length of the left wing of the th fracture, m |
| : | Length of the right wing of the th fracture, m |
| : | Length of the discrete segment (), m |
| : | Dimensionless length of the discrete segment (), dimensionless |
| : | -coordinate of the intersection between the th fracture and -axis, m |
| : | Distance between and , |
| : | Shape factor of matrix blocks, m−2 |
| : | Shape factor of vug system, m−2 |
| : | Oil density in natural fracture system, kg/m3 |
| : | Oil density at the reference pressure , kg/m3 |
| : | Porosity of natural fracture system, fraction |
| : | Porosity of matrix system, fraction |
| : | Porosity of vug system, fraction |
| : | Porosity of natural fracture system at the reference pressure , fraction |
| : | Porosity of matrix system at the reference pressure , fraction |
| : | Porosity of vug system at the reference pressure , fraction |
| : | Storativity ratio of natural fracture system, dimensionless |
| : | Storativity ratio of matrix system, dimensionless |
| : | Storativity ratio of vug system, dimensionless |
| : | Interporosity flow coefficient between matrix and natural fractures, dimensionless |
| : | Interporosity flow coefficient between vugs and natural fractures, dimensionless |
| : | Permeability modulus, Pa−1 |
| : | Dimensionless permeability modulus, dimensionless |
| : | Oil viscosity, Pa·s |
| : | Circumference ratio, 3.1415926…, dimensionless. |
| D: | Dimensionless. |
| : | Laplace transformation. |
Conflict of Interests
The authors declare that there is no conflict of interests related to this paper.
Acknowledgments
Support for this work is provided by the National Natural Science Foundation of China (Grants nos. 51404206 and 51304165), the National Science Fund for Distinguished Young Scholars of China (Grant no. 51125019), the Project of National First-Level Discipline in Oil and Gas Engineering from China’s Central Government for Development of Local Colleges and Universities, and SWPU Science & Technology Fund (Project no. 2013XJZ005). The authors also would like to thank the reviewers and editors for their helpful comments.