domain: chim sde math

text:
Alessio Fumagalli
email: alessio.fumagalli1984@gmail.com

Francesco Saverio Patacchini

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

Keywords: Mathematics Subject Classification. 76S05, 35A01, 35A15, 65N50, 65N55 Adaptive constitutive law, elliptic equation, fractured porous media, non-linearity, variational formulation

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 [START_REF] Morales | A Darcy-Brinkman model of fractures in porous media[END_REF][START_REF] Rybak | A dimensionally reduced stokes-Darcy model for fluid flow in fractured porous media[END_REF]. 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 [START_REF] Frih | Modeling fractures as interfaces: a model for Forchheimer fractures[END_REF][START_REF] Knabner | Mathematical analysis of a discrete fracture model coupling Darcy flow in the matrix with Darcy-Forchheimer flow in the fracture[END_REF] 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 [START_REF] Zeng | A criterion for non-Darcy flow in porous media[END_REF]), 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 [START_REF] Oden | Variational Methods in Theoretical Mechanics[END_REF][START_REF] Spena | A potential formulation of non-linear models of flow through anisotropic porous media[END_REF]. For related time-dependent problems involving the minimization of a free energy giving rise to a gradient flow structure, we refer the reader to [START_REF] Both | Free energy diminishing discretization of darcy-forchheimer flow in poroelastic media[END_REF][START_REF] Peletier | Variational modelling: energies, gradient flows, and large deviations[END_REF]. Moreover, for a thorough overview of analytical aspects of porous medium equations, with and without time dependence, see [START_REF] Vazquez | The Porous Medium Equation[END_REF], and for an example of existence and uniqueness theory of non-linear Darcy-like equations, see [START_REF] Antontsev | A model porous medium equation with variable exponent of nonlinearity: existence, uniqueness and localization properties of solutions[END_REF]. 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 [START_REF] Keilegavlen | Porepy: an open-source software for simulation of multiphysics processes in fractured porous media[END_REF], 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

Λ(𝑢) = -∇𝑝 + 𝑓 , (2.1) div𝑢 = 𝑞, (2.2)
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:

{ 𝑢 • 𝑛 = 𝑢 0 on Σ v , 𝑝 = 𝑝 0 on Σ p , (2.3)
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

Λ S (𝑢) = -Δ𝑢, Λ D (𝑢) = 𝐾 -1 𝑢, Λ B (𝑢) = -Δ𝑢 + 𝐾 -1 𝑢, Λ DF (𝑢) = 𝐾 -1 𝑢 + 𝛼‖𝑢‖ 𝑟-2 𝑢, (2.4)
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

Λ(𝑢) = { Λ 1 (𝑢) wherever ‖𝑢‖ < ū, Λ 2 (𝑢) wherever ‖𝑢‖ > ū. (2.5)
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:

Problem 2.1 (Strong form). Find 𝑢 : Ω → R 𝑑 and 𝑝 : Ω → R such that ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ div𝑢 = 𝑞 on Ω, Λ(𝑢) = -∇𝑝 + 𝑓 on Ω, 𝑢 • 𝑛 = 𝑢 0 on Σ v , 𝑝 = 𝑝 0 on Σ p ,
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 [START_REF] Granero-Belinchón | Well-posedness and decay to equilibrium for the Muskat problem with discontinuous permeability[END_REF], by Darcy-Stokes interface coupling problems in [START_REF] Angot | Well-posed Stokes/Brinkman and Stokes/Darcy coupling revisited with new jump interface conditions[END_REF][START_REF] Bernardi | Coupling Darcy and Stokes equations for porous media with cracks[END_REF] or fracture barrier problems in [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF]. 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 [START_REF] Bulíček | A generalization of the Darcy-Forchheimer equation involving an implicit, pressuredependent relation between the drag force and the velocity[END_REF] 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

|Ω| ∫ Ω 𝑝 = 𝜛, (2.6)
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

⎧ ⎪ ⎨ ⎪ ⎩ div𝑢 = 𝑞 on Ω, Λ(𝑢) = -∇𝑝 + 𝑓 on Ω, 𝑢 • 𝑛 = 𝑢 0 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 𝑢(𝑥) = ∫ 𝑥 0 𝑞(𝑦) d𝑦 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

𝑞(𝑥) = ⎧ ⎪ ⎨ ⎪ ⎩ 4
for all 𝑥 ∈ (0, 0.25), 0 for all 𝑥 ∈ [0.25, 0.75), -4 for all 𝑥 ∈ [0.75, 1), 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 𝑞 such as to satisfy the compatibility condition

∫ Ω 𝑞 = ∫ 𝜕Ω 𝑢 0 = 0.)

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 ℐ:

𝑛 ∑ 𝑖=1 𝑢 𝑖 • 𝑛 𝑖 = 0 and 𝑝 𝑖 = 𝑝 𝑗 , ∀ 𝑖, 𝑗 = 1, . . . , 𝑛, (2.7)
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 [START_REF] Alboin | Domain decomposition for some transmission problems in flow in porous media[END_REF][START_REF] Amir | Décomposition de domaine et préconditionnement pour un modèle 3D en milieu poreux fracturé[END_REF][START_REF] Berrone | A PDE-constrained optimization formulation for discrete fracture network flows[END_REF][START_REF] Berrone | On simulations of discrete fracture network flows with an optimization-based extended finite element method[END_REF]. 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 [START_REF] Ahmed | A reduced fracture model for two-phase flow with different rock types[END_REF][START_REF] Cancès | Finite volume scheme for two-phase flows in heterogeneous porous media involving capillary pressure discontinuities[END_REF][START_REF] Johannes Van Duijn | Effective equations for two-phase flow in porous media: the effect of trapping on the microscale[END_REF][START_REF] Kumar | Formal upscaling and numerical validation of unsaturated flow models in fractured porous media[END_REF].

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

Ω 1 (𝑢) = {𝑥 ∈ Ω | ‖𝑢(𝑥)‖ < 1}, Ω 2 (𝑢) = {𝑥 ∈ Ω | ‖𝑢(𝑥)‖ > 1}, Γ(𝑢) = Ω ∖ (Ω 1 (𝑢) ∪ Ω 2 (𝑢)) = {𝑥 ∈ Ω | ‖𝑢(𝑥)‖ = 1}, ( 3.1)
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:

Ω 1 (𝑢) = 𝑢 -1 (𝐵 1 (0)), Ω 2 (𝑢) = 𝑢 -1 (𝐸 1 (0)), Γ(𝑢) = 𝑢 -1 (𝑆 1 (0)),
where 𝐵 1 (0) and 𝑆 1 (0) stand respectively for the unit open ball and unit sphere in R 𝑑 centered at the origin 0 and 𝐸 1 (0) = R 𝑑 ∖ (𝐵 1 (0) ∪ 𝑆 1 (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

Λ 1 (𝑢) = 𝜑 1 ( ‖𝑢‖ 2 ) 𝑢 𝜒 Ω1(𝑢) and Λ 2 (𝑢) = 𝜑 2 ( ‖𝑢‖ 2 ) 𝑢 𝜒 Ω2(𝑢) , (3.2)
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

𝑐𝑎 𝑟-2 2 𝜑 2 (𝑎) 𝐶 ( 1 + 𝑎 𝑟- 2 2
) for all 𝑎 1. (

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

Λ 1 (𝑢) = 𝜆 1 𝑢 𝜒 Ω1(𝑢) and Λ 2 (𝑢) = 𝜆 2 𝑢 𝜒 Ω2(𝑢) , that is, 𝜑 1 ≡ 𝜆 1 and 𝜑 2 ≡ 𝜆 2 , or Λ 1 (𝑢) = 𝜆 1 𝑢 𝜒 Ω1(𝑢) and Λ 2 (𝑢) = (𝜆 21 + 𝜆 22 ‖𝑢‖)𝑢 𝜒 Ω2(𝑢) . 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 [START_REF] Spena | A potential formulation of non-linear models of flow through anisotropic porous media[END_REF], are for instance given by the following:

Λ 1 (𝑢) = 𝜑 1 (𝜆 1 𝑢 • 𝑢)𝜆 1 𝑢 𝜒 Ω1(𝑢)
and analogously for Λ 2 , where 𝜆 1 ∈ R 𝑑×𝑑 is a symmetric positive definite matrix encoding a tensor permeability. Even more general forms are envisageable:

Λ 1 (𝑢) = 𝑛 ∑ 𝑗=1 𝜑 1,𝑗 (𝜆 1,𝑗 𝑢 • 𝑢)𝜆 1,𝑗 𝑢 𝜒 Ω1(𝑢) ,
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

Λ(𝑢) = { Λ 1 (𝑢) + Λ 2 (𝑢) + ℎ𝜒 Γ(𝑢) } ℎ∈𝐿 𝑠 (Ω) , (3.5)
that is, Λ(𝑢) is the set of all the functions of the form Λ 1 (𝑢) + Λ 2 (𝑢) + ℎ𝜒 Γ(𝑢) , where ℎ ∈ 𝐿 𝑠 (Ω). 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

Λ(𝑢) = { 𝛽(𝜆𝑢 • 𝑢)𝜆𝑢 + Λ 1 (𝑢) + Λ 2 (𝑢) + ℎ𝜒 Γ(𝑢) } ℎ∈𝐿 𝑠 (Ω) ,
where 𝜆 ∈ R 𝑑×𝑑 is symmetric positive definite and 𝛽 : [0, ∞) → [0, ∞) is a continuous, increasing function such that 𝛽 + 𝜑 2 satisfies (3.3) in place of 𝜑 2 .

Weak formulation

We fix

𝑞 ∈ 𝐿 𝑟 (Ω), 𝑢 0 ∈ 𝐿 𝑟 (Σ v ), 𝑝 0 ∈ 𝑊 1 𝑟 ,𝑠 (Σ p ) and 𝑓 ∈ 𝐿 𝑠 (Ω). Consider the space  𝑊 1,𝑠 (Ω) = ⎧ ⎪ ⎨ ⎪ ⎩ { 𝜉 ∈ 𝑊 1,𝑠 (Ω) | 1 |Ω| ∫ Ω 𝜉 = 𝜛 } if Vol 𝑑-1 (Σ p ) = 0, { 𝜉 ∈ 𝑊 1,𝑠 (Ω) | 𝜉 = 𝑝 0 on Σ p } if Vol 𝑑-1 (Σ p ) > 0,
where we recall that Vol 𝑑-1 is the (𝑑 -1)-dimensional Lebesgue measure. We give a first weak formulation of Problem 2.1:

Problem 3.3 (Weak form I). Find (𝑢, 𝑝) ∈ 𝐿 𝑟 (Ω) × ̃ 𝑊 1,𝑠 (Ω) so that there exists Λ 𝑢 ∈ Λ(𝑢) such that ∫ Ω Λ 𝑢 • 𝜙 = - ∫ Ω ∇𝑝 • 𝜙 + ∫ Ω 𝑓 • 𝜙 ∀𝜙 ∈ 𝐿 𝑟 (Ω), ∫ Ω ∇𝜓 • 𝑢 = - ∫ Ω 𝑞𝜓 + ∫ Σv 𝑢 0 𝜓 ∀𝜓 ∈ ̃ 𝑊 1,𝑠 (Ω),
where we recall that Λ is given in (3.5).

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

𝑊 1,𝑠 0 (Ω) = ⎧ ⎨ ⎩ { 𝜉 ∈ 𝑊 1,𝑠 (Ω) | ∫ Ω 𝜉 = 0 } if Vol 𝑑-1 (Σ p ) = 0, { 𝜉 ∈ 𝑊 1,𝑠 (Ω) | 𝜉 = 0 on Σ p } if Vol 𝑑-1 (Σ p ) > 0,
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:

Problem 3.4 (Weak form II). Find (𝑢, 𝑝) ∈ 𝐿 𝑟 (Ω) × 𝑊 1,𝑠 0 (Ω) so that there exists Λ 𝑢 ∈ Λ(𝑢) such that ∫ Ω Λ 𝑢 • 𝜙 = - ∫ Ω ∇𝑝 • 𝜙 + ∫ Ω 𝑓 0 • 𝜙 ∀𝜙 ∈ 𝐿 𝑟 (Ω), ∫ Ω ∇𝜓 • 𝑢 = - ∫ Ω 𝑞𝜓 + ∫ Σv 𝑢 0 𝜓 ∀𝜓 ∈ 𝑊 1,𝑠 0 (Ω).
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 .

Thanks to Riesz's representation theorem, we shall equivalently manipulate functions in 𝐿 𝑠 (Ω) as elements of (𝐿 𝑟 (Ω)) * , the topological dual of 𝐿 𝑟 (Ω). In particular, this means that, given 𝑢 ∈ 𝐿 𝑟 (Ω), operators Λ 𝑢 in Λ(𝑢) will also be seen as maps from 𝐿 𝑟 (Ω) to (𝐿 𝑟 (Ω)) * ∼ 𝐿 𝑠 (Ω). Denoting by ⟨•, •⟩ the dual mapping on 𝐿 𝑠 (Ω) × 𝐿 𝑟 (Ω), the weak formulation in Problem 3.4 can be equivalently rewritten as follows: Problem 3.5 (Weak form III). Find (𝑢, 𝑝) ∈ 𝐿 𝑟 (Ω) × 𝑊 1,𝑠 0 (Ω) so that there exists Λ 𝑢 ∈ Λ(𝑢) such that

⟨Λ 𝑢 , 𝜙⟩ = -⟨∇𝑝, 𝜙⟩ + ⟨𝑓 0 , 𝜙⟩ ∀𝜙 ∈ 𝐿 𝑟 (Ω), ⟨∇𝜓, 𝑢⟩ = - ∫ Ω 𝑞𝜓 + ∫ Σv 𝑢 0 𝜓 ∀𝜓 ∈ 𝑊 1,𝑠 0 (Ω).
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 𝐿 𝑠 (Ω) 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,[START_REF] Salas | An analysis of a mixed finite element method for a Darcy-Forchheimer model[END_REF]);

(2) we derive an energetic formulation whose minimizers are solutions to this reduced problem on the velocity field (cf. Problem 4.9 and [START_REF] Spena | A potential formulation of non-linear models of flow through anisotropic porous media[END_REF]); (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 .

Reduction to a problem on the velocity field

Define the maps 𝐵 : 𝐿 𝑟 (Ω) → (𝑊 1,𝑠 0 (Ω)) * and 𝐵 * :

𝑊 1,𝑠 0 (Ω) → (𝐿 𝑟 (Ω)) * by 𝐵(𝜙)(𝜓) = 𝐵 * (𝜓)(𝜙) = ⟨∇𝜓, 𝜙⟩ for all 𝜙 ∈ 𝐿 𝑟 (Ω) and 𝜓 ∈ 𝑊 1,𝑠 0 (Ω), (4.1)
Let us introduce the set

𝑉 = { 𝜙 ∈ 𝐿 𝑟 (Ω) | ∀ 𝜓 ∈ 𝑊 1,𝑠 0 (Ω), ⟨∇𝜓, 𝜙⟩ = 0 } ,
which satisfies 𝑉 = Ker(𝐵). We write 𝑉 ⊥ ⊂ 𝐿 𝑠 (Ω) the polar subspace of 𝑉 , that is,

𝑉 ⊥ = {𝑔 ∈ 𝐿 𝑠 (Ω) | ∀ 𝜙 ∈ 𝑉, ⟨𝑔, 𝜙⟩ = 0}
, and 𝐿 𝑟 (Ω)/𝑉 ⊂ 𝐿 𝑟 (Ω) the quotient space of 𝐿 𝑟 (Ω) by 𝑉 .

The following three lemmas are inspired from their equivalents in [START_REF] Audu | On the well-posedness of generalized Darcy-Forchheimer equation[END_REF].

Lemma 4.1. The maps 𝐵 and 𝐵 * in (4.1) are isomorphisms from 𝐿 𝑟 (Ω)/𝑉 to (𝑊 1,𝑠 0 (Ω)) * and from 𝑊 1,𝑠 0 (Ω) to 𝑉 ⊥ , respectively. Proof. By Lemma 2.1 of [START_REF] Audu | On the well-posedness of generalized Darcy-Forchheimer equation[END_REF], it suffices to show that there exists 𝛾 > 0 such that inf

𝜓∈𝑊 1,𝑠 0 (Ω) 𝜓̸ =0 sup 𝜙∈𝐿 𝑟 (Ω) 𝜙̸ =0 ⟨∇𝜓, 𝜙⟩ ‖𝜓‖ 𝑊 1,𝑠 0 (Ω) ‖𝜙‖ 𝐿 𝑟 (Ω)
𝛾.

To this end, let 𝜓 ∈ 𝑊 1,𝑠 0 (Ω) with 𝜓 ̸ = 0. Then, by identifying 𝐿 𝑠 (Ω) with (𝐿 𝑟 (Ω)) * we get

‖𝜓‖ 𝑊 1,𝑠 0 (Ω) = ‖∇𝜓‖ 𝐿 𝑠 (Ω) = sup 𝜙∈𝐿 𝑟 (Ω) 𝜙̸ =0 ⟨∇𝜓, 𝜙⟩ ‖𝜙‖ 𝐿 𝑟 (Ω) , so that 1 = ‖𝜓‖ 𝑊 1,𝑠 0 (Ω) ‖𝜓‖ 𝑊 1,𝑠 0 (Ω) = sup 𝜙∈𝐿 𝑟 (Ω) 𝜙̸ =0 ⟨∇𝜓, 𝜙⟩ ‖𝜓‖ 𝑊 1,𝑠 0 (Ω) ‖𝜙‖ 𝐿 𝑟 (Ω)
• Taking above the infimum over all 𝜓 ∈ 𝑊 1,𝑠 0 (Ω) with 𝜓 ̸ = 0 ends the proof. Now, via the following two lemmas, we simplify our weak formulation in Problem 3.5 into a problem restricted to 𝑉 (cf. Problem 4.3 below). Lemma 4.2. There exists a unique û ∈ 𝐿 𝑟 (Ω)/𝑉 such that

⟨∇𝜓, û⟩ = - ∫ Ω 𝑞𝜓 + ∫ Σv 𝑢 0 𝜓 for all 𝜓 ∈ 𝑊 1,𝑠 0 (Ω).
Proof. For all 𝜓 ∈ 𝑊 1,𝑠 0 (Ω), let

𝐹 (𝜓) = - ∫ Ω 𝑞𝜓 + ∫ Σv 𝑢 0 𝜓,
and bound 𝐹 (𝜓) as follows:

|𝐹 (𝜓)| ‖𝑞‖ 𝐿 𝑟 (Ω) ‖𝜓‖ 𝐿 𝑠 (Ω) + ‖𝑢 0 ‖ 𝐿 𝑟 (Σv) ‖𝜓‖ 𝐿 𝑠 (𝜕Ω) ( 𝐶‖𝑞‖ 𝐿 𝑟 (Ω) + ‖𝑢 0 ‖ 𝐿 𝑟 (Σv) ) ‖𝜓‖ 𝑊 1,𝑠 0 (𝜕Ω)
, where we used boundedness of the trace operator on 𝜓 and where 𝐶 > 0 is a constant stemming from Poincar's inequality. Recalling that indeed 𝑞 and 𝑢 0 are assumed to belong to 𝐿 𝑟 (Ω) and 𝐿 𝑟 (Σ v ), respectively, the map 𝐹 is then a bounded linear map from 𝑊 1,𝑠 0 (Ω) into R, i.e., 𝐹 ∈ (𝑊 1,𝑠 0 (Ω)) * . By Lemma 4.1, we thus know there exists a unique û ∈ 𝐿 𝑟 (Ω)/𝑉 such that

⟨∇𝜓, û⟩ = 𝐵(û)(𝜓) = 𝐹 (𝜓) = - ∫ Ω 𝑞𝜓 + ∫ Σv 𝑢 0 𝜓 for all 𝜓 ∈ 𝑊 1,𝑠 0 (Ω),
which is the desired result. Furthermore, for all 𝜓 ∈ 𝑊 1,𝑠 0 (Ω) we find ⟨∇𝜓, 𝑣⟩ = ⟨∇𝜓, 𝑢⟩ -⟨∇𝜓, û⟩ = 0, so that 𝑣 ∈ 𝑉 and 𝑣 satisfies Problem 4. [START_REF] Amir | Décomposition de domaine et préconditionnement pour un modèle 3D en milieu poreux fracturé[END_REF].

Suppose now that 𝑣 ∈ 𝑉 satisfies Problem 4.3. Then, Λ 𝑣 -𝑓 0 ∈ 𝑉 ⊥ . By Lemma 4.1 we know 𝐵 * is an isomorphism from 𝑊 1,𝑠 0 (Ω) to 𝑉 ⊥ , so that there exists a unique 𝑝 ∈ 𝑊

Energetic formulation

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

𝑔 𝑣 ∈ 𝜕ℱ(𝑣) ⇐⇒ 𝑔 𝑣 ∈ 𝑉 * and ∀ 𝜙 ∈ 𝑉, lim inf 𝛿→0 + ℱ(𝑣 + 𝛿𝜙) -ℱ(𝑣) 𝛿 ⟨𝑔 𝑣 , 𝜙⟩.
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

Ψ(𝑎) = { Φ 1 (𝑎) for all 𝑎 1, Φ 2 (𝑎) for all 𝑎 > 1, (4.2)
where

Φ 1 , Φ 2 : [0, ∞) → R
𝛿→0 + 𝒟(𝑣 + 𝛿𝜙) -𝒟(𝑣) 𝛿 ⟨𝑔 𝑣 , 𝜙⟩.
Write 𝑉 1 the subset of 𝑉 consisting of the functions which are supported in Ω 1 (𝑣 + û). Then, for all 𝜙 ∈ 𝑉 1 , since -𝜙 ∈ 𝑉 1 as well, we get lim inf

𝛿→0 + 1 𝛿 ∫ Ω1(𝑣+ û) ( Ψ ( ‖𝑣 + 𝛿𝜙 + û‖ 2 ) -Φ 1 ( ‖𝑣 + û‖ 2 )) ⟨𝑔 𝑣 , 𝜙⟩ and lim inf 𝛿→0 + 1 𝛿 ∫ Ω1(𝑣+ û) ( Ψ ( ‖𝑣 -𝛿𝜙 + û‖ 2 ) -Φ 1 ( ‖𝑣 + û‖ 2 )) -⟨𝑔 𝑣 , 𝜙⟩,
which, by Lebesgue's dominated convergence theorem and the differentiability of Φ 1 , yield

∫ Ω Λ 1 (𝑣 + û) • 𝜙 = ∫ Ω1(𝑣+ û) 𝜑 1 ( ‖𝑣 + û‖ 2 ) (𝑣 + û) • 𝜙 ⟨𝑔 𝑣 , 𝜙⟩ and - ∫ Ω Λ 1 (𝑣 + û) • 𝜙 = - ∫ Ω1(𝑣+ û) 𝜑 1 ( ‖𝑣 + û‖ 2 ) (𝑣 + û) • 𝜙 -⟨𝑔 𝑣 , 𝜙⟩.
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

⟨𝑔 𝑣 , 𝜙⟩ = ⟨Λ 2 (𝑣 + û), 𝜙⟩ for all 𝜙 ∈ 𝑉 2 .
Thus the function defined by

Λ 𝑣 = Λ 1 (𝑣 + û) + Λ 2 (𝑣 + û) + 𝑔 𝑣 𝜒 Γ(𝑣+ û)
is such that ⟨Λ 𝑣 , 𝜙⟩ = ⟨𝑔 𝑣 , 𝜙⟩ for all 𝜙 ∈ 𝑉 . Therefore, from Problem 4.7 we obtain ⟨Λ 𝑣 , 𝜙⟩ = ⟨𝑓 0 , 𝜙⟩ for all 𝜙 ∈ 𝑉 .

Moreover Λ 𝑣 ∈ Λ(𝑣 + û), where we recall Λ is in (3.5). Hence 𝑣 is solution to Problem 4.3.

We now define the energy ℰ : 𝑉 → R associated to 𝒟 by ℰ(𝑣) = 𝒟(𝑣) -⟨𝑓 0 , 𝑣 + û⟩ for all 𝑣 ∈ 𝑉 . Problem 4.9 (Energy minimization). Find 𝑣 ∈ 𝑉 such that 𝑣 ∈ 𝑀 ℰ .

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.

Case 𝜆 1 𝜆 2

The result we wish to show here is the following: for all 𝑎 1, 𝜑 2 (𝑎) for all 𝑎 > 1.

Theorem 4.
Since 𝜑 1 and 𝜑 2 are increasing on [0, 1] and [1, ∞), respectively, and 𝜑 1 (1) = 𝜆 1 and 𝜑 2 (1) = 𝜆 2 , the derivatives Φ′ 1 and Φ′ 2 are also increasing so that Φ1 and Φ2 are convex. Furthermore, because 𝜆 1 𝜆 2 and 𝜑 2 is increasing, we have Φ′ 2 (𝑎) Φ′ 1 (𝑎) for all 𝑎 > 1; thus, the fact that Φ 1 (1) = Φ 2 (1) (and so Φ1 (1) = Φ2 (1)) yields Φ1 (𝑎) Φ2 (𝑎) for all 𝑎 > 1. Similarly, we get that Φ1 (𝑎) Φ2 (𝑎) for all 𝑎 < 1. Let 𝑎, 𝑏 ∈ [0, ∞) and 𝑡 ∈

[0, 1]. If 𝑎, 𝑏 1, then (1 -𝑡)𝑎 + 𝑡𝑏 1 and Ψ((1 -𝑡)𝑎 + 𝑡𝑏) = Φ 1 ((1 -𝑡)𝑎 + 𝑡𝑏) = Φ1 ((1 -𝑡)𝑎 + 𝑡𝑏) (1 -𝑡) Φ1 (𝑎) + 𝑡 Φ1 (𝑏) = (1 -𝑡)Ψ(𝑎) + 𝑡Ψ(𝑏),
and similarly if 𝑎, 𝑏 > 1. If now 𝑎 1, 𝑏 > 1 and (1 -𝑡)𝑎 + 𝑡𝑏 1, then we have

Ψ((1 -𝑡)𝑎 + 𝑡𝑏) = Φ 1 ((1 -𝑡)𝑎 + 𝑡𝑏) = Φ1 ((1 -𝑡)𝑎 + 𝑡𝑏) (1 -𝑡) Φ1 (𝑎) + 𝑡 Φ1 (𝑏) (1 -𝑡) Φ1 (𝑎) + 𝑡 Φ2 (𝑏) = (1 -𝑡)Ψ(𝑎) + 𝑡Ψ(𝑏),
and similarly if (1 -𝑡)𝑎 + 𝑡𝑏 > 1 or 𝑎 > 1 and 𝑏 1. In all cases, we see that Ψ((1 -𝑡)𝑎 + 𝑡𝑏) (1 -𝑡)Ψ(𝑎) + 𝑡Ψ(𝑏), so that Ψ is convex. Note that one could reach the same conclusion using that Ψ(𝑎) = max( Φ1 (𝑎), Φ2 (𝑎)) for all 𝑎 0.

Since 𝜑 1 and 𝜑 2 are positive on (0, 1) and [1, ∞), respectively, we get that Ψ is increasing. Therefore, the function Ψ ∘ ‖•‖ 2 is strictly convex.

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

𝐾 > ℰ(𝑣 𝑛 ) = ∫ Ω1(𝑣𝑛+ û) Φ 1 ( ‖𝑣 𝑛 + û‖ 2 ) + ∫ Ω2(𝑣𝑛+ û) Φ 2 ( ‖𝑣 𝑛 + û‖ 2 ) - ∫ Ω 𝑓 0 • (𝑣 𝑛 + û) Φ 1 (0)|Ω| + 𝑐 𝑟 ∫ Ω2(𝑣𝑛+ û) ‖𝑣 𝑛 + û‖ 𝑟 -‖𝑓 0 ‖ 𝑠 𝐿 𝑠 (Ω) ‖𝑣 𝑛 + û‖ 𝑟 𝐿 𝑟 (Ω) ( Φ 1 (0) - 𝑐 𝑟 ) |Ω| + 𝑐 𝑟 ‖𝑣 𝑛 + û‖ 𝑟 𝐿 𝑟 (Ω) -‖𝑓 0 ‖ 𝐿 𝑠 (Ω) ‖𝑣 𝑛 + û‖ 𝐿 𝑟 (Ω) ,
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 𝑤∈𝑉 ℰ(𝑤) = lim inf 𝑛→∞ ℰ(𝑣 𝑛 ) ℰ(𝑣) 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

𝒟((1 -𝑡)𝑣 + 𝑡𝑤) = ∫ Ω Ψ ( ‖(1 -𝑡)𝑣 + 𝑡𝑤 + û‖ 2 ) = ∫ Ω Ψ ( ‖(1 -𝑡)(𝑣 + û) + 𝑡(𝑤 + û)‖ 2 ) = (1 -𝑡) ∫ Ω∖𝐴 Ψ ( ‖𝑣 + û‖ 2 ) + 𝑡 ∫ Ω∖𝐴 Ψ ( ‖𝑤 + û‖ 2 ) + ∫ 𝐴 Ψ ( ‖(1 -𝑡)(𝑣 + û) + 𝑡(𝑤 + û)‖ 2 ) < (1 -𝑡) ∫ Ω∖𝐴 Ψ ( ‖𝑣 + û‖ 2 ) + 𝑡 ∫ Ω∖𝐴 Ψ ( ‖𝑤 + û‖ 2 ) + (1 -𝑡) ∫ 𝐴 Ψ ( ‖𝑣 + û‖ 2 ) + 𝑡 ∫ 𝐴 Ψ ( ‖𝑤 + û‖ 2 ) = (1 -𝑡) ∫ Ω Ψ ( ‖𝑣 + û‖ 2 ) + 𝑡 ∫ Ω Ψ ( ‖𝑤 + û‖ 2 ) = (1 -𝑡)𝒟(𝑣) + 𝑡𝒟(𝑤),
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.

Proof. Without loss of generality, take Ω = (0, 1).

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

𝜓(𝑥) = ∫ 𝑥 0 𝜂 for all 𝑥 ∈ (0, 1)
satisfies 𝜓 ′ = 𝜂 and 𝜓(0) = 𝜓(1) = 0. Let now 𝑣 ∈ 𝑉 and define 𝜂 ∈ 𝑉 as

𝜂(𝑥) = 𝑣(𝑥) - ∫ 1 0
𝑣 for all 𝑥 ∈ (0, 1), so that obviously ∫ 1 0 𝜂 = 0. Write 𝜓 ∈ 𝑊 1,𝑠 0 ((0, 1)) a primitive of 𝜂 and compute, for all 𝜙 ∈ 𝑉 ,

∫ 1 0 ( 𝑣 - ∫ 1 0 𝑣 ) 𝜙 = ∫ 1 0 𝜂𝜙 = ∫ 1 0 𝜓 ′ 𝜙 = 0.
Thus 𝑣 = ∫ 1 0 𝑣 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:

𝐸(𝛼) = ∫ 1 0 Ψ ( (𝛼 + û) 2 ) -𝛼 ∫ 1 0 𝑓 0 for all 𝛼 ∈ 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 inf

𝛽∈R 𝐸(𝛽) = lim inf 𝑛→∞ 𝐸(𝛼 𝑛 ) 𝐸(𝛼) inf 𝛽∈R 𝐸(𝛽),
showing that 𝛼 is a global minimizer of 𝐸 and the constant function 𝛼 belongs to 𝑀 ℰ .

Case Vol 0 (Σ v ) > 0. Note that here Σ p ∈ {∅}∪{{0}}∪{{1}}. 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 𝜓 ′ 𝜙 = 0, 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.

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 𝜓(𝑥) = ∫ 𝑥 0 𝑣 - ∫ 1 0 (∫ 𝑦 0 𝑣 ) d𝑦 if Σ p = ∅, 𝜓(𝑥) = ∫ 𝑥 0 𝑣 - ∫ 1 0 𝑣 if Σ p = {1}, 𝜓(𝑥) = ∫ 𝑥 0 𝑣 if Σ p = {0},

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 [START_REF] Keilegavlen | Porepy: an open-source software for simulation of multiphysics processes in fractured porous media[END_REF].

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

𝑢 (𝑖) , 𝑝 (𝑖) ) ← PDE(𝒞 (𝑖-1) ℎ ); for 𝐸 = (𝑒1, 𝑒2) ∈ Ω (𝑖-1) ℎ do (𝑐1, 𝑐2) ← ( ‖𝑢 (𝑖) (𝑒1)‖ < 1, ‖𝑢 (𝑖) (𝑒2)‖ < 1 ) ; if 𝑐1 ̸ = 𝑐2 then 𝒞 (𝑖) ℎ in 𝐸 ← Transition zone(𝑢 (𝑖) ); end end 𝑑 ← Distance(Γ (𝑖) ℎ , Γ (𝑖-1) ℎ ); 𝑖 ← 𝑖 + 1; while 𝑑 > 𝜖Ω;

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 [START_REF] List | A study on iterative methods for solving Richards' equation[END_REF][START_REF] Pop | Mixed finite elements for the Richards' equation: linearization procedure[END_REF][START_REF] Radu | A robust, mass conservative scheme for two-phase flow in porous media including Hölder continuous nonlinearities[END_REF] 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 [START_REF] Formaggia | A reduced model for Darcy's problem in networks of fractures[END_REF], 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 [START_REF] Raviart | A mixed finite element method for second order elliptic problems[END_REF][START_REF] Roberts | Mixed and hybrid methods[END_REF]. 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 [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF][START_REF] Formaggia | A reduced model for Darcy's problem in networks of fractures[END_REF], 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 [START_REF] Flemisch | Benchmarks for singlephase flow in fractured porous media[END_REF] in Section 6.2.3.

These examples were developed with the open source library PorePy [START_REF] Keilegavlen | Porepy: an open-source software for simulation of multiphysics processes in fractured porous media[END_REF]. The associated scripts are freely accessible. PorePy uses Gmsh [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF] 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

𝜇𝑘 -1 𝑢 = -∇𝑝 and 𝜇𝑘 -1 𝑢 + 𝜌𝛽‖𝑢‖𝑢 = -∇𝑝,
where 𝜇 is the viscosity of the liquid in [Pa s], 𝑘 the scalar intrinsic rock permeability in [m 2 ], 𝜌 the liquid density [kg m -3 ] and 𝛽 the so-called Forchheimer coefficient in [m -1 ]. We consider water as liquid so 𝜇 = 10 -3 [Pa s] and 𝜌 = 10 3 [kg m -3 ]. Following [START_REF] Zeng | A criterion for non-Darcy flow in porous media[END_REF], we can define the Forchheimer number 𝐹 𝑜 to discriminate when non-Darcy effects should be included in the model or not, namely, ]. We can estimate a fluid velocity to be in the order of 10 -4 and with this value the Forchheimer number is estimated as

𝐹 𝑜 = 𝑘𝛽𝜌‖𝑢‖ 𝜇 ; for 𝐹 𝑜 < 0.
𝐹 𝑜 (𝑥) = { 0.5 for 𝑥 50[m], 0.05 for 𝑥 > 50[m].
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

Ω 1 = {𝑥 ∈ Ω | ‖𝑢(𝑥)‖ < 𝑢}, Ω 2 = {𝑥 ∈ Ω | ‖𝑢(𝑥)‖ > 𝑢}, Γ = {𝑥 ∈ Ω | ‖𝑢(𝑥)‖ = 𝑢}.
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 Ω (0) 1 = Ω, 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 𝐾 -1 𝑢 = -∇𝑝 + 𝑓 with 𝐾 = 𝑘𝐼, where 𝐼 is the identity matrix and (With the notation used previously, 𝑘 𝑖 = 1/𝜆 𝑖 , 𝑖 ∈ {1, 2}.) Specifically, we have

𝑘 = { 𝑘 1 = 1 in Ω 1 , 𝑘 2 = 10 in Ω 2 . ( 6
Λ(𝑢) = { 𝑢 in Ω 1 , 0.1𝑢 in Ω 2 . ( 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 Ω 1 and a Darcy-Forchheimer flow in Ω 2 . Specifically, we have

Λ(𝑢) = { 𝑢 in Ω 1 , (0.01 + 3‖𝑢‖)𝑢 in Ω 2 . (6.3)
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

𝑞(𝑥) = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if 𝑥 1 0.3, -1 if 0.3 < 𝑥 1 < 0.7, 1 if 𝑥 1 0.7,
and 𝑓 = [5 • 10 -2 , 0, 0] ⊤ .

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 (that is, 𝜆 2 ) for the linear example in Section 6.2.1. The red asterisks mean that the maximum number of iterations is 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 and4 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. 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

err p = ‖𝑝 ref -𝑝‖ ‖𝑝 ref ‖ and err u = ‖𝑢 ref -𝑢‖ ‖𝑢 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.

Crossing domains

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

𝑞 horiz (𝑥) = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if 𝑥 1 0.3 -1 if 0.3 < 𝑥 1 < 0.7 1 if 𝑥 1 0.7 and 𝑞 vert (𝑥) = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if 𝑥 2 0.3 -1 if 0.3 < 𝑥 2 < 0.7 1 if 𝑥 2 0.7
.

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

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 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 Ω, while on the bottom for the vertical one. 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 [START_REF] Ahmed | A reduced fracture model for two-phase flow with different rock types[END_REF][START_REF] Audu | On the well-posedness of generalized Darcy-Forchheimer equation[END_REF][START_REF] Berrone | A PDE-constrained optimization formulation for discrete fracture network flows[END_REF][START_REF] Berrone | On simulations of discrete fracture network flows with an optimization-based extended finite element method[END_REF]. 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 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. 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 [START_REF] Flemisch | Benchmarks for singlephase flow in fractured porous media[END_REF]. 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 [START_REF] Ahmed | A reduced fracture model for two-phase flow with different rock types[END_REF][START_REF] Amir | Décomposition de domaine et préconditionnement pour un modèle 3D en milieu poreux fracturé[END_REF][START_REF] Angot | Well-posed Stokes/Brinkman and Stokes/Darcy coupling revisited with new jump interface conditions[END_REF], 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. 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.

Problem 4 . 3 ( 2 . 4 . 4 .

 43244 Weak form restricted to 𝑉 ). Find 𝑣 ∈ 𝑉 such that there exists Λ 𝑣 ∈ Λ(𝑣 + û) satisfying⟨Λ 𝑣 , 𝜙⟩ = ⟨𝑓 0 , 𝜙⟩ for all 𝜙 ∈ 𝑉 ,where û is given by Lemma 4.Lemma Problem 3.5 is equivalent to Problem 4.3.Proof. We first suppose that (𝑢, 𝑝) is solution to Problem 3.5. Then we decompose 𝑢 as 𝑢 = (𝑢-û)+ û =: 𝑣+ û. By Problem 3.5, there exists Λ 𝑣 ∈ Λ(𝑢) = Λ(𝑣 + û) such that ⟨Λ 𝑣 , 𝜙⟩ = ⟨𝑓 0 , 𝜙⟩ for all 𝜙 ∈ 𝑉 .

  1,𝑠 0 (Ω) with ⟨Λ 𝑣 -𝑓 0 , 𝜙⟩ = 𝐵 * (-𝑝)(𝜙) = -⟨∇𝑝, 𝜙⟩ for all 𝜙 ∈ 𝐿 𝑟 (Ω). Furthermore, writing 𝑢 = 𝑣 + û we get ⟨∇𝜓, 𝑢⟩ = ⟨∇𝜓, 𝑣⟩ + ⟨∇𝜓, û⟩ = ⟨∇𝜓, û⟩ = -∫ Ω 𝑞𝜓 + ∫ Σv 𝑢 0 𝜓, since 𝑣 ∈ 𝑉 . We thus have that (𝑢, 𝑝) satisfies Problem 3.5, with 𝑝 unique.

Definition 4 . 5 (

 45 Fréchet subdifferential and strong local minimizer). Let ℱ : 𝑉 → R. For all 𝑣 ∈ 𝑉 we define the (Fréchet) subdifferential 𝜕ℱ(𝑣) of ℱ at 𝑣 by:

(4. 4 )

 4 Let us write 𝑀 ℰ ⊂ 𝑉 the, possibly empty, set of local minimizers of ℰ. Consider the following minimization problem.

satisfies

  𝜓 ′ = 𝜂, and ∫ 1 0 𝜓 = 0 if Σ p = ∅ and 𝜓(𝑥) = 0 for 𝑥 ∈ Σ p otherwise. Let 𝑣 ∈ 𝑉 , write 𝜓 a primitive of 𝑣 and compute, for all 𝜙 ∈ 𝑉 ,

Algorithm 1 :

 1 Transition zone tracking algorithm. Data: For 𝑖 = 1: the configuration 𝒞 (𝑖-1) ℎ and the tolerances (𝜖Γ, 𝜖Ω); Result: For 𝑖 = end: the configuration 𝒞 (𝑖) ℎ ; do (

  1 a pure Darcy model can be considered, while for 𝐹 𝑜 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 𝑘(𝑥) = { 5 • 10 -11 [m -2 ] for 𝑥 50, 5 • 10 -12 [m -2 ] for 𝑥 > 50. A pressure gradient is set on the boundary with 𝑝(𝑥 = 0[m]) = 10 6 [Pa] and 𝑝(𝑥 = 100[m]) = 2 • 10 6 [Pa], a scalar source term 𝑞 = 10 -6 [s -1 ] is set in the interval (30[m], 70[m]) and the Forchheimer coefficient is taken as 𝛽 = 10 8 [m -1

. 1 )Figure 1 .

 11 Figure 1. Solutions for the problem of Section On the left, using the proposed adaptive model; in the center, with the fixed flow regime; on the right, the error between them.

Figure 2 .

 2 Figure2. 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 Ω that belongs to Ω 1 or not. The pressure profile is amplified by a factor of 10.

Figure 3 .

 3 Figure 3. Number of iterations of Algorithm 1 for different values of the inverse of 𝑘 2 (that is, 𝜆 2 ) for the linear example in Section 6.2.1. The red asterisks mean that the maximum number of iterations is

Figure 4 .

 4 Figure 4. The two solutions computed by the algorithm for 𝑘 2 = 0.25; see the linear case in Section 6.2.1.

Figure 5 .

 5 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 Ω that belongs to Ω 1 or not. The pressure profile is amplified by a factor of 100.

Figure 8 .

 8 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 Ω 𝑠 .

Figure 11 .

 11 Figure 11. Zoom around the short fracture branches Ω 𝑠 of the "if Ω 1 " condition for different iterations, respectively[START_REF] Ahmed | A reduced fracture model for two-phase flow with different rock types[END_REF][START_REF] Amir | Décomposition de domaine et préconditionnement pour un modèle 3D en milieu poreux fracturé[END_REF][START_REF] Angot | Well-posed Stokes/Brinkman and Stokes/Darcy coupling revisited with new jump interface conditions[END_REF], for the problem of Section 6.2.3.

  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: Problem 4.7 (Dissipation form). Find 𝑣 ∈ 𝑉 such that there exists 𝑔 𝑣 ∈ 𝜕𝒟(𝑣) satisfying ⟨𝑔 𝑣 , 𝜙⟩ = ⟨𝑓 0 , 𝜙⟩ ∀𝜙 ∈ 𝑉. Lemma 4.8. Any solution to Problem 4.7 is also solution to Problem 4.3.

	are primitives of 𝜑1 2 and 𝜑2 2 (cf. (3.2)) such that Φ 1 (1) = Φ 2 (1) = 0. We define the
	dissipation 𝒟 : 𝑉 → R by
	𝒟(𝑣) =	∫	Ψ	(	‖𝑣 + û‖ 2	)	for all 𝑣 ∈ 𝑉 .	(4.3)
		Ω
	The following result holds:

Proof. Suppose that 𝑣 ∈ 𝑉 is solution to Problem 4.7. Then we can pick 𝑔 𝑣 ∈ 𝜕𝒟(𝑣) ⊂ 𝐿 𝑠 (Ω) so that Problem 4.7 holds. By definition, for all 𝜙 ∈ 𝑉 we must have lim inf

  10 (Existence and uniqueness when 𝜆 1 𝜆 2 ). Suppose that 𝜆 1 𝜆 2 . Then, Problem 4.9 has a unique solution.

	Proof. Let us prove that the integrand Ψ ∘ ‖•‖	2 (cf. (4.2)) of 𝒟 is strictly convex. Define the functions
	Φ1 , Φ2 : [0, ∞) → R by
	Φ1 (𝑎) =	{	Φ 1 (𝑎)	for all 𝑎 1,
	Then Φ1 and Φ2 are differentiable with
	Φ′ 1 (𝑎) =	1 2	{ 𝜑 1 (𝑎) for all 𝑎 1, 𝜆 1 for all 𝑎 > 1,	and Φ′ 2 (𝑎) =	2 1	{	𝜆 2

𝜆1 2 (𝑎 -1) for all 𝑎 > 1, and Φ2 (𝑎) = { 𝜆2 2 (𝑎 -1) for all 𝑎 1, Φ 2 (𝑎) for all 𝑎 > 1.

Table 1 .

 1 Numbers of iterations and errors computed for the example in Section 6.2.1, nonlinear case.

	𝜖nl	itout itin errp	erru
	10 -1	2	3	1.2 • 10 -5	8.7 • 10 -6
	10 -2	2	4	1.8 • 10 -7	1.3 • 10 -7
	10 -3	2	5	2.8 • 10 -9	2.2 • 10 -9
	10 -4	2	6	4.4 • 10 -11 2.1 • 10 -11
	10 -5	2	7	2.4 • 10 -12 1.0 • 10 -11
	10 -8	2	11	0	0
	10 -12 2	15	-	-

https://www. edpsciences.org/en


Arnaud Disant
email: arnaud.disant@ucd.ie

Frederic Dias
email: 2frederic.dias@ucd.ie

MICROWAVE PROPAGATION IN MARITIME ENVIRONMENT

Keywords: Real-time data, WiFi at sea, microwaves, antennas, data transmission

A simple question that arises when dealing with maritime communications is: How would one offload large quantities of data from sea to shore and vice versa if one cannot use conventional solutions such as satellite communications or cellular data? In this note we describe the first prototype solutions that were produced. They gave rise to SeaFi, which has become an enabling technology in the context of oceanic and coastal research. Measurements at sea are key to scientific research project such as for exemple HIGHWAVE, a recently started ERC Advanced Grant project that relies partly on the possibility to measure breaking waves in realtime. Since the bottleneck is the real-time transmission of data, transferring measurement data at sea using SeaFi instead of using conventional satellite communications or cellular data connection, quickly became an evidence. To assess the possibilities, a series of offshore experiments were performed from May to July 2018. Those experiments led on June 6 th , 2018 to a world record for the longest wireless microwave transmission at sea, between a moving ship and a lighthouse, using the SeaFi communication system. In addition, the proposed solution also promises a better future for lighthouses around the world that are now gradually falling into disuse. This breakthrough in maritime telecommunications could change the way scientific researchers retrieve data in real time at sea in the coming years.

INTRODUCTION

Whereas some data buoys or oil and gas rigs can be connected by optical fiber laid out on the seafloor, there are realistically only two options at present when dealing with large datasets on a ship/deep water data buoy at sea: satellite or cellular. But neither was designed to deal with large scale marine datasets. Integration of satellite communications on a buoy or an Unmanned Surface Vehicle (USV) is costly. Cellular data is only sporadically available at sea. SeaFi, with wireless microwave transmission, offers a new third option specifically designed for large-scale marine data environments. The large-scale offshore experiment described in this paper, which took place on May 5 th and June 6 th , 2018, is a proof of concept.

In telecommunications, wide bandwidth data transmissions conveying multiple signals are also referred to as broadband. Data transmission refers to a digital bitstream over a point-to-point or even point-to-multipoint communications channel. The carrier can be a coaxial cable, optical fibre, radio or a twisted pair of copper cables. In wireless communication channels, the data are represented as an electromagnetic signal, such as an electrical voltage, radio wave, microwave, or infrared signal.

The first SeaFi prototype, developed in 2012, simply aimed to connect a local service vessel (tug) to the internet using conventional off the shelf WiFi equipment. It quickly proved not to be suited to the maritime environment. It was too fragile and had a poor tilting tolerance at sea. The prototype evolved into a rugged seaworthy piece of equipment designed for ships and data buoys.

With support from the Irish Naval Service, sea trials involving Offshore Patrol Vessels (OPVs) and Rigid Inflatable Boats (RIBs) showed greater potential. This is when we decided to build our own equipment. The system would be drawing from Wireless Internet Service Providers (WISP) standards but, additionally, require more robust antennas on shore. By 2014, the concept of a Wireless Maritime Area Network (WMAN) was born and, with it, a new way of building smarter cost-effective ports and intelligent coastlines, as well as offshore renewable and salvage operations, without compromising throughput, security or data transfers.

Beyond scientific data collection, the applications for SeaFi can also be of societal nature. In 2017 an experiment was conducted aboard the MV Stena Horizon to demonstrate the superiority of WMAN in a port. It has since caught the attention of the maritime community in the UK and Ireland. Besides improving security and business efficiency and providing much lower communication costs, SeaFi offers the possibility of greatly improving links of sea workers with their families during long absences.

The goal of the large-scale offshore experiment was to assess the resilience of a wireless microwave link while determining the furthest practical maritime broadband transmission without a satellite or cellular network, using the SeaFi wireless maritime communication system between a moving ship and a lighthouse. The SeaFi experiment of June 2018 demonstrated the superiority in terms of reliability and running costs of the SeaFi wireless microwave link over conventional technologies such as satellite or cellular data communications. This experiment brought permanent access to the ship equipment from the shore side (remote connection). Using an Offshore Supply Vessel (OSV) operating on the same route, instead of an OPV, allowed controlled observations depending on:

1.

Coastal In this note, we first describe the various communication methodologies. We then provide an account of the experiment. Finally, we discuss the outcomes of the experiment.

COMMUNICATION METHODOLOGIES

Cellular networks and satellite communications

Invented to satisfy the growing demand for data from mobile phones, the 3G/4G cellular network is expensive to put in place, and so is primarily used in larger population areas. The availability of 3G/4G at sea is purely coincidental and unreliable, as it was never designed for this type of use.

So, when networks become overloaded (saturated networks phenomenon), the units most remote from the mast are generally the first to be disconnected from the network. Cellular was originally invented for mobility with station roaming from cell to cell, for example a phone in a car travelling on the motorway.

Cellular networks data capacity is often increased by splitting cells to reduce the size, and thus increase the number of cells serving a given area. This is done by deploying more radio towers/antennas and shrinking the reach of each tower by reducing the radiated power of its radio transmissions. This allows the radio spectrum to be reused for multiple simultaneous transmissions within the geographic area. It also increases the amount of traffic that a Hz of spectrum can carry within an overall geographic area (measured by bps/km 2 ). But while effective at deepening wireless network capacity, this method is also expensive -requiring the construction of extra towers/antennas, deploying more radios and base station equipment, as well as extending additional backhaul links to link new towers back into the mobile operator's core network [START_REF] Clarke | Expanding mobile wireless capacity: The challenges presented by technology and economics[END_REF].

Specifically designed for maritime use, the Iridium satellite system is theoretically based on 80 satellites. The first generation of satellites that was put in orbit between 1997 and 2002 (at a cost of $5 billion) is now reaching end of life in orbit with the new generation satellites taking over (Iridium Satellite Constellation, 2019). Technically, a satellite is a WiFi system where the base station access point (the satellite) is 700 km away from the vessels (the station) allowing communication under the horizon. The satellite concept designed when smartphones did not exist must now face the bottlenecked era of the all-connected Internet of Things (IoT). More recent projects such as Amazon, Starlink and OneWeb are promoting Low Orbit Satellites, which together with Nanosatellites (satellite weighing less than 10 kg), seem to be the most popular commercial solution in the next decade [START_REF] Boiardt | Low earth orbit nanosatellite communications using Iridium's network[END_REF].

SeaFi in a nutshell

SeaFi facilitates the creation of private networks in ports and coastal areas by enabling microwave connections without relying on conventional telecommunications infrastructures (Figure 1). These networks can be used to connect, for example, lighthouses, maritime wind turbines, offshore drilling platforms, sea tugs, gantry cranes, data collection buoys and wave energy converters. But

SeaFi can also work in ship-to-ship mode, and be combined with satellite communications to connect, for example, ships during salvage operations or ships with oil and gas platforms.

short-wave radio band, which is much lower in frequency. SeaFi is operating on the 2.4 Gigahertz (GHz) frequency, using interconnected coastal stations (connection access points) and embarked stations at sea. This frequency is gradually being abandoned on shore because of its channels limitation and hardly used at sea. Fourteen channels are designated in the 2.4 GHz range, spaced 5 MHz apart from each other. SeaFi is based on industry standards like the 802.11 wireless standard that was developed by the Institute of Electrical and Electronic Engineers, while its antennas are designed to withhold the harshest conditions of maritime environment, including humidity, salinity, rain fade, rolling, pitching and wind.

THE EXPERIMENT

The purpose of the experiment

The main purpose of the experiment was to assess the longest possible distance at which the signal from SeaFi's bespoke maritime WiFi service is still reliable. Lining up irrefutable evidence of an invisible point at sea where intangible events occurred on a ship is not easy. A large-scale fieldtesting environment involving a series of coastal stations was setup on the south coast of Ireland.

The experiments reported in this document focus essentially on the coastal station located on Roches Point lighthouse 25 m above sea level and one ship station located south of the lighthouse on board the OSV Ocean Spey.

The coastal station technical layout

The lighthouse of Roches Point was first connected to the internet in 2013. The internet connection is provided by SeaFi Port of Cork WMAN (Figure 1). The WMAN network is composed of a series of 10 interconnected coastal stations. The main junction point with the internet is in a secured bastion on Mitchell's Fort on Spike Island exactly 2.75 nautical miles (NM) away. The C180 radio transmission equipment is located on the south side of Roches Point lighthouse below the lantern to avoid any obstruction to the light quadrants 90°-210° (120°).

The ship station technical layout

The ship station was at the base of the mast in a Faraday cage-type compartment such that the only way in or out for the radio signal was the SeaFi S360 antenna, positioned approximately 25 m above sea level (Figure 2). The SeaFi ship station was composed of a firewall, a radio and an omnidirectional antenna. The smart cabinet enclosing all the transmission was not radio transparent and located in a metal fireproof bunker (battery compartment with no windows). This cabinet was closed and did not let radio waves through. As such, the equipment couldn't receive the signal of a GPS nor any cellular data for that matter.

Before the experiment

To define the furthest point at sea where the internet connection is still useable, while a ship is moving at a speed of 10 knots and radio waves are travelling at the speed of light, meaningful tests had to be designed. The OSV Ocean Spey was connected to SeaFi Port of Cork WMAN at 14:45 on May 9th, 2018 for the shakedown cruise. The fine tuning of the antenna was done at sea. The software had been pre-configured in the workshop. Witnesses were carefully appointed, to encompass a wide array of skills from mathematics to technology. A few hours before the experiment on June 6th, the witnesses inspected the SeaFi ship station on board the OSV Ocean Spey in the Ringaskiddy Deep Water Berth. The bridge and the mast of the ship were searched for any trace of internet connectivity other than the SeaFi radio link.

During the experiment

The observations of the witnesses were recorded independently from each other, which allowed cross-verification of observations. At exactly 14:00 on June 6th, At exactly 16:31, Roches Point issued one last message to the ship station. At 16:32, the SeaFi ship station on the Ocean Spey went silent after 2.5 hours of uninterrupted transmission. The OSV Ocean Spey was identified at this time by the Irish Naval Service vessel LÉ James Joyce, which of 10 -5 cm and gamma rays on the order of 10 -10 cm or shorter, all propagate at the same velocity through a common medium.

For two high-orbit satellites, for example, microwave transmission consists in pointing the antennas towards each other and launching the wave into space. Since the wave obeys the wellknown Friis equation [START_REF] Friis | A Note on a Simple Transmission Formula[END_REF][START_REF] Kopp | Microwave and Millimetric Wave Propagation[END_REF], its intensity is inversely proportional to the square of the distance. On earth, more parameters come into play to affect the signal. Low-frequency radio waves travel easily through brick and stone and very low frequency (VLF) waves even penetrate seawater. As the frequency rises, absorption effects become more important. At microwave or higher frequencies, absorption by molecular resonances in the atmosphere (mostly from water and oxygen) is a major factor in radio propagation at sea [START_REF] Liebe | An updated model for millimeter wave propagation in moist air[END_REF]. For example, in the 58-60 GHz band, there is a major absorption peak which makes this band useless for longdistance use (Federal Communications Commission, 1997). While all visible light, most radio light, and some infrared light passes through the atmosphere, the atmosphere blocks most ultraviolet light (UV) and all X-rays and gamma rays from reaching the surface of Earth. Heavy rain and falling snow also affect microwave absorption. On earth, radio waves are subject to distortion, diffraction, refraction and reflection. At sea level this is made even more complex by natural phenomena such as water absorption and sea surface reflection. [START_REF] Algor | Radio-wave propagation measurements over sea water[END_REF] Another parameter to take into consideration is the Fresnel space factor [START_REF] Guo | Fresnel Zone Antennas[END_REF].

The electromagnetic penetration is infinite in a perfectly transparent substance and zero in a reflective material (for example, metals). At the microwave frequency (2.45 GHz, λ = 12.24 cm), most energy is absorbed by water. This became more obvious as the ship reached the horizon. The line of sight (LoS) horizon for our antennas located at 25 m at sea and 25 m on land was approximately 36 km. As the distance between the two antennas increased, the Fresnel space diameter increased too, causing the signal to hit the sea. The loss of signal occurred at 19.4 NM (approximately 35.9 km). At this distance the Fresnel space radius was over 32.7 m, causing possible signal absorption by sea water too.

Figure 1 :

 1 Figure 1 : Wireless Maritime Area Network (WMAN) representing Irish Navy vessel Offshore Patrol Vessel (OPV) acknowledging transmission from Offshore Supply Vessel (OSV) using SeaFi.

Figure 2 :

 2 Figure 2 : The photo shows the OSV Ocean Spey passing in front of Roches Point lighthouse on May 9th, 2018, during the shakedown cruise. All the evidence used in this note is dated and all photos are also geolocated.

Figure 3 :

 3 Figure 3 : SeaFi operator consoles showing GPS positions and OSV moving at a speed of 9.6 kts with SeaFi enabled IP cameras, emails, instant messaging connection. OSV Ocean Spey SeaFi operator also performed a Skype call with several newspaper magazines in Ireland and UK.

  2018, the OSV Ocean Spey left port, setting course south of Roches Point lighthouse in Cork Harbor. For several hours, every 15 minutes, the operator on the bridge of the Ocean Spey sent emails to the witnesses and made Skype calls to recipients in Cork and London. The witnesses confirmed in real time that the emails originated from the SeaFi ship station. Each email sent to shore from the ship was delivered without satellite nor cellular connection. Some witnesses with ICT expertise performed tests such as trace route, IP lookup, after reading email headers information to verify the unique public IP address of the Port of Cork WMAN. Since the email test could be sent remotely, the team on board the Ocean Spey also decided to do undisputable real time Skype calls with witnesses in Ireland and UK. The calls were recorded and can be viewed in a video available on YouTube.

Figure 4 :

 4 Figure 4 : Emails sent are tagged with the public IP address of the sender, in the case of Port of Cork WMAN, the public IP address assigned by the Internet Service Provider is easily traceable.

  measured propagation loss at three frequencies(30, 140, and 412 MHz) over various sea-water paths out to 40 NM between a moderately elevated shore-based receiving site and a floating transmitter platform subject to wave motion. Experiments were conducted in order to determine the frequency region most conducive to trans horizon communications over a sea-water path of about 20 to 50 NM, considering various practical factors. A special requirement of the study was that one terminal be floating essentially at the water surface, and that the other terminal have a limited elevation ranging between 50 and 100 ft. Our own measurements of attenuation are shown in Figure3. An unexplained spike in the noise occurs at a distance of about 17 km (9.2 NM).Since most microwave links are used for terrestrial point-to-point or point-to-multipoint, manufacturers have failed to design an efficient alignment procedure for microwave equipment used at sea. It is almost impossible to line up a microwave antenna on a ship or a buoy because of the three-dimensional motion. For this reason, SeaFi S360 ship stations are different in their design from coastal base station antennas. For example they can compensate for pitch and roll, and they are designed to withstand the pressure of breaking waves.

Figure 6 :

 6 Figure 6 : Fresnel space factor and attenuation curves

ACKNOWLEDGEMENTS

The authors thank all the witnesses: Lt Cdr Martin Brett, Karl Grabe, Dr Sean McSweeney, Theophile Buyssens, Capt. Bill Kavanagh, Sean O'Callaghan, Capt. Mark Sloan and James Darwin. The contribution of Theophile Buyssens to prepare the curves of Figure 3 and the video log (SeaFi Offshore Experiment) posted on YouTube is also acknowledged. All experiments were performed while the second author was preparing the proposal of HIGHWAVE. Since then, HIGHWAVE has been approved for funding. Therefore, we wish to acknowledge the European Research Council (ERC-2019-AdG 833125-HIGHWAVE). SeaFi is an important part of HIGHWAVE. The second author also acknowledges the support of Science Foundation Ireland (SFI) under grant number SFI/12/RC/2302.

James Joyce confirmed the reporting of the position of the Ocean Spey at Lat 51 30.05N Lon 08 03.70W -19.4 NM (or 35.9 km) from Roches Point at the time the last transmission was made.

RESULTS

Although the main objective of the experiment was to assess the longest distance at which the microwave signal was still useable for data transmission, the signal attenuation was also studied.

All channels could be used for communications reliably to 22 NM or more with low transmitter power. It is evident from the following data, however, that there were some shortcomings due to setting up the antenna on the upper deck (Figure 6).

All electromagnetic radiation propagates at a velocity of 3 x 10 8 m/s in vacuum. Although wavelengths of radio waves are on the order of 10 3 m while those of light waves are on the order