{'fraction_non_alphanumeric': 0.031884517766497464, 'fraction_numerical': 0.007336611675126904, 'mean_word_length': 4.797011494252874, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiscale Extended Finite Element Method for Deformable Fractured Porous Media

Deformable fractured porous media appear in many geoscience applications. While the extended finite element (XFEM) has been successfully developed within the computational mechanics community for accurate modeling of the deformation, its application in natural geoscientific applications is not straightforward. This is mainly due to the fact that subsurface formations are heterogeneous and span large length scales with many fractures at different scales. In this work, we propose a novel multiscale formulation for XFEM, based on locally computed basis functions. The local multiscale basis functions capture the heterogeneity and discontinuities introduced by fractures. Local boundary conditions are set to follow a reduced-dimensional system, in order to preserve the accuracy of the basis functions. Using these multiscale bases, a multiscale coarse-scale system is then governed algebraically and solved, in which no enrichment due to the fractures exist. Such formulation allows for significant computational cost reduction, at the same time, it preserves the accuracy of the discrete displacement vector space. The coarse-scale solution is finally interpolated back to the fine scale system, using the same multiscale basis functions. The proposed multiscale XFEM (MS-XFEM) is also integrated within a two-stage algebraic iterative solver, through which error reduction to any desired level can be achieved. Several proof-of-concept numerical tests are presented to assess the performance of the developed method. It is shown that the MS-XFEM is accurate, when compared with the fine-scale reference XFEM solutions. At the same time, it is significantly more efficient than the XFEM on fine-scale resolution. As such, it develops the first scalable XFEM method for large-scale heavily fractured porous media.


Introduction
Subsurface geological formations are often highly heterogeneous and heavily fractured at multiple scales. Heterogeneity of the deformation properties (e.g. elasticity coefficients) can be of several orders of magnitude which occurs at fine scale (cm) resolution. The reservoirs also span large scales, in the order of kms. Numerical simulation of mechanical deformation for such complex systems is necessary to optimise the geo-engineering operations [1,2], and assess their safety and manage the associated risks (e.g. fracture propagation, fault slip and induced seismicity). Though being crucially important, simulation of these systems are beyond the scope of classical numerical schemes.
Presence of highly heterogeneous coefficients with high resolution within largescale domains has been systematically addressed in the computational geoscience community through the development of multiscale finite element and finite volume methods [3][4][5][6]. Recent developments also include mechanical deformation coupled with fluid pore pressure dynamics [7][8][9][10][11]. In presence of many fractures, however, the complexity of the computational model increases significantly. As such, development of a robust multiscale strategy for deformation of heavily fractured porous media, which also allows for convergent systematic error reduction [12][13][14], is of high interest in the geoscience community.
The presence of fractures within the computational domain can be included explicitly by two approaches of (1) unstructured grid and (2) immersed or embedded methods. The unstructured grid approach [15][16][17] generates a discrete computational domain in which fractures are always at the interfaces of elements. This allows for convenient treatment of their effect, however, at the cost of complex meshing. The complex mesh generation for three-dimensional (3D) large scale domains with many fractures is challenging, specially when fractures dynamically extend their geometries. On the other hand, the immersed or embedded approach allows for independent grids for matrix block and fractures, by introducing enrichment of the discrete connectivity (for flow) and shape functions (for mechanics) [18][19][20][21][22]. These enriched formulations are aimed at representing discontinuities within the overlapping matrix cell, without any adjustment nor refinement of the grid [23]. The enrichment strategy for modeling deformation using finite-element schemes in presence of fractured media are referred to as 'extended finite element (XFEM)' methods.
XFEM enriches the partition of unity (PoU) [24] by introducing additional degrees of freedom (DOF) at the existing element nodes. There exists sets of enriched functions to capture the jump discontinuity in the displacement field, when the fracture element cuts through the entire cell, and the tip when a fracture ends within the domain of an element (i.e. its tip is inside an element). [25][26][27][28][29]. For these jump and tip scenarios, additional shape functions are introduced which are multiplied by the original shape functions and supplement the discrete displacement approximation space.
When it comes to geoscience applications, the XFEM is not an attractive method, due to its excessive additional degrees of freedom to capture the many fractures. As such one has to develop a scalable approach, in order to allow for accurate yet efficient application of XFEM to simulate deformation in geological formations. This paper develops a multiscale XFEM (referred to as MS-XFEM) which offers a scalable efficient strategy to model large-scale fractured systems. MS-XFEM imposes a coarse mesh on the given fine-scale mesh. The main novel idea behind MS-XFEM is to use XFEM to computationally solve for local coarse-scale (multiscale) basis functions. These basis functions capture the fractures and coefficient heterogeneity within each coarse element. The solving strategy of these local coarse-scale basis functions can be either geometric or algebraic [13,30,31]. We prefer algebraic construction, as it allows for blackbox integration of the method within any existing XFEM simulator. Once the basis functions are solved, they will be clustered in the matrix of Prolongation (P), which maps the coarse-scale solution to the fine-scale one. Note that there will be no additional multiscale basis functions due to jump or tips, and that only 4 multiscale basis functions per element exist for 2D structured grids (8 in 3D) in each direction (x, y, and z).
The fine-scale XFEM system is then mapped to the coarse grid by using the Restriction (R) operator, which is defined based on the FEM, as the transpose of the prolongation operator. The approximate fine-scale solution is finally obtained after mapping the coarse-scale solution to the fine scale, by using the prolongation operator.
The approximate solution of MS-XFEM can be found acceptable for many applications, however error control and reduction to any desired level is necessary to preserve its applicability for challenging cases. As such, the MS-XFEM is integrated within the two-stage iterative solver in which the MS-XFEM is paired with an efficient iterative smoother (here ILU(0)) to reduce the error [32,33]. One can also use the Krylov subspace methods (e.g. GMRES) to enhance the convergence, which stays outside the scope of this paper.
Several proof-of-concept numerical tests are presented to assess the accuracy of the presented MS-XFEM without and with iterative improvements. The test cases include large deformations which may not be realistic in geoscience applications, but important to be studied in order to quantify the errors in large deformation scenarios. From these results it becomes clear that the MS-XFEM, despite using no enriched basis functions at coarse scale, presents an efficient and accurate formulation to study deformation of fractured geological media.
The structure of this paper is set as the following. Next, the governing equations and the fine scale XFEM method are introduced. Then, the MS-XFEM method is presented in detail, with emphasis on the construction of local multiscale basis function and the approximate fine scale solution. Then, different numerical test cases are presented. Finally, concluding remarks are discussed.

Governing Equations and Fine-scale XFEM System
Consider the domain Ω bounded by Γ as shown in figure 1. Prescribed displacements or Dirichlet boundary condition are imposed on Γ u , while tractions are imposed on Γ t . The crack surface Γ c (lines in 2-D and surfaces in 3-D) is assumed to be traction-free.
where σ is the stress tensor and u is the displacement field over the whole domain. − → n is the normal vector pointing outside the domain [34,35].
The constitutive law with linear elasticity assumption reads where, ∇ s denotes the symmetrical operator and C is the property tensor defined as with λ and µ denoting the Lame's parameters [36,37].
The strain tensor ε is expressed as where, ∇ denotes the gradient operator.
Substituting Eqs. (5) and (6)  Equation (7) is then solved for computational domains with cracks (representing faults and fractures). This is done by the extended finite element (XFEM) method, which is briefly revisited in the next section.

Extended Finite Element Method (XFEM)
The FEM with smooth shape functions N i provides an approximate numerical solution to Eq. (7) for displacement unknowns, i.e., This formula can be used for computational domains without discontinuity. The FEM approximation is insufficient to capture discontinuities imposed by the existence of the fractures and faults. As such, the XFEM method introduces two sets of enrichment to the original FEM in order to allow it to capture the discontinuities without adapting the grid. These enrichment sets are associated with the body and tip of the fractures and faults. The body is enriched by jump functions, and the tip by tip enrichment functions [25]. Below brief descriptions of these two enrichment functions are provided.

Jump enrichment
The jump enrichment represents the discontinuity involved in the displacement field across the fracture and fault main body. The jump enrichment is often chosen as the step or Heaviside function, which can be expressed as Note that Ω + and Ω − zones are determined based on the normal vector pointing out of the fracture curve. For line fractures, the direction can be any side, as long as all discrete elements use the same + and -sides for a fracture.

Tip enrichment
The tip enrichment represents the discontinuity of the displacement field near the fracture tip. This type of enrichment function, denoted by F l , is based on the auxiliary displacement field near the fracture tip and contains four functions, i.e., These four functions around the fracture tip inside the element are plotted in Figure 2. The red segment, shown on the base of the plots, represents the fracture which ends in the element. Note that only the √ rsin( θ 2 ) contains a discontinuity around the fracture tip, while other functions are smooth.

Enrichment mechanism
To decide whether the node is enriched or not, the node location related to the fracture is the key factor. The sketch of the enrichment mechanism is shown in Figure 3. More precisely, in this figure, the tip and jump enriched nodes are highlighted in red and black, respectively. Figure 3: Enrichment mechanism: node I and J will be enriched using tip and jump functions.

XFEM linear system
The XFEM approximates the continuum displacement field u at fine-scale mesh resolution h by u h which is defined as The first term in the right-hand-side (RHS) of Eq. (10) is the contribution of the classical finite element method. This term captures the smooth deformation, using classical shape functions. The second term, however, represents the contribution of the jump enrichment. Note that the jump enrichment is modeled by the weighted Heaviside functions, with weights being the classical shape functions. There will be as many jump enrichment functions as the number of fractures inside an element. Finally, the third term in the RHS is the contribution of the fracture tips. Note that if several fracture tips end up in an element, there will be 4 additional DOFs per tip per direction in that element.
The resulting linear system entails the nodal displacement unknowns u, as well as the jump level a and tip weight b per fracture (and fault). The augmented Compared to the classical FEM, there exist several additional blocks involved in the stiffness matrix, due to the existence of the discontinuities. The advantage of XFEM is that it does not rely on complex mesh geometry, instead, it allows fractures to overlap with the matrix elements. On the other hand, for geoscientific fractured systems, the additional DOFs due to the enrichment procedure results in excessive computational costs. This imposes a significant challenge for the XFEM application in geoscience applications. In this paper, we develop a scalable multiscale procedure which constructs a coarse-scale system based on locally supported basis functions. The method is described in the next section.

Multiscale Extended Finite Element Method (MS-XFEM)
A multiscale formulation provides an approximate solution u h to the fine-scale XFEM deformation u h through where N H i are the coarse-scale (multiscale) basis functions and u H i are the coarse-scale nodal displacements at coarse mesh Ω H . Note that this multiscale formulation does not include any enrichment functions. Instead, all enrichment functions are incorporated in the construction of accurate coarse-scale basis functions N H . This allows for significant computational complexity reduction, and makes the entire formulation attractive for field-scale geoscientific applications.
Next, construction of the coarse-scale system and the basis functions will be presented.

Coarse scale linear system
MS-XFEM solves the linear deformation system on a coarse mesh, imposed on a given fine-scale mesh, as shown in figure 4. The coarsening ratio is defined as the ratio between the coarse mesh size and fine-scale mesh size. The multiscale formula (12) can be algebraically expressed as where P is the matrix of basis functions (i.e., prolongation operator) and d H is the coarse-scale deformation vector for u H unknowns. Algebraic formulation allows for convenient implementation of the proposed MS-XFEM, and its integration as a black-box tool for any given classical XFEM solver. Therefore, the remainder of the article will be devoted to the formulation.
The coarse-scale solution d H needs to be found by solving a coarse-scale system. To construct the coarse-scale system and solve for d H , one has to restrict (map) the fine-scale linear system(K h d h = f h ) to the coarse-scale, i.e., Here, R is the restriction operator with the size of Ω H × Ω h+j+t , where Ω h+j+t is the size of the fine-scale enriched XFEM system including jump and tip enrichment. Prolongation operator P has the dimension of Ω h+j+t × Ω H . This results in the coarse-scale system matrix K H size of Ω H × Ω H .
The finite-element-based restriction function is introduced as the transpose of the prolongation matrix, i.e., R = P T .
Therefore, the coarse-scale matrix K H is symmetric-positive-definite (SPD), if K h is SPD. Once the coarse-scale system is solved on Ω H space for d H , one can find the approximate fine-scale solution using Eq. (13). Overall, the multiscale procedure can be summarised as finding an approximate solution d h according to Next, the prolongation operator P, i.e., the basis functions are explained in detail. Once P is known, all terms in Eq. (16) are defined.

Construction of multiscale basis functions
To obtain the basis functions, the governing equation (7) without any source term using XFEM, i.e., (11) needs to be solved in each coarse element ω H .
This can be expressed as solving subject to local boundary conditions. Here, we develop a reduced-dimensional equilibrium equation to solve for the boundary cells [5,38], i.e., Here, Γ H denotes the boundary cells of the coarse element Ω H . In addition, ∇ S denotes the reduced dimensional divergence and symmetrical gradient operators, which act parallel to the direction of the local domain boundary. For 2D geometries, the reduced-dimensional boundary condition represents the 1D (rod) deformation model along the coarse element edges. Note that the local basis functions involve transverse equilibrium, i.e., therefore the prolongation matrix P reads  Note that, as shown in Fig. 5, the boundary problem is solved for both edges which have the node H at one of their end values. More precisely, e.g., to find the basis function in x-direction for node H, we set the value of u x (H) = 1 at the location H. This causes extension of the horizontal boundary cells and bending of the vertical boundary.
Once the boundary values are found, the internal cells are solved subjected to Dirichlet values for the boundary cells. An illustration of a basis function obtained using this algorithm is presented in figure 6. Note that the illustrated basis function captures the fractures, because of the XFEM enrichment procedure. The basis function N H i will be stored in the column i of the prolongation operator P. Once all basis functions are found, the operator P is also known and one can proceed with the multiscale procedure as explained before.
Next, we explain how the basis functions can be algebraically computed based on the given XFEM fine-scale system. This crucial step allows for convenient integration of our multiscale method into a given XFEM simulator.

Algebraic construction of multiscale basis functions
The basis function formulation (17) subjected to the local boundary condition (18) can be constructed and solved purely algebraically. This is important, since it allows for convenient integration of the devised multiscale method into any existing XFEM simulator.
Consider the coarse cell (local domain) as shown in figure 7. The cells are split into 3 categories of internal, edge and vertex (node), depending on their locations [13]. To develop the basis functions, first the fine-scale stiffness matrix K h is permuted, such that the terms for vertex, then edge and finally the internal cells appear. The permutation operator T as such reorders K h into K v such that Here, I represents the internal nodes, E represents the edge nodes and V represents the vertex nodes. The permuted linear system, therefore, reads Note that the permuted system collects all entries of the XFEM discrete system belonging to I, E, and V cells. Therefore, the XFEM enrichment entries due to tips and jumps are within their corresponding I, E, and V entries.
The reduced-dimensional boundary condition is now being imposed by replacing the 2D equation for E by a 1D XFEM discrete system. This leads the entrȳ K V I to vanish, as there will be no connectivity between the edge and internal cells for the edge cells. This 1D edge equations can then be expressed as Knowing that the solution at vertex cells will be obtained from the coarse-scale system, the reorder fine-scale system matrix can now be reduced to Note that the equations for basis functions do not have any source terms in their right-hand-side. The upper-triangular matrix of Eq. (23) can be easily inverted to give the prolongation operator, i.e., given the coarse nodes solutions d V , one can obtain the solution at the edge via similarly, the solution at the internal cells reads Note that P E and P I are the sub-matrices of the prolongation operator, i.e., Here, I V V is the diagonal identity matrix equal to the size of the vertex nodes.
After defining the prolongation operator algebraically, based on the entries of the 2D XFEM (for internal cells) and 1D XFEM (for edge cells), one can find the multiscale solution.

Iterative multiscale procedure (iMS-XFEM)
The multiscale solutions with the accurate XFEM basis functions can be used to provide an approximate efficient solution for many practical applications. However, it is important to control the error and reduce it to any desired tolerance [12] if needed. As such, the MS-XFEM is paired with a fine-scale smoother (here, ILU(0) [32]) to allow for error reduction. Note that this iterative procedure can also be used within a GMRES iterative loop [39] to enhance convergence rates. The study of the most efficient iterative strategy to reduce the error is outside the scope of the current paper. The iterative procedure then reads • Construct the P and R operators Smoothing stage (apply n s times ILU(0)): δd ν+1 = (M ns ILU(0) ) −1 r ν+1/2 Note that n s is defined by user.

Numerical Test Cases
In this section several test cases are considered to investigate the performance of MS-XFEM both as approximate solver and integrated within the iterative error reduction procedure.

Test case 1: Single fracture in a heterogeneous domain
In the first test, a square 2D domain of L × L with L = 10[m] is considered, which contains a single horizontal fracture in its centre, as shown in figure 8a.  A basis function for a fractured local domain is illustrated in figure 10. Note that the discontinuity is captured by the basis functions, since XFEM is used to solve for it. The effect of the coarsening ratio is shown in figure 11. The error e in this figure is computed using where N is the number of fine-scale mesh nodes. u i,M S and u i,f denote the proloned MS-XFEM solution displacement field and fine-scale solution displacement field, respectively. The MS-XFEM errors are due to the local boundary conditions used to calculate the basis functions, and also because no additional enrichment functions are imposed at coarse scale. Note this means that for heterogeneous domains there can be a finer resolution for coarse cells at which the local boundary conditions impose more errors compared with coarser resolutions. In spite of this, figure 11 clearly shows, for this example, a decaying trend of the error with respect to the finer coarse mesh is observed.

Iterative MS-XFEM
As discussed in section 3.4, one can pair the MS-XFEM in an iterative strategy in which the error is reduced to any desired level [13]. From figure 9 (c) that with 3 times of fine scale smoothers applied in the second stage after 3 iterations the MS-XFEM result has been improved compared to the fine-scale result and the error is decreased significantly. Results of the iterative MS-XFEM procedure (iMS-XFEM) are shown in figure 12. Different smoothing steps per iteration values n s are used. Note than neither GMRES [39] nor any other iterative convergence enhancing procedure is used here. Clearly, one can reduce the multiscale errors to the machine accuracy by applying iMS-XFEM iterations. In particular, for practical applications, one can stop iterations after a few counts, once the error norm is below the level of uncertainty τ within the parameters of the problem.
In which τ is chosen as 10 −10 in here. In figure 12 it is shown that convergence is achieved with n s rounds of the fine scale smoother involved in the second stage.

Test case 2: heterogeneous reservoir with multiple fractures
The second test case is set to model deformation in a heterogeneous reservoir with more fractures. The size and the heterogeneous properties of this test case are the same as those in test case 1. Here, more fractures are considered. In addition, compared to test case 1, the east and west boundaries are also observing distributed loads, as shown in figure 13.  An example of two basis functions for this test case is shown in figure 15. In the local plot, it is illustrated how 1 (15a) and 2 (15b) fractures are captured by the basis functions. The iMS-XFEM procedure, as explained before, is now applied to reduce the multiscale errors to machine precision. Still in figure 14 (c) the result quality is improved a lot after 3 iterations with 3 times fine-scale smoothers applied in the second stage. Note that since no GMRES nor a complete smoother is used, but the incomplete smoother ILU(0) for its efficiency. Results are shown in figure 16.

Conclusion
A multiscale procedure for XFEM is proposed to model deformation of geological heterogeneous fractured fields. The method resolves the discontinuities through local multiscale basis functions, which are computed using XFEM subjected to local boundary conditions. The coarse-scale system is obtained by using the basis functions, algebraically, which does not have any additional enrichment functions, in contrast to the local basis function systems. This procedure makes the MS-XFEM very efficient. Also, by combining it with a fine-scale smoother, an iterative MS-XFEM (iMS-XFEM) procedure is developed, which allows to reduce the error to any desired level of accuracy.
Two heterogeneous test cases were studied as proof-of-concept, to investigate the performance of the MS-XFEM. It was shown that MS-XFEM results in acceptable solutions, when no iterations are imposed. By applying iterations, one can further improve the results. For practical applications, when parameters are uncertain, only a few iterations can be applied to maintain (and control) the MS-XFEM quality of the solution.