Generalized Multiscale Finite Element Method for Elastic Wave Propagation in the Frequency Domain

: In this work, we consider elastic wave propagation in fractured media. The mathematical model is described by the Helmholtz problem related to wave propagation with speciﬁc interface conditions (Linear Slip Model, LSM) on the fracture in the frequency domain. For the numerical solution, we construct a ﬁne grid that resolves all fracture interfaces on the grid level and construct approximation using a ﬁnite element method. We use a discontinuous Galerkin method for the approximation by space that helps to weakly impose interface conditions on fractures. Such approximation leads to a large system of equations and is computationally expensive. In this work, we construct a coarse grid approximation for an effective solution using the Generalized Multiscale Finite Element Method (GMsFEM). We construct and compare two types of the multiscale methods—Continuous Galerkin Generalized Multiscale Finite Element Method (CG-GMsFEM) and Discontinuous Galerkin Generalized Multiscale Finite Element Method (DG-GMsFEM). Multiscale basis functions are constructed by solving local spectral problems in each local domains to extract dominant modes of the local solution. In CG-GMsFEM, we construct continuous multiscale basis functions that are deﬁned in the local domains associated with the coarse grid node and contain four coarse grid cells for the structured quadratic coarse grid. The multiscale basis functions in DG-GMsFEM are discontinuous and deﬁned in each coarse grid cell. The results of the numerical solution for the two-dimensional Helmholtz equation are presented for CG-GMsFEM and DG-GMsFEM for different numbers of multiscale basis functions.


Introduction
Understanding the complex processes in fractured media is necessary in many real world applications, for example, in exploring and developing hydrocarbon reservoirs. Fractured media are characterized by a complex fracture distribution and presence of the fractures at multiple scales [1,2]. Moreover, fractures are characterized as thin and long interfaces. Numerical methods for the solution of the problems in fractured media are being actively developed in recent years. In Reference [3] the authors developed finite difference schemes for modeling wave propagation in a fractured medium. Fractured media are derived when the fracture interfaces are aligned with the boundaries of the finite-difference grid. For flow problems in fractured media, the multiscale finite volume method on unstructured grids for discrete fracture model (MS-DFM) is presented in Reference [4]. The use of unstructured grids requires rather complicated grid generators and significantly increases the run time, also, its use is undesirable when the number of fractures is large. In Reference [5] the author proposed an approach called Chimera Grid Approach. This technique avoids the use of curvilinear or unstructured grids, replacing them with a set of arbitrary rotated rectangular grids. In Reference [6], the authors proposed an approach to modeling the propagation of seismic waves in a medium containing subvertical fractured inhomogeneities using the grid-characteristic method on structural grids without constructing additional grids. To simulate wave propagation in fractured media, the discontinuous Galerkin (DG) method is widely used. The main disadvantage of the DG method is the excessive use of memory due to a large number of degrees of freedom. But in Reference [7] the authors propose an Enriched Galerkin finite element method (EGM) for elastic wave propagation. EGM formulated by enriching the conforming continuous Galerkin finite element method with piecewise constant functions.
In this work, we consider elastic wave propagation in fractured media. The Linear Slip Model (LSM) is used to simulate the influence of the fracture on the distribution of the seismic field, where the stress components are proportional to the displacement [1,2]. LSM requires constructing a fine grid that resolves all fracture interfaces on the grid level, which leads to a very large system of equations. To construct a suitable mathematical model for seismic wave propagation, multiscale methods or homogenization techniques are used. Homogenization or effective media theories replace a heterogeneous rock with an equivalent homogeneous rock, where the corresponding properties are calculated under the static-deformation regimes. An overview and comparison of the existing effective media theories for fractured formations is presented in Reference [8]. The authors consider a three-dimensional formulation and use direct finite-element simulations in order to calculate the effective properties of the fractured media. Direct computational studies provide an independent verification of theoretical predictions and can be used to examine realistic fracture models that violate the conventional assumptions, where fractures have notoriously irregular shape and can form interconnected networks. On the other hand, this approach is computationally expensive for complex fracture distributions. Moreover, the errors of the homogenization methods depend on the coarse-grid mesh.
In the multiscale methods, a central research problem is constructing the coarse grid approximation for faster results, where multiscale basis functions are computed on a fine grid to capture the influence of fractures and other heterogeneity [9][10][11][12][13][14]. We consider seismic waves in fractured media and construct multiscale basis functions for coarse grid simulations in the two-dimensional formulation [15][16][17][18][19]. These multiscale basis functions can capture the influence of the fractures on a coarse grid and reduce the number of unknowns in calculations. Multiscale basis functions are constructed using the solution of the local spectral problem on a fine grid in each coarse grid cell. In our previous work, we presented GMsFEM for elastic wave propagation in time domain in fractured media [15,20].
In this paper, we construct two types of multiscale spaces to solve the Helmholtz equation for fractured media. The numerical implementation is based on the global projection approach to construct a coarse-grid approximation. The local spectral problems that we used for selecting the dominant modes are frequency-independent. The presented method has much in common with classical homogenization methods and, therefore, the applicability is limited by low frequencies.
We present numerical results and consider several test cases with different fracture configurations. Numerical simulations show that the method is accurate for complex cases and can significantly reduce the dimension of the system.
The work is organized as follows-in Section 2, we provide the mathematical model of the Helmholtz problem in fractured media and construct a fine grid approximation using the finite element method in Section 3. In Section 4, we present a coarse grid approximation using the Generalized Multiscale Finite Element Method (GMsFEM) algorithm, where we describe constructing of the two types of multiscale basis functions-(1) Continuous Galerkin Generalized Multiscale Finite Element Method (CG-GMsFEM) with continuous multiscale basis functions and (2) Discontinuous Galerkin Generalized Multiscale Finite Element Method (DG-GMsFEM) with discontinuous multiscale basis functions. Finally, numerical results are presented in Section 5.

Problem Formulation
We consider a problem in a fractured medium. Let Ω ⊂ R d is the bounded computational domain and γ ⊂ R d−1 is the fracture interface, d = 2 (see Figure 1 for illustration). The elastic wave propagation is described by Helmholtz equation in the computational domain Ω [21,22] where ω is the frequency, σ is the stress tensor, ρ is the density and f is the source function.
We have linear elastic stress-strain constitutive relation where I is the unit tensor, λ and µ are the Lame parameters. In this paper, we focus on the two-dimentional formulation of the problem. Since the problem is considered in fractured media, for numerical simulations of the elastic wave equation, the linear-slip model (LSM) is applied on the fracture interface γ [23,24]. Namely, we assume the fractures have a vanishing width across which the traction is considered continuous. Following the linear-slip model, we have a linear relation between the traction vector and the magnitude of the discontinuity in the displacement field as follows where [u] is the jump of the displacement field at the fracture, σn is the traction vector at the surface of the fracture γ and Z is the fracture compliance matrix. In the two-dimensional isotropic case, the compliance matrix is diagonal and positive definite where z 1 and z 2 are the normal and tangential compliance. In the numerical simulations we will use a scalar fracture compliance, with z = z 1 = z 2 and Z = zI, where I is the identity matrix [23,24]. In the computations we use the first order absorbing boundary condition, because the energy of waves needs to be absorbed in artificial boundaries in order to avoid spurious reflections caused by the finite computational domain [25] where A = n 1 n 2 −n 2 n 1 c p 0 0 c s n 1 −n 2 n 2 n 1 .
Here n = (n 1 , n 2 ) is the outward normal to the boundary and c p = (λ + 2µ)/ρ, c s = µ/ρ are the S-and P-wave velocities. We will impose this absorbing condition in the variational formulation of the problem.

Fine Grid Approximation
We use the Interior Penalty Discontinuous Galerkin (IPDG) finite element method for the fine-grid approximation. IPDG method allows for discontinuities in the displacement field to simulate fractures with the linear-slip model [26][27][28]. We construct a fine grid for fracture domain Ω that resolve fracture interface γ on the grid level.
Let T h be a triangulation of the computational domain Ω and h be the fine mesh size. We define Γ h as the set of all the interior faces between the elements T h , Γ c ⊂ Γ h be the subset of all faces, where the displacement field is continuous and Γ γ ⊂ Γ h be the subset of facet that represent fractures, Γ b is a subset of faces on the boundary and Γ h = Γ c ∪ Γ γ ∪ Γ b . Let e ∈ Γ h be the edge between the elements ι 1 and ι 2 , then the average {·} and jump [·] of a vector function u on e are given by The variational formulation of the elastic wave equation using the interior penalty discontinuous Galerkin method in fractured media is defined as follows: where ζ is the penalty parameter, τ(u) = σn is the traction vector, V is a suitable finite dimension function space. Here u is the complex valued function and u = ∑ j u j φ j , φ j are linear basis functions for the fine scale approximation. By definition of the following bilinear and linear forms we obtain the following formulation We can write the complex valued problem in matrix form where K h is the stiffness matrix, M h is the mass matrix and B h is the boundary mass matrix.

Multiscale Method on the Coarse Grid
In this section, we describe construction of the multiscale basis functions and coarse grid approximation [10,15,17,[29][30][31]. For construction of the coarse grid approximation, we will use a projection approach and define a projection matrix R. We consider two types of the multiscale basis functions: (1) CG-GMsFEM with continuous multiscale basis functions and (2) DG-GMsFEM with discontinuous multiscale basis functions in coarse cells. By the algorithmic point of view, the differences between two considered methods are concentrated on the definition of the local domains and construction of the multiscale basis functions. After construction, we collect multiscale basis functions into the projection matrix R. Therefore, we will define construction of the two projection matrices R CG and R DG .
The main steps of the multiscale method is similar for both approaches (CG-GMsFEM and DG-GMsFEM). We have following steps in the multiscale computational algorithm: 1.
construction of the coarse grid and local domains; 2.
construction of the multiscale basis functions by the solution of the local eigenvalue problem in each local domains; 3.
construction of the projection matrix R (from fine grid to coarse grid) using computed multiscale basis functions; 4.
construction of the fine grid system and projection to the coarse grid using matrix R; 5.
solution of the reduced order model and reconstruction of the fine grid solution.
We start with definition of the coarse grid. Let T H be the coarse grid partitioning of the domain Ω, T H = ∪ N c i=1 K i with mesh sizes H h > 0, where K i is the coarse cell and N c is number of coarse cells. We note that, the coarse grid edges are conforming with the fine grid. In this work, we use a structured coarse grid with quadrilateral coarse cells. In general, shape of the coarse grid can be complex (unstructured coarse grid) and can be constructed, for example, using traditional mesh partitioning.
Next, we define local domains. We have • for CG-GMsFEM, the local domain ω i are defined as a coarse neighborhood that contains four coarse grid quadrilateral cells around coarse grid node, where i = 1, . . . , N v and N v is the number of the coarse grid nodes (see Figure 2 for illustration); • for DG-GMsFEM, the local domain is the coarse grid cell K i , where i = 1, . . . , N c and N c is the number of the coarse grid cells.
After that, we solve local eigenvalue problems in each local domains in order to construct the multiscale basis functions and using them we define a projection matrices R CG and R DG . The form of the eigenvalue problems will be presented later for CG-GMsFEM and DG-GMsFEM.
Finally, the multiscale reduced order discrete system can be calculated by projecting the fine-scale matrices onto the coarse grid with the global projection matrix R = R CG or R = R DG where U H is the complex valued solution on multiscale space, M H and K H are the coarse-scale mass and stiffness matrices, and B H is the coarse-scale boundary mass matrix After calculation of the coarse-scale solution U H , we can recover the fine-scale solution by U ms = R T U H . Next, we consider construction of multiscale basis functions based on two approaches-continuous Galerkin method and discontinuous Galerkin method.

Multiscale Basis Functions for CG-GMsFEM
In CG-GMsFEM, the local domain ω i is obtained by the combining all the coarse cells around one vertex of the coarse grid. To construct the multiscale basis functions, we start with the solution of the following local eigenvalue problem in ω i where where ψ ω i j = χ i φ ω i j (j = 1, . . . , M) and χ i is the piece-wise bi-linear shape function on the coarse grid that equals to 1 at the coarse vertex x i , equals to 0 at all other coarse vertices.
The global projection matrix defined as follows

Multiscale Basis Functions for DG-GMsFEM
In contrast to the continuous Galerkin method, the discontinuous Galerkin approach treats local domain as a coarse cells K i . For DG-GMsFEM, we construct two local multiscale spaces (boundary and interior) by solution of the local eigenvalue problems on each coarse grid cell K i ∈ T H . Boundary and interior basis functions are both defined in the local domain K i up to fine grid resolution, φ

o (indices b and o refer to the boundary and internal basis, respectively). Multiscale basis functions differ by the definition of the local spectral problem.
To construct the boundary multiscale basis functions, we solve following spectral problem in K i : where To construct a boundary multiscale basis functions, we select the first M b eigenvectors φ , and define following local projection matrix The interior multiscale basis functions are defined to capture interior eigenmodes for K i and use following spectral problem to identify the important modes where with homogeneous Dirichlet boundary conditions. We select the first M o eigenvectors φ The global projection matrix defined as follows The complete stability and convergence analysis of the GMsFEM was presented in References [17,30]. We expect a similar behavior to the presented algorithms. The error of both multiscale methods depends on the coarse-scale mesh size (H) and the number of multiscale basis functions (Λ = 1/λ M+1 ) (see Theorem 1 and 3 in Reference [17]).

Numerical Results
In this section, we present the numerical results for the coarse-scale approximation using multiscale basis functions. We consider Continuous Galerkin (CG) and Discontinuous Galerkin (DG) methods in GMsFEM. The basis functions of the multiscale space are constructed following the procedure described above.
For numerical simulations, we use the following parameters. Computational domain is presented in Figure 3 and have dimensions Ω = [0, L x ] × [0, L y ] with L x = L y = 500 [m]. In Figure 4, we show a fine grid (green color) and coarse grid (blue color) for two test geometries with different length of fractures: • Geometry 1: Fine grid contains 16,077 vertices and 31,752 triangle elements.
Note that, the fine grids are unstructured grids that resolve the fractures. Coarse grid is uniform and contains 441 vertices and 400 rectangular elements.  To compare a results, we calculate following relative errors in % where u ms and u are multiscale solutions and reference solutions (fine grid solution). In Figures 5 and 6, we present real parts Re(u) of solution for Geometry 1 and 2 with f 0 = 15, where multiscale solver constructed using DG-GMsFEM. Reference (fine-grid) solution is shown on the top, and multiscale is on bottom. From the figures we can see that the length of fractures affect the solution. If we take a longer crack length, then the waves are reflected more strongly. The fine grid system that used to calculate reference solution has size DOF f = 190,512 for Geometry 1 and   Relative errors are shown in Tables 1 and 2 for two approaches proposed above and for two geometries (Figure 3). DOF c is number of unknowns in multiscale solver, M is the number of multiscale basis functions for CG-GMsFEM, and M o and M b are the numbers of multiscale interior and boundary basis functions for DG-GMsFEM, respectively. We observe that the larger number of multiscale basis functions can reduce errors.

Conclusions
We construct reduced order model using Generalized Multiscale Continuous Galerkin Finite Element Method (CG-GMsFEM) and Generalized Multiscale Discontinuous Galerkin Finite Element Method (DG-GMsFEM). We present numerical results for domain with fractures for two approaches to demonstrate efficiency of presented methods. Numerical results show that the presented multiscale methods give a good approximation of the solution and reduce the size of the system for both types of methods (CG-and DG-GMsFEM). We observe a good reduction of the L 2 errors for both methods when the number of the multiscale basis functions increases. In order to further error reduction, we will add an online residual-based multiscale basis function and oversampling techniques in future works. Moreover, more complex mathematical models and multiscale methods for their solutions will be considered in future works.