# Numerical validation of an adaptive model for the determination of nonlinear-flow regions in highly heterogeneous porous media

Alessio Fumagalli, Francesco Saverio Patacchini

► To cite this version:

Alessio Fumagalli, Francesco Saverio Patacchini. Numerical validation of an adaptive model for the determination of nonlinear-flow regions in highly heterogeneous porous media. Mathematics and Computers in Simulation, 2025, 230, pp.306-322. 10.1016/j.matcom.2024.10.024. hal-04929252

HAL Id: hal-04929252

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

Submitted on 4 Feb 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>

ELSEVIER

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

Original articles

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

# Numerical validation of an adaptive model for the determination of nonlinear-flow regions in highly heterogeneous porous media

Alessio Fumagalli<sup>a,*</sup>, Francesco Saverio Patacchini<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Porous medium

Heterogeneous medium

Adaptive constitutive law

Numerical validation

## ABSTRACT

An adaptive model for the description of flows in highly heterogeneous porous media is developed in Fumagalli and Patacchini (2022, 2023). There, depending on the magnitude of the fluid's velocity, the constitutive law linking velocity and pressure gradient is selected between two possible options, one better adapted to slow motion and the other to fast motion. We propose here to validate further this adaptive approach by means of more extensive numerical experiments, including a three-dimensional case, as well as to use such approach to determine a partition of the domain into slow- and fast-flow regions.

## 1. Introduction

In the context of flow modeling in porous media, Darcy's law is omnipresent. It is a linear relationship between velocity and pressure gradient, and allows for tensor permeabilities. Its validity is supported both theoretically, via a derivation from Navier-Stokes equations, and empirically, through countless experiments. When velocities deviate from those typically found in porous media, as a result of particularly heterogeneous permeabilities, Darcy's law can, however, be challenged. This was noted experimentally in many studies [1-3]. Indeed, in highly permeable areas of a porous medium, such as in fractures, macropores or conduits, the velocity may get high enough to generate nonnegligible inertial and frictional effects which translate into a nonlinear relation with the pressure gradient. Experimentally, this deviation is noticed as an overestimation of the velocity by Darcy's law compared to that empirically measured. Indeed, from an energetic viewpoint, the effect of the nonlinearity is to penalize high velocities by giving them a high energetic cost. The nonlinearity can also be retrieved theoretically, from averaging and homogenization methods, without, however, providing an explicit formulation of this nonlinearity [4-8]. Instead, these methods yield the nonlinearity as a perturbation, or correction, of Darcy's law which needs to be fitted empirically. One possible expression for the nonlinearity, which we assume in the present work, is polynomial, meaning that the pressure gradient is a polynomial function of the velocity. This approach includes the widely used Darcy-Forchheimer law, which is a quadratic correction.

Numerically, the nonlinearity thus introduced increases precision but slows down simulations. Nevertheless, the portions of the porous medium where the nonlinear correction has a significant impact are typically small compared to the size of the whole domain. Therefore, the use of the nonlinear perturbation of Darcy's law can be limited to a small subdomain, while the classical Darcy law can be kept in most of the medium. As a result, a-priori knowledge of where these nonlinear regions lie should help improve accuracy without significant computational cost. Indeed, if a domain partition can be devised subdividing the medium into a slow region and a fast region, then a domain-decomposition algorithm [9-11] can be run over this partition where Darcy's law is solved in the slow region and a nonlinear law, such as the Darcy-Forchheimer law, is solved in the fast region. We expect this

<footnote><sup>a</sup> MOX Laboratory, Department of Mathematics, Politecnico di Milano, p.za Leonardo da Vinci 32, Milano 20133, Italy</footnote>

<footnote><sup>b</sup> IFP Energies nouvelles, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: alessio.fumagalli@polimi.it (A. Fumagalli).</footnote>
approach to be more precise than solving Darcy's law globally and faster than solving the nonlinear law globally. The crucial point here is to get the knowledge a priori of where each region lies; we aim to answer this point in the present paper, while we leave the domain-decomposition step for future investigation. Methods a posteriori are also envisageable [12,13], though not explored here.

To locate the linear and nonlinear regions, we propose to use the adaptive model introduced in [14,15]. There, the authors derive a model where the constitutive law linking velocity and pressure gradient switches from linear to nonlinear according to a fixed threshold on the seepage-flux magnitude: when the magnitude is below the threshold, the linear law is solved, and when it is above, the nonlinear law is solved instead. As a result, a law which naturally adapts itself given the flux speed is obtained. Although well-posedness of the adaptive model is proved in [14,15], the adaptive law is discontinuous in speed, leading to a jump of pressure gradient at the transition zones between the two laws. Consequently, the adaptive model is rather unsuited for numerical simulation. To circumvent this, the authors in [15] tackle this discontinuity via a regularization argument. More precisely, an energetic, or variational, formulation of the discontinuous problem is derived and smoothed by convolving it with a Gaussian kernel. Then, as an optimization problem, the smooth constitutive model ensues and is shown to converge to the original adaptive model as the regularization parameter vanishes. This regularized law is precisely that which we propose to use in this paper to determine the partition of the medium into slow and fast regions. The result of this task is a labeling of each mesh cell as belonging to one or the other subdomain and strongly depends on the tolerance we give ourselves on the error allowed when solving Darcy's law rather than a nonlinear one; this tolerance is encoded in the very definition of the flux threshold.

This approach may seem to yield a slower model than a globally nonlinear one, since indeed the determination of the regions is carried out by solving another nonlinear, albeit smooth, problem, that is, the regularized problem. Overall then, solving the regularized model to obtain the regions and subsequently implementing a domain decomposition according to the obtained regions should be slower than simply solving the nonlinear model globally. However, we hope that our two-step approach (solving the regularized model followed by a domain decomposition) opens the possibility in a future work to introduce machine-learning techniques to determine the regions in a computationally efficient way, as started in [16].

The technique introduced here has potential applications in fractured porous media. Indeed, fractures can be modeled as lower-dimensional objects due to their thinness, yet they can act as highly permeable pathways for fluid flow. An example of this, where the model within the fractures is predefined and distinct from the surrounding medium, is discussed in [11]. Another possible application involves elastic models, where the appropriate model (e.g., elastic or hyperelastic) could be automatically selected based on the stress-displacement regime experienced by the material, as described in [17].

The paper is organized as follows: in Section 2, we summarize the various nonlinear laws considered by our model, we derive an expression for the flux threshold depending on the prescribed tolerance and we recall the adaptive and regularized models of [14,15]; in Section 3, we show the validity of the approach by experimenting on four different test cases in two and three space dimensions, the first two cases being inspired by an application in landfill management [18] and the remaining two cases being taken from the SPE10 benchmark reservoir scenario [19]. The codes used to obtain our results are written in Python and can be found at [20], while the underlying discretization and equation resolution are performed using the open source Python simulation tool PorePy [21,22].

## 2. Model

We place ourselves in the modeling context of [15]; for completeness, we summarize it here. In addition, we introduce some notions and quantities of particular interest in the present paper.

### 2.1. Notation

We denote by $\Omega \subset \mathbb{R}^d$ the porous medium, which we assume to be open, bounded and with Lipschitz boundary $\partial\Omega$; we write $n$ the outward normal unit vector of $\partial\Omega$. Also, we write $\mathbb{R}_+ = [0, \infty)$ and $\mathbb{N}_0 = \mathbb{N} \cup \{0\}$.

The unknowns of the problems discussed throughout the paper are the fluid's pressure $p: \Omega \to \mathbb{R}$ and the seepage flux $u: \Omega \to \mathbb{R}^d$ defined by $u = \rho\phi\nu$, where $\phi$ is the medium's porosity and $\rho$ and $\nu$ are the fluid's density and velocity. Furthermore, we write $\nu$ and $k$ the fluid's kinematic viscosity and the medium's scalar permeability. We suppose that $\rho, \nu, k: \Omega \to (0, \infty)$ and $\phi: \Omega \to (0, 1)$ are space-dependent knowns of the problem.

### 2.2. Conservation equations

We study the stationary flow of the fluid through the porous medium. The conservation of mass then reads

$$ \mathrm{div} \, u = q \quad \text{in } \Omega, $$

where $q: \Omega \to \mathbb{R}$ is a known fluid mass source, and the conservation of momentum is given by

$$ -\nabla p + f = \lambda[u] \quad \text{in } \Omega, \qquad (2.1) $$

where $f: \Omega \to \mathbb{R}^d$ is a known vector of external body forces, possibly including gravity. The conservation of momentum shows that the pressure gradient and the external forces balance the drag force $-\lambda[u]$ undergone by the fluid. We refer to (2.1) as the seepage, or constitutive, law and to $\lambda$ as the drag operator.
## 2.3. Drag operator

We let $\lambda$ be of the form

$$ \lambda[u] : \Omega \ni x \mapsto \alpha(x, \|u(x)\|)d(x)u(x), \qquad (2.2) $$

where $d: \Omega \to (0, \infty)$ is the medium's *drag*, that is, $d = v/k$, and $\alpha: \Omega \times \mathbb{R}_+ \to \mathbb{R}_+$ is a space- and flux-dependent function (the product $\alpha v$ possibly being interpreted as the kinematic viscosity perceived by the fluid in motion). This assumption physically means that $\lambda[u]$ is colinear with $u$ and has magnitude directly dependent on the position in $\Omega$ and on the Euclidean norm $\|u\|$. Moreover, it does not exhibit any anisotropy, be it from the permeability or other quantities. Extensions to tensor permeabilities are envisageable [15,23,24] and one possibility is in fact mentioned in the context of our three-dimensional test case (cf. Section 3.3), but further investigation is required. Regarding anisotropies coming from physical quantities other than permeability, there does not seem to be any consensus on how to incorporate them [23], so that we do not consider this possibility in this paper.

Volume-averaging and homogenization methods yield relevant information on the expression of the function $\alpha$ directly from the Navier–Stokes equations [4–8]. By considering the interaction between micro- and macroscopic scales, these methods can be used to motivate the following choice:

$$ \alpha = 1 + \eta, $$

where $\eta: \Omega \times \mathbb{R}_+ \to \mathbb{R}_+$ is referred to as the *inertial correction* and behaves linearly about 0 with respect to its second variable [8]; in particular, $\eta(x, 0) = 0$ for all $x \in \Omega$. With this information, we take $\lambda$ in (2.2) as

$$ \lambda[u] : \Omega \ni x \mapsto (1 + \eta(x, \|u(x)\|))d(x)u(x). \qquad (2.3) $$

Additional theoretical arguments require that $\eta$ be monotonously nondecreasing and coercive in its second variable [14,15], which intuitively implies that the viscosity perceived by the fluid cannot decrease with its speed.

## 2.4. Inertial correction

Unfortunately, volume-averaging and homogenization techniques do not give a unique and explicit expression of the inertial correction, so that assumptions on its form and experimental fitting of coefficients are needed.

### 2.4.1. Assumption

We assume that the inertial correction in (2.3) is such that $\eta(x, 0) = 0$ and $\eta(x, \cdot)$ is entire for all $x \in \Omega$, i.e., that it can be expanded as

$$ \eta(x, a) = \sum_{i=1}^{\infty} \eta_i(x) a^i \quad \text{for all } x \in \Omega \text{ and } a \ge 0, \qquad (2.4) $$

where, for all $i \in \mathbb{N}$, the empirical function $\eta_i: \Omega \to \mathbb{R}_+$ is defined by

$$ \eta_i(x) = \frac{1}{i!} \frac{\partial^i \eta}{\partial a^i}(x, 0). $$

A less general expression for the inertial correction is found in [23, Equation (5)].

Note that the entirety of $\eta$ sets the constraint that

$$ \lim_{i \to \infty} \sqrt[i]{\eta_i(x)} = 0 \quad \text{for all } x \in \Omega. $$

Further note that the nonnegativity of the functions $\eta_i$ is not a mathematically, nor physically, necessary condition but is assumed here for ease of presentation.

### 2.4.2. Coefficients

For all $i \in \mathbb{N}$, we introduce the following dimensionless function $c_i: \Omega \to \mathbb{R}_+$, referred to as *inertial coefficient*:

$$ c_i = \left( \frac{\mu}{\sqrt{k}} \right)^i \eta_i, $$

where $\mu := \rho v$ is the fluid's dynamic viscosity. With this notation, the inertial correction from (2.4) satisfies

$$ \eta(x, a) = \sum_{i=1}^{\infty} c_i(x) \left( \left( \sqrt{k}/\mu \right) (x) a \right)^i \quad \text{for all and } x \in \Omega \text{ and } a \ge 0. \qquad (2.5) $$

The nonzero inertial coefficients in (2.5) are precisely the quantities to be chosen experimentally to get an expression for the inertial correction.
### 2.4.3. Reynolds number

In the present porous-medium context [23], given a flux $u$, the Reynolds number $Re[u] : \Omega \to \mathbb{R}_+$ is set by

$$Re[u] = \frac{\sqrt{K}}{\mu} \|u\|.$$

With this notation, together with (2.5), one can rewrite the seepage law (2.1) as

$$-\nabla p + f = \left( 1 + \sum_{i=1}^{\infty} c_i \mathrm{Re}[u]^i \right) \mathbf{d}u \quad \text{in } \Omega. \qquad (2.6)$$

## 2.5. Truncated seepage laws

In practice, one cannot experimentally tune all the nonzero inertial coefficients since, a priori, these are infinitely many. Therefore, we propose to truncate the law in (2.6), i.e., to keep only a finite number of coefficients. Another way of reducing the number of coefficients is presented in Section 2.8.

We let $I \subset \mathbb{N}$ be the set of indices we want to keep in the truncated law, so that (2.6) rewrites as

$$-\nabla p + f = \left( 1 + \sum_{i=1}^{\infty} c_i \mathrm{Re}[u]^i \right) \mathbf{d}u \quad \text{in } \Omega. \qquad (2.7)$$

In the limit case $I = \mathbb{N}$, we recover the untruncated law (2.6). By considering finite index sets, a wide range of seepage laws can be obtained, many of which are already proposed in various theoretical and numerical studies [25–29]. There does not seem to be a general agreement on which is best.

### 2.5.1. Darcy’s law

Taking $I = \emptyset$ in (2.7) gives the widely used Darcy law:

$$-\nabla p + f = \mathbf{d}u \quad \text{in } \Omega.$$

Note that, here, no coefficient needs fitting so that this law is in fact not empirical.

### 2.5.2. Single-power law

For some $m \in \mathbb{N}$, we take $I = \{m\}$ in (2.7) to get

$$-\nabla p + f = (1 + c_F \mathrm{Re}[u]^m) \mathbf{d}u \quad \text{in } \Omega,$$

where the coefficient $c_m$ is preferably denoted by $c_F : \Omega \to \mathbb{R}$, which is the Forchheimer coefficient (or Ergün number [30,31]), containing information on porosity and tortuosity [6]. When $m = 1$, this is called the Darcy–Forchheimer law. Note that this formulation stays meaningful for any positive, noninteger value of $m$.

### 2.5.3. Power-expansion law

By setting $I = \{1, 2, \dots, m\}$ in (2.7) for some $m \in \mathbb{N}$, we yield

$$-\nabla p + f = (1 + c_1 \mathrm{Re}[u] + c_2 \mathrm{Re}[u]^2 + \dots + c_m \mathrm{Re}[u]^m) \mathbf{d}u \quad \text{in } \Omega,$$

which, as with the single-power law, includes the Darcy–Forchheimer law. Single-power and power-expansion laws may thus be called generalized Forchheimer laws.

## 2.6. Comparing seepage laws

Given a flux $u$ and two drag operators $\lambda$ and $\hat{\lambda}$, the resulting forces $-\nabla p + f$ and $-\nabla \hat{p} + f$ are obtained by the laws (cf. (2.1))

$$\begin{cases} -\nabla p + f = \lambda[u], \\ -\nabla \hat{p} + f = \hat{\lambda}[u], \end{cases} \quad \text{in } \Omega. \qquad (2.8)$$

The drag operator $\hat{\lambda}$ is referred to as the reference operator.

### 2.6.1. Untruncated error

Subtracting the two equations in (2.8) and using the untruncated law (2.6) with inertial coefficients $c := (c_i)_{i \in \mathbb{N}}$ and $\hat{c} := (\hat{c}_i)_{i \in \mathbb{N}}$, the local error $d_{c,\hat{c}}[u] : \Omega \to \mathbb{R}_+$ of the first law relative to the second is defined by

$$d_{c,\hat{c}}[u] = \frac{\left| \sum_{i=1}^{\infty} (c_i - \hat{c}_i) \mathrm{Re}[u]^i \right|}{1 + \sum_{i=1}^{\infty} \hat{c}_i \mathrm{Re}[u]^i}.$$
### 2.6.2. Truncated error

For truncated laws, subtracting again the two equations in (2.8) and using this time (2.7) with coefficients $c := (c_i)_{i \in \mathbb{N}}$ and $\hat{c} := (\hat{c}_i)_{i \in \mathbb{N}}$ and index sets $I$ and $\hat{I}$, we define the local error $d_{c,I;\hat{c},\hat{I}}[u] : \Omega \to \mathbb{R}_+$ as

$$d_{c,I;\hat{c},\hat{I}}[u] = \frac{|\sum_{i \in \hat{I}} c_i \operatorname{Re}[u]^i - \sum_{i \in \hat{I}} \hat{c}_i \operatorname{Re}[u]^i|}{1 + \sum_{i \in \hat{I}} \hat{c}_i \operatorname{Re}[u]^i}.$$

**Example 2.1 (Equal Coefficients and Darcy–Forchheimer).** Of particular interest to us is the case $c = \hat{c}$ and $I = \emptyset$. In this setting, the above error becomes

$$d_{\hat{c},\hat{I}}[u] := d_{\hat{c},\emptyset;\hat{c},\hat{I}}[u] = \frac{\sum_{i \in \hat{I}} \hat{c}_i \operatorname{Re}[u]^i}{1 + \sum_{i \in \hat{I}} \hat{c}_i \operatorname{Re}[u]^i} < 1.$$

More specifically, we shall consider $\hat{I} = \{m\}$, with $m \in \mathbb{N}$ and $\hat{c}_m = c_F$, that is, the case when the first law is Darcy's and the reference law is a single-power law; then,

$$d_{c_F,m}[u] := d_{\hat{c},\emptyset;\hat{c},\{m\}}[u] = \frac{c_F \operatorname{Re}[u]^m}{1 + c_F \operatorname{Re}[u]^m} =: \frac{\operatorname{Fo}_m[u]^m}{1 + \operatorname{Fo}_m[u]^m}, \quad (2.9)$$

where the $m$th Forchheimer number $\operatorname{Fo}_m[u] : \Omega \to \mathbb{R}_+$ is defined by

$$\operatorname{Fo}_m[u] = c_F^{1/m} \frac{\sqrt{\kappa}}{\mu} \|u\| = c_F^{1/m} \operatorname{Re}[u].$$

The dimensionless number $\operatorname{Fo}_m[u]^m$ quantifies the weight of the nonlinear term in the single-power law relative to 1, i.e., relative to the linear term. The first Forchheimer number, resulting from considering the Darcy–Forchheimer law as reference, is simply referred to as the Forchheimer number in the literature [32,33].

### 2.6.3. Flux threshold and subdomains

In the setting given by Example 2.1, there exists a function $\Delta_{\hat{c},\hat{I}} : \Omega \times \mathbb{R}_+ \to [0, 1)$ such that, for all $x \in \Omega$, there holds

$$d_{\hat{c},\hat{I}}[u](x) = \Delta_{\hat{c},\hat{I}}(x, \|u(x)\|),$$

and $\Delta_{\hat{c},\hat{I}}(x, \cdot)$ is an increasing bijection. Let $\delta \in [0, 1)$, referred to as error tolerance.

**Threshold.** By bijectivity, there is a unique $w_\delta : \Omega \to \mathbb{R}_+$ so that $\Delta_{\hat{c},\hat{I}}(x, w_\delta(x)) = \delta$ for all $x \in \Omega$. We then set the flux threshold $\bar{u}_\delta$ as

$$\bar{u}_\delta = \inf_{\Omega} w_\delta \in \mathbb{R}_+. \quad (2.10)$$

When the reference law is a single-power law with exponent $\hat{m} \in \mathbb{N}$, then, from (2.9), one finds

$$w_\delta = \operatorname{Fo}_\delta \frac{\mu}{c_F^{1/m} \sqrt{\kappa}},$$

where $\operatorname{Fo}_\delta := \sqrt[m]{\delta/(1-\delta)}$ is a critical Forchheimer number, which yields

$$\bar{u}_\delta = \operatorname{Fo}_\delta \inf_{\Omega} \frac{\mu}{c_F^{1/m} \sqrt{\kappa}}.$$

In particular, when $m=1$, i.e., the reference law is Darcy–Forchheimer, we get

$$w_\delta = \operatorname{Fo}_\delta \frac{\mu}{c_F \sqrt{\kappa}} \quad \text{and} \quad \bar{u}_\delta = \operatorname{Fo}_\delta \inf_{\Omega} \frac{\mu}{c_F \sqrt{\kappa}} = \frac{\delta}{1-\delta} \inf_{\Omega} \frac{\mu}{c_F \sqrt{\kappa}}.$$

**Subdomains.** We then define the subdomains

$$\Omega_\delta[u] = \{x \in \Omega \mid \|u(x)\| < \bar{u}_\delta\} \quad \text{and} \quad \hat{\Omega}_\delta[u] = \{x \in \Omega \mid \|u(x)\| > \bar{u}_\delta\},$$

which we respectively refer to as the slow and fast subdomains.

If $x \in \Omega_\delta[u]$, we know that solving at $x$ the law given by $\lambda$ (cf. first equation in (2.8)) does not generate a local error greater than $\delta$ relative to solving at $x$ the law stemming from $\hat{\lambda}$ (cf. second equation in (2.8)). Indeed, suppose $x \in \Omega_\delta[u]$; then,

$$\|u(x)\| < \bar{u}_\delta \le w_\delta(x),$$

so that

$$d_{\hat{c},\hat{I}}[u](x) = \Delta_{\hat{c},\hat{I}}(x, \|u(x)\|) < \Delta_{\hat{c},\hat{I}}(x, w_\delta(x)) = \delta.$$

Hence, by solving the seepage law for $\lambda$ in $\Omega_\delta[u]$ and that for $\hat{\lambda}$ in $\hat{\Omega}_\delta[u]$, we make sure that the local error stays below the error tolerance $\delta$. Also note that if $x \in \Omega$ and $\operatorname{Fo}_m[u](x) > w_\delta$, then $x \in \hat{\Omega}_\delta[u]$; this is why $\operatorname{Fo}_\delta$ is seen as a critical Forchheimer number above which we are certain that nonlinear effects are nonnegligible.
## 2.7. Final model

Suppose we are in the framework of Section 2.6 with $\delta \in [0, 1)$.

Furthermore, let $\Sigma_v, \Sigma_p \subset \partial\Omega$ be relatively open in $\partial\Omega$ (i.e., $\Sigma_v$ and $\Sigma_p$ are each the intersection of an open subset of $\mathbb{R}^d$ with $\partial\Omega$) and such that $\partial\Omega = \Sigma_v \cup \Sigma_p$ and $\Sigma_v \cap \Sigma_p = \emptyset$. We also let $u_0 : \Sigma_v \to \mathbb{R}$ and $p_0 : \Sigma_p \to \mathbb{R}$ be functions setting the conditions on the boundary for the flux and pressure. For simplicity, a map on $\Omega$ and its trace on $\partial\Omega$ are denoted by the same symbol.

### 2.7.1. Adaptive setting

We consider the problem where the seepage law for $\lambda$ applies in the slow subdomain and that for $\hat{\lambda}$ in the fast one. In practice, the reference drag operator $\hat{\lambda}$ is numerically more costly to handle, so that by only solving it in the fast subdomain, rather than the entire domain, we reduce the overall computational cost. This discontinuous problem, where the discontinuity in flux arises from the jump occurring between the seepage laws, is referred to as the *adaptive problem* and is as follows:

**Problem 2.2 (Adaptive).** Find $u : \Omega \to \mathbb{R}^d$ and $p : \Omega \to \mathbb{R}$ such that

$$ 
\left\{ 
\begin{array}{ll} 
\text{div } u = q & \text{in } \Omega, \\ 
-\nabla p + f = \lambda[u] & \text{in } \Omega_{\delta}[u], \\ 
-\nabla p + f = \hat{\lambda}[u] & \text{in } \hat{\Omega}_{\delta}[u], \\ 
u \cdot n = u_0 & \text{on } \Sigma_v, \\ 
p = p_0 & \text{on } \Sigma_p. 
\end{array} 
\right. 
\qquad (2.11) 
$$

The rigorous setting in which to solve such problem is described in [15], from where well-posedness can be deduced via a straightforward generalization of [15, Assumption 4.1 and Corollary 4.4] to space-dependent physical parameters, such as $\mu, \nu, \kappa$ and the inertial coefficients, which we further assume to be bounded below and above by strictly positive constants in $\Omega$. In [15], the appropriate functional framework is developed and a multivalued formulation is derived to tackle the flux discontinuity at the transition zone between the laws.

**Remark 2.3 (Average Pressure).** If the boundary piece $\Sigma_p$ verifies $\text{Vol}^{d-1}(\Sigma_p) = 0$, where $\text{Vol}^{d-1}$ is the $(d-1)$-dimensional Lebesgue measure, to ensure uniqueness of the pressure satisfying Problem 2.2, one imposes a constraint on the average of $p$:

$$ \frac{1}{|\Omega|} \int_{\Omega} p = \bar{p}, \qquad (2.12) $$

for a given $\bar{p} \in \mathbb{R}$. Tacitly, we thus require (2.12) in (2.11) whenever $\text{Vol}^{d-1}(\Sigma_p) = 0$.

### 2.7.2. Regularized setting

Numerically, it is not desirable to solve a discontinuous problem such as Problem 2.2. Indeed, since we use an explicit fixed-point algorithm, convergence is not guaranteed without continuity. We therefore follow the approach developed once again in [15], which consists in averaging the two drag operators $\lambda$ and $\hat{\lambda}$ via a convolution to obtain a smooth law. This introduces an additional parameter, $\epsilon > 0$, which vanishes as the regularization weakens. The details of this regularization are found in [15] and the corresponding Python code at [20]. Still, for completeness, we give here an informal derivation of this regularized law and the respective regularized problem.

The authors in [15] exploit the variational, or energetic, formulation of the adaptive model to regularize the law via a convolution of the underlying energy functional. It is shown there that a pair $(u, p)$ is solution to Problem 2.2 if and only if it is a saddle point of the energy functional $\mathcal{E}$ defined by

$$ \mathcal{E}(\varphi, \psi) = \int_{\Omega} f \cdot \varphi - \int_{\Omega} q\psi + \int_{\Sigma_v} u_0\psi - \int_{\Omega} \nabla\psi \cdot \varphi - D(\varphi), $$

where $D$ is called the *dissipation* and is minimized on a certain restricted function space by a certain component of the solution $u$ (cf. [14,15]). It encodes the information of the two chosen drag operators $\lambda$ and $\hat{\lambda}$ and is of the form

$$ D(\varphi) = \frac{1}{2} \int_{\Omega} \mathrm{d}1 (\Psi(\|\varphi\|^2) \mathbb{1}_{\Omega_{\delta}[\varphi]} + \hat{\Psi}(\|\varphi\|^2) \mathbb{1}_{\hat{\Omega}_{\delta}[\varphi]}), $$

where $\mathbb{1}_A$ stands for the indicator function of set $A$, and $\Psi, \hat{\Psi} : \mathbb{R}_+ \to \mathbb{R}_+$ satisfy

$$ \lambda[\varphi] = \Psi'(\|\varphi\|^2)\mathrm{d}1\varphi \quad \text{and} \quad \hat{\lambda}[\varphi] = \hat{\Psi}'(\|\varphi\|^2)\mathrm{d}1\varphi. $$

We note that the discontinuity of the adaptive law is reflected in the discontinuity of the integrand of $D$ with respect to the variable $\|\varphi\|$. In fact, $D$ is not differentiable and merely has a subdifferential, which is equal to the right-hand side of the adaptive law in Problem 2.2. To make the dissipation smooth, the integrand of $D$ is convolved with a Gaussian kernel $G_\epsilon$ of mean 0 and standard deviation some fixed $\epsilon > 0$. The convolution is applied on the space of flux magnitudes and yields the regularized dissipation $D_\epsilon$, defined as

$$ D_{\epsilon}(\varphi) = \frac{1}{2} \int_{\Omega} \mathrm{d}1 G_{\epsilon} * \left[ (\Psi(\|\varphi\|^2) \mathbb{1}_{\Omega_{\delta}[\varphi]} + \hat{\Psi}(\|\varphi\|^2) \mathbb{1}_{\hat{\Omega}_{\delta}[\varphi]}) \right], $$
and the corresponding energy $\mathcal{E}_\epsilon$, given by

$$ \mathcal{E}_\epsilon(\varphi, \psi) = \int_{\Omega} \mathbf{f} \cdot \varphi - \int_{\Omega} q\psi + \int_{\Sigma_v} u_0\psi - \int_{\Omega} \nabla\psi \cdot \varphi - D_\epsilon(\varphi). $$

The regularized problem then consists in finding a pair $(u_\epsilon, p_\epsilon)$ which is a saddle point of $\mathcal{E}_\epsilon$. Equivalently, after differentiation of the regularized dissipation, it consists in letting the seepage law be

$$ -\nabla p_\epsilon + \mathbf{f} = \lambda_{\delta,\epsilon}[u_\epsilon], $$

where the regularized drag operator $\lambda_{\delta,\epsilon}$ is given by

$$ \lambda_{\delta,\epsilon}[u_\epsilon] = \mathfrak{d} G'_\epsilon * \left[ \left( \Psi(\|u_\epsilon\|^2) \mathbb{1}_{\Omega_\delta[u_\epsilon]} + \hat{\Psi}(\|u_\epsilon\|^2) \mathbb{1}_{\hat{\Omega}_\delta[u_\epsilon]} \right) \right]. $$

This drag operator contains the information of the threshold flux $\bar{u}_\delta$ via the indicator functions and the resulting regularized problem is the following:

**Problem 2.4 (Regularized).** Find $u_\epsilon : \Omega \to \mathbb{R}^d$ and $p_\epsilon : \Omega \to \mathbb{R}$ such that

$$ \left\{ \begin{array}{ll} \operatorname{div} u_\epsilon = q & \text{in } \Omega, \\ -\nabla p_\epsilon + \mathbf{f} = \lambda_{\delta,\epsilon}[u_\epsilon] & \text{in } \Omega, \\ u_\epsilon \cdot \mathbf{n} = u_0 & \text{on } \Sigma_v, \\ p_\epsilon = p_0 & \text{on } \Sigma_p. \end{array} \right. $$

Problem 2.4 is indeed the problem we validate numerically in Section 3 as a tool to generate partitions of porous media into fast and slow subdomains. Thanks to [15, Corollary 5.3 and Theorem 5.4], it is well posed and converges in a weak sense to the adaptive problem as $\epsilon \to 0$. Thus, the slow and fast subdomains $\Omega_\delta[u_\epsilon]$ and $\hat{\Omega}_\delta[u_\epsilon]$ it generates differ from those stemming from the adaptive model, namely, $\Omega_\delta[u_\epsilon]$ and $\hat{\Omega}_\delta[u_\epsilon]$, but are expected, though not yet proven, to be not too different for a reasonably small choice of $\epsilon$, such as $\epsilon = 0.1$, which is the value we take for the numerical experiments below.

## 2.8. Remark: exponential seepage laws

As an intriguing observation, we show here an alternative way of reducing the number of inertial coefficients in the untruncated law (2.6), which yields laws different from those derived in Section 2.5 and, to our knowledge, not found in the current literature, although we do not explore them further in the present paper.

### 2.8.1. Derivation

The idea is to reduce the number of coefficients to a single one by requiring that there be a coefficient $c : \Omega \to \mathbb{R}$ satisfying

$$ i! c_i = c^i, \quad i \in \mathbb{N}. $$

Moreover, we give ourselves a truncating parameter $m \in \mathbb{N}_0 \cup \{\infty\}$.

From (2.5), for all $x \in \Omega$ and $a \ge 0$, we yield

$$ \begin{aligned} 1 + \eta(\mathbf{x}, a) &= \sum_{i=0}^{m} \frac{1}{i!} \left( \left( c\sqrt{\kappa}/\mu \right)(\mathbf{x})a \right)^i \\ &= G_m \left( \left( c\sqrt{\kappa}/\mu \right)(\mathbf{x})a \right) \exp \left( \left( c\sqrt{\kappa}/\mu \right)(\mathbf{x})a \right), \end{aligned} $$

where $G_m : \mathbb{R}_+ \to (0, 1]$ is defined, for all $b \in \mathbb{R}_+$, as

$$ G_m(b) = \begin{cases} \frac{\Gamma(m+1, b)}{\Gamma(m+1, 0)} & \text{if } m \in \mathbb{N}_0, \\ 1 & \text{if } m = \infty, \end{cases} $$

where $\Gamma$ is the upper incomplete gamma function [34]. Thus, (2.6) gives

$$ -\nabla p + \mathbf{f} = G_m(c \operatorname{Re}[\mathbf{u}]) \exp(c \operatorname{Re}[\mathbf{u}]) \mathfrak{d} \mathbf{u} \quad \text{in } \Omega. \quad (2.13) $$

We may call this law full exponential law if $m = \infty$ and truncated exponential law if $m \in \mathbb{N}_0$. This is again a generalized Forchheimer law since, from it, we retrieve the Darcy–Forchheimer law if $m = 1$ and $c = c_F$, and indeed Darcy’s law if $m = 0$. However, when $m > 1$, we yield laws different from those presented in Section 2.5.

Mathematically, this formulation makes sense for some nonpositive and noninteger values of $m$ since the gamma function is well defined on some of those values. Further note that $G_m \uparrow G_\infty$ pointwise as $m \to \infty$, so that indeed the truncated exponential law approaches the full law as the truncating parameter $m$ grows. Note, however, that well-posedness for such exponential laws is not guaranteed by the arguments in [15] since operator boundedness does not hold in this case.
### 2.8.2. Local error

Following the reasoning of Section 2.6.2 for exponential laws, using (2.13) with coefficients $c$ and $\hat{c}$ and truncating parameters $m$ and $\hat{m}$, the local error $d_{c,m;\hat{c},\hat{m}}[u] : \Omega \to \mathbb{R}_+$ may be defined as

$$d_{c,m;\hat{c},\hat{m}}[u] = \frac{|G_m(c \operatorname{Re}[u]) \exp(c \operatorname{Re}[u]) - G_{\hat{m}}(\hat{c} \operatorname{Re}[u]) \exp(\hat{c} \operatorname{Re}[u])|}{G_{\hat{m}}(\hat{c} \operatorname{Re}[u]) \exp(\hat{c} \operatorname{Re}[u])}.$$

**Example 2.5 (Equal Coefficients and Full Exponential).** The particular case $c = \hat{c}$ and $\hat{m} > m$ simplifies into

$$d_{m,c;\hat{m}}[u] := d_{c,m;c,\hat{m}}[u] = 1 - \frac{G_m(c \operatorname{Re}[u])}{G_{\hat{m}}(c \operatorname{Re}[u])} < 1,$$

which, when $\hat{m} = \infty$, further cuts down to

$$d_{m;c}[u] := d_{c,m;c,\infty}[u] = 1 - G_m(c \operatorname{Re}[u]). \qquad (2.14)$$

The considerations in Section 2.6.3 still hold in the exponential setting given in Example 2.5. Indeed, one can find $\Delta_{m;c,\hat{m}} : \Omega \times \mathbb{R}_+ \to [0, 1)$ so that, for all $x \in \Omega$, there holds

$$d_{m,c;\hat{m}}[u](x) = \Delta_{m;c,\hat{m}}(x, \|u(x)\|),$$

where $\Delta_{m;c,\hat{m}}(x, \cdot)$ is an increasing bijection. Using the same definition as in (2.10), when the reference law is the full exponential, from (2.14), one gets

$$w_{\delta} = G_{m}^{-1}(1 - \delta) \frac{\mu}{c \sqrt{\kappa}},$$

which yields

$$\bar{u}_{\delta} = G_{m}^{-1}(1 - \delta) \inf_{\Omega} \frac{\mu}{c \sqrt{\kappa}}.$$

The inverse $G_m^{-1}$ can, in general, be computed numerically; but, if, for instance $m = 0$, that is, the drag operator $\lambda$ is Darcy's, then $G_0^{-1} = -\ln$, so that

$$\bar{u}_{\delta} = -\ln(1 - \delta) \inf_{\Omega} \frac{\mu}{c \sqrt{\kappa}}.$$

## 3. Numerical validation

In this section, we present four test cases to validate the performance and quality of our approach, that is, the results of the numerical resolution of Problem 2.4 using Raviart-Thomas finite elements and an explicit fixed-point algorithm, as described in detail in [15]. Note that, in this section, as opposed to the rest of the paper, we refer to the model in Problem 2.4 as the *adaptive* model, rather than *regularized*; this is to keep stressing the fact that there always is an underlying adaptability between the nonlinear and linear laws in the numerically solved model.

In Section 3.1, two examples are proposed with two or more highly conductive channels and a background heterogeneity, as inspired by [18]. The examples in Sections 3.2 and 3.3 are based on the SPE10 benchmark study [19], where we consider, respectively, a single layer and multiple layers.

### 3.1. Highly permeable channels

This test case is based on the examples reported in [18] for a Darcy-Forchheimer system, and adapted here for our purposes. We consider in particular two configurations depending on the number of highly conductive channels that are present in the problem. In the first case, reported in Section 3.1.1, two vertical channels are present, while in Section 3.1.2 a network of vertical and horizontal channels is considered.

For both tests, we consider the following physical data for the fluid: dynamic viscosity $\mu = 0.001[\text{Pa s}]$ and density $\rho = 998.0[\text{kg m}^{-3}]$. For the porous medium, the permeability is reported in Fig. 1 for each configuration. To pass from permeability to porosity and vice-versa when needed, we use the Kozeny-Carman law:

$$\kappa = \kappa_{\text{ref}} \frac{(1 - \phi_{\text{ref}})^2}{\phi_{\text{ref}}^3} \frac{\phi^3}{(1 - \phi)^2},$$

where $\kappa_{\text{ref}} = 1.0152441851 \cdot 10^{-9}[\text{m}^2]$ and $\phi_{\text{ref}} = 35[\%]$ (cf. [18]).

Unless stated otherwise, we fix the error tolerance $\delta = 0.1[-]$, yielding the critical Forchheimer number $F_{0,\delta} \approx 0.11[-]$ (cf. Section 2.6); moreover, as chosen in [18], we let $m = 1$, corresponding to the classical Darcy-Forchheimer model, with coefficient $c_F = 0.55[-]$ (cf. Section 2.5).

For each case, we consider two scenarios for different values of the top inflow magnitude: Scenario a, with $u_{\text{in}} = 7 \cdot 10^{-3}[\text{kg m}^{-2} \text{s}^{-1}]$; Scenario b, with $u_{\text{in}} = 7 \cdot 10^{-2}[\text{kg m}^{-2} \text{s}^{-1}]$. On the left and right parts of the boundary, we set zero flux, while on the bottom part, we impose the hydrostatic pressure $p_{\text{atm}} + \rho g h$, where $p_{\text{atm}} = 1.01325 \cdot 10^5[\text{Pa}]$, $h = 10[\text{m}]$ and $g = 9.81[\text{m s}^{-2}]$.
<figure><img src="image_5.png" /><figcaption>Fig. 1. On the left, the log-permeability for the example of Section 3.1.1; on the right, for Section 3.1.2.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 2. Solutions obtained to the adaptive model for the example in Section 3.1.1, in the top row for Scenario a and in the bottom row for Scenario b. The pressure is multiplied by 1000, and the velocity arrows inside the channels are scaled by 1/6 and 1/21 in comparison with those outside, respectively for each scenario.</figcaption></figure>

### 3.1.1. Two vertical channels

We consider here the case where only two highly conductive channels are present in the domain, spanning the domain from top to bottom. Their permeability is roughly two orders of magnitude higher than that of the surrounding porous medium, and they are aligned with the flow direction given by the top boundary condition (cf. left in Fig. 1). Depending on the flux intensity, Scenario a or Scenario b, we expect to observe a Darcy–Forchheimer, i.e., nonlinear, flow developing inside these channels and a Darcy, i.e., linear, flow outside.

In Fig. 2, we report the solution obtained for both scenarios. We notice that for Scenario b, as expected, the flow in the two channels is higher than in Scenario a, resulting in a higher Forchheimer number so that more cells are selected in the Darcy–Forchheimer region, or subdomain (identified in the legend of the figure with the value 1). For Scenario b, both channels along with some extra cells are selected to be in the Darcy–Forchheimer region. For Scenario a, only cells in the channels are selected; more precisely, all cells are selected in the left channel, possibly because it is narrower than the right one and thus develops a higher flow rate.

Table 1 shows the $L^2$ errors computed relative to a reference solution, namely, the solution to the global Darcy–Forchheimer model, that is, the Darcy–Forchheimer model applied in the whole domain, without any adaptivity involved. The errors are rather small, although we see that if a global Darcy model is used, we make an error that reaches several orders of magnitude more than with the adaptive model, especially for Scenario b, even if the cells selected in the Darcy–Forchheimer subdomain in the adaptive model are many fewer than those in the Darcy subdomain. In both scenarios, the flux $u$ is affected by higher errors than the pressure $p$.
<caption>Table 1</caption>

L² errors of the adaptive and global Darcy solutions relative to the global Darcy–Forchheimer solutions, broken down in the Darcy–Forchheimer, Darcy and whole domains, for both scenarios of the example in Section 3.1.1.

<table><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Forchheimer region</th><th colspan="2">Darcy region</th><th colspan="2">Whole domain</th></tr><tr><th>Darcy</th><th>Adaptive</th><th>Darcy</th><th>Adaptive</th><th>Darcy</th><th>Adaptive</th></tr></thead><tbody><tr><td rowspan="2">Sce. a</td><td><i>err</i><sub>p</sub></td><td>5.01e-10</td><td>3.64e-12</td><td>4.31e-10</td><td>8.36e-11</td><td>4.33e-10</td><td>8.28e-11</td></tr><tr><td><i>err</i><sub>u</sub></td><td>6.89e-06</td><td>1.47e-06</td><td>4.95e-06</td><td>2.61e-06</td><td>6.28e-06</td><td>1.95e-06</td></tr><tr><td rowspan="2">Sce. b</td><td><i>err</i><sub>p</sub></td><td>3.98e-08</td><td>1.70e-12</td><td>4.35e-08</td><td>9.68e-12</td><td>4.33e-08</td><td>9.37e-12</td></tr><tr><td><i>err</i><sub>u</sub></td><td>6.20e-05</td><td>3.43e-08</td><td>1.27e-04</td><td>3.68e-07</td><td>6.28e-05</td><td>4.83e-08</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 3. Darcy (red) and Darcy–Forchheimer (blue) regions depending on the error tolerance $\delta$ decreasing from left to right according to $\delta = 0.05$, $\delta = 0.0125$ and $\delta = 0.003125$, in the top row for Scenario a and in the bottom row for Scenario b, for the example in Section 3.1.1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

In Fig. 3, we decrease the error tolerance $\delta$ for both scenarios. We notice that by setting lower and lower values of $\delta$, coherently, the Darcy–Forchheimer cells become more and more numerous. First, the two channels become entirely of Darcy–Forchheimer type; then, the top boundary condition and the heterogeneities of the background permeability start to play a role in determining where the rest of the Darcy–Forchheimer subdomain develops and thus which cells should be selected. By taking $\delta$ even smaller, we expect that all the cells will be selected to be of Darcy–Forchheimer type (see Fig. 3).

### 3.1.2. Network of channels

In this example, we consider the same background permeability as in the previous case, on which we place a network of intersecting horizontal and vertical channels with same width. All these channels are completely immersed in the domain without touching the boundary and they have permeabilities much higher than the surrounding medium. We expect, compared with the previous case, a more intricate flow distribution and thus a complex subdivision of cells into linear and nonlinear regions.

In Fig. 4, we report the solutions for Scenario a and Scenario b by, recall, fixing the value of the error tolerance $\delta = 0.1$. As expected, $||u||$ for Scenario b is much higher than for Scenario a, which results in a higher Forchheimer number and consequently in more cells marked as Darcy–Forchheimer. We notice that, for Scenario a, the Darcy–Forchheimer cells lie mostly inside the network, with only a few more cells leaking out of it and connecting with the bottom boundary. In contrast, for Scenario b, two rather large regions, about the top and bottom of the domain, are marked as Darcy–Forchheimer and connect with the network.

Table 2 shows the relative errors made by the global Darcy model and the adaptive model with respect to the global Darcy–Forchheimer model, considered thus to be the reference. As for the two-channel case, we notice that the adaptive scheme is more accurate, by several order of magnitude, than the global Darcy model. In particular, the detrimental impact of the error made in the cells marked as Darcy–Forchheimer brings down the error in the whole domain, even if these cells are many fewer than the Darcy cells. The adaptive approach gives thus rather accurate results.

In Fig. 5, we show the region split by decreasing the error tolerance $\delta$. As expected, the cells marked as Darcy–Forchheimer grow in number by filling up the zones near the top and bottom parts of the domain and then the areas in between the highly conductive channels.
<figure><img src="image_8.png" /><figcaption>Fig. 4. Solutions obtained to the adaptive model for the example in Section 3.1.2, in the top row for Scenario a and in the bottom row for Scenario b. The pressure is multiplied by 1000, and the velocity arrows inside the channels are scaled by 1/11 and 1/120 in comparison with those outside, respectively for each scenario.</figcaption></figure>

<caption>Table 2</caption>

$L^2$ errors of the adaptive and global Darcy solutions relative to the global Darcy–Forchheimer solutions, broken down in the Darcy–Forchheimer, Darcy and whole domains, for both scenarios of the example in Section 3.1.2.

<table><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Forchheimer region</th><th colspan="2">Darcy region</th><th colspan="2">Whole domain</th></tr><tr><th>Darcy</th><th>Adaptive</th><th>Darcy</th><th>Adaptive</th><th>Darcy</th><th>Adaptive</th></tr></thead><tbody><tr><td rowspan="2">Sce. a</td><td>err<sub>p</sub></td><td>6.08e-10</td><td>5.56e-13</td><td>6.03e-10</td><td>6.67e-13</td><td>6.03e-10</td><td>6.60e-13</td></tr><tr><td>err<sub>u</sub></td><td>1.43e-04</td><td>8.02e-08</td><td>4.76e-04</td><td>5.27e-07</td><td>1.56e-04</td><td>1.07e-07</td></tr><tr><td rowspan="2">Sce. b</td><td>err<sub>p</sub></td><td>6.53e-08</td><td>1.85e-13</td><td>5.75e-08</td><td>1.24e-13</td><td>6.02e-08</td><td>1.47e-13</td></tr><tr><td>err<sub>u</sub></td><td>1.55e-03</td><td>7.54e-10</td><td>9.02e-04</td><td>3.59e-08</td><td>1.55e-03</td><td>9.09e-10</td></tr></tbody></table>

### 3.1.3. Behavior with error tolerance

For both the two-channel and network examples, we plot in Fig. 6 the errors obtained on the adaptive model relative to the global Darcy–Forchheimer model as a function of the error tolerance $\delta$. The dashed lines for the relative errors, corresponding to Scenario b, drop to zero starting from a certain value of $\delta$, since, from that value on, all cells in the adaptive approach are classified as Darcy–Forchheimer. This phenomenon is not seen in Scenario a since the intensity of the top condition is one order of magnitude smaller than in the other scenario. In all examples and scenarios, the trend of the error is as expected: by decreasing $\delta$, the relative errors decrease and, simultaneously, the numbers of cells belonging to the Darcy–Forchheimer subdomain increase.

## 3.2. One layer from SPE10

In this example, we consider a single layer, namely, Layer 35, of the SPE10 benchmark study [19] represented as a bidimensional domain. For the fluid, we consider the following physical parameters: the dynamic viscosity $\mu = 3 \cdot 10^{-4}[\text{Pa} \cdot \text{s}]$ and the density $\rho = 1025[\text{kg m}^{-3}]$. The permeability is given by the benchmark study and it is represented in Fig. 7, where we notice some highly permeable paths surrounded by areas of low permeability.

As in the previous examples, the error tolerance is by default set to $\delta = 0.1$; the nonlinear exponent $m = 1$, so that we are considering the Darcy–Forchheimer law, and the Forchheimer coefficient is set as $c_F = 0.55[-]$.

We consider zero flux all around the boundary; moreover, we set five wells: one injector in the center of the domain and four identical producers at the corners. We consider different scenarios depending on the flux imposed in the injection well: in Scenario a, we set $u_{in} = 10[\text{kg s}^{-1}]$; in Scenario b, $u_{in} = 50[\text{kg s}^{-1}]$; in Scenario c, $u_{in} = 200[\text{kg s}^{-1}]$.

We represent in Fig. 8 the solution to the adaptive model obtained for the different scenarios. We notice that by increasing the inflow at the injection well, the cells marked as Darcy–Forchheimer become more numerous and, in particular, are preferably selected along the high-permeability channels connecting some of the wells. Even if the two wells on the bottom left and top right are not in highly permeable areas, for Scenario c, some Darcy–Forchheimer cells start to be selected around these two wells and, in particular for the second one, in the highly permeable region separating the well itself from the central injector.
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Darcy (red) and Darcy–Forchheimer (blue) regions depending on the error tolerance $\delta$ decreasing from left to right according to $\delta = 0.05$, $\delta = 0.0125$ and $\delta = 0.003125$, in the top row for Scenario a and in the bottom row for Scenario b, for the example in Section 3.1.2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Pressure (in red) and flux (in blue) $L^2$ relative errors by taking smaller values of the error tolerance $\delta$. The solid lines refer to Scenario a and the dashed lines to Scenario b for both the examples of Sections 3.1.1 (two-channel) and 3.1.2 (network). In green, the ratio between the number of Darcy–Forchheimer cells and the total number of cells. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Fig. 9 shows the Darcy–Forchheimer cells with decreasing error tolerance $\delta$. We notice that the number of cells belonging to the Darcy–Forchheimer region increases by first following highly permeable paths from the injection well to the production wells, and then accepting also regions of lower permeability. In the last figure, the error tolerance is so small that most of the cells are categorized as Darcy–Forchheimer. We expect that by reducing $\delta$ even further, we would get all the cells into the nonlinear subdomain.

In Fig. 10, we show the $L^2$ errors in pressure and flux for the adaptive model relative to the global Darcy–Forchheimer model, as well as the number of Darcy–Forchheimer cells selected by the adaptive model, all as functions of the error tolerance $\delta$. As for the examples of Section 3.1, by lowering the error tolerance, the adaptive approach tends to have a solution that is more similar to the global Darcy–Forchheimer model, and thus yields a smaller error. In fact, when the percentage of the Darcy–Forchheimer cells reaches 1, the adaptive model is equivalent to the global Darcy–Forchheimer one and thus has an error of exactly zero.
<figure><img src="image_11.png" /><figcaption>Fig. 7. Logarithm of permeability for the example in Section 3.2.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 8. Solutions obtained to the adaptive model for the example in Section 3.2, in the top row for Scenario a, in the middle row for Scenario b and in the bottom row for Scenario c. The velocity arrows are scaled by 1/92, 1/23 and 1/4, respectively for each scenario.</figcaption></figure>

### 3.3. Three-dimensional test case

In this last test case, we consider a tridimensional portion of the SPE10 benchmark [19], composed of the first 15 layers. However, to keep the computational cost limited, we select only the first 10 cells in the x direction and the first 20 in the y direction. As done in the previous example, we consider five wells: one in the center of the domain, the injector, and four in the corners, the producers. Here, the wells are represented as monodimensional objects and a Peaceman model [35,36] is used for the coupling with the surrounding porous medium. The wells are placed in the middle of each of their cells and are nonmatching with the
<figure><img src="image_13.png" /><figcaption>Fig. 9. Darcy (red) and Darcy–Forchheimer (blue) regions depending on the error tolerance $\delta$ decreasing from left to right according to $\delta = 0.05$, $\delta = 0.0125$ and $\delta = 0.003125$, in the top row for Scenario a and in the bottom row for Scenario c, for the example in Section 3.2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 10. Pressure (in red) and flux (in blue) $L^2$ relative errors by taking smaller values of the error tolerance $\delta$. The solid lines refer to Scenario a, the dashed lines to Scenario b, and the dashed-and-dotted lines to Scenario c, for the example of Section 3.2. In green, the ratio between the number of Darcy–Forchheimer cells and the total number of cells. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 11. Log-permeability field for the example in Section 3.3 in the horizontal (left) and vertical (right) directions.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. Darcy (opaqued red) and Darcy–Forchheimer (blue) regions for different well depths for the example in Section 3.3. From left to right: Scenario a (3.5 cells deep), Scenario b (7.5) and Scenario c (9.5). The velocity arrows for Scenario a are halved compared to the other two scenarios. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

porous-medium grid. By keeping the injection flow rate fixed, equal to 10[kg s⁻¹], we change the depth of the wells and see the impact on the selection procedure of Darcy and Darcy–Forchheimer cells by the adaptive model.

Although the models presented in Section 2 only consider a scalar permeability K, we wish here to take into account the anisotropy given by the SPE10 case study in the vertical direction. To do so, we alter the Darcy–Forchheimer law slightly:

$$-\nabla p + f = (1 + c_F \operatorname{Re}[u]) \nu K^{-1} u \quad \text{in } \Omega,$$

where K is now the permeability tensor and the Reynolds number is updated into

$$\operatorname{Re}[u] = \frac{\det(K)^{1/6}}{\mu} \|u\|.$$

We refer the reader to [23,24] for a justification of such a law. Note, however, that this form does not fit in the context of [14,15], where the adaptive model is derived and shown to be well posed, since, for that, the norm of the flux in the Reynolds number would need to be weighed by K⁻¹. We leave the precise numerical treatment of the tensor model for future investigation and accept for now the validity, in our context, of the above three-dimensional version of the Darcy–Forchheimer law. We have in fact, for the SPE10 study, a permeability tensor which is diagonal with equal values in the x and y directions and a different value in the z direction. A representation of the permeability field is reported in Fig. 11.

In Fig. 12, we show a comparison of the obtained adaptive solution by making the wells deeper: 3.5 cells for Scenario a, 7.5 cells for Scenario b, and 9.5 cells for Scenario c. We notice that, by considering the wells deeper, we get more Darcy–Forchheimer cells since indeed more regions of high permeability are thus connected. For the same reason, we observe that the number of Darcy–Forchheimer cells on the top layer reduces in favor of deeper connections.

## 4. Conclusion

In conclusion, our study presents a strategy for an automatic selection of linear- and nonlinear-flow regions in highly heterogeneous porous media, particularly answering the challenges raised by deviations from Darcy’s law in highly permeable areas. Through a combination of theoretical derivations, empirical observations, and numerical simulations, we have demonstrated the significance of nonlinear corrections to Darcy’s law by accurately capturing complex flow behaviors, especially in regions with high permeability.

Moreover, our approach, by using an adaptive regularized model based on a physically motivated flux threshold, offers a promising strategy for identifying and distinguishing between linear and nonlinear regimes within porous media. A local tolerance on the error, from which the threshold value is derived, has been introduced to switch between the linear and nonlinear laws effectively based on local physical properties. The experimental validation we have conducted across several numerical test cases shows the effectiveness of our methodology, although, as mentioned in the introduction, more investigation is needed to combine our strategy with a domain-decomposition technique and thus yield an complete, accurate and computationally efficient approach for flows in highly heterogeneous porous media.
Further research directions, each presenting its own challenges, include extending the proposed scheme to tensor-valued permeability (a tentative approach being used in Section 3.3 herein), incorporating space derivatives of the velocity into the drag operator (thus including Brinkman's model), and addressing time-dependent compressible flow. Numerically, it would also be valuable to explore more complex and realistic scenarios to assess the practical effectiveness of our methodology.

## CRediT authorship contribution statement

**Alessio Fumagalli:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.
**Francesco Saverio Patacchini:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.

## Acknowledgments

The present research is part of the activities of “Dipartimento di Eccellenza 2023–2027”, Italian Minister of University and Research (MUR), grant Dipartimento di Eccellenza 2023–2027.

## References

[1] W. Sobieski, A. Trykozko, Darcy's and Forchheimer's laws in practice. Part 1. The experiment, Tech. Sci. 17 (4) (2014) 321–335.

[2] W. Sobieski, A. Trykozko, Darcy's and Forchheimer's laws in practice. Part 2. The numerical model, Tech. Sci. 17 (4) (2014) 337–350.

[3] T. Zhang, Y. Zhao, Q. Gan, L. Yuan, G. Zhu, Y. Cai, B. Cao, Experimental investigation of forchheimer coefficients for non-Darcy flow in conglomerate-confined aquifer, Geofluids (2018) 1–21, http://dx.doi.org/10.3390/w14040546.

[4] E. Marušić-Paloka, A. Mikelić, The derivation of a nonlinear filtration law including the inertia effects via homogenization, Nonlinear Anal. Theory Methods Appl. 42 (1) (2000) 97–137, http://dx.doi.org/10.1016/S0362-546X(98)00346-0.

[5] S.P. Neuman, Theoretical derivation of Darcy's law, Acta Mech. 25 (3–4) (1977) 153–170, http://dx.doi.org/10.1007/BF01376989.

[6] D. Ruth, H. Ma, On the derivation of the forchheimer equation by means of the averaging theorem, Transp. Porous Media 7 (3) (1992) 255–264, http://dx.doi.org/10.1007/BF01063962.

[7] S. Whitaker, Flow in porous media I: A theoretical derivation of Darcy's law, Transp. Porous Media 1 (1) (1986) 3–25, http://dx.doi.org/10.1007/BF01036523.

[8] S. Whitaker, The Forchheimer equation: A theoretical development, Transp. Porous Media 25 (1) (1996) 27–61, http://dx.doi.org/10.1007/BF00141261.

[9] M.S. Espedal, K.J. Hersvik, B.G. Ersland, Domain decomposition methods for flow in heterogeneous porous media, Contemp. Math. 218 (1998) 104–120.

[10] J.O. Skogestad, E. Keilegavlen, J.M. Nordbotten, Domain decomposition strategies for nonlinear flow problems in porous media, J. Comput. Phys. 234 (2013) 439–451, http://dx.doi.org/10.1016/j.jcp.2012.10.001.

[11] E. Ahmed, A. Fumagalli, A. Budiša, A multiscale flux basis for mortar mixed discretizations of reduced Darcy–Forchheimer fracture models, Comput. Methods Appl. Mech. Engrg. 354 (2019) 16–36, http://dx.doi.org/10.1016/j.cma.2019.05.034, URL https://www.sciencedirect.com/science/article/pii/S0045782519303044.

[12] F. Févotte, A. Rappaport, M. Vohralík, Adaptive regularization, discretization, and linearization for nonsmooth problems based on primal–dual gap estimators, Comput. Methods Appl. Mech. Engrg. 418 (2024) 116558, http://dx.doi.org/10.1016/j.cma.2023.116558.

[13] D.A. Di Pietro, M. Vohralík, S. Yousef, Adaptive regularization, linearization, and discretization and a posteriori error control for the two-phase Stefan problem, Math. Comp. 84 (2015) 153–186.

[14] A. Fumagalli, F.S. Patacchini, Model adaptation for non-linear elliptic equations in mixed form: existence of solutions and numerical strategies, ESAIM Math. Model. Numer. Anal. 56 (2) (2022) 565–592, http://dx.doi.org/10.1051/m2an/2022016, URL https://www.esaim-m2an.org/articles/m2an/abs/2022/02/m2an210134/m2an210134.html.

[15] A. Fumagalli, F.S. Patacchini, Well-posedness and variational numerical scheme for an adaptive model in highly heterogeneous porous media, J. Comput. Phys. 475 (2023) 111844, http://dx.doi.org/10.1016/j.jcp.2022.111844, URL https://www.sciencedirect.com/science/article/pii/S002199912200907X.

[16] C. Giovannini, A. Fumagalli, F. Patacchini, Predicting nonlinear-flow regions in highly heterogeneous porous media using adaptive constitutive laws and neural networks, 2024, ArXiv preprint 2408.00533.

[17] O. Coussy, Poromechanics, Wiley, 2004.

[18] R. Winter, A. Valsamidou, H. Class, B. Flemisch, A Study on Darcy versus Forchheimer Models for Flow through Heterogeneous Landfills including Macropores, Water 14 (4) (2022) http://dx.doi.org/10.3390/w14040546, URL https://www.mdpi.com/2073-4441/14/4/546.

[19] M.A. Christie, M.J. Blunt, SPE-66599-MS, Soc. Pet. Eng., Houston, Texas, 2001, p. 13, http://dx.doi.org/10.2118/66599-MS.

[20] A. Fumagalli, F.S. Patacchini, https://github.com/compgeo-mox/validate_adaptative,

[21] E. Keilegavlen, R. Berge, A. Fumagalli, M. Starnoni, I. Stefansson, J. Varela, I. Berre, PorePy: An open-source software for simulation of multiphysics processes in fractured porous media, Comput. Geosci. (2020) http://dx.doi.org/10.1007/s10596-020-10002-5, URL https://link.springer.com/article/10.1007%2Fs10596-020-10002-5.

[22] https://github.com/pmgbergen/porepy,

[23] F. Russo Spena, A. Vacca, A potential formulation of non-linear models of flow through anisotropic porous media, Transp. Porous Media 45 (2001) 405–421.

[24] P.M. Knupp, J.L. Lage, Generalization of the Forchheimer-extended Darcy flow model to the tensor premeability case via a variational principle, J. Fluid Mech. 299 (1995) 97–104, http://dx.doi.org/10.1017/S0022112095003430.

[25] E. Aulisa, L. Bloshanskaya, L. Hoang, A. Ibragimov, Analysis of generalized Forchheimer flows of compressible fluids in porous media, J. Math. Phys. 50 (10) (2009) 103102, http://dx.doi.org/10.1063/1.3204977.

[26] L. Hoang, A. Ibragimov, Structural stability of generalized Forchheimer equations for compressible fluids in porous media, Nonlinearity 24 (1) (2011) 1–41, http://dx.doi.org/10.1088/0951-7715/24/1/001.

[27] V.A. Kovtunenko, Mixed variational problem for a generalized Darcy–Forchheimer model driven by hydraulic fracture, J. Vib. Test. Dyn. Syst. 7 (1) (2023) 15–21, http://dx.doi.org/10.5890/JVTSD.2023.03.003.

[28] J.D. Audu, F.A. Fairag, S.A. Messaoudi, On the well-posedness of generalized Darcy–Forchheimer equation, Bound. Value Probl. 2018 (123) (2018) http://dx.doi.org/10.1186/s13661-018-1043-6.
[29] M. Firdaouss, J.-L. Guermond, P. Le Quéré, Nonlinear corrections to Darcy's law at low Reynolds numbers, J. Fluid Mech. 343 (1997) 331–350, http://dx.doi.org/10.1017/S0022112097005843.

[30] A.P. Baitharu, S. Sahoo, G.C. Dash, Numerical approach to non-Darcy mixed convective flow of non-Newtonian fluid on a vertical surface with varying surface temperature and heat source, Karbala Int. J. Mod. Sci. 6 (3) (2020).

[31] A. Barhoi, C. Gohain, G.C. Hazarika, M. Clement Joe Anand, Numerical solution of skin friction in porous medium using Brinkman Forchheimer model under fuzzy environment, Int. J. Appl. Comput. preprint (2022).

[32] P. Macini, E. Mesini, R. Viola, Laboratory measurements of non-Darcy flow coefficients in natural and artificial unconsolidated porous media, J. Pet. Sci. Eng. 77 (3) (2011) 365–374, http://dx.doi.org/10.1016/j.petrol.2011.04.016.

[33] Z. Zeng, R. Grigg, A criterion for non-Darcy flow in porous media, Transp. Porous Media 63 (2006) 57–69.

[34] M.A. Chaudhry, S.M. Zubair, On a Class of Incomplete Gamma Functions with Applications, Chapman & Hall, 2001.

[35] D.W. Peaceman, Interpretation of well-block pressures in numerical reservoir simulation (includes associated paper 6988), Soc. Pet. Eng. J. 18 (03) (1978) 183–194, http://dx.doi.org/10.2118/6893-PA.

[36] D.W. Peaceman, Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability, Soc. Pet. Eng. J. 23 (03) (1983) 531–543, http://dx.doi.org/10.2118/10528-PA.