MODEL ADAPTATION FOR NON-LINEAR ELLIPTIC EQUATIONS IN MIXED FORM: EXISTENCE OF SOLUTIONS AND NUMERICAL STRATEGIES

. 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.


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.
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 (multiphysics) 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 ( = 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.

Proposed model
We identify the porous medium with an open, bounded, connected set Ω ⊂ R with Lipschitz boundary Ω. We suppose Ω is filled with a fluid of constant density . A general constitutive relation coupled with the mass conservation in Ω reads where : Ω → R is the unknown velocity of the fluid, : Ω → R is the unknown pressure of the fluid, : Ω → R 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 is a given body force (e.g., gravity). Equation (2.1) is a law between the velocity field and the pressure field via some operator Λ. Let Σ v , Σ p ⊂ Ω be relatively open in Ω (i.e., Σ v and Σ p are each the intersection of an open subset of R with Ω) and such that Ω = Σ v ∪ Σ p and Σ v ∩ Σ p = ∅. To (2.2) we add the following boundary conditions: where is the outward normal unit vector of Ω. Here, 0 : Σ v → R and 0 : Σ p → R are given functions setting the conditions on the boundary on and , respectively. We are denoting maps on Ω and their traces on the boundary Ω by the same notation.

Velocity-pressure constitutive law
Examples of laws Λ relating velocity and pressure via (2.1) are where : Ω → R × is the permeability tensor and 0 and ∈ [2, ∞) some parameters. The notation ‖·‖ stands for the Euclidean norm on R . Choosing Λ S gives Stokes' equation, Λ D gives Darcy's equation, Λ B gives Brinkman's equation and Λ DF gives the generalized Darcy-Forchheimer equation (which simplifies into the classical Darcy-Forchheimer when = 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 Λ D may be preferred, whereas where it is high the Darcy-Forchheimer law Λ DF might be a better choice. For this reason, in this paper we consider the case where we need to choose from two laws Λ 1 and Λ 2 according to some threshold speed¯> 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 Being able to impose a law at the transition zone {‖ ‖ =¯} 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: where Λ is any law such as in (2.5).
Clearly, when Λ 1 ̸ = Λ 2 , this choice of adaptable law introduces a discontinuity at ‖ ‖ =¯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 Ω 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 Σ p is such that Vol −1 (Σ p ) = 0, where Vol −1 stands for the ( − 1)-dimensional Lebesgue measure, in order to ensure uniqueness of the pressure field satisfying Problem 2.1 when Λ is a classical continuous law, as opposed to (2.5), one imposes a constraint on the average of the pressure field: 1 |Ω| for some ∈ R, normally assumed null. We shall therefore impose (2.6) whenever Vol −1 (Σ p ) = 0, so that our problem in this case becomes: find functions and such that on Ω, under the constraint (2.6), where Λ 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 −1 (Σ 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 Λ 1 and Λ 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 = 1, Ω = (0, 1), Vol 0 (Σ p ) = 0, 0 ≡ 0 and ≡ 0, and set the threshold speed to¯= 1. Note that in this case, the velocity part of a strong solution ( , ) to Problem 2.1 is entirely determined independently of the chosen laws; indeed is the function defined by for all ∈ (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 is equal to 1 on a subset of (0, 1) with non-zero Lebesgue measure. Indeed, it is enough to pick

Multiple sub-domains
In the case Ω 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 Ω we introduce ⊂ Ω to be a sub-domain, with ∋ their total number. Clearly, given two distinct and , with ̸ = , we have˚∩˚= ∅ and also that Ω = ∪ =1 . We consider 2 sub-domains that meet at ℐ whose closure ℐ = ∩ =1 . To complete model in Problem 2.1 we impose the following classical conditions on ℐ: where with a superscript we indicate the corresponding object restricted to , and is a unit vector tangent to and pointing to ℐ, in the mono-dimensional case, and in addition normal to , in the multidimensional 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 Ω 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 ℐ, used for unsaturated and two-phase flow models, see for example [1,13,19,22].

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¯to be equal to 1. For a given Lebesgue measurable field : Ω → R , we write where Ω 1 ( ), Ω 2 ( ) and Γ( ) are what we have already respectively referred to as the low-speed region, highspeed region and transition zone (associated with ). Note that thanks to the measurability assumption on , these three subsets are also measurable. Obviously the family := {Ω 1 ( ), Ω 2 ( ), Γ( )} forms a partition of Ω, and we will refer to 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: where 1 (0) and 1 (0) stand respectively for the unit open ball and unit sphere in R centered at the origin 0 ). Note that if is not continuous, then Ω 1 ( ) and Ω 2 ( ) may not be open. For all ∈ [1, ∞), ∈ (0, ∞) and ⊂ R measurable we will denote by ( ) and , ( ) the Lebesgue space of measurable functions on with integrable th power and the th-order Sobolev space associated to ( ); we will also write ( ) in place of ( ( )) . As usual, in these spaces, equality between two functions is always intended in the almost-everywhere sense.

Assumptions on the velocity-pressure laws
In the following, the operator laws Λ 1 and Λ 2 are assumed to be of the form where is the characteristic function of any set ⊂ R . The functions 1 , 2 : [0, ∞) → [0, ∞) are continuous and increasing on [0, 1] and [1, ∞), respectively. Furthermore, 2 satisfies the following assumption: there exist 2 and , > 0 such that A recurrent notation we will use is 1 := 1 (1) and 2 := 2 (1), (3.4) and will call the difference 2 − 1 the transition zone inverse permeability jump. Symmetrically, whenever 1 := 1/ 1 and 2 := 1/ 2 are considered (cf. in particular Sects. 5 and 6), the difference 2 − 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 2 is increasing on [1, ∞), one must have 2 and 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 Λ D and Λ DF (cf. (2.4)), as desired in the first place. Indeed, one is allowed to consider that is, 1 ≡ 1 and 2 ≡ 2 , or that is, 1 ≡ 1 and 2 ( ) = 21 + 22 √ , where 1 , 2 , 21 and 22 are positive scalars. In the former case we would require = 2, whereas in the latter = 3.
Remark 3.1. This setting where Λ 1 and Λ 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: and analogously for Λ 2 , where 1 ∈ R × is a symmetric positive definite matrix encoding a tensor permeability. Even more general forms are envisageable: where the 1, and 1, are different law functions and permeability tensors and the dot · stands for the Euclidean inner product in R . 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 = −1 the conjugate exponent of . Thanks to the continuity of 1 and the right-hand inequality in (3.3), we see that the operators Λ 1 and Λ 2 map (Ω) into (Ω). Because we allow for any law on the transition zone, this invites us to consider the multi-valued setting where our combined law Λ maps (Ω) into the power set 2 (Ω) and is given by We highlight the fact that the transition zone law ℎ in the definition above really is a function in (Ω) and not a trace on Γ( ). In this way, ℎ Γ( ) is simply the restriction of an (Ω) function to Γ( ). This is also coherent with our energetic approach below, where we will link to Λ the subdifferential of a functional defined on (Ω), whose topological dual is indeed (Ω).
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 Λ in (3.5) is replaced by where ∈ R × is symmetric positive definite and : [0, ∞) → [0, ∞) is a continuous, increasing function such that + 2 satisfies (3.3) in place of 2 .
We set 0 = in the case Vol −1 (Σ p ) = 0 and 0 = + ∇( 0 ) in the case Vol −1 (Σ p ) > 0, with : 1 , (Σ p ) → 1, (Ω) any extension operator being right-inverse of the 1, (Ω)-trace operator, and by defining the Sobolev space We endow 1, 0 (Ω) with the norm ‖ ‖ 1, 0 (Ω) = ‖∇ ‖ (Ω) for all ∈ 1, 0 (Ω). 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: We emphasize here that the well-posedness of Problem 3.4 is not affected by the choice of the extension in the definition of 0 .
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 (Ω) function as transition zone law (cf. ℎ 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.

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, 2 − 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 1 2 from the non-convex case 1 > 2 (cf. Thms. 4.10 and 4.11). We are only able to treat the one-dimensional setting = 1 when 1 > 2 .
Proof. By Lemma 2.1 of [6], it suffices to show that there exists > 0 such that .

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 ℱ : → R. For all ∈ we define the (Fréchet) subdifferential ℱ( ) of ℱ at by: We say that ∈ is a (strong) local minimizer of ℱ if there exists > 0 such that for all ∈ we have ℱ( + ) ℱ( ) for all ∈ [0, ).
Remark 4.6. The important property of the subdifferential to keep in mind here is that if ∈ is a local minimizer of a functional ℱ : → R, then 0 (Ω) ∈ ℱ( ). When ℱ is convex, the reverse of this statement is also true: if 0 (Ω) ∈ ℱ( ), then is a local minimizer of ℱ.
We write Ψ : [0, ∞) → R the function given by where Φ 1 , Φ 2 : [0, ∞) → R are primitives of 1 2 and 2 2 (cf. (3.2)) such that Φ 1 (1) = Φ 2 (1) = 0. We define the dissipation : → R by Thanks to our assumptions on 1 and 2 we easily get that the domain of is indeed all of , and it is thus a well-defined functional from into R. Consider the following problem: The following result holds: Write 1 the subset of consisting of the functions which are supported in Ω 1 ( +ˆ). Then, for all ∈ 1 , since − ∈ 1 as well, we get which, by Lebesgue's dominated convergence theorem and the differentiability of Φ 1 , yield All in all we get ⟨ , ⟩ = ⟨Λ 1 ( +ˆ), ⟩ for all ∈ 1 .
Similarly, denoting by 2 the subset of consisting of the functions which are supported in Ω 2 ( +ˆ), we get Thus the function defined by is such that ⟨Λ , ⟩ = ⟨ , ⟩ for all ∈ . Therefore, from Problem 4.7 we obtain Moreover Λ ∈ Λ( +ˆ), where we recall Λ is in (3.5). Hence is solution to Problem 4.3.
We now define the energy ℰ : → R associated to by Let us write ℰ ⊂ the, possibly empty, set of local minimizers of ℰ. Consider the following minimization problem.
By Remark 4.6, any solution to Problem 4.9 is also solution to Problem 4.7; if ℰ 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.
We now want to use the direct method of the calculus of variations to show that ℰ has in fact a unique global (and thus local) minimizer. Let ( ) ∈N ⊂ be a minimizing sequence for ℰ. Then we know there exists ∈ N large enough and > 0 such that ℰ( ) < for all > . Without loss of generality we can therefore assume that the sequence (ℰ( )) ∈N is bounded by some constant > 0. Hence, thanks to the left-hand inequality in (3.3), for all ∈ N we have which shows that the sequence (‖ ‖ (Ω) ) ∈N is bounded. Thus we can extract a subsequence from ( ) ∈N , still denoted ( ) ∈N which converges weakly to some ∈ (Ω). Since further is weakly closed we in fact have ∈ . Because Ψ ∘ ‖·‖ 2 is convex, the dissipation is weakly lower semi-continuous and so is the energy ℰ. We therefore yield inf so that ℰ( ) = inf ∈ ℰ( ) and is a global minimizer of ℰ and so ∈ ℰ . To show that ℰ is a singleton it is enough to prove that ℰ is strictly convex. Let , ∈ and ∈ (0, 1). Suppose furthermore that ̸ = and write ⊂ Ω the set where and are different; the Lebesgue measure of is therefore positive. Note that we must have Ψ( ) ̸ = 0 for all ̸ = 1. Using the strict convexity of Ψ ∘ ‖·‖ 2 we therefore get so that the dissipation is strictly convex. Consequently, the energy ℰ is also strictly convex and the proof is over.

Case 1 > 2
This case is more difficult to tackle than the case 1 2 . Indeed, we lose the convexity of the integrand Ψ ∘ ‖·‖ 2 (cf. proof of Thm. 4.10), which by Tonelli's theorem of functional analysis means that ℰ 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 ℰ.
To simplify our task at this point, we shall restrict to the one-dimensional case (i.e., = 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 depending on the boundary conditions, as will be clear from the proof of the following theorem: Theorem 4.11 (Existence and uniqueness when 1 > 2 ). Let = 1 and suppose that 1 > 2 . If Vol 0 (Σ v ) = 0, then Problem 4.9 has a solution. If instead Vol 0 (Σ v ) > 0, then Problem 4.9 has a unique solution.
Case Vol 0 (Σ v ) = 0. First note that here Σ p = {0, 1}. Let us characterize in this case. To this end, note that any ∈ with ∫︀ 1 0 = 0 has a primitive in 1, 0 ((0, 1)); indeed, the function : (0, 1) → R defined by and is constant (almost everywhere). Since constant functions clearly belong to , this shows that is in fact the set of all constant functions on (0, 1) and we can identify with R.
This in particular means that the energy ℰ defined in (4.4) can be identified with the following function : R → R: where we recall that Ψ is given in (4.2) andˆ:=ˆis as in Lemma 4.2. We can use the direct method of the calculus of variations on in the Euclidean topology in R. Let ( ) ∈N ⊂ R be a minimizing sequence for . Using a similar calculation as in the proof of Theorem 4.10, thanks to (3.3) we can show that ( ) ∈N is bounded in R. Therefore, there exists a subsequence of ( ) ∈N , still denoted by ( ) ∈N , converging to some ∈ R. By continuity of Ψ and Fatou's lemma we get that is lower semi-continuous on R, so that showing that is a global minimizer of and the constant function belongs to ℰ .
Similarly to the previous case, let us characterize . Note that any ∈ has a primitive in 1, 0 ((0, 1)); indeed, the function : (0, 1) → R defined for all ∈ (0, 1) by so that = 0. This shows that = {0}, i.e., contains only the zero function on (0, 1). Trivially then, the zero function is the unique global minimizer of the energy ℰ in (4.4) and ℰ is a singleton. Note that in this case the energy is only trivially convex.

Numerical approximation
In this part we present the numerical approximation for the considered problem when = 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 Ω ℎ composed of non-overlapping segments ∈ Ω ℎ that approximate Ω; we clearly have Ω ℎ = ∪ ∈Ω ℎ . At the discrete level we can define the approximate configuration ℎ which is given by the set ℎ = {Ω 1,ℎ , Ω 2,ℎ , Γ ℎ }, where respectively Ω 1,ℎ ⊂ Ω ℎ and Ω 2,ℎ ⊂ Ω ℎ are the approximations of the domains Ω 1 and Ω 2 . Γ ℎ is the approximation of the transition zone Γ. Note that for simplicity we are dropping the dependence of these regions on the velocity field. For each element we name 1 and 2 its two extremal vertices and ℎ its length. We let ℎ = max ∈Ω ℎ ℎ 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 ( ℎ , ℎ ) the approximation of ( , ) for a given mesh.

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 a tentative configuration ( ) ℎ approximates ℎ . For a given configuration ( −1) ℎ , the proposed algorithm solves the differential problem obtaining a new velocity field ( ) . To speed up the computation, we evaluate condition (2.5) only at the extremities of each grid element . If we obtain opposite values, we can thus determine the position of a transition zone in the considered element up to a given tolerance Γ . This part can be coded with an embarrassingly parallel workload, speeding up the algorithm. Having checked all the elements and computed the new configuration ( ) ℎ , the algorithm resets with the new configuration as starting point. The exit strategy considers the position of Γ ℎ for two successive iterations: if their distance is smaller than a given threshold Ω , 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 ; this is a direct implication of the fact that Ω ( ) 1,ℎ or Ω ( ) 2,ℎ 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 Ω is comparable with the mesh size given at the outset of the problem. Also, the value of Γ 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 Γ.

Discretization in space
Following the algorithm discussed before, we assume that a configuration ( ) ℎ is given. In the case of nonlinear constitutive relations Λ, 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 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 Λ.
Following [15], to solve the problem we consider the mixed finite element approximation of lowest-order degree. For a given mesh Ω ℎ we have ( ℎ , ℎ ) ∈ (P 0 (Ω ℎ ), RT 0 (Ω ℎ )), 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 ℎ goes to zero.

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.

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 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 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 = 50[m] while for the former we set = 10 −4 [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 ≈ 48.2[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.

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 . The threshold on the velocity norm is set as = 0.15, so that we have If not otherwise specified, we consider a mesh size equal to ℎ = 5 · 10 −2 as well as Ω = ℎ; 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 Γ = 10 −10 . In all the cases, the initial configuration is chosen to be Ω (With the notation used previously, = 1/ , ∈ {1, 2}.) Specifically, we have 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 Ω 1 and a Darcy-Forchheimer flow in Ω 2 . Specifically, we have If not specified, we consider 50 as maximum for the number of iterations of the non-linear solver with tolerance nl equal to 10 −4 .

Single domain
In this first case, we consider Ω = (0, 1). Boundary conditions are set to zero for the pressure. In the sequel, we consider a linear and non-linear law relationship between and . The scalar and vector source terms are set equal to Linear case. In this part, we consider the linear case with Λ 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 Ω 1 " legend represents is a binary outcome saying if the region is Ω 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 Γ, which might change during the determination of the final configuration. By changing the  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 2 while fixing 1 = 1. By increasing the maximum number of iterations to 1000, Figure 3 on the left shows the impact of changing 2 on the number of iterations for the Algorithm 1. We see that for 2 1 the solution is always computable while we cannot draw the same conclusion for 2 < 1 .
An example for 2 = 0.25 is shown in Figure 4, where we notice that "if Ω 1 " and "condition" are perfectly flipped between any two successive iterations. The algorithm "jumps" between two states and thus does not converge. By setting = 0 and (1) = 0.2, the latter to avoid that the algorithm converges in 1 iteration for each value of 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 2 > 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 2 1 .
Non-linear case. In this part, we consider the non-linear case with Λ 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 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 ( ref , ref ) with tolerance 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 and err u for computed as 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.
Linear case. In this part, we consider the linear case with Λ as in (6.2). Figure 6 shows the graphical representation of the solution for both the horizontal and vertical part of Ω for all the iterations of Algorithm 1. We notice the influence of the crossing through a velocity jump in Ω horiz , while the pressure profile is continuous as condition (2.7) imposes. Also in this case, by changing the initial condition we obtain the same final outcome, where all the domain becomes Ω 2,ℎ .
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 Λ 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 which presents a jump only in Ω horiz . Nevertheless, condition (2.7) is respected at the intersection since a velocity jump is also present in Ω vert . The representation of only ‖ ‖ 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.

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 Ω 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. 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: Ω 2,ℎ , 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 in Ω 1 , (0.01 + 0.25‖ ‖) in Ω 2 .
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 Ω 2,ℎ , in the longest fracture branches and the low-speed region Ω 1,ℎ mainly in the fracture branches at the outflow. There is a transition zone from Ω 2,ℎ to Ω 1,ℎ that mainly takes place in Ω . A zoom-in is reported in Figure 11 to better clarify the evolution of Ω 2,ℎ and Ω 1,ℎ 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.

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 9. Solution for different iterations, respectively (0, 1, 4), for the linear problem of Section 6.2.3. On the top the condition "if Ω 1 " is represented with red indicating "true" and blue "false". On the bottom the norm of the velocity. 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 Ω 1 " is represented with red indicating "true" and blue "false". On the bottom the norm of the velocity. Figure 11. Zoom around the short fracture branches Ω of the "if Ω 1 " condition for different iterations, respectively (1,3,4), for the problem of Section 6.2.3. 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 > 1 and the transition zone inverse permeability jump is negative, and the extension of the existence results to tensor permeabilities when > 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 the development of the tracking algorithm for > 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.