Online Coupled Generalized Multiscale Finite Element Method for the Poroelasticity Problem in Fractured and Heterogeneous Media

In this paper, we consider the poroelasticity problem in fractured and heterogeneous media. The mathematical model contains a coupled system of equations for fluid pressures and displacements in heterogeneous media. Due to scale disparity, many approaches have been developed for solving detailed fine-grid problems on a coarse grid. However, some approaches can lack good accuracy on a coarse grid and some corrections for coarse-grid solutions are needed. In this paper, we present a coarse-grid approximation based on the generalized multiscale finite element method (GMsFEM). We present the construction of the offline and online multiscale basis functions. The offline multiscale basis functions are precomputed for the given heterogeneity and fracture network geometry, where for the construction, we solve a local spectral problem and use the dominant eigenvectors (appropriately defined) to construct multiscale basis functions. To construct the online basis functions, we use current information about the local residual and solve coupled poroelasticity problems in local domains. The online basis functions are used to enrich the offline multiscale space and rapidly reduce the error using residual information. Only with appropriate offline coarse-grid spaces can one guarantee a fast convergence of online methods. We present numerical results for poroelasticity problems in fractured and heterogeneous media. We investigate the influence of the number of offline and online basis functions on the relative errors between the multiscale solution and the reference (fine-scale) solution.


Introduction
The poroelasticity problem has applications in many different engineering disciplines such as petroleum engineering, agricultural science, and biomedicine [1][2][3]. The mathematical model of the poroelasticity problem is described by a coupled system of equations for fluid pressure and displacements known as the Biot model [4][5][6][7]. The poroelasticity equations describe the fluid flow in the pores and the medium's stress-strain state [8][9][10]. Numerical solutions of flow and mechanics systems have often been considered separately using splitting schemes [11][12][13][14]. However, a fully coupled scheme may be preferable.
Heterogeneous and fractured microstructures and high-contrast physical properties are the key characteristics of porous media and modern composite materials. It is well known that almost all materials used in modern life and industry, both produced and found in nature, are heterogeneous and multicomponent. These materials have rich and complex internal structures. Relevant examples can be given from all science fields, such as heterogeneous (composite) solids, mixtures, multicomponent liquids, soils and rocks, and biological tissues. The internal structure or microstructure plays a crucial role in understanding and controlling the macroscopic (continuum) behavior of such materials [15][16][17]. In the reservoir simulation, the fluid flow and mechanics in the fractured porous media play an essential role. Fracture networks have a complex structure and have a significant impact on the flow processes in poroelastic media [18][19][20][21]. There are various fracture models, each of which has benefits and drawbacks and application areas [22][23][24]. If a medium consists of a relatively low number of dominant fractures, one can use the discrete fracture model (DFM) [24][25][26][27]. This model requires the construction of conforming grids with fracture geometry. The presence of fracture networks and heterogeneous properties leads to complications to the effective simulation due to the complexity of scales. In numerical simulations of the poroelasticity problems, these multiscale properties should be accurately accounted for by high-resolution simulations and require many computational resources, since the fracture and heterogeneity resolution in approximation increases the number of unknowns.
Multiscale methods are widely used to construct a reduced-order model [28][29][30]. The generalized multiscale finite element method (GMsFEM) for solving flow problems in fractured media is considered in [31,32]. In [33,34], we considered the coarse-grid problem construction for the poroelasticity problem in heterogeneous media. The multiscale finite volume method for the solution of the poroelasticity problem is presented in [9,[35][36][37][38]. We also note that many global model reduction techniques have successfully been applied to solve challenging problems [39,40]. In this paper, our main focus is on local model reduction techniques; in particular, on online approaches. The online generalized multiscale finite element method is presented in [41][42][43]. A new method for solving flow and poroelasticity problems in fractured and heterogeneous porous media is presented in [18,[44][45][46][47].
This work presents the online coupled generalized multiscale finite element method to solve the poroelasticity problem in fractured and heterogeneous media. For fine-grid approximation, we apply the finite element method with the discrete fracture model. The coarse-grid approximation is based on the generalized multiscale finite element method (GMsFEM). We present the construction of the offline and online multiscale basis functions. The offline multiscale basis functions are precomputed for the given heterogeneity and fracture network geometry. The construction of the offline basis functions is given separately for pressure and displacement fields and involves the solution of local spectral problems on the snapshot space for flow and mechanics. We use current information about the local residual and solve coupled poroelasticity problems in local domains to construct the online basis functions. The offline basis functions are selected using eigenvectors that correspond to dominant eigenvalues. The latter introduces a cut-off in the spectrum. We choose dominant modes based on the smallest eigenvalues following [31], as the convergence rate (offline and online) controlled by the reciprocal of the smallest eigenvalue of the corresponding eigenvector is not selected in the construction of multiscale spaces.
The online basis functions are used to enrich the offline multiscale space and reduce error using residual information. The online GMsFEM was introduced in [41] for steadystate flow equations. In this work, the authors show that only with an appropriate offline space does the online method converge (if we iterate on the residual), and the number of iterations is independent of the contrast. This does not hold if the offline space does not contain a sufficient number of basis functions. The online convergence is related to the cutoff in the spectral problem. The extensions of these ideas to more complex time-dependent poroelasticity problems in heterogeneous and fractured media are considered in this paper. In particular, we update online basis functions at certain time instants and study how the offline space affects the accuracy of the method. Numerical results show the efficiency of the presented method due to the coupled construction of the online residual-based multiscale basis functions.
We present numerical results. In our tests, we consider heterogeneous permeability fields and fractures. We use different numbers of offline basis functions and enrich the offline space with online multiscale basis functions at certain time instants. Our numerical results show that with fewer updates and online basis functions, the accuracy of multiscale simulations can be improved.
The work is organized as follows. Section 2 presents a mathematical model and a fine-grid approximation of the poroelasticity problem in a fractured and heterogeneous medium. In Section 3, we present a coarse grid approximation using the generalized multiscale finite element method, where we describe the calculation of the de-coupled offline multiscale basis functions for flow and mechanics and describe the construction of the coarse-grid system. Then, we present an algorithm for constructing the coupled online multiscale basis functions in Section 4. Numerical results for two-dimensional model poroelasticity problems in the fractured and heterogeneous medium are presented in Section 5. Finally, we present conclusions.

Problem Formulation and Approximation on a Fine Grid
Let Ω ⊂ R d be the domain for the matrix and γ ⊂ R d−1 be the domain for the fracture network, where d denotes a geometrical dimension of the problem and d = 2 for a two-dimensional case. The mathematical model is described by a coupled system of equations for pressure in a porous matrix p m ∈ Ω, pressure in a fracture network p f ∈ γ, and displacements u ∈ Ω [46]: where k α = κ α /ν, (α = m, f ), κ m and κ f are the matrix and fracture permeabilities, f m and f f refer to the matrix and fracture source terms, α m and α f are the matrix and fracture Biot coefficients, M m and M f are the matrix and fracture Biot moduli, r m f and r f m are the transfer terms between porous matrix and fractures, ν is the fluid viscosity, and σ is the stress tensor. The relation between the stress and strain tensors is given as where ε is the strain tensor, and λ and µ are the Lamé coefficients. We consider a system of Equation (1) with the following initial conditions: and boundary conditions where r is a mass transfer coefficient, s is an external pressure, Γ L ∪ Γ R ∪ Γ T ∪ Γ B = ∂Ω; Γ L , Γ R , Γ B , Γ T denote left, right, bottom, and top boundaries, respectively; and Γ γ = Γ L ∩ ∂γ (see Figure 7 for illustration).
Variational formulation. For approximating the system of Equation (1) with boundary conditions (3), we use a finite element method. We define the following functional spaces The variational formulation of the poroelasticity problem in fractured media can be written as follows.
where the bilinear and linear forms are as follows: where Π m f is the projection operator from W m to W f defined as [48] for details). Discrete form. Next, we construct an unstructured fine grid T h that explicitly resolves the fracture geometry and use the discrete fracture model [19,24]. For the approximation by time, we use an implicit discretization of the equations with a time step τ. We have the following discrete system in matrix form on the fine grid for where andy h is the solution from the previous time step.
Here, In this work, we use a modified DFM approach to simplify the matrix construction and consider the case when α f = 0. The nodes of the fracture element and the matrix element are coincident at the interface because the unstructured fine grid T h explicitly resolves the fracture geometry. For simplicity, we suppose p h = p h m = p h f on the fractures and p h = p h m on the porous matrix. Using the superposition principle (by summing the equations for fractures and matrix nodes), we obtain the following coupled system of equations for in a grid construction on a system of fractures, it is important to consider the intersections of fractures and the intersections of the intersection segments. The grid-building algorithm consists of the following steps: • The sequential division of the intersection segments and the fracture boundaries by points into segments of a given size, which are subsequently used as mesh faces; • Uniform filling of fracture regions by points with a restriction on the proximity of internal points to the boundary and points on intersection segments; • Construction of triangulation for each fracture according to the given frameworks (points, segments).

Multiscale Method for a Coarse-Grid Approximation on Offline Space
We use the generalized multiscale finite element method (GMsFEM) to construct the coarse-grid approximation of the poroelasticity problem in a fractured and heterogeneous medium. In this computational algorithm, the first four steps are offline (preprocessing) steps for a given fracture geometry and heterogeneity.
Offline stage • Coarse grid and local domain construction. • The solution of the local problems with different boundary conditions to construct a snapshot space in each local domain.
• An offline space construction via the solution of the local spectral problems on the snapshot space. • Generation of the projection matrix.
Next, we move on to the online stage that includes the construction of the coarse-grid system using the precalculated projection matrix with offline multiscale basis functions.
Online stage • Construction of the coarse grid system using the projection matrix. • Solving the problem on the coarse grid at the current time step. • Moving to the next time step.
In this work, we use a continuous Galerkin approximation on the coarse grid. We define local domain ω l (see Figure 5) for multiscale basis functions as a combination of the several coarse-grid cells that share the same coarse-grid node (l = 1, ..., N H v , N H v is the number of coarse-grid vertices), see Figure 1. We denote the coarse grid with T H . We start with the construction of the offline space, where the generation of the offline basis functions for displacements and pressure are given separately. The offline basis construction contains two steps: (1) snapshot space construction and (2) solution of the local spectral problem on snapshot space.
It is possible to speed up the construction of multiscale basis functions. First, we need to construct local domains. Then, we construct multiscale basis functions for each domain, dividing them among processors. Parallelization works seamlessly because local domains are independent of each other.
Multiscale basis functions for pressure. To construct a snapshot space, we solve the following local problem in domain where with γ ω l = γ ∩ ω l and δ j i being the Kronecker delta function for j = 1, .., N is the number of fine grid nodes on the ∂ω l ).
We define a snapshot space for pressure in fractured media as follows: Next, we solve the following local spectral problem on the snapshot space where A l,snap p Here for matrices, we have We choose eigenvectors (see Figures 2 and 3)φ l,j (j = 1, .., M l,p ) corresponding to the first smallest M l,p eigenvalues. The eigenvalues are ordered in increasing order as are the corresponding eigenfunctions. The most dominant mode is represented by the eigenfunctions with the smallest eigenvalue. We will choose the eigenfunctions corresponding to small eigenvalues because the convergence rate of GMsFEM is a reciprocal of the smallest eigenvalue for which the corresponding eigenvector is not included in multiscale space. Furthermore, we multiply to the linear partition of unity functions χ l for obtaining conforming basis functions W ms = span{φ l,j , l = 1, ..., N H v , j = 1, ..., M l,p }, where φ l,j = χ lφl,j . Multiscale basis functions for displacements. For construction of a snapshot space, we solve the following local problem in domain ω l . Find Ψ l,j ∈ V h (ω l ) such that where For the construction of the multiscale basis, we solve the following local spectral problem on the snapshot space: We choose eigenvectors (see Figure 4) Φ Finally, we obtain the following reduced-order model whereC H = RCR T ,Ã H = RÃR T , andF H = RF. After obtaining a coarse-scale solution, we can reconstruct a fine-scale solution The coarse-scale system size is N H = ∑

Online Enrichment of Multiscale Space
To improve the accuracy for the presented multiscale approximation, we present the construction of the online residual-based multiscale basis functions. We compute the coupled online basis functions after solving the coarse-scale system on the offline space using residual information on the online stage.
Online stage • Construction of the coarse-grid system using the projection matrix. • Solving the problem on the coarse grid at the current time step. • Multiscale space enrichment by calculation of the online basis functions. We enrich the offline space and update the projection matrix using the obtained online basis functions. After that, we repeatedly solve the current time step problem on the coarse grid to update the solution. • Moving to the next time step.
To construct the local residual-based online multiscale basis functions (see Figures 5 and 6), we solve the following local problems in each ω l . Find (Υ l,p with Here, we have the following bilinear forms and the right-hand side is based on the local residual information where the subscript ms denotes the multiscale solution.
where k is the number of online iterations for the current time step. We will update online basis functions for some time steps. Next, we present an algorithm for the multiscale method with the online residualbased multiscale basis functions. Let R o f f be the projection matrix constructed using offline multiscale basis functions

Multiscale algorithm with online enrichment
• Define projection matrix R =Ř for the current time step n with R = R o f f for n = 0 (Ř is the projection matrix from the previous time step). • Construct and solve the coarse-scale problem at the current time step.
and y ms,k = (R k ) T y H,k . For the projection matrix, we have R k = R o f f for k = 0 and In general, we can add online basis functions adaptively only for some local domains with a large residual [41,42]. Next, we present numerical results for heterogeneous and fractured poroelastic media.

Numerical Results
This section presents numerical results for coupled flow and mechanics in heterogeneous and fractured media. Calculations are performed on a computation domain, which is shown in Figure 7 for Ω = [0, 50] 2 . The fine grid consists of 12,699 vertices, 24,996 cells, and 37,694 facets, and the coarse grid contains 121 vertices, 100 cells, and 220 facets. The fine grid is conformed with the coarse-grid edges and resolves fracture geometry at the grid level.
We consider two cases For initial conditions, we set p 0 = 0 and u 0 = 0. For boundary conditions, we set (3) with r = 10 4 and s = 1.
We consider an isotropic permeability field and elasticity tensor with where ν = 0.3 is Poisson's ratio and E is Young's modulus. Heterogeneous Young's modulus and heterogeneous permeability are presented in Figure 8. For constructing the unstructured grid that conformed with discrete fractures and coarse grid facets, we used GMSH software [49]. The implementation of the fine grid solver and multiscale method was based on the FEniCS library [50].  To compare the results, we use the fine-grid solution as a reference solution and calculate the relative L 2 norm and H 1 semi-norm of errors between the multiscale solution and the reference solution (p ms , u ms ) is the multiscale solution, and (p, u) is the reference (fine grid) solution. We note that, since matrix R needs to be re-created, when we add the online basis functions, we add/update them for some selected time steps, in our case 5 and 10. We note that when we add new online basis functions at some selected time step based on current residuals, we remove previously calculated online basis functions and keep them until we update the online basis functions. More online updates initially help to reduce the error due to the initial condition. We use multiscale basis functions from the offline space as initial basis functions.
In Figures 9 and 10, we present the distributions of pressure and displacement in X and Y directions at the final time for Cases 1 and 2, respectively. The first row represents pressure distribution (p), and the second and third rows are displacements in the X and Y directions (u x and u y ). In each figure, we have three columns: the first column represents the reference (fine grid) solution, the second column is the multiscale solution using M o f f offline basis functions for pressure and displacements (M on = 0), and the third column is the multiscale solution using one online residual-based basis function (M on = 1). In Case 1, the size of the system for the reference solution on the fine grid is DOF h = 12,944. For the multiscale solution, using 2 offline basis functions, we have 0.605% and 15.28% of L 2 and H 1 errors with DOF H = 726 (5.608% from DOF h ) for pressure. After adding one online basis function at every 5 time steps, we reduce the errors to 0.052% and 2.354% for L 2 and the H 1 norm with DOF H = 1089 for pressure. In Case 2, the size of the system for the reference solution on the fine grid is DOF h = 12,944, and we have 3.297% and 22.77% of L 2 and H 1 errors using 4 offline basis functions with DOF H = 1452 (11.217% from DOF h ) for pressure. By adding one online basis function at every 5 time steps, we have 1.910% and 16.95% of L 2 and H 1 errors with DOF H = 1815 for pressure.
In Tables 1 and 2, we present the relative L 2 and H 1 errors without online basis and after adding 1 or 2 online basis functions for Case 1 and Case 2. We have calculated errors for different numbers of offline multiscale basis functions (M o f f ) and considered updating the online basis functions at every 5 and 10 time steps. For the poroelasticity problem in fractured media in Case 2, we obtain large errors with M o f f = 1, 2, and online enrichment does not work. This happens because we did not take a sufficient number of the offline basis functions (continuums) to capture local flows between the porous matrix and fractures [19].
For M o f f = 4, 8, 12, the presented multiscale method provides good results for fractured media in Case 2. We observe that most of the errors vanished in both cases after using just 8 offline multiscale basis functions. With a sufficient number of offline basis functions by adding only one online basis function, we obtain a significant error reduction with a small system size. We observe that the errors decrease as the number of online basis functions increases, as expected.
We present relative L 2 error dynamics in Figures 11 and 12 for Case 1 and Figures 13 and 14 for Case 2. For Case 1, we take two offline multiscale basis functions and four offline multiscale basis functions for Case 2. If we add online multiscale basis functions, the errors reduce significantly. We observe that the errors reduce over time. From the presented results, we see that updating the online basis functions at every 10 time steps is enough to obtain a good solution for both cases.

Conclusions
We have presented a residual-based online generalized multiscale finite element method for the poroelasticity problem in fractured and heterogeneous media. As a reference solution to compare the results of the presented multiscale method, we have used a simulation on the fine grid. For fine-scale approximation, we have used the finite element method on unstructured grids that resolves fractures on the grid level, where the discrete fracture model is used to represent flow in fractures. We have presented a coarse-grid approximation technique based on the generalized multiscale finite element method. The construction of the two types (offline and online) of the multiscale basis functions was presented. The online basis functions were constructed locally using current residual information by solving the coupled poroelasticity equations. Numerical results for a two-dimensional formulation of the poroelasticity problem were presented for two cases: Case 1-heterogeneous media and Case 2-fractured heterogeneous media. We presented relative errors between the multiscale solution and the reference (fine-scale) solution for a different number of offline and online multiscale basis functions. We showed that online basis functions can significantly reduce errors for pressure and displacements due to coupled construction and using the current residual information. Cases with online basis functions updated at every 5 and 10 time steps were considered. Numerical results demonstrate the efficiency of the presented algorithm with online enrichment of the multiscale space for solving poroelasticity problems in fractured and heterogeneous media.

Conflicts of Interest:
The authors declare no conflict of interest.