A multiscale flux basis for mortar mixed discretizations of reduced Darcy-Forchheimer fracture models

In this paper, a multiscale flux basis algorithm is developed to efficiently solve a flow problem in fractured porous media. Here, we take into account a mixed-dimensional setting of the discrete fracture matrix model, where the fracture network is represented as lower-dimensional object. We assume the linear Darcy model in the rock matrix and the non-linear Forchheimer model in the fractures. In our formulation, we are able to reformulate the matrix-fracture problem to only the fracture network problem and, therefore, significantly reduce the computational cost. The resulting problem is then a non-linear interface problem that can be solved using a fixed-point or Newton-Krylov methods, which in each iteration require several solves of Robin problems in the surrounding rock matrices. To achieve this, the flux exchange (a linear Robin-to-Neumann co-dimensional mapping) between the porous medium and the fracture network is done offline by pre-computing a multiscale flux basis that consists of the flux response from each degree of freedom on the fracture network. This delivers a conserve for the basis that handles the solutions in the rock matrices for each degree of freedom in the fractures pressure space. Then, any Robin sub-domain problems are replaced by linear combinations of the multiscale flux basis during the interface iteration. The proposed approach is, thus, agnostic to the physical model in the fracture network. Numerical experiments demonstrate the computational gains of pre-computing the flux exchange between the porous medium and the fracture network against standard non-linear domain decomposition approaches.


Introduction
Using the techniques of domain decomposition [18], a first reduced model has been proposed for flow in a porous medium with a fracture in which the flow in the fracture is governed by the Darcy-Forchheimer's law while that in the surrounding matrix is governed by Darcy's law.
We consider here the generalized model given in [28], for which we let Ω to be a bounded domain in R d , d = 2, 3, with boundary Γ = ∂Ω, and we let γ ⊂ Ω be a (d − 1)-dimensional surface that divide Ω into two sub-domains: Ω = Ω 1 ∪ Ω 2 ∪ γ, where γ = ∂Ω 1 ∩ ∂Ω 2 and Γ i = ∂Ω i ∩ ∂Ω, i = 1, 2. The reduced model problem as presented in [28] is as follows: 1a) for i = 1, 2, together with ∇ τ · u γ = f γ + (u 1 · n 1 + u 2 · n 2 ) in γ, (1.2b) p γ = 0 in ∂γ, (1.2c) and subject to the following interface conditions −u i · n i + α γ p i = α γ p γ on γ, (1.3) for i = 1, 2. Here, ∇ τ denotes the (d − 1)-dimensional gradient operator in the plane of γ, the coefficient K i , i = 1, 2, is the hydraulic conductivity tensor in the sub-domain Ω i , and K γ is the hydraulic conductivity tensor in the fracture, I ∈ R d×d is the identity matrix, n i is the outward unit normal vector to ∂Ω i , and β γ is a non-negative scalar known as the Forchheimer coefficient. In (1.3), the coefficient α γ is a function proportional to the normal component of the permeability of the physical fracture and inversely proportional to the fracture width/aperture. We refer to [19] for a more detailed model description. For illustration purposes, we give a simple graphical example of a fractured porous medium in Figure 1. 3) can be seen as a Robin boundary condition for the sub-domain Ω i with a dependence on the pressure on the fracture γ. Clearly, if β γ = 0, then (1.2) is reduced to a linear Darcy flow in the fracture. The homogeneous Dirichlet boundary conditions (1.1c) and (1.2c) are considered merely for simplicity. The functions f i ∈ L 2 (Ω i ), i = 1, 2 and f γ ∈ L 2 (γ) are source terms in the matrix and in the fracture, respectively.
The mixed-dimensional problem (1.1)-(1.3) is an alternative to the possibility of using a very fine grid in the physical fracture and a necessarily much coarser grid away from the fracture. This idea was developed in [2] for highly permeable fractures and in [5] for fractures that may be highly permeable or nearly impermeable. We also refer to [30,34,14] for similar models. For all of the above models, where the linear Darcy's law is used as the constitutive law for flow in the fractures as well as in the surrounding domains, there are interactions between fractures and surrounding domains. This coupling is ensured using Robin type conditions as in [33], delivering discontinuous normal velocity and pressure across the fractures. Particularly, for fractures with large enough permeability, Darcy's law is replaced by Darcy-Forchheimer's law as established in [28], which complicates the coupling with the surrounding medium.
Several numerical schemes have been developed for fracture models, such as a cell-centered finite volume scheme in [24], an extended finite element method in [10], a mimetic finite difference [6] and a block-centered finite difference method in [31]. The aforementioned numerical approaches solve coupled fracture models directly. However, different equations defined in different regions are varied in type, such as coupling linear and non-linear systems, and often interface conditions involve new variables in different domains, which results in very complex algebraic structures. Particularly, several papers deal with the analysis and implementation of mixed methods applied to the above model problem in the linear case, on conforming and non-conforming grids [33,9,17,38,39]. In [18], the model problem (1.1)-(1.3) was solved using domain decomposition techniques based on mixed finite element methods (see [2] for the linear counterpart).
The purpose of this paper is to propose an efficient domain decomposition method to solve (1.1)-(1.3) based on the multiscale mortar mixed finite element method (MMMFEM) [22]. The method reformulates (1.1)-(1.3) into an interface problem by eliminating the sub-domain variables. The resulting interface problem is a superposition of a non-linear operator handling the flow on the fracture and a linear operator presenting the flux contribution from the sub-domains. When applying the MMMFEM, an outer-inner iterative algorithm like, the Newton-GMRes (or any Krylov solver) method or fixed-point-GMRes method, is used to solve the interface problem. As an example, if a fixed-point method (outer) is adopted, the linearized interface equation for the interface update can be solved with a domain decomposition algorithm (inner), in which at each iteration sub-domain solves, together with inter-processor communication, are required. The main issue of this outer-inner algorithm is that it leads to an excessive calculation from the sub-domains, as the dominant computational cost is measured by the number of sub-domain solves.
The new implementation recasts this algorithm by distinguishing the linear and non-linear contributions in the overall calculation and employing the multiscale flux basis functions from [23] for the linear part of the problem, before the non-linear interface iterations begin. The fact that the non-linearity in (1.1)-(1.3) is only within the fracture, we can adopt the notion that sub-domain problems can be expressed as a superposition of multiscale basis functions. In our terminology the mortar variable considered in [22] becomes the fracture pressure, these multiscale flux basis with respect to the fracture pressure can be computed by solving a fixed number of Robin sub-domain problems, that is equal to the number of fracture pressure degrees of freedom per sub-domain. Furthermore, this is done in parallel without any inter-processor communication.
An inexpensive linear combination of the multiscale flux basis functions then circumvents the need to solve any sub-domain problems in the inner domain decomposition iterations. This procedure can be enhanced by applying interface preconditioners as in [4,33,25] and by using a posteriori error estimates of [36] to adaptively refine the mesh grids. This calculation made in an offline step typically spares numerous unnecessary sub-domain solves. Precisely, in the original implementations, the number of sub-domain solves is approximately equal to where N lin is the number of iterations of the linearization procedure, and N k dd denotes the number of domain decomposition iterations (GMRes or any Krylov solver). For the new implementation, the number of sub-domains solves will be reduced if N lin k=1 N k dd exceeds the maximum number of fracture pressure degrees of freedom on any sub-domain.
This step of freezing the contributions on the flow from the rock matrices can be easily coded, cheaply evaluated, and efficiently used in practical simulations, i.e, it permits reusing the same basis functions to extend (1.1)-(1.3) to simulate various linear and non-linear models for flow in the fracture, such as generalized Forchheimer's laws: as well as exploring the fracture and barrier cases and comparing in a cheap way various non-linear solvers to (1.1)-(1.3). Crucially, the present approach can naturally be integrated into discrete fracture networks (DFNs) models [38,39,20,16], which in contrast to discrete fracture models (DFMs), do not consider the flow in the surrounding sub-domains, but handle both a large number of fractures and a complex interconnecting network of these fractures. For the presenting setting, we allow for the discretization of (1.1)-(1.3) by different numerical methods applied separately in the surrounding sub-domains and in the fracture. We allow for the cases where the grids of the porous sub-domains do not match along the fracture, where different mortar grid elements are used. We also investigate the case where the permeability in the fracture K γ is much lower than the permeability in the surrounding matrix K.
The library PorePy [27] has been used and extended to cover the numerical schemes and examples introduced in this article. The main contribution to the library is the implementation of the multiscale and domain decomposition frameworks. Even if we focus on lowest-order Raviart-Thomas-Nédélec finite elements, our implementation is agnostic with respect to the numerical scheme. The example presented are also available in the GitHub repository.
This paper is organized as follows: Firstly, the variational formulation of the problem and the MMM-FEM approximation are given in Section 2. Therefrom, the reduction of the original problem into nonlinear interface problem is introduced. The linearization-domain-decomposition procedures are formulated in Section 3. Section 4 describes the implementation based on the multiscale flux basis. We show that structurally the same implementation can be extended for more complex intersecting fractures model. Finally, we showcase the performance of our method on several numerical examples in Section 5 and draw the conclusions in Section 6.
2 Non-linear domain decomposition method As explained earlier, it is natural to solve the mixed-dimensional problem (1.1)-(1.3) using domain decomposition techniques. To this aim, we introduce the weak spaces in each sub-domain Ω i , i = 1, 2, and define their global versions by Equivalently, we introduce the weak spaces on the fracture γ, i.e, Following [18,33], a mixed-dimensional weak form of (1.1)-(1.3) asks for (u, p) ∈ V × M and (u γ , p γ ) ∈ V γ × M γ such that, for each i ∈ {1, 2}, where we introduced the functions K and f in Ω 1 ∪ Ω 2 such that K i = K| Ωi , and f i = f | Ωi , i = 1, 2. The jump · is defined by u · n := u 1 · n 1 + u 2 · n 2 , with n i the outer unit normal vector of Ω i on γ, for i = 1, 2. Finally, the non-linear term is defined as The reader is referred to [28] for proof of the existence and uniqueness of a solution to the variational formulation (2.1).

The discrete problem
Let T h,i be a partition of the sub-domain Ω i into either d-dimensional simplicial and/or rectangular elements. We also let T h,γ to be a partition of the fracture γ into (d − 1)-dimensional simplicial and/or rectangular elements. Note that, for both partitions, general elements can be treated via sub-meshes, see [37] and the references therein. Moreover, we assume that each partition is conforming within each sub-domain as well as in the fracture. The meshes T h,i , i = 1, 2, are allowed to be non-conforming on the fracture-interface γ, but also different from T h,γ . We then set T h = ∪ 2 i=1 T h,i and denote by h the maximal element diameter in T h . For the scalar unknowns, we introduce the approximation spaces M h = M h,1 × M h,2 and M h,γ , where M h,i , i = 1, 2, respectively M h,γ , is the space of piecewise constant functions associated with T h,i , i = 1, 2, respectively T h,γ . For the vector unknowns, we introduce the approximation spaces V h = V h,1 × V h,2 and V h,γ , where V h,i , i = 1, 2 and V h,γ , are the lowest-order Raviart-Thomas-Nédélec finite elements spaces associated with T h,i , i = 1, 2 and T h,γ , respectively. Clearly, in contrast to what is done in [22,12], we use the same order of the polynomials for the interface-pressure and the normal traces of the sub-domain velocities on the interface. The discrete mixed-dimensional finite element approximation of (2.1) is as follows: The next step in formulating a multiscale flux basis algorithm to solve (2.2) is to adopt domain decomposition techniques to reduce the global mixed-dimensional problem to an interface problem posed only on the fracture [3].

Reduction to interface problem
We introduce the discrete (linear) Robin-to-Neumann operator S RtN i , i = 1, 2: is the solution of the sub-domain problems with source term f , homogeneous Dirichlet boundary condition on ∂Ω, and λ as a Robin boundary condition along the fracture γ, i.e, for i = 1, 2, Then we set With these notations, we can see that solving (2.2) is equivalent to solving the following non-linear mixed interface problem: or equivalently The above distinction is classical in domain decomposition techniques in which we split the sub-domain problems into two families of local problems on each Ω i : one is with zero source and specified Robin value on the fracture-interface, and the other is with zero Robin value on the fracture-interface and specified source. In compact form, the mixed interface Darcy-Forchheimer problem (2.5) can be rewritten as This system is a non-linear mixed interface problem [13] that can be solved iteratively by using fixed point iterations or via a Newton-Krylov method. To present the two approaches, let us first consider the linear context, i.e, suppose the operator K −1 (·) is linear. Then (2.7) is the system associated to the linear mixed Darcy problem on the fracture that can be solved using a Krylov type method, such as GMRes or MINRes. Given an initial guess w as an approximate solution to (2.7), where A γ is the associated stiffness matrix of the linear system, b γ is the right-hand side, and Clearly, each GMRes iteration needs to evaluate the action of the Robin-to-Neumann type operator S γ via (2.6), representing physically the contributions on the flow from the rock matrices, i.e, to solve one Robin sub-domain problem per sub-domain. Thus the GMRes algorithm is implemented in the matrix-free context [21,23,41]. One can easily observe that the evaluation of S γ dominates the total computational costs in (2.8). In practice, this step is done in parallel and involves inter-processor communication across the fractureinterface [21]. To present the evaluating algorithm of S γ , we let D T h,i : V h,i ·n i | γ → M h,γ be the L 2 -orthogonal projection from the normal trace of the velocity space onto the mortar space normal trace of the velocity space in sub-domain Ω i , i = 1, 2, onto the pressure space on the fracture M h,γ . We then summarize the evaluation of the interface operator by the following steps:

Non-linear interface iterations
In this section, we form two linearization-domain-decomposition algorithms to solve the mixed interface Darcy-Forchheimer problem (2.5). For the linearization (outer) of (2.5), a first algorithm is based on a fixed-point method is presented and the second one is based on Newton-GMRes method [32,35]. For the solver of the inner systems (domain decomposition systems), both methods uses the GMRes method (2.8) to solve the reduced mixed interface problems. Note that the two approaches have competitive performance for such nonlinear model problems and they lead to different applications of the multiscale flux basis functions of Section 4.

Method 1: fixed-point-GMRes
We consider first a standard fixed-point approach to solve the interface Darcy-Forchheimer problem (2.5) (see [35]). Given an initial value u h,γ , being the solution of a linear Darcy, for k = 1, 2, . . . , until convergence, find (u This process is linear and can be solved using GMRes method (2.8), where each iteration needs to set up the action of the Robin-to-Neumann operator S γ using Algorithm 2.1. The above fixed-point-GMRes algorithm is iterated until a fixed-point residual tolerance reaches some prescribed value. The result of this procedure is then used to generate the solution in the sub-domains via for i = 1, 2, requiring two additional sub-domain solves, and where p h,γ indicates the fracture pressure at convergence.
Remark 3.1 (An alternative to (3.1)). A well-known drawback of GMRes algorithm for solving the interfacefracture problem (3.1) is that the number of iterations depends essentially on the number of sub-domain solves. A preconditioner is then necessary to reduce the iterations number to a reasonable level. To this aim, it is possible to reformulate (3.1) into a primal problem: at the iteration k ≥ 1, by solving for the sole scalar unknown p which can be discretized with a cell-centered finite volume method, leading to a symmetric and positive definite system that can be solved with a CG method. The CG method can be equipped with a preconditioner being the inverse of the discrete counterpart of the operator −∇ τ · [−K(p (k−1) )∇ τ ] (see [4,1] for more details).
Remark 3.2 (The total computational costs). The total computational costs in the inner-outer iterative approach (3.1) is dominated by the number of sub-domain solves required. Precisely, the total number of subdomain solves is given by where N lin is the number of iterations of the fixed-point procedure as outer-loop algorithm, and N k dd denotes the number of inner loop domain decomposition iterations (GMRes) at the fixed-point iteration k ≥ 1.

Method 2: Newton-GMRes
In the second approach, we propose Newton's method to solve the interface Darcy-Forchheimer problem (2.5). For simplicity of notation, we introduce the following The non-linear variational form (2.4) may be rewritten in the following canonical form: where F γ is the residual expression from the mixed system given as follows: In the next step, we calculate the Jacobian given by J γ (u h,γ in the directions of δu h,γ and δp h,γ , respectively. This is can be formally obtained by computing

This definition yields
where ⊗ denotes the standard tensor product. In each Newton iteration, we solve the following linear variational problem: find δu h,γ ∈ V h,γ and δp h,γ ∈ M h,γ , such that In compact form, the linear system for the Newton step has the following mixed structure (3.7) The interface system (3.7) is then solved with the GMRes iterations (2.8). On each GMRes iteration, we need to evaluate the action of the Robin-to-Neumann operator S γ using Algorithm 2.1. The solution of the interface problem is therefore obtained in an iterative fashion using the following update equations until the Newton residual reaches some prescribed tolerance: The result of this iterative approach is then used to infer the solution in the sub-domains using (3.2), which needs two additional sub-domain solves. Remark 3.3 (An alternative to (3.6)). For the mixed Jacobian problem in the fracture (3.6), it is possible to adopt the idea introduced in Remark 3.2 to reduce the computational cost by reformulating (3.6) into a cellcentered finite volume problem with the pressure step δp h,γ as the sole variable. The resulting system is also symmetric definite and positive and can be solved with the CG method equipped with a local preconditioner.

Outer-inner interface iterations with multiscale flux basis
As noticed previously, the dominant computational cost in the above linearization-domain-decomposition procedures comes from the sub-domain solves to evaluate the action of S γ using Algorithm 2.1 (step 2(b)). We recall that the number of sub-domain solves required by each method is approximately equal to where N lin is the number of iterations of the linearization procedure, and N k dd denotes the number of domain decomposition iterations (GMRes or any Krylov solver). Even though all sub-domain solves can be computed in parallel, this still be very costly; first, as the non-linear interface solver may converge very slowly and, second, that at each linearization iteration the condition number of the linearized interface problem ((3.1) for Method 1 and (3.6) for Method 2) is large due to a highly refined mesh.
One way to reduce the computational costs, is to employ the multiscale flux basis, following [23]. The motivation of these techniques in this work stems from eliminating the dependency between the total number of solves and the employed outer-inner procedure on the interface-fracture. This is easily achieved by precomputing and storing the flux sub-domain responses, called multiscale flux basis, associated with each fracture pressure degree of freedom on each sub-domain.
The multiscale flux basis requires solving a fixed number of linear sub-domain solves and permits retrieving the action of S γ on M h,γ by simply taking a linear combination of multiscale flux basis functions. As a result, the number of sub-domains solves is now independent of the used linearization procedure as well as of the used solver for the inner domain decomposition systems. In practice, the number of sub-domains solves will be reduced if

Multiscale flux basis
Following [23], we define (Φ h,γ ) N h,γ =1 to be the set of basis functions on the interface pressure space M h,γ , where N h,γ is the number of pressure degrees of freedom on sub-domain γ. As a result, on the fractureinterface, we have We compute the multiscale flux basis functions corresponding to (Φ h,γ ) N h,γ =1 using the following algorithm: 3. Form the multiscale flux basis for sub-domain Ω i , i.e, Once the multiscale flux basis functions are constructed for each sub-domain, the action of interface operator S RtN i , and then also the action of S γ via (2.6), is replaced by a linear combination of the multiscale flux basis functions Ψ h,γ,i . Specifically, for an interface datum ϕ h,γ ∈ M h,γ , we have ϕ h,γ = N h,γ =1 ϕ h,γ Φ h,γ , and for i = 1, 2, Remark 4.2. We observe that each fracture pressure basis function Φ h,γ on the fracture-interface corresponds to exactly two different multiscale flux basis functions, one for Ω 1 and one for Ω 2 . For the case of a fractures network, say γ = ∪ i =j γ ij , where γ ij is the fracture between the sub-domain Ω i and Ω j , the previous basis reconstruction is then applied independently on each fracture.

Application on intersecting fractures model: solving the DFNs system
In this part, we first introduce and describe the case of intersecting fractures, and then we provide our amendments to the previous algorithms.

Mathematical model
For the sake of simplicity, we consider the Darcy-Forchheimer model in a two-dimensional geological domain made up with three sub-domains Ω i , i = 1, 2, 3, physically subdivided by fractures γ i,j , 1 ≤ i < j ≤ 3. The rock matrix is now defined as Ω = each fracture γ i,j , with unknowns denoted by (u γi,j , p γi,j ). See Figure 3 (left) as an example.
They are coupled using the Robin boundary conditions given by for 1 ≤ i < j ≤ 3, where the coefficient α i,j can now be different in each fracture. To close the system, we need to impose transmission conditions between the fractures at the (d − 2)-dimensional interface T . On the intersection T , we set − u γi,j · n i,j + α γi,j p γi,j = α γi,j p T on T, (4.2a) 1≤i<j≤3 u γi,j · n i,j = 0 on T, (4.2b) for 1 ≤ i < j ≤ 3 where n i,j is the outer unit normal vector to ∂γ i,j .
For the partition of the sub-domain Ω i , 1 ≤ i ≤ 3, and the fractures γ i,j , 1 ≤ i < j ≤ 3, we extend the notation introduced in Subsection 2.1. We let T h,i be a partition of the sub-domain Ω i into 2-dimensional simplicial elements and let T h,γi,j to be a partition of the fracture γ i,j into 1-dimensional simplicial elements. Again, the meshes T h,i , 1 ≤ i ≤ 3, are allowed to be non-conforming on the fractures γ i,j , 1 ≤ i < j ≤ 3, but also different from those used in γ i,j , 1 ≤ i < j ≤ 3 (see Figure 2 for more details). We also extend the same notation for the approximation spaces in the sub-domains and in the fractures, and additionally we let M h,T be the space endowed with constant functions on T .

Domain decomposition formulation
The extension of the reduced interface problem (2.5) to the present intersecting fractures setting is as follows: On each fracture, the Robin-to-Neumann operator S γi,j and the linear functional g γi,j , 1 ≤ i < j ≤ 3, are now given by The above problem can be seen as a DFNs system on the set of fractures, and as a domain decomposition problem between the 1-dimensional fractures γ i,j , 1 ≤ i < j ≤ 3, cf. [38,39,20] for more details.

Iterative procedure
We propose to solve the non-linear domain decomposition problem (4.3) using the fixed-point approach in Subsection 3.1. This iterative process is now equipped with the multiscale flux basis of Section 4 to lessen the interface iterations. To this aim, we introduce Applying the fixed-point algorithm on the set of interface Darcy-Forchheimer equations (4.3) can be interpreted as follows: at the iteration k ≥ 1, we solve using GMRes method until a fixed tolerance is reached. Again, the evaluation of S γ in each interface GMRes iteration dominates the total computational costs of this outer-inner procedure. Note that each inner iteration also requires the evaluation of the Dirichlet-to-Neumann operator S T , which requires solves in the fractures. The complete algorithm when equipped with multiscale flux basis is now given by the following algorithm.
While i ≤ 3. 5. Given an initial guess u ii. Compute the jump across all the fractures:

Numerical examples
In this section, we validate the model and analysis presented in the previous parts by means of numerical test cases. We have chosen three examples designed to show how the proposed linearization-domaindecomposition approaches equipped with multiscale flux basis behaves vs the standard ones in various physical and geometrical situations. To compare these approaches, the main criteria considers the number of solutions of the higher-dimensional sub-problems since it constitutes the major computational cost. We consider solving the problem in the network of fractures as negligible. Since each of the higher-dimensional sub-problem is linear and will be solved many times, we consider an LU-factorization of the system matrix and a forward-backward substitution algorithm to compute the numerical solution. It results in a computational cost reduced to O(n 2 ) flops each time, where n is the size of the matrix. For bigger systems, an iterative scheme is preferable. We use the PorePy [27] library, which is a simulation tool for fractured and deformable porous media written in Python. PorePy uses SciPy [26] as default sparse linear algebra. All the examples are reported in the GitHub repository of PorePy, we want to stress again that even if we focus on lowest-order Raviart-Thomas-Nédélec finite elements, our implementation is agnostic with respect to the numerical scheme.
For the multiscale flux basis scheme presented in Section 4, for a fixed rock matrix grid and normal fracture permeability it is possible to compute once all the basis functions. The results in the next parts should be read under this important property of the method, thus in many cases only a pure fracture network will be solved at a negligible computational cost. The multiscale basis functions are computed and stored in an offline phase prior the simulation (called online).
Unless otherwise noted, the tolerance for the relative residual in the inner GMRes algorithm is taken to be 10 −6 . The same tolerance is chosen for the outer Newton/fixed-point algorithm. We consider an LU-factorization of the fracture network matrix [8,11] as the preconditioner of the GMRes method.
In the examples, we use the abbreviation MS when the linearization-domain-decomposition approach is equipped with multiscale flux basis techniques, and DD for the corresponding classical approach.

Problem setting
To validate the performance of the two proposed algorithms, we consider the first problem presented in the benchmark study [15]. The unit square domain Ω, depicted in Figure 4, has unitary permeability of the rock matrix and it is divided into 10 sub-domains by a set of fractures with fixed aperture equal to 10 −4 . At the boundary, we impose zero flux condition on the top and bottom, unitary pressure on the right, and Contrary to what has been done in the benchmark paper, we consider four different scenarios for the fracture permeabilities, by having high or low values in the tangential and normal parts. Thus, we have the case (i) with high permeable fractures, case (ii) has low permeable fractures, while cases (iii) and (iv) have mixed high and low permeability in normal and tangential directions. See Table 1 for a summary of the fracture permeability in each case. Case (i) and (ii) have the same permeabilities used in the benchmark paper [15].
In the following examples, we consider the maximal number of rock matrix solves to be 10 4 , and we mark with ∞ if this is exceeded.  5.2 Darcy model: β γ = 0 The first example considers the Forchheimer coefficient set to zero, thus the problem becoming linear. The results for different level of discretization are reported in Table 2. We indicate by level 1 a grid with 110 triangles and 26 segments, level 2 with 1544 triangles and 84 segments, and level 3 with 3906 triangles and 138 segments.  Table 2: Total number of the higher-dimensional problem solves for the case study of example in Subsection 5.2. For each level of refinement cases marked in † share the same multiscale flux basis, which can be constructed only once. The same is valid for §. Table 2 shows the results of this example for the physical considerations of Table 1. We notice that for high permeable fractures (case (i) and (iii)), the standard domain decomposition method performs better than our method with multiscale flux basis, while the opposite occurs for low permeable fractures. A possible explanation is related to the ratio between normal and tangential permeability. The normal permeability determines how strong the flux exchange is between the rock matrix and the fractures (thus, the communications at each DD iteration), while for small values of the tangential permeability the fractures are more influenced by the surrounding rock matrices. The opposite occurs in the case of high tangential permeability. Additionally, the choice of the preconditioner for DD slightly goes in favor of high permeable fractures due to the dominating role of the fracture flow in the system. We also recall that the number of higher-dimensional problem solves does not depend on the number of outer-inner interface iterations, but only on the number of local mortar degrees of freedom on the fractures network. A further important result in this experiments, is that case (i) and (iv) share the same value of α γ , thus the multiscale flux basis are computed only once per level of refinement. The same applies to case (ii) and (iii). As a result, the developed method is globally more efficient than the classical approach. That is, the results in Table 2 shows a reduction of the number of the higher-dimensional problem solves from 195 to 56 for level 1, from 212 to 186 for level 2, and from 312 to 280 for level 3. Note that the two methods produce the same solution for all the cases, within the same relative convergence tolerance. The numerical solution for all cases is reported in Figure 5.
The next series of numerical experiments aims at assessing the stability of the domain decomposition approach with respect to GMRes tolerance. The multiscale flux basis approach provides the extra flexibility to do such analysis with negligible costs, by reusing the stored multiscale flux basis used for the results of Table 2 but now with different tolerance for GMRes. Further, this set of test cases aims assessing how the overall gain for an entire simulation in terms of number of higher-dimensional problem solves can be appreciated or depreciated with more or less stringent stopping criteria for GMRes; this is a preparatory step to address the complete approaches of Section 3 for the full nonlinear problem, which requires several solves of linear Darcy problems, for which one should formulate the stopping criteria very carefully. In Table 2, we have considered the relative residual to be below 10 −6 , while in Table 3 we present the results in the case of 10 −4 and 10 −8 . Based on the results of Table 3, we can conclude that even with less stringent criterion, a considerable gain in terms of number of higher-dimensional problem solves can be achieved. We also see that all the results are free of oscillations and neither the fracture, barrier, or the very different tangential and normal permeabilities pose any problems for the domain decomposition approach. Based on the above results and in what follows we consider 10 −6 as tolerance for the GMRes algorithm.

Forchheimer model
In this second example we consider case (i) and (iii) for the fracture permeabilities since the Forchheimer model requires high permeable fractures. In this problem, we fix the computational grid level 2 of Table 2 and we change the value of β γ in order to compare the performances of Method 1 and Method 2 with and without multiscale flux basis. The Forchheimer coefficient here varies as {1, 10 2 , 10 4 , 10 6 }. These values are reasonable since in our model we do not explicitly scale β γ by the aperture, as done in [19,29]. Therefore, the last two values are more realistic. The stopping criteria for both methods is based on the relative residual criteria (4.6) with a threshold fixed as 10 −6 . The initial guess is taken by solving the linear Darcy by taking β γ equal to zero.  For Method 1, the number of higher-dimensional problem solves is reported in Table 4. As expected, Method 1 equipped with multiscale flux basis (MS) performs all the higher-dimensional problem solves in the offline phase, thus the outer-inner interface iterations for the resulting fracture network problem do not influence the total computational costs. On the contrary, the computational costs of the classical approach (DD) is influenced by the non-linearity, by varying β γ , as well as by the ratio of the normal and tangential permeabilities, by varying K γ and α γ . Particularly, the total gain of the new approach is more significant when the non-linear effects becomes more important (by increasing the value of β γ ). Furthermore, for the entire simulation of each case of Table 4, the multiscale flux basis are computed only once. As a conclusion, the entire simulation of case (i) required for Method 1 1600 higher-dimensional problem solves, while for Method 1 with multiscale flux basis, this number is reduced by 95%. For case (iii), we reduce the computational costs by 88%.
The numerical solution for two values of β γ are reported in Figure 6 for both cases. Despite the different values of β γ , we notice that the graphical results are very similar in the case of low α γ . While for high value of α γ , the resulting apparent permeability given by K γ (1 + K −1 γ β γ |u γ |) −1 decreases (for a fixed |u γ |) and the fractures are less prone to be the main path for the flow. Also as stated previously, since we do not explicitly scale β γ by the aperture, values of β γ > 10 4 are more likely for real applications.
For Method 2, involving Newton's method for the linearization step, the number of higher-dimensional problem solves are reported in Table 5. As expected, Method 2 is more efficient than Method 1 in terms of On the left, we consider value of the Forchheimer coefficient equal to β γ = 1 and on the right a high value β γ = 10 6 . In all the cases, the velocity is represented by arrows (purple for the fractures) proportional to its magnitude.  the number of higher-dimensional problem solves required, regardless of using multiscale flux basis in the domain decomposition algorithm. Again, the number of solves for the classical approach(DD) depends on the used parameters. This table demonstrates (as shown with Method 1) that as the value of β γ is increased, there is a point after which Method 2 with multiscale flux basis is more efficient than without multiscale flux basis. In that case, the gain in the number of solves becomes more significant when decreasing the value of α γ . Note that, in practice, the simulations for Method 2 with multiscale flux basis are performed with negligible computational costs as we reused the flux basis inherited from Method 1. This point together with the fact that the total number of solves required by the entire simulation of case (i) is now reduced by 57% as well as that of case (i) is reduced by 80% showcase the performance of Method 2 with multiscale flux basis.
To sum up, equipping Method 1 and 2 with multiscale flux basis leads to powerful tools to solve complex fracture network with important savings in terms of the number of higher-dimensional problem solves. Note that, as known, one limitation of Method 2 involving Newton method is that a good initial value is usually required to obtain a solution. A good combination of both methods can also be used, in which one can perform first some fixed-point iterations and then switch to Newton method. Concerning the computational costs, let us point out that the fixed-point algorithm of Method 1 requires at each iteration the assembly of the matrix corresponding to the linearization of the Darcy-Forchheimer equations and the solution of a linear system. The Newton method in Method 2 is slightly more expensive since one has to assemble two matrices at each iteration and to update the right-hand side.

Heterogeneous Forchheimer model
In this example we assign to the two biggest fractures (one horizontal and one vertical) high permeability while to the others low permeability. For the highly permeable fractures we adopt the physical parameters of case (i), while for those with lower permeabilities, the physical parameters corresponding to case (ii) together with zero Forchheimer coefficient. In this case, we want to test the applicability of Method 1 with and without multiscale flux basis on highly heterogeneous setting for both the permeability and the flow models. We consider level 2 for the computation and, subsequently, the local grids of the low permeable fractures are coarsened by having half of the original number of elements resulting in 60 mortar unknowns instead of 84.
As usually, we compare the method with and without multiscale flux basis in terms of the number of higher-dimensional problem solves. The results are represented in  can see that the classical approach is outperformed with the approach equipped with multiscale flux basis, particularly, the total computational costs is drastically reduced when the non-linear effects becomes more important. The entire simulation of Table 6 required 2708 higher-dimensional problem solves for the classical approach while the same approach equipped with multiscale flux basis required 62 solves. The overall gain is then of 94% which can also be appreciated for level 3. Similar conclusions as above can be drawn for Method 2, namely in terms of reduction of the solves (not shown). An example of solution is given in Figure  7.

Generalized Forchheimer model
As stated previously, another advantage distinguishes our approach is that it can integrate easily more complex problems. Here, we apply our procedure to a more general model describing the pressure-flow , we have an exponential relation between the fluid viscosity and the pressure. We consider problem (2.1) where the non-linear term is as follows where ζ being a model parameter. Thus, the non-linear effects are now dependent on both the pressure and the velocity. For a more detailed discussion we refer to [40]. For the present setting, the fracture permeabilities are set as in case (i) and (iii) of Table 1.
For the discretization of the mixed geometry, we consider level 2. We use Method 1 with and without multiscale basis functions. Also, it was not necessary to recompute the basis functions, since we can reuse the stored multiscale flux basis from the previous test case and solve then only on the fracture network the above more complex Darcy-Forchheimer model. The total number the higher-dimensional problem solves for β γ = 20 and ζ ∈ (0.5, 5, 7.5) is reported in Table 7. As expected, for such a strong non-linearity, the results case (i) case (iii) ζ MS DD MS DD 0.5 86 † (5) 71 (5) 86 § (4) 176 (4) 5 86 † (4) 648 (4) 86 § (6) ∞ 7.5 86 † (3) 5317 (3) 86 § (4) ∞ Table 7: Total number of the higher-dimensional problem solves required by Method 1 for the case study in Subsection 5.5. The number of the fixed-point iterations are in brackets. Within each case the construction of the multiscale basis is done only once, we mark by † (respectively §) common computations.
shows that a considerable gain in terms of higher-dimensional problem solves can be achieved. Particularly, for large values of ζ the classical approach becomes uncompetitive to the new approach. In Figure 8 we report the solution for ζ = 5.

Conclusions
In this work, we have presented a strategy to speed up the computation of a Darcy-Forchheimer model for flow and pressure in fractured porous media by means of multiscale flux basis, that represent the interdimensional flux exchange. The scheme transforms a computationally expensive discrete fracture model to a more affordable discrete fracture network, where in the latter only a co-dimensional problem is solved. The multiscale flux basis are computed in an offline stage of the simulation and, despite the particular choice done in this paper, are completely agnostic to the model in the fracture network. The numerical results show the speed-up gain compared to a more classical linearization-domain-decomposition approaches, where solves in both the matrix and the fracture network are required along the entire outer-inner iterative method. Crucially, an important number of the outer-inner interface iterations may be spared.
With the proposed approach we are able to predict the computational effort needed to solve the problem since it is directly related to the number of mortar grids in the fracture network. Furthermore, the multiscale flux basis can be reused when the fracture network geometry, rock matrix properties, and normal permeability are fixed. Theoretical findings and numerical results show the validity of the proposed approach and of its aforementioned properties.
Even if not explicitly considered in this work, it is possible to further increase the efficiency of the proposed scheme by the following two steps. First, compute a multiscale flux basis only in the related connected part of the rock matrix. Second, use an adaptive stopping criteria for the inner-outer iterative method based on a posteriori error estimates. These enhancements are a part of future work along with the extension in three-dimensions.