# Model Adaptation for Non Linear Elliptic Equations in Mixed Form Existence of Solutions and Numerical Strategies

Alessio Fumagalli, Francesco Saverio Patacchini

► To cite this version:

Alessio Fumagalli, Francesco Saverio Patacchini. Model Adaptation for Non Linear Elliptic Equations in Mixed Form Existence of Solutions and Numerical Strategies. ESAIM: Mathematical Modelling and Numerical Analysis, 2022, 56 (2), pp.565 - 592. 10.1051/m2an/2022016. hal-03625900

HAL Id: hal-03625900

https://ifp.hal.science/hal-03625900

Submitted on 31 Mar 2022

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
# MODEL ADAPTATION FOR NON-LINEAR ELLIPTIC EQUATIONS IN MIXED FORM: EXISTENCE OF SOLUTIONS AND NUMERICAL STRATEGIES

ALESSIO FUMAGALLI<fnref n="1" />,*, AND FRANCESCO SAVERIO PATACCHINI<fnref n="2" />

**Abstract.** Depending on the physical and geometrical properties of a given porous medium, fluid flow can behave differently, going from a slow Darcian regime to more complicated Brinkman or even Forchheimer regimes for high velocity. The main problem is to determine where in the medium one regime is more adequate than others. In order to determine the low-speed and high-speed regions, this work proposes an adaptive strategy which is based on selecting the appropriate constitutive law linking velocity and pressure according to a threshold criterion on the magnitude of the fluid velocity itself. Both theoretical and numerical aspects are considered and investigated, showing the potentiality of the proposed approach. From the analytical viewpoint, we show existence of weak solutions to such model under reasonable hypotheses on the constitutive laws. To this end, we use a variational approach identifying solutions with minimizers of an underlying energy functional. From the numerical viewpoint, we propose a one-dimensional algorithm which tracks the transition zone between the low- and high-speed regions. By running numerical experiments using this algorithm, we illustrate some interesting behaviors of our adaptive model on academic cases and on small networks of intersecting fractures.

**Mathematics Subject Classification.** 76S05, 35A01, 35A15, 65N50, 65N55.

Received July 16, 2021. Accepted February 1, 2022.

## 1. INTRODUCTION

Porous media are found in many applications related to the exploitation of the underground for energy management and extraction. The void space between the rock grains is filled by one or more fluids which move at relatively high of low speed, depending on many factors. Because of subsequent mineralization by chemical reactions, the void space of a porous medium may be partially or completely filled by additional material, substantially altering its physical properties and changing the fluid circulation. Additional rock deformation may change even more the hydraulic properties of the medium by, e.g., compacting the grain and leaving less free space for the fluid. A noticeable example of this is a fractured porous medium since fractures may have substantially different properties from the surrounding rock matrix and thus be more or less favourable to fluid circulation.

<footnote>*Keywords and phrases.* Adaptive constitutive law, elliptic equation, fractured porous media, non-linearity, variational formulation.</footnote>

<footnote n="1">¹ Department of Mathematics, Politecnico di Milano, p.za Leonardo da Vinci 32, Milano 20133, Italy.</footnote>

<footnote n="2">² IFP Energies nouvelles, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France.</footnote>

<footnote>*Corresponding author: alessio.fumagalli1984@gmail.com</footnote>
Depending on several aspects, e.g., the micro-structure and hydraulic aperture in the case of a fracture, the flow can be classified into different regimes corresponding to increasing flow rates and thus increasing mathematical and numerical difficulties. In many applications for low flow rates, caused by a combination of filling materials and a relatively low-pressure gradient, Darcy flow can be considered. Most of the research so far has been focused on this flow model; however, its validity is questionable and, to a large extent, it is insufficient for real problems. For increased flow rates, for example when a fracture is a variably open and narrow channel or the packaging of grains is too coarse, viscous effects become important and a Brinkman or Stokes equation is more coherent as a model to describe the flow; see [25,32]. Finally for high velocities, because of inertial effects, experiments show deviations from the previous models, which indicates the need for a non-linear correction term. For a fractured porous medium, the authors in [16,21] proposed a reduced model based on a Darcy–Forchheimer flow to capture this phenomenon. This effect is even more evident in large objects like faults, spanning several hundreds of meters.

Depending on its nature, a porous medium may exhibit all flow regimes in different regions separated by transition zones and coupled by suitable conditions. The transition zones can be located using the Forchheimer and Reynolds numbers (see [36]), which depend on the water velocity. This is the setting for the present paper, where the positions of the transition zones are not fixed at the outset of the problem: we obtain a (multi-physics) non-linear, free-boundary problem on the porous medium. We present a mathematical framework that is able to adapt the constitutive law in accordance with the flow regime. As a simplification, we assume that only two laws, i.e., two speed regimes, can be prescribed, leaving the case of multiple laws as a future work. A theoretical analysis is developed along with a numerical algorithm that tracks the transition zone. Theoretically, we show existence of solutions to our problem *via* the minimization of an underlying energy in the case when no derivatives are involved in the laws (thus excluding Stokes' and Brinkman's cases, which would require a different analysis left here for later work). This underlying energy is derived as it is classical in variational methods for PDEs; see for instance [26,34]. For related time-dependent problems involving the minimization of a free energy giving rise to a gradient flow structure, we refer the reader to [11,27]. Moreover, for a thorough overview of analytical aspects of porous medium equations, with and without time dependence, see [35], and for an example of existence and uniqueness theory of non-linear Darcy-like equations, see [5]. We do not show uniqueness since we only have that minimizers of the energy form a subset of the solution set of the original problem. We show that the convex nature of the problem is strongly intertwined with the “direction” of the jump between the high-speed and low-speed laws; indeed, if the permeability increases from low- to high-speed the problem is non-convex, otherwise it is convex. In the former case we are forced to restrict our existence result to one space dimension ($d = 1$). Numerically, the examples we give, first in an academic setting made of only one-dimensional domains and then in a simplified fracture network, illustrate the quality and applicability of the proposed algorithm, especially in the non-convex case. In the convex case, which is the analytically more favorable one, the algorithm often features oscillations between configurations which prevents it from converging. For the numerical simulations we used the library PorePy [20], a simulation tool written in Python for fractured and deformable porous media which is freely available on GitHub along with the numerical tests proposed in this work.

The paper is organized as follows. In Section 2 we describe our model and introduce the equations as well as some notation. In Section 3 we give the rigorous mathematical setting, along with the assumptions on the constitutive laws and the weak formulation of the problem. Section 4 is dedicated to our results on the existence of solutions and their proofs. Section 5 introduces the discrete formulation of the problem and a suitable algorithm to solve it. Numerical examples are reported in Section 6 for increasing geometrical and physical complexity. The work finishes with conclusions in Section 7.

## 2. PROPOSED MODEL

We identify the porous medium with an open, bounded, connected set $\Omega \subset \mathbb{R}^d$ with Lipschitz boundary $\partial\Omega$. We suppose $\Omega$ is filled with a fluid of constant density $\rho$. A general constitutive relation coupled with the mass
conservation in $\Omega$ reads

$$ \mathbf{\Lambda}(\mathbf{u}) = -\nabla p + \mathbf{f}, \qquad (2.1) $$

$$ \mathrm{div}\mathbf{u} = q, \qquad (2.2) $$

where $\mathbf{u}: \Omega \to \mathbb{R}^d$ is the unknown velocity of the fluid, $p: \Omega \to \mathbb{R}$ is the unknown pressure of the fluid, $q: \Omega \to \mathbb{R}$ is a given source term allowing, for instance, for fluid mass to be exchanged between fractures or wells and the surrounding porous medium, or rock matrix, and $\mathbf{f}$ is a given body force (e.g., gravity). Equation (2.1) is a law between the velocity field $\mathbf{u}$ and the pressure field $p$ via some operator $\mathbf{\Lambda}$. 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 = \bar{\Sigma}_v \cup \bar{\Sigma}_p$ and $\Sigma_v \cap \Sigma_p = \emptyset$. To (2.2) we add the following boundary conditions:

$$ \left\{ \begin{array}{ll} \mathbf{u} \cdot \mathbf{n} = u_0 & \text{on } \Sigma_v, \\ p = p_0 & \text{on } \Sigma_p, \end{array} \right. \qquad (2.3) $$

where $\mathbf{n}$ is the outward normal unit vector of $\partial\Omega$. Here, $u_0: \Sigma_v \to \mathbb{R}$ and $p_0: \Sigma_p \to \mathbb{R}$ are given functions setting the conditions on the boundary on $\mathbf{u}$ and $p$, respectively. We are denoting maps on $\Omega$ and their traces on the boundary $\partial\Omega$ by the same notation.

## 2.1. Velocity-pressure constitutive law

Examples of laws $\mathbf{\Lambda}$ relating velocity and pressure via (2.1) are

$$ \begin{aligned} \mathbf{\Lambda}_S(\mathbf{u}) &= -\mathbf{\Delta}\mathbf{u}, & \mathbf{\Lambda}_D(\mathbf{u}) &= \mathbf{K}^{-1}\mathbf{u}, \\ \mathbf{\Lambda}_B(\mathbf{u}) &= -\mathbf{\Delta}\mathbf{u} + \mathbf{K}^{-1}\mathbf{u}, & \mathbf{\Lambda}_{DF}(\mathbf{u}) &= \mathbf{K}^{-1}\mathbf{u} + \alpha\|\mathbf{u}\|^{r-2}\mathbf{u}, \end{aligned} \qquad (2.4) $$

where $\mathbf{K}: \Omega \to \mathbb{R}^{d \times d}$ is the permeability tensor and $\alpha \ge 0$ and $r \in [2, \infty)$ some parameters. The notation $\|\cdot\|$ stands for the Euclidean norm on $\mathbb{R}^d$. Choosing $\mathbf{\Lambda}_S$ gives Stokes' equation, $\mathbf{\Lambda}_D$ gives Darcy's equation, $\mathbf{\Lambda}_B$ gives Brinkman's equation and $\mathbf{\Lambda}_{DF}$ gives the generalized Darcy-Forchheimer equation (which simplifies into the classical Darcy-Forchheimer when $r = 3$).

To the authors' knowledge, a known issue that has not yet found a documented answer is the case when one needs to choose a combination of laws such as those given as examples in (2.4), rather than a single one, i.e., when one needs to couple different velocity-pressure laws according to some validity criterion which selects the better-adapted law. As already mentioned in the introduction, this criterion should depend on the speed regime (or Reynolds number) of the flow. For example, where the Reynolds number is low Darcy's law $\mathbf{\Lambda}_D$ may be preferred, whereas where it is high the Darcy-Forchheimer law $\mathbf{\Lambda}_{DF}$ might be a better choice. For this reason, in this paper we consider the case where we need to choose from two laws $\mathbf{\Lambda}_1$ and $\mathbf{\Lambda}_2$ according to some threshold speed $\bar{u} > 0$. We expect that generalizing our results to more than two laws (e.g., having a low-speed regime, a transitional regime and a high-speed regime) should not be difficult. In this setting, the law operator in (2.1) takes the form

$$ \mathbf{\Lambda}(\mathbf{u}) = \begin{cases} \mathbf{\Lambda}_1(\mathbf{u}) & \text{wherever } \|\mathbf{u}\| < \bar{u}, \\ \mathbf{\Lambda}_2(\mathbf{u}) & \text{wherever } \|\mathbf{u}\| > \bar{u}. \end{cases} \qquad (2.5) $$

Being able to impose a law at the transition zone $\{\|\mathbf{u}\| = \bar{u}\}$ is out of the scope of this paper – we will therefore consider our problem solved whenever we find a pressure field and a velocity field such that (2.2), (2.1) and (2.5) hold (see Rem. 3.6 for a consequence of this choice on the problem of uniqueness of solutions). This summarizes as in the following problem:

**Problem 2.1 (Strong form).** Find $\mathbf{u}: \Omega \to \mathbb{R}^d$ and $p: \Omega \to \mathbb{R}$ such that

$$ \left\{ \begin{array}{ll} \mathrm{div}\mathbf{u} = q & \text{on } \Omega, \\ \mathbf{\Lambda}(\mathbf{u}) = -\nabla p + \mathbf{f} & \text{on } \Omega, \\ \mathbf{u} \cdot \mathbf{n} = u_0 & \text{on } \Sigma_v, \\ p = p_0 & \text{on } \Sigma_p, \end{array} \right. $$
where $\Lambda$ is any law such as in (2.5).

Clearly, when $\Lambda_1 \neq \Lambda_2$, this choice of adaptable law introduces a discontinuity at $||u|| = \bar{u}$ which we shall treat carefully when studying existence of solutions. Note that this discontinuity is not one coming from an a priori physical subdivision of the domain, where different parts of $\Omega$ would be assigned different permeabilities, but rather one stemming a posteriori from the velocity field itself. This makes our setting contrast with that of inhomogeneous porous media, such as described by the Muskat problem in [18], by Darcy-Stokes interface coupling problems in [4, 7] or fracture barrier problems in [24]. In Section 3.1 we will complete the strong formulation in Problem 2.1, which is somewhat still formal since it lacks a law at the transition zone; indeed, we will introduce a multi-valued weak setting so as to be able to treat the transition zone without imposing any given law on it. We refer the reader to [12] for a multi-valued monotone operator approach for pressure-dependent permeabilities; note that the operator therein is continuous in velocity.

**Remark 2.2.** Whenever the boundary piece $\Sigma_p$ is such that $Vol^{d-1}(\Sigma_p) = 0$, where $Vol^{d-1}$ stands for the $(d-1)$-dimensional Lebesgue measure, in order to ensure uniqueness of the pressure field satisfying Problem 2.1 when $\Lambda$ is a classical continuous law, as opposed to (2.5), one imposes a constraint on the average of the pressure field:

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

for some $\bar{\varpi} \in \mathbb{R}$, normally assumed null. We shall therefore impose (2.6) whenever $Vol^{d-1}(\Sigma_p) = 0$, so that our problem in this case becomes: find functions $u$ and $p$ such that

$$ \left\{ \begin{array}{ll} \mathrm{div}u = q & \text{on } \Omega, \\ \Lambda(u) = -\nabla p + f & \text{on } \Omega, \\ u \cdot n = u_0 & \text{on } \partial\Omega, \end{array} \right. $$

under the constraint (2.6), where $\Lambda$ is now given in (2.5). For ease of discussion, however, we will often omit (2.6) and only refer to Problem 2.1 as being our problem, even when $Vol^{d-1}(\Sigma_p) = 0$; this average condition will nevertheless be naturally encoded in our weak formulation.

As we will see, we focus in this paper on Darcy-like operator laws, in the sense that they involve no derivatives of the velocity field, so that we need to exclude Stokes' and Brinkman's equations as admissible examples. Under some additional conditions (depending in particular on the "direction" of the jump between laws $\Lambda_1$ and $\Lambda_2$ at the transition zone), we show existence of solutions via the study of an energetic formulation of Problem 2.1. Indeed, we are able to define an energy functional on the space of velocity fields whose minimizers satisfy (some weak version of) Problem 2.1. Although we have uniqueness for this energetic formulation in some conditions (see Sect. 4), this property does not transfer to Problem 2.1 – this is related to the treatment of the transition zone, and one approach to circumvent this issue would be to show that the transition zone must have zero Lebesgue measure so that it would not play a role in defining weak solutions. Remark 2.3 below shows, unfortunately, that this is not the case in general. We will present in Remark 3.6 another potential approach to tackle uniqueness, but which we do not explore further in this paper.

**Remark 2.3.** It is fair at this point to ask whether a transition zone can ever be of non-zero Lebesgue measure and therefore be relevant in a weak setting. We give here a simple, but non-trivial, one-dimensional example illustrating that this can happen. Suppose $d = 1$, $\Omega = (0,1)$, $Vol^0(\Sigma_p) = 0$, $u_0 \equiv 0$ and $f \equiv 0$, and set the threshold speed to $\bar{u} = 1$. Note that in this case, the velocity part of a strong solution $(u,p)$ to Problem 2.1 is entirely determined independently of the chosen laws; indeed $u$ is the function defined by

$$ u(x) = \int_{0}^{x} q(y) \, dy \quad \text{for all } x \in (0, 1). $$
We therefore see that it is possible to choose a source term, which we recall in the fracture setting describes the fluid exchange between fracture and rock matrix, such that $u$ is equal to 1 on a subset of $(0, 1)$ with non-zero Lebesgue measure. Indeed, it is enough to pick

$$q(x) = \begin{cases} 4 & \text{for all } x \in (0, 0.25), \\ 0 & \text{for all } x \in [0.25, 0.75), \\ -4 & \text{for all } x \in [0.75, 1), \end{cases}$$

since this gives that the transition zone is $[0.25, 0.75]$. We thus see that a simple case can lead to transition zones of full dimension. (Note that in this strong solution setting with no-flux boundary conditions, we had to choose $q$ such as to satisfy the compatibility condition $\int_{\Omega} q = \int_{\partial\Omega} u_0 = 0$.)

## 2.2. Multiple sub-domains

In the case $\Omega$ is composed of multiple sub-domains (e.g., intersecting fracture branches), forming thus a network of equi-dimensional objects, we can extend the previous model by including suitable conditions at their intersections. Given $\Omega$ we introduce $\omega^i \subset \Omega$ to be a sub-domain, with $n_\omega \ni i$ their total number. Clearly, given two distinct $\omega^i$ and $\omega^j$, with $i \neq j$, we have $\dot{\omega}^i \cap \dot{\omega}^j = \emptyset$ and also that $\bar{\Omega} = \cup_{i=1}^{n_\omega} \bar{\omega}^i$.

We consider $2 \le n \le n_\omega$ sub-domains that meet at $\mathcal{I}$ whose closure $\bar{\mathcal{I}} = \cap_{i=1}^n \bar{\omega}^i$. To complete model in Problem 2.1 we impose the following classical conditions on $\mathcal{I}$:

$$\sum_{i=1}^{n} \mathbf{u}^{i} \cdot \mathbf{n}^{i} = 0 \quad \text{and} \quad p^{i} = p^{j}, \quad \forall i, j = 1, \dots, n, \qquad (2.7)$$

where with a superscript $i$ we indicate the corresponding object restricted to $\omega^i$, and $\mathbf{n}^i$ is a unit vector tangent to $\omega^i$ and pointing to $\mathcal{I}$, in the mono-dimensional case, and in addition normal to $\partial\omega^i$, in the multi-dimensional case. This condition is frequently used, see for instance [2, 3, 8, 9]. The first condition in (2.7) is a direct consequence of the conservation of mass at the intersection, while the second can be derived from each constitutive relation of the form (2.1). These conditions do not put any additional difficulties in the analysis and are therefore considered only in the numerical examples, while in the analysis we will keep assuming that $\Omega$ is a single domain.

It is possible to consider more complex conditions such as the so-called non-linear transmission conditions, which assume non-linear relations for the pressure. For more general conditions on $\mathcal{I}$, used for unsaturated and two-phase flow models, see for example [1, 13, 19, 22].

## 3. MATHEMATICAL SETTING

We give in this section the rigorous mathematical setting. In particular we introduce the assumptions on the underlying law operators as well as the weak formulation of our problem. We shall use the convention to use boldfaced symbols for vectors and vector-valued functions.

From now on, without loss of generality we take $\bar{u}$ to be equal to 1. For a given Lebesgue measurable field $u: \Omega \to \mathbb{R}^d$, we write

$$\Omega_1(\mathbf{u}) = \{\mathbf{x} \in \Omega \mid \|\mathbf{u}(\mathbf{x})\| < 1\}, \quad \Omega_2(\mathbf{u}) = \{\mathbf{x} \in \Omega \mid \|\mathbf{u}(\mathbf{x})\| > 1\}, \\ \Gamma(\mathbf{u}) = \Omega \setminus (\Omega_1(\mathbf{u}) \cup \Omega_2(\mathbf{u})) = \{\mathbf{x} \in \Omega \mid \|\mathbf{u}(\mathbf{x})\| = 1\}, \qquad (3.1)$$

where $\Omega_1(\mathbf{u})$, $\Omega_2(\mathbf{u})$ and $\Gamma(\mathbf{u})$ are what we have already respectively referred to as the low-speed region, high-speed region and transition zone (associated with $\mathbf{u}$). Note that thanks to the measurability assumption on $\mathbf{u}$, these three subsets are also measurable. Obviously the family $\mathcal{C} := \{\Omega_1(\mathbf{u}), \Omega_2(\mathbf{u}), \Gamma(\mathbf{u})\}$ forms a partition of
$\Omega$, and we will refer to $\mathcal{C}$ as the configuration of the problem, especially for the numerics in Sections 5 and 6. We can rewrite these sets in a more compact form:

$$ \Omega_1(\boldsymbol{u}) = \boldsymbol{u}^{-1}(B_1(\boldsymbol{0})), \quad \Omega_2(\boldsymbol{u}) = \boldsymbol{u}^{-1}(E_1(\boldsymbol{0})), \quad \Gamma(\boldsymbol{u}) = \boldsymbol{u}^{-1}(S_1(\boldsymbol{0})), $$

where $B_1(\boldsymbol{0})$ and $S_1(\boldsymbol{0})$ stand respectively for the unit open ball and unit sphere in $\mathbb{R}^d$ centered at the origin $\boldsymbol{0}$ and $E_1(\boldsymbol{0}) = \mathbb{R}^d \setminus (B_1(\boldsymbol{0}) \cup S_1(\boldsymbol{0}))$. Note that if $\boldsymbol{u}$ is not continuous, then $\Omega_1(\boldsymbol{u})$ and $\Omega_2(\boldsymbol{u})$ may not be open.

For all $n \in [1, \infty)$, $m \in (0, \infty)$ and $A \subset \mathbb{R}^d$ measurable we will denote by $L^n(A)$ and $W^{m,n}(A)$ the Lebesgue space of measurable functions on $A$ with integrable $n$th power and the $m$th-order Sobolev space associated to $L^n(A)$; we will also write $L^n(A)$ in place of $(L^n(A))^d$. As usual, in these spaces, equality between two functions is always intended in the almost-everywhere sense.

### 3.1. Assumptions on the velocity-pressure laws

In the following, the operator laws $\mathbf{\Lambda}_1$ and $\mathbf{\Lambda}_2$ are assumed to be of the form

$$ \mathbf{\Lambda}_1(\boldsymbol{u}) = \phi_1(\|\boldsymbol{u}\|^2)\boldsymbol{u} \chi_{\Omega_1(\boldsymbol{u})} \quad \text{and} \quad \mathbf{\Lambda}_2(\boldsymbol{u}) = \phi_2(\|\boldsymbol{u}\|^2)\boldsymbol{u} \chi_{\Omega_2(\boldsymbol{u})}, \qquad (3.2) $$

where $\chi_A$ is the characteristic function of any set $A \subset \mathbb{R}^d$. The functions $\phi_1, \phi_2: [0, \infty) \to [0, \infty)$ are continuous and increasing on $[0, 1]$ and $[1, \infty)$, respectively. Furthermore, $\phi_2$ satisfies the following assumption: there exist $r \ge 2$ and $c, C > 0$ such that

$$ ca^{\frac{r-2}{2}} \le \phi_2(a) \le C \left(1 + a^{\frac{r-2}{2}}\right) \quad \text{for all } a \ge 1. \qquad (3.3) $$

A recurrent notation we will use is

$$ \lambda_1 := \phi_1(1) \quad \text{and} \quad \lambda_2 := \phi_2(1), \qquad (3.4) $$

and will call the difference $\lambda_2 - \lambda_1$ the transition zone inverse permeability jump. Symmetrically, whenever $k_1 := 1/\lambda_1$ and $k_2 := 1/\lambda_2$ are considered (cf. in particular Sects. 5 and 6), the difference $k_2 - k_1$ will be called the transition zone permeability jump. We will see that the sign of this jump is an important feature which determines the convexity of the energy functional underlying the problem. Note that because $\phi_2$ is increasing on $[1, \infty)$, one must have $c \le \lambda_2$ and $C \ge \lambda_2/2$ in (3.3).

Interesting examples that fall into the above requests, in particular satisfying (3.2) with (3.3), are combinations of scalar versions of the Darcy and Darcy-Forchheimer laws $\mathbf{\Lambda}_D$ and $\mathbf{\Lambda}_{DF}$ (cf. (2.4)), as desired in the first place. Indeed, one is allowed to consider

$$ \mathbf{\Lambda}_1(\boldsymbol{u}) = \lambda_1 \boldsymbol{u} \chi_{\Omega_1(\boldsymbol{u})} \quad \text{and} \quad \mathbf{\Lambda}_2(\boldsymbol{u}) = \lambda_2 \boldsymbol{u} \chi_{\Omega_2(\boldsymbol{u})}, $$

that is, $\phi_1 \equiv \lambda_1$ and $\phi_2 \equiv \lambda_2$, or

$$ \mathbf{\Lambda}_1(\boldsymbol{u}) = \lambda_1 \boldsymbol{u} \chi_{\Omega_1(\boldsymbol{u})} \quad \text{and} \quad \mathbf{\Lambda}_2(\boldsymbol{u}) = (\lambda_{21} + \lambda_{22}\|\boldsymbol{u}\|)\boldsymbol{u} \chi_{\Omega_2(\boldsymbol{u})}. $$

that is, $\phi_1 \equiv \lambda_1$ and $\phi_2(a) = \lambda_{21} + \lambda_{22}\sqrt{a}$, where $\lambda_1, \lambda_2, \lambda_{21}$ and $\lambda_{22}$ are positive scalars. In the former case we would require $r = 2$, whereas in the latter $r = 3$.

**Remark 3.1.** This setting where $\mathbf{\Lambda}_1$ and $\mathbf{\Lambda}_2$ are as in (3.2) physically restricts us to scalar permeabilities. More general laws including tensor permeabilities, as motivated in [34], are for instance given by the following:

$$ \mathbf{\Lambda}_1(\boldsymbol{u}) = \phi_1(\lambda_1 \boldsymbol{u} \cdot \boldsymbol{u}) \lambda_1 \boldsymbol{u} \chi_{\Omega_1(\boldsymbol{u})} $$
and analogously for $\Lambda_2$, where $\lambda_1 \in \mathbb{R}^{d \times d}$ is a symmetric positive definite matrix encoding a tensor permeability. Even more general forms are envisageable:

$$ \mathbf{\Lambda}_1(\mathbf{u}) = \sum_{j=1}^{n} \phi_{1,j}(\mathbf{\Lambda}_1 \mathbf{u} \cdot \mathbf{u}) \mathbf{\Lambda}_1 \mathbf{u} \chi_{\Omega_1}(\mathbf{u}), $$

where the $\phi_{1,j}$ and $\mathbf{\Lambda}_1$ are n different law functions and permeability tensors and the dot $\cdot$ stands for the Euclidean inner product in $\mathbb{R}^d$. We leave these general laws to a future work. We claim that the techniques we use in the present paper for scalar laws should extend to tensor laws without too much difficulty.

We denote by $s = \frac{r}{r-1}$ the conjugate exponent of $r$. Thanks to the continuity of $\phi_1$ and the right-hand inequality in (3.3), we see that the operators $\mathbf{\Lambda}_1$ and $\mathbf{\Lambda}_2$ map $L^r(\Omega)$ into $L^s(\Omega)$. Because we allow for any law on the transition zone, this invites us to consider the multi-valued setting where our combined law $\mathbf{\Lambda}$ maps $L^r(\Omega)$ into the power set $2^{L^s(\Omega)}$ and is given by

$$ \mathbf{\Lambda}(\mathbf{u}) = \{\mathbf{\Lambda}_1(\mathbf{u}) + \mathbf{\Lambda}_2(\mathbf{u}) + h\chi_{\Gamma(\mathbf{u})}\}_{\mathbf{h} \in L^s(\Omega)}, \quad (3.5) $$

that is, $\mathbf{\Lambda}(\mathbf{u})$ is the set of all the functions of the form $\mathbf{\Lambda}_1(\mathbf{u}) + \mathbf{\Lambda}_2(\mathbf{u}) + h\chi_{\Gamma(\mathbf{u})}$, where $h \in L^s(\Omega)$. We highlight the fact that the transition zone law $h$ in the definition above really is a function in $L^s(\Omega)$ and not a trace on $\Gamma(\mathbf{u})$. In this way, $h\chi_{\Gamma(\mathbf{u})}$ is simply the restriction of an $L^s(\Omega)$ function to $\Gamma(\mathbf{u})$. This is also coherent with our energetic approach below, where we will link to $\mathbf{\Lambda}$ the subdifferential of a functional defined on $L^r(\Omega)$, whose topological dual is indeed $L^s(\Omega)$.

**Remark 3.2.** The existence results that we present in this paper still apply if we consider some “background” law with tensor permeability. Indeed, our proofs remain essentially untouched if the combined law $\mathbf{\Lambda}$ in (3.5) is replaced by

$$ \mathbf{\Lambda}(\mathbf{u}) = \{\beta(\mathbf{\lambda}\mathbf{u} \cdot \mathbf{u})\mathbf{\lambda}\mathbf{u} + \mathbf{\Lambda}_1(\mathbf{u}) + \mathbf{\Lambda}_2(\mathbf{u}) + h\chi_{\Gamma(\mathbf{u})}\}_{\mathbf{h} \in L^s(\Omega)}, $$

where $\mathbf{\lambda} \in \mathbb{R}^{d \times d}$ is symmetric positive definite and $\beta: [0, \infty) \to [0, \infty)$ is a continuous, increasing function such that $\beta + \phi_2$ satisfies (3.3) in place of $\phi_2$.

## 3.2. Weak formulation

We fix $q \in L^r(\Omega)$, $u_0 \in L^r(\Sigma_v)$, $p_0 \in W^{\frac{1}{r},s}(\Sigma_p)$ and $f \in L^s(\Omega)$. Consider the space

$$ \widetilde{W}^{1,s}(\Omega) = \begin{cases} \{\xi \in W^{1,s}(\Omega) \mid \frac{1}{|\Omega|} \int_{\Omega} \xi = \varpi\} & \text{if } \mathrm{Vol}^{d-1}(\Sigma_p) = 0, \\ \{\xi \in W^{1,s}(\Omega) \mid \xi = p_0 \text{ on } \Sigma_p\} & \text{if } \mathrm{Vol}^{d-1}(\Sigma_p) > 0, \end{cases} $$

where we recall that $\mathrm{Vol}^{d-1}$ is the $(d-1)$-dimensional Lebesgue measure. We give a first weak formulation of Problem 2.1:

**Problem 3.3 (Weak form I).** Find $(\mathbf{u}, p) \in L^r(\Omega) \times \widetilde{W}^{1,s}(\Omega)$ so that there exists $\mathbf{\Lambda}_u \in \mathbf{\Lambda}(\mathbf{u})$ such that

$$ \int_{\Omega} \mathbf{\Lambda}_u \cdot \varphi = - \int_{\Omega} \nabla p \cdot \varphi + \int_{\Omega} f \cdot \varphi \quad \forall \varphi \in L^r(\Omega), $$

$$ \int_{\Omega} \nabla \psi \cdot \mathbf{u} = - \int_{\Omega} q\psi + \int_{\Sigma_v} u_0 \psi \quad \forall \psi \in \widetilde{W}^{1,s}(\Omega), $$

where we recall that $\mathbf{\Lambda}$ is given in (3.5).
We set $f_0 = f$ in the case $Vol^{d-1}(\Sigma_p) = 0$ and $f_0 = f + \nabla(Ep_0)$ in the case $Vol^{d-1}(\Sigma_p) > 0$, with $E: W^{1,s}(\Sigma_p) \to W^{1,s}(\Omega)$ any extension operator being right-inverse of the $W^{1,s}(\Omega)$-trace operator, and by defining the Sobolev space

$$W_0^{1,s}(\Omega) = \begin{cases} \{\xi \in W^{1,s}(\Omega) \mid \int_{\Omega} \xi = 0\} & \text{if } \mathrm{Vol}^{d-1}(\Sigma_p) = 0, \\ \{\xi \in W^{1,s}(\Omega) \mid \xi = 0 \text{ on } \Sigma_p\} & \text{if } \mathrm{Vol}^{d-1}(\Sigma_p) > 0, \end{cases}$$

We endow $W_0^{1,s}(\Omega)$ with the norm $\|\psi\|_{W_0^{1,s}(\Omega)} = \|\nabla\psi\|_{L^s(\Omega)}$ for all $\psi \in W_0^{1,s}(\Omega)$. From the linearity in Problem 3.3 with respect to the pressure field, the formulation in Problem 3.3 is equivalent to the following one:

**Problem 3.4 (Weak form II).** Find $(u, p) \in L^r(\Omega) \times W_0^{1,s}(\Omega)$ so that there exists $\Lambda_u \in \Lambda(u)$ such that

$$\int_{\Omega} \Lambda_u \cdot \varphi = - \int_{\Omega} \nabla p \cdot \varphi + \int_{\Omega} f_0 \cdot \varphi \quad \forall \varphi \in L^r(\Omega),$$

$$\int_{\Omega} \nabla \psi \cdot u = - \int_{\Omega} q\psi + \int_{\Sigma_v} u_0 \psi \quad \forall \psi \in W_0^{1,s}(\Omega).$$

We emphasize here that the well-posedness of Problem 3.4 is not affected by the choice of the extension $E$ in the definition of $f_0$.

Thanks to Riesz's representation theorem, we shall equivalently manipulate functions in $L^s(\Omega)$ as elements of $(L^r(\Omega))^*$, the topological dual of $L^r(\Omega)$. In particular, this means that, given $u \in L^r(\Omega)$, operators $\Lambda_u$ in $\Lambda(u)$ will also be seen as maps from $L^r(\Omega)$ to $(L^r(\Omega))^* \sim L^s(\Omega)$. Denoting by $\langle \cdot, \cdot \rangle$ the dual mapping on $L^s(\Omega) \times L^r(\Omega)$, the weak formulation in Problem 3.4 can be equivalently rewritten as follows:

**Problem 3.5 (Weak form III).** Find $(u, p) \in L^r(\Omega) \times W_0^{1,s}(\Omega)$ so that there exists $\Lambda_u \in \Lambda(u)$ such that

$$\langle \Lambda_u, \varphi \rangle = - \langle \nabla p, \varphi \rangle + \langle f_0, \varphi \rangle \quad \forall \varphi \in L^r(\Omega),$$

$$\langle \nabla \psi, u \rangle = - \int_{\Omega} q\psi + \int_{\Sigma_v} u_0 \psi \quad \forall \psi \in W_0^{1,s}(\Omega).$$

Our goal is now to show existence for this weak formulation.

**Remark 3.6.** We repeat that uniqueness is a tougher question which we shall not answer in the present paper, but which will be part of a future work. Nonetheless, we give here a potential approach to handle it. The core point lies in the way our model is posed, in the fact that we are admitting any $L^s(\Omega)$ function as transition zone law (cf. h in (3.5)). This facilitates our task of finding solutions via the minimization of the underlying energy, since the subdifferential of the dissipation associated with the energy necessarily gives subsets of the law in (3.5) (as it is shown in the proof of Lem. 4.8), i.e., subsets of the admissible laws. However, at this point one could still have solutions with admissible transition zone laws which are not part of the minimizers of the underlying energy – this is why we cannot deduce uniqueness in our setting, even if the minimizers are unique and the energy is convex. An idea, then, would be to characterize the subderivatives at the transition zone and then update the transition zone law in (3.5) in a way as to match exactly the subdifferential of the dissipation.

## 4. EXISTENCE

We state and prove here our results on existence for Problem 3.5. As already mentioned, we will see that our results depend on the sign of the transition zone inverse permeability jump, $\lambda_2 - \lambda_1$. The strategy is the following:

(1) we reduce Problem 3.4 on the velocity and pressure fields into an equivalent problem on the velocity field only (cf. Problem 4.3 and [6, 33]);
(2) we derive an energetic formulation whose minimizers are solutions to this reduced problem on the velocity field (cf. Problem 4.9 and [34]);

(3) we study the existence of minimizers for this energetic formulation distinguishing the convex case $\lambda_1 \le \lambda_2$ from the non-convex case $\lambda_1 > \lambda_2$ (cf. Thms. 4.10 and 4.11). We are only able to treat the one-dimensional setting $d = 1$ when $\lambda_1 > \lambda_2$.

## 4.1. Reduction to a problem on the velocity field

Define the maps $B: L^r(\Omega) \to (W_0^{1,s}(\Omega))^*$ and $B^*: W_0^{1,s}(\Omega) \to (L^r(\Omega))^*$ by

$$B(\varphi)(\psi) = B^*(\psi)(\varphi) = \langle \nabla \psi, \varphi \rangle \quad \text{for all } \varphi \in L^r(\Omega) \text{ and } \psi \in W_0^{1,s}(\Omega), \quad (4.1)$$

Let us introduce the set

$$V = \{ \varphi \in L^r(\Omega) \mid \forall \psi \in W_0^{1,s}(\Omega), \langle \nabla \psi, \varphi \rangle = 0 \},$$

which satisfies $V = \text{Ker}(B)$. We write $V^\perp \subset L^s(\Omega)$ the polar subspace of $V$, that is, $V^\perp = \{g \in L^s(\Omega) \mid \forall \varphi \in V, \langle g, \varphi \rangle = 0\}$, and $L^r(\Omega)/V \subset L^r(\Omega)$ the quotient space of $L^r(\Omega)$ by $V$.

The following three lemmas are inspired from their equivalents in [6].

**Lemma 4.1.** The maps $B$ and $B^*$ in (4.1) are isomorphisms from $L^r(\Omega)/V$ to $(W_0^{1,s}(\Omega))^*$ and from $W_0^{1,s}(\Omega)$ to $V^\perp$, respectively.

*Proof.* By Lemma 2.1 of [6], it suffices to show that there exists $\gamma > 0$ such that

$$\inf_{\substack{\psi \in W_0^{1,s}(\Omega) \\ \varphi \neq 0}} \sup_{\substack{\varphi \in L^r(\Omega) \\ \varphi \neq 0}} \frac{\langle \nabla \psi, \varphi \rangle}{\|\psi\|_{W_0^{1,s}(\Omega)} \|\varphi\|_{L^r(\Omega)}} \ge \gamma.$$

To this end, let $\psi \in W_0^{1,s}(\Omega)$ with $\psi \neq 0$. Then, by identifying $L^s(\Omega)$ with $(L^r(\Omega))^*$ we get

$$\|\psi\|_{W_0^{1,s}(\Omega)} = \|\nabla \psi\|_{L^s(\Omega)} = \sup_{\substack{\varphi \in L^r(\Omega) \\ \varphi \neq 0}} \frac{\langle \nabla \psi, \varphi \rangle}{\|\varphi\|_{L^r(\Omega)}},$$

so that

$$1 = \frac{\|\psi\|_{W_0^{1,s}(\Omega)}}{\|\psi\|_{W_0^{1,s}(\Omega)}} = \sup_{\substack{\varphi \in L^r(\Omega) \\ \varphi \neq 0}} \frac{\langle \nabla \psi, \varphi \rangle}{\|\psi\|_{W_0^{1,s}(\Omega)} \|\varphi\|_{L^r(\Omega)}}.$$

Taking above the infimum over all $\psi \in W_0^{1,s}(\Omega)$ with $\psi \neq 0$ ends the proof. $\square$

Now, via the following two lemmas, we simplify our weak formulation in Problem 3.5 into a problem restricted to $V$ (cf. Problem 4.3 below).

**Lemma 4.2.** There exists a unique $\hat{u} \in L^r(\Omega)/V$ such that

$$\langle \nabla \psi, \hat{u} \rangle = - \int_{\Omega} q \psi + \int_{\Sigma_v} u_0 \psi \quad \text{for all } \psi \in W_0^{1,s}(\Omega).$$

*Proof.* For all $\psi \in W_0^{1,s}(\Omega)$, let

$$F(\psi) = - \int_{\Omega} q \psi + \int_{\Sigma_v} u_0 \psi,$$

and bound $F(\psi)$ as follows:

$$|F(\psi)| \le \|q\|_{L^r(\Omega)} \|\psi\|_{L^s(\Omega)} + \|u_0\|_{L^r(\Sigma_v)} \|\psi\|_{L^s(\partial\Omega)} \le (C \|q\|_{L^r(\Omega)} + \|u_0\|_{L^r(\Sigma_v)}) \|\psi\|_{W_0^{1,s}(\partial\Omega)},$$

where we used boundedness of the trace operator on $\psi$ and where $C > 0$ is a constant stemming from Poincar's
inequality. Recalling that indeed $q$ and $u_0$ are assumed to belong to $L^r(\Omega)$ and $L^r(\Sigma_v)$, respectively, the map $F$ is then a bounded linear map from $W_0^{1,s}(\Omega)$ into $\mathbb{R}$, i.e., $F \in (W_0^{1,s}(\Omega))^*$. By Lemma 4.1, we thus know there exists a unique $\hat{u} \in L^r(\Omega)/V$ such that

$$ \langle \nabla \psi, \hat{u} \rangle = B(\hat{u})(\psi) = F(\psi) = - \int_{\Omega} q\psi + \int_{\Sigma_v} u_0\psi \quad \text{for all } \psi \in W_0^{1,s}(\Omega), $$

which is the desired result.

**Problem 4.3** (Weak form restricted to $V$). Find $v \in V$ such that there exists $\Lambda_v \in \Lambda(v + \hat{u})$ satisfying

$$ \langle \Lambda_v, \varphi \rangle = \langle f_0, \varphi \rangle \quad \text{for all } \varphi \in V, $$

where $\hat{u}$ is given by Lemma 4.2.

**Lemma 4.4.** *Problem 3.5 is equivalent to Problem 4.3.*

*Proof.* We first suppose that $(u, p)$ is solution to Problem 3.5. Then we decompose $u$ as $u = (u - \hat{u}) + \hat{u} =: v + \hat{u}$. By Problem 3.5, there exists $\Lambda_v \in \Lambda(u) = \Lambda(v + \hat{u})$ such that

$$ \langle \Lambda_v, \varphi \rangle = \langle f_0, \varphi \rangle \quad \text{for all } \varphi \in V. $$

Furthermore, for all $\psi \in W_0^{1,s}(\Omega)$ we find $\langle \nabla \psi, v \rangle = \langle \nabla \psi, u \rangle - \langle \nabla \psi, \hat{u} \rangle = 0$, so that $v \in V$ and $v$ satisfies Problem 4.3.

Suppose now that $v \in V$ satisfies Problem 4.3. Then, $\Lambda_v - f_0 \in V^\perp$. By Lemma 4.1 we know $B^*$ is an isomorphism from $W_0^{1,s}(\Omega)$ to $V^\perp$, so that there exists a unique $p \in W_0^{1,s}(\Omega)$ with

$$ \langle \Lambda_v - f_0, \varphi \rangle = B^*(-p)(\varphi) = -\langle \nabla p, \varphi \rangle \quad \text{for all } \varphi \in L^r(\Omega). $$

Furthermore, writing $u = v + \hat{u}$ we get

$$ \langle \nabla \psi, u \rangle = \langle \nabla \psi, v \rangle + \langle \nabla \psi, \hat{u} \rangle = \langle \nabla \psi, \hat{u} \rangle = - \int_{\Omega} q\psi + \int_{\Sigma_v} u_0\psi, $$

since $v \in V$. We thus have that $(u, p)$ satisfies Problem 3.5, with $p$ unique.

## 4.2. Energetic formulation

We first give the definition of Fréchet subdifferential and strong local minimizer in our setting:

**Definition 4.5** (Fréchet subdifferential and strong local minimizer). Let $\mathcal{F}: V \to \mathbb{R}$. For all $v \in V$ we define the (Fréchet) subdifferential $\partial \mathcal{F}(v)$ of $\mathcal{F}$ at $v$ by:

$$ g_v \in \partial \mathcal{F}(v) \iff g_v \in V^* \text{ and } \forall \varphi \in V, \liminf_{\delta \to 0^+} \frac{\mathcal{F}(v + \delta\varphi) - \mathcal{F}(v)}{\delta} \geq \langle g_v, \varphi \rangle. $$

We say that $v \in V$ is a (strong) local minimizer of $\mathcal{F}$ if there exists $\eta > 0$ such that for all $\varphi \in V$ we have $\mathcal{F}(v + \delta\varphi) \geq \mathcal{F}(v)$ for all $\delta \in [0, \eta)$.

**Remark 4.6.** The important property of the subdifferential to keep in mind here is that if $v \in V$ is a local minimizer of a functional $\mathcal{F}: V \to \mathbb{R}$, then $0_{L^s(\Omega)} \in \partial \mathcal{F}(v)$. When $\mathcal{F}$ is convex, the reverse of this statement is also true: if $0_{L^s(\Omega)} \in \partial \mathcal{F}(v)$, then $V$ is a local minimizer of $\mathcal{F}$.
We write $\Psi: [0, \infty) \to \mathbb{R}$ the function given by

$$ \Psi(a) = \begin{cases} \Phi_1(a) & \text{for all } a \le 1, \\ \Phi_2(a) & \text{for all } a > 1, \end{cases} \qquad (4.2) $$

where $\Phi_1, \Phi_2: [0, \infty) \to \mathbb{R}$ are primitives of $\frac{\phi_1}{2}$ and $\frac{\phi_2}{2}$ (cf. (3.2)) such that $\Phi_1(1) = \Phi_2(1) = 0$. We define the dissipation $\mathcal{D}: V \to \mathbb{R}$ by

$$ \mathcal{D}(v) = \int_{\Omega} \Psi(\|v + \hat{u}\|^2) \quad \text{for all } v \in V. \qquad (4.3) $$

Thanks to our assumptions on $\phi_1$ and $\phi_2$ we easily get that the domain of $\mathcal{D}$ is indeed all of $V$, and it is thus a well-defined functional from $V$ into $\mathbb{R}$. Consider the following problem:

**Problem 4.7 (Dissipation form).** Find $v \in V$ such that there exists $g_v \in \partial\mathcal{D}(v)$ satisfying

$$ \langle g_v, \varphi \rangle = \langle f_0, \varphi \rangle \quad \forall \varphi \in V. $$

The following result holds:

**Lemma 4.8.** Any solution to Problem 4.7 is also solution to Problem 4.3.

*Proof.* Suppose that $v \in V$ is solution to Problem 4.7. Then we can pick $g_v \in \partial\mathcal{D}(v) \subset L^s(\Omega)$ so that Problem 4.7 holds. By definition, for all $\varphi \in V$ we must have

$$ \liminf_{\delta \to 0^+} \frac{\mathcal{D}(v + \delta\varphi) - \mathcal{D}(v)}{\delta} \ge \langle g_v, \varphi \rangle. $$

Write $V_1$ the subset of $V$ consisting of the functions which are supported in $\Omega_1(v + \hat{u})$. Then, for all $\varphi \in V_1$, since $-\varphi \in V_1$ as well, we get

$$ \liminf_{\delta \to 0^+} \frac{1}{\delta} \int_{\Omega_1(v+\hat{u})} \left( \Psi(\|v + \delta\varphi + \hat{u}\|^2) - \Phi_1(\|v + \hat{u}\|^2) \right) \ge \langle g_v, \varphi \rangle $$

and

$$ \liminf_{\delta \to 0^+} \frac{1}{\delta} \int_{\Omega_1(v+\hat{u})} \left( \Psi(\|v - \delta\varphi + \hat{u}\|^2) - \Phi_1(\|v + \hat{u}\|^2) \right) \ge -\langle g_v, \varphi \rangle, $$

which, by Lebesgue's dominated convergence theorem and the differentiability of $\Phi_1$, yield

$$ \int_{\Omega} \Lambda_1(v + \hat{u}) \cdot \varphi = \int_{\Omega_1(v+\hat{u})} \phi_1(\|v + \hat{u}\|^2) (v + \hat{u}) \cdot \varphi \ge \langle g_v, \varphi \rangle $$

and

$$ -\int_{\Omega} \Lambda_1(v + \hat{u}) \cdot \varphi = -\int_{\Omega_1(v+\hat{u})} \phi_1(\|v + \hat{u}\|^2) (v + \hat{u}) \cdot \varphi \ge -\langle g_v, \varphi \rangle. $$

All in all we get

$$ \langle g_v, \varphi \rangle = \langle \Lambda_1(v + \hat{u}), \varphi \rangle \quad \text{for all } \varphi \in V_1. $$

Similarly, denoting by $V_2$ the subset of $V$ consisting of the functions which are supported in $\Omega_2(v + \hat{u})$, we get

$$ \langle g_v, \varphi \rangle = \langle \Lambda_2(v + \hat{u}), \varphi \rangle \quad \text{for all } \varphi \in V_2. $$

Thus the function defined by

$$ \Lambda_v = \Lambda_1(v + \hat{u}) + \Lambda_2(v + \hat{u}) + g_v \chi_{\Gamma(v+\hat{u})} $$

is such that $\langle \Lambda_v, \varphi \rangle = \langle g_v, \varphi \rangle$ for all $\varphi \in V$. Therefore, from Problem 4.7 we obtain

$$ \langle \Lambda_v, \varphi \rangle = \langle f_0, \varphi \rangle \quad \text{for all } \varphi \in V. $$

Moreover $\Lambda_v \in \Lambda(v + \hat{u})$, where we recall $\Lambda$ is in (3.5). Hence $v$ is solution to Problem 4.3. $\square$
We now define the energy $\mathcal{E}: V \to \mathbb{R}$ associated to $\mathcal{D}$ by

$$ \mathcal{E}(v) = \mathcal{D}(v) - \langle \boldsymbol{f}_0, v + \hat{u} \rangle \quad \text{for all } v \in V. \qquad (4.4) $$

Let us write $M_{\mathcal{E}} \subset V$ the, possibly empty, set of local minimizers of $\mathcal{E}$. Consider the following minimization problem.

**Problem 4.9 (Energy minimization).** Find $v \in V$ such that $v \in M_{\mathcal{E}}$.

By Remark 4.6, any solution to Problem 4.9 is also solution to Problem 4.7; if $\mathcal{E}$ is convex these problems are actually equivalent. By Lemma 4.8 it therefore suffices to find a solution to Problem 4.9 in order to get the desired existence result on our original Problem 3.5. The rest of this section will thus be solely dedicated to solving Problem 4.9.

### 4.3. Case $\lambda_1 \le \lambda_2$

The result we wish to show here is the following:

**Theorem 4.10 (Existence and uniqueness when $\lambda_1 \le \lambda_2$).** Suppose that $\lambda_1 \le \lambda_2$. Then, Problem 4.9 has a unique solution.

*Proof.* Let us prove that the integrand $\Psi \circ \|\cdot\|^2$ (cf. (4.2)) of $\mathcal{D}$ is strictly convex. Define the functions $\bar{\Phi}_1, \bar{\Phi}_2: [0, \infty) \to \mathbb{R}$ by

$$ \bar{\Phi}_1(a) = \begin{cases} \Phi_1(a) & \text{for all } a \le 1, \\ \frac{\lambda_1}{2}(a-1) & \text{for all } a > 1, \end{cases} \quad \text{and} \quad \bar{\Phi}_2(a) = \begin{cases} \frac{\lambda_2}{2}(a-1) & \text{for all } a \le 1, \\ \Phi_2(a) & \text{for all } a > 1. \end{cases} $$

Then $\bar{\Phi}_1$ and $\bar{\Phi}_2$ are differentiable with

$$ \bar{\Phi}'_1(a) = \frac{1}{2} \begin{cases} \phi_1(a) & \text{for all } a \le 1, \\ \lambda_1 & \text{for all } a > 1, \end{cases} \quad \text{and} \quad \bar{\Phi}'_2(a) = \frac{1}{2} \begin{cases} \lambda_2 & \text{for all } a \le 1, \\ \phi_2(a) & \text{for all } a > 1. \end{cases} $$

Since $\phi_1$ and $\phi_2$ are increasing on $[0, 1]$ and $[1, \infty)$, respectively, and $\phi_1(1) = \lambda_1$ and $\phi_2(1) = \lambda_2$, the derivatives $\bar{\Phi}'_1$ and $\bar{\Phi}'_2$ are also increasing so that $\bar{\Phi}_1$ and $\bar{\Phi}_2$ are convex. Furthermore, because $\lambda_1 \le \lambda_2$ and $\phi_2$ is increasing, we have $\bar{\Phi}'_2(a) \ge \bar{\Phi}'_1(a)$ for all $a > 1$; thus, the fact that $\Phi_1(1) = \Phi_2(1)$ (and so $\bar{\Phi}_1(1) = \bar{\Phi}_2(1)$) yields $\bar{\Phi}_1(a) \le \bar{\Phi}_2(a)$ for all $a > 1$. Similarly, we get that $\bar{\Phi}_1(a) \ge \bar{\Phi}_2(a)$ for all $a < 1$. Let $a, b \in [0, \infty)$ and $t \in [0, 1]$. If $a, b \le 1$, then $(1-t)a + tb \le 1$ and

$$ \begin{aligned} \Psi((1-t)a + tb) &= \Phi_1((1-t)a + tb) = \bar{\Phi}_1((1-t)a + tb) \\ &\le (1-t)\bar{\Phi}_1(a) + t\bar{\Phi}_1(b) = (1-t)\Psi(a) + t\Psi(b), \end{aligned} $$

and similarly if $a, b > 1$. If now $a \le 1$, $b > 1$ and $(1-t)a + tb \le 1$, then we have

$$ \begin{aligned} \Psi((1-t)a + tb) &= \Phi_1((1-t)a + tb) = \bar{\Phi}_1((1-t)a + tb) \\ &\le (1-t)\bar{\Phi}_1(a) + t\bar{\Phi}_1(b) \le (1-t)\bar{\Phi}_1(a) + t\bar{\Phi}_2(b) = (1-t)\Psi(a) + t\Psi(b), \end{aligned} $$

and similarly if $(1-t)a + tb > 1$ or $a > 1$ and $b \le 1$. In all cases, we see that

$$ \Psi((1-t)a + tb) \le (1-t)\Psi(a) + t\Psi(b), $$

so that $\Psi$ is convex. Note that one could reach the same conclusion using that

$$ \Psi(a) = \max(\bar{\Phi}_1(a), \bar{\Phi}_2(a)) \quad \text{for all } a \ge 0. $$
Since $\phi_1$ and $\phi_2$ are positive on $(0, 1)$ and $[1, \infty)$, respectively, we get that $\Psi$ is increasing. Therefore, the function $\Psi \circ \|\cdot\|^2$ is strictly convex.

We now want to use the direct method of the calculus of variations to show that $\mathcal{E}$ has in fact a unique global (and thus local) minimizer. Let $(v_n)_{n \in \mathbb{N}} \subset V$ be a minimizing sequence for $\mathcal{E}$. Then we know there exists $N \in \mathbb{N}$ large enough and $K > 0$ such that $\mathcal{E}(v_n) < K$ for all $n > N$. Without loss of generality we can therefore assume that the sequence $(\mathcal{E}(v_n))_{n \in \mathbb{N}}$ is bounded by some constant $K > 0$. Hence, thanks to the left-hand inequality in (3.3), for all $n \in \mathbb{N}$ we have

$$ 
\begin{aligned}
K > \mathcal{E}(v_n) &= \int_{\Omega_1(v_n+\hat{u})} \Phi_1(\|v_n + \hat{u}\|^2) + \int_{\Omega_2(v_n+\hat{u})} \Phi_2(\|v_n + \hat{u}\|^2) - \int_{\Omega} f_0 \cdot (v_n + \hat{u}) \\
&\geq \Phi_1(0)|\Omega| + \frac{c}{r} \int_{\Omega_2(v_n+\hat{u})} \|v_n + \hat{u}\|^r - \|f_0\|_{L^s(\Omega)} \|v_n + \hat{u}\|_{L^r(\Omega)}^r \\
&\geq \left(\Phi_1(0) - \frac{c}{r}\right)|\Omega| + \frac{c}{r} \|v_n + \hat{u}\|_{L^r(\Omega)}^r - \|f_0\|_{L^s(\Omega)} \|v_n + \hat{u}\|_{L^r(\Omega)}^r,
\end{aligned}
$$

which shows that the sequence $(\|v_n\|_{L^r(\Omega)})_{n \in \mathbb{N}}$ is bounded. Thus we can extract a subsequence from $(v_n)_{n \in \mathbb{N}}$, still denoted $(v_n)_{n \in \mathbb{N}}$ which converges weakly to some $v \in L^r(\Omega)$. Since further $V$ is weakly closed we in fact have $v \in V$. Because $\Psi \circ \|\cdot\|^2$ is convex, the dissipation $\mathcal{D}$ is weakly lower semi-continuous and so is the energy $\mathcal{E}$. We therefore yield

$$ \inf_{w \in V} \mathcal{E}(w) = \liminf_{n \to \infty} \mathcal{E}(v_n) \geq \mathcal{E}(v) \geq \inf_{w \in V} \mathcal{E}(w), $$

so that $\mathcal{E}(v) = \inf_{w \in V} \mathcal{E}(w)$ and $v$ is a global minimizer of $\mathcal{E}$ and so $v \in M_{\mathcal{E}}$.

To show that $M_{\mathcal{E}}$ is a singleton it is enough to prove that $\mathcal{E}$ is strictly convex. Let $v, w \in V$ and $t \in (0, 1)$. Suppose furthermore that $v \neq w$ and write $A \subset \Omega$ the set where $v$ and $w$ are different; the Lebesgue measure of $A$ is therefore positive. Note that we must have $\Psi(a) \neq 0$ for all $a \neq 1$. Using the strict convexity of $\Psi \circ \|\cdot\|^2$ we therefore get

$$ 
\begin{aligned}
\mathcal{D}((1-t)v + tw) &= \int_{\Omega} \Psi(\|(1-t)v + tw + \hat{u}\|^2) = \int_{\Omega} \Psi(\|(1-t)(v + \hat{u}) + t(w + \hat{u})\|^2) \\
&= (1-t) \int_{\Omega \setminus A} \Psi(\|v + \hat{u}\|^2) + t \int_{\Omega \setminus A} \Psi(\|w + \hat{u}\|^2) \\
&\quad + \int_{A} \Psi(\|(1-t)(v + \hat{u}) + t(w + \hat{u})\|^2) \\
&< (1-t) \int_{\Omega \setminus A} \Psi(\|v + \hat{u}\|^2) + t \int_{\Omega \setminus A} \Psi(\|w + \hat{u}\|^2) \\
&\quad + (1-t) \int_{A} \Psi(\|v + \hat{u}\|^2) + t \int_{A} \Psi(\|w + \hat{u}\|^2) \\
&= (1-t) \int_{\Omega} \Psi(\|v + \hat{u}\|^2) + t \int_{\Omega} \Psi(\|w + \hat{u}\|^2) = (1-t)\mathcal{D}(v) + t\mathcal{D}(w),
\end{aligned}
$$

so that the dissipation $\mathcal{D}$ is strictly convex. Consequently, the energy $\mathcal{E}$ is also strictly convex and the proof is over. $\square$

### 4.4. Case $\lambda_1 > \lambda_2$

This case is more difficult to tackle than the case $\lambda_1 \le \lambda_2$. Indeed, we lose the convexity of the integrand $\Psi \circ \|\cdot\|^2$ (cf. proof of Thm. 4.10), which by Tonelli's theorem of functional analysis means that $\mathcal{E}$ is not weakly lower semi-continuous. As a consequence we cannot use the direct method of the calculus of variations in the weak topology to deduce the existence of a minimizer of $\mathcal{E}$.
To simplify our task at this point, we shall restrict to the one-dimensional case (i.e., $d = 1$) and leave the higher-dimensional case for future investigation. The results we present here therefore apply to the numerical experiments we present below. The one-dimensional case is simpler since we can easily characterize the space $V$ depending on the boundary conditions, as will be clear from the proof of the following theorem:

**Theorem 4.11** (Existence and uniqueness when $\lambda_1 > \lambda_2$). Let $d = 1$ and suppose that $\lambda_1 > \lambda_2$. If $\mathrm{Vol}^0(\Sigma_v) = 0$, then Problem 4.9 has a solution. If instead $\mathrm{Vol}^0(\Sigma_v) > 0$, then Problem 4.9 has a unique solution.

*Proof.* Without loss of generality, take $\Omega = (0, 1)$.

**Case** $\mathrm{Vol}^0(\Sigma_v) = 0$. First note that here $\Sigma_p = \{0, 1\}$. Let us characterize $V$ in this case. To this end, note that any $\eta \in V$ with $\int_0^1 \eta = 0$ has a primitive in $W_0^{1,s}((0, 1))$; indeed, the function $\psi: (0, 1) \to \mathbb{R}$ defined by

$$ \psi(x) = \int_0^x \eta \quad \text{for all } x \in (0, 1) $$

satisfies $\psi' = \eta$ and $\psi(0) = \psi(1) = 0$. Let now $v \in V$ and define $\eta \in V$ as

$$ \eta(x) = v(x) - \int_0^1 v \quad \text{for all } x \in (0, 1), $$

so that obviously $\int_0^1 \eta = 0$. Write $\psi \in W_0^{1,s}((0, 1))$ a primitive of $\eta$ and compute, for all $\varphi \in V$,

$$ \int_0^1 \left( v - \int_0^1 v \right) \varphi = \int_0^1 \eta \varphi = \int_0^1 \psi' \varphi = 0. $$

Thus $v = \int_0^1 v$ and $v$ is constant (almost everywhere). Since constant functions clearly belong to $V$, this shows that $V$ is in fact the set of all constant functions on $(0, 1)$ and we can identify $V$ with $\mathbb{R}$.

This in particular means that the energy $\mathcal{E}$ defined in (4.4) can be identified with the following function $E: \mathbb{R} \to \mathbb{R}$:

$$ E(\alpha) = \int_0^1 \Psi((\alpha + \hat{u})^2) - \alpha \int_0^1 f_0 \quad \text{for all } \alpha \in \mathbb{R}, $$

where we recall that $\Psi$ is given in (4.2) and $\hat{u} := \hat{u}$ is as in Lemma 4.2. We can use the direct method of the calculus of variations on $E$ in the Euclidean topology in $\mathbb{R}$. Let $(\alpha_n)_{n \in \mathbb{N}} \subset \mathbb{R}$ be a minimizing sequence for $E$. Using a similar calculation as in the proof of Theorem 4.10, thanks to (3.3) we can show that $(\alpha_n)_{n \in \mathbb{N}}$ is bounded in $\mathbb{R}$. Therefore, there exists a subsequence of $(\alpha_n)_{n \in \mathbb{N}}$, still denoted by $(\alpha_n)_{n \in \mathbb{N}}$, converging to some $\alpha \in \mathbb{R}$. By continuity of $\Psi$ and Fatou's lemma we get that $E$ is lower semi-continuous on $\mathbb{R}$, so that

$$ \inf_{\beta \in \mathbb{R}} E(\beta) = \liminf_{n \to \infty} E(\alpha_n) \ge E(\alpha) \ge \inf_{\beta \in \mathbb{R}} E(\beta), $$

showing that $\alpha$ is a global minimizer of $E$ and the constant function $\alpha$ belongs to $M_{\mathcal{E}}$.

**Case** $\mathrm{Vol}^0(\Sigma_v) > 0$. Note that here $\Sigma_p \in \{\emptyset\} \cup \{\{0\}\} \cup \{\{1\}\}$. Similarly to the previous case, let us characterize $V$. Note that any $v \in V$ has a primitive in $W_0^{1,s}((0, 1))$; indeed, the function $\psi: (0, 1) \to \mathbb{R}$ defined for all $x \in (0, 1)$ by

$$ \left\{ \begin{array}{ll} \psi(x) = \int_0^x v - \int_0^1 \left( \int_0^y v \right) dy & \text{if } \Sigma_p = \emptyset, \\ \psi(x) = \int_0^x v - \int_0^1 v & \text{if } \Sigma_p = \{1\}, \\ \psi(x) = \int_0^x v & \text{if } \Sigma_p = \{0\}, \end{array} \right. $$
satisfies $\psi' = \eta$, and $\int_0^1 \psi = 0$ if $\Sigma_p = \emptyset$ and $\psi(x) = 0$ for $x \in \Sigma_p$ otherwise. Let $v \in V$, write $\psi$ a primitive of $v$ and compute, for all $\varphi \in V$,

$$ \int_0^1 v\varphi = \int_0^1 \psi'\varphi = 0, $$

so that $v = 0$. This shows that $V = \{0\}$, i.e., $V$ contains only the zero function on $(0, 1)$. Trivially then, the zero function is the unique global minimizer of the energy $\mathcal{E}$ in (4.4) and $M_{\mathcal{E}}$ is a singleton. Note that in this case the energy is only trivially convex. $\square$

# 5. NUMERICAL APPROXIMATION

In this part we present the numerical approximation for the considered problem when $d = 1$; we will still keep boldfaced notation for vectors and vector-valued functions for coherence with most of the previous sections. In particular, in Section 5.1 the algorithm to track the transition zone is described. Later in Section 5.2 we describe the discretization adopted for a known transition zone. As mentioned in the introduction, for the implementation we have used the flexible framework of the PorePy library; see [20].

We introduce a mesh $\Omega_h$ composed of non-overlapping segments $E \in \Omega_h$ that approximate $\Omega$; we clearly have $\bar{\Omega}_h = \cup_{E \in \Omega_h} \bar{E}$. At the discrete level we can define the approximate configuration $C_h$ which is given by the set $C_h = \{\Omega_{1,h}, \Omega_{2,h}, \Gamma_h\}$, where respectively $\Omega_{1,h} \subset \Omega_h$ and $\Omega_{2,h} \subset \Omega_h$ are the approximations of the domains $\Omega_1$ and $\Omega_2$. $\Gamma_h$ is the approximation of the transition zone $\Gamma$. Note that for simplicity we are dropping the dependence of these regions on the velocity field. For each element $E$ we name $e_1$ and $e_2$ its two extremal vertices and $h_E$ its length. We let $h = \max_{E \in \Omega_h} h_E$ be the mesh size. In the case of a fracture network where the fracture are represented by segments, we suppose that each intersection is respected by the grid. We indicate by $(u_h, p_h)$ the approximation of $(u, p)$ for a given mesh.

## 5.1. Transition zone tracking algorithm

We suppose that the equations are discretized with a numerical scheme that gives an accurate enough velocity field. We consider an iterative scheme such that for each step $i$ a tentative configuration $C_h^{(i)}$ approximates $C_h$. For a given configuration $C_h^{(i-1)}$, the proposed algorithm solves the differential problem obtaining a new velocity field $u^{(i)}$. To speed up the computation, we evaluate condition (2.5) only at the extremities of each grid element $E$. If we obtain opposite values, we can thus determine the position of a transition zone in the considered element up to a given tolerance $\epsilon_\Gamma$. This part can be coded with an embarrassingly parallel workload, speeding up the algorithm. Having checked all the elements and computed the new configuration $C_h^{(i)}$, the algorithm resets with the new configuration as starting point. The exit strategy considers the position of $\Gamma_h$ for two successive iterations: if their distance is smaller than a given threshold $\epsilon_\Omega$, then a stable configuration is reached and the algorithm ends. We summarize in Algorithm 1 the implemented scheme.

We see that with this algorithm we cannot locate multiple transition zones in a single element $E$; this is a direct implication of the fact that $\Omega_{1,h}^{(i)}$ or $\Omega_{2,h}^{(i)}$ would be smaller than the mesh size in this case. The algorithm is also unable to track transitions zone which are full-dimensional (in this case of space dimension 1) since we choose to place only one point as the transition zone within an element whose vertices have velocities around the threshold speed.

To avoid unnecessary loops due to the chosen accuracy, numerical experiments showed that a good value of $\epsilon_\Omega$ is comparable with the mesh size given at the outset of the problem. Also, the value of $\epsilon_\Gamma$ is set quite small to determine the transition zones with high precision. Numerical examples when the transition zone inverse permeability jump is non-positive indicated convergence of the proposed scheme, independently from the starting configuration, thus suggesting that the problem has a unique solution.

**Remark 5.1.** Since the law on the transition zone is not defined, we can simply assume conservation of fluxes and continuity of pressure across $\Gamma$.
**Algorithm 1:** Transition zone tracking algorithm.

**Data:** For $i = 1$: the configuration $C_h^{(i-1)}$ and the tolerances $(\epsilon_\Gamma, \epsilon_\Omega)$;

**Result:** For $i = \text{end}$: the configuration $C_h^{(i)}$;

do

$$ (\mathbf{u}^{(i)}, p^{(i)}) \leftarrow \text{PDE}(C_h^{(i-1)}); $$

for $E = (e_1, e_2) \in \Omega_h^{(i-1)}$ do

$$ (c_1, c_2) \leftarrow (\|\mathbf{u}^{(i)}(e_1)\| < 1, \|\mathbf{u}^{(i)}(e_2)\| < 1); $$

if $c_1 \neq c_2$ then

$$ C_h^{(i)} \text{ in } E \leftarrow \text{Transition_zone}(\mathbf{u}^{(i)}); $$

end

end

$$ d \leftarrow \text{Distance}(\Gamma_h^{(i)}, \Gamma_h^{(i-1)}); $$

$$ i \leftarrow i + 1; $$

while $d > \epsilon_\Omega$;

## 5.2. Discretization in space

Following the algorithm discussed before, we assume that a configuration $C_h^{(i)}$ is given. In the case of non-linear constitutive relations $\Lambda$, an iterative scheme can be used, e.g., fixed-point or Newton or L-scheme; for the last one, see [23, 28, 29] applied to Richards' equation and two-phase flow in porous media. The solution is computed up to a tolerance $\epsilon_{nl}$. In our implementation we have considered the fixed-point iteration scheme. For this, to discuss the numerical approximation, we assume thus a linear constitutive relation $\Lambda$.

Following [15], to solve the problem we consider the mixed finite element approximation of lowest-order degree. For a given mesh $\Omega_h$ we have $(\mathbf{u}_h, p_h) \in (\mathbb{P}_0(\Omega_h), \mathbb{RT}_0(\Omega_h))$, where the first is the space of constant piecewise polynomial and the second the Raviart-Thomas space; see [30, 31]. This pair of discrete spaces is stable and gives a good approximation for the velocity field, essential for our purposes. The resulting discrete problem is well-posed, see also [10, 15], and the discrete solution converges to the exact one as $h$ goes to zero.

## 6. NUMERICAL EXAMPLES

In this part, we present numerical evidence for the quality and effectiveness of the previously introduced framework. Particularly for increasing geometrical and physical complexity. We consider four cases, starting from a realistic example comparing the proposed model with a fixed flow regime in Section 6.1. Then we study the cases of a single domain in Section 6.2.1, two crossing domains in Section 6.2.2, and the fracture network of Benchmark 1 from [14] in Section 6.2.3.

These examples were developed with the open source library PorePy [20]. The associated scripts are freely accessible. PorePy uses Gmsh [17] to construct the grids.

### 6.1. Comparison with a fixed flow regime model

We consider a realistic setting and compare the proposed approach with an *a-priori* choice of the flow regimes based on physical considerations. The problem design is a rather simple one, where we need to guess the position of the transition zone. Nevertheless, we will see in the next examples that, for a general problem, this might be a challenging question, thus supporting the need of an adaptive strategy.

First, we rewrite the Darcy and Darcy-Forchheimer equations as

$$ \mu k^{-1} \mathbf{u} = -\nabla p \quad \text{and} \quad \mu k^{-1} \mathbf{u} + \rho \beta \|\mathbf{u}\| \mathbf{u} = -\nabla p, $$
where $\mu$ is the viscosity of the liquid in [Pa s], $k$ the scalar intrinsic rock permeability in [m$^2$], $\rho$ the liquid density [kg m$^{-3}$] and $\beta$ the so-called Forchheimer coefficient in [m$^{-1}$]. We consider water as liquid so $\mu = 10^{-3}$[Pa s] and $\rho = 10^3$[kg m$^{-3}$].

Following [36], we can define the Forchheimer number $F_o$ to discriminate when non-Darcy effects should be included in the model or not, namely,

$$F_o = \frac{k\beta\rho\|\mathbf{u}\|}{\mu};$$

for $F_o < 0.1$ a pure Darcy model can be considered, while for $F_o \ge 0.1$, because of inertial effects, the Darcy–Forchheimer model should be used instead. The domain is of length 100[m] with permeability given by

$$k(x) = \begin{cases} 5 \cdot 10^{-11}[\text{m}^{-2}] & \text{for } x \le 50, \\ 5 \cdot 10^{-12}[\text{m}^{-2}] & \text{for } x > 50. \end{cases}$$

A pressure gradient is set on the boundary with $p(x = 0[\text{m}]) = 10^6[\text{Pa}]$ and $p(x = 100[\text{m}]) = 2 \cdot 10^6[\text{Pa}]$, a scalar source term $q = 10^{-6}[\text{s}^{-1}]$ is set in the interval (30[\text{m}], 70[\text{m}]) and the Forchheimer coefficient is taken as $\beta = 10^8[\text{m}^{-1}]$. We can estimate a fluid velocity to be in the order of $10^{-4}$ and with this value the Forchheimer number is estimated as

$$F_o(x) = \begin{cases} 0.5 & \text{for } x \le 50[\text{m}], \\ 0.05 & \text{for } x > 50[\text{m}]. \end{cases}$$

Thus the Darcy–Forchheimer model should be assumed in the first part of the domain and the Darcy model in the second. This is of course an estimate since $F_o$ depends on the fluid velocity which is an unknown, so that we cannot guarantee that this choice correspond to the actual position of the transition zone.

We apply the proposed model and compare the result with the a-priori set flow regime. For the latter the transition is at $x = 50[\text{m}]$ while for the former we set $\bar{u} = 10^{-4}[\text{m s}^{-1}]$. Figure 1 compares the two solutions obtained.

We notice that, even if the difference are not macroscopic, the transition zone between the two regimes is located at $x \approx 48.2[\text{m}^{-1}]$, showing that the actual Darcy model should be considered before what is assumed in the fixed flow regime setting. We have about 4% relative error in the determination of the transition zone. Also, the relative errors for the pressure and velocity show a difference mostly located at the centre of the domain.

Even in this simple setting, we see the importance of accurately choosing the appropriate regions for each flow regime. For more complex cases, their a-priori dislocation can be even more challenging.

## 6.2. The adaptive model

In all remaining cases, linear and non-linear velocity-pressure relations are considered as well as the impact of a possible vector source term $f$. The threshold on the velocity norm is set as $\bar{u} = 0.15$, so that we have

$$\Omega_1 = \{\mathbf{x} \in \Omega \mid \|\mathbf{u}(\mathbf{x})\| < \bar{u}\}, \quad \Omega_2 = \{\mathbf{x} \in \Omega \mid \|\mathbf{u}(\mathbf{x})\| > \bar{u}\}, \quad \Gamma = \{\mathbf{x} \in \Omega \mid \|\mathbf{u}(\mathbf{x})\| = \bar{u}\}.$$

If not otherwise specified, we consider a mesh size equal to $h = 5 \cdot 10^{-2}$ as well as $\epsilon_\Omega = h$; the maximum number of iterations to reach a stable configuration with Algorithm 1 is set to 50; and to track the transition zone we set $\epsilon_\Gamma = 10^{-10}$. In all the cases, the initial configuration is chosen to be $\Omega_1^{(0)} = \Omega$, i.e., the whole domain coincides with the low-speed region.

Our simulations, unless mentioned otherwise, are based on the following combinations of laws. In the linear case, we consider a combination of classical Darcy laws between the velocity and pressure, namely $K^{-1}\mathbf{u} = -\nabla p + \mathbf{f}$ with $K = kI$, where $I$ is the identity matrix and

$$k = \begin{cases} k_1 = 1 & \text{in } \Omega_1, \\ k_2 = 10 & \text{in } \Omega_2. \end{cases} \qquad (6.1)$$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Solutions for the problem of Section 6.1. On the left, using the proposed adaptive model; in the center, with the fixed flow regime; on the right, the error between them.</figcaption></figure>

(With the notation used previously, $k_i = 1/\lambda_i$, $i \in \{1, 2\}$.) Specifically, we have

$$ \mathbf{\Lambda}(\mathbf{u}) = \begin{cases} \mathbf{u} & \text{in } \Omega_1, \\ 0.1\mathbf{u} & \text{in } \Omega_2. \end{cases} \qquad (6.2) $$

When the non-linear case is studied, a non-linear and heterogeneous relationship between the velocity and the pressure is set. We assume a linear Darcy flow in $\Omega_1$ and a Darcy-Forchheimer flow in $\Omega_2$. Specifically, we have

$$ \mathbf{\Lambda}(\mathbf{u}) = \begin{cases} \mathbf{u} & \text{in } \Omega_1, \\ (0.01 + 3\|\mathbf{u}\|)\mathbf{u} & \text{in } \Omega_2. \end{cases} \qquad (6.3) $$

If not specified, we consider 50 as maximum for the number of iterations of the non-linear solver with tolerance $\epsilon_{nl}$ equal to $10^{-4}$.

### 6.2.1. Single domain

In this first case, we consider $\Omega = (0, 1)$. Boundary conditions are set to zero for the pressure. In the sequel, we consider a linear and non-linear law relationship between $\mathbf{u}$ and $p$. The scalar and vector source terms are set equal to

$$ q(\mathbf{x}) = \begin{cases} 1 & \text{if } x_1 \le 0.3, \\ -1 & \text{if } 0.3 < x_1 < 0.7, \\ 1 & \text{if } x_1 \ge 0.7, \end{cases} \quad \text{and} \quad \mathbf{f} = [5 \cdot 10^{-2}, 0, 0]^\top. $$

**Linear case.** In this part, we consider the linear case with $\mathbf{\Lambda}$ as in (6.2). The numerical solution is reported in Figure 2 along with some snapshots of the tentative solutions from Algorithm 1. What the “if $\Omega_1$” legend represents is a binary outcome saying if the region is $\Omega_1$ or not. The “condition” legend represents the configuration at the previous algorithm iteration. The first figure gives the initial condition imposed, the second the configuration after 3 steps and the third the final solution (at iteration 6). We notice the creation of multiple transitions zone $\Gamma$, which might change during the determination of the final configuration. By changing the
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Solution for different iterations for the linear problem of Section 6.2.1. From the left, at iterations 0, 2 and 6. The green line represents if the portion of $\Omega$ that belongs to $\Omega_1$ or not. The pressure profile is amplified by a factor of 10.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Number of iterations of Algorithm 1 for different values of the inverse of $k_2$ (that is, $\lambda_2$) for the linear example in Section 6.2.1. The red asterisks mean that the maximum number of iterations is reached.</figcaption></figure>

initial condition we get to the same stationary solution, and by refining the grid we obtain a stable outcome similar to the one presented in Figure 2.

In Figures 3 and 4 we study a case different from (6.2), where we test our algorithm for various permeabilities $k_2$ while fixing $k_1 = 1$. By increasing the maximum number of iterations to 1000, Figure 3 on the left shows the impact of changing $k_2$ on the number of iterations for the Algorithm 1. We see that for $k_2 \ge k_1$ the solution is always computable while we cannot draw the same conclusion for $k_2 < k_1$.

An example for $k_2 = 0.25$ is shown in Figure 4, where we notice that “if $\Omega_1$” and “condition” are perfectly flipped between any two successive iterations. The algorithm “jumps” between two states and thus does not converge.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. The two solutions computed by the algorithm for $k_2 = 0.25$; see the linear case in Section 6.2.1.</figcaption></figure>

By setting $f = 0$ and $p(1) = 0.2$, the latter to avoid that the algorithm converges in 1 iteration for each value of $k_2$, we can do the same analysis and obtain the plot in Figure 3 on the right. We deduce that the presence of the vector source has an impact on the computability of the solution when $k_2 > k_1$, that is, when the transition zone permeability jump is positive.

We can conclude that, even in this simple setting, the obtained numerical evidence is interesting and gives a valid support to the developed theory, at least when $k_2 \ge k_1$.

**Non-linear case.** In this part, we consider the non-linear case with $\Lambda$ as in (6.3). Figure 5 shows the numerical solution obtained at different iteration steps of Algorithm 1. The stable solution is reached very quickly and only two iterations of the outer scheme are needed. We see the effect on the pressure of the non-linear law, which changes shape between the initial configuration and iteration 1. By changing the parameters, it is possible to show that the obtained solution is independent from the initial condition and stable with respect the grid refinement, once the mesh size is small enough to separate close transition zones. In this case 2 iterations are needed to reach the stable solution, the first requires only 1 iteration and the second 6 iterations.

We consider now the effect of the tolerance $\epsilon_{nl}$ imposed in the non-linear solver, in particular its effect on the number of iterations and resulting error. By keeping fixed the spatial discretization, we compute a reference solution ($u_{ref}, p_{ref}$) with tolerance $\epsilon_{nl} = 10^{-12}$. We report in Table 1 the comparison with higher tolerances. In all cases the first iteration requires only two non-linear cycles, since the initial configuration has only a linear problem. At the second iteration the non-linear steps depend on the chosen tolerance, this value is reported in the table. In particular, we notice that both errors $err_p$ for $p$ and $err_u$ for $u$ computed as

$$err_p = \frac{\|p_{ref} - p\|}{\|p_{ref}\|} \quad \text{and} \quad err_u = \frac{\|u_{ref} - u\|}{\|u_{ref}\|}$$

have a monotone decay. The norms in the previous expression are the Euclidean norms of the solution vector. The error is rather small and decays very quickly, reaching zero for the non-linear tolerance equal to $10^{-3}$. The outer iterations are not influenced by this parameters, probably due to the small errors obtained in all the cases.

Also in this case we can conclude that, even in this simple setting, the obtained numerical evidence is interesting and gives a valid support to the developed theory.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Solution for different iterations for the non-linear problem of Section 6.2.1. From the left, at iterations 0 and 1. The green line represents if the portion of the $\Omega$ that belongs to $\Omega_1$ or not. The pressure profile is amplified by a factor of 100.</figcaption></figure>

<table><caption>TABLE 1. Numbers of iterations and errors computed for the example in Section 6.2.1, non-linear case.</caption><thead><tr><th>ε<sub>n1</sub></th><th>it<sub>out</sub></th><th>it<sub>in</sub></th><th>err<sub>p</sub></th><th>err<sub>u</sub></th></tr></thead><tbody><tr><td>10<sup>-1</sup></td><td>2</td><td>3</td><td>1.2 · 10<sup>-5</sup></td><td>8.7 · 10<sup>-6</sup></td></tr><tr><td>10<sup>-2</sup></td><td>2</td><td>4</td><td>1.8 · 10<sup>-7</sup></td><td>1.3 · 10<sup>-7</sup></td></tr><tr><td>10<sup>-3</sup></td><td>2</td><td>5</td><td>2.8 · 10<sup>-9</sup></td><td>2.2 · 10<sup>-9</sup></td></tr><tr><td>10<sup>-4</sup></td><td>2</td><td>6</td><td>4.4 · 10<sup>-11</sup></td><td>2.1 · 10<sup>-11</sup></td></tr><tr><td>10<sup>-5</sup></td><td>2</td><td>7</td><td>2.4 · 10<sup>-12</sup></td><td>1.0 · 10<sup>-11</sup></td></tr><tr><td>10<sup>-8</sup></td><td>2</td><td>11</td><td>0</td><td>0</td></tr><tr><td>10<sup>-12</sup></td><td>2</td><td>15</td><td>—</td><td>—</td></tr></tbody></table>

### 6.2.2. Crossing domains

In this second case, we consider the domain made of two crossing mono-dimensional domains, e.g., crossing fractures. We set $\Omega = \Omega_{horiz} \cup \Omega_{vert}$ where $\Omega_{horiz} = (0, 1) \times \{0.5\}$ and $\Omega_{vert} = \{0.5\} \times (0, 1)$; both $\Omega_{horiz}$ and $\Omega_{vert}$ are identifiable with $(0, 1)$. We consider a zero vector source term $f$ and a scalar source term given on $\Omega_{horiz}$ and $\Omega_{vert}$ respectively by

$$ q_{\text{horiz}}(\mathbf{x}) = \begin{cases} 1 & \text{if } x_1 \le 0.3 \\ -1 & \text{if } 0.3 < x_1 < 0.7 \\ 1 & \text{if } x_1 \ge 0.7 \end{cases} \quad \text{and} \quad q_{\text{vert}}(\mathbf{x}) = \begin{cases} 1 & \text{if } x_2 \le 0.3 \\ -1 & \text{if } 0.3 < x_2 < 0.7 \\ 1 & \text{if } x_2 \ge 0.7 \end{cases} $$

On the boundary we set $p(0, 0.5) = 0$ and $p(0.5, 0) = p(1, 0.5) = p(0.5, 1) = 0.1$.

**Linear case.** In this part, we consider the linear case with $\Lambda$ as in (6.2). Figure 6 shows the graphical representation of the solution for both the horizontal and vertical part of $\Omega$ for all the iterations of Algorithm 1. We notice the influence of the crossing through a velocity jump in $\Omega_{horiz}$, while the pressure profile is continuous as
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Solution for different iterations for the linear problem of Section 6.2.2. The pressure profile is amplified by a factor of 4. On the top for the horizontal part of $\Omega$, while on the bottom for the vertical one.</figcaption></figure>

condition (2.7) imposes. Also in this case, by changing the initial condition we obtain the same final outcome, where all the domain becomes $\Omega_{2,h}$.

Also with this more complex case, the obtained numerical evidence is insightful and shows good properties for the developed approximation framework which is apparently applicable to this case.

**Non-linear case.** In this part, we consider the non-linear case with $\Lambda$ as in (6.3). The obtained numerical solution is reported in Figure 7. The scheme takes 5 iterations to converge with increasing number of non-linear solver iterations as (1, 6, 8, 9). The obtained solution shows that the high-speed model, being Darcy-Forchheimer, is more proper to describe most of the problem leaving the slow Darcian regime in the vicinity of the boundary where the non-zero pressure condition is imposed. It is important to note that the plots show the norm of $u$ which presents a jump only in $\Omega_{\text{horiz}}$. Nevertheless, condition (2.7) is respected at the intersection since a
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Solution for different iterations for the non-linear problem of Section 6.2.2. The pressure profile is amplified by a factor of 8. On the top for the horizontal part of Ω, while on the bottom for the vertical one.</figcaption></figure>

velocity jump is also present in $\Omega_{\text{vert}}$. The representation of only $||u||$ hides this details. By changing the initial condition or refining the mesh, we obtain again the same final outcome.

We can conclude that also in presence of a non-linear and heterogeneous law the proposed framework works properly on this case.

### 6.2.3. Multiple fracture network

We consider now a complex fracture network, with geometry taken from Benchmark 1 of [14]. It is composed of 6 intersecting fractures as Figure 8 shows, along with the set of boundary conditions.

We denote by $\Omega_s$ the set of smaller fracture branches, which will be useful in the following. Null and unitary vector and scalar sources are considered, respectively. As done with the previous examples, we will consider the linear and non-linear case in the next Sections.
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Representation of the fracture network for the examples in Section 6.2.3. We have reported the boundary conditions as well as the portion of the network with small branches in red $\Omega_s$.</figcaption></figure>

**Linear case.** We consider here the same linear relations as in Section 6.2.1; see (6.2). The obtained solution is represented in Figure 9. The scheme converges after 4 iterations of Algorithm 1. We see an interesting result: $\Omega_{2,h}$, which is the high-velocity region, is automatically positioned on the main pathways between the inflow and outflow parts of the network. Again, the algorithm showed robustness when changing the initial configuration and refining the grid.

Even for this complex configuration, the proposed algorithm is capable to compute a reasonable solution with a limited cost.

**Non-linear case.** We consider in this part the non-linear case, where the constitutive law combination is given as

$$ \mathbf{\Lambda}(\mathbf{u}) = \begin{cases} \mathbf{u} & \text{in } \Omega_1, \\ (0.01 + 0.25\|\mathbf{u}\|)\mathbf{u} & \text{in } \Omega_2. \end{cases} $$

Algorithm 1 takes 4 steps to reach the final configuration with an average of 18 iterations of the non-linear solver for each step. Figure 10 shows the obtained numerical solution for different iterations. The obtained solution is again insightful, positioning the high-speed region, given by $\Omega_{2,h}$, in the longest fracture branches and the low-speed region $\Omega_{1,h}$ mainly in the fracture branches at the outflow. There is a transition zone from $\Omega_{2,h}$ to $\Omega_{1,h}$ that mainly takes place in $\Omega_s$.

A zoom-in is reported in Figure 11 to better clarify the evolution of $\Omega_{2,h}$ and $\Omega_{1,h}$ at the small fracture branches.

This final example shows a very interesting and physically sound final configuration, which might have been hard to predict without the framework introduced in this work. All these examples showed the applicability and importance of the model adaptation and support the presented Algorithm 1 to be a valid approach for its solution.

## 7. CONCLUSION

In this work we introduced a new model for a porous medium that is able to adapt the constitutive relation between velocity and pressure depending on the local magnitude of the fluid velocity, which is part of the unknowns.

We presented a mathematical formulation for it, and with an energy argument we were able to show that under reasonable hypotheses on the constitutive law the problem has a solution. When the transition zone
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Solution for different iterations, respectively (0, 1, 4), for the linear problem of Section 6.2.3. On the top the condition “if $\Omega_1$” is represented with red indicating “true” and blue “false”. On the bottom the norm of the velocity.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 10. Solution for different iterations, respectively (1, 3, 4), for the problem of Section 6.2.3, non-linear case. On the top the condition “if $\Omega_1$” is represented with red indicating “true” and blue “false”. On the bottom the norm of the velocity.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Zoom around the short fracture branches $\Omega_s$ of the “if $\Omega_1$” condition for different iterations, respectively (1, 3, 4), for the problem of Section 6.2.3.</figcaption></figure>

<table><caption>TABLE 2. Summary of results according to the transition zone inverse permeability jump $\lambda_{21} := \lambda_2 - \lambda_1$.</caption><thead><tr><th></th><th>λ<sub>21</sub> > 0</th><th>λ<sub>21</sub> = 0</th><th>λ<sub>21</sub> < 0</th></tr></thead><tbody><tr><td>Existence of solutions</td><td>Yes</td><td>Yes</td><td>Yes if d = 1</td></tr><tr><td>Convexity of energy</td><td>Yes</td><td>Yes</td><td>No*</td></tr><tr><td>Convergence of algorithm (d = 1)</td><td>No</td><td>Yes</td><td>Yes</td></tr></tbody></table>

Notes. (*)Unless in trivial case $d = 1$ and $\Sigma_v \neq \emptyset$.

inverse permeability jump is non-negative, the problem is convex and we could show existence in any space dimension; when it is negative, however, the problem becomes non-convex and we had to restrict our proof of existence to one space dimension. We also introduced a discrete algorithm that, for a given problem, tracks the low- and high-speed regions as well at the transition zone separating them. We considered various constitutive relations for distinct parts of the network, such as the classical Darcy law and the non-linear Darcy-Forchheimer law. Several numerical examples showed the validity of the proposed approach by increasing the geometrical and physical complexity of the problem. We noticed that when the transition zone inverse permeability jump is positive, the algorithm seems not to converge and oscillate indefinitely between two configurations. In the complementary non-positive case, the algorithm seems to behave and converge nicely. We summarize these results in Table 2.

From a modeling point of view, as mentioned in the introduction, future extensions will be the inclusion of transmission conditions with the rock matrix and the possibility of having more than two constitutive laws for the problem. Additionally, derivative-dependent laws such as Stokes' and Brinkman's will be considered.

From an analytical point of view, open questions include the existence of solutions when $d > 1$ and the transition zone inverse permeability jump is negative, and the extension of the existence results to tensor permeabilities when $d > 1$; see Remark 3.1. In addition, characterizing admissible constitutive laws on the transition zone (rather than leaving the choice as a free parameter of the model, as we did here) as well as determining the Hausdorff dimension of the transition zone are interesting questions that would help us study uniqueness of solutions; recall Remarks 2.3 and 3.6.

From a numerical point of view, further extensions will be the development of the tracking algorithm for $d > 1$ and the proof of its convergence when the transition zone inverse permeability jump is non-positive, and the development of an alternative algorithm when the jump is positive.
REFERENCES

[1] E. Ahmed, J. Jaffré and J.E. Roberts, A reduced fracture model for two-phase flow with different rock types. *Math. Comput. Simul.* **137** (2017) 49–70.

[2] C. Alboin, J. Jaffré, J.E. Roberts, X. Wang and C. Serres, Domain decomposition for some transmission problems in flow in porous media. In: Numerical Treatment of Multiphase Flows in Porous Media (Beijing, 1999). Vol. 552 of *Lecture Notes in Phys*. Springer, Berlin (2000) 22–34.

[3] L. Amir, M. Kern, V. Martin and J.E. Roberts, Décomposition de domaine et préconditionnement pour un modèle 3D en milieu poreux fracturé. In: *Proceeding of JANO 8, 8t Conference on Numerical Analysis and Optimization (2005).*

[4] P. Angot, Well-posed Stokes/Brinkman and Stokes/Darcy coupling revisited with new jump interface conditions. *ESAIM: M2AN* **52** (2018) 1875–1911.

[5] S.N. Antontsev and S.I. Shmarev, A model porous medium equation with variable exponent of nonlinearity: existence, uniqueness and localization properties of solutions. *Nonlinear Anal. Theory Methods App.* **60** (2005) 515–545.

[6] J. Audu, F. Fairag and S. Messaoudi, On the well-posedness of generalized Darcy–Forchheimer equation. *Boundary Value Prob.* **2018** (2018) 1–17.

[7] C. Bernardi, F. Hecht and O. Pironneau, Coupling Darcy and Stokes equations for porous media with cracks. *ESAIM: M2AN* **39** (2005) 7–35.

[8] S. Berrone, S. Pieraccini and S. Scialò, A PDE-constrained optimization formulation for discrete fracture network flows. *SIAM J. Sci. Comput.* **35** (2013) B487–B510.

[9] S. Berrone, S. Pieraccini and S. Scialò, On simulations of discrete fracture network flows with an optimization-based extended finite element method. *SIAM J. Sci. Comput.* **35** (2013) 908–935.

[10] D. Boffi, F. Brezzi and M. Fortin, Mixed Finite Element Methods and Applications. Springer Series in Computational Mathematics. Springer, Berlin, Heidelberg (2013).

[11] J.W. Both, J.M. Nordbotten and F.A. Radu, Free energy diminishing discretization of darcy-forchheimer flow in poroelastic media. In: Finite Volumes for Complex Applications IX – Methods, Theoretical Aspects, Examples, edited by R. Klöfkorn, E. Keilegavlen, F. A. Radu, and J. Fuhrmann. Springer International Publishing, Cham (2020) 203–211.

[12] M. Bulíček, J. Málek, and J. Žabenský, A generalization of the Darcy–Forchheimer equation involving an implicit, pressure-dependent relation between the drag force and the velocity. *J. Math. Anal. App.* **424** (2015) 785–801.

[13] C. Cancès, Finite volume scheme for two-phase flows in heterogeneous porous media involving capillary pressure discontinuities. *ESAIM: M2AN* **43** (2009) 973–1001.

[14] B. Flemisch, I. Berre, W. Boon, A. Fumagalli, N. Schwenck, A. Scotti, I. Stefansson and A. Tatomir, Benchmarks for single-phase flow in fractured porous media. *Adv. Water Res.* **111** (2018) 239–258.

[15] L. Formaggia, A. Fumagalli, A. Scotti and P. Ruffo, A reduced model for Darcy’s problem in networks of fractures. *ESAIM: M2AN* **48** (2014) 1089–1116.

[16] N. Frih, J.E. Roberts and A. Saada, Modeling fractures as interfaces: a model for Forchheimer fractures. *Comput. Geosci.* **12** (2008) 91–104.

[17] C. Geuzaine and J.-F. Remacle, Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities. *Int. J. Numer. Methods Eng.* **79** (2009) 1309–1331.

[18] R. Granero-Belinchón and S. Shkoller, Well-posedness and decay to equilibrium for the Muskat problem with discontinuous permeability. *Trans. Am. Math. Soc.* **372** (2019) 2255–2286.

[19] C. Johannes van Duijn, H. Eichel, R. Helmig and I.S. Pop, Effective equations for two-phase flow in porous media: the effect of trapping on the microscale. *Transp. Porous Media* **69** (2007) 411–428.

[20] E. Keilegavlen, R. Berge, A. Fumagalli, M. Starnoni, I. Stefansson, J. Varela and I. Berre, Porepy: an open-source software for simulation of multiphysics processes in fractured porous media. *Comput. Geosci.* **25** (2021) 243–265.

[21] P. Knabner and J.E. Roberts, Mathematical analysis of a discrete fracture model coupling Darcy flow in the matrix with Darcy–Forchheimer flow in the fracture. *ESAIM: M2AN* **48** (2014) 1451–1472.

[22] K. Kumar, F. List, I.S. Pop and F.A. Radu, Formal upscaling and numerical validation of unsaturated flow models in fractured porous media. *J. Comput. Phys.* **407** (2020) 109138.

[23] F. List and F.A. Radu, A study on iterative methods for solving Richards’ equation. *Comput. Geosci.* **20** (2016) 341–353.

[24] V. Martin, J. Jaffré and J.E. Roberts, Modeling fractures and barriers as interfaces for flow in porous media. *SIAM J. Sci. Comput.* **26** (2005) 1667–1691.

[25] F.A. Morales and R.E. Showalter, A Darcy–Brinkman model of fractures in porous media. *J. Math. Anal. App.* **452** (2017) 1332–1358.

[26] J.T. Oden and J.N. Reddy, Variational Methods in Theoretical Mechanics. Springer Berlin Heidelberg (1976).

[27] M.A. Peletier, Variational modelling: energies, gradient flows, and large deviations. Preprint arXiv:1402.1990 (2014).

[28] I.S. Pop, F. Radu and P. Knabner, Mixed finite elements for the Richards’ equation: linearization procedure. *J. Comput. Appl. Math.* **168** (2004) 365–373.

[29] F.A. Radu, K. Kumar, J.M. Nordbotten and I.S. Pop, A robust, mass conservative scheme for two-phase flow in porous media including Hölder continuous nonlinearities. *IMA J. Numer. Anal.* **38** (2017) 884–920.

[30] P.-A. Raviart and J.-M. Thomas, A mixed finite element method for second order elliptic problems. *Lecture Notes Math.* **606** (1977) 292–315.
[31] J.E. Roberts and J.-M. Thomas, Mixed and hybrid methods. In: Handbook of Numerical Analysis. Vol. II, *Handb. Numer. Anal., II*. North-Holland, Amsterdam (1991) 523–639.

[32] I. Rybak and S. Metzger, A dimensionally reduced stokes-Darcy model for fluid flow in fractured porous media. *App. Math. Comput.* **384** (2020) 125260.

[33] J.J. Salas, H. López and B. Molina, An analysis of a mixed finite element method for a Darcy–Forchheimer model. *Math. Comput. Modell.* **57** (2013) 2325–2338.

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

[35] J.L. Vazquez, The Porous Medium Equation. Oxford University Press (2006).

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

Subscribe to Open (S2O)
A fair and sustainable open access model

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

This journal is currently published in open access under a **Subscribe-to-Open model (S2O)**. S2O is a transformative model that aims to move subscription journals to open access. Open access is the free, immediate, online availability of research articles combined with the rights to use these articles fully in the digital environment. We are thankful to our subscribers and sponsors for making it possible to publish this journal in open access, free of charge for authors.

**Please help to maintain this journal in open access!**

Check that your library subscribes to the journal, or make a personal donation to the S2O programme, by contacting
*subscribers@edpsciences.org*

More information, including a list of sponsors and a financial transparency report, available at: https://www.edpsciences.org/en/maths-s2o-programme