# Expanded mixed multiscale finite element methods and their applications for flows in porous media

## Abstract

We develop a family of expanded mixed Multiscale Finite Element Methods (MsFEMs) and their hybridizations for second-order elliptic equations. This formulation expands the standard mixed Multiscale Finite Element formulation in the sense that four unknowns (hybrid formulation) are solved simultaneously: pressure, gradient of pressure, velocity, and Lagrange multipliers. We use multiscale basis functions for both the velocity and the gradient of pressure. In the expanded mixed MsFEM framework, we consider both separable and non-separable spatial scales. Specifically, we analyze the methods in three categories: periodic separable scales, convergent separable scales, and a continuum of scales. When there is no scale separation, using some global information can significantly improve the accuracy of the expanded mixed MsFEMs. We present a rigorous convergence analysis of these methods that includes both conforming and nonconforming formulations. Numerical results are presented for various multiscale models of flow in porous media with shale barriers that illustrate the efficacy of the proposed family of expanded mixed MsFEMs.

xpanded mixed multiscale finite element methods, non-separable scales, hybridization, two-phase flows

65N99, 65N30, 34E13

## 1 Introduction

There are many fundamental and practical problems involving a wide range of length scales. Typical examples include highly heterogeneous porous media and composite materials with fine micro-structures. In practice, these scales are too fine to treat with direct numerical simulation as the computational cost far exceeds the capabilities of computers for the foreseeable future. As a result, it is a significant challenge, both theoretically and numerically, to treat problems with multiple-scales effectively. A variety of numerical algorithms, varying from upscaling to multiscale methods, have been developed to capture the influence of fine-scale spatial heterogeneity on coarse-scale properties of the solution. In most upscaling methods the coarse-scale model is developed by numerically homogenizing (or averaging) the parameters of the fine-scale model, such as permeability. The form of the coarse-scale model is typically assumed to be the same as the fine-scale model. Simulations are performed using the coarse-scale model, and coarse-scale quantities of interest are readily computed. However, with this parameter upscaling, fine-scale properties of the solution can no longer be recovered. In contrast, multiscale methods carry fine-scale information throughout the simulation, and the coarse-scale equations are generally not expressed analytically, but rather formed and solved numerically.

Various numerical multiscale approaches have been developed during the past decade. A Multiscale Finite Element Method (MsFEM) was introduced in [23] and takes its origin from the pioneering work [9]. Its main idea is to incorporate fine-scale information into the finite element basis functions and capture their effect on coarse scales via finite element formulations. In many cases, the multiscale basis functions can be pre-computed and used repeatedly in subsequent computations with different source terms, boundary conditions and even slightly different coefficients. In some situations the bases can be updated adaptively. This leads to a large computational saving in upscaling multi-phase flows where the pressure equation needs to be solved many times dynamically. There are a number of other multiscale numerical methods with different approaches, such as residual free bubbles [11], the variational multiscale method [24], two-scale conservative subgrid upscaling [4] and the heterogeneous multiscale method [17]. Arbogast et al. [6] used a domain decomposition approach and variational mixed formulation to develop a multiscale mortar mixed MsFEM. Jenny et al. [25] have used a set of multiscale basis functions similar to [23] to develop a Multiscale Finite Volume Method (MsFV). A Multilevel Multiscale Method [28] was proposed in the framework of the Mimetic Finite Difference Method that efficiently evolves a hierarchy of coarse-scale models. In recent years, multiscale methods (e.g., [10, 15, 18]) have been developed to address high-contrast in multiscale coefficients.

The mixed multiscale method was first developed by Arbogast et.al. [4] as a locally conservative subgrid upscaling method, and was later analyzed and extended in [5] using a variational multiscale formulation. Chen and Hou developed a local multiscale basis equation for velocity and combined it with a mixed finite element formulation to propose a mixed MsFEM [14]. The mixed MsFEMs retain local conservation of mass and have been found to be useful for solving flow equations in heterogeneous porous media and other applications. However, in many practical situations, the permeability may be very low, and even vanish in some regions of the domain (e.g., permeability in shale). In this case its inverse is not readily available for use in the standard mixed MsFEM formulation, and hence, the direct application of these methods to practical problems may be restricted. To overcome this limitation, we propose a family of expanded mixed MsFEMs for these multiscale applications.

Expanded Mixed Finite Element Methods (MFEMs) were studied in classic finite element spaces in the past. Chen [13] developed and analyzed expanded MFEMs for second-elliptic equations, and obtained optimal error estimates for linear elliptic equations and certain nonlinear equations in standard finite element spaces. Woodward et al. [32] performed an error analysis of an expanded MFEM using the lowest-order Raviart-Thomas space for Richards equation. In [7] Arbogast et al. established the connection between the expanded MFEM and a certain cell-centered finite difference methods. A dual-dual formulation was introduced in [21] for the expanded MFEM with Raviart-Thomas spaces. The expanded formulation extends the standard mixed formulation in the sense that three variables are explicitly approximated, namely, the scalar unknown (e.g., pressure), its gradient and the velocity. The expanded MFEMs are suitable for the cases where the tensor coefficient of the underlying partial differential equations is small and even partially vanishing inside cells, which may be viewed as an extreme case of high-contrast of coefficients.

Porous media formations are often created by complex geological processes and may contain materials with a widely varying ability to transmit fluids. The permeability of the porous media may change dramatically, from almost impermeable barriers to highly permeable channels. This high-contrast of the permeability creates significant challenges for the simulation of flow through porous media. If the simulation fails to capture the influence of these barriers (e.g., shales) and channels (e.g., sand lenses), errors in the quantities of interest will likely be unacceptably large. Moreover, if the low permeability regions are thin with a small inter-barrier spacing, a fully resolved fine-grid discretization may be difficult to obtain and costly to apply. Thus, a multiscale simulation on a coarse grid may be necessary for these situations. To simulate the flows on a coarse grid when the permeability is almost vanishing in some regions, a reduced-contrast approximation was developed in MsFEM [16] to lower the variance of the coefficients. In [2], Aarnes et al. proposed an automatic approach to detect the barrier and iteratively split the coarse cells around barriers to obtain improved solutions in mixed MsFEM. For these barrier situations, the proposed expanded mixed MsFEM can provide an easy and direct approach to perform the computation.

The purpose of this work is to develop a framework of expanded mixed MsFEMs and rigorously analyze the convergence for different multiscale cases. The work enriches the studies on mixed MsFEMs from previous works [1, 2, 14]. The multiscale phenomena can be roughly classified into two categories: separable and non-separable spatial scales. In the case of separable scales one can localize the computation of multiscale basis functions. However, these approaches usually produce resonance errors Strategies exist to reduce the resonance errors, such as the oversampling technique introduced in [23]. Recently, a new technique was proposed in [22], where a zero-th order term is artificially added to the standard multiscale basis equation [23] to make the associated Green’s function decay exponentially, and consequently, the resonance error can be reduced. However, these strategies for reducing the resonance errors usually result in a nonconforming FEM, and moreover, they do not remedy the poor performance observed when the there is no scale separation. Instead, when media exhibit strong non-separable scales, some global information is needed for representing non-local effects. If the global information captures all relevant information from different scales, then resonance errors are removed and approximation accuracy is significantly improved [1, 10, 30]. If global information (or field) is used, we refer to the multiscale methods as global multiscale methods; otherwise, we refer to them as local multiscale methods. We consider the expanded mixed MsFEM for both separable and non-separable scales. Multiscale basis functions are employed for both the velocity and the gradient of the scalar. We present convergence analysis for three typical multiscale cases: periodic highly-oscillatory separable scales, convergent separable scales, and a continuum of scales. For the expanded mixed MsFEM, we consider both conforming multiscale bases and nonconforming/oversampling multiscale bases. The computation of expanded mixed MsFEM is very similar to the standard mixed MsFEM, but the expanded mixed MsFEM gives more information about the solution and provides an accurate approximation of the scalar field’s gradient as well. This gradient is often needed in applications, for example, the pressure gradient is used for flows with gravity in computing upstream directions.

We consider the proposed expanded MsFEMs in the non-hybrid (standard) form and the corresponding hybrid form. The hybridization of the mixed formulation is of increasing importance to the study of mixed methods. Hybridization was initially devised by Fraejis de Veubeke [19] as an efficient implementation technique for mixed finite elements. Specifically, hybridization localizes the mass matrix on each cell, and hence, enables the local elimination of the velocity to obtain a sparse positive definite system. Later hybridization was applied to produce a better approximation of the scalar unknown over each cell. For example, Arnold and Brezzi [8] showed that using the Lagrange multiplier unknowns introduced by hybridization in a local post-processing procedure improved the accuracy of the scalar unknown. In the framework of expanded mixed MsFEMs, we first show the equivalence between each expanded mixed MsFEM and its hybridization, and then analyze the convergence of these methods. In particular, we obtain the convergence rates of the Lagrange multipliers for each of the three multiscale cases.

The rest of the paper is organized as follows. Section 2 is devoted to formulating the expanded mixed formulation for a model elliptic equation in an abstract framework. In Section 3, we present multiscale basis functions and formulations of expanded mixed MsFEM for separable scales and continuum scales. In Section 4, we present convergence analysis for both expanded mixed MsFEMs and their hybridizations. Local multiscale methods and global multiscale methods are analyzed for the expanded mixed MsFEMs. In Section 5, the expanded mixed MsFEMs are applied to various multiscale models of flow in porous media to demonstrate their efficacy. Finally, some comments and conclusions are made.

## 2 Background

In this section we highlight the formulation of the expanded MFEM, and summarize important notation. We focus on the dual-dual formulation of the continuum problem because it most naturally supports our analysis. We review important results of the corresponding discretization as well.

### 2.1 The continuum expanded mixed formulation

Let be a bounded domain in , , with Lipschitz boundary . For a subdomain , and , and denote the usual Sobolev space and Legesgue space, respectively. The norm and seminorm of are denoted by and , respectively. When , is written as with norm and seminorm . The norm of is denoted by . We also use the following spaces

For a normed space on , denotes the underlying norm on . We denote to be the restriction of to subdomain . In the paper, is the usual inner product, and is the inner product.

We consider the following elliptic problem

(1) |

We may assume is symmetric positive definite. Let and . The equation (1) can be rewritten as

(2) |

Here usually refers to the fluid velocity. We define the following function spaces for solutions:

The expanded mixed formulation of (1) reads: find in , satisfying

(3) |

We have the following theorem, which addresses the relation between the solution of (1) and the solution of (3). {theorem} (See Theorem 3.5 in [13].) If is the solution of (3), then is the solution of (1), satisfying and . Conversely, if is the solution of (1), then (3) has the solution satisfying and .

For theoretical analysis, we introduce the operators associated with (3). Let and , where the notation denotes the dual of a space. These are defined, respectively, by

from which we define by

where is the transpose operator of . Let denote the null operator. Then can be rewritten as

Let and . They are defined, respectively, by

Using this operator notation, the expanded mixed formulation (3) is equivalent to finding in such that

(4) |

Equation (4) is called a dual-dual mixed formulation of (3) [20], since the operator itself has a similar dual-type structure. Because the operator satisfies a continuous inf-sup condition [12] and the operator satisfies a continuous inf-sup condition on , the system (4) has a unique solution in . Moreover, there exists a , independent of the solution, such that [21]

(5) |

where

### 2.2 The discrete expanded mixed formulation

Let , and be approximations of operators , and , respectively. Let be an approximation of . Then the numerical formulation of (4) reads

(6) |

Let be a numerical inner product. Let ,
and be finite dimensional approximation of ,
and , respectively. Set .
Then we have the following abstract result for well-posedness of
(6).
{lemma}[20]
Assume that

(1) there exists a positive constant independent of such
that for any ,

(7) |

(2) there exists a positive constant independent of such that for any

(8) |

(3) there exist positive constant and constant independent of such that

Then there exists a unique solution of (6).

We note that should be an element broken norm, provided that is a nonconforming space in , i.e, . We have the following Strang-type lemma for the problem (6). {lemma}[20] Let and be the unique solutions of (4) and (6), respectively. Then there exists , where , , and are defined in Lemma 2.2, such that

(9) |

Again, in Lemma 2.2 should be an element broken norm, provided that is a nonconforming space in . The second term on the right hand side of (9) is the so-called consistency error. It is easy to check that

Consequently, if , then

which is the case for a conforming expanded mixed FEM.

As a general remark, the generic constant is assumed to be independent of the mesh diameter throughout the paper.

## 3 Formulations of expanded mixed MsFEMs

In this section, we describe how to construct multiscale basis functions for expanded mixed MsFEMs.

Let be a quasi-uniform partition of and be a representative coarse mesh with . Let . Let be a classic mixed finite element space pair such as the Raviart-Thomas, Brezzi-Douglas-Marini or Brezzi-Douglas-Fortin-Marini spaces (cf. [12]). In this section, we present multiscale basis functions for standard local expanded mixed MsFEM, oversampling expanded mixed MsFEM, and global expanded mixed MsFEM.

### 3.1 Local expanded mixed MsFEM

Let (the finite element velocity space localized on ) and be a velocity basis function defined on . Following the mixed MsFEMs developed in [14], we define the multiscale basis equation for the local expanded mixed MsFEM in the following way:

(10) |

where is the outward unit normal to . We can obtain different mixed MsFEM basis functions corresponding to different choices of . Chen and Hou [14] choose to be the lowest Raviart-Thomas basis. The weak expanded mixed formulation of (10) is to find such that

(11) |

with on . The function is a basis function for the gradient variable and is a basis function for velocity . We define finite element spaces for the local expanded mixed MsFEM as follows:

The finite element space for pressure in the expanded mixed MsFEM is , which is the same as in the classic mixed finite element pairs. On each coarse element , the velocity basis functions are associated with its faces. Therefore, velocity basis functions of have the support of two coarse elements sharing a common interface. However, the basis functions of are supported in only one coarse element. Each basis function of is associated with a velocity basis function via , but it is restricted to only one coarse element K.

The local expanded mixed MsFEM formulation of the global problem (3) reads: find such that

(12) |

subject to the global boundary conditions.

It is well known [13] that the linear system of equations in the expanded mixed FEM produces an indefinite matrix, which is a considerable source of trouble when solving the linear system. We can introduce Lagrange multipliers on interfaces of cells and localize mass matrix to each element to obtain a sparse symmetric positive definite system by elimination, which is suitable for many linear solvers. This will give rise to a hybrid formulation. To this end we require some further notations.

For the hybrid formulation we define the MsFEM velocity space to be , whose normal components are not necessarily continuous on , the collection interior interfaces of . They are defined as follows:

We note that basis functions in do not require normal continuity. The basis function in is supported on a coarse grid, and this is different from the basis function in whose support is two coarse cells sharing the common face. We note that

The finite element space for the Lagrange multiplier is defined as

Define the operator by

(13) |

The hybridization of local expanded mixed MsFEM formulation of (12) reads: find such that

(14) |

### 3.2 Oversampling expanded mixed MsFEM

If the equation (10) is solved in a block larger than , and the interior information of the solution is taken to construct the basis functions in , then this results in the oversampling technique introduced in [23, 14]. We note that oversampling MsFEM is a modified local MsFEM. Using the oversampling can reduce resonance error.

Here we follow the outline in [14] to present the oversampling multiscale basis functions. Let be the solution of the equation

(15) |

Let and be indices of faces of . We define

(16) |

where the constants are chosen such that

(17) |

We define

We introduce an operator whose local form is defined by

(18) |

Then we define finite element spaces for the oversampling expanded mixed MsFEM as follows:

Here, is to impose some intrinsic continuity of the normal components of velocity multiscale basis functions across the interfaces . However, and so the oversampling multiscale method is nonconforming. The expanded oversampling mixed MsFEM formulation of the global problem (3) reads: find such that

(19) |

subject to the global boundary conditions.

### 3.3 Global expanded mixed MsFEM

To remove the resonance error and substantially improve the accuracy, we can use global information to construct the multiscale basis functions. Suppose that there exist global fields that capture the non-local features of the solution of equation (1). Following the global mixed MsFEM framework proposed in [1], we define the multiscale basis equations for the global expanded mixed MsFEM by

(21) |

where and are indices of the global information. References [1, 10, 30] provide some options for the global fields. To reduce the computational cost of using global fields, we can pre-compute these fields at some intermediate scale [26].

The finite element spaces for the global expanded mixed MsFEM are defined by

Consequently the global expanded mixed MsFEM formulation of the global problem (3) reads: find such that

(22) |

subject to the global boundary conditions.

We define the MsFEM velocity space and the finite element space of the Lagrange multipliers for the hybrid formulation of (22) by

The functions in may not have normal continuity. The relation between and is expressed by the identity .

Let the operator be defined in a similar way to (13). Then the hybridization of the global expanded mixed MsFEM formulation of (22) reads: find such that

(23) |

We note that (14), (20) and (23) can be rewritten as the operator formulation

(24) |

This is actually the hybrid approximation of the operator equation (6).

## 4 Convergence analysis

We will focus on the analysis of the local expanded mixed MsFEM. The analysis of the oversampling expanded mixed MsFEM and the global expanded mixed MsFEM is very similar to the local expanded mixed MsFEM but requires additional notation and definitions. We sketch the analysis of the oversampling expanded mixed MsFEM and global expanded mixed MsFEM and present their convergence results. For the local expanded mixed MsFEM, we consider micro-scales that can be treated in terms of periodic homogenization and convergent homogenization.

For the analysis, we will focus on the case when the velocity space on , , is the lowest-order Raviart-Thomas space () and . We use to build the source term and boundary conditions for the MsFEM basis defined in equation (10).

### 4.1 Inf-sup condition for expanded mixed MsFEMs

In this subsection, we discuss the inf-sup condition associated with problem (12). To simplify the presentation, we define three discrete operators associated with (12). Specifically, let , and , then the operators are defined as

By Lemma 2.2, we require that the operator is bounded and positive for well-posedness of problem (12). If satisfies the assumption

(25) |

then the operator is bounded and positive. However, we can weaken the assumption (25); for example, may be partially vanishing inside fine cells and positive and bounded everywhere else. Then the operator is still positive and bounded. For simplicity of presentation, we use assumption (25) for the analysis.

We have the following lemma for the discrete inf-sup condition of problem (12). {lemma} For any , there exists a positive constant independent of such that

(26) |

Let