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 kilometers.Numerical simulation of mechanical de-formation for such complex systems is necessary to optimize the geo-engineering operations [1,2], and assess their safety and manage the associated risks (e.g.fracture propagation, fault slip and induced seismicity).However, classical numerical schemes become too computationally costly to be applied to these systems.It is therefore necessary to develop scalable simulation methods, which preserve the heterogeneity and complex deformation physics of a fractured medium without compromising the accuracy of the simulations.
The presence of fractures within computational domains can be treated explicitly by two approaches of (1) unstructured grid and (2) immersed or embedded methods.The unstructured grid approach [3][4][5][6][7] 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 mesh geometries.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) [8][9][10][11][12][13][14].These enriched formulations are aimed at representing discontinuities within the overlapping matrix cell, without any adjustment nor refinement of the grid [15].The enrichment strategies for modeling deformation of fractured media, using finite-element schemes, are referred to as 'extended finite element (XFEM)' methods.
XFEM enriches the partition of unity (PoU) [16] by introducing additional degrees of freedom (DOFs) 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) [17][18][19][20][21][22].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.
Presence of highly heterogeneous coefficients with high resolution within large-scale domains has been systematically addressed in the computational geoscience community through the development of multiscale finite element and finite volume methods [23][24][25][26].Recent developments also include mechanical deformation coupled with fluid pore pressure dynamics [27][28][29][30][31], in which local multiscale basis functions were computed using either linear or reduced dimensional boundary conditions.In presence of fractures, specially many fractures, however, 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 [32][33][34], is of high interest in the geoscience community.
When it comes to geoscience applications, the classical XFEM is not an attractive method, due to its excessive additional degrees of freedom to capture the many fractures.
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 enriched 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 [33,45,46].We prefer algebraic construction, as it allows for black-box integration of the method within any existing XFEM simulator.Once the basis functions are solved, they will be clustered in the prolongation matrix (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) matrix, which is defined based on the FEM, as the transpose of P. The approximate fine-scale solution is finally obtained after mapping the coarse-scale solution to the fine scale grid, by using the prolongation operator.
The approximate solution of MS-XFEM is 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 [47,48].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 Fig. 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.
The momentum balance equations and boundary conditions read where σ is the stress tensor and u is the displacement field over the whole domain.− → n is the normal vector pointing outside the domain [49,50].
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 [5,51].
The strain tensor ε is expressed as where, ∇ denotes the gradient operator.
Substituting Eqs. ( 5) and ( 6) in the governing Eq. ( 1) results in a 2nd order Partial Differential Equation (PDE) for displacement field u Eq. ( 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)
FEM employs smooth shape functions N i for each node i ∈ I , where I is the set of all nodes, to provide an approximate numerical solution to Eq. (7) for displacement unknowns, i.e., u = i∈I u i N i .This smooth 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 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 [17].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.,

)}.
Here, θ ranges from [−π , π], which is defined based on the local polar coordinates of a point with respect to the fracture tip.An illustration is given in Fig. 2.

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 Fig. 3.More precisely, in this figure, the tip and jump enriched nodes are highlighted in circular point and square point, respectively.

XFEM linear system
The XFEM approximates the displacement field u at fine-scale resolution h by u h which is defined as where N, H and F l represent, respectively, the classical FEM shape functions, the Heaviside function and the tip enrichment functions.The fine-scale mesh has h nodes.Moreover, u denotes the standard DOFs associated to the classical finite element method.The first term in the right-hand-side (RHS) of Eq. ( 9) 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 XFEM linear system 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 Fig. 4. The coarsening ratio is defined as the ratio between the coarse mesh size and fine-scale mesh size.
The multiscale formula (11) can be algebraically expressed as where P is the matrix of basis functions (i.e., prolongation operator) and d H is the coarse-scale algebraic deformation vector, corresponding to the physical node-based 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., (R K h P) 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.(12).Overall, the multiscale procedure can be summarized 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. ( 15) are defined.

Construction of multiscale basis functions
To obtain the basis functions, the governing Eq. ( 7) without any source term using XFEM, i.e., Eq. ( 10) 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 [25,52], i.e., Here, H denotes the boundary cells of the coarse element H .In addition, ∇ S denotes the reduced dimensional divergence Fig. 5 shows an example of a local system to be solved for basis functions belonging to the highlighted node H in x and y directions.Note that the Dirichlet value of 1 is set at H for each directional basis functions, while all other 3 coarse mesh nodes are set to 0. Furthermore, 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.
The sub-block matrices P xx and P xy refer to the deformation along x-axis and y-axis, respectively, in response to the unit horizontal displacement at the coarse node H , as shown in Fig. 5(a).Similarly, P yy and P yx refer to the deformation along y-axis and x-axis, respectively, in response to the unit vertical displacement at the coarse node H , as shown in Fig. 5(a).
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 Fig. 6.It is worth to be emphasized that the local boundary, when crossed by fractures, will be enriched by jump enrichment functions.Specially, there exists a chance that fracture tips end exactly at local boundaries.In that case, since the reduced-dimensional boundary condition solves 1D problems along the boundary, the tip enrichment is replaced by the jump enrichment, and consequently the original fine-scale XFEM stiffness matrix is also modified in the same manner.
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 (16) subjected to the local boundary condition ( 17) 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 Fig. 7.The cells are split into 3 categories of internal, edge and vertex (node), depending on their locations [33].
Note that the vertex nodes are in fact the coarse mesh nodes, where the coarse-scale solution will be computed.The basis functions are needed to interpolate the solution between the vertex cells through the edge and internal cells.
To develop the basis functions, first the fine-scale stiffness matrix K h is permuted, such that the terms for internal, then edge, and finally the vertex 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 entry K E 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 consider any source terms in their right-hand-side, and that we have replaced the equations for the vertex cells by the coarse-scale system (13).
Once the coarse-scale solutions d V are found, the upper-triangular matrix of Eq. ( 22) can be conveniently inverted to obtain the prolongation operator.More precisely, 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 [32] if needed.As such, the MS-XFEM is paired with a fine-scale smoother (here, ILU(0) [47]) to allow for error reduction.Note that this iterative procedure can also be used within a GMRES iterative loop [53] 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 reads Smoothing stage: iterate from j = 1 till j = n s times [2a] Apply ILU(0): δd ν+1/2+1/2( j/n s ) = M −1 ILU(0) r ν+1/2+1/2( j−1)/n s [2b] Update residual r ν+1/2+1/2( j/n s ) Here, the approximate ILU(0) smoother operator M ILU(0) is applied for n s times on the updated residual vector.

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 center, as shown in Fig. 8(a).The fine-scale mesh consists of 40 × 40 cells, while the MS-XFEM contains only 5 × 5 coarse grids.This results in a coarsening ratio of 8 ×8.The heterogeneous Young's modulus distribution is shown in Fig. 8(b), while the Possion's ratio is assumed to be constant 0.2 everywhere.The fracture tip coordinates are shown in Fig. 8(a).The Dirichlet boundary condition is set at the south face, while the north boundary is under distributed upward load with q = 5 × 10 5 [N/m] magnitude.The east and west boundaries are set as stress-free.
Results are shown in Fig. 9.It is clear that the results of MS-XFEM on only 5 grid cells is in reasonable agreement with that of the fine-scale XFEM solver using a 40 × 40 mesh.Note that no enrichment for the MS-XFEM is used, and the basis functions are computed using the XFEM method on local domains.
The error at each node i is defined as the difference between fine-scale XFEM (i.e., u i, f ) and MS-XFEM (i.e., u i,M S ) solutions, i.e., The stress field σ yy , obtained from the displacement field, is also plotted in Fig. 11      A basis function for a fractured local domain is illustrated in Fig. 12.Note that the discontinuity is captured by the basis functions, since XFEM is used to solve for it.
Coarsening ratio indicates the number of fine-cells in each coarse cell.The effect of the coarsening ratio is shown in Fig. 13.To do this study, the fine-scale mesh is kept constant, while the coarsening ratio is changed.This results in different coarse-scale system sizes to solve for the same fine-scale system.The normalized 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 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, Fig. 13 shows, for this example, a decaying trend of the error with respect to the finer coarse mesh is observed.
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 [33].From Fig. 14 that with 5 times of fine scale smoothers applied in the second stage after 5 iterations the multiscale solution error is decreased significantly.Also as shown in Fig. 11(a) and (c) after 5 iterations the stress fields of iterative MS-XFEM (iMS-XFEM) and fine-scale XFEM match well.Convergence history of the iMS-XFEM is shown in Fig. 15.The blue circle represents the normalized error of the result without the iterative strategy, while the red circle represents the normalized error after 5 iterations within each 5 smoothing operations are applied.These 2 circles correspond to the results shown in Fig. 14 and Fig. 10.Different smoothing steps per iteration values n s can be used.Note that neither GMRES [53] 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   Here, τ is chosen as 10 −10 to confirm convergence can be achieved to machine accuracy.Fig. 15 shows the convergence history for different smoothing step counts n s per iteration.

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 Fig. 16.
Simulation results for both fine-scale XFEM and MS-XFEM are shown in Fig. 17.Multiscale solutions are obtained by 8 × 8 coarsening ratio, applied to the fine-scale mesh with 40 × 40 fine cells.This leads to significant cost reduction, as the multiscale coarse system has only 5 × 5 cells with no additional degrees of freedom due to enrichments.The black lines on Fig. 17(b) represent the coarse scale mesh.It is clear that the MS-XFEM (without any iterations) result is significantly different than the fine-scale one in the zones of large deformations.Stress fields for all three cases of fine-scale XFEM, MS-XFEM and iMS-XFEM are shown in Fig. 18.The iMS-XFEM results are shown after 10 iterations with n s = 5.Note that after iterations, multiscale stress distributions agree with the fine-scale one.Additionally, one can employ adaptive coarsening ratios, so to minimize the need to use an iterative strategy [54].ence around other three fractures have been improved significantly.Note that no GMRES nor an unconditionally-convergent smoother is used, but the ILU(0) is used as the second stage smoother.ILU(0) provides an approximate (incomplete) decomposition of the fine scale system for maintaining the efficiency [53].As shown in Fig. 21, a higher number of smoothing steps would result in faster convergence.However, it comes with additional computational costs.Similar as the extensive studies for multiscale simulation of flow in porous media (see e.g.[55][56][57]), for deformation simulation, a CPU-based study should be performed to obtain the optimum combination of coarse-grid resolution, type and count of smoothing steps.

Computational efficiency analysis
Appropriate CPU-based speedup analyses is out of the scope of this paper.However, to provide an estimate of the computational efficiency of the proposed MS-XFEM strategy, compared with a fine-scale XFEM, an operator-cost-based analysis is presented below.Then, a scalability test with respect to the performance of error reduction strategy with increasing density of the fractures is also studied and presented.

Approximate speedup based on an operation-based analysis
The approximate operation-based approach builds on the assumption that solving a linear system of size ζ costs O (ζ m ) operations, i.e., C = αζ m , where C is the cost function.As such, the approximate cost function for fine-scale XFEM system (i.e., C f ) with N f nodes is Here, ψ represents the number of extra degrees of freedom due to XFEM enrichments.also since there are x and y displacement directions for each node, a factor 2 is used.
MS-XFEM (with and without iterations) provides an approximate solution to the fine-scale system in a procedure that includes three main steps: (1) calculate basis function, (2) solve coarse scale system, and (3) iteratively improve the results.Consider the coarsening ratio of , thus, the coarse girds size N c is  Assuming that γ enrichments per each basis function is considered, basis functions for each coarse node of a coarse cell costs where the factor 2 represents displacements in x and y directions.Note all N c coarse nodes will have basis functions.The coarse scale system entails no extra degree of freedom due to enrichments.As such, its cost C c can be approximated as Again the factor 2 is for the 2D (x, y) unknowns.The cost of the linear smoother I LU (0) is approximated by setting m = 1 and a constant is β, i.e., where the n s is the number of smoothing steps per iMS-XFEM iterations.
If no iterative strategy is applied, then, the MS-XFEM computational cost C M S can be approximated as where If the iterative strategy is applied for N it times, then the iMS-XFEM computational cost C iM S can be approximated as Note that if N it = 0, the cost function C iM S becomes identical with C M S .The speedup factor φ is now introduced as the cost ratios between iMS-XFEM and fine-scale XFEM, i.e., φ = This function results in the speedup for MS-XFEM (no iterations) when N it = 0. To provide a better sense about the speedup, two tests with different fine scale grid numbers of N f = 10 4 and N f = 10 8 are considered.For the case with N f = 10 4 , coarsening ratios of ∈ {2 × 2, 5 × 5, 10 × 10, 20 × 20, 50 × 50} are considered, while for the case with N f = 10 8 , ∈ {10 × 10, 20 × 20, 50 × 50, 100 × 100, 1000 × 1000] are employed.It is assumed that α/β = 1, meaning the scalability factor of the smoother is close to a linear solver.In addition, we set m = 1.3, as the exponent of linear solver computational complexity.The same factor was used in the literature [58].Furthermore, extra unknowns due to XFEM enrichments are considered to be 10% of fine-scale grid cells, i.e., ψ = 0.1N f .The basis functions are then considered to have these extra DOFs equally distributed, which leads to γ ≈ ψ N c . With these settings, one can study the overall approximate speedup of the proposed MS-XFEM, based on an operational-based analyses.For that MS-XFEM (no iterations), iMS-XFEM with N it = 2 and N it = 5 are considered.For both iMS-XFEM cases, n s = 8.The speedup factor φ for all tests is plotted in Fig. 22.
From Fig. 22, it is clear that the MS-XFEM provides a significant speedup compared with fine-scale XFEM.The speedup becomes bigger when the problem size is bigger.Also, from the presented operational-based analyses, one can consider an  optimum coarsening ratio, in which the multiscale procedure would perform optimal.Note that, similar as for multiphase flow simulations [57,59], a CPU-based analyses is needed to draw more conclusive remarks about the real speedup of the MS-XFEM.

Scalability test
In this test the effect of the density of fractures on convergence rate of the iMS-XFEM is investigated.This is specially relevant to geoscientific applications.A fixed domain size of 10 × 10 m 2 is considered.As it can be seen in Fig. 24, the convergence rate of the developed iMS-XFEM is not much affected by the increasing fracture density.This proof-of-the concept analysis indicates the applicability of the devised method for real-field applications.

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.Through an iterative procedure, only the coarse-scale system and smoothing steps are repeated; while basis functions are not updated.
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.The operational-based computational speedup indicated that the developed multiscale procedure is promising for real-field applications.Further research is needed to test the method for complex systems with 3D heterogeneous fractured media, specially coupled with flow and heat transfer.Additionally, to minimize the need for an iterative strategy and adaptive coarse grid can be employed [46,60,61].We expect the speedup increase when time-dependent problems are solved, as the basis functions would need to be infrequently updated, only when local properties change.Nevertheless, the real speedup will be found when the method is implemented in compilable language.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 3 .
Fig.3.Enrichment mechanism: node I and J will be enriched using tip and jump functions.
Furthermore, a denotes the extra DOFs associated to the jump enriched node.For the jump enriched nodes, in the 2D domains, each node would contain 2 extra DOFs.Similarly, b indicates the extra DOFs associated to the tip enrichment, which adds four extra DOFs per direction (8 in total in a 2D domain) for each tip inside an element.

Fig. 4 .
Fig.4.Illustration of the multiscale mesh imposed on the given fine-scale mesh, with a coarsening ratio of 3 × 3.
and symmetrical gradient operators, which act parallel to the direction of the local domain boundary.Moreover, C r is the reduced-dimensional (here, 1D) average elasticity tensor along the boundary of the coarse cell.More precisely, for boundaries parallel to x-direction, the only nonzero entries are C x r (1, 1) = λ and C x r (3, 3) = μ; while its nonzero entries for the boundaries parallel to y-direction are C y r (2, 2) = λ and C y r (3, 3) = μ.Here, the λ and μ are averaged elasticity properties of the adjacent 2D cells belonging to the boundary 1D edges.For 2D geometries, the reduced-dimensional boundary condition represents the 1D (rod) deformation model along the coarse element edges.The local deformation in response to horizontal and vertical unity nodal displacements will be in general two-dimensional.Therefore, the prolongation matrix P will have off-block-diagonal entries (i.e., P xy = 0 and P yx = 0) and reads

Fig. 5 .
Fig. 5. Illustration of the multiscale local basis functions, constructed using XFEM for the node H in x direction (a) and y direction (b).

Fig. 6 .
Fig. 6.Illustration of a basis function that captures the discontinuity of a fracture.Model size is 10 × 10 m 2 .Yellow segment represents the fracture extending from (3, 3) to (7, 3) coordinates.The shown basis function belongs to the coarse node located at the coordinate (4, 4).(For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Illustration of the 3 categories of Internal, Edge, and Vertex cells, corresponding to the position of each fine cell within a coarse element.

Fig. 10
Fig.10illustrates the error distribution for Test Case 1.Note that the zones of high errors are mainly near the fracture.

Fig. 8 .
Fig. 8. Test case 1: (a) illustration of the model setup, (b) heterogeneous Young's modulus distribution.Note the units are SI.The black lines represent the coarse mesh.

Fig. 9 .
Fig. 9. Test Case 1: displacement field for a heterogeneous fractured reservoir using (a) fine scale XFEM and (b) MS-XFEM.Black lines represent the coarse scale mesh.

Fig. 10 .
Fig. 10.Test Case 1: displacement difference between fine scale XFEM and MS-XFEM in (a) x direction (b) y direction.Black lines represent coarse scale mesh.

Fig. 12 .
Fig. 12. Test case: basis functions of single fracture test case.Single discontinuity is captured by axial equilibrium and transverse equilibrium solutions.This basis function is centered at coordinate (4, 6).

Fig. 14 .
Fig. 14.Test Case 1: displacement difference between fine scale XFEM and MS-XFEM in (a) x direction (b) y direction after 5 iterations, within which 5 smoothing steps were employed.Solid black lines represent the coarse mesh.a few iteration counts, once the error norm falls below a certain tolerance τ , considering the level of uncertainty within the parameters of the problem, i.e., e i τ , i = x, y.

Fig. 15 .
Fig. 15.Test case 1: iteration history of iMS-XFEM procedure with different number of smoothing per step.Errors for displacement in x (a) and y (b) directions are shown.The blue and the red circles represent the normalized errors corresponding to the results shown in Fig. 10 and Fig. 14, respectively.

Fig. 16 .
Fig. 16.Test case 2: multiple fractures within a heterogeneous reservoir under tension stress across three boundaries.

Fig. 17 .
Fig. 17.Test Case 2: displacement field for a heterogeneous media with multiple fractures for (a) fine scale XFEM and (b) MS-XFEM without iterative improvements.Note the 40 × 40 fine-scale XFEM system is solved by applying only 5 × 5 coarse cells, using MS-XFEM.Black lines represent the coarse scale mesh lines.An example of two basis functions for this test case is shown in Fig. 19.It illustrates, in 2 different views, how two fractures are captured by the local basis functions.The iMS-XFEM procedure, as explained before, is now applied to reduce the multiscale errors to machine precision.Fig.20(c) and (d) present improved results compared to 20(a) and (b), after 10 iterations with 5 times fine-scale smoother applied in each iteration.The large deformation of the most left fracture still requires more iterations however the differ-

Fig. 18 .
Fig. 18.Test Case 2: stress field σ yy comparison between (a) fine scale XFEM and (b) MS-XFEM (c) MS-XFEM after 10 iterations with 5 rounds of fine scale smoother applied in the second stage.

Fig. 19 .
Fig. 19.Illustration of the basis functions P yy where two discontinuities are captured in a local coarse cell.

Fig. 20 .
Fig. 20.Test Case 2: displacement difference between fine scale XFEM and MS-XFEM (a) x direction (b) y direction (c) x direction after 10 iterations (d) y direction after 10 iterations.5 rounds of fine scale smoother are used in the second stage.Black lines represent the coarse mesh.

Fig. 21 .
Fig. 21.Iteration history of iMS-XFEM procedure with different number of smoothing per step.Errors for displacement in x (a) and y (b) direction are shown.N c = N f .

Fig. 22 .
Fig. 22. Speedup factor for different coarsening ratios and different iteration numbers for fine-scale gird of N f = 10 8 (a) and N f = 10 4 (b).The extra DOFs due to XFEM enrichment is ψ = 0.1N f is 10 7 .For both cases the number of smoothing steps per iteration is set to n s = 8.

Fig. 24 .
Fig. 24.Iteration history of iMS-XFEM procedure with different numbers of extra DOFs per step for displacement in x (a) and y (b) directions.F1 to F7 phases represent the increasing fracture density as show in Fig. 23.

Fine and coarse grids are 40 × 40 and 5 ,
respectively.As shown in Fig.23, starting from phase F 1 to phase F 7 one fracture is added in the domain.For convenience, homogeneous Young's modulus of E = 10 7 Pa and the Possion's ratio of 0.2 are considered.The boundary conditions are the same as in test case 2.