Algebraic multiscale method for flow in heterogeneous porous media with embedded discrete fractures (F-AMS)

TL;DR

This paper presents an algebraic multiscale method (F-AMS) for simulating flow in heterogeneous porous media with embedded fractures, offering scalable, accurate, and mass-conservative solutions adaptable to complex fracture networks.

Contribution

The novel F-AMS framework constructs independent multiscale coarse grids for matrix and fractures, with flexible coupling strategies and proven scalability for complex heterogeneous media.

Findings

F-AMS is insensitive to high contrast in conductivities.

The method achieves efficient convergence with reduced CPU times.

It provides accurate, mass-conservative flux reconstructions after iterations.

Abstract

This paper introduces an Algebraic MultiScale method for simulation of flow in heterogeneous porous media with embedded discrete Fractures (F-AMS). First, multiscale coarse grids are independently constructed for both porous matrix and fracture networks. Then, a map between coarse- and fine-scale is obtained by algebraically computing basis functions with local support. In order to extend the localization assumption to the fractured media, four types of basis functions are investigated: (1) Decoupled-AMS, in which the two media are completely decoupled, (2) Frac-AMS and (3) Rock-AMS, which take into account only one-way transmissibilities, and (4) Coupled-AMS, in which the matrix and fracture interpolators are fully coupled. In order to ensure scalability, the F-AMS framework permits full flexibility in terms of the resolution of the fracture coarse grids. Numerical results are presented for two- and three-dimensional heterogeneous test cases. During these experiments, the performance of F-AMS, paired with ILU(0) as second-stage smoother in a convergent iterative procedure, is studied by monitoring CPU times and convergence rates. Finally, in order to investigate the scalability of the method, an extensive benchmark study is conducted, where a commercial algebraic multigrid solver is used as reference. The results show that, given an appropriate coarsening strategy, F-AMS is insensitive to severe fracture and matrix conductivity contrasts, as well as the length of the fracture networks. Its unique feature is that a fine-scale mass conservative flux field can be reconstructed after any iteration, providing efficient approximate solutions in time-dependent simulations.

Tables

Table 1. Table 1: F-AMS solution algorithm

Repeat the following steps until convergence to the desired accuracy is reached:

1. Initialize: $p^{\nu }\leftarrow p^{\nu +1}$

2. Update linear system entries: ${𝑨}^{\nu }$ and $q^{\nu }$ in (17)

3. Update residual: $r^{\nu }=q^{\nu }-{𝑨}^{\nu }{p}^{\nu }$

4. Compute (or adaptively update) $𝒫$: follow either coupling strategy from

Subsection 3.1.

5. Multiscale Stage: $\delta p^{\nu +1/2}=𝓟{\left(𝓡{𝑨}^{\nu }𝓟\right)}^{-1}𝓡r^{\nu }$

6. Update residual $r^{\nu +1/2}=r^{\nu }-{𝑨}^{\nu }\delta p^{\nu +1/2}$

7. Smoothing Stage: $\delta p^{\nu +2/2}={𝑴}_S^{-1}{r}^{\nu +1/2}$

8. Update solution: $p^{\nu +1}=\left(p^{\nu }+\delta p^{\nu +1/2}+\delta p^{\nu +2/2}\right)$

Table 2. Table 2: Distance-based fracture coarsening algorithm.

Repeat the following for each fracture network, $f_i$, which has $N_{cells}$ fine cells:

1. Choose $d_{min}$, the minimum distance between two fracture coarse nodes.

2. Initialize three empty queues: $Q_{vertex},Q_1,Q_2$.

3. Initialize two vectors of length $N_{cells}$: $level$ (set to $\infty$) and $primal$ (set to $0$).

4. $N_{primal}\leftarrow 0$

5. Choose a starting cell from $f_i$ and add it to $Q_{vertex}$.

Repeat until $Q_{vertex}$ is empty:

6. $vertex\leftarrow$ extract top of $Q_{vertex}$

7. $N_{primal}\leftarrow N_{primal}+1$

8. Create primal block number $N_{primal}$ with $vertex$ as its coarse node.

9. $primal\left[vertex\right]\leftarrow N_{primal}$ and $level\left[vertex\right]\leftarrow 0$.

10. Add $vertex$ to $Q_1$.

For $dist$ from $1$ up to and including $d_{min}$:

Repeat until $Q_1$ is empty:

11. $cell\leftarrow$ extract top of $Q_1$

For each neighbour of $cell$, $neig{h}_j$, with $level\left[neig{h}_j\right]>dist$:

12. remove $neig{h}_j$ from $Q_{vertex}$, if it is present (i.e. $level\left[neig{h}_j\right]$ = $d_{min}$)

13. $level\left[neig{h}_j\right]\leftarrow dist$

14. $primal\left[neig{h}_j\right]\leftarrow N_{primal}$

15. add $neig{h}_j$ to $Q_2$.

16. swap $Q_1$ and $Q_2$.

17. empty $Q_1$ into $Q_{vertex}$.

At this point, $primal\left[i\right]$ gives the index of fine cell $i$’s primal block, while $level\left[i\right]$ is

the distance from cell $i$ to the nearest vertex. The fine cells which were not marked

as vertices will form edges on the dual-coarse grid.

Table 3. Table 3: Performance of F-AMS during the scale sensitivity test cases, when using different coarsening factors. The middle row for each test (shown in bold) is the optimum configuration, whose results were presented in the body of the manuscript.

Scale	Coarsening ratio (matrix, fracs)	Total CPU time (sec)	# iterations
${32}^3$	$2\times 2\times 2$, $2\times 2$	$1.904$	$11$
	$\mathrm{𝟔}\times \mathrm{𝟔}\times \mathrm{𝟔}$, $\mathrm{𝟔}\times \mathrm{𝟔}$	$0.330$	$\mathrm{𝟐𝟑}$
	$8\times 8\times 8$, $8\times 8$	$0.351$	$29$
${64}^3$	$4\times 4\times 4$, $4\times 4$	$3.270$	$13$
	$\mathrm{𝟔}\times \mathrm{𝟔}\times \mathrm{𝟔}$, $\mathrm{𝟔}\times \mathrm{𝟔}$	$2.186$	$\mathrm{𝟏𝟗}$
	$9\times 9\times 9$, $9\times 9$	$2.371$	$27$
${128}^3$	$6\times 6\times 6$, $6\times 6$	$21.790$	$17$
	$\mathrm{𝟖}\times \mathrm{𝟖}\times \mathrm{𝟖}$, $\mathrm{𝟖}\times \mathrm{𝟖}$	$17.620$	$\mathrm{𝟐𝟐}$
	$11\times 11\times 11$, $11\times 11$	$21.350$	$31$
${256}^3$	$8\times 8\times 8$, $8\times 8$	$164.600$	$18$
	$\mathrm{𝟏𝟎}\times \mathrm{𝟏𝟎}\times \mathrm{𝟏𝟎}$, $\mathrm{𝟏𝟎}\times \mathrm{𝟏𝟎}$	$150.400$	$\mathrm{𝟐𝟑}$
	$17\times 17\times 17$, $17\times 17$	$252.100$	$40$

Equations

$$\displaystyle\bigg{[}\frac{\partial(\phi\rho)}{\partial t}-\nabla\cdot(\rho\lambda\cdot\nabla p)\bigg{]}^{m}$$

$$\displaystyle\bigg{[}\frac{\partial(\phi\rho)}{\partial t}-\nabla\cdot(\rho\lambda\cdot\nabla p)\bigg{]}^{f}$$

$$q_V^{mw}=PI{\lambda }^m(p^w-p^m)/V\equiv {\beta }^m(p^w-p^m),$$

$$q_V^{mf}=CI{\lambda }^{f-m}(p^f-p^m)/V\equiv {\eta }^m(p^f-p^m)$$

$$q_A^{fm}=CI{\lambda }^{f-m}(p^m-p^f)/A\equiv {\eta }^f(p^m-p^f),$$

$${\int }_V{q}_V^{mf}dV=-{\int }_A{q}_A^{fm}dA.$$

$$C{I}_{ij}=\frac{A_{ij}}{⟨d{⟩}_{ij}},$$

$$\displaystyle\bm{A}^{\nu}p^{\nu+1}\equiv{\left[\begin{array}[]{ccc}\bm{A}^{mm}&\bm{A}^{mf}&\bm{A}^{mw}\\ \bm{A}^{fm}&\bm{A}^{ff}&\bm{A}^{fw}\\ \bm{A}^{wm}&\bm{A}^{wf}&\bm{A}^{ww}\end{array}\right]^{\nu}}\left[\begin{array}[]{c}p^{m}\\ p^{f}\\ p^{w}\end{array}\right]^{\nu+1}=\left[\begin{array}[]{c}q^{m}\\ q^{f}\\ q^{w}\end{array}\right]^{\nu}\equiv q^{\nu},$$

$$p^m\approx p^{\prime m}=\sum _{i=1}^{N_{cm}}{\Phi }_i^{mm}{\breve{p}}_i^m+\sum _{i=1}^{N_f}\sum _{j=1}^{N_{c{f}_i}}{\Phi }_j^{m{f}_i}{\breve{p}}_j^{f_i}+\sum _{k=1}^{N_w}{\Phi }_j^{mw}{\breve{p}}_k^w,$$

$$p^f\approx p^{\prime f}=\sum _{i=1}^{N_{cm}}{\Phi }_i^{fm}{\breve{p}}_i^m+\sum _{i=1}^{N_f}\sum _{j=1}^{N_{c{f}_i}}{\Phi }_j^{f{f}_i}{\breve{p}}_j^{f_i}+\sum _{k=1}^{N_w}{\Phi }_j^{fw}{\breve{p}}_k^w$$

$$p^{w_i}=p^{\prime w_i}={\breve{p}}^{w_i}\forall i\in \{1,\cdots ,N_w\}.$$

$$p^m\approx {p^{\prime }}^m={\mathcal{P}}^m\breve{p}\equiv [{\mathcal{P}}^{mm}{\mathcal{P}}^{mf}{\mathcal{P}}^{mw}][{\breve{p}}^m{\breve{p}}^f{\breve{p}}^w{]}^T$$

$$p^f\approx {p^{\prime }}^f={\mathcal{P}}^f\breve{p}\equiv [{\mathcal{P}}^{fm}{\mathcal{P}}^{ff}{\mathcal{P}}^{fw}][{\breve{p}}^m{\breve{p}}^f{\breve{p}}^w{]}^T,$$

$$\displaystyle\bm{\mathcal{P}}^{m}=\left[\begin{array}[]{ccc|ccccccc|ccc}\vdots&\cdots&\vdots&\vdots&\cdots&\vdots&\cdots&\vdots&\cdots&\vdots&\vdots&\cdots&\vdots\\ \Phi^{mm}_{1}&\cdots&\Phi^{mm}_{N_{cm}}&\Phi^{mf_{1}}_{1}&\cdots&\Phi^{mf_{1}}_{N_{cf_{1}}}&\cdots&\Phi^{mf_{N_{f}}}_{1}&\cdots&\Phi^{mf_{N_{f}}}_{N_{cf_{N_{f}}}}&\Phi^{mw}_{1}&\cdots&\Phi^{mw}_{N_{w}}\\ \vdots&\cdots&\vdots&\vdots&\cdots&\vdots&\cdots&\vdots&\cdots&\vdots&\vdots&\cdots&\vdots\end{array}\right].$$

$$\displaystyle\bm{\mathcal{P}}^{f}=\left[\begin{array}[]{ccc|ccccccc|ccc}\vdots&\cdots&\vdots&\vdots&\cdots&\vdots&\cdots&\vdots&\cdots&\vdots&\vdots&\cdots&\vdots\\ \Phi^{fm}_{1}&\cdots&\Phi^{fm}_{N_{cm}}&\Phi^{ff_{1}}_{1}&\cdots&\Phi^{ff_{1}}_{N_{cf_{1}}}&\cdots&\Phi^{ff_{N_{f}}}_{1}&\cdots&\Phi^{ff_{N_{f}}}_{N_{cf_{N_{f}}}}&\Phi^{fw}_{1}&\cdots&\Phi^{fw}_{N_{w}}\\ \vdots&\cdots&\vdots&\vdots&\cdots&\vdots&\cdots&\vdots&\cdots&\vdots&\vdots&\cdots&\vdots\end{array}\right].$$

$$\displaystyle\bm{\mathcal{P}}=\left[\begin{array}[]{c}\bm{\mathcal{P}}^{m}\\ \bm{\mathcal{P}}^{f}\\ \bm{\mathcal{P}}^{w}\end{array}\right]=\left[\begin{array}[]{ccc}\bm{\mathcal{P}}^{mm}&\bm{\mathcal{P}}^{mf}&\bm{\mathcal{P}}^{mw}\\ \bm{\mathcal{P}}^{fm}&\bm{\mathcal{P}}^{ff}&\bm{\mathcal{P}}^{fw}\\ \bm{\mathcal{P}}^{wm}&\bm{\mathcal{P}}^{wf}&\bm{\mathcal{P}}^{ww}\\ \end{array}\right],$$

$${\delta }_{ij}=\{\begin{array}{ll}1,& \text{ if }i=j\\ 0,& \text{ if }i\ne j\end{array}.$$

$$-\nabla \cdot ({\lambda }^{\ast }\cdot \nabla {\Phi }^{\ast \bullet })+\sum _{j\in {\text{conn}}_{mf}^{\ast }}{\eta }_j^{\ast }\xi ({\Phi }^{\ast \bullet })+\sum _{j\in {\text{perf}}_w^{\ast }}{\beta }_j^{\ast }({\Phi }^{\ast \bullet }-{\Phi }^{w\bullet })=0,$$

$$\xi ({\Phi }^{\ast \bullet })=0\forall {\Phi }^{\ast \bullet }\in \{{\Phi }^{mm},{\Phi }^{mf},{\Phi }^{mw},{\Phi }^{fm},{\Phi }^{ff},{\Phi }^{fw}\},$$

$$\xi ({\Phi }^{ff})=0$$

$$\xi ({\Phi }^{mf})={\Phi }^{mf}-{\Phi }^{ff}.$$

$$\xi ({\Phi }^{mm})={\Phi }^{mm}.$$

$$\xi ({\Phi }^{mm})=0.$$

$$\xi ({\Phi }^{fm})={\Phi }^{fm}-{\Phi }^{mm}.$$

$$\xi ({\Phi }^{ff})={\Phi }^{ff}.$$

$$\xi ({\Phi }^{mm})={\Phi }^{mm}-{\Phi }^{fm},$$

$$\xi ({\Phi }^{ff})={\Phi }^{ff}-{\Phi }^{mf}$$

$$T_{ratio}=\frac{⟨T{⟩}_{frac}}{⟨T{⟩}_{rock}},$$

$$\displaystyle\underbrace{\left[\begin{array}[]{cc}\bm{A}^{mm}&\bm{A}^{mf}\\ \bm{A}^{fm}&\bm{A}^{ff}\end{array}\right]}_{\bm{A}}\underbrace{\left[\begin{array}[]{c}p^{m}\\ p^{f}\end{array}\right]}_{p}=\underbrace{\left[\begin{array}[]{c}q^{m}\\ q^{f}\end{array}\right]}_{q}.$$

$$\displaystyle{\bm{\wp}}\equiv{\left[\begin{array}[]{cc}\bm{\wp}^{m}&0\\ 0&\bm{\wp}^{f}\end{array}\right]},$$

$$\displaystyle\underbrace{\left[\begin{array}[]{cccc|ccc}\bm{A}^{{I^{m}}{I^{m}}}&\bm{A}^{{I^{m}}{F^{m}}}&0&0&\bm{A}^{{I^{m}}{F^{f}}}&\bm{A}^{{I^{m}}{E^{f}}}&\bm{A}^{{I^{m}}{V^{f}}}\\ \bm{A}^{{F^{m}}{I^{m}}}&\bm{A}^{{F^{m}}{F^{m}}}&\bm{A}^{{F^{m}}{E^{m}}}&0&\bm{A}^{{F^{m}}{F^{f}}}&\bm{A}^{{F^{m}}{E^{f}}}&\bm{A}^{{F^{m}}{V^{f}}}\\ 0&\bm{A}^{{E^{m}}{F^{m}}}&\bm{A}^{{E^{m}}{E^{m}}}&\bm{A}^{{E^{m}}{V^{m}}}&\bm{A}^{{E^{m}}{F^{f}}}&\bm{A}^{{E^{m}}{E^{f}}}&\bm{A}^{{E^{m}}{V^{f}}}\\ 0&0&\bm{A}^{{V^{m}}{E^{m}}}&\bm{A}^{{V^{m}}{V^{m}}}&\bm{A}^{{V^{m}}{F^{f}}}&\bm{A}^{{V^{m}}{E^{f}}}&\bm{A}^{{V^{m}}{V^{f}}}\\ \hline\cr\bm{A}^{{F^{f}}{I^{m}}}&\bm{A}^{{F^{f}}{F^{m}}}&\bm{A}^{{F^{f}}{E^{m}}}&\bm{A}^{{F^{f}}{V^{m}}}&\bm{A}^{{F^{f}}{F^{f}}}&\bm{A}^{{F^{f}}{E^{f}}}&0\\ \bm{A}^{{E^{f}}{I^{m}}}&\bm{A}^{{E^{f}}{F^{m}}}&\bm{A}^{{E^{f}}{E^{m}}}&\bm{A}^{{E^{f}}{V^{m}}}&\bm{A}^{{E^{f}}{F^{f}}}&\bm{A}^{{E^{f}}{E^{f}}}&\bm{A}^{{E^{f}}{V^{f}}}\\ \bm{A}^{{V^{f}}{I^{m}}}&\bm{A}^{{V^{f}}{F^{m}}}&\bm{A}^{{V^{f}}{E^{m}}}&\bm{A}^{{V^{f}}{V^{m}}}&0&\bm{A}^{{V^{f}}{E^{f}}}&\bm{A}^{{V^{f}}{V^{f}}}\end{array}\right]}_{\bm{\wp A\wp^{T}}}\underbrace{\left[\begin{array}[]{c}p^{I^{m}}\\ p^{F^{m}}\\ p^{E^{m}}\\ p^{V^{m}}\\ \hline\cr p^{F^{f}}\\ p^{E^{f}}\\ p^{V^{f}}\\ \end{array}\right]}_{\bm{\wp}p}=\underbrace{\left[\begin{array}[]{c}q^{I^{m}}\\ q^{F^{m}}\\ q^{E^{m}}\\ q^{V^{m}}\\ \hline\cr q^{F^{f}}\\ q^{E^{f}}\\ q^{V^{f}}\\ \end{array}\right]}_{\bm{\wp}q}.$$

Full text

††thanks: m.tene@tudelft.nl††thanks: malkobaisi@pi.ac.ae††thanks: h.hajibeygi@tudelft.nl

1 Introduction

In many geoscience applications, including hydrocarbon production and geothermal energy exploitation, the target formations are naturally fractured. These formations are often highly heterogeneous matrix rock is crossed by several networks of lower-dimensional highly-conductive fractures at multiple length scales [1]. This raises important challenges for flow simulation, motivating the development of advanced modelling and numerical solution techniques.

Among the proposed methods, the hierarchical fracture modelling approach allows for avoiding complexities associated with the discretization and dynamic nature of fracture geometries [2]. In this approach, small-scale fractures (below the matrix grid resolution) are homogenized within the matrix rock, altering its effective permeability [3]. The remaining fractures are then represented as explicit control volumes [4, 5]. If the fracture and matrix grids are generated independently, then the formation is said to be discretized according to the Embedded Discrete Fracture Model (EDFM) [6]. Alternatively, the fracture cells can be constrained to lie at the interfaces of matrix cells, i.e. by employing Discrete Fracture Modelling methods (DFM), which require unstructured grids [7]. Both DFM and EDFM have been applied to reservoirs with complex fracture geometries [8] and fluid physics [9, 10]. Recent developments include higher-order approximation schemes within finite-volume [11] and finite-element [12] methods.

Note that the total number of degrees of freedom (DOF), even after homogenizing small-scale fractures, is beyond the scope of classical simulation methods. This motivates the development of an efficient multiscale method for heterogeneous fractured porous media.

Multiscale finite element (MSFE) [13] and finite volume (MSFV) methods [14] have been introduced and evolved mainly for heterogeneous, but non-fractured, porous media (see [15] for a comparison). Recent developments include efficient solution of the pressure equation for multi-component displacements within sequentially- [16, 17] and fully-implicit [18, 19] frameworks. They have also been extended to capture complex wells [20, 21] and to the transport equations [22]. In addition, enriched multiscale methods are targeted at media with high conductivity contrasts [23, 24, 25], with modifications to maintain their monotonicity [26, 27]. More recently, a multiscale formulation has been devised to support unstructured grids [28].

The first application of MSFV methods to fractured reservoirs was developed on the basis of coupling the matrix pressure to the average pressure in each fracture network, at coarse scale [5]. Based on their results for two-dimensional (2D) problems, this method was efficient for media with highly conductive fracture networks which span short spatial length scales (relative to that of the domain). However, convergence was observed to degrade for test cases with significant variations in the pressure distribution along the fracture network. In combination with streamline-based mixed formulations, multiscale methods have also been employed to 2D fractured reservoirs [29]. More recently, a multiscale approach was developed for 2D reservoirs which assigned one coarse node at each fracture intersection only, with no coarse nodes in the matrix [30]. Note that none of these methods include 3D heterogeneous reservoirs nor has their performance been benchmarked against a commercial linear solver. More importantly, the literature is lacking a multiscale method which allows for flexible coarse grids inside the matrix as well as its embedded fractures and, thus, able to accommodate heterogeneous cases with fracture networks of different length scales.

This paper presents the development of an Algebraic MultiScale method for heterogeneous Fractured porous media (F-AMS) using EDFM. Given a partition of the fine-scale cells into primal and dual coarse blocks for both the matrix and fracture networks, the algorithm algebraically constructs the multiscale prolongation (mapping coarse- to fine-scale) and restriction (mapping fine- to coarse-scale) operators. The prolongation operator columns are the local basis functions, solved on dual-coarse cells, for both matrix and fractures.

F-AMS supports four different matrix-fracture coupling strategies, at the coarse-scale. First, the Decoupled-AMS basis functions are defined by neglecting the contribution of a medium’s coarse solution (e.g., fractures) in the interpolated solution in the other (e.g., matrix), thus preserving sparsity in the resulting coarse-scale system. Then, two semi-coupled (one-way) strategies, Rock-AMS and Frac-AMS, are considered. The Rock-AMS approach constructs a prolongation operator in which the matrix coarse solutions also contribute in computing the interpolated fine-scale solution in neighbouring fractures. Similarly, Frac-AMS considers the influence of the fracture coarse solution when interpolating the pressure inside the surrounding porous rock. Finally, the fully coupled strategy, Coupled-AMS, is devised, where coarse-scale solutions from both media play a role in finding the fine-scale solution of each other. This last approach, although allowing for full fracture-matrix coupling, leads to a dense coarse-scale system with additional overhead during the associated algebraic (matrix-vector, matrix-matrix) operations. As such, for practical applications, the Coupled-AMS prolongation operator may require tuning via truncation, where values below a specified threshold are algebraically deleted, followed by a rescaling step, to maintain partition of unity. This option is also investigated in the paper.

To summarize, F-AMS allows for arbitrary coarse grid resolutions in both fractures and matrix, as well as all possible coarse-scale coupling between them. Furthermore, once these coarse grids are defined, the F-AMS procedure is formulated and implemented in algebraic form, in line with the previously published formulations of incompressible (AMS) [31] and compressible (C-AMS) [32] flows. In the limit, if the Frac-AMS coupling strategy and only one coarse node per fracture network is used, F-AMS automatically reduces to the method described in [5]. However, this setup proves inefficient for many of the test cases in this paper. From a bottom-up perspective, F-AMS extends the AMS prolongation operator, as previously described in [31], with additional columns. Some of these columns correspond to the enriched fracture coarse domain, as explained above. The remainder represent local well basis functions for Peaceman wells [20].

In order to test F-AMS method, a proof-of-concept implementation is developed with a focus on reservoirs defined on 3D structured grids with embedded vertical fracture plates (for the challenges associated with unstructured multiscale simulation see [28]). For the presented experiments, a finite-element (FE) restriction operator is employed to obtain a symmetric-positive-definite coarse system. If approximate (non-converged) F-AMS solutions are used, in the last iteration step, a finite-volume (FV) restriction operator is employed followed by mass-conservative reconstruction of fine-scale flux field. The performance of this in-house object oriented serial-processing simulator was measured based on CPU times, as well as convergence rates. Numerical test cases are considered in order to study the effect of the different components of the algorithm, namely, the coarsening ratios and basis function coupling strategies. F-AMS is developed as an accurate multiscale (approximate) pressure solver. In order to assess its scalability, however, its performance (based on CPU times) is compared against the commercial Algebraic MultiGrid (AMG) solver, SAMG [33], as a preconditioner. Results of these systematic studies show that only few DOF per fracture network are necessary to obtain a good trade-off between convergence rate and computational expense. In conclusion, F-AMS is found efficient and scalable for solving flow in heterogeneous and naturally fractured porous media. Its development marks an important step forward towards the integration multiscale methods as “black-box” pressure solvers within existing reservoir simulators, with the possibility of extension to more complex physics and scenarios.

The paper is structured as follows. First, the EDFM fine-scale discrete system is described in Section 2. Then, the components of the F-AMS algorithm are detailed in Section 3. Section 4 consists of numerical results for both 2D and 3D test cases. Finally, conclusions and remarks make the subject of Section 5.

2 Fine-scale discretized system

The mass-conservation equations for single-phase flow in fractured media, using Darcy’s law, can be written as

[TABLE]

for the matrix (superscript m) and

[TABLE]

for the fracture (superscript f) spatial domains. Here, the phase mobility $\lambda=k^{*}/\mu$ consists of fluid viscosity ($\mu$) and effective fracture ($k^{*}=k^{f}$) or rock ($k^{*}=k^{m}$) permeability. Note that the latter can also account for the homogenized small-scale fractures, as described in the hierarchical fracture model [4, 3]. Also, $\rho$ and $\phi$ are the fluid density and rock porosity, respectively. Here, for simplicity, fracture porosity is always considered to be $1$. The $q^{mw}$ and $q^{fw}$ denote the matrix and fracture external source terms, respectively, i.e. from injection/production wells. For a perforated matrix control volume $V$, it reads

[TABLE]

where $\beta^{m}=PI~{}\lambda^{m}/V$, $q^{mw}_{tot}=\int_{V}q^{mw}_{V}dV$ is the total injection rate, and $PI$ is the productivity index [34]. Similarly, the fracture-matrix coupling terms are modeled such that for a matrix volume $V$ intersecting with a fracture surface $A$ one obtains

[TABLE]

and

[TABLE]

where $\eta^{m}=CI~{}\lambda^{f-m}/V$ and $\eta^{f}=CI~{}\lambda^{f-m}/A$. This ensures the total flux between a fracture element of area $A$ and a matrix element of volume $V$ is equal [5], i.e.,

[TABLE]

The $\lambda^{f-m}$ is the effective mobility at the interface between the fractures and their surrounding matrix. The $CI$ is the connectivity index, defined on a discrete system as

[TABLE]

where $A_{ij}$ is the area fraction of fracture element $i$ overlapping with the matrix element $j$, and $\langle d\rangle_{ij}$ is the average distance of the two elements [5]. Finally, the $Q^{m}$ and $Q^{f}$ terms describe other external source terms for matrix and fractures (e.g., gravity terms).

These equations are to be solved for matrix and fracture pressures, $p^{m}$ and $p^{f}$, on the matrix $\Omega^{m}$ and fracture $\Omega^{f}$ domains, as depicted in Fig. 1. Note that a fracture network can consist of several fractures, which are represented in a lower dimensional space, i.e. $\Omega^{f}\subset R^{n-1}$, than the matrix (reservoir rock) $\Omega^{m}\subset R^{n}$. The main advantage of this type of formulation is that the matrix and fracture grids are independent and, thus, can be freely adapted to accommodate the appropriate physics for each medium. This is especially important in highly fractured reservoirs or when fractures are generated (and closed) during simulation, e.g., in geothermal formations [35].

The incompressible single-phase pressure solution obtained using the EDFM approach for a 2D fractured reservoir model, shown in Fig. 1(a), is provided in Fig. 1(b). Two pressure-constrained wells are placed on the East and West boundaries, and the reservoir rock is homogeneous.

When non-linearities are present (e.g., compressible flows), the flow equations need to first be linearized, i.e.,

[TABLE]

and then iteratively solved, in a Newton-Raphson loop, until the converged solution is achieved. Note that this system (17) shows an implicit treatment of the coupling between fracture and matrix through the $\bm{A}^{fm}$ entries, and that $\bm{A}$ can be non-symmetric, due to the compressibility effects [32, 36].

Developing an efficient solution strategy for the linearized system (17) is quite challenging for several reasons. On the one hand, the size of this system can exceed several millions of unknowns for realistic test cases. On the other hand, the value of the condition numbers for the system matrix is worsened by high contrasts between reservoir properties (matrix permeability is highly heterogeneous over large scales, fractures are typically much more conductive than the matrix, etc.).

Clearly a classical upscaling method cannot be employed here due to the highly resolved fractures, which play an important role in mass transport. This creates a niche for conservative multiscale methods, which have the important advantage of solving coarse-scale problems while honouring fine-scale data [37, 38] in an iterative error reduction strategy [39, 40, 41] which allows for mass-conservative flux reconstruction at any stage [42]. Next, the development of the F-AMS method is presented.

3 Algebraic Multiscale Method for Heterogeneous Porous Media with Embedded Discrete Fractures (F-AMS)

This section describes the F-AMS procedure, an efficient multiscale solution strategy for Eq. (17). Given a computational domain with $N_{f}$ fracture networks and $N_{w}$ wells, F-AMS first superimposes two coarse grids on top of both the matrix and fracture domains. The primal-coarse grid is a non-overlapping decomposition of the domain, inside which a fine-scale cell is selected as coarse node (Fig. 2(a) for 2D, and 2(c) and 2(e) for 3D cases). By connecting the coarse nodes, a secondary overlapping coarse grid is obtained, which is called the dual coarse grid (Figs. 2(b), 2(d) and 2(f)). There exist $N_{cm}$ and $N_{dm}$ matrix primal-coarse and dual-coarse blocks and, similarly, each fracture network $f_{i}$ contains $N_{cf_{i}}$ and $N_{df_{i}}$ fracture primal-coarse and dual-coarse blocks. Note that $N_{w}$ (injection or production) wells exist in the domain, as driving forces for the flow.

F-AMS approximates the solution to Eq. (17), $p$, as a superposition of coarse-scale solutions ($\breve{p}$) using locally computed basis functions ($\Phi$), i.e.

[TABLE]

for the matrix and

[TABLE]

for the fractures, respectively. The basis functions associated with matrix coarse cells (i.e., $\Phi^{m*}$) are $\Phi_{i}^{mm}$ for matrix-matrix effects, $\Phi_{j}^{mf}$ for the matrix-fracture coupling, and $\Phi_{j}^{mw}$ matrix-well interactions. These basis functions (interpolators) are employed in order to capture the effects of all the important factors (matrix, fractures, and wells) in the construction of a good approximation for the matrix pressure field $p^{m}$. Similarly for fractures, $\Phi^{f*}$ consists of the contributions from the matrix $\Phi^{fm}$, fractures $\Phi^{ff}$, and wells $\Phi^{fw}$, if present.

One of the novel aspects of this work is that the pressure field inside fractures, $p^{f}$, is included explicitly in the multiscale formulation (Eq. (19)). This means that the fracture grid cells are also decomposed into primal and dual coarse blocks (Fig. 2), similar to the matrix. Their solutions are also mapped to the coarse scale and back to the original resolution, again, similar to the matrix. More specifically, each fracture network $f_{i}$ is decomposed into $N_{cf_{i}}$ primal-coarse grid blocks, for which sets of basis functions are calculated. One could employ the same formulation for wells, i.e., discretize them into several fine-scale cells which can then be coarsened on the superimposed multiscale grids. However, for the sake of simplicity, in the experiments presented in this work, each well is assigned a single fine-scale DOF, which is mapped to the coarse-scale using the identity restriction operator, i.e.

[TABLE]

In algebraic notation, the superpositions (18) and (19) can be expressed as

[TABLE]

and

[TABLE]

respectively.

The basis functions are assembled in the columns of the multiscale prolongation operator, $\bm{\mathcal{P}}$, with the dimension of $N_{fine}\times N_{coarse}$, where $N_{fine}$ and $N_{coarse}$ are the total number of fine- and coarse-scale control volumes, respectively. The part of $\bm{\mathcal{P}}$ corresponding to the matrix fine-cells is defined as

[TABLE]

Notice the three sub-blocks which represent matrix-matrix, matrix-fracture, and matrix-well coupling. Similarly, the prolongation operator for fractures can be stated as

[TABLE]

Algebraically, the complete F-AMS prolongation operator reads

[TABLE]

where $\bm{\mathcal{P}}^{wm}$ and $\bm{\mathcal{P}}^{wf}$ are set to zero, while $\bm{\mathcal{P}}^{ww}$ is the identity matrix.

Note that the prolongation operator, as described in Eq. (37), allows full flexibility in consideration of the fracture-matrix coupling in the interpolated solution, i.e. via the values in $\bm{\mathcal{P}}^{mf}$ and $\bm{\mathcal{P}}^{fm}$. This leads to the definition of four operators, differentiated by the coupling strategy they employ:

Decoupled-AMS: $\bm{\mathcal{P}}^{mf}=0$ and $\bm{\mathcal{P}}^{fm}=0$ 2. 2.

Frac-AMS: only $\bm{\mathcal{P}}^{fm}=0$. 3. 3.

Rock-AMS: only $\bm{\mathcal{P}}^{mf}=0$. 4. 4.

Coupled-AMS: $\bm{\mathcal{P}}^{mf}$ and $\bm{\mathcal{P}}^{fm}$ both non-zero.

The first option, i.e., Decoupled-AMS, constructs the most sparse $\bm{\mathcal{P}}$ and thus has an efficient setup phase. The fourth option, i.e., Coupled-AMS, can lead to more accurate multiscale simulations, however, it can severely increase the density of the multiscale operators. In such cases, one may be able to obtain a trade-off between the quality of the prolongation operator and its sparsity via truncation, followed by a rescaling of the rows to ensure partition of unity. A CPU-based study considering the overhead introduced by the density of $\bm{\mathcal{P}}$ is presented in Section 4.

In order to construct the coarse-scale system, F-AMS also needs the specification of a restriction operator, which is a map from fine- to coarse-scale (dimension $N_{coarse}\times N_{fine}$). Due to its algebraic formulation, F-AMS can accommodate multiscale finite volumes (MSFV), multiscale finite elements (MSFE) or even a hybrid multiscale finite elements and volumes restriction (MSMIX). More specifically, the MSMIX employs a FV-based restriction for part of the domain (e.g., fractures or wells), and FE for the rest (e.g., matrix rock). It is important to note that, after any MSFV stage, it is possible to construct a mass-conservative flux field for both matrix and fractures. As such, in multiphase simulations, if iterations for pressure Eq. (17) are stopped before full convergence is achieved, MSFV needs to be employed before solving transport equations. On the other hand, MSFE leads to a symmetric-positive-definite (SPD) coarse system if the fine-scale system matrix is also SPD [31, 32], and is the option used during all numerical experiments presented in this work. Note that MSMIX can be tweaked to achieve the desired compromise between MSFV and MSFE.

In the following sub-sections, first, the formulation of the local basis functions is explained. Then, the F-AMS system, and finally the simulation strategy is described in detail.

3.1 Basis function formulations

As stated before, F-AMS constructs a non-overlapping partition on the given fine-scale computational domain for both matrix rock and fracture cells, i.e. the primal-coarse grid. Then, by connecting the coarse nodes, the overlapping decomposition of the domain, i.e., dual-coarse grid, is obtained. Following the original description of the MSFV basis functions [14, 42] and its algebraic description [43, 31, 32], local basis functions are calculated for each coarse node $i$, corresponding to each dual block, by respecting the wirebasket hierarchy [44]. First the pressure in the coarse nodes (also called vertices, shown in orange in Fig. 2(b)) is set to

[TABLE]

Then the dual blocks in the neighbourhood of node $i$ are resolved, in sequence, as follows: first all the 1D dual blocks (or edges, shown in green in Fig. 2(b)), followed by the 2D (or faces, shown in blue in Fig. 2(b)) and, finally, if applicable, the 3D dual blocks (or interiors, shown in purple in Fig. 2(d)). The fact that each dual block (e.g. edge) neglects the transmissibilities to neighbouring cells belonging to blocks of inferior rank in the wirebasket hierarchy (i.e. faces and interiors), constitutes the localization assumption [31], which ensures that each basis function has a limited support.

The pressure values obtained in the manner described above, for each coarse node $i$, are assembled in column $i$ of the prolongation operator to form basis function $\Phi_{i}$.

It is important to note that by having independent fine-scale grids for each media, a matrix cell (say from a face block) can be directly connected to fracture cells belonging to dual blocks of any type (vertex, edge or face). This is an important difference from non-fractured media, where the two-point flux approximation (TPFA) stencil on structured grids ensured that any dual block would have connections only to blocks of directly superior or inferior rank in the wirebasket hierarchy (e.g., the external neighbours of a face cell are either edge or interior cells). As such, the multscale localization assumption needs to be extended to account for the connection between the two media. In the F-AMS framework, this leads to definition of basis functions which account for different degrees of coupling between the matrix and its perforating fractures.

In order to provide a compact definition of the various basis functions (matrix, fracture, well), paired with one of the four different coupling strategies considered (Decoupled-AMS, Rock-AMS, Frac-AMS, Coupled-AMS), the following “skeleton” is introduced

[TABLE]

which is solved for all basis functions $\forall\Phi^{*\bullet}\in\{\Phi^{mm},\Phi^{mf},\Phi^{mw},\Phi^{fm},\Phi^{ff},\Phi^{fw}\}$, subject to the localization assumption within each domain. Recall that $\bm{\mathcal{P}}^{wm}$ and $\bm{\mathcal{P}}^{wf}$ are zero, while $\bm{\mathcal{P}}^{ww}$ is the identity matrix. The skeleton expression (39) is based on the incompressible pressure equation, since it was found computationally efficient, even when compressibility is involved [32]. In its definition, $\text{perf}_{w}^{*}$ represents the set of (matrix or fracture) cells perforated by wells. Moreover, $\text{conn}_{mf}^{*}$ is the set of all cells with cross-media connectivities from the corresponding (matrix or fracture) domain. Finally, the $\xi(\Phi^{*\bullet})$ function gives the type of matrix-fracture coupling captured by the basis function, and will be specified separately for each strategy, as follows:

Decoupled-AMS: all basis functions have no-flow boundary conditions between the matrix and fracture domains, i.e.,

[TABLE]

which means that the fracture-matrix coupling term is completely omitted in Eq (39). The prolongation is solved algebraically, as described in Appendix A.1. Figure 3 illustrates the step-by-step procedure for a fracture and matrix basis function belonging to the 2D reservoir from Fig. 1. Note that the support of each of the interpolators is restricted to their containing medium. Finally, the Decoupled-AMS approach can be seen as applying the original AMS to separate sub-domains (i.e., matrix and fractures), having them coupled only at the coarse-scale system.

Frac-AMS: the fracture basis functions in the fracture domain, $\Phi^{ff}$, are first computed subject to no-flow conditions towards the matrix, the same as in the Decoupled-AMS (Fig. 3(d) and 3(e)), i.e. by substituting

[TABLE]

in Eq. (39). Then, the obtained values are fixed and used as Dirichlet boundary conditions while computing $\Phi^{mf}$, for which the matrix-fracture transmissibility is taken into account, i.e.,

[TABLE]

On the other hand, the matrix basis functions, $\Phi^{mm}$, are solved by setting $\Phi^{fm}=0$ as Dirichlet boundary condition, i.e.,

[TABLE]

This procedure is performed algebraically, as described in Appendix A.2. Note from Fig. 4 that, after this computation, the fracture functions have non-zero values in the matrix, while the support of the matrix basis functions is restricted to the rock domain.

Rock-AMS: First, $\Phi^{mm}$ is computed with no-flow to the fractures, as with Decoupled-AMS (Figs. 3(a)-3(c)), i.e.,

[TABLE]

Then, the values are fixed and used as Dirichlet boundaries while solving $\Phi^{fm}$, for which the fracture-matrix connections are taken into account in Eq. 39, i.e.,

[TABLE]

For the fracture functions, $\Phi^{mf}=0$ which is used as Dirichlet condition to compute $\Phi^{ff}$, i.e.,

[TABLE]

Appendix A.3 presents the algebraic procedure corresponding to this coupling strategy. Note from Fig. 5 that, in Rock-AMS, the matrix basis functions have non-zero values in the fractures, while the opposite does not hold.

Coupled-AMS: In order to preserve the two-way coupling between fractures and matrix, adjacent dual blocks of the same type are merged (e.g. fracture edges with the matrix edges they perforate), as shown in Fig. 6. To clarify, two blocks are considered adjacent if there is a non-zero transmissibility between a cell from one of them and a cell from the other.

On the resulting dual-coarse grid, computation of the basis functions follow the usual wirebasket hierarchy, with full consideration of the coupling, i.e., using

[TABLE]

for matrix and fracture cells, respectively (see Appendix A.4 for the detailed algebraic procedure). Figure 7 presents an illustration of the Coupled-AMS basis functions. Note that fracture functions have non-zero values in the matrix, and vice versa, the matrix basis functions also have non-zero values inside fractures.

It is worth mentioning that, by construction, all four coupling strategies result in basis functions which form a partition of unity.

The consideration of wells is similar to what was described in the literature for 2D problems [20], but extended here to 3D problems and integrated within the F-AMS framework. Each well is represented as a single coarse node and a well function is computed locally. The resulting values are assembled in the designated column of the prolongation operator. Note that, even for well perforations in the matrix, the corresponding well functions can have non-zero values inside fractures as well if either Coupled-AMS or Rock-AMS are employed.

As mentioned before, the basis functions have local support. However, in the Coupled-AMS case, the merging of dual blocks (Fig. 6) can increase this support substantially in cases with a high density of interconnected (or long) fractures. This, in turn, can lead to a dense prolongation operator, with a severe impact on computational performance. To overcome this, one can impose a limiting criterion on the merging of the dual blocks, for example a maximum number of fine grid cells. Alternatively, or in combination with the previous method, one can discard the non-zeros from $\bm{\mathcal{P}}$ which lie below a specified threshold, and rescale the rows accordingly to preserve the partition of unity. The latter choice is studied in detail in Section 4.

It is important to note that F-AMS basis functions are computed at the beginning of a time-dependent simulation, and adaptively updated only if the fine-scale properties change beyond a threshold value [40]. Next, the F-AMS solution algorithm is described.

3.2 The F-AMS solution algorithm

In addition to the prolongation, in order to obtain the coarse-scale pressure system, the restriction operator $\bm{\mathcal{R}}$ (i.e., map from fine to coarse scale) is now defined. As previously described, F-AMS identifies three distinct types of features in the domain, i.e., matrix, fractures, and wells, therefore, $\bm{\mathcal{R}}$ can be defined in a much more general form than in the previous studies [43, 31, 32].

The first option is to apply a FV-based restriction to all domains, i.e. the MSFV restriction operator, $\bm{\mathcal{R}}^{FV}$, where the entry at row $i$ and column $j$ is $1$ only if the fine-scale cell $j$ (either from the matrix, fractures, or wells) belongs to primal-coarse block $i$. MSFV ensures mass conservation, at the coarse-scale, thus allowing the reconstruction of a fine-scale conservative flux field. However, in previous works [31, 27], it has been found sensitive to the heterogeneity contrast in the domain. Alternatively, one can construct a Galerkin-FE-based restriction operator, as $\bm{\mathcal{R}}^{FE}=\bm{\mathcal{P}}^{T}$, traditionally called MSFE, which leads to a symmetric-positive-definite (SPD) coarse linear system, if the fine-scale system (17) is also SPD. Finally, one can consider a third option, where some of the features (e.g., matrix) are restricted according to FE, while, for the rest, FV is used. This will be referred to as the MSMIX restriction operator. Note that, for multiphase flow problems, if the pressure system is not solved to machine accuracy, a final iteration with $\bm{\mathcal{R}}^{FV}$ needs to be employed, followed by a conservative flux reconstruction stage, in order to facilitate the solution of the transport equations [36, 22].

Although F-AMS can be used as a single-sweep multiscale solver, where the approximate solution, $p^{\prime}$, is used with no iterations, previous studies have shown that an iterative procedure is needed for highly-heterogeneous reservoirs [37]. To this end, one needs to pair the F-AMS multiscale step with a fine-scale smoother, which ensures error reduction to any desired level. The F-AMS algorithm can now be compiled, as shown in Table 1.

The smoothing operator, $\bm{M}^{-1}_{S}$, approximates the inverse of the complete fine-scale linear operator, $\bm{A}^{\nu}$, via ILU(0) decomposition [45]. Note that the contrast between the matrix and fracture transmissibility values is usually severe, leading to a high condition number in Eq. (17). In such cases, F-AMS can easily be extended to include another smoothing stage, which employs iterations on the sub-block systems corresponding to each media, i.e. $\bm{A}^{mm}$ and $\bm{A}^{ff}$ from (17). A detailed study of the impact of such a smoothing stage is beyond the scope of this paper, and makes the object of future research.

Next, numerical results are presented in order to study the effect of each component on the performance of the F-AMS algorithm. Then, the scalability of F-AMS is studied, as linear solver, in a CPU benchmark, where the SAMG commercial solver [33] is used as reference.

4 Numerical Results

The aim of this section is to investigate the performance of F-AMS while simulating flow through fractured porous media. First, a 2D reservoir with heterogeneous matrix rock and a relatively complex fracture network is considered. A distance-based graph algorithm is introduced, in order to automate the fracture coarsening process. The convergence behaviour of F-AMS is studied, considering each of the four different coupling strategies introduced in the previous section. Then, the same fracture network is extruded along the Z axis and embedded in a 3D heterogeneous domain, for which simulations are performed considering different coarsening strategies, fracture-matrix conductivity ratios, fracture densities and domain sizes. CPU times are measured in detail for both the setup and solution stages in all test cases, and compared to those obtained using the industrial-grade SAMG solver [33]. Finally, the same 3D reservoir is used to investigate the effect of heterogeneous fracture conductivities, spanning several orders of magnitude.

During the upcoming simulations, special attention is given to the conductivity contrast between the matrix and the fracture domains. The transmissibility ratio $T_{ratio}$ is introduced as

[TABLE]

i.e. the ratio between the average fracture $\langle T\rangle_{frac}$ and matrix $\langle T\rangle_{rock}$ transmissibility values, respectively.

It is important to note that, in all test cases, F-AMS employs a FE restriction operator. Furthermore, the coarse-scale linear system and the basis functions in each dual block are all solved using a direct solver, based on LU decomposition, from the PETSc package [46].

For some experiments, a detailed breakdown of the CPU time spent in each stage of the F-AMS algorithm is presented. In the legends of the corresponding bar plots, the “Initialization” refers to the time spent on allocation of memory for the various data structures, the “Operators” represents the computation of basis functions and construction of the restriction and prolongation matrices (Step 4 in Table 1). Also, “Fine linsys. constr.” denotes computation of the transmissibility values and fine-scale linear system assembly. In addition, the matrix multiplications resulting in the coarse-scale system are labelled as “Coarse linsys, constr.”, while “Solution” stands for the solution of the coarse system followed by the interpolation (Step 5 in Table 1). Finally, “Smoother” accounts for Step 7.

4.1 Distance-based fracture coarsening

A fracture network can be represented as a graph, in which the fracture lines (plates) are the arcs, while the nodes are the locations at which fractures intersect. This leads to a quasi-unstructured grid, where the complexity mostly revolves around the representation of the intersections. In this work, each intersection is assigned a pressure value, which is explicitly represented in the fine-scale linear system (17) via an equation describing the conservation of flux coming from/going into the fracture control volumes it connects.

As previously described, F-AMS requires primal- and dual-coarse grids in the (quasi-unstructured) fracture domain. In order to hide this complexity from the user, this paper introduces a distance-based algorithm for the automatic coarsening of fracture networks, as described in Table 2. The network’s graph is traversed in a breadth-first order such that a distance of at least $d_{min}$ cells is guaranteed between each pair of resulting coarse nodes. As such, $d_{min}$ can be seen as a fracture coarsening factor. Note that choosing $d_{min}=\infty$ results in a single coarse node, as shown in Fig. 8(a). Moreover, Fig. 8(b) depicts the result of the coarsening algorithm for $d_{min}=20$ cells on a fairly complex fracture network.

4.2 F-AMS convergence

The fracture network from Fig. 8 was embedded into a heterogeneous (patchy) matrix rock with two pressure-constrained wells added on the West and East boundaries. This 2D test case, depicted in Fig. 9, was used to study the convergence properties of F-AMS, with the four coupling strategies presented before.

In order to test the accuracy of the basis functions as pressure interpolators, F-AMS was stopped after Step 5 of its first iteration (see Table 1). The solution is shown in Figs. 10 and 11. Using a single fracture coarse DOF leads to a poor approximation of the pressure distribution, especially for Decoupled- and Frac-AMS. Figure 10(b) depicts the results of Frac-AMS with 1 fracture coarse DOF. In this setup, F-AMS treats fractures similar to [5]. Note that, due to the large length scale of the network, having a constant interpolator for the pressure along the network results in an initial solution which lacks a lot of the fine-scale features. In contrast, the Rock-AMS and Coupled-AMS (Figs. 10(c) and 10(d)) place a lot more emphasis on the matrix basis functions and, since, in this test case, the rock heterogeneity is the main source of approximation error, their results are more accurate. It may seem unexpected that the Coupled-AMS performs slightly worse than Rock-AMS. This can be attributed to the fact that a single fracture DOF is not sufficient to accurately capture the pressure distribution in the large fracture network, especially under the localization assumption. However, when only few additional coarse DOF are added in the fracture domain (as shown in Fig. 11), the situation improves dramatically for Decoupled-AMS, Frac-AMS and Coupled-AMS. Note that Rock-AMS, on the other hand, is insensitive to this change.

For the results in Fig. 12, as well as all subsequent experiments in this paper, F-AMS was iterated until converged to a residual 2-norm of $10^{-6}$. It is clear that Rock-AMS shows a good convergence rate on this particular 2D test case, regardless of the fracture coarsening factor. Also, the other strategies reach a similar behaviour when the fracture network is enriched with only few additional coarse-scale DOF.

4.3 Basis function truncation

In order to get an idea of the performance of F-AMS on realistic fractured reservoirs, Fig. 13 introduces a 3D scenario, where the fracture network from Fig. 9 was extruded and discretized along the Z axis. Two pressure-constrained horizontal wells are placed on opposite edges of the domain boundary. Figure 13(d) shows the fine-scale pressure solution obtained on the heterogeneous (patchy) matrix permeability field shown in Fig. 13(b). Note that, even though the matrix-fracture conductivity contrast is of only two orders of magnitude, this is enough to make the pressure distribution in the fracture network insensitive to the matrix heterogeneity (see the approximately constant pressure in Fig. 13(c)).

The procedure described in Table 2 is followed to determine the fracture coarse nodes along the projection of the network on the X-Y plane. Then, the resulting coarse grids are extruded along the Z-axis uniformly, with the vertical distribution of the coarse nodes honouring the user-specified coarsening ratio (see Fig. 2).

The increased number of cells in both the matrix ($64\times 64\times 64$), as well as the fracture network ($575\times 64$), compared to the 2D case, can lead to a larger number of non-zeros for the basis functions which take into account the coupling between the two media. This can be particularly severe for Coupled-AMS (Fig. 7), since the high density of fractures can cause a large number of dual blocks to be merged. The resulting basis functions have a wider support and can potentially lead to more accurate interpolations, however, the added density to the prolongation operator will also increase the computational effort necessary to construct and solve the coarse-scale linear system (i.e., Step 5 in Table 1).

One can limit the density of $\bm{\mathcal{P}}$ by truncating basis function values below a specified threshold, $\alpha\in[0,1)$. However, in order to preserve partition of unity, the affected rows in $\bm{\mathcal{P}}$ need to be rescaled by dividing the remaining values by the row sum. Figure 14 shows the CPU time spent by F-AMS on the 3D test case, while varying the value of $\alpha$. Notice that the very restrictive value of $\alpha=10^{-1}$ leads to an increase in the number of overall iterations, because the smoother needs to compensate for the induced inaccuracy of the basis functions. However, starting with $\alpha=10^{-2}$, the convergence is no longer much affected and the algorithm gains efficiency from the reduced number of FLOPS necessary to perform $\bm{\mathcal{R}}\bm{A}^{\nu}\bm{\mathcal{P}}$ and invert the result. The truncation has the biggest impact on the Coupled-AMS strategy, which experiences a speed-up factor of 2, compared to the un-truncated case (last bar in Fig. 14, where $\alpha=0$). Figure 15 shows that, for this coupling strategy, when only a single DOF is considered for fracture network, the support of basis functions can be as big as the span of the fracture network. Also, this figure shows that after the truncation stage the locality of the basis function support can be maintained. As conclusion to this study, the subsequent experiments will use a value of $\alpha=10^{-2}$, regardless of the choice of basis function coupling strategy.

4.4 Sensitivity to the coarsening factor

The sensitivity of F-AMS to the number of coarse DOF in the fracture network, as well as the matrix coarsening ratio is studied for the 3D test case shown in Fig. 13. The coarsening factor is defined as the average number of fine cells contained in one (matrix or fracture) primal-coarse block, along each axis. Recall that, along the fracture length, this is given by the $d_{min}$ parameter in Table 2. The experiment designed for this purpose is focused on “isotropic” coarsening factors, mainly due to the point-wise nature of ILU(0), which was chosen as global smoothing stage for the implementation of F-AMS used here (Table 1).

Figure 16 shows the F-AMS CPU times obtained with three different coarsening factors for the matrix, as well as the fracture network. From this figure, the Coupled-AMS is found to automatically adapt itself to the coarsening ratio. Its convergence rate is surpassed only in cases where there is a large discrepancy between the rock and fracture coarsening ratios. However, this comes with the additional computational cost of having basis functions with wider support. Still, due to the truncation factor $\alpha=10^{-2}$, the construction and solution of coarse system for Coupled-AMS is not significantly higher than that of the alternative strategies. In addition, in each row, the optimum simulation results are obtained when fracture and matrix coarsening ratios are similar.

Based on this study and unless otherwise stated, the coarsening ratio of $8$ in each direction for both matrix and fracture media is employed in the experiments presented in the following subsections. Note that this option leads to more efficient coarse-scale systems than the alternative option of using coarsening ratios of $4$. In addition, the number of linear iterations can be significantly reduced when F-AMS is employed as preconditioner for GMRES [45].

4.5 Sensitivity to the transmissibility ratio

The next set of experiments aims to investigate the sensitivity of F-AMS to the conductivity contrast between the matrix and fractures. The transmissibility ratio $T_{ratio}$, as defined in Eq. (48), is varied over several orders of magnitude (Fig. 17) while measuring CPU times and number of linear iterations performed by the Richardson loop. A coarsening factor of $8$ was chosen for both media, based on the results from the previous subsection. As the network becomes more conductive, the influence of the matrix heterogeneity on the fracture pressure decreases (see Figs. 13(c) and 13(d)). As such, Rock-AMS exhibits a degradation in performance for higher $T_{ratio}$ (Fig. 17), while the reverse is true for Frac-AMS (Fig. 17). On the other hand, by automatically adapting to the change, the Coupled-AMS strategy remains relatively insensitive to $T_{ratio}$, as shown in Fig. 17. Finally, Decoupled-AMS requires the most number of iterations when the fracture and matrix transmissibility values are close (Fig. 17), since in this case the two-way coupling between the media is the most pronounced.

Note that the solver could not converge to the chosen tolerance, of $10^{-6}$ residual norm, when F-AMS was iterated in a Richardson’s loop for $T_{ratio}\geq 10^{6}$ and the same holds for SAMG [33]. However, as shown in Fig. 18, both methods converge successfully when employed as preconditioners for GMRES [45].

4.6 CPU benchmark study

This final subsection presents the results of a benchmark study between F-AMS and SAMG [33] on 3D heterogeneous fractured reservoirs. Both methods are employed as preconditioners to GMRES [45] and iterated until converged with a residual 2-norm below $10^{-6}$.

Unlike the Richardson loop, similar performance was observed for all experiments when F-AMS is used as preconditioner to GMRES, regardless of which coupling strategy was chosen. Therefore, the presentation of the results is restricted to Decoupled-AMS, for conciseness. Note that the Coupled-AMS strategy provides a more general framework, however, with a more complex implementation.

At each GMRES iteration, SAMG employs a single V-cycle. It is important to note that SAMG is a commercial black-box package. Thus, it is not possible to measure its CPU breakdown as accurately as for F-AMS. Instead, the time spent on its first iteration is considered as “Initialization”, while subsequent iterations were labelled as “Solution”. Finally, for both SAMG and F-AMS, the setup and construction of the operators are performed only once, at the beginning of the iteration procedure.

This study is aimed only to demonstrate the scalability of the F-AMS method. Note that a unique advantage of F-AMS over SAMG is that a fine-scale mass conservative flux field can be reconstructed after any iteration stage, once the coarse-scale system with $\bm{\mathcal{R}}^{FV}$ restriction operator is solved.

4.6.1 Transmissibility contrast

The test case from Fig. 13 is used, with different values of fracture-matrix transmissibility contrasts, i.e., $T_{ratio}$ in Eq. (48). As can be seen in Fig. 19, the number of iterations and the CPU time for both F-AMS and SAMG is insensitive to the contrast. This is a significant achievement for F-AMS, compared to [5].

4.6.2 Fracture density

What follows is a study of the scalability of F-AMS when faced with a dynamic fracture network, where the number of fracture plates is increased step by step. The 3D fracture map shown in Fig. 13(a) is considered, where the network is now created through 4 phases, as presented in Fig. 20. Note that, as new fractures are added, not only the number of DOF increases, but also the pressure variation along the network can increase. The detailed description of the CPU times obtained using F-AMS in these four cases are depicted in Fig. 21. It is clear that, by maintaining the prescribed fracture coarsening factor of $8\times 8$, F-AMS maintains virtually the same convergence rate. The slight increase in CPU time is mainly due to computation of extra fracture basis functions, as well as the construction and solution of a slightly larger coarse-scale linear system.

In consequence, by having multiple coarse-scale DOF in each fracture network, F-AMS can automatically scale with fracture length and density. This is in contrast to [5], where the use of a single fracture basis function would lead to a drastic deterioration of the multiscale convergence for test cases containing fracture networks with large length scales.

4.6.3 Domain scale

The scalability of F-AMS, benchmarked with SAMG, is investigated for heterogeneous (patchy) reservoir of increasing size. To this end, both the matrix and fracture fine-scale grid resolution is varied from $32^{3}$ matrix and $320\times 32$ fractures (smallest) up to $256^{3}$ matrix and $2117\times 256$ fracture cells (see Fig. 22). The transmissibility ratio between the two media is $T_{ratio}=10^{2}$. Figure 23 shows the obtained CPU times. F-AMS and SAMG both maintain their convergence rates and experience a similar level of scalability, in terms of CPU time, i.e., they grow linearly with the problem size. During these experiments, and as previously studied in Subsection 4.4, F-AMS was found very sensitive to the coarsening strategy used. The reported results use the optimum coarsening ratio found during repeated experiments (some of which are detailed in Table 3). Note that SAMG uses adaptive coarsening at each coarse level in its V-cycles.

4.6.4 Heterogeneous fractures

Finally, the sensitivity of the F-AMS method (benchmarked with SAMG) to strongly heterogeneous fracture properties is investigated. The permeability of of the $127$ fracture plates from Fig. 13(a) is randomly perturbed to a span $6$ orders of magnitude. Figure 24 shows that, if an appropriate coarsening ratio is chosen –in this case $6\times 6\times 8$ in the matrix and $4\times 8$ in the fractures–, then the F-AMS and SAMG performances are comparable. In addition to the coarsening ratio, multiscale methods are also sensitive to the heterogeneity contrasts (here, in both fractures and matrix). Improvements can be achieved by adapting the coarse grid geometry to follow the fracture and matrix conductivity distribution, or by enriching the prolongation operator with additional basis functions [23, 24, 25]. These are subjects of future studies.

The results of the experiments presented in this section show that the performance and scalability of F-AMS is comparable to that of SAMG. As such, even in its current proof of concept implementation, F-AMS is found a promising multiscale method for naturally fractured porous media. Note that, for realistic multiphase test cases, simulations can be further accelerated by employing only few iterations of F-AMS, followed by a flux reconstruction stage, leading to efficient approximate solutions [47].

5 Conclusions

In this paper, a novel general multiscale framework, F-AMS, was devised for efficient and accurate simulation of flow through heterogeneous porous media with embedded fractures of various length scales. For the first time in the multiscale community, the possibility to prescribe an arbitrary coarse grid in each fracture network was presented. Then, for each coarse node (from both matrix and fractures), a locally-supported basis function was defined, by considering one of the four cross-media coupling strategies (Decoupled-AMS, Frac-AMS, Rock-AMS or Coupled-AMS). All of these flexibilities allow the user to tweak the trade-off between the computational budget of the setup stage and the convergence rate.

Aligned with the EDFM approach of having independent grids for fracture and matrix [2], this work also introduced a distance-based automatic coarsening algorithm for the fracture domain. This allows the user to specify the desired (uniform) coarsening factor for the fracture domain, in similar fashion to the matrix. In addition, the effect of truncating small non-zeros from the prolongation operator was studied, in order to maintain efficiency, especially for the Coupled-AMS strategy. For all test cases considered, the truncation value of $\alpha=10^{-2}$ was found optimum.

The numerical results illustrate that F-AMS (similar as SAMG) is insensitive to the fracture-matrix conductivity contrast, and - importantly - shares the same scalability with respect to the fracture density, domain scale and heterogeneous properties. However, the performance of F-AMS was found to dramatically degrade if a sub-optimal coarsening strategy is chosen for challenging test cases. The method can be further extended to address this, e.g. by considering enhanced coarsening strategies, different choices for the smoother or employment of enriched prolongation operators [23, 24].

Another important finding of this paper was that all basis function coupling strategies perform similar when F-AMS is used as preconditioner to GMRES. This recommends the Decoupled-AMS approach for commercial reservoir simulation, due to its convenient implementation, an attractive feature for real-field applications.

In summary, it is concluded that F-AMS is an important multiscale development for flow in heterogeneous media with embedded fractures. It was shown that only few fracture coarse nodes are required to deliver good approximate pressure solutions, at the original fine-scale resolution. Future developments of F-AMS will include consideration of complex physics such as stress-dependent fracture networks and also capillarity and gravity effects for multiphase flow simulations.

Acknowledgements

Financial support of PI/ADNOC is acknowledged. Also, many thanks are due towards the members of the DARSim research group of TU Delft, for the useful discussions during the development of F-AMS method.

Appendix A Algebraic Formulation of the F-AMS prolongation operators

Consider the fine-scale system for matrix and fractures, $\bm{A}p=q$, i.e.,

[TABLE]

The permutation operator $\bm{\wp}$ containing matrix and fracture permutation block operators $\bm{\wp}^{m}$ and $\bm{\wp}^{f}$, respectively,

[TABLE]

is defined such that it reorders the linear system (55) based on the wirebasket ordering [48, 44, 31] of Internal (I), Face (F), Edge (E) and Vertex (V) for both matrix (superscript m) and fracture (superscript f) unknowns, i.e.,

[TABLE]

Note that fractures have only Face, Edge, and Vertex cells, since they are represented in a lower-dimensional space than the matrix. Also, according to the two-point flux approximation (TPFA) stencil for structured grids, $\bm{A}^{{I^{m}}{E^{m}}}$, $\bm{A}^{{E^{m}}{I^{m}}}$, $\bm{A}^{{I^{m}}{V^{m}}}$, $\bm{A}^{{V^{m}}{I^{m}}}$, $\bm{A}^{{V^{m}}{F^{m}}}$, $\bm{A}^{{F^{m}}{V^{m}}}$, $\bm{A}^{{F^{f}}{V^{f}}}$, $\bm{A}^{{V^{f}}{F^{f}}}$ are zero. More importantly, for media with embedded fractures, the coupling off-diagonal blocks $\bm{A}^{mf}$ and $\bm{A}^{fm}$ are full, i.e., each matrix cell may overlap with fracture cells of any type (F, E, or V). This is the main reason behind the consideration of the four types of basis functions, each with a different level of matrix-fracture coupling, as previously discussed in this paper.

The algebraic construction of the prolongation operator for each strategy is described next.

A.1 Decoupled-AMS

In the Decoupled-AMS prolongation operator, the matrix-fracture coupling terms are completely neglected. To this end, all off-diagonal block matrix entries (belonging to $\bm{A}^{mf}$ and $\bm{A}^{fm}$) are set to zero. In addition, similar to the AMS [31] and C-AMS [32] methods, the linear system is further simplified to account for the localization boundary condition within each medium (by neglecting connectivity between each cell and its lower-ranked neighbours in the wirebasket hierarchy). This leads to the following approximate linear system:

[TABLE]

Here, the diagonal blocks marked as $\overline{\bm{A}}$ indicate that the matrix-matrix and fracture-fracture transmissibilities, neglected due to this localization assumption ($\bm{A}^{{F^{m}}{I^{m}}}$, $\bm{A}^{{E^{m}}{F^{m}}}$, $\bm{A}^{{V^{m}}{E^{m}}}$ and $\bm{A}^{{E^{f}}{F^{f}}}$, $\bm{A}^{{V^{f}}{E^{f}}}$, respectively) have also been removed from the diagonal term. At the same time, the notation $\underline{\bm{A}}$ indicates diagonal blocks where the neglected matrix-fracture transmissibilities have been removed from the diagonal term. Finally,

[TABLE]

are the components of the coarse-scale system.

After solving for the coarse-scale pressures, $\breve{p}={\breve{\bm{A}}}^{-1}\breve{q}$, the approximate system can be inverted algebraically, due to its upper-triangular structure. Consequently, the prolongation operator, which satisfies $p^{\prime}=\bm{\mathcal{P}}\breve{p}$, reads

[TABLE]

where $\bm{I}$ is the identity matrix and the transpose operator ${\bm{\wp}^{T}}$ back-transforms the wirebasket ordering into the natural ordering. Also, $\bm{\mathcal{P}}^{E^{m}V^{m}}$, $\bm{\mathcal{P}}^{F^{m}V^{m}}$ and $\bm{\mathcal{P}}^{E^{f}V^{f}}$ are sub-blocks of the prolongation with the corresponding rows and columns given in the superscripts. For example, $\bm{\mathcal{P}}^{E^{m}V^{m}}=-{\overline{\underline{\bm{A}}}^{{E^{m}}{E^{m}}}}^{-1}\bm{A}^{{E^{m}}{V^{m}}}$ and $\bm{\mathcal{P}}^{F^{m}V^{m}}=-{\overline{\underline{\bm{A}}}^{{F^{m}}{F^{m}}}}^{-1}\bm{A}^{{F^{m}}{E^{m}}}\bm{\mathcal{P}}^{E^{m}V^{m}}$. Note that, once computed, the higher-rank sub-blocks of $\bm{\mathcal{P}}$ become boundary conditions for the values of basis functions in lower-rank cells, in accordance to the localization assumption (e.g. the values obtained for matrix edges, $\bm{\mathcal{P}}^{E_{m}V_{m}}$, are used to compute the prolongation in adjacent faces, $\bm{\mathcal{P}}^{F^{m}V^{m}}$).

A.2 Frac-AMS

The Frac-AMS approach considers the effect of the $\bm{A}^{mf}$ transmissibilities when computing basis functions. This leads to the following approximate system operator

[TABLE]

where $\bm{A}^{{F^{m}}{I^{m}}}$, $\bm{A}^{{E^{m}}{F^{m}}}$, $\bm{A}^{{V^{m}}{E^{m}}}$, $\bm{A}^{{V^{m}}{F^{f}}}$ and $\bm{A}^{{V^{m}}{E^{f}}}$ are set to zero due to localization boundary condition corresponding to Frac-AMS coupling for matrix, while, at the same time, the $\bm{A}^{{E^{f}}{F^{f}}}$, $\bm{A}^{{V^{f}}{E^{f}}}$ are zero in the fracture equations.

The Frac-AMS prolongation operator reads

[TABLE]

where the superscript f (e.g. from $\bm{\mathcal{P}}^{fV^{f}}$) corresponds to all the fracture cells, regardless of their containing dual block. Similar as in the previous case, the sub-blocks of prolongation operator $\bm{\mathcal{P}}^{F^{m}V^{m}}$, $\bm{\mathcal{P}}^{F^{m}V^{f}}$, $\bm{\mathcal{P}}^{fV^{f}}$, $\bm{\mathcal{P}}^{E^{m}V^{m}}$, $\bm{\mathcal{P}}^{E^{m}V^{f}}$, $\bm{\mathcal{P}}^{fV^{f}}$ and $\bm{\mathcal{P}}^{E^{f}V^{f}}$ represent the corresponding rows and columns given in their superscripts. For example, $\bm{\mathcal{P}}^{E^{m}V^{f}}=-{\overline{\bm{A}}^{{E^{m}}{E^{m}}}}^{-1}\bm{A}^{{E^{m}}{f}}\bm{\mathcal{P}}^{fV^{f}}$ and, specially,

[TABLE]

A.3 Rock-AMS

For Rock-AMS, the $\bm{A}^{mf}$ transmissibilities are set to zero,

[TABLE]

where, the localization boundary condition was also applied.

The Rock-AMS prolongation operator reads

[TABLE]

where the superscript m (e.g. from $\bm{\mathcal{P}}^{mV^{m}}$) corresponds to all the matrix cells, regardless of their containing dual block. The sub-blocks of the prolongation operator, i.e.,$\bm{\mathcal{P}}^{F^{m}V^{m}}$, $\bm{\mathcal{P}}^{E^{m}V^{m}}$, $\bm{\mathcal{P}}^{E^{f}V^{m}}$, $\bm{\mathcal{P}}^{mV^{m}}$ and $\bm{\mathcal{P}}^{E^{f}V^{f}}$, are defined similar to the previous cases in the sense that they represent the corresponding rows and columns given in their superscripts. For example, $\bm{\mathcal{P}}^{E^{m}V^{m}}=-{\overline{\underline{\bm{A}}}^{{E^{m}}{E^{m}}}}^{-1}\bm{A}^{{E^{m}}{V^{m}}}$ and, specially,

[TABLE]

A.4 Coupled-AMS

In the Coupled-AMS approach, all adjacent Face and Edge blocks are merged between the media (Fig. 6), i.e. $F=F^{m}\bigcup F^{f}$ and $E=E^{m}\bigcup E^{f}$. Also, let $V$ denote the set of coarse nodes, irrespective of their location. In this new setting, the $\bm{\wp}$-reordered approximate linear system is defined as

[TABLE]

where the localization boundary condition was appropriately employed.

Then, the Coupled-AMS prolongation operator reads

[TABLE]

Appendix B Coarsening ratios used during scale sensitivity test

The scale sensitivity experiment revealed that the performance of F-AMS is highly dependent on the coarsening ratios used. Only the optimum configuration was featured in the plots presented in the body of the manuscript. Table 3 lists experimental results obtained when using primal grids with more refined and more coarse blocks, respectively, for comparison purposes.