A Generalized Multiscale Finite Element Method for Poroelasticity Problems I: Linear Problems ^{†}^{†}thanks: This work was supported by RFBR (project N 153120856)
Abstract
In this paper, we consider the numerical solution of poroelasticity problems that are of Biot type and develop a general algorithm for solving coupled systems. We discuss the challenges associated with mechanics and flow problems in heterogeneous media. The two primary issues being the multiscale nature of the media and the solutions of the fluid and mechanics variables traditionally developed with separate grids and methods. For the numerical solution we develop and implement a Generalized Multiscale Finite Element Method (GMsFEM) that solves problem on a coarse grid by constructing local multiscale basis functions. The procedure begins with construction of multiscale bases for both displacement and pressure in each coarse block. Using a snapshot space and local spectral problems, we construct a basis of reduced dimension. Finally, after multiplying by a multiscale partitions of unity, the multiscale basis is constructed in the offline phase and the coarse grid problem then can be solved for arbitrary forcing and boundary conditions. We implement this algorithm on two heterogenous media and compute error between the multiscale solution with the finescale solutions. Randomized oversampling and forcing strategies are also tested.
1 Introduction
Problems of mechanics and flow in porous media have wide ranging applications in many areas of science and engineering. Particularly in geomechanical modeling and its applications to reservoir engineering for enhanced production and environmental safety due to overburden subsidence and compaction [16, 17]. One of the key challenges is the multiscale nature of the geomechanical problems. Heterogeneity of reservoir properties should be accurately accounted in the geomechanical model, and this requires a high resolution solve that adds many degrees of freedom that can be computationally costly. Moroever, there are disparate scales between the often relatively thin reservoir structure and the large overburden surrounding the reservoir that adds more complexity to the simulation. Therefore, we propose a multiscale method to attempt overcome some of these challenges.
The basic mathematical structure of the poroelasticity models are usually coupled equations for pressure and displacements known as Biot type models [21]. For pressure, or flow equations, we have the parabolic equation Darcy equation with a time dependent coupling to volumetric strain. The stress equation is the quasistatic elasticity equations with a coupling to the pressure gradients as a forcing. Poroelastic models of this type have been explored in the petroleum engineering literature in the context of geomechanics for some time [18, 19, 9, 12] to name just a few. There are noted issues that arise. The first being heterogeneities of the reservoir and surrounding media add many complications to the effective simulation due to complexity of scales. Moreover, development of flow and mechanics simulation were often considered separately. Progress was made on this problem by considering various coupling strategies [19]. However, in the instance that the physics is not well understood a fully coupled scheme may be desired. This separation of development from flow and mechanics methods adds the complication of the computational grids not being the same in each regime. Some effort has been made in the improvements of gridding techniques between geomechanical and flow calculations [20] and references therein.
As briefly noted before, typically for numerical solution of such coupled systems time splitting schemes are often used. Various splitting techniques for poroelastic equations have been explored and analyzed in the context of reservoir geomechanics in [10, 11]. Also, in the context of poroelasticty and thermoelastic equations, various splitting techniques have been analyzed and implemented [15]. The primary splitting techniques are the undrained, fixedstress, and fully implicit. Due to observed better errors, we will primarily consider the less computationally costly fixedstress splitting and the more robust, yet with a loss in some matrix sparsity, fully implicit coupled approach.
Once the equations have been split in time we wish to resolve in space and will utilize a multiscale method. There are many very effective multiscale frameworks that have been developed in recent years. There are rigorous approaches based on homogenization of partial differential equations, where effective equations are derived based finescale equations at the microstructure level [8, 7]. However, these approaches may have limited computational use and more practical multiscale methods are used. Examples include the Heterogeneous Multiscale Method (HMM), where macroscale equations on coarsegrids are solved while the effective coefficients on the finescale are resolved at each coarse grid nodes [22, 23]. An approach based on the Variational Multiscale Method (see [24]), where coarsegrid quasiinterpolation operators are used to build an orthogonal splitting into a multiscale space and a finescale space [25]. Finescale space corrections are then localized to create a computationally efficient scheme. In this paper, we will use the Generalized Multiscale Finite Element Method framework, which is a generalization of the multiscale finite element method [2].
To efficiently solve these splitting schemes and overcome some of the challenges of heterogenous reservoir properties and gridding issues between mechanics and flow, we will develop a Generalized Multiscale Finite Element Method (GMsFEM) [1]. Our GMsFEM has the advantage of being able to capture small scale features from the heterogeneities into coarsegrid basis functions and offline spaces, as well as having a unified computational grids for both mechanics and flow solves. The offline multiscale basis construction may proceed in both fluid and mechanics in parallel and both constructions are comparable. We proceed by first generating a coarsegrid and in each grid block a local static problem with varying boundary conditions is solved to construct the snapshot spaces. We then perform a dimension reduction of the snapshot space by solving auxiliary eigenvalue problems. Taking the corresponding smallest eigenpairs, and multiplying by a multiscale partition of unity we are able to construct our offline basis. In this greatly reduced dimension offline basis, the online solutions may be calculated for pressure and displacements for any viable boundary condition or forcing.
The work is organized as follows. In Section 2 we provide the mathematical background of the poroelasticity problem. We will introduce the Biot type model and highlight where the heterogeneities primarily occur. In our formulation, the computational domain will be entirely inside of the fluid filled, or reservoir, region. However, coupling to regimes of pure elasticity to model the overburden are of course possible. In Section 3, to outline the difficulties in full direct numerical simulation we introduce the finescale discretizations using coupled and splitted schemes. Once we split the porooelastic system we will be able to apply our multiscale method. In Section 4, we present our GMsFEM algorithm and outline its construction procedure. We will use the offline multiscale basis functions to calculate accurately pressure and displacements, at a reduced dimension and computational cost in the online phase. Finally, numerical implementations are presented in Section 5. Using the GMsFEM, we compare the multiscale solution to finescale solutions and give error estimates. We will present two different examples with varying coefficients. Additionally, we will implement and discuss different strategies with oversampling and randomized forcings to construct the multiscale spaces.
2 Problem formulation
We denote our computational domain to be a bounded Lipschitz region. We consider linear poroelasticity problem where we wish to find a pressure and displacements satisfying
(1a)  
(1b) 
with initial condition for pressure We write the boundary of the domain into four sections . We suppose the following boundary conditions on each portion
and
Here the primary sources of the heterogeneities in the physical properties arise from , the stress tensor and , the permeability. We denote to be the Biot modulus, is the fluid viscosity, and is the BiotWillis fluidsolid coupling coefficient. Here, is a source term representing injection or production processes and is the unit normal to the boundary. Body forces, such as gravity, are neglected. In the case of a linear elastic stressstrain constitutive relation we have that the stress tensor and symmetric strain gradient may be expressed as
where , are Lame coefficients, is the identity tensor. In the case where the media has heterogeneous material properties the coefficients and may be highly variable.
The above poroelasticity problem (1a), assuming a linear elastic stressstrain relation, can be written in operator matrix form:
(2)  
(3) 
where
and and are gradient and divergence operators and .
3 FineScale Discretization
We will now present splitting methods for the above system in the context of solving the finescale approximation. This will highlight the areas where we would like to utilize a multiscale method when solving in the spatial variables due to the degrees of freedom required in resolving the system. For approximating the numerical solution to (1) on finescale grid we use a standard finite element method. We begin by giving the corresponding variational form of the continuous problem written as
(4)  
(5) 
for , where
and the test spaces with homogeneous boundary conditions are given by
Here for bilinear and linear forms we have define
and
Here under the integrand denotes the standard inner product. In Section 5, we will discretize the spaces using a finescale standard FEM and denote them and , being the finegrid size. The FEM using these spaces will serve as a reference solution for our GMsFEM outlined in Section 4.
To solve the above system we first discretize in time. This discretization leads to several possible couplings between timesteps and the two equations of prorelasticity. We proceed by giving the coupled and socalled fixedstress splitting [10, 15]. The standard fully implicit finitedifference scheme, or coupled scheme, can be used for the timediscretization and is given by
(6a)  
(6b) 
with , , where , , and . For time discretization we can apply many different splitting techniques which often occur in the literature.
Another we shall consider here is the fixedstress splitting scheme
(7a)  
(7b) 
where the variational form is rewritten with
and is the drained modulus
where is the Poisson ratio and is the elastic modulus. When we utlize the fixedstress splitting scheme, first we solve pressure equation for given data at the previous timesteps. Then, passing this new pressure information, we return to the quasistatic stress equation and calculate displacements at .
4 GMsFEM for Poroelasticity
In the GMsFEM presented here, we will focus on the development in the fixedstress splitting (7). We will however give numerical examples from both coupling strategies. The fixedstress splitting decouples the flow and mechanics equations. We will first present the offline multiscale basis construction in the fluid or pressure solve then its construction in the mechanics or displacement calculation step. In this algorithm, due to the heterogeneities arising primarily from the permeability and the stress tensor , we will solve local problems in each of the relevant portions of the variational form to construct the offline multiscale spaces. We now outline the general procedure of the GMsFEM algorithm.
The overall finescale model equations will be solved on a finegrid using spaces and , and will act as our reference solutions. Once the finegrid is established we must introduce the concepts of coarsegrids and their relationships. To this end, let be a standard conforming partition of the computational domain into finite elements. We refer to this partition as the coarsegrid and assume that each coarse element is partitioned into a connected union of fine grid blocks. The fine grid partition will be denoted by , and is by definition a refinement of the coarse grid . We use , where is the number of coarse nodes, to denote the vertices of the coarse mesh , and define the neighborhood of the node by
See Figure 1 for an illustration of neighborhoods and elements subordinated to the coarse discretization. We emphasize that the use of is to denote a coarse neighborhood, and we use to denote a coarse element throughout the paper.
Boadly speaking, the GMsFEM algorithm consist of several steps:

Step 1: Generate the coarsegrid, .

Step 2: Construct the snapshot space, used to compute an offline space, by solving many local problems on the finegrid.

Step 3: Construct a small dimensional offline space by performing dimension reduction in the space of local snapshots.

Step 4: Use small dimensional offline space to find the solution of a coarsegrid problem for any force term and/or boundary condition.
As noted previously, because coupled system of equations for poroelasticity problems can be solved using splitting scheme, we can construct multisclate basis functions for pressure and displacements separately. We begin by considering the pressure solve, then, the displacement solve.
4.1 Pressure Solve
Recall, for the numerical solution of pressure equation on coarse grid we consider the continuous Galerkin (CG) formulation (7b) given by
(8) 
where is used to denote the space spanned by multiscale basis functions , each of which is supported in . The index represents the numbering of these multiscale basis functions. We will now show how to construct the offline multiscale space . In turn, the CG solution of the form
will be sought.
We begin by construction of a snapshot space . We use harmonic extensions
(9) 
Here are defined by , where denotes the finegrid boundary node on . For simplicity, we will omit the index when there is no ambiguity.
Let be the number of functions in the snapshot space in the region , and define
for each coarse subdomain . We denote the corresponding matrix of snapshot functions to be
To construct the offline space , we perform a dimension reduction of the space of snapshots by using an auxiliary spectral decomposition. More precisely, we solve the eigenvalue problem in the space of snapshots:
(10) 
where
where and denote fine scale matrices
Here, are finescale basis functions.
We then choose the smallest eigenvalues from Eq. (10) and form the corresponding eigenvectors in the space of snapshots by setting
for , where are the coordinates of the vector . We denote the span of this reduced space as .
For construction of the offline space, to ensure the functions we construct form an conforming basis, we define multiscale partition of unity functions
(11)  
for all . Here is a continuous on and is linear on each edge of . We could choose to also be selected shape function, Neumann conditions, or boundary conditions on larger domains in the context of oversampling.
Finally, we multiply the partition of unity functions by the eigenfunctions in the offline space to construct the resulting basis functions
(12) 
where denotes the number of offline eigenvectors that are chosen for each coarse node . We note that the construction in Eq. (12) yields continuous basis functions due to the multiplication of offline eigenvectors with the initial (continuous) partition of unity. Next, we define the continuous Galerkin spectral multiscale space as
(13) 
Using a single index notation, we may write , where denotes the total number of basis functions in the spaces , for .
Denote the matrix
where are used to denote the nodal values of each basis function defined on the fine grid. Then, the variational form in (8) yields the following linear algebraic system
(14) 
where
Here, being the operator corresponding right hand side data from the previous time step and denotes the coarsescale nodal values of the discrete CG solution. We also note that the operator matrix may be analogously used in order to project coarse scale solutions onto the fine grid
4.2 Displacement Solve
We now suppose that we have solved for the finegrid pressure by the GMsFEM pressure solve in the previous section. We must now solve the mechanics equations (7a). Since the construction of the multiscale offline space remains very similar in this setting, we will be a bit more brief on its construction. Recall, for discretization of the displacements equation we rewrite equation as follows
(15) 
where . The corresponding continuous Galerkin (CG) formulation for displacements equations is given by:
(16) 
where , where are finescale basis functions, and we construct the multiscale offline space .
For construction of multiscale basis functions for displacements we use similar algorithm that we used for pressure. For construction of a snapshot space we solve following problem in
(17) 
Let be the number of functions in the snapshot space in the region , and define
for each coarse subdomain . Note we are using the same notation but with different harmonic extensions. We denote the corresponding matrix of snapshot functions, again with similar notation, to be
Again, we perform a dimension reduction of the space of snapshots by using an auxiliary spectral decomposition. We solve the eigenvalue problem in the space of snapshots
(18) 
where where
where and denote fine scale matrices
and
Here, are finescale basis functions.
We then choose the smallest eigenvalues from Eq. (18) and form the corresponding eigenvectors in the space of snapshots by setting
for , where are the coordinates of the vector . We denote the span of this reduced space as .
For construciton of multiscale partition of unity functions for the mechanics solve, we proceed as before and solve for all
(19)  
(20) 
Here is a continuous function on and is linear on each edge of . Finally, we multiply the partition of unity functions by the eigenfunctions in the offline space to construct the resulting basis functions
(21) 
where denotes the number of offline eigenvectors that are chosen for each coarse node . Next, we define the spectral multiscale space as
(22) 
Using a single index notation, we may write , where denotes the total number of basis functions in the space , for all .
And after construction we denote the matrix
where are used to denote the nodal values of each basis function defined on the fine grid. Then, the variational form in (16) yields the following linear algebraic system
(23) 
where , and
5 Numerical Examples
In this section, we present numerical examples to demonstrate the performance of the GMsFEM for computing the solution of the poroelasticity problem in heterogenous domains. Although we presented the algorithm in the fixedstress splitting, we are able to apply the same offline spaces as their construction remains the same in the fully coupled setting. However, in the coupled setting the equations (14) and (23) will no longer be decoupled and must be solved simultaneously.
We will implement a single complicated geometry with contrasting parameter values. We provide two cases one with lower contrast in elastic properties and another with higher contrast. We present the algorithm applied to these heterogenous coefficients in both the fully coupled and fixed stress time splittings. We give the errors with varying multiscale basis functions and over time. We then will apply the GMsFEM method with oversampling and with snapshots with randomized boundary conditions to obtain good accuracy, while having to solve fewer snapshot solutions. The effects of higher contrast in properties will also be discussed.
5.1 GMsFEM Implementation
First, we take the computational domain as a unit square , and set the source term in (1). We utilize heterogeneous coefficients that have different values in two subdomains. We denote each region as subdomain 1 and 2, and use following coefficients: for the Biot modulus we take and for permeability in the two regions. For fluid viscosity we take and fluidsolid coupling constant . For the elastic properties, we present results for two test cases. In Case 1, the elastic modulus is given by in each respective subdomain and in Case 2, we have . The Poisson’s ratio is , and these can be related to the parameters and , for via the relation
in each subdomain. The subdomains for coefficients shown in Fig. 2, where the background media in red is the subdomain 1, and isolated particles and strips in blue are the subdomain 2.
As we have chosen we must use boundary conditions to force flow and mechanics. In these tests, we use following boundary conditions:
and
and finally,
Here and are left and right boundaries, and are top and bottom boundaries respectively. We set and to drive the flow, and thus, the mechanics.
In Fig. 3 we show the coarse and fine grids. The coarse grid consists of 36 nodes and 50 triangle cells, and the fine mesh consists of 3721 nodes and 7200 triangle cells. The number of time steps is and the maximal time being set at . As an initial condition for pressure we use . The reference solution computed by using a standard FEM (linear basis functions for pressure and displacements) on the fine grid and using a fully coupled scheme. The pressure and the displacement fields for Case on the finegrid are presented on the left column of Fig. 4  5.
We test the fully coupled and fixedstress splitting schemes. The errors will be measured in weighted and weighted norm and seminorm for pressure
and for displacements
Here and are finescale and coarsescale using GMsFEM solutions, respectively for pressure and displacements.
Recall, we will use a few multiscale basis functions per each coarse node , and these number of coarse basis defines the problem size (dimension of offline spaces, and ). We suppose that in each patch we take the same number of multiscale basis functions for pressure, , for . Similarly for displacements we take , for . Varying the basis functions in both pressure and displacement multiscale spaces we recorded the errors at the final times.
In Tables 1 and 2, we present the weighted and errors for Case 1 and Case 2 of the coefficients in geometry Fig. 2 using the fully coupled scheme. We compare these to a finescale solution space with dimension 11163. In these tables, and are number of multiscale basis functions for each neighborhoods, the second column show the dimension of the offline space, the next two columns present the weighted and errors for pressure and last two columns show the weighted and errors for displacements. We see that the errors in pressure remain similar in both cases because the permeability parameters remain the same and the change is in elastic properties between scenarios. In Case 2, pictured in Table 2, we see great errors in displacements throughout when compared to Case 1 in Table 1 because the elastic properties in Case 2 have several orders of higher contrast.
In a similar setting, we consider the fixedstress splitting. For Case 1 we present the results in Table 3, the errors are very similar compared to the corresponding fully coupled scheme. This may be because we are comparing a finescale fully coupled scheme to a multiscale fully coupled scheme and similarly, a finescale splitting scheme to a multiscale splitting scheme and the errors do not differ very much between the two schemes here. For Case 2 we present the errors in Table 4 and again see that the errors are higher when compared to the lower contrast scenario. Comparing these results with the Case 2 using the fully coupled scheme, presented in Table 2, we see that both the pressure errors and displacement errors are much greater in this sequential coupling. This disparity is particularly striking when few multiscale basis functions are used.
We also include plots over time of the error with respect to number of basis functions used. We present the results from the fully coupled scheme. In Fig. 6 and 7 we show errors over time for and multiscale basis functions for each Thus, the dimensions of offline spaces are 432, 864, 1296 and 1728, respectively. We observe that errors decrease as we increase the dimension of the offline space as expected. We observe the errors in Fig. 6 are generally better than the errors Fig. 7, again, due to the lower contrast in Case 1. We see that in both cases most of the error vanished after the use of just 8 multiscale basis functions. In general, the error remains stable in time with a slight decrease over time.
Pressure errors  Displacements errors  
dim()  
2  216  0.06  0.08  0.06  0.13 
2  360  0.06  0.08  0.06  0.12 
4  432  0.01  0.01  0.04  0.11 
2  648  0.06  0.08  0.02  0.06 
4  720  0.01  0.01  0.01  0.03 
8  864  0.0003  0.002  0.002  0.03 
2  936  0.06  0.08  0.02  0.05 
4  1008  0.01  0.01  0.01  0.02 
8  1152  0.0003  0.002  0.0009  0.01 
12  1296  0.0001  0.001  0.0009  0.01 
2  1224  0.06  0.08  0.02  0.05 
4  1296  0.01  0.01  0.01  0.01 
8  1440  0.0003  0.002  0.0008  0.01 
12  1584  0.0001  0.001  0.0007  0.01 
16  1728  0.0001  0.0007  0.0007  0.01 
Pressure errors  Displacements errors  
dim()  
2  216  0.06  0.08  0.25  0.26 
2  360  0.06  0.08  0.22  0.24 
4  432  0.02  0.01  0.19  0.24 
2  648  0.06  0.08  0.08  0.13 
4  720  0.02  0.01  0.01  0.08 
8  864  0.001  0.002  0.01  0.08 
2  936  0.06  0.08  0.07  0.11 
4  1008  0.02  0.01  0.02  0.04 
8  1152  0.0003  0.002  0.004  0.03 
12  1296  0.0001  0.001  0.004  0.03 
2  1224  0.06  0.08  0.07  0.11 
4  1296  0.02  0.01  0.02  0.03 
8  1440  0.0003  0.002  0.001  0.02 
12  1584  0.0001  0.001  0.001  0.02 
16  1728  0.0001  0.0006  0.001  0.02 
Pressure errors  Displacements errors  
dim()  
2  216  0.06  0.08  0.06  0.13 
2  360  0.06  0.08  0.06  0.12 
4  432  0.01  0.01  0.04  0.11 
2  648  0.06  0.08  0.02  0.06 
4  720  0.01  0.01  0.01  0.03 
8  864  0.0003  0.002  0.002  0.03 
2  936  0.06  0.08  0.02  0.05 
4  1008  0.01  0.01  0.01  0.02 
8  1152  0.0003  0.002  0.0009  0.01 
12  1296  0.0001  0.001  0.0009  0.01 
2  1224  0.06  0.08  0.02  0.05 
4  1296  0.01  0.01  0.01  0.01 
8  1440  0.0003  0.002  0.0008  0.01 
12  1584  0.0001  0.001  0.0007  0.01 
16  1728  0.0001  0.0007  0.0007  0.01 
Pressure errors  Displacements errors  
dim()  
2  216  0.30  0.26  0.45  0.46 
2  360  0.30  0.26  0.42  0.45 
4  432  0.01  0.01  0.33  0.38 
2  648  0.30  0.25  0.36  0.48 
4  720  0.006  0.01  0.04  0.15 
8  864  0.001  0.006  0.04  0.15 
2  936  0.30  0.25  0.37  0.50 
4  1008  0.006  0.01  0.007  0.06 
8  1152  0.002  0.006  0.007  0.06 
12  1296  0.001  0.004  0.007  0.06 
2  1224  0.30  0.25  0.38  0.50 
4  1296  0.006  0.01  0.003  0.03 
8  1440  0.002  0.006  0.002  0.02 
12  1584  0.001  0.004  0.002  0.02 
16  1728  0.0009  0.003  0.002  0.02 
5.2 GMsFEM with Randomized Oversampling
In this section we consider the oversampling randomized algorithm proposed in [5]. In this algorithm, instead of solving harmonic extensions (9 and 17) for each fine grid node on the boundary, we solve a small number of harmonic extension local problems with random boundary conditions. More precisely, we let
where are independent identical distributed standard Gaussian random vectors on the fine grid nodes of the boundary. The advantage of this algorithm lies in the fact that a much fewer number of snapshot basis functions are calculated, while maintaining accuracy. In addition, we will use an oversampling strategy. This is done to reduce the mismatching effects of boundary conditions imposed artificially in the construction of snapshot basis functions. We will denote the extended coarse grid neighborhood for , by . Here for example, , would mean the coarse grid neighborhood plus all 1 layer of adjacent fine grida of , and so on.
In Fig. 8 and Fig. 9, we show the weighted and errors over time for Case 1 and 2 using the randomized GMsFEM with oversampling using different numbers of multiscale basis functions. The oversampled region is chosen, that is, the oversampled region contains an extra 4 fine grid cells layers around . Here, we use only the fully coupled scheme. We use a snapshot ratio of between the standard number of snapshots and the randomized algorithm. Comparing results from Fig. 8, the randomized algorithm, to Fig. 6, the standard GMsFEM, we observe that the randomized algorithm is slightly less accurate but at the advantage of having less snapshot solutions required.
In Table 5 and 6 we investigate the effect of the oversampling as we increase the number of fine grid extensions for and . We present the data of the randomized snapshots for last time step. We see that oversampling can help to improve the results initially, but the improvements level off as large oversampling domains do not give significant improvement in the solution accuracy. Again the effects of the high contrast of Case 2 can be seen in the data as the oversampling performs slightly worse than in the lower contrast regime.
pressure errors  displacements errors  
without oversampling,  
4  0.05  0.04  0.16  0.22 
8  0.05  0.03  0.16  0.21 
12  0.02  0.02  0.16  0.21 
16  0.004  0.01  0.15  0.21 
with oversampling,  
4  0.05  0.03  0.14  0.21 
8  0.04  0.03  0.12  0.19 
12  0.007  0.01  0.09  0.17 
16  0.002  0.009  0.08  0.16 
with oversampling,  
4  0.05  0.03  0.09  0.17 
8  0.04  0.02  0.06  0.14 
12  0.006  0.01  0.04  0.11 
16  0.001  0.008  0.02  0.08 
with oversampling,  
4  0.05  0.03  0.09  0.17 
8  0.04  0.02  0.06  0.13 
12  0.009  0.01  0.02  0.09 
16  0.002  0.007  0.02  0.07 
pressure errors  displacements errors  
without oversampling,  
4  0.05  0.04  0.36  0.31 
8  0.05  0.03  0.35  0.31 
12  0.02  0.02  0.34  0.31 
16  0.006  0.01  0.34  0.31 
with oversampling,  
4  0.05  0.03  0.33  0.31 
8  0.04  0.03  0.30  0.30 
12  0.01  0.01  0.25  0.27 
16  0.009  0.009  0.22  0.25 
with oversampling,  
4  0.05  0.03  0.28  0.29 
8  0.03  0.02  0.19  0.24 
12  0.007  0.01  0.11  0.20 
16  0.002  0.008  0.07  0.15 
with oversampling,  
4  0.05  0.03  0.24  0.27 
8  0.04  0.02  0.14  0.22 
12  0.01  0.01  0.07  0.17 
16  0.002  0.007  0.06  0.14 
6 Conclusion
Simulating poroelasticity is difficult due the complex heterogeneities and because of the complexity of gridding the flow and mechanics regimes in such media. Therefore, in this paper we developed a Generalized Multiscale Finite Element Method for a linear poroelastic media. We first presented the general poroelasticity framework of Biot and its subsequent solution by fixed stress time splitting methods. Although fully coupled schemes are considered numerically, this splitting lays the framework for the application of the GMsFEM to the decoupled poroelastic equations. We then outline the construction of the multiscale spaces in both fluid and mechanics regimes. The algorithm is then implemented on a single geometry with two different cases of elastic parameters. We show the errors relative to the fine scale solution over time and with varying multiscale basis functions. Finally, we implemented oversampling strategies and randomized boundary conditions when solving for the snapshot space. As in cases of reservoir compaction, the permeability may depend on pressure resulting in a nonlinear relation. In future studies, we will develop a GMsFEM for such nonlinear poroelastic problems.
References
 [1] Y. Efendiev, J. Galvis and T. Hou Generalized Multiscale Finite Element Methods. Journal of Computational Physics, V. 251, pp.116135, 2013.
 [2] Y. Efendiev and T. Hou Multiscale Finite Element Methods: Theory and Applications. Springer, New York, Surveys and Tutorials in the Applied Mathematical Sciences, V. 4, 2009.
 [3] Y. Efendiev, J. Galvis, G. Li and M. Presho Generalized Multiscale Finite Element Methods. Oversampling strategies. International Journal for Multiscale Computational Engineering, 12(6), 2014.
 [4] E. Chung, Y. Efendiev and S. Fu Generalized multiscale finite element method for elasticity equations. To appear in International Journal on Geomathematics, 2014.
 [5] V. Calo, Y. Efendiev, J. Galvis, and G. Li Randomized oversampling for generalized multiscale finite element methods. http://arxiv.org/pdf/1409.7114.pdf, 2014.
 [6] E. Chung, Y. Efendiev and C. Lee Mixed generalized multiscale finite element methods and applications. To appear in Multicale Model. Simul, 2014.
 [7] D. L. Brown, P. Popov, and Y. Efendiev Effective equations for fluidstructure interaction with applications to poroelasticity. Applicable Analysis, 93(4):771790, 2014.
 [8] D. L. Brown, P. Popov, and Y. Efendiev On homogenization of stokes flow in slowly varying media with applications to fluidstructure interaction. GEM  International Journal on Geomathematics, 2(2):281305, 2011.
 [9] A. Mikelic and M. F. Wheeler, Convergence of iterative coupling for coupled flow and geomechanics, Springer, Computational Geosciences, pp. 1–7, 2013.
 [10] J. Kim, H. A. Tchelepi and R. Juanes Stability, accuracy, and efficiency of sequential methods for coupled flow and geomechanics, SPE Journal, 16(2), pp. 249262, 2011.
 [11] J. Kim Sequential methods for coupled geomechanics and multiphase flow. Stanford University, PhD thesis, 2010.
 [12] S. E. Minkoff, C. M. Stone, S. Bryant, M. Peszynska and M. F. Wheeler Coupled fluid flow and geomechanical deformation modeling. Elsevier, Journal of Petroleum Science and Engineering, V. 28, N. 1, pp. 37–56, 2003.
 [13] B. Jha and R. Juanes A locally conservative finite element framework for the simulation of coupled flow and reservoir geomechanics. Springer, Acta Geotechnica, V. 2, N. 3, pp. 139–153, 2007.
 [14] F. Armero and J. C. Simo A new unconditionally stable fractional step method for nonlinear coupled thermomechanical problems. Wiley Online Library, International Journal for numerical methods in Engineering, V. 35, N. 4, pp. 737–766, 1992.
 [15] A. E. Kolesov, P. N. Vabishchevich and M. V. Vasilyeva Splitting schemes for poroelasticity and thermoelasticity problems. Computers & Mathematics with Applications, 67(12), pp. 21852198, 2014.
 [16] C. M. Sayers and P. Schitjens An introduction to reservoir geomechanics. Society of Exploration Geophysicists, The Leading Edge, 26(5), pp. 597601, 2007.
 [17] M. Zoback Reservoir geomechanics. Cambridge University Press, 2010.
 [18] A. Settari and D. A. Wlaters Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction. Society of Petroleum Engineers, Spe Journal, 6(03), pp. 334342, 2001.
 [19] A. Settari and F.M. Mourits A coupled reservoir and geomechanical simulation system. Society of Petroleum Engineers, Spe Journal, 3(03), pp. 219226, 1998.
 [20] V.K. Shrivastava, D. Tran, L.X. Nghiem, and B. F. Kohse Use of Improved Gridding Technique in Coupled Geomechanics and Compositional Reservoir Flow Simulation. Society of Petroleum Engineers, SPE Middle East Oil and Gas Show and Conference, 2011.
 [21] M.A. Biot General theory of three dimensional consolidation. Journal of Applied Physics, 12, pp. 155164, 1941.
 [22] W. E and B. Engquist The heterogeneous multiscale methods. Communications in Mathematical Sciences, 1(1), pp. 87132, 2003.
 [23] A. Abdulle, W. E, B. Engquist, and E. VandenEijnden The heterogeneous multiscale methods. Acta Numerica, 21, pp. 187, 2012.
 [24] T. J. R. Hughes and G. Sangalli Variational multiscale analysis: the finescale Green’s function, projection, optimization, localization, and stabilized methods. SIAM Journal on Numerical Analysis, 45(2), pp. 539557, 2007.
 [25] A. Målqvist and D. Peterseim Localization of elliptic multiscale problems. Mathematics of Computation, 83(290), pp. 25832603, 2014.