Multiscale Simulations for Coupled Flow and Transport Using the Generalized Multiscale Finite Element Method

1 Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong, China; E-Mail: eric.t.chung@gmail.com 2 Center for Numerical Porous Media (NumPor), King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia 3 Department of Mathematics, Texas A&M University, College Station, TX 77843, USA; E-Mails: sidnet123@gmail.com (W.T.L.); ustcrenjun@gmail.com (J.R.)


Introduction
Many porous media problems occur over multiple scales.Because of scale disparity, some type of model reduction or averaging are often employed.Typical flow-based upscaling approaches ( [1][2][3]) compute the permeabilities on a coarse (computational) grid and solve the flow and transport equations together on a coarse grid.These approaches are easy to implement with a minimal code modification; however, these approaches can result in large errors for complex heterogeneities because limited information is used in upscaling.In this paper, we consider a multiscale method for a coupled flow and transport system, where the flow equation for the pressure field is described by a steady state elliptic equation (derived by Darcy law) and the transport equation for the concentration field is described by a convection-dominated parabolic equation.
Our approach derives its foundation from rigorous multiscale techniques developed recently .In our approach, multiscale basis functions are computed and used to solve flow equations.The main contributions of this paper are (1) to use multiscale basis functions for both flow and transport equations and (2) design of a novel mixed multiscale methods for convection-dominated transport equations within the context of Generalized Multiscale Finite Element Method (GMsFEM).The GMsFEM is introduced in [16], where the objective is to construct multiscale basis functions in each coarse domain.The main distinction of this approach from some of existing multiscale methods is that the GMsFEM systematically constructs multiscale basis functions using local spectral problems and the local snapshot spaces.Both of these concepts (local spectral problem and local snapshot spaces) are important in a design of the method.
In our proposed approach, we use a mixed formulation for both the flow and the transport equations.The use of a mixed formulation for the flow equation guarantees the mass conservation [7,[26][27][28][29][30].The use of a mixed formulation is important for the transport equation, guarantees the mass conservation and helps with the stabilization of it.In particular, by adding more flux basis functions, we can achieve a better stability based on our numerical studies [31].In the mixed formulation for a coupled system, we first define snapshot spaces for the appropriately defined fluxes in the flow and the concentration equations.For the pressure and the concentration fields, we use piecewise constant basis functions.The snapshot space represents the solution space in each coarse region and provides a mass conservative approach.In the snapshot space, we perform a local spectral decomposition to identify multiscale basis functions.The functions corresponding to dominant modes are used to construct coarse spaces for the fluxes.In particular, we use only a few degrees of freedom (maximum two) in each coarse region.For the test functions for the convection-dominated transport, we use the solutions of adjoint problem.As a result, we have a full coarse-grid approximation for the coupled flow and transport system.
We would like to discuss a relation between our proposed approaches and the flow-based upscaling [2,3].In the latter, the upscaled permeabilities are computed and, then, the coupled flow and transport equation is solved on a coarse grid.The computations of upscaled permeabilities can be compared to the computations of multiscale basis functions, which are performed locally in each coarse region by solving local problems.The resulting coarse-grid problem in the upscaled model is similar to our resulting coarse-grid system.There are several main advantages of the proposed method: (1) it can use more multiscale basis functions in each region and achieve a desired accuracy; (2) the number of multiscale basis functions can be adaptively selected in space and time; (3) our approach can easily reconstruct the fine-scale distributions of the solution; (4) our approaches can handle any coarse-grid geometries.
We consider several numerical examples by considering two different hydraulic conductivity fields.In our numerical examples, we use a relatively coarse grid with the coarsening factor of 20.We study the numerical errors in the concentration field and depict concentration profiles at different time instants.Our numerical results show that the use of only one basis function per coarse edge (which is similar to numerical upscaling of hydraulic conductivity) does not produce an accurate solution, while if we use two basis functions, we can obtain a much more accurate approximation of the solution.These preliminary numerical results show that our fully coarse-grid model for the coupled flow and transport equations can accurately predict the solution using only a few degrees of freedom in each coarse region.
Novelty.The novelty of our work can be summarized as follows: • The development of Generalized Multiscale Finite Element Method for a mixed formulation of the convection-diffusion equation using different test and trial spaces.
• The development and implementation of a mass conservative mixed GMsFEM for the coupled flow and transport equations.
The paper is organized as follows.In the next section, we present preliminaries and describe coarse and fine grids as well as a mixed formulation.Section 3 is devoted to the construction of multiscale basis functions.Numerical results are presented in Section 4. Finally, we present conclusions.

Preliminaries and Notations
We consider a coupled flow and transport system that arises in many hydrological applications Here Ω is the computational domain, T > 0 is a fixed time, κ is a hydraulic conductivity, c is the concentration, p is the pressure, D is a diffusivity of the medium, and f c , f q are source terms.The equations are subject to the boundary conditions c(x, t) = 0 on ∂Ω × [0, T ], p(x) = g on ∂Ω and the initial condition c(x, 0) = 0 in Ω.In order to have a mass conservation, we re-write the system in a mixed formulation.∂c ∂t ( Here q is an auxiliary variable, we call it flux or flux-concentration.It is well known that the Galerkin method inherits the stability of the continuous problem, and it yields to spurious oscillation when the convection coefficient is larger than the diffusive one (v D).The mixed finite element methods can be used to achieve a mass conservation.In this paper, we aim to construct multiscale basis functions for the flux q and the velocity v.As for the concentration c and the pressure p, we choose piecewise constant functions on a coarse grid.
Next, we introduce the discretized form of the Problem (2).We first comment on some notations of the coarse and the fine mesh.We assume that Ω is partitioned in a usual way into a union of rectangles (in 2-D) or hexadedrals (in 3-D) denoted by T H , where H is the coarse mesh size.The fine mesh is obtained by splitting each coarse block in T H into a connected union of fine-grid blocks, and we denote it by T h with mesh size h, h < H.In addition, we denote by E H the set of all edges/faces of elements in T H , i.e., E H = {E i : 1 ≤ i ≤ N e }, where N e is the number of coarse edges/faces.Note that each E i ∈ E H is the union of some fine-grid edges/faces.That is, E i = {e j : e j is on E i }, where e j are the fine-grid edges/faces lying on E i .We also define the coarse neighborhood of the edge/face E i by See Figure 1 for an illustration of a coarse neighborhood associated with one coarse edge, where the coarse mesh and the fine mesh are represented by black lines and gray lines, respectively.

A fine cell in
An illustration of a coarse neighborhood ω i (in green) corresponding to the coarse edge E i .
Next, we define the coarse trial spaces for q, v, c, and p.The concentration c and the pressure p are approximated in the space of piecewise constant functions with respect to the coarse grid T H .We denote this space by Q H , and As for the flux-concentration q and the velocity v, the basis functions are associated with coarse edges and are supported in the corresponding coarse neighborhoods.In particular, to obtain the basis functions for a coarse edge E i , we will solve a local problem in the coarse neighborhood ω i .Let Φ q i and Φ v i be the set of multiscale basis functions for q and v with respect to E i , respectively.Then, we define the multiscale solution space for q and v as Furthermore, we discretize the time interval [0, T ] uniformly by the points where τ is the time step and N = T /τ .Using the above spaces, the GMsFEM approximation of (2) can be given as follows.For all n ≥ 1, find subject to the boundary condition p H = g.Here, we let , where z = c H or q H .For θ = 0, 1/2 and 1, we get the Backward Euler, Crank-Nicolson and Forward Euler Methods respectively.Note that W q H is the testing space for q n H .We emphasize that in general, W q H = V q H , which allows a better stability.The construction of V q H , V v H , and W q H will be discussed in the next section.In addition, to get the fine-scale solution, we use the standard lowest-order Raviart-Thomas space for the approximation of (2) on the fine grid T h .The fine-scale solution (c h , q h ) will be considered as a reference solution for the comparison with the multiscale solution (c H , q H ).

Generalized Multiscale Finite Element Method
In this section, we will discuss the construction of the vector solution spaces V q H , V v H , and the testing space W q H .Note that solving q H in (6) requires v H . Therefore, we will first compute v H using the last two equations of (6).
We briefly introduce the basic idea of the GMsFEM.In the framework of GMsFEM, it follows an offline-online procedure to generate the multiscale solution space.At the offline stage, we first construct the snapshot space obtained by solutions of local problems with all possible boundary conditions, which are resolvable by the fine-grid.Then, we choose a subspace of the snapshot space by selecting dominant modes via some local spectral decompositions.The resulting reduced dimensional space, called the offline space, will be used at the online stage to get a multiscale solution.Note that the snapshot space can be taken as the solution space to obtain a solution with a good accuracy.However, it is expensive to use the snapshot space due to its high dimensionality.Thus, there is a need to perform a space reduction within the snapshot space, and choose only few dominant modes to form the offline space.We remark that for parameter dependent problems (in which the coefficients depend on some extra parameters), one needs to construct an online space using the offline space and an appropriate spectral problem.

Multiscale Solution Space V v H
We form the multiscale solution space V v H for v H .In order to build V v H , we first construct the snapshot space, denoted by V v snap .For each coarse-grid edge E i ∈ E H , a local problem is defined on each of its neighboring coarse cells.Because the snapshot space is supposed to capture the multiscale information within this local domain, the corresponding boundary conditions of the local problem should cover all possibilities.Formally, the local problem is defined as the following: For each E i ∈ E H , on its neighboring coarse cell K = K r or K s (see Figure 2), find Here, j varies over all fine-grid edges on E i .For the convenience, we use the same notations as the solution to denote local snapshot fields.Here α i j is a constant defined on K, n is a fixed unit-normal vector on ∂K, and δ i j (x) is defined on E i with respect to the fine-grid edges.Note that E i = {e i j : e i j in E i }.We require δ i j has value 1 on e i j and 0 on the other fine-grid edges, 1 ≤ j ≤ J i , where J i is the number of the fine-grid edges on E i .That is, 1, on e i j , 0, on the other fine-grid edges.
Neighboring coarse cells K r and K s corresponding to the coarse-grid edge E i .The fine-grid edge e i j and e i j are in red.
Note that the constant α i j in ( 7) is chosen so that the compatibility condition The above local Problem (7) can be solved numerically on the underlying fine grid of K by the lowest-order Raviart-Thomas element method.
After solving the local problems (on K r and K s ) with respect to E i , a local snapshot V v,i snap is generated by the solutions of (7).That is Then, the snapshot space V v snap is formed as We re-numerate the snapshot functions by a single index to create the snapshot matrix where J i is the dimension of the snapshot space V v snap .Note that the matrix R v snap maps from the coarse space V v snap to the fine space.The next step is to construct the offline space V v off for v H .For a coarse-grid edge E i , we define the following eigenvalue problem in V v,i snap : where and The eigenvalue Problem (9) will produce a set of eigenvectors for generating the local offline space V v,i off by first selecting L i off eigenvectors corresponding to the smallest L i off eigenvalues, and then Then, the local offline space with respect to E i is formed as and the offline space V v off is formed as We re-numerate the basis functions by a single index to create the offline matrix where We take the multiscale solution space as Then, v H is approximated by solving the problem: subject to the boundary condition p H = g.
Remark 1.The design of a "good" eigenvalue problem is a key ingredient of the GMsFEM.It plays a critical role in reducing the dimension of the coarse solution space.After applying such eigenvalue problem, a few dominant modes are selected to capture the multiscale behavior of the multiscale solution up to some certain extent.
Remark 2. Generally, the offline space should contain a constant basis function.Therefore, the local snapshot space V v,i snap , as a solution space of the eigenvalue Problem (9), can be decomposed into We note that the element with {(1, 1, • • • , 1)} in the snapshot space corresponds to the local solution with the unit normal flux on the edge and the elements of V v,i 2 have average zero normal flux traces.We eliminate the constant boundary condition basis function in the snapshot space and put it, (1, 1, 1 , as the first basis function of the offline space, and identify the other dominant modes in the space V v,i 2 .

Multiscale Solution Space V q H
We construct the solution space V q H and the testing space W q H for q H .They follow a similar procedure to the one for v H .In order to build V q H , the snapshot space V q snap is to be constructed.We define a local problem corresponding to each coarse-grid edge where K = K r or K s depicted in Figure 2. The local snapshot space V q,i snap , the snapshot space V q snap , and the snapshot matrix R q snap are formed as respectively.
To construct the testing space W q H , a snapshot space W q snap is also needed.We solve the adjoint problem of (12) w i q,j − D∇p i j = 0 in K, Then, the snapshot space W q snap is formed by the snapshot functions ω i q,j for all E i ∈ T H .
The local offline space V q,i off for u H is obtained by solving the following eigenvalue problem in V q,i snap : where Here, we let D = 1+|v H | D , where |v H | denotes the vector l 2 norm of v H .Then, the local offline space V q,i off , the offline space V q off , and the offline matrix R q off are formed as: respectively.Here φ i q,k = is the dominant mode in the space V q,i snap .
We let V q H = V q off and define the testing space as where n is the unit-normal vector of the coarse-grid edge E. Now, we have the multiscale solution space V q H and the testing space W q H .The GMsFEM solution (c H , q H ) can be computed through (6).
Remark 3. We remark that the fine-scale velocity field is constructed and used in computing multiscale basis functions for the concentration.In future, we plan to study joint multiscale basis construction procedures, which can compute multiscale basis functions for the concentration without solving the velocity field.Because the concentration field is a nonlinear function of the velocity field, an online procedure for computing multiscale basis functions for the concentration will be employed.

Numerical Results
We present a representative set of numerical experiments that demonstrate the performance of the mixed GMsFEM for approximating the coupled flow and transport Equation (2).The computational domain Ω = [0, 1] × [0, 1].We consider two different permeability fields κ, as depicted in Figures 3a,b.The source term f c to be used are plotted in Figures 3c,d, and f q is chosen to be 0. In all experiments below, the fine grid T h is fixed to be 200 × 200 uniform mesh, i.e., h = 1/200.The coarse grid consists of 20 × 20 uniform meshes, i.e., H = 1/20.As for the time discretization, the time step is chosen to be τ = 1 × 10 −4 and the Crank-Nicolson scheme is used.We denote c h as the average value of c h in each coarse cell, i.e., c h = 1 |K| K c h , for every K ∈ T H .The relative numerical errors are measured in L 2 norm for the concentration.In our numerical experiments, we approximate the velocity v H and the flux q H (at any time instant) in the offline space V v H and V q H in which each coarse edge is equipped with L i off basis functions.We choose at most L i off = 2 and compare concentration profiles with L i off = 1.Note that there are H/h basis functions associated with each coarse edge.We only select up to the two dominant modes per edge to form the offline space.
The first experiment is for the Problem (2), where the diffusivity D = 1, the source term f c = f 1 , the permeability field κ = κ 1 , and the boundary condition g = 7xy.The boundary condition represents the flow from lower left corner to upper right corner.The Peclet number in the entire domain is 180.18, while the coarse-grid Peclet numbers vary with the maximum one being 11.01.In Figure 4, we depict the averaged fine-scale solutions and the GMsFEM solutions at the time T = 0.015, 0.04, and 0.1.We choose the averaged fine-scale solution since piecewise constant basis functions are used in the approximation of the concentration on a coarse grid.The averaged fine-scale solution is depicted on the top of Figure 4.In the middle row, we depict the GMsFEM solution, where only one basis function per edge is used for the approximation of the velocity and one basis function is used for the approximation of the flux-concentration field.The L 2 errors in the concentration field are 21.22%, 30.16%, and 27.19% at T = 0.015, 0.04, and 0.1.At the bottom row of Figure 4, we depict the concentrations computed using two basis functions for the velocity and two basis functions for the flux-concentration field.As we observe, the results look much better when two basis functions are used.For example, at T = 0.1, we can observe that the concentration profile is similar to the averaged fine-scale concentration profile.
The L 2 errors in the concentration field are 8.1%, 11.1%, and 8.44% at T = 0.015, 0.04, and 0.1.Note that since we use 20 × 20 coarse grid and piecewise constant basis functions for the concentration, one expects a first order convergence with respect to the coarse mesh (i.e., about 5% error in our case).
In our second set of numerical examples, we consider the permeability κ = κ 2 , a source term f c = f 2 , and a boundary condition g = 8(x − y).We let the diffusivity to be D = 1.In this case, the problem becomes convection-dominated, where the coarse-grid Peclet number has a maximum value around 8.1.The numerical results are shown in Figure 5, where the top row shows the averaged fine-scale concentration profiles, the middle row shows the GMsFEM solution with one velocity basis per edge and one flux-concentration basis function, and the bottom row shows the GMsFEM solution, which uses two velocity basis functions and two flux-concentration basis functions.In all numerical results, we use piecewise constant pressure and piecewise constant concentration basis functions, as mentioned earlier.Thus, for the concentration field, there is an irreducible error of order H (5% in our case).First, we study the numerical results using only one multiscale basis function (the middle row of Figure 5).The L 2 errors in the concentration field are 77.6%,82.53%, and 77.96% at T = 0.01, 0.02, and 0.1.We observe from the figure that the results with one basis function per edge do not capture the averaged fine-scale solution behavior correctly.On the other hand, we observe a good agreement when two multiscale basis functions are used (the bottom row of Figure 5).We observe that the main features of the concentration profile are captured.The L 2 errors in the concentration field are 14.9%, 18.24%, and 15.71% at T = 0.01, 0.02, and 0.1.
In our numerical simulations, we limit ourselves to piecewise constant (lowest order) concentration basis functions.This allows a low computational cost and avoids a design of more complex concentration basis functions, which can satisfy inf-sup conditions (see e.g., [32]).Since we use lowest order basis functions for the concentration, the accuracy is limited to O(H).This is reasonable as in large-scale simulations (due to very detailed κ), the coarse grid sizes can be taken sufficiently small.In general, the numerical results shown above can be improved using more basis functions for the flux-concentration and designing multiscale basis functions for the concentration field.One of our main objectives is to show that the results with one basis function per edge (which is similar to flow-based upscaling) can be substantially improved if one uses more multiscale basis functions, and we also provide a procedure for computing the multiscale basis functions and suitable global formulations.

Discussions on Time Dependent Velocity
In our examples above, the velocity field was time-independent, and thus, the coupling between the flow and transport is weak.In above examples, we derived a multiscale coupling framework and demonstrated a robust framework for multiscale convection-diffusion equation.When the velocity field becomes time-dependent, one needs to re-calculate multiscale basis functions as velocity varies.There are several procedures, which we will mention, after a brief discussion on time varying velocity field scenarios.First, we remark on several scenarios, when the velocity field can depend on time (we call it "time-dependent" in our further discussions, even though the time may enter implicitly).The velocity field can depend on time via time-varying concentration field as in miscible or immiscible flows (e.g., [33] and the referenes therein) or the velocity field can, in addition, have a time-dependent component due to compressibility.We plan to address these application problems in our future works.
Computations of multiscale basis functions for the time-dependent velocity field can be performed offline [16,34] or using the residuals in the online stage [12,13].For the offline computations, one needs a richer class of snapshot vectors, which can span possible velocity fields and perform the same local spectral decompositions and global coupling framework, as we discussed above."Rich" snapshot spaces are discussed in [16] for parameter-dependent problem (the time can be treated as a parameter).Another approach for identifying multiscale basis functions, when the snapshot space is very large, is the use of l 1 -minimization techniques as proposed in [34].All these approaches can be used adaptively in space and time.Finally, one can use residual-based online basis functions [12,13], where one additional multiscale basis function is computed adaptively using the residual information.The success of these approaches depends on the offline space as it is shown, thus a good offline space is needed to achieve a high accuracy.
In some time-dependent velocity cases, one can get a good accuracy without modifying multiscale basis functions.For simplicity, we consider a "miscible" type flow and transport, which are described by We take ζ(c) = (c 2 + 0.2(1 − c) 2 ).We use the same setup as in the first example.The source term is the same, that is, f c = f 1 and the boundary condition is g = 28xy.The numerical results are shown in Figure 6.The first row shows the averaged fine-scale concentration profiles and the second row shows the GMsFEM solution with two velocity basis and two flux-concentration basis functions.The L 2 errors in the concentration field are 5.81%, 8.63%, and 7.72% at T = 0.015, 0.04, and 0.1, respectively.These errors are comparable to the case with time-independent velocity field.We note that the are close to the irreducible error due to the coarse-mesh size.

Conclusions
In this paper, we develop a Generalized Multiscale Finite Element Method (GMsFEM) for coupled flow and transport equations.The transport equation is convection dominated and we choose an appropriate test space to achieve a stability and improve numerical results.We study a mixed formulation for both flow and transport, which guarantees mass conservation.The multiscale spaces for the flux and the velocity fields are constructed by appropriately choosing the snapshot spaces and performing local spectral decompositions.The design of the local test and trial multiscale spaces is one of our novel contributions.From the numerical experiments, we observe that one needs to have more than one multiscale basis function to capture the concentration profile on the coarse grid.In particular, using only two multiscale basis functions for the pressure equation and two multiscale basis functions for the concentration equation, we can achieve good accuracy.Our future studies will include a systematic design of multiscale basis functions for the concentration that can provide a stable discretization.

f 2 Figure 3 .
Figure 3. Two permeability fields (a,b) and two source terms (c,d).

Figure 4 .
Figure 4. First row: Fine-scale solutions c h .Second row: GMsFEM solutions c H (one multiscale basis function per edge for both velocity and flux-concentration).Third row: GMsFEM solutions c H (two multiscale basis function per edge for both velocity and flux-concentration).From left column to right column: T = 0.015, 0.04, 0.1.We use κ = κ 1 , f c = f 1 , H = 1/20.

Figure 5 .
Figure 5. First row: Fine-scale solutions c h .Second row: GMsFEM solutions c H (one multiscale basis function per edge for both velocity and flux-concentration).Third row: GMsFEM solutions c H (two multiscale basis function per edge for both velocity and flux-concentration).From left column to right column: T = 0.01, 0.02, 0.1.We use κ = κ 2 , f c = f 2 , H = 1/20.