Abstract
The homogenization of matrix and short fractures is one of the conventional approaches to deal with a plenty of fractures in different scales. However, the accuracy of this approach is still a question when long fractures and short fractures are distributed in the homogenized model. This paper describes a new hybrid method in which the long fractures will be modeled explicitly by the embedded fracture continuum approach and short fractures are considered through the homogenized technique. The author used this hybrid method to demonstrate the effect of fractures which are intersected to the well on the oil production rate as well as the elapsed time of a fractured reservoir in a depletion process. The advantages of the new hybrid method are easy assembling of numerous fractures into the model and incorporation of the complex fracture behaviour into the model.
1. Introduction
Prediction of the oil production rate of fractured reservoirs is still a challenge due to the complexity of fracture distributions as well as their behaviours. The permeability of a fracture is usually much larger than that of the rock matrix. Hence, the fluid will flow mostly through the fracture network if fractures are connected. This implies that the fluid transport in the naturally fractured reservoirs will be governed by the fracture connectivity and their distribution [1, 2]. There are three approaches commonly used in modeling naturally fractured formations: (1) continuum models using effective properties of representative element volume, (2) discrete fractured models in which fractures are introduced explicitly, and (3) hybrid models that combine the above two models.
The first approach to model a fractured reservoir in the literature is the homogenized approach [3–5]. The fact of using the equivalent medium in the modeling means that the porous space and permeability of the rock matrix as well as those of the fracture network are gathered and characterized through only a single porosity and a single permeability represented by the overall (effective) permeability. Thus, the case study (known hereafter as the homogenized model) considers only the effect of the single porosity and single permeability. The advantage and limitation of using the effective single continuum approach were discussed in [3, 4] by comparing with the traditional dual continuum model in [6]. These authors discussed that this approach can capture the baffled flow effect, while the dual continuum approach is better equipped to enhance flow through the fracture system which is however limited for the models they used, in the regular and orthogonal fracture geometry. The unsteady-state flow effects between fracture and matrix (matrix-fracture transfer) are not considered in the effective single continuum approach which however can be more easily addressed in the dual continuum model by using the matrix-fracture transfer functions (sugar cube-like matrix and fracture geometry). Thus, as detailed in [4], the effective single continuum modeling approach can be applied when the unsteady-state between matrix and surrounding fractures is not a dominant feature.
The second approach in this study is that fractures are explicitly introduced in the model. The question of what fracture should be considered explicitly in a model is an open question. In fact, in every modeling approach (and this is valuable for all models not only for fractures), the discontinuities of lower scale than that of observation are neglected. For example, when a rock sample is tested in the laboratory, the pores or cracks not detectable by naked observation are neglected. When cracks in site are observed, the cracks under a given value are simply neglected, and the media under this scale are considered as continuum. Likewise, in our model, only the long fractures with length superior to the representative element volume (REV) size will be considered explicitly. The choice of only long fractures modeled explicitly in this work comes from the fact that long factures play an important role in fractured reservoir behaviour [7]. This concept can be explained by the fact that, in reality, the fractured rock mass is highly heterogeneous, and the assumption of using the uniform aperture of the fracture networks to simplify the modeling by using the equivalent medium through the upscaling approach can over/underestimate the problem. In such heterogeneous media, there is a high probability of the existence of several initial fractures/faults (or even induced fractures during drilling) with a large aperture. The consideration of these fractures in the classical upscaling approach can violate the notion of REV. The scenario proposed in this work can be an alternative approach to simplify the problem and matches well with the idea proposed by some other scholars [2, 4]. In their work, Lee et al. [4] proposed a hierarchical approach for modeling fluid flow in a naturally fractured reservoir with multiple length-scale fractures. They classified the fractures as short, medium, and long fractures. The short and medium fractures were associated with the matrix through the homogenization technique to define the effective porous medium, while the long fractures (considered as the major conduits) are separated and modeled explicitly in the model. The flexibility and performance of this hierarchical approach were demonstrated by these authors in modeling the complex fractured rock masses.
In this paper, the embedded fracture continuum (EFC) approach [7–9] will be used to model a hypothetical reservoir whose natural fractures coincide with that of the Sellafield site [10]. The objective is to illustrate the effects of intersected fractures on the oil production rate of a fractured reservoir. The choice to use the data of Sellafield is motivated only by the completeness of the data. To demonstrate the effect of intersected fractures on the oil production rate of a fractured reservoir, two scenarios are studied: the first scenario is the homogenized approach as described above; the other is explicitly a fracture approach as explained in the preceding paragraph. The second scenario will be separated in two cases which aim to investigate in more detail the effect of these fractures on a hypothetical well behaviour (well production notably). In the first case, we consider that at least one of these long fractures intersects the well (known hereafter as the conductive fracture model), and in the other case, any long fractures intersect the well (known as the nonconductive fracture model).
The configuration of problems considered here is quite similar as the one studied in the work of [2]. In that work, the author would like to elucidate the influence of well and long fracture intersection on the productivity of the well. However, in comparison with that contribution [2], there are some essential differences in the present study. In effect, to consider the communication between porous matrix and fractures, as well as fractures and well, Li and Lee [2] used the transport index between matrix-fracture transfer and the fracture transmissibility to the well. The later parameter is determined by using directly Peaceman’s productivity index (PI) which is largely used in practice for the well drilled in the medium without fracture intersection. The author assumed that the pressure drop along the fracture inside the well block is negligible, and the productivities from the fracture and well are superposed. This assumption means that the fractures connected to the well will become part of extended well geometry. However, as the fracture surface is much larger than the well surface, the production from the fracture surface will be much larger than that from the well surface. Inversely, in the present work, we interest only the production calculated from the well surface. This well surface is explicitly considered in our model through which the production rate (fluid flux) will be calculated directly while the drop of pressure in the fractures is decided by the transient flow in their corresponding fracture cells.
The implementation of these fractures in the model at the present stage seems easy by using the EFC approach based on the fracture cell concept [7–9]. Note that each fracture cell in the model represents a porous medium which has its own porosity and permeability and is different from those of the fractured matrix. These models represent in effect the double porosity and double permeability medium to distinguish with the initial scenario based on the single effective medium. Otherwise, as mentioned above, in all calculations, we neglected the unsteady-state feature between matrix and fractures in this fractured matrix. This EFC approach was implemented with the help of an open source code deal.II [7–9, 11, 12].
2. Equivalent Elastic Properties and Permeability of the Fracture Cell in the Embedded Fracture Continuum Approach
2.1. Equivalent Elastic Properties of the Fracture Cell
The first work in this approach was reported by Figueiredo et al. [13, 14] in which they considered that the elastic properties of all fracture cells generated in the whole medium are isotropic and are characterized by the same Young’s modulus calculated from the fracture cell intersected by a horizontal fracture. Otherwise, Poisson’s ratio is equal to one of the intact rock:where Em, , kn, and h indicate, respectively, Young’s modulus, Poisson’s ratio of the intact matrix, the normal stiffness of fracture, and the element size.
More precisely, in these last works, one does not distinguish fracture cells intersected by one or many fractures (Figure 1(a)). It also means that this approximation is only applicable for fractures owing to the same mechanical properties.
(a)
(b)
(c)
In this work, another isotropic approximation of the equivalent poroelastic properties of the fracture cell is proposed. Difference with the previous approximation (equation (1)) is that the isotropic approximation introduced here allows modeling the fractures having different mechanical properties by accounting for the contribution of each fracture in the isotropic properties of the fracture cell. More precisely, for the fracture cell of type I intersected by an inclined fracture (Figure 1(b)), their approximated isotropic moduli are proposed as follows:where d is defined as follows:
Similarly, for the fracture cells of type II intersected by more than one fracture (Figure 1(c)), the effective properties are obtained using the hypothesis of noninteraction leading to the summation of all fracture contributions:where i and j denote ith and jth fractures.
2.2. Equivalent Permeability of the Fracture Cell
Based on the hypothesis that the flow in the fractured element is principally conducted through the fracture inside the domain and flow in fracture is assumed to be laminar between two infinite parallel smooth plates (cubic law), the permeability is defined as , where is the aperture of fracture. By equivalence of the flux flow through the fracture cell and one through the fracture, we have the formula to determine the isotropic equivalent permeability of the fracture cell of type I [7–9, 15]:
Correspondingly, for the fracture cell of type II, the contribution of multiple fractures on the isotropic permeability of the fracture cell is added up using the following general formula:
3. Application of the Embedded Fracture Continuum Approach to Investigate the Effect of Intersected Fractures on the Oil Production Rate of Fractured Reservoirs
The study consists of simulating the production of a vertical well (primary depletion) with 1 m of diameter. To this end, the 2D plane strain assumption is adopted. Otherwise, to simplify the problem, only a quarter of the model will be considered (Figure 2). Regarding the dimension of the well, the saturated reservoir of 200 m width and 200 m length (with respect to the center of the well) was selected. The initial pore pressure and total isotropic horizontal stress are assumed to be equal to 10.0 MPa and 24.46 MPa, respectively, which corresponds to the values measured at the depth of 725 m in the Sellafield site [7–9]. Otherwise, concerning the boundary conditions, we consider that the top and right boundaries are closed and fixed while a zero-constant pressure is imposed on the wall of the well.
(a)
(b)
(c)
As mentioned in Introduction, the simulation will be conducted with different scenarios. As the first scenario, we replace the fractured medium by the homogeneous equivalent poroelastic medium [7, 9]. Recall that, for Sellafield, the lower cutoff length for crack distribution is 0.5 m [10]. These fractures are embedded in the fractured matrix with a uniform aperture of 65 µm (Case 1). In the second scenario, the fractures are explicitly introduced in the model. The second scenario will be separated in two cases (Case 2 and Case 3) which aim to investigate in more detail the effect of these fractures on a hypothetical well behaviour (well production notably). As illustrated in Figure 2, no long fractures intersect the well in Case 2, and it is called the nonconductive fracture model, while in Case 3, at least one long fracture intersects the well, and it is hereafter called the conductive fracture model.
For these calculations, the following poroelastic properties of the equivalent porous matrix and long fractures are chosen:
For the equivalent fractured matrix, the effective properties calculated in [7, 9] will be used. As observed from these last calculations, the anisotropic degree of this fractured rock is low; thus, we can consider it as an isotropic material. In this case, the Biot modulus and the Biot coefficient can be determined from the following equations:where is the bulk modulus of skeleton of the rock matrix (Ks = 271 GPa) and is the porosity of the homogenized porous matrix (the initial porosity of the intact matrix is taken as ). The viscosity and compressibility of fluid (oil) are 0.01(Pa·s) and 5.0 × 10−10 (Pa−1), respectively. The effective properties of REV are summarized in Table 1.
For long fractures, as an example, only 6 fractures longer than 100 m with 1 mm of width are modeled explicitly (Figure 2). The properties of fractures like the normal stiffness and shear stiffness are taken to and , respectively, at the initial state. The corresponding permeability of these fractures based on the cubic law is Kf = 8.33 10−8 (m2).
In the following, for each model, the fluid flow simulation will be conducted. In this type of simulation, we consider that the fluid flows toward the well (drawdown effect) and hence well production (primary depletion) is from the contribution of both the permeable fractured matrix and long fractures. The mechanical effect will be ignored in the first subcase by considering the constant aperture and hydraulic properties of long fractures and will be then considered (but in an implicit manner) in the second subcase by varying the aperture and hydraulic properties of long fractures during the transient flow. More precisely, in the latter subcase, closure of fractures (updated aperture and permeability of fractures) will be taken into account as a consequence of the mechanics effect. From several contributions [16], it showed that geomechanics can impact the fluid flow notably through the closure/opening of fractures which in turn induces a variation of the fractures’ permeability. To account for this effect, many authors proposed to use an empirical stress-dependent function of the fractures’ aperture. Among different functions presented in the literature, the following relationship proposed in [17, 18] (called the Bandis and Barton model) is the most commonly used:
Following this model, the closure/opening of fractures is principally affected by the effective normal stress acting on their surfaces (shear stress is neglected). In equation (8), and are the aperture and normal stiffness of fracture at the initial state, respectively; is the actual fracture aperture calculated from the variation of the effective normal stress ( and indicate the actual and initial effective normal stress acting on the fracture plane, respectively). Otherwise, through the definition of the effective stress (we can consider the simplest case b = 1), if the variation of the total stress is negligible, the evolution of fracture aperture can be defined as a function of the pressure drawdown :
This last equation will be used in our flow simulations to activate the fracture aperture corresponding to the actual pressure () during the transient fluid flow. For the illustration purpose, this empirical relation is presented in Figure 3 by using the properties of the long fractures as mentioned above. Note that, in our code, this pressure drawdown is calculated as the mean value of pressure drawdown in each fracture cell; it means that the fractures’ aperture and permeability will be updated at the end of each time step for each segment of fracture embedded in the fracture cell.
The production rate and cumulative production of the well for the three considered models are presented in Figures 4 and 5 respectively. As the first observation, we can note that there is a significant difference for the third case when the long fractures intersected with the well keeping their initial aperture during the simulation. With respect to the other case, the production rate of this last case is much higher (about twenty times) in the first days of production and drops quickly after one month. This statement is well illustrated in Figure 5 with a linear variation of the cumulative production which attains finally its asymptotic value after this mentioned period. The difference in the production rate is however much smaller (about four times) when the aperture of the conductive fractures is updated during the simulation by considering the geomechanical effects described through the Bandis and Barton model. A rapid drop of pressure in the conductive fractures in the first days induced a closure of aperture and hence induced the decrease in the fractures’ permeability. Therefore, the production rate and the cumulative production decrease versus the elapsed time. The production is over at about 11 months (Table 2). The results obtained from the homogenized model or the nonconductive fracture model present no significant difference regardless of the aperture of fractures. The similar rate and cumulative production were stated in the whole simulation time with only moderate difference at the end days of the production process. The influence of the fracture closure as highlighted in Figures 4 and 5 is really negligible. The production in the nonconductive models is finished at about 17 months, while for the single effective medium scenario, it is about 19 months. For more information about the well behaviours in these studied scenarios, we presented in Figure 6 the distribution of pore pressure at 10 hours of fluid flow simulation. A drop of pressure in the fracture cells of the corresponding fracture intersected to the well is captured. Comparing with the case that the fracture aperture is updated, this drop is much more pronounced if the fracture aperture is kept constant during the simulation. It is important to note also that this drop is not instantaneous or uniform in the fractures but mostly concentrated around the well (about 30 m to 40 m with respect to the center of the well). It seems that the drop along the conductive fracture could play an important role on the productivity of well and cannot be neglected as assumed in the work in [2]. Besides, due to the rapid drawdown of pressure, the aperture of the conductive fractures decreases quickly which is about half of the initial aperture of 1 mm after 10 h and attains its maximal closure after 1 month (the corresponding final aperture is about 0.2 mm) (Figure 7). Two points (point A located at (38.3 m, 32.14 m) and point B located at (25.0 m, 43.30 m), as sketched in Figure 2), which are the same distance to the center of the wellbore, are selected to observe. The purpose of studying these two points is to compare the pore pressure behaviour between one point belonging to the long conductive fracture and the other belonging to the fractured matrix. All these observations are confirmed in Figures 8 and 9 in which the variation of pore pressure and fractures’ aperture and permeability at two points (point A and point B) with respect to elapsed time are highlighted. Quite similar results of pore pressure evolution were noted for all studied cases except for the conductive pressure case with constant aperture. The controlled point A which belongs to the long conductive fracture networks presents a drop of pressure at the very early instant (about 0.1 h). We also observed the drop of pressure at the controlled point B which belongs to the fractured matrix beginning later at about 100 h.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
The results obtained in this part confirm the conclusion in [2] which showed that the intersection of fractures and well plays a significant role on the well productivity. In this context, one observed a high production rate in the first stage which is followed by a quick drop in the second stage. This is also illustrated by the fact that the cumulative production will attain rapidly the maximum in comparison with the other cases of the equivalent single medium or nonintersected well and fractures. The geomechanics effect described in an implicit manner by using a pressure-dependent function will reduce this production rate and delay the primary depletion time. As mentioned in the first paragraph of this section, the two scenarios considered here represent in effect the single porosity/single permeability medium and double porosity/double permeability. Thus, thanks for these simulations [2]; we also highlighted the effect of these two concepts on the behaviour of the fractured reservoir.
4. Conclusions
The paper presented an effective hybrid simulation method to model oil production of naturally fractured reservoirs that include homogenized media of short fractures and explicated long fracture networks. A depletion simulation was calculated on the basis of this hybrid method. In this method, fractures are classified into “short” and “long” fractures on the basis of fracture length relative to the block size. The contributions from short fractures are calculated by a numerical method to effective properties of representative element volume of fractured reservoirs. Based on this new hybrid method, we studied the oil production rate as well as cumulative production of a hypothetical reservoir at the Sellafield site. A high production rate can be obtained early, and a sharp production drop follows to establish a stable rate if a fracture intersects with the well. The oil production rate is still much higher than the one in the case of a reservoir without intersection of the production well and fractures in the stability state.
In addition, the efficiency of the new hybrid approach on modeling a complex fracture network was demonstrated through this paper. This hybrid method reduces effectively the total number of explicated fractures in the model in comparison with the discrete fracture model. Thus, the computational time of this hybrid method is effective in the actual site. Moreover, the new hybrid method provides an easy way to develop nonlinear fractures as the Bandis and Barton model in this article.
The hybrid method of this paper does not work for the sugar cube approximation for matrix and fracture configurations as in the conventional dual porosity-permeability model. More precisely, the distribution of natural fractures can be directly modeled through homogenization of short fractures and networking of explicated long fractures. As a result, this hybrid method can be computationally much more efficient than the dual permeability and porosity model.
Data Availability
The numerical data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This paper is based on the results of PhD thesis obtained at the University of Orleans (France). The author would like to thank Prof. Dashnor HOXHA and Dr. Duc Phi DO of the University of Orleans (France) and his colleagues of the University of Transport and Communications (Vietnam) for their supports on this work. This work was supported by the Ministry of Education and Training of Socialist Republic of Vietnam (the No. 911 scholarship project).