Generalized Multiscale Multicontinuum Model for Fractured Vuggy Carbonate Reservoirs

Generalized Multiscale Multicontinuum Model for Fractured Vuggy Carbonate Reservoirs

Min Wang Department of Mathematics, Texas A&M University, College Station, TX 77843, USA (    Siu Wun Cheung Department of Mathematics, Texas A&M University, College Station, TX 77843, USA (    Eric T. Chung Department of Mathematics, The Chinese University of Hong Kong, Shatin, New Territories, Hong Kong SAR, China (    Maria Vasilyeva Institute for Scientific Computation, Texas A&M University, College Station, TX 77843 & Department of Computational Technologies, North-Eastern Federal University, Yakutsk, Republic of Sakha (Yakutia), Russia, 677980. (    Yuhe Wang Department of Petroleum Engineering, Texas A&M University at Qatar, Doha, Qatar (

Simulating flow in a highly heterogeneous reservoir with multiscale characteristics could be considerably demanding. To tackle this problem, we propose a numerical scheme coupling the Generalized Multiscale Finite Element Method (GMsFEM) with a triple-continuum model aimed at a faster simulator framework that can explicitly represent the interactions among different continua. To further enrich the descriptive ability of our proposed model, we combine the Discrete Fracture Model (DFM) to model the local effects of discrete fractures. In the proposed model, GMsFEM, as an advanced model reduction technique, enables capturing the multiscale flow dynamics. This is accomplished by systematically generating an approximation space through solving a series of local snapshot and spectral problems. The resulting eigen-functions can pass the local features to the global level when acting as basis functions in coarse problems. Our goal in this paper is to further improve the accuracy of flow simulation in complicated reservoirs especially for the case when multiple discrete fractures located in single coarse neighborhood and multiscale finite element methods fail. Together with a detailed description of the model, several numerical experiments are conducted to confirm the success of our proposed method. A rigid proof is also given in the aspect of numerical analysis.

1 Introduction

Simulating fluid flow in a fractured vuggy carbonate reservoir has always been of great interest and challenge to both academia and the petroleum industry. A large share of worldwide hydrocarbon resources are stored in such reservoirs. For this reason, it is a long-lasting quest of many scientists and engineers to accurately model the underground flow dynamics for a better understanding of fractured vuggy carbonates. Nevertheless, the simulation of flow within such heterogeneous formations is notorious for natural coexistence of different continua, and distinct scales of underlying porous media. Not only the sizes of vugs are in different orders of magnitude, the fracture configurations are also highly diverse. On top of all, the connectivity of vugs and fractures are complicated, further exacerbating the fluid flow modeling in such media. Mass would transfer among different scales and continua in different forms, which makes flow modeling a task full of challenges.

Among many tools that have been developed to address the characteristics of such problems, an easy and convenient way to explore the dynamics is to use a fine scale simulation. A fine partition is required to divide the domain into local pieces and then a global description of the flow is obtained by puttting the local solutions together. A common fine scale simulations can be conducted under frameworks such as the Finite Element Method (FEM) [4], Finite Volume Method (FVM) or Finite Difference Method (FDM).

Yet, due to high complexity of the media, it is impossible by nature for a fine simulation to simultaneously resolve all-scale media effects. One would either lose important small-scale information or run into a system that is extremely expensive to solve. Thus, a model-reduction scheme is necessary to design a practical numerical method for flow simulation in complicated porous media such as fractured vuggy carbonates.

Homogenization is commonly used in many model reduction methods [18, 3]. With homogenization, the domain of interest is partitioned into many coarse blocks, and the effective properties are calculated for each coarse block aiming at homogenizing local heterogeneous media using information in finer scale within the block. These pre-computed properties can thus catch and average fine scale characteristics and further calibrate the coarse solution accordingly [8]. This up-scaling scheme has been proved to be quite effective for simulations on media with scale separation or periodicity [13], but it fails to model the interaction between matrix, fractures and vugs.

That is exactly what motivated the multi-continuum model. To represent the flow between different continua, each continua in the domain of interest is now considered as a system that expands the whole domain. For example, a fractured and vuggy heterogeneous carbonate reservoir can now be described as three parallel continuum: fractures, matrix and vugs. Different continua coexists at every spot of the domain while they macroscopically interact with each other. The interaction between different continua is coupled based on the mass conservation law. For each continua, both intra and inter flow transfers are modeled. The first multi-continuum model, the dual porosity model (DPM), is proposed by Barenblatt for flow through naturally fissured rock [5]. In his work, two continuum were proposed to represent low and high porosity continua, respectively [23, 21]. And later, a third continua was introduced by researchers to generalize the application of such an idea [20, 27, 25, 24]. Thanks to the simplicity and flexibility of this method, it has been widely adopted in many fields [17]. However, the limitation of the multi-continuum model is also prominent, as it assumes all continuum are connected globally. This assumption will only be valid when each continua has merely global effects. For well-developed fractures, one continua is sufficient to represent its effects, yet such continua fails to model, for example, the independent long fractures as theses fractures may have local responses that can not be globalized.

Such ”dicontinuum,” like discrete fractures, are pretty common in reservoir modeling . This inspires us to also consider a discrete model. Ideally, the discrete models should be able to describe medias such as fractures that mainly contribute to local flow transfers. One common choice is the Embedded Fracture Model (EDFM). Discrete fractures and rock volume cells are considered as two separate systems, and the flow transfer between them is modeled element by element [15, 16, 6]. EDFM is computationally effective in the sense its discretization of rock (matrix) is independent of the spatial distribution of fractures. Yet, its effectiveness is sensitive to conductivity contrast. Another classic “dicontinuum” approach is the Discrete Fracture Model (DFM). DFM, on the other hand, can also be used to describe cross-flow between the fractures and matrix which are spatially adjacent to each other [19, 14, 1, 2]. Unlike developed fractures, each discrete fracture is described using a dimension element instead of a -dimensional continua. Researches have shown that using discrete networks can accurately resolve the local characteristics of the media with discrete fractures regardless of the conductivity contrast. Discontinuous models for vugs were also developed which take both free and porous flow into consideration [26].

In order to overcome the constraints of the homogenization scheme, and enrich the heterogeneous information reserved from the local fine-sacle region, the Multiscale Finite Element Method (MsFEM) was developed [10]. Like homogenization, the domain of interest is first partitioned with a coarse mesh. Each local block is then further partitioned with a finer mesh. With this two-level mesh, MsFEM can build a series of basis functions in each local region. These basis functions are obtained by solving a series of local problems which incorporate the heterogeneous characteristics of the local coarse region. The basis functions for all coarse regions will work together as a new global approximation space to replace the standard one used in FEM. With such setting, the size of the resulting numerical system is reduced significantly. MsFEM can be easily coupled with DFM when modeling flow in a fractured reservoir. However, more in-depth investigation indicates that when there are more than one independent fracture network in a coarse local region, MsFEM is not able to restore the local dynamics accurately [27].

The Generalized Multiscale Finite Element Method (GMsFEM) was later developed to tackle such a problem. It has been proven that GMsFEM can strengthen the ability of MsFEM on solving multi-scale problems. By conducting a spectral decomposition over the local snapshot space, GMsFEM can identify basis functions corresponding to dominant modes of local heterogeneous regions [9, 10]. This makes automatic enrichment of the multiscale space possible [11, 22]. In the paper [7], a one-one correspondence between GMsFEM basis functions and high-conductivity networks was presented. This property of basis makes GMsFEM a necessity when dealing with practical examples like carbonate reservoir simulations, as many fracture channels may coexist in a single local region. By nature of MsFEM and GMsFEM, the production of basis functions are independently conducted within each coarse neighborhood. And we can employ parallel computing to speed up the multiscale basis generation process.

With a GMsFEM framework, a multicontinuum model, and the discrete fracture network, we can inherit the merits of the three, and couple them together to improve the capability as well as accuracy of our simulation. In this paper, we use this fully coupled system to describe the flow in fractured and vuggy heterogeneous reservoirs. Conforming unstructured mesh is used to surrender the random discrete fracture networks. Fractures are treated hierarchically. Highly developed fractures with only global effects are modeled as a fracture continuum, while fractures that have local effects are embedded as discrete fracture networks. For independent vugs, a continuum is used to represent their effects with specific configurations such that no intra-flow is considered. GMsFEM allows us to consistently develop an approximation space that contains prominent sub-grid scale heterogeneous background information based on the multi-continuum and DFM.

The paper is organized as follows: in Section 2, the problem under discussion is clarified, followed by Section 3 which briefly reviews the multi-continuum model. In Section 4, a step-by-step illustration on GMsFEM together with a priori error estimate is provided. The details of time discretization of our problem is also discussed in this section. In Section 5, we present multiple numerical results to varify the effectiveness of porposed methods. Lastly, this paper is concluded by Section 6.

2 Preliminaries

In this paper, we consider a 2-dimensional flow problem in a multiscale porous media. We assume that the dynamic of flow is governed by the Darcy’s Law. For simplicity, we ignore the gravity and the capillary pressure effects.

2.1 Equation for slightly compressive flow in porous media

In specific, we consider the following equation for slightly compressive flow in heterogeneous porous media,


Here, is the computational domain. is compressibility and is viscosity of the liquid. is the formation volume factor (FVF) at reference pressure and is a FVF at reservoir condition. They are used to quantify compressibility of the target liquid. represents porosity of the fracture vuggy media, while is a permeability function that bears multiscale features in the media (See Figure 1 for an illustration). The solution to be sought is pressure , given a production rate .

Figure 1: Permeability field with multiscale features

Limiting our interests to slightly compressible liquid, we can further employ the simplified correlation between the formation volume factor and the pressure


to rewrite (1) and get


In the following sections, we will derive our method based on (3) along with Dirichlet or Neumann boundary conditions on

Throughout this paper, we assume , and are constants. (2) can then be reformulated as



is a constant, while

is a field map in .

2.2 Fine-scale spatial discretization

For flow in a fractured and vuggy media, the multiscale flow problem described in (4) becomes more complicated as the fractures and vugs have very different hydraulic properties from its background matrix. They can bring in extra transfer and storage mechanics to the flow. The fractures amplify the complexity of modeling as they can have a wide range of scales and topology. In order to delicately model the their effects on flow, we apply a hierarchical approach. Fractures that have only local effects can be resolved by fine mesh. Thus, the computational domain can be partitioned into


where and correspond to matrix and fracture regions respectively. is a 1-dimensional region that represents one resolved fracture with an aperture . Those fractures that have global effects and are not resolved by mesh can later be handled by representing them as one continua. So are the effects of vugs, which will be discussed in details in next section. For a fine-scale approximation of , we discretize the PDE on a fine grid, and apply Finite Element Method as well as DFM. Specifically, all integrations in the weak form of (4), will now be taken separately in both and with distinct hydraulic parameters. To compromise arbitrary fractures , one need to adopt an unstructured fine-scale mesh. The resulting semi-discrete numerical system is


Here, is a standard FEM basis function. means taking directional derivative along the degenerated fracture . Note that the aperture effects are considered in .

are again constants.

3 Multi-continuum Model

To explicitly represent the global effects of unresolved fractures, vugs and matrix, we introduce the multi-continuum methods. We consider the media as a coupled system of three parallel continua: matrix, unresolved fractures(usually natural fractures), and vugs. They coexist everywhere in our computational domain, while they interact with each other via mass transfer (see Figure 2 for an illustration).

Figure 2: Illustration of triple-continuum model

Without of loss of generality, we assume that all continuum interact with each other. If we denote the flow pressure for continua as , and write the interaction between continua and as , we can then establish a system of PDE following (4) to describe the flow mechanism in each continua.



Here can be , ,or which stands for matrix, unresolved fractures and vugs respectively. We further assume that there is no intra-flow inside the vugs and all vugs only act as a storage in this system. That is to say, we only consider the case when all vugs are independent from each other. Mass transfer due to inflow of liquid along vugs can be disregarded in any element of the domain. Therefore, mass balance equation for vugs can be written as


The term represents the mass transfer from continua to continua . This transfer can be modeled using [23, 21]

where . Here is a shape factor, and is taken as harmonic mean of the permeability and .

With (7) and (8), we can derive the weak formulation of our proposed triple-continuum system. For matrix and unresolved fracture, we have


For vugs, we have


Here, is any testing function in the same space as . We mention that equation of , and are coupled through term , thus this coupled system should be solved on a Cartesian product space . In our proposed approach, we take for all continuum .

To express effects of both unresolved and resolved fractures on flow dynamics, we manage to incorporate DFM when solving this multicontinuum equation system (9)–(10). Like what we have in (6), we assume corresponds to a 1-D domain that serves as a resolved fracture region. All integrations on is thus rewritten as . For example, (9) can be rewritten as


after applying DFM to its original form. Similarly, incorporating DFM in (10) yields


The fine-scale FEM solution should be sought in , where , where the is a conforming finite element space of the continuum on a fine partition of domain. We also remark that the shape factor is taken to be proportional to .

For the purpose of simpler notations in the analysis presented in Appendix A, we rewrite the derived system (11)– (12) in a more general -continuum setting. First, we denote the Sobolev space . On each continuum, given a fixed , we define bilinear forms:


Given a fixed , we further define the following coupled bilinear forms on


Then the weak formulation (11)– (12) can be written as: find , where , such that for all , where ,


with continua and the continuum indices representing the matrix, fracture and vug components in order.


In order to reduce the computational cost, we would like to solve the equation system (7) and (8) on coarse mesh. However, permeability coefficient is heterogeneous in space, thus a standard FEM solution on coarse mesh will be inaccurate as it loses subgrid information. Therefore, we use GMsFEM to construct multiscale basis that contains local heterogeneous permeability information. By replacing the standard FEM basis with GMsFEM basis, we are able to obtain a better accuracy and sustain an affordable computational cost.

In this section, we briefly review the procedure for GMsFEM. Roughly speaking, the construction of GMsFEM basis consists of two stages: solving snapshot problems and conducting spectral decomposition. Both steps are conducted locally.

4.1 Coarse and Fine Mesh

We first introduce the coarse grid with mesh size . And each coarse block in can be denoted as . can be further refined by an unstructured fine mesh with mesh size . See Figure 3 for an illustration. We assume is fine enough to resolve all underlying fine-scale properties of . Let be the set of all coarse nodes of , where is the total number of coarse nodes. We then define the coarse neighborhood of node as

Figure 3: Coarse & Fine Mesh. Left: coarse mesh with discrete fractures; Upper Right: A coarse neighborhood; Lower Right: A coarse block with a discrete fracture and refined mesh.

4.2 Snapshot Space

A snapshot space is an auxiliary space constructed within each coarse neighborhood . We omit the subscript for simplicity. There are a few different ways of constructing snapshot space [9]. In this paper, we take solutions to the following harmonic extension problems as snapshot basis functions of three coninuum. The snapshot space is exactly the span of all such basis functions.

The snapshot problems are designed analogue to the steady state equation of (7) and (8). We consider a coupled snapshot system in a coarse neighborhood , in which we find such that


is defined on all fine-scale nodes of . If the set represents all fine-scale nodes on boundary, we have

Here, is standard basis in . So far, we have constructed the local snapshot space as:

The global snapshot space is defined as the sum of all local snapshot spaces, i.e.

Remark When solving local snapshot problem (16) on the fine mesh within , one should also apply the idea of DFM and replace all integral by and all coefficient correspondingly.

4.3 Spectral Problem

To further reduce the dimension of resulting system, we conduct a spectral decomposition on . Such decomposition will automatically detect the dominant modes. More precisely, we sought eigenpairs for the following local spectral problem



The form of and are inspired by analysis which will be demonstrated in next section along with Appendix A. We sort the eigenvalues of (17) in ascending order, and we take the first eigenfunctions . Then the -th multiscale basis function in is defined by

where is a partition of unity function for coarse grid on a coarse neighborhood . By multiplying , we obtained a set of conforming multiscale basis functions supported in . Using the multiscale basis functions for all coarse regions , we construct the multiscale space

We remark that , where is the standard FEM approximation space on . When the multiscale space is established, we can find a coarse-scale solution on with less computational effort.

Once the multscale space is constructed, the GMsFEM solution is given by: find , where , such that for all , where ,


4.4 A-priori error estimates

In this section, we present some a-priori error estimates of the semi-discrete problem. The proofs of these estimates will be left to Appendix A.

We suppose the field has a upper bound and a lower bound on . We further assume that the fields and has a uniform upper bound and a uniform lower bound , i.e.


Next, we introduce some metrics on . The bilinear form can further induce a norm

We also define a norm by




The first theorem provides an estimate of the error between the weak solution and the multiscale solution by the projection error of onto the multiscale space in various metrics.

Theorem 1.

Let be the weak solution in (15) and be the multiscale numerical solution in (18). Assume and . Then we have


In light of Theorem 1, we have to establish an estimate of the projection error of onto the multiscale space in various metrics on the right hand side of (22), in order to complete the convergence analysis. With the assumption that the irreducible error between the Sobolev space and the snapshot space is small, which holds when a sufficiently large number of snapshot solutions is taken, we define an approximation of in the snapshot space by


and provide an estimate of the projection error of onto the snapshot space .

Theorem 2.

Let and be reference solution and snapshot projection of as defined in (15) and (23). Then we have



4.5 An implementation view

In this section, we derive the fully discrete system and present the implementation details. We adopt the implicit Euler scheme for time discretization to the semi-discrete GMsFEM system (18). Suppose the time domain is partitioned into equal subintervals of length , and denote the -th time instant by . Using backward difference, the fully discrete GMsFEM scheme is to, successively for find such that


where the subscript denotes the evaluation of a time-dependent function at the time instant and an initial condition is given. At each time instant , (25) gives rise to a nonlinear algebraic system in the coefficients with respect to the multiscale basis functions. With a sufficiently small time step size, we can adopt a direct linearization approach by replacing the field by and derive


Alternatively, we can use an iterative approach. More precisely, we can construct a sequence whose fixed point is the solution and truncate the successive iterations when a stopping criterion is fulfilled. In this case, we start with an initial guess and solve for


We remark that it is equivalent to the linearization approach if we stop after one iteration.

5 Numerical Results

In this section, we apply our proposed methods to a realistic fractured and vuggy reservoir. All three continuum have heterogeneous permeability background (see Figure 1 for the permeability of matrix) and discrete fracture networks are embeded in this reservoir like in Figure 4. An unstructured fine mesh is used to resolve the discrete fractures networks(see Figure 5). The descriptive parameters of this reservoir are listed in Figure 1 and Table 4. All numerical results are implemented using FEniCS Library.

Figure 4: Idealized discrete fracture network(DFN)
Figure 5: Unstructured fine mesh

Numerical experiments are conducted from different aspects. Performance are compared between MsFEM and GMsFEM, nonzero source term and nonzero mixed boundary condition. We also discuss the impact of the number of basis function selected to the solution accuracy. We remark that all examples are conduced using direct linearization approch as the iterative approach do not significantly improve the results for our problem, which indicates that the nonlinearity in our problem is not very strong.

5.1 Comparison of MsFEM and GMsFEM

In this subsection, we discuss the necessity to apply GMsFEM. From Figure 6, we can tell that, even with similar number of degrees of freedom, the MsFEM is not able to resolve the true solution, thus GMsFEM must be applied to generate meaningful results. This is especially true when there are multiple discrete fracture networks coexist in a single coarse neighborhood. Many numerical experiments have shown that MsFEM basis functions are not able to handle homogeneous background and multiple discrete fracture networks simultaneously. Figure 6 shows the solution we obtained using MsFEM and GMsFEM respectively when a single source is placed at the bottom left corner. The error of MsFEM solution can be as large as .

Figure 6: Comparison between GMsFEM and MsFEM solution with heterogeneous background and discrete fracture network. Left: GMsFEM solution with DOF=2646. Right MsFEM solution with DOF =2400.

5.2 GMsFEM solution for different boundary condition and source

In this subsection, we demonstrate the performance of our proposed triple continuum GMsFEM solution to problem (7) and (8) ,where lagging coefficient scheme is used to linearize the problem.

Different boundary condition and source term settings are tested for coupled GMsFEM approach.

Figure 7: Triple-Continuum, heterogeneous background flow simulation matrix with top and bottom nonzero Dirichlet boundary condition. Zero Neumann boundary is applied to left and right boundary. First row: fine-scale reference solution, DOF = 80229. Second row: Coupled coarse-scale GMsFEM solution with 8 basis, DOF = 3528.
Number of Basis Day 1 Day 10 Day 20
2 17.21 27.22 66.44
4 14.88 17.27 43.65
8 4.72 11.86 13.31
16 4.24 12.05 12.58
Table 1: relative errors(%) of numerical results for mixed boundary condition. Nonzero Dirichlet boundary condition is imposed on top and bottom boundary. Zero Neumann boundary is applied to left and right boundary.
Figure 8: Illustration of error trend with time for different number of basis for dirichlet boundary condition case
Figure 9: Flow simulation results for a triple continuum heterogeneous background matrix with no flow boundary condition. Injector located at bottom left corner. First row : Fine-scale reference solution. Second row: Coupled coarse-scale GMsFEM solution. DOF is same as in Figure 7.
Number of Basis Day 1 Day 10 Day 20
2 15.79 10.05 11.42
4 5.48 5.89 8.53
8 2.84 6.20 8.51
16 1.12 6.30 8.49
Table 2: relative errors(%) of numerical results for zero Neumann boundary condition.
Figure 10: Illustration of error trend with time for different number of basis for single source

From both error tables and solution figures , we come to the conclusion that: 1) For nonzero mixed boundary condition case, the GMsFEM solution can obtain a good result when using 8 basis or more. 2) For zero Neumann boundary and single point source term case, the coupled approach can obtain good approximation of fine-scale solution with 4 basis or more. 3) For both cases, the coupled approach can give us an acceptable solution.

From Figure 10 , Figure 8, Table 2 and Table 1, we can tell that the error of solution decrease when we increase the number of eigen-functions used.

6 Conclusion

In this paper, we proposed a triple continuum GMsFEM method as a fast solver of flow problems in heterogeneous domain. A fractured and vuggy reservoir is modeled as a coupled system of three continuum. Fractures are treated hierarchically, fractures with only global effects are considered as a continua, while the ones have local effects are represented as discrete fracture networks using DFM. The system coupling DFM and three continuum are discretized spatially following the Generalized Multiscale Finite Element Method(GMsFEM) for accurate and fast solution. Coupled assembling is provided to construct GMsFEM multiscale space. The convergence of our proposed method is proved strictly following mild assumptions. Later, the performance is tested using multiple examples with different settings. We conclude that GMsFEM is necessary for complicated discrete fracture networks, the proposed approach can provide competitive approximation for both mixed boundary conditions and a single source case. The number of basis are also discussed and chosen. From the numerical exmaples, we can see that selecting enough number of basis is crutial to the accuracy of our proposed method.

In short, we claim that our proposed method can accomplish the flow simulation task with both accuracy and efficiency. Nevertheless, we notice that our proposed method is only good for the case when a clear description of the discrete fracture networks is known. For reservoirs containing uncertainties, further exploration is desired. Besides, for vugs with turbulent flow inside, one will end up with a coupled PDE system containing Navier-Stokes system. Future investigations are required to expand our work to such cases.

Appendix A Proofs of error estimates

In this section, we present the proofs of the error estimates in Theorem 1 and Theorem 2.

a.1 Proof of Theorem 1


Using (15) and (18), we have

Let and take , we have

From this equation ,we can further get the following by the definition of and the bounded condition (19) of ,

By Cauchy-Schwarz Inequality, this implies


The last term on the right-hand side of (28) can be written as


Following [12], we employ generalized Holder’s Inequality and the definition of to obtain


Further, with Ladyzhenskaya’s Inequality, there exists some constant such that


There also exist some constant such that

For the fracture part, we have


To sum up, we have for any ,


for some constant . Plug back to (28), and notice that we can use Young’s Inequality to derive