# A posteriori-based, local multilevel mesh refinement for the Darcy porous media flow problem

Soleiman Yousef

► To cite this version:

Soleiman Yousef. A posteriori-based, local multilevel mesh refinement for the Darcy porous media flow problem. Journal of Computational and Applied Mathematics, 2025, 454, pp.116162. 10.1016/j.cam.2024.116162. hal-04916336

HAL Id: hal-04916336

https://ifp.hal.science/hal-04916336v1

Submitted on 28 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# A posteriori-based, local multilevel mesh refinement for the Darcy porous media flow problem

Soleiman Yousef

IFP Energies nouvelles, 1 & 4 av. Bois Préau, 92852 Rueil-Malmaison, France

## ARTICLE INFO

**Keywords:**

Adaptive mesh refinement
A posteriori error estimates
Flux interface correction
Darcy flow problem

## ABSTRACT

In this work we develop an a posteriori-based adaptive local multilevel mesh refinement strategy for the Darcy porous media flow problem discretized with the finite volume method. Our approach is based on coupling the FIC "flux interface correction" method [1] and the equilibrated flux a posteriori error estimate [2]. The choice of the FIC method is motivated by the fact that it is well-adapted to problems discretized in conservative form since that the FIC provides corrections by balancing fluxes computed from both coarse and fin grids across the interface. The equilibrated flux a posteriori error estimate provides a guaranteed upper bound on the total error in the fluxes and takes advantage of the conservative scheme in such a way that the estimate can be easily coded and cheaply evaluated. We use this a posteriori estimate to detect automatically the zone of interest where the error is dominant in order to defined the local sub-grids. Numerous numerical experiments on practical problems illustrate the performance of our methodology. A comparison with a classical *h*-adaptive method is also carried out in this work.

## 1. Introduction

In numerical simulations, many physical models require, for a sufficiently faithful representation of the associated data, a very large numerical grid. Despite the performance of machines and the progress made in the development of high performance computing algorithms, the large size of the grid makes the computation time of a simulation excessively long. For example, geological models are generally defined on a grid containing several hundred thousand or even several millions of cells, depending on the geographical area. These models are then enriched with physical properties (porosity, permeability, relative permeability and capillary pressure laws, etc.) and thermodynamic properties. These different properties are generally known approximately hence the need to perform some preliminary simulations in order to finalize the construction of this model. For a single simulation, we then have a chain of processes making it very time-consuming. It is therefore crucial to improve the runtime performance of the simulators.

Adaptive mesh refinement AMR appears to be a judicious solution for saving, dynamically, both memory resources and computation times, while guaranteeing a given numerical accuracy of the obtained results. In all AMR techniques, the ultimate goal is to get as close as possible to the exact solution, which remains unknown. However, two main families of refinement methods can be distinguished : First, the adaptive methods: which consist in locally enriching the discretization either by the number of elements (*h*-refinement strategy [3]), or by increasing the order of consistency (analogous to the *p*-refinement strategy [4]) or by moving mesh nodes towards areas of interest (*s*-refinement strategy [5]), or by overlay finer meshes on the coarse mesh in areas of interest to form a composite mesh (*r*-refinement strategy [6]). Combining these different methods is also possible. The most popular

<footnote>E-mail address: soleiman.yousef@ifpen.fr.</footnote>
are hp-adaptive methods [7,8], but the literature also includes hr-adaptive methods [9]. Second, the local multilevel mesh refinement methods [10], where the initial computational domain is always covered by the coarsest mesh, and the problem is solved on several fine grid levels (local patches) covering only part of the domain. Such methods avoid the generation of non-conforming meshes due to refinement, and always work on meshes of good quality. Additionally, solving sub-problems on local patches allows us to have systems to solve with a moderate number of degrees of freedom. From an implementation point of view, local multilevel methods reveal great potential in parallel, as they have the advantage of distributing computations over patches in an almost individual way, while limiting the communications required on distributed architectures, which is an interesting property for implementing dynamic load-balancing mechanisms with low communication.

Several variants of local multilevel methods have been proposed in the literature that share the communication step between the coarse grid and local patches. This is an extension step, where the coarse solution is used to define the boundary conditions of the local fine problems. However, these variants differ by the correction step of the coarse solution is corrected using the fine solutions. Among these methods we find the LDC (Local defect correction) method [11], which is well suited to finite element discretization and seems to be the most popular, the FAC (Fast Adaptive Composite) method [12], which uses an intermediate composite problem between two grid levels, and finally the FIC (Flux interface correction) method [1], which consists in ensuring the conservation of flows between sub-levels, lends itself well to conservative methods, such as finite volume methods [13,14].

Whatever the refinement strategy chosen, its success depends crucially on the tool used to indicate zones to refine where the error is dominant. Using a posteriori error estimators to drive adaptive methods has always been a concern for numerical engineers. In contrast, for local multilevel refinement methods, potential refinement zones are quite often defined a priori. Nevertheless, recent work on the combination of ZZ (Zienkiewicz and Zhu) type a posteriori error estimators with the LDC method [15,16] can be found on mechanical problems with finite element discretizations. In this work, we rather explore the FIC method, which is less studied than the LDC method, but has the advantage to ensure the conservation of physical quantities which is an essential property in the design of many numerical schemes. We also aim to combine this method with appropriate error estimators.

A posteriori error estimation for the discretization of partial differential equations is a vast field of research. In this work, we are particularly interested in estimators based on equilibrated flux, cf. [17–19]. These estimators have the advantage of providing an upper bound on the error that is fully computable, with modest evaluation cost. They also satisfy all the conditions of local efficiency, robustness and sometimes asymptotic accuracy. In this respect, the work of [2] has shown that a posteriori error estimators are an extremely effective tool for steering adaptive methods in geoscience applications. This efficiency is explained by the mathematical rigor with which they are developed, unlike other more heuristic tools [20]. In this work we then combine the equilibrated flux a posteriori error estimate of [2] with the FIC “flux interface correction” method [1] in order to develop an a posteriori-based adaptive local multilevel mesh refinement strategy.

This contribution is organized as follows. Section 2 presents the model problem, a single phase Darcy flow, and the corresponding cell-centered finite volume discretization. In Section 3 we introduce the flux interface correction method and propose an a posteriori-based local multilevel mesh refinement strategy. Finally, in Section 4 we illustrate our theoretical analysis by numerical results on a simple example then on a two-phase flow model and perform a comparison with an h-refinement strategy.

## 2. Single phase Darcy flow and its finite volume discretization

In this section, we introduce a model problem and its finite volume discretization.

### 2.1. Model problem

Let $\Omega \in \mathbb{R}^d, d \ge 1$, be a bounded connected polytopal domain and denote by $\Gamma$ the Lipschitz-continuous boundary $\partial\Omega$ of $\Omega$. Let $p$ be the pressure head, we then define the Darcy velocity

$$ \mathbf{u} := -\underline{\mathbf{K}}\nabla p, \qquad (2.1) $$

with $\underline{\mathbf{K}}$ a symmetric, positive definite, and bounded tensor with values in $\mathbb{R}^{d \times d}$, which represents the permeability of the medium. Given $f \in L^2(\Omega)$, we consider the problem that consists in seeking $p : \Omega \to \mathbb{R}$ such that

$$ -\nabla \cdot (\underline{\mathbf{K}} \nabla p) = f \quad \text{in } \Omega, \qquad (2.2a) $$

$$ p = p_D \quad \text{on } \Gamma. \qquad (2.2b) $$

We denote by $H_{p_D}^1(\Omega)$ the subspace of $H^1(\Omega)$ of functions whose trace on the boundary $\Gamma$ is a given function $p_D$. Then, the weak formulation of (2.2) reads : Find $p \in H_{p_D}^1(\Omega)$ with $p = p_0 + p_D, p_0 \in H_0^1(\Omega)$, such that

$$ (\underline{\mathbf{K}}\nabla p_0, \nabla v) = (f, v) - (\underline{\mathbf{K}}\nabla p_D, \nabla v) \quad \forall v \in H_0^1(\Omega). \qquad (2.3) $$

By the Riesz representation theorem in $H_0^1(\Omega)$ for the inner product $\langle u, v \rangle := (\underline{\mathbf{K}}\nabla u, \nabla v)$, this defines $p_0 \in H_0^1(\Omega)$ uniquely. Then, with $p = p_0 + p_D$, this gives us the existence and uniqueness of the weak solution of the problem (2.3). In what follows, we use the energy norm: for $v \in [L^2(\Omega)]^d$,

$$ \|v\|_{\underline{\mathbf{K}}^{-\frac{1}{2}}; \Omega} := \left\{ \int_{\Omega} \left| \underline{\mathbf{K}}^{-\frac{1}{2}}(x)v(x) \right|^2 dx \right\}^{\frac{1}{2}}. $$
## 2.2. Mesh description

Let $\mathcal{T}_h$ be a mesh defined as a finite collection of control volume $K$ such that $\cup_{K \in \mathcal{T}_h} K = \bar{\Omega}$. We denote by $|K|$ the $d$-dimensional Lebesgue measure of $|K|$, and we assume that $|K| > 0$. We denote by $\partial K = \bar{K} \setminus K$ the boundary of $K$ and by $\mathbf{x}_K$ the barycenter of $K$. Let $\mathcal{E}_h$ a finite family of disjoint subsets of $\bar{\Omega}$ (the faces of the mesh). We assume that for any $K \in \mathcal{T}_h$, there exists a subset $\mathcal{E}_K$ of $\mathcal{E}_h$ such that $\partial K = \cup_{\sigma \in \mathcal{E}_K} \bar{\sigma}$. Then, for any $\sigma \in \mathcal{E}_h$, we denote by $\mathcal{T}_\sigma = \{K \in \mathcal{T}_h | \sigma \in \mathcal{E}_K\}$ (the neighbors of $\sigma$), and assume that either $\mathcal{T}_\sigma$ has exactly one element, and then $\sigma \subset \partial\Omega$ (the set of these faces, called boundary faces, is denoted $\mathcal{E}_h^{\text{ext}}$), or $\mathcal{T}_\sigma$ has exactly two elements (the set of these faces, called interior faces, is denoted $\mathcal{E}_h^{\text{int}}$). Suppose that each boundary face lies entirely in $\Gamma$ and denote by $\mathcal{E}_h^D$ the ensemble of these faces. For all $K \in \mathcal{T}_h$ and all $\sigma \in \mathcal{E}_K$, we denote by $\mathbf{n}_{K,\sigma}$ the unit normal vector to $\sigma$ outward to $K$. For any $\sigma \in \mathcal{E}_h$, we denote by $\mathbf{x}_\sigma$ its barycenter. In what follows we let, for all $K \in \mathcal{T}_h$ and all $\sigma \in \mathcal{E}_K$, $d_{K,\sigma} := \text{dist}(\mathbf{x}_K, \mathbf{x}_\sigma)$ and $\mathbf{n}_K := (\mathbf{n}_{K,\sigma}) \cdot \mathbf{n}_{K,\sigma}$.

## 2.3. Cell-centered finite volume scheme

The cell centered finite volume scheme [21] is a two-point flux approximation scheme (TPFA scheme), where the unknowns of the model are discretized using one value per cell. Then, system (2.2) is discretized as follows: find $p_h := (p_K)_{K \in \mathcal{T}_h} \in \mathbb{P}_0(\mathcal{T}_h)$ such that

$$ \sum_{\sigma \in \mathcal{E}_K} \mathcal{T}_{K,\sigma}(p_h) = (f, 1)_K, \quad \forall K \in \mathcal{T}_h. \qquad (2.4) $$

where $\mathcal{T}_{K,\sigma}(p_h)$ is the two-point finite volume approximation of the normal component of the average velocity $\mathbf{u}$ on $\sigma$:

$$ \mathcal{T}_{K,\sigma}(p_h) := \begin{cases} |\sigma|\alpha_K (p_K - p_\sigma), & \text{if } \sigma \in \mathcal{E}_h^D, \\ |\sigma|\frac{\alpha_K \alpha_L}{\alpha_K + \alpha_L} (p_K - p_L), & \text{if } \sigma = K|L \in \mathcal{E}_h^{\text{int}}, \end{cases} \qquad (2.5) $$

with $\alpha_K := \frac{\mathbf{n}_K^\sigma}{d_{K,\sigma}}$, $\forall K \in \mathcal{T}_h$.

**Remark 2.1.** Note that, by construction of the cell centered finite volume scheme, if we combine the flux expression (2.5) with the flux conservation propriety: $\mathcal{T}_{K,\sigma}(p_h) + \mathcal{T}_{L,\sigma}(p_h) = 0, \forall \sigma = K|L \in \mathcal{E}_h^{\text{int}}$, then we can evaluate a value of the pressure per face

$$ p_{\sigma} = \frac{\alpha_{K} p_{K} + \alpha_{L} p_{L}}{\alpha_{K} + \alpha_{L}}, \quad \forall \sigma \in \mathcal{E}_{h}^{\text{int}}. $$

For the exterior faces, we use the boundary condition (2.2b) together with (2.5) to write

$$ p_{\sigma} = p_{D}(\mathbf{x}_{\sigma}) \quad \forall \sigma \in \mathcal{E}_{h}^{D}. $$

# 3. Coupling the flux interface correction method and the equilibrated flux a posteriori estimate

In this section, our goal is to propose an a posteriori-based local multilevel mesh refinement strategy. Accordingly, we first introduce the Flux Interface Correction method. After, we explain how to use the equilibrated flux a posteriori error estimate in order to guide the FIC method. Finally, we discuss about patch generations algorithm which take place as a crucial step for the efficiency of the proposed refinement algorithm.

## 3.1. The flux interface correction method

The flux interface correction method [1] is a multilevel mesh refinement technique that consists in generating, from an initial coarse grid, several local levels of finer subgrids. On each subgrid an independent local problem is solved which is the same main problem described in Section 2.1. We rely the different levels of subgrids through prolongation and restriction operators.

### 3.1.1. The prolongation step

In order to solve a local problem on a fine subgrid, boundary conditions need to be prescribed. We impose inhomogeneous Dirichlet boundary conditions to the pressure as follows:

* If a part of the boundaries coincide with boundary of the initial coarse grid we keep the boundary conditions given by the physical problem data.

* Otherwise, we impose new inhomogeneous boundary conditions via a prolongation operator $P_{l-1,l}$: Based on the description of the cell centered finite volume scheme (2.4)–(2.5) a value of the pressure at each face of the boundary $p_{\sigma l}$, $\sigma^l \in \gamma^l$ must be imposed. We use the advantage of the finite volume method that gives us a value of the pressure at each face of the coarser grid $p_{\sigma l-1}$, $\sigma^{l-1} \in \mathcal{E}_h^{l-1}$, see Remark 2.1. Then, the prolongation operator $P_{l-1,l}$ can be a linear interpolation operator

$$ p_{\sigma l} = \frac{p_{K^{l-1}} + p_{\sigma l-1}}{2}. $$
### 3.1.2. The restriction step

Before we resolve the problem on the coarse grid at level $l-1$ (iteration $k+1$) we need to add some corrective terms issued from the solution obtained on the fine local subgrids at level $l$ (iteration $k$). We first introduce the restriction operator $R^{l,l-1}$

$$R_{l,l-1}^k(p_{hl})(K^{l-1}) := \sum_{\sigma^l \in \partial K^{l-1}; K^l \in K^{l-1}} F_{K^l, \sigma^l}^k(p_{hl})$$

Following [1], we then define the flux correction residual on each coarse control volume $K^{l-1}$ by:

$$r_{l-1}^{k+1}(p_{hl})(K^{l-1}) := \frac{1}{|K^{l-1}|} \left\{ R_{l,l-1}^k(p_{hl})(K^{l-1}) - \sum_{\sigma^{l-1} \in \partial K^{l-1}} F_{K^{l-1}, \sigma^{l-1}}^k(p_{hl}) \right\}$$

This correction terms are then added to the right-hand side of the coarse grid problem :

$$f_{l-1}^{k+1} = f_{l-1}^{k} + \sum_{K^{l-1} \in \mathcal{T}_{hl}^{l-1}} \mathbb{1}_{K^{l-1}} r_{l-1}^{k+1}(p_{hl})(K^{l-1}), \quad k \ge 0,$$

where $\mathbb{1}_{K^{l-1}}$ denotes the characteristic function of the control volume $K^{l-1}$.

### 3.2. Equilibrated flux a posteriori estimates

In order to identify the different zones where it is necessary to solve on a finer subgrid we make use of the equilibrated flux a posteriori error estimates. As shown in [2], the development of the a posteriori error estimate, for single-phase Darcy flow, consists in considering an $H^1$-conforming potential reconstruction $s_h$ and an $H(\text{div}, \Omega)$-conforming flux reconstruction $u_h$. This two ingredients, essential for deriving a posteriori error estimates, are fictitious in the sense that they do not need to be constructed in practice in order to evaluate the resulting a posteriori error estimate, as we will just need for the evaluation some degrees of freedom related to the finite volume discrete solution and the discrete flux, see [2, Section 3].

Following [2, Theorem 3.12], we obtain a guaranteed a posteriori error estimate of the form

$$\|\mathbf{u} - \mathbf{u}_h\|_{K^{-\frac{1}{2}}; \Omega} \le \left\{ \sum_{K \in \mathcal{T}_h} \eta_K^2 \right\}^{1/2},$$

where the local estimator $\eta_K$ merely consists in some local matrix-vector multiplications, where, on each mesh element, the matrices are either directly inherited from the given numerical method, or easily constructed from the element geometry, while the vectors are the flux and potential values on the given element. The use of this kind of a posteriori error estimates, for a large family of spatial discretizations, has already proved to be efficient in practice by providing a valuable decision tool to drive adaptive algorithm. For more details, we refer to [2,20,22].

### 3.3. Patch generation algorithms

A crucial step in the proposed methodology rely on having a tool which allows, in an efficient way, to generate local patches from the set of cells flagged to be refined by the a posteriori estimate. In this work, we are based on the Berger-Rigoustos algorithm proposed initially in [23] then analyzed and improved by Liven in [24]. We restrict the algorithm with the objective to define rectangular boxes that cover the flagged cells. The generated patches are non-overlapping and respect two main criteria : a minimum size and a maximum size which are a user-defined numbers fixed generally in with respect to the coarse mesh size. For more details about the Berger-Rigoustos algorithm see [23,24].

<figure><img src="image_4.png" /><figcaption>Fig. 1. Cell to refine (left), patches (right).</figcaption></figure>

In Fig. 1, we test the patch generation algorithm on an example where the flagged cells are in separated zones. We observe that the algorithm is able to define, in an efficient way, patches that covers exactly the flagged cells. Fig. 2, shows more complicated example where we have a continuous flagged zone. Then, we have many configurations of the constructed patches by tuning the max size of the generated patches. A big max size number may include extra non-flagged cells while a small max size number may generate a lot of patches. From a parallel point of view, having more small patches could be more efficient than having less big patches that include unnecessary cells.
<figure><img src="image_5.png" /><figcaption>Fig. 2. Cell to refine (left), patches : size max 4 (middle), size max 6 (right).</figcaption></figure>

### 3.4. An a posteriori-based, multilevel local mesh refinement algorithm

We are now able to propose our a posteriori-based, multilevel local mesh refinement algorithm. From an initial problem like (2.2) we compute an initial coarse solution then using the a posteriori error estimate of Section 3.2 we flag the cells where the error is dominant in respect with a refinement criteria. The flagged cells is then used as an input of the patch generation algorithms of Section 3.3. Once the patches are generated we refine them locally before moving to the next step which consists in applying the flux interface correction method of Section 3.1. The FIC method is applied in an iterative way until we observe no significant improvement in the solution measured by computing the $L^2$-norm of the difference between the last two solutions or until we reach a user-defined maximum number of iterations. This a posteriori-based multilevel local mesh refinement algorithm is illustrated in Fig. 3.

<figure><img src="image_6.png" /><figcaption>Fig. 3. An a posteriori-based, multilevel local mesh refinement algorithm.</figcaption></figure>

Note that with the ability to provide a guaranteed upper bounds for a large family of spatial discretization (see Section 3.2), the Algorithm 3 is very generic and can be applied with different discretization methods and different choices of the local multilevel refinement strategy. For example, instead of choosing Finite Volume as the discretization method and Flux Interface Correction (FIC) as the local multilevel mesh refinement strategy, a Finite Element method accompanied by a Local Defect Correction (LDC) method can fit perfectly into our proposed algorithm. In this case, the algorithm is simply adapted by using the finite element method to solve the global coarse problem and the various local fine problems and by employing the restriction step of the LDC instead of the FIC method.

**Remark 3.1 (Unsteady problem).** Algorithm 3 is illustrated for a single resolution. For a non-stationary problems, we need then to apply it at each time step. This will lead us to recompute the estimators, identify new flagged cells and recompute the patches distribution at each time step. As discussed is Section 3.2 and shown in [2,20] the cost of the estimators evaluation is low compared with the total computation CPU time. Moreover, the patch generation algorithm is near instantaneous so that the cost of repeating the proposed a posteriori-based, multilevel local mesh refinement algorithm will mainly depends on the number of iterations of the FIC method. In practice, we observe that we can simply take $k = 1$ in Fig. 3.
<figure><img src="image_7.png" /><figcaption>Fig. 4. Initial solution (left), error estimate (right).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Generated patch (left), patch solution (right).</figcaption></figure>

# 4. Numerical experiments

In this section, we validate our proposed adaptive algorithm on a simple example of a single phase flow problem before testing on a more practical example issued from a two-phase flow model.

## 4.1. Single phase flow model

For the first test, we consider a simple two-dimensional problem of form (2.2) in the domain $\Omega = (0, 1) \times (0, 1)$ with $u_D = 0$, $K = \text{Id}$, and the load term function $f$ given such that the (regular) exact solution reads $p(x) = 2^{4\alpha} x^\alpha (1-x)^\alpha y^\alpha (1-y)^\alpha$ with $\alpha = 100$.

The problem is discretized with the cell-centered finite volume scheme of Section 2.3. The a posteriori error estimate in this case indicates one zone which is centered at the middle of the domain, see the right part of Fig. 4. This leads the adaptive algorithm 3 to designate one local patch. The solution on this local patch, needed to correct the initial solution, is shown in the right part of Fig. 5. Note that in the right part of Fig. 5 the zone outsize the centered local patch is not considered in the fine resolution so that uniform refinement process is applied just on a patch of $7 \times 7$ elements. With this approach we then avoid the generation of nonconforming meshes due to the refinement and then to solve the problem on the local patch we still apply the cell-centered finite volume scheme.

## 4.2. Two-phase flow model

We consider an immiscible two-phase model where the phases are water and gas. Each of them is composed of only one component, water component and carbon dioxide component, respectively. As we have an immiscible model we will use the same index to represent the phase and its component. Consequently, we use the indices w,g to represent the water and the gas, respectively.

We have the following relation between the phase pressures

$$P_g = P_w + P_{g,w} \qquad (4.1)$$
where $P_{c,g,w}$ is the capillary pressure. Consequently, the unknowns of the model are the water pressure $P_w$ and the saturations $S_p$, with $p \in \{g, w\}$. The system of governing equations is given by

$$ 
\begin{align} 
\phi \partial_t (S_w) + \nabla \cdot \left( \rho_w \frac{k_{r,w}}{\mu_w} \mathbf{v}_w \right) &= q_w, \\ 
\phi \partial_t (S_g) + \nabla \cdot \left( \rho_g \frac{k_{r,g}}{\mu_g} \mathbf{v}_g \right) &= q_g, \\ 
S_w + S_g &= 1. 
\end{align} 
\tag{4.2} $$

with no-flow boundary conditions. For all $p \in \{g, w\}$, $v_p := v(P_p)$ is the phase velocity given by Darcy's law,

$$ v(P_p) = -\underline{\mathbf{K}} (\nabla P_p - \rho_p g) $$

where $g = -g\nabla z$ denotes the gravity vector acting in the negative $z$ direction and $g$ its Euclidean norm.

### 4.2.1. Numerical setting

We consider the injection of CO₂ into a reservoir initially saturated with water. The spatial domain is $\Omega = (0, 1010)m \times (0, 1010)m$ discretized by a grid of 51 × 51 uniform cells and the process is simulated for $t_F = 5$ days with an initial time step $\tau^0 = 1 \times 10^3$ s. We adopt a two-spot pattern, where we have one injection well located at the lower left corner with pressure fixed to $P_{inj} = 5.6 \times 10^7$ Pa and a production well at the upper right corner with pressure fixed to $P_{pro} = 2.7 \times 10^7$ Pa. The required physical properties are chosen as follows:

* the porosity $\phi = 0.1$,

* The permeability field $\underline{\mathbf{K}} = 1 \times 10^{-13} m^2$,

* the water molar density $\rho_w = 1025$ kg m⁻³ and the gas molar density $\rho_g = 1.89$ kg m⁻³,

* the water phase viscosity $\mu_w = 1 \times 10^{-3}$ Pa s and the gas phase viscosity $\mu_g = 1.42 \times 10^{-5}$ Pa s,

* the relative permeability

$$ k_{r,p}(S_p) = \begin{cases} 1 & \text{if } S_p \ge 1, \\ \frac{S_p - S_p^{\text{res}}}{1 - S_p^{\text{res}}} & \text{if } S_p^{\text{res}} < S_p < 1, \\ 0 & \text{if } S_p \le S_p^{\text{res}}, \end{cases} $$

where the residual saturations are respectively given by $S_w^{\text{res}} = 0.1$ and $S_g^{\text{res}} = 0.1$.

For the phase pressures, we choose $P_w$ as the reference pressure $P$, i.e. $P_{c,w} = 0$, and follow the Brooks–Corey model [25] for the gas phase capillary pressure law, i.e.,

$$ P_g = P + P_{cg}(S_g), \quad P_{cg}(S_g) = P_e \cdot (S_e)^m, \quad S_e = 1 - \frac{S_g - S_g^{\text{res}}}{1 - S_w^{\text{res}} - S_g^{\text{res}}}, $$

with $P_e = 8.73 \times 10^5$ Pa, $m = -\frac{1}{2.89}$, and $S_w^{\text{res}}, S_g^{\text{res}}$ the residual saturations defined previously.

We consider a cell-centered finite volume scheme in space and the backward Euler scheme in time with the Newton method for the resulting nonlinear systems. Then, we use the BICGStab method to solve the outcome linear system at each Newton iteration. For more details about the numerical resolution see [26].

### 4.2.2. A posteriori-based, multilevel local mesh refinement

Considering the two-phase flow model, we deal with an unsteady nonlinear coupled problems, then for the refinement process, we need to apply Algorithm 3 at each time step of the resolution. For the error estimate, as mentioned in Section 3.2 a guaranteed equilibrated flux a posteriori error estimate is already available for the two-phase flow model [2,20,22]. Here, we rely on our recent work of [2]. The set of cells flagged by the error estimate is then used to generate locales patches as described in Section 3.3. Moreover, for the FIC correction method, we follow the prolongation step with Dirichlet boundary conditions and the restriction step as described in Section 3.1.

Fig. 6 depicts, at a fixed time step, the evolution of the approximate water saturation (left) and the corresponding error estimator (right). We see that the estimator detects well the error following the saturation front. The cells where the error is dominant is then flagged and added to the zone to refine. Then the patch generation algorithm, in this case, gives eight locales patches tagged in the left part of Fig. 7 with different colors in order to visually identify them. Afterwards, the two phase flow problem is solved on each of the eight local patches refined in an uniform way with two level of refinement. Note that, to solve on the local patches we still apply, for the space discretization, the two-point flux approximation scheme as our algorithm avoids the generation of nonconforming meshes due to the refinement. Algorithm 3, is then repeated at each time step. Note that, since the saturation front is moving during the time simulation, the zone to refine is then updated at each time step. Therefore, we need to update the generation of the local patches at each time step also. However, as mentioned in Section 3.4 the cost of the estimators evaluation is low compared with the total computation CPU time and the patch generation algorithm is near instantaneous.
<figure><img src="image_9.png" /><figcaption>Fig. 6. Approximate solution, error estimate.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Generated patches and patches solutions.</figcaption></figure>

The right part of Fig. 7 shows, at the same time step, the resulting local solutions of the two-phase flow problem projected on each of the initial patches inside the global mesh in such a way we can compare the result in the right part of Fig. 7 with the approximate solution in the left part of Fig. 6 for the local refinement zone.

As explained previously and detailed in Section 3.4, at each time step, the different local solutions are used to correct the global coarse solution following the restriction step of Section 3.1. Note that, for this practical example we observe that one iteration of the FIC method is sufficient to have a significant improvement in the solution (corresponding to $k = 1$ in Algorithm 3).

### 4.2.3. A comparison with an h-refinement strategy

We consider here an h-refinement strategy to adapt the mesh. For the decision tool, as for the multilevel mesh refinement, we rely on the equilibrated flux a posteriori error estimate of Section 3.2. We propose then the following mesh adaptive algorithm:

**Algorithm 4.1 (h-refinement algorithm).**

Fix the fractions of cells to refine, $\zeta_{ref}$, and to derefine, $\zeta_{deref}$

while $t^n \le t_F$ do {Time loop}

Solve the system

Compute the local spatial estimators.

Refine the cells $K \in \mathcal{T}_h^n$ such that $\eta_K^n \ge \zeta_{ref} \max_{L \in \mathcal{T}_h^n} \{\eta_L^n\}$.

Derefine the cells $K \in \mathcal{T}_h^n$ such that $\eta_K^n \le \zeta_{deref} \max_{L \in \mathcal{T}_h^n} \{\eta_L^n\}$.

Update data

end while
<figure><img src="image_11.png" /><figcaption>Fig. 8. *h*-adaptivity.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Patch-based refinement v.s. *h*-adaptive refinement.</figcaption></figure>

By applying Algorithm 4.1, we generate non-conforming meshes so the two-point flux approximation scheme will be non-consistent. To overcome this problem, we choose as a space discretization method, when considering the *h*-refinement strategy, the multi-point flux approximation L-method which is more expensive than the TPFA but seem to be a good choice for adaptive mesh refinement, see [27].

Fig. 8 shows the evolution of the approximate water saturation at two different simulation time while using algorithm 4.1. We remark that the refinement follows the front of the water saturation as time evolves, and then the derefinement process is effected in the zones abandoned by the water front.

The efficiency of the two a posteriori-based AMR method (the multilevel local mesh refinement and the *h*-refinement strategy) can be appreciated in Fig. 9 which illustrates the cumulated rate of water recovered from the production well during the simulation; we compare here the result on a uniform fine mesh, a uniform coarse mesh, an adapted mesh with the a posteriori-based multilevel local mesh refinement strategy and an adapted mesh with the a posteriori-based *h*-refinement strategy. We observe that applying the refinement strategy gives us a very good precision whether we opt for the multilevel local mesh refinement or the *h*-refinement strategy.

Note, however, that when the *h*-refinement strategy is applied, it was required to choose a more advanced discretization scheme to deal with the nonconforming meshes. Something that is not necessary with the local multilevel mesh refinement method since it avoids the generation of nonconforming meshes and always works on meshes of good quality. This gives our proposed a posteriori-based, local multilevel mesh refinement strategy a major advantage from a practical point of view as it allows us in this work to stay with the cell-centered finite volume scheme which is very popular in reservoir engineering because it is easy to code and it meets the industrial constraints of robustness and low computational cost. Conversely, when using the *h*-adaptive method as we generate nonconforming meshes then the cell-centered finite volume scheme is no longer consistent and we needed to code another
more advanced scheme adapted to nonconforming meshes (in this case we used the multi-point flux approximation L-method). The
advantage of the local multilevel method in parallel computation still to be investigated in a future work.

5. Conclusion

This paper presents an a posteriori-based adaptive local multilevel mesh refinement strategy. The a posteriori error estimate is used to efficiently identify the different zones where it is necessary to solve on a finer subgrid. On each subgrid an independent local problem is solved which is the same main problem where we rely the different levels of subgrids through prolongation and restriction operators. The efficiency of the method in providing an accurate solution, while working on a smaller grid size compared to a uniform fine mesh, is illustrated and a comparison with a classical *h*-refinement strategy is also carried out. As mentioned in Section 4, the main advantage of our method over the *h*-adaptive method is that we always work with matching meshes (the coarse grid and the uniformly refined local patches). This gives the possibility to use a simple discretization scheme with low computational cost. In contrast, with the *h*-adaptive method we generate nonconforming meshes and this leads to the use of a more “expensive” scheme adapted to nonconforming meshes. As a perspective, we still need to study the potential of the proposed method in parallel, in particular on load balancing mechanisms.

Data availability

The data that has been used is confidential.

References

[1] P. Angot, J.P. Caltagirone, K. Khadra, Une méthode adaptative de raffinement local : la correction du flux à l'interface, C. R. Math. Acad. Sci. Paris 315 (1) (1992) 739-745.

[2] M. Vohralík, S. Yousef, A simple a posteriori estimate on general polytopal meshes with applications to complex porous media flows, Comput. Methods Appl. Mech. Engrg. 331 (2018) 728-760.

[3] L. Demkowicz, P. Devloo, J. Oden, On an h-type mesh-refinement strategy based on minimization of interpolation errors, Comput. Methods Appl. Mech. Engrg. 53 (1) (1985) 67-89.

[4] R.R. Salagame, A.D. Belegundu, A simple p-adaptive refinement procedure for structural shape optimization, Finite Elem. Anal. Des. 24 (3) (1997) 133-155.

[5] Z. Yue, D.H. Robbins Jr., Adaptive superposition of finite element meshes in elastodynamic problems, Internat. J. Numer. Methods Engrg. 63 (11) (2005) 1604-1635.

[6] W. Cao, W. Huang, R.D. Russell, Anr-adaptive finite element method based upon moving mesh PDEs, J. Comput. Phys. 149 (2) (1999) 221-244.

[7] I. Babuška, B. Guo, The h, p and h-p version of the finite element method; basis theory and applications, Adv. Eng. Softw. 15 (3) (1992) 159-174.

[8] P. Daniel, A. Ern, I. Smears, M. Vohralík, An adaptive hp-refinement strategy with computable guaranteed bound on the error reduction factor, Comput. Math. Appl. 76 (5) (2018) 967-983.

[9] M. Piggott, C. Pain, G. Gorman, P. Power, A. Goddard, h, r, and hr adaptivity with applications in numerical ocean modelling, Ocean Model. 10 (1) (2005) 95-113, The Second International Workshop on Unstructured Mesh Numerical Modelling of Coastal, Shelf and Ocean Flows.

[10] A. Brandt, Multi-level adaptive solutions to boundary-value problems, Math. Comp. 31 (138) (1977) 333-390.

[11] W. Hackbusch, Local defect correction method and domain decomposition techniques, in: Defect Correction Methods: Theory and Applications, Springer Vienna, Vienna, 1984, pp. 89-113.

[12] S. McCormick, Fast adaptive composite grid (FAC) methods: Theory for the variational case, in: Defect Correction Methods: Theory and Applications, Springer Vienna, Vienna, 1984, pp. 115-121.

[13] O. Pironneau, A. Lozinski, A. Perronnet, F. Hecht, Numerical zoom for multiscale problems with an application to flows through porous media, Discrete Contin. Dyn. Syst. 23 (1-2) (2009) 265-280.

[14] I. Ramière, P. Angot, M. Belliard, A general fictitious domain method with immersed jumps and multilevel nested structured meshes, J. Comput. Phys. 225 (2) (2007) 1347-1387.

[15] L. Barbié, I. Ramière, F. Lebon, Strategies involving the local defect correction multi-level refinement method for solving three-dimensional linear elastic problems, Comput. Struct. 130 (2014) 73-90.

[16] H. Liu, I. Ramière, F. Lebon, On the coupling of local multilevel mesh refinement and ZZ methods for unilateral frictional contact problems in elastostatics, Comput. Methods Appl. Mech. Engrg. 323 (2017) 1-26.

[17] W. Prager, J.L. Synge, Approximations in elasticity based on the concept of function space, Quart. Appl. Math. 5 (1947) 241-269.

[18] P. Ladevèze, J.-P. Pelle, Mastering Calculations in Linear and Nonlinear Mechanics, in: Mechanical Engineering Series, Springer-Verlag, New York, 2005, p. xii+413, Translated from the 2001 French original by Theofanis Strouboulis.

[19] A. Ern, M. Vohralík, A posteriori error estimation based on potential and flux reconstruction for the heat equation, SIAM J. Numer. Anal. 48 (1) (2010) 198-223.

[20] J.-M. Gratien, O. Ricois, S. Yousef, Reservoir simulator runtime enhancement based on a posteriori error estimation techniques, Oil Gas Sci. Technol. - Rev. IFP Energies Nouvelles 71 (5) (2016) 59.

[21] R. Eymard, T. Gallouët, R. Herbin, Finite volume methods, in: Handbook of Numerical Analysis, Vol. VII, North-Holland, Amsterdam, 2000, pp. 713-1020.

[22] G. Mallik, M. Vohralík, S. Yousef, Goal-oriented a posteriori error estimation for conforming and nonconforming approximations with inexact solvers, J. Comput. Appl. Math. 366 (2020) 112367.

[23] M. Berger, I. Rigoutsos, An algorithm for point clustering and grid generation, IEEE Trans. Syst. Man Cybern. 21 (5) (1991) 1278-1286.

[24] O.E. Livne, Minimum and Maximum Patch Size Clustering on a Single Refinement Level, Tech. rep. UUSCI-2006-002, Scientific Computing and Imaging Institute, University of Utah, 2006.

[25] R.J. Brooks, A.T. Corey, Hydraulic Properties of Porous Media, Hydrology Paper 3, Colorado State University, Fort Collins, 1964.

[26] J.-M. Gratien, O. Ricois, S. Yousef, Reservoir simulator runtime enhancement based on a posteriori error estimation techniques, Oil Gas Sci. Technol. Rev. IFP Energies Nouvelles 71 (5) (2016) 59.

[27] M. Wolff, Y. Cao, B. Flemisch, R. Helmig, B. Wohlmuth, Multi-point flux approximation L-method in 3D: numerical convergence and application to two-phase flow through porous media, Radon Ser. Comput. Appl. Math. de Gruyter 12 (2013) 39-80.