Hichem Gargoubi
email: hichem.gargoubi@ipeit.rnu.tn

Sayed Kossentini
email: kossentinisayed@hotmail.fr

On Hyperbolic Integers

Keywords: Hyperbolic numbers, Hyperbolic Gaussian integers, f -ring, -group, Riesz space, Lattice. Mathematics Subject Classification : 11R04, 06F25

. We introduce the ring of hyperbolic integers Z h as a sub f -ring of the ring Z D of integers of D. We prove that Z h is the unique, up to ring isomorphism, Archimedean f -ring of quadratic integers. Our study focuses on arithmetic properties of Z h related to its lattice-ordered structure. We show that many of basic properties of the ring of integers Z such as primes, unique factorization theorem and the notions of floor and ceiling functions can be extended to Z h . A surprising fact is that prime numbers seen as hyperbolic integers are semiprimes. We also obtain some properties of hyperbolic Gaussian integers . As an application, we discuss the Dirichlet divisor problem using hyperbolic intervals.

Introduction

In order to solve problems concerning certain classes of integers, number theorists of the XIX century are led to study generalizations of the usual arithmetic of the natural numbers in more general settings. In these analogous of Z, concepts like unique factorization into prime elements, Euclidean division and modular arithmetic are developed. One can cite, for instance, Gaussian integers Z[i] and Kummer's cyclotomic integers Z[exp(2πi/n)] whose interest came about little by little.

In a process of generalization of the above constructions, Dedekind [START_REF] Dedekind | Theory of algebraic integers[END_REF] introduced the notion of ring of integers O K of a number field K. The ring of integers Z is the simplest ring of integers. Namely, Z = O Q where Q is the field of rational numbers. The ring of Gaussian integers Z[i] is the ring of integers of the number field Q(i). And, the ring of cyclotomic integers Z[exp(2πi/n)] is the ring of integers of the cyclotomic field Q(exp(2πi/n)).

Another structure of the ring of integers Z is useful, it corresponds to its order structure : it's clear that one has to take into account the relation between its divisibility and its order properties as the existence of a unique positive gcd and the existence of unique factorization into product of positive primes. In this paper Z is seen as an Archimedean f -ring. Indeed, the general notion of f -algebra is simultaneously a Riesz space (or vector lattice) and an associative real algebra that fulfills certain "positivity" conditions. A typical example of f -algebras is the linear space of real valued continuous functions on a topological space. Obviously, the fundamental example of Archimedean f -algebras is the field of real numbers.

The purpose of this paper is to answer the following general question: can we extend the order structure of Z to some of its generalizations in a way compatible with its arithmetic characteristics ?

The first element of answer is given by Theorem 3.1. We prove that there is no analogous of positivity in the ring O K . Actually, we prove a more general result concerning ring extensions of Q. This leads us to consider non division extension of reals numbers. More precisely, we consider the ring of hyperbolic numbers D = z = x + jy : x, y ∈ R, j / ∈ R; j 2 = 1 .

(1.1)

Hyperbolic numbers (also called duplex numbers) are an extension of real numbers defined in the same way that complex numbers C but with an imaginary unit j satisfying j 2 = 1 (instead of i 2 = -1). It's clear that D is not a division algebra. However, it enjoys an important order structure which makes it into the unique (up to an isomorphism) two dimensional Archimedean f -algebra. Therefore, basic notions of real analysis as sign, absolute value, Archimedean and Dedekind completeness are extended to hyperbolic numbers [START_REF] Gargoubi | f -Algebra Srtucture On Hyperbolic Numbers[END_REF]. Note that complex numbers and hyperbolic numbers are the only real commutative Clifford algebras:

D ∼ = Cl R (1, 0) and C ∼ = Cl R (0, 1)
The notion of partial order on D stimulates many authors and leads to interesting applications in different areas of mathematics. Alpay et al. [START_REF] Alpay | Basics of Functional Analysis with Bicomplex Scalars, and Bicomplex Schur Analysis[END_REF] investigated the D-normed bicomplex modulus. In probability theory it is shown in [START_REF] Alpay | Klomogrov's axioms for probabilities with values in hyperbolic numbers[END_REF] that Kolmogrov's axioms and Bays theorems holds in the context of D-valued probabilities. Kumar et al. [START_REF] Kumar | Orlicz Spaces with Bicomplex Scalars[END_REF] introduced the notion of D-valued measure on a sigma algebra. As an Application to fractal geometry, a concept of Cantor sets in hyperbolic numbers was developed by Balankin et al. [START_REF] Balankin | Cantor-type sets in hyperbolic numbers[END_REF] and Téllez-Sánchez et al. [START_REF] Téllez-Sánchez | More about Contor like sets in hyperbolic numbers[END_REF]. Recently, the authors of the present paper used in [START_REF] Gargoubi | Bicomplex numbers as a normal complexified f -algebra[END_REF] lattice-theoretical results to go further in the development of the theory of bicomplex zeta function. Further applications can be found in [START_REF] Kumar | Bicomplex Weighted Hardy Spaces and Bicomplex C * -algebra[END_REF][START_REF] Kumar | Topological Bicomplex modules[END_REF][START_REF] Luna-Elizarraras | On linear functionals and Hahn-Banach theorems for hyperbolic and bicomplex modules[END_REF][START_REF] Luna-Elizarraras | On the Laurent series for bicomplex holomorphic functions: Complex Var[END_REF][START_REF] Rochon | Bicomplex quantum mechanics II. The Hilbert Space[END_REF].

The aim of this paper is to investigate a lattice-ordered ring ( -ring) of algebraic integers where we are able to generalize many of the basic divisibility and order properties of Z. Our main result (Theorem 3.2) is the existence of a unique, up to ring isomorphism, Archimedean f -ring of quadratic integers called the ring of hyperbolic integers denoted Z h . Namely, Z h is the ring of integers of the extension Q(j) = {α + jβ ; α, β ∈ Q}.

The present paper is organized in the following way: In section 2, we recall some notions and terminology concerning -groups, f -rings and f -algebras and present basic notions and properties of hyperbolic numbers that will be used throughout this article. Section 3 introduces the lattice-ordered ring Z h of hyperbolic algebraic integers, and various of its properties are established. We introduce the notions of hyperbolic floor and ceiling functions which generalize that of real numbers. Section 4 and 5 are devoted to the divisibility in Z h . Many of the basic concepts of the arithmetic of Z are extended to Z h : the Euclidean division, the existence of a unique positive gcd, the existence of a unique factorization into a product of positive primes. In section 6, we establish some properties of the hyperbolic Gaussian integers as a subring of Z h . As an application, we discuss the Dirichlet divisor problem using hyperbolic intervals.

Preliminaries

In this section, we recall basic facts which we will use throughout this paper.

Basic lattice concepts

We recall some notions and terminology about -groups, f -rings [START_REF] Bigard | Groupes et Anneaux Réticulés[END_REF][START_REF] Steinberg | Lattice-ordered Rings and Modules[END_REF] and f -algebras [START_REF] Zaanen | Riesz Spaces II[END_REF]. Let G be a group which is also a partially ordered set. The group operation is denoted additively even if G is not assumed abelian, and so the identity element and the inverse of a ∈ G are denoted by 0 and -a, respectively. G is called a partially ordered group if the partial order ≤ satisfies: for any a, b ∈ G,

a ≤ b ⇒ a + c ≤ b + c and c + a ≤ c + b for all c ∈ G.
In the partially ordered group G, an element a is called positive if a ≥ 0. The set G + of all positive elements is called the positive cone of G. The partially ordered group G is said to be Archimedean if for each nonzero a in G the set {na : n ∈ Z} has no upper bound in G; equivalently, a, b ∈ G + and na ≤ b for all n ∈ N ⇒ a = 0. The partially ordered group G is a lattice-ordered group (an -group) if the partial order is a lattice order (i.e., the supremum a ∨ b and the infimum a ∧ b exist in G for all a, b ∈ G). Every element a in an -group G can be written as a = a + -a -, where a + = a ∨ 0 and a -= -a ∨ 0. The absolute value of a is defined as |a| = a ∨ (-a) = a + + a -. Any Archimedean -group is abelian. A real vectorial space V is said to be a vector lattice or Riesz space if V as a group is an -group satisfying the property: for any a, b ∈ V, a ≤ b ⇒ αa ≤ αb for all α ∈ R + . A ring R is called an f -ring if R is an -group and for any a, b ∈ R + ,

ab ∈ R + and a ∧ b = 0 ⇒ ac ∧ b = ca ∧ b = 0 for all c ∈ R + .
An f -ring is Archimedean if its underlying group is Archimedean. For each element a in an f -ring we have a + a -= 0. Two f -rings R and S are called -isomorphic if there exists a ring isomorphism ϕ from R to S satisfying ϕ(a∨b) = ϕ(a)∨ϕ(a) and ϕ(a∧b) = ϕ(a)∧ϕ(a) for all a, b ∈ R. An associative real algebra is an f -algebra if it is an f -ring and its underlying group is a vector lattice. We give above some examples of partially order groups, -groups and f -rings. 1. The additive group G = Z, Q, or R is an Archimedean totally -group with the usual order between real numbers, and |x| = x ∨ (-x) = max{x, -x} for all x ∈ G. 2. Let (G, P ) be the partially order group R × R with positive cone P .

• If P = {(x, y) : x > 0 or x = 0 and y ≥ 0}, then (G, P ) is a totally ordered group which is not Archimedean since for any n ∈ N, n(0, 1) ≤ (1, 0). The absolute value in (G, P ) is given by |(x, y)| = max{(x, y), (-x, -y)}.

• If P = {(x, y) : x > 0 and y > 0 or (x, y) = (0, 0)}, then (G, P ) is an Archimedean partially ordered group but not an -group.

• If P = {(x, y) : x ≥ 0 and y ≥ 0}, then (G, P ) is an Archimedean -group, and |(x, y)| = (|x|, |y|).

Let Z[ε]

= {z = x + εy : x, y ∈ Z, ε / ∈ R; ε 2 = 0} be the ring of dual Gaussian integers. Let z ∈ Z[ε] belongs to the positive cone P of Z[ε] if Re(z) > 0 or Re(z) = 0 and Im(z) ≥ 0. Then (Z[ε], P ) is a totally f -ring but not Archimedean, and the absolute value is given by.

|z| = max{z, -z} = z if z ≥ 0 -z if z ≤ 0 4. The ring Z[i]
of Gaussian integers cannot be made into an f -ring since i 2 = -1 (the squares in any f -ring are positive). Nevertheless, there is a partial order on Z[i] that makes it into an Archimedean -group; namely, z ≤ w in Z[i] if and only if Re(z) ≤ Re(w) and Im(z) ≤ Im(w). Thus, for any z ∈ Z[i], the absolute value is

|z| = |Re(z)| + i|Im(z)|.
It is worth noticing that both |z| and the modulus √ z z are generalizations of the usual absolute value in the sens that they coincide on real numbers. However, the first one belongs to Z[i] as an -group, and the second is a positive real number which represents the euclidean distance from z to the origin 0 as a lattice point.

Hyperbolic numbers

We recall basic properties of hyperbolic numbers equipped by its natural Archimedean f -algebra structure (see [START_REF] Gargoubi | f -Algebra Srtucture On Hyperbolic Numbers[END_REF]). Hyperbolic numbers defined by (1.1) is a commutative ring with group of units defined by

D * = z ∈ D : z h = 0 ,
where z h := z z = x 2 -y 2 denotes the hyperbolic square-norm of z = x + jy, (x, y ∈ R) and z is the conjugate of z given by exchanging y ←→ -y. The hyperbolic plane has an important basis {e 1 , e 2 } where

e 1 = 1 + j 2 , e 2 = 1 -j 2 .
e 1 and e 2 are mutually complementary idempotent zero divisors, i.e., e 2 1 = e 1 ; e 2 2 = e 2 ; e 1 + e 2 = 1; e 1 e 2 = 0.

(2.1)

In this basis, each hyperbolic number z can be written as

z = π 1 (z)e 1 + π 2 (z)e 2 , (2.2)
where the maps π 1 , π 2 : D → R are a pair of surjective ring homomorphisms defined by

π 1 (x + jy) = x + y and π 2 (x + jy) = x -y.
From the representation (2.2), called the spectral decomposition [START_REF] Sobczyk | The hyperbolic number plane[END_REF], algebraic operations correspond to coordinate-wise operations, the square norm of z is the product π 1 (z)π 2 (z) and its conjugation is given by exchanging π 1 (z) ↔ π 2 (z). Moreover, we can define a partial order ≤ on D that makes it into Archimedean f -algebra, where

z, w ∈ D; z ≤ w if and only if π k (z) ≤ π k (w), k = 1, 2. (2.3)
From this ordering the lattice operation are

z ∨ w = max{π 1 (z), π 2 (z)}e 1 + max{π 1 (z), π 2 (z)}e 2 , z ∧ w = min{π 1 (z), π 2 (z)}e 1 + min{π 1 (z), π 2 (z)}e 2 .
Moreover, z ∨ w and z ∧ w can be expressed as a I(D)-combination of z and w, where I(D) means the set of all idempotent elements of D. More precisely, this property is formulated in the following result.

Proposition 2.1 (Proposition 3.1 in [START_REF] Gargoubi | f -Algebra Srtucture On Hyperbolic Numbers[END_REF]). For any z, w ∈ D there exist unique u, v ∈ I(D)

satisfying uv = 0 and u + v = 1 such that z ∨ w = uz + vw and z ∧ w = vz + uw.
The Riesz space D is Dedekind complete i.e. every nonempty subset A that is bounded above (resp. below), has a supremum sup A (resp. infimum inf A), and

sup A = sup π 1 (A)e 1 + sup π 2 (A)e 2 , (2.4
)

inf A = inf π 1 (A)e 1 + inf π 2 (A)e 2 .
(2.5)

In the ring of hyperbolic numbers there is a multiplicative group S called group of signs given by

S = 1, -1, j, -j ∼ = Z/2Z × Z/2Z. (2.6)
Theorem 2.1 (Theorem 5.1 in [START_REF] Gargoubi | f -Algebra Srtucture On Hyperbolic Numbers[END_REF]). Let z ∈ D, then there exists an element ε ∈ S such that εz ≥ 0.

If z h = 0 then ε is unique, called sign of z, denoted sgn(z) and given by

sgn(z) = z |z| .
The f -algebra D under the norm

z R := min α ∈ R + : α ≥ |z| = |z| ∨ |z| for all z ∈ D,
is a unital Banach lattice algebra, i.e. the norm . R satisfies the properties:

(i) |z| ≤ |w| implies z R ≤ w R ; (ii) zw R ≤ z R w R and 1 R = 1. As |z| ∨ |z| = max{|π 1 (z)|, |π 2 (z)
|}, then using the standard basis {1, j} an explicit expression of z R is given by the formula

x + jy R = max{|x + y|, |x -y|} for all x, y ∈ R.
The Banach algebra structure allows us to define the exponential of z, for any hyperbolic number z, as

e z := ∞ n=0 z n n! = e π 1 (z) e 1 + e π 2 (z) e 2 .
The hyperbolic exponential function exp is a group homomorphism from the additive group D to the multiplicative group D * . But unlike complex numbers, the hyperbolic exponential function is one-to-one with exp

(D) = D + * = {z = z 1 e 1 + z 2 e 2 : z 1 , z 2 ∈ R * + }. Therefore, exp : D -→ D +
* is a group isomorphism, and this leads to define the hyperbolic logarithm function as the inverse isomorphism ln = exp -1 . Given z ∈ D + * and α ∈ D. Define the hyperbolic exponentiation z α as z α := e α ln(z) = e π 1 (α) ln(π 1 (z)) e 1 + e π 2 (α) ln(π 2 (z)) e 2 .

Write z = z 1 e 1 + z 2 e 2 , α = α 1 e 1 + α 2 e 2 . Then from the above formula we obtain

z α := z α 1 1 e 1 + z α 2 2 e 2 .
(2.7)

Finally, let us mention that we will use the following notation: for any z, w ∈ D we write

z ≺ w if and only if w -z ∈ D + * . Therefore, if z, w ∈ R then z < w in R if and only if z ≺ w in D.
3 Hyperbolic integers

Basic definitions and properties

Let R be a ring extension of Q with degree n, i.e., a unital commutative ring in which its underling group is a Q-vector space with dimension n. Therefore, each α ∈ R is an algebraic number ( i.e., a root of a polynomial P ∈ Z[X]) since it satisfies an equation

a n α n + a n-1 α n-1 + • • • + a 0 = 0, a i ∈ Z.
If in addition a n = 1, α is said to be an algebraic integer. The set O R of all algebraic integers of R is a ring [21, Chap 2 ] called the ring of integers of R. In particular, a finite field extension K of Q is usually refereed as a number field. The ring O K of its integers is Notherian but in general is not a unique factorization domain see [START_REF] Ireland | A Classical Introduction to Modern Number Theory[END_REF]Chap 12]. In lattice-ordered rings framework we obtain the following result. Theorem 3.1. Let R be a ring extension of Q with degree = 1. If R is an integral domain then its ring of integers O R cannot be made into Archimedean f -ring.

Proof. Let R be an integral ring extension of Q with degree d > 1. Suppose that O R can be made into Archimedean f -ring. Let x ∈ O R then x + = 0 or x -= 0 since x + x -= 0 and O R is an integral domain. This means that x ≥ 0 or x ≤ 0, since x can be written as x = x + -x -with x + , x -≥ 0. It follows that for any u, v ∈ O R , either u -v ≥ 0 or u -v ≤ 0, and hence u ≥ v or u ≤ v; that is, O R is totally ordered group. Therefore, from Hölder theorem [5, Theorem 2.6.3], O R is isomorphic to a subgroup of R. This yields a contradiction since dim Q R = 1. This completes the proof.

The aim of this paper is the characterization of all Archimedean f -rings of quadratic integers. From Theorem 3.1, these rings must be non integral domains. In fact, we show that, up to an -isomorphism, there is only one Archimedean f -rings of quadratic integers., namely the ring of integers of the extension Q(j) = {α + jβ ; α, β ∈ Q}. Theorem 3.2 (Hyperbolic integers). The ring of integers of Q(j) is given by

Z h := Ze 1 + Ze 2 ,
and it is the unique, up to order an ring isomorphism, Archimedean f -ring of quadratic integers called the ring of hyperbolic integers.

Proof. The Z-module Z h = Ze 1 + Ze 2 is a subring and sublattice ( closed under ∨ and ∧) of the Archimedean f -algebra D. Therefore it's an Archimedean f -ring under the partial order induced from D, and for any u = ne 1 + me 2 and v = pe 1 + qe 2 we have

u ∨ v = max{n, p}e 1 + max{m, q}e 2 and u ∧ v = min{n, p}e 1 + min{m, q}e 2 . (3.1)
We prove now that Z h is the ring of integers of Q(j). To do this, we will use the decomposition Q(j) = Qe 1 + Qe 2 that follows from the identities: e 1 + e 2 = 1 and e 1 -e 2 = j.

Let υ = αe 1 + βe 2 ∈ Qe 1 + Qe 2 then υ ∈ O Q(j) if and only if there exists (a, b) ∈ Z 2 such that υ 2 + aυ + b = 0, i.
α = -a + √ a 2 -4b 2 = n 2 and β = -a - √ a 2 -4b 2 = m 2 , (3.2)
where n, m, α + β and αβ are integers, so that n+m 2 , nm 4 ∈ Z. This holds only if n and m are even, i.e., from (3.2) α, β ∈ Z.

Hence Z h = O Q(j) . Let O R be the ring of integers of a quadratic ring extension R of Q. Suppose that O R is an Archimeadean f -ring. Then it contains an element b with u 1 = b + = 0 and u 2 = b -= 0.
Otherwise, v + = 0 or v -= 0 for every v ∈ O R , and so it's a totally ordered Archimedean group which implies (by Hölder theorem [5, Theorem 2.6.3]), that O R is isomorphic to a subgroup of R, which is a contradiction. As O R is an f -ring we must have u 1 u 2 = b + b -= 0 which implies that {u 1 , u 2 } are linearly independent over Q. Indeed, αu 1 +βu 2 = 0 implies αu 2 1 = 0 and βu 2 2 = 0 and so α = β = 0 since any unital Archimedean f -ring is semiprime (i.e., 0 is the only nilpotent element). Let

α 1 , α 2 ∈ Q be such that 1 = α 1 u 1 + α 2 u 2 . Thus (v 1 , v 2 ) := (α 1 u 1 , α 2 u 2
) is a basis of R and satisfying the properties

1 = v 1 + v 2 , v 1 v 2 = 0 =⇒ v 2 1 = v 1 , v 2 2 = v 2 . (3.3) Let now υ = αv 1 + βv 2 ∈ R = Qv 1 + Qv 2 then from (3.
3), a similar reasoning to that of

Q(j) = Qe 1 + Qe 2 shows that υ ∈ O R if and only α, β ∈ Z, that is, O R = Zv 1 + Zv 2 .
We claim that it is an Archimedean f -ring under the partial order (nv 1 + mv 2 ) ≤ (pv 1 + qv 2 ) if and only if n ≤ p and m ≤ q in Z. It is clear that the positive cone is closed under multiplication. Also, O R is an -group and for any two elements

u = nv 1 + mv 2 and v = pv 1 + qv 2 , u ∨ v = max{n, p}v 1 + max{m, q}v 2 and u ∧ v = min{n, p}v 1 + min{m, q}v 2 . (3.4)
The Archimedean property follows from that of Z with the usual order. For the f -ring property, we will prove that for any positive elements u, v and w we have wu

∧ v = 0 whenever u ∧ v = 0. Write u = n 1 v 1 + n 2 v 2 , v = m 1 v 1 + m 2 v 2 and w = p 1 v 1 + p 2 v 2 with n 1 , n 2 , m 1 , m 2 , p 1 , p 2 ∈ Z + . Let c 1 = min{p 1 , p 2 ,
ϕ : nv 1 + mv 2 → ne 1 + me 2 .
Clearly ϕ is bijective and ( by (3.1), (3.4)), it preserves ∨ and ∧. Moreover it follows from (2.1) and (3.3) that ϕ is also a ring homomorphism. Therefore ϕ is an -isomorphism between the two f -rings.

D = e 1 Z R ⊕ e 2 Z R
. Indeed, an hyperbolic number α is an algebraic integer if and only if α 1 = π 1 (α) and α 2 = π 2 (α) are real algebraic integers. From this point of view the ring of hyperbolic integers Z h can be seen as the smallest (with respect to inclusion) sub f -ring of Z D containing {e 1 , e 2 }.

As for integers of a quadratic field Q( √ d), every hyperbolic integer a is the root of a monic polynomial P ∈ Z[X] given by P (X) = X 2 -2Re(a)X + a h , where a h = aā and Re(a) is the real part of a. However, Z h has zero divisors which are the multiples ne of e ∈ {e 1 , e 2 } with n ∈ Z \ {0}. For the units of Z h we have.

Proposition 3.1. The units of Z h coincides with the group of signs of D (2.6); that is

U (Z h ) = S = 1, -1, j, -j .
Proof. The units of Z h are characterized by all υ ∈ Z h such that υ h = ±1 since the square norm . h is multiplicative and υ satisfies: In the hyperbolic plane D ≡ R 2 , Z h is a 'square' full lattice [START_REF] Martinet | Perfect Lattices in Euclidean Spaces[END_REF] with the fundamental parallelepiped P = {z ∈ D : 0 ≤ z ≺ 1} and minimal elements ±e 1 , ±e 2 . (ii) If A is bounded from below and closed under ∧, then it has a smallest element.

Proof. (i) If

A is finite, it is clear that max A = a∈A a. Otherwise, A is a countable set which means that A can be viewed as a sequence (a n ). Let

z n = a 0 ∨ • • • ∨ a n for n = 0, 1, • • • . Then, (z n ) n≥0
is an increasing sequence of A which is bounded above. It follows that (z n ) become constant, i.e., there exists an integer N ∈ N such that z n = z N for all n ≥ N . Hence max A = z N .

(ii) We apply (i) for -A and use the duality formula inf A = -sup(-A).

Ideals of Z h

In this subsection, we establish some properties involving ideals of Z h .

Proposition 3.3. For every ideal I in the ring Z h there exists a unique positive element g I such that I = g I Z h . Moreover, I is a sublattice of Z h .

Proof. Let I be an ideal of the ring Z h . Since for every k = 1, 2 the map π k is a surjective ring homomorphism from Z h to Z, then π k (I) is an ideal of the principal ideal domain Z. Therefore, there is a unique positive integer n k such that π k (I) = n k Z. Thus, the element g I = n 1 e 1 + n 2 e 2 generates I and it is the only positive one. It follows from Proposition 2.1 that I is a sublattice of Z h .

Recall that an -subgroup C (i.e., subgroup and sublattice) of an -group G is said to be convex if 0 ≤ a ≤ b in G and b ∈ C imply a ∈ C. An -ideal of an -ring R is a convex -subgroup of R that is also an ideal of R. The following characterizes -ideals of the -ring Z h . Proposition 3.4. An ideal of Z h is an -ideal if and only if it is generated by an idempotent element.

Proof. Let I be an ideal of Z h with positive generator g I . Proposition 3.3 means that I is an -subgroup of Z h . So, I is an -ideal if and only if I is convex. Suppose that g I is an idempotent element, i.e., g I ∈ {0, 1, e 1 , e 2 }. It is obvious that I is convex if g I = 0 or g I = 1. Assume that g I = e ∈ {e 1 , e 2 } that means I = eZ h = eZ. Let a, b ∈ Z h be such that 0 ≤ a ≤ b and b ∈ I. Then a = αe for some real α, since eR is an order ideal of D (see Theorem 3.5 in [START_REF] Gargoubi | f -Algebra Srtucture On Hyperbolic Numbers[END_REF]). We must also have α ∈ Z, because a ∈ Z h . Hence, a ∈ I and this proves that I is convex. Conversely, assume that I is convex. We have to prove that the generator g I of I is an idempotent element. We distinguish two cases: (i) if g I h = 0 the case g I = 0 is trivial. Suppose that g I = ne for some e ∈ {e 1 , e 2 } and some integer n ≥ 1. Then, we have 0 ≤ e ≤ ne with ne ∈ I which implies that e ∈ I = neZ. Thus, nk = 1 for some k ∈ Z. This yields that n = 1, and hence g I = e.

(ii) if g I h = 0, then 0 ≤ 1 ≤ g I with g I ∈ I which means that 1 ∈ I, i.e., I = Z h . Therefore, g I = 1.

Proposition 3.5. Let I be an ideal of Z h with the positive generator g I = n 1 e 1 + n 2 e 2 . Then,

Z h /I Z/n 1 Z × Z/n 2 Z.
In particular, Z h /I Z if and only if I is a nontrivial -ideal.

Proof. Let I be an ideal in Z h with the positive generator g I = n 1 e 1 + n 2 e 2 . One can easily see that the mapping

Z h /I ȧ → π 1 (a), π 2 (a) ∈ Z/n 1 Z × Z/n 2 Z, establishes an isomorphism of Z h /I with Z/n 1 Z × Z/n 2 Z.
In particular, from Proposition 3.4, Z h /I Z if and only if I is a nontrivial -ideal. By (2.4) and (2.5), we derive that

Hyperbolic floor and ceiling functions

z D = π 1 (z) e 1 + π 2 (z) e 2 and z D = π 1 (z) e 1 + π 2 (z) e 2 . (3.5) Therefore z -1 ≺ z D ≤ z ≤ z D ≺ z + 1 for all z ∈ D. (3.6)
Let as consider a hyperbolic closed interval (see [START_REF] Balankin | Cantor-type sets in hyperbolic numbers[END_REF] and [START_REF] Téllez-Sánchez | More about Contor like sets in hyperbolic numbers[END_REF] ) defined by As on real numbers, the functions . D and . D allows one to determine the number N Z h (I) of hyperbolic integers in a hyperbolic interval I by considering the four types above.

Proposition 3.6. Let α, β ∈ D be such that α ≤ β then

N Z h ([α, β] D ) = β D -α D + 1 h , N Z h ([α, β) D ) = β D -α D h , N Z h ((α, β] D ) = β D -α D h , N Z h ((α, β) D ) = β D -α D -1 h .
Proof. Let us denote by

I 1 = [α, β] D , I 2 = [α, β) D , I 3 = (α, β] D and I 4 = (α, β) D . Then the sets Z h ∩ I k are bijectively mapped onto Z 2 ∩ ϕ(I k ) via the map from D to R 2 defined by ϕ(z) = (π 1 (z), π 2 (z)). Thus, for k = 1, . . . , 4, N Z h (I k ) = #Z 2 ∩ R k where R k = ϕ(I k ) are the rectangles R 1 = [π 1 (α), π 1 (β)] × [π 2 (α), π 2 (β)] , R 2 = [π 1 (α), π 1 (β)) × [π 2 (α), π 2 (β)) R 3 = (π 1 (α), π 1 (β)] × (π 2 (α), π 2 (β)] , R 4 = (π 1 (α), π 1 (β)) × (π 2 (α), π 2 (β)) .
Therefore,

N Z h (I 1 ) = ( π 1 (β) -π 1 (α) + 1) ( π 2 (β) -π 2 (α) + 1) , N Z h (I 2 ) = ( π 1 (β) -π 1 (α) ) ( π 2 (β) -π 2 (α) ) , N Z h (I 3 ) = ( π 1 (β) -π 1 (α) ) ( π 2 (β) -π 2 (α) ) , N Z h (I 4 ) = ( π 1 (β) -π 1 (α) -1) ( π 2 (β) -π 2 (α) -1) .
Finally, the results follows from (3.5) and the propriety z h = π 1 (z)π 2 (z).

Divisibility

First properties

Divisibility in Z h is defined in a natural way: we say b divides a , or a is a multiple of b (and write b|a) if a = bc for some c ∈ Z h . In this case, we call b a divisor of a. (iii) a|b and b|a if and only if a = εb for some unit ε, i.e., (by Proposition 5.1 in [START_REF] Gargoubi | f -Algebra Srtucture On Hyperbolic Numbers[END_REF]) if and only if |a| = |b|.

Hyperbolic euclidean division and congruence

Theorem 4.1. Let a, b ∈ Z h with b h = 0, then there exist unique q, r ∈ Z h such that

a = bq + r, 0 ≤ r ≺ |b|.
The hyperbolic integers q and r are called, respectively, the quotient and the remainder of the division of a by b.

Proof. We consider first the uniqueness. Assume that

a = bq 1 + r 1 = bq 2 + r 2 , 0 ≤ r 1 , r 2 ≺ |b|. Then, 0 ≤ |q 1 -q 2 | = |r 1 -r 2 | |b| ≺ 1.
Hence, q 1 = q 2 which implies r 1 = r 2 . Consider now the existence. Put As for integers, congruences in Z h are defined using divisibility. Since congruence modulo 0 means equality and a|b if and only if |a| b, we usually assume the modulus is a nonzero positive element of Z h .

(ii) if also a, b, v ∈ Z then a ≡ b mod v in Z h if and only if a ≡ b mod v in Z; (iii) a ≡ b mod v in Z h if and only if π k (a) ≡ π k (b) mod π k (v) in Z for k = 1, 2; (iv) a ≡ b mod v if and only if ā ≡ b mod v. Proof. Straightforward. Let v ∈ Z h , v
0. Then (by proposition 3.5), the number of the residue classes modulo v is its square norm v h . For instance, the four binary classes are the set

Z h /2Z h = { 0, 1, e 1 , e 2 }.
(4.2)

Positive gcd and positive lcm

According to Proposition 3.3, for every a, b ∈ Z h the ideals aZ h + bZ h and aZ h ∩ bZ h are generated by a unique positive element. This justifies the following result. As an immediate consequence of Theorem 4.2 we have the following properties of the gcd Z h and lcm Z h which are extension of the corresponding ones in Z. Proposition 4.3. For a, b ∈ Z h we have

(i) gcd Z h (|a|, |b|) = gcd Z h (a, b) and lcm Z h (|a|, |b|) = lcm Z h (a, b); (ii) gcd Z h (a, b) lcm Z h (a, b) = |ab|; (iii) gcd Z h (a, b) = gcd Z h (ā, b) and lcm Z h (a, b) = lcm Z h (ā, b); (iv) gcd Z h (a ∨ b, a ∧ b) = gcd Z h (a, b) and lcm Z h (a ∨ b, a ∧ b) = lcm Z h (a, b),

Primes and irreducibles in Z h

In this section we characterize all prime and irreducible elements of Z h . We also extend the unique factorization theorem of Z to Z h . The set of all prime numbers: 2, 3, 5, 7, 11, • • • will be denoted by P. For basic notions and terminology about prime and irreducible elements we refer to [16, Chap II].

Characterization

Theorem 5.1 (Hyperbolic primes). The following statements are satisfied.

(i) Prime elements of Z h are u = εv where ε is a unit and v ∈ {e 1 , e 2 , pe 1 + e 2 , e 1 + pe 2 : p ∈ P}.

(ii) Irreducible elements of Z h are u = εv where ε is a unit and v is in (ii) It is well known that Gaussian primes are, up to units, prime numbers p with p ≡ 3 mod 4 or Gaussian integers z = a + ib with norm z z = a 2 + b 2 is a prime number. Similarly, (non-zero-divisor) hyperbolic primes are, up to units, hyperbolic integers υ with square norm υ h = υ ῡ = p, where p is a prime number. However, for each prime number p one has the decomposition p = p e 1 p e 2 . It is quite remarkable to see that p viewed as an hyperbolic integer is in fact "semiprime".

P h := {pe 1 + e 2 ,

Unique factorization theorem

The fundamental theorem of arithmetic states that every nonzero integers n can be written uniquely in the form

n = ε p∈P p vn(p) ,
where ε is a unit ( ε = sgn(n)) and v n : P -→ N with v n (p) = 0 for a finite number of p.

Using the hyperbolic exponentiation (2.7), the following statements show that this property can be generalized to hyperbolic integers.

Theorem 5.2. Every a ∈ Z h with a h = 0 can be written uniquely in the form

a = ε p∈P p va(p) ,
where ε is a unit and v a : P -→ Z + h with v a (p) = 0 for a finite number of p. Proof. Let a ∈ Z h with a h = 0. By Theorem 2.1, there is a unique unit ε ∈ S such that a = ε|a|. Let n 1 , n 2 ∈ N be such that |a| = n 1 e 1 +n 2 e 2 . Then, n 1 , n 2 = 0. By fundamental theorem of arithmetic, for every n ∈ N there is a unique application Theorem 5.3 (Unique factorization theorem). Every a ∈ Z h with a h = 0 can be written uniquely in the form a = ε

u∈P h u va(u) ,
where ε is a unit and v a : P h -→ N with v a (u) = 0 for almost all u.

Proof. The proof follows immediately from Theorem 5.2 by observing that hyperbolic primes are p e 1 and p e 2 with p ∈ P ( Theorem 5.1).

Hyperbolic Gaussian integers

By analogy to complex numbers, the hyperbolic Gaussian integers or more simply the h-Gaussian integers (also called split Gaussian integers [START_REF] Cifuente | The Ring of Integers in the Canonical Structures of the Planes[END_REF]) are the set

G h = Z[j] := x + jy : x, y ∈ Z .
It is a subring of Z h with zero divisors that are the multiples n(1 ± j), n ∈ Z \ {0}, and units that are 1, -1, j, and -j. From the four binary classes (4.2) of hyperbolic integers we have the following characterization of h-Gaussian integers. Observing that [-ρ, ρ] D is the closed ball B R (0, ρ) in (D, . R ) where . R is the lattice norm (2.2), then Eq (7.2) yields that there exists a hyperbolic integer a = ne 1 + me 2 0 such that a R = max{n, m} > ρ and a h = nm ≤ ρ. Therefore, from the identity nm = max{n, m} min{n, m} one has that ρ ≥ nm > ρ. Which is a contraction. Hence,

Z h ∩ D + (ρ) = Z h ∩ D + (ρ).
Let n be an integer Proof. From proposition 7.3 ∆(ρ) is given by ∆ -(ρ), and the error is the order of δ(ρ). Thus, the proof follows from the definition of δ by observing that d(n) = O(n ) for every > 0.

  e., from (2.1) (α 2 + aα + b)e 1 + (β 2 + aβ + b)e 2 = 0. This means that α and β are the roots of x 2 + ax + b. Then up to a permutation of the roots one has

1} and c 2 = max{p 1 , p 2 ,

 212 1}. Thus from (3.4), we get c 1 (u ∧ v) ≤ wu ∧ v ≤ c 2 (u ∧ v), and hence u ∧ v = 0 implies wu ∧ v = 0. It remains to prove that the f -rings O R and Z h are -isomorphic. Define the mapping

Remark 3 . 1 .

 31 Using notations of [21, Chap 2 ] one can define the set Z D as the integral closure of Z in D i.e. the ring of algebraic integers of D. From Proposition 2.1, It is clear that Z D is a sublattice of D, under the induced partial order (2.3), and then an Archimedean f -ring. One can easy check that Z

υ 2 +

 2 av + b = 0, where a, b ∈ Z with a = 2Re(υ) and b = υ h . Write υ = ne 1 + me 2 , then υ h = nm = ±1 if and only if (n, m) ∈ {±(1, 1), ±(1, -1)}. As e 1 + e 2 = 1 and e 1 -e 2 = j, we get υ = ±1, ±j.

Figure 1 :

 1 Figure 1: Hyperbolic integers with fundamental parallelepiped P .

  (i) If A is bounded from above and closed under ∨, then it has a largest element.

Definition 3 . 1 .

 31 Let z ∈ D. Then the sets E + (z) := {k ∈ Z h : k ≤ z} and E -(z) := {k ∈ Z h : k ≥ z} are two nonempty sublattices of Z h . Thus from Proposition 3.2, the notions of floor . and ceiling . functions on real numbers can be extended to the hyperbolic numbers in the following way. The functions . D and . D from D to Z h defined byz D := max k ∈ Z h : k ≤ z , z D := min k ∈ Z h : k ≥ zare called, respectively, hyperbolic floor function and hyperbolic ceiling function.

Figure 2 :

 2 Figure 2: Non degenerate hyperbolic closed interval [α, β] D . similar way replacing ≤ by ≺ in left-right and left/right, respectively. However, all these intervals are empty if α -β h = 0. One has (α, β) D = [α, β] D \ (∂ α ∪ ∂ β ), (α, β] D =

Proposition 4 . 1 .

 41 For a, b ∈ Z h we have (i) a|b in Z h implies a h b h ∈ Z; (ii) a|b and b h = 0 implies |a| ≤ |b|; (iii) a|b and b|a if and only if |a| = |b|. Proof. (i) If a|b then b = ac which implies, by multiplicity of . h , that b h = a h c h . Hence, a h b h ∈ Z. (ii) If a|b and b h = 0, then b = ca for some c ∈ Z h with c h = 0 that means |c| ≥ 1. Therefore, |b| -|a| = (|c| -1)|a| ≥ 0.

D

  and r = a -bq, where ε = sgn(b) = |b| b and . D is the hyperbolic floor function (3.1). Then, we have q, r ∈ Z h and a = bq + r. It remains to prove that 0 ≤ r ≺ |b|. From (3.6) one has by -εb = -|b| and use r = a -bq to get the desired inequality.

Definition 4 . 1 .

 41 Let a, b, v ∈ Z h , we write a ≡ b mod v if and only if v|(a -b).

Proposition 4 . 2 .

 42 For a, b, c, v ∈ Z h one has (i) a ≡ b mod v and c ≡ d mod v implies a + c ≡ b + d mod v and ac ≡ bd mod v;

Theorem 4 . 2 .

 42 Every a, b ∈ Z h have a unique positive greatest common divisor gcd Z h (a, b) and a unique positive latest common multiple lcm Z h (a, b). Moreover, gcd Z h (a, b) = gcd(π 1 (a), π 1 (b))e 1 + gcd(π 2 (a), π 2 (b))e 2 , lcm Z h (a, b) = lcm(π 1 (a), π 1 (b))e 1 + lcm(π 2 (a), π 2 (b))e 2 .

Remark 4 . 1 .

 41 In view of Proposition 4.1, the quazi-order on Z + h defined by a ≤ b if and only if a|b, is a partial order, and under such order, Z + h is a lattice ordered multiplicative monoid with a ∧ b = gcd Z h (a, b) and a ∨ b = lcm Z h (a, b) for all a, b ∈ Z + h .

  µ n : P -→ N with µ a (n) = 0 for almost all p such that n = p∈P p µn(p) . Let v a : P -→ Z + h be the function defined by v a (p) = µ n 1 (p)e 1 + µ n 2 (p)e 2 . Therefore |a| = p∈P p va(p) .

DefineD

  (ρ) = D(ρ) ∩ [-ρ, ρ] D . Geometrically, D (ρ) is the square [-ρ, ρ] D if ρ ≤ 1, and D (ρ) [-ρ, ρ] D if ρ > 1 as represented in Figure 3.

Figure 3 : 1 2 j ln ρ . Proposition 7 . 1 .

 3171 Figure 3: A representation of D (ρ) for ρ > 1 with τ = √ ρe 1 2 j ln ρ .

≥ 2 .

 2 Define ξ k , λ k , µ k and η k such thatξ k = √ ρe jk 2n ln ρ , λ k = ξ k ∧ √ ρ, for k = 0, • • • , n; η k = ξ k ∨ ξ k-1 , µ k = ξ k ∧ ξ k-1 , for k = 1, • • • , n.Thus, and referring to Figure 4, proposition 7.1 yields that D - n (ρ) ≤ D(ρ) ≤ D + n (ρ), (7.3)

Figure 4

 4 Figure 4: .

Proposition 7 . 2 . 1 2 e 1 , λ k = ρ 1 2 e 1 + ρ n-k 2n e 2 , ξ k = ρ n+k 2n e 1 + ρ n-k 2n e 2 , η k = ρ n+k 2n e 1 N 2 ,N= ρ 1 2 2 + 2 ρ -ρ 1 2 -.Theorem 7 . 1 .

 72111221212271 where D - n (ρ) = N Z h (0, √ ρ] D + 2N Z h ((α, τ ] D ) + 2 k=n k=2 N Z h ((λ k , ξ k-1 ] D ) , (7.4)D + n (ρ) = N Z h (0, √ ρ] D + 2N Z h ((α, τ ] D ) + 2 k=n k=1 N Z h ((λ k , η k ] D ) .(7.5)Let ∆ - n (ρ) and ∆ + n (ρ) be such that∆ - n (ρ) = D - n (ρ) -(ρ ln ρ + (2γ -1)ρ), ∆ + n (ρ) = D + n (ρ) -(ρ ln ρ + (2γ -1)ρ). Let δ(ρ) be a function defined for ρ > 1 by δ(ρ) = 0 if ρ / ∈ N and δ(ρ) = d(ρ) + χ N ( √ ρ),otherwise. One has lim supn→∞ (∆ + n (ρ) -∆ - n (ρ)) ≤ δ(ρ).Proof. We have∆ + n (ρ) -∆ - n (ρ) = D + n (ρ) -D - n (ρ) = 2 k=n k=1 N Z h ((µ k , η k ] D ) .Denote by J = Z h ∩(D (ρ)\γ * ρ ) where γ * ρ is the image of γ ρ (t) = √ ρe jt defined in 0, 1 2 ln ρ . Since J is a nonempty finite set one obtainsd γ * ρ = min h∈J d(h, γ * ρ ) = min h∈J inf ξ∈γ * ρ h -ξ > 0.Let N be an integer such that2 √ ρ sinh( ln ρ 4N ) cosh(ln ρ) < d γ * ρ .From the above inequality we have∆ - n (ρ) = ψ(ρ) + 2 k=n k=2 N Z h ((λ k , ξ k-1 ] D ) , ψ(ρ) = N Z h (0, √ ρ] D + 2N Z h ((α, τ ] D ) -(ρ ln ρ -(2γ -1)ρ).Straightforward calculations giveτ = ρe 1 + e 2 , α = ρ Z h ((α, τ ] D ) = ρ -ρ 1 Z h ((λ k , ξ k-1 ] (ρ ln ρ -(2γ -1)ρ).Finally ∆ -(ρ) = ψ(ρ) + 2 lim inf For every real > 0 we have ∆(ρ) = ψ(ρ) + 2 lim inf

  e 1 + pe 2 : p ∈ P}.where P h is defined as the set of hyperbolic prime numbers (or hyperbolic primes) , Proof. Let v = ne 1 + me 2 be a nonzero and nonunit positive element of Z h . (i) By Proposition 3.5, v is prime if and only if Z/n Z × Z/m Z is an integral domain, i.e, if and only if one of Z/n Z or Z/m Z is an integral domain and the other one is zero. It follows that (n, m) ∈ {(1, 0), (0, 1), (p, 1), (1, p) : p ∈ P}. Which means that v ∈ {e 1 , e 2 , pe 1 + e 2 , e 1 + pe 2 : p ∈ P}.(ii) If v is irreducible, then vZ h is a maximal ideal, since by Proposition 3.3, every ideal in the ring Z h is principal. Therefore, vZ h is a prime ideal which means that v is prime. So from (i), either v = e 1 or v = e 2 or v = pe 1 + e 2 or v = e 1 + pe 2 for some prime number p. We will prove that each non-zero-divisor prime element v is irreducible, since the atoms e 1 and e 2 are not. Let a, b ∈ Z h be such that v = ab. Then, taking the norm . h , we obtain p = a h b h . Since p is irreducible in Z it follows that a h = ±1 or b h = ±1. Hence, either a ∈ S or b ∈ S. Remark 5.1. (i) Theorem (5.1) shows that hyperbolic primes are positive non-zerodivisor prime elements of Z h and they are the hyperbolic integers of the form pe 1 +e 2 or e 1 +pe 2 with p ∈ P. Write p = pe 1 +pe 2 , then using the hyperbolic exponentiation (2.7) we have pe 1 + e 2 = p e 1 and e 1 + pe 2 = p e 2 . Thus hyperbolic primes are the set P h := {p e : (p, e) ∈ P × {e 1 , e 2 }}.

Acknowledgment. We are grateful to Mounir Bezzarga for enlightening discussions and fruitful comments.

Theorem 6.1. Let a ∈ Z h , then a ∈ G h if and only if either a ≡ 0 mod 2 or a ≡ 1 mod 2.

Proof. Let a = ne 1 + me 2 = n+m 2 + j n-m 2 ∈ Z h with n, m ∈ Z. So, a ∈ G h if and only if n ≡ m mod 2, i.e., if and only if either a ≡ 0 mod 2 or a ≡ 1 mod 2.

In view of Theorem 2.1 and by the fact that the units of G h are the set S one can see that G h is closed under absolute value. But, it is not an -subgroup of Z h , since 0 ∨ j = e 1 / ∈ G h . However, the next result gives under which condition the supremum of two incomparable (with respect to the order induced by

Proof. Let a, b ∈ G h be incomparable. Then, from Proposition 2.1 we can write

Proposition 6.2. Every a ∈ G h with a h = 0 can be written uniquely in the form

where ε a is a unit, v a : P -→ Z + h with v a (p) = 0 for almost all p and ν ∈ Z + h is such that ν = 0 if a ≡ 1 mod 2 and ν 0 if a ≡ 0 mod 2.

Proof. From Theorem 5.2, a can be uniquely expressed in the form

where ε a is a unit and P va -→ Z + h with v a (p) = 0 for almost all p and ν = v a (2) ∈ Z + h . Therefore, a ≡ 2 ν mod 2 since, p va(p) ≡ 1 mod 2 for p = 2 . It follows that ν = 0 if a ≡ 1 mod 2 and ν 0 if a ≡ 0 mod 2.

Dirichlet divisor problem

The Dirichlet divisor problem, arises from estimating D(ρ) := n≤ρ d(n), where d(n) is the number of positive divisors of n. A well known result is D(ρ) = ρ ln ρ+(2γ-1)ρ+∆(ρ), where γ is Euler's constant and ∆(ρ) is the error term. The Dirichlet divisor problem asks for the correct order of magnitude of ∆(ρ) as ρ -→ ∞ (see e.g. [START_REF] Krätzel | Lattice points[END_REF]Chap 5]). From a geometrical point of view D(ρ) is equal to the number of lattice points in the first quadrant under the hyperbola xy = ρ. Thus, this is equivalent to determine the number of hyperbolic integers a 0 with a h ≤ ρ, i.e.,

Then, for every n ≥ N and for every k

It follows from the inclusion above that

Thus, it follows from

This completes the proof.

-ρ


Alexandre Ern
email: ern@cermics.enpc.fr

Igor Mozolevski
email: igor.mozolevski@mtm.ufsc.br

Luciane Schuh
email: lucianes@ime.usp.br

Discontinuous Galerkin

Keywords: two-phase flows, heterogeneous porous media, discontinuous capillary pressure, discontinuous Galerkin, interface condition, velocity reconstruction, weighted averages, secondary oil recovery PACS: 02.70.Dh, 47.56.+r, 02.60.Cb

come

2000 MSC: 65M60, 65N30, 76T99, 76S05

Introduction

Two-phase immiscible incompressible flows through porous media are important in many applications, e.g., secondary oil recovery in petroleum engineering [START_REF] Aziz | Petroleum Resevoir Simulation[END_REF]. Swift changes of rock properties in an oil reservoir lead to capillary pressure discontinuities which can considerably reduce the recovery factor and even cause oil trapping; see, e.g., [START_REF] Bertsch | Analysis of oil trapping in porous media flow[END_REF][START_REF] Ersland | Numerical methods for flow in a porous medium with internal boundaries[END_REF][START_REF] Van Duijn | The effect of capillary forces on immiscible two-phase flow in heterogeneous porous media[END_REF]. For such applications, the development of numerical methods that can cope satisfactorily with two-phase flows in heterogeneous porous media with discontinuous capillary pressures is a relevant and challenging problem.

The physical principles governing two-phase immiscible incompressible flows in porous media and the corresponding mathematical models can be found, e.g., in [START_REF] Bear | Dynamic of Fluids in Porous Media[END_REF][START_REF] Chavent | Mathematical Models and Finite Elements for Reservoir Simulation[END_REF][START_REF] Chen | Computational Methods for Multiphase Flows in Porous Media[END_REF][START_REF] Helmig | Multiphase Flow and Transport Processes in the Subsurface[END_REF]. The governing equations generally consist of an elliptic Darcy-type equation for the pressure and a nonlinear degenerate parabolic equation with a nonlinear advective term for the saturation. The dependent variables are either the total or the wetting-phase pressure and either the wetting or the non-wetting phase saturation. The governing equations are often solved sequentially: at each time step, the pressure equation is solved first, then the total velocity is computed from the pressure gradient, and finally, the saturation equation is advanced in time using the total velocity in the advective term. This approach substantially alleviates the nonlinearities in the coupled system of governing equations. Furthermore, working with the global pressure instead of the wetting-phase pressure offers the advantage of less strongly coupled equations, since the capillary pressure no longer enters explicitly the pressure equation.

Discontinuous Galerkin methods (DG) offer an attractive tool to approximate (in space) two-phase flows in porous media. Advantages of such methods include their local (elementwise) conservation property, their flexibility in using non-matching meshes and variable polynomial degrees, and their potential to capture shocks sharply while maintaining high accuracy in smooth regions. A mathematical analysis of DG methods applied to Friedrichs' systems encompassing both elliptic and hyperbolic PDE's can be found in [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF]. Different versions of DG methods have been considered for the sequential pressure-saturation formulation of two-phase flows. In [START_REF] Bastian | Discontinuous Galerkin methods twophase flow in porous media[END_REF], the nonsymmetric interior penalty DG method was applied to the pressure equation and to the diffusive term in the saturation equation, upwinding was used for the advective term, and the total velocity was recovered using the postprocessing proposed in [START_REF] Bastian | Superconvergence and H(div) projection for discontinuous Galerkin methods[END_REF]. Similar techniques together with further variants of DG methods for the diffusive terms were investigated numerically in [START_REF] Klieber | Adaptive simulations of two-phase flow by discontinuous Galerkin methods[END_REF] in the context of adaptive methods. In [START_REF] Eslinger | Discontinuous Galerkin Finite Element Methods apllied to two-phase, air-water flow problems[END_REF], a compressible air-water two-phase problem was studied numerically using various variants of DG methods for the pressure equation and the so-called Local DG method for the saturation equation. In general, when approximating rough solutions by DG methods, spurious oscillations triggered by Gibbs-type phenomena occur and, in the context of two-phase flows, can become significant. To control them, slope limiters can be used, as advocated for instance in [START_REF] Bastian | Discontinuous Galerkin methods twophase flow in porous media[END_REF]. However, slope limiters are fairly complicated to be designed for unstructured meshes and there is no analysis available in two or three space dimensions. One possibility proposed and assessed numerically in [START_REF] Epshteyn | Fully implicit discontinuous finite element methods for two-phase flow[END_REF] to avoid spurious oscillations is to consider a fully implicit DG method for the complete coupled system of pressure and saturation equations, but this entails substantial computational costs.

All of the above DG methods for two-phase flows cannot cope satisfactorily with capillary pressure discontinuities produced at the interface between two subdomains with different rock properties. Indeed, such discontinuities require a nonlinear interface condition for the pressure and the saturation, and both variables can exhibit a nonzero jump at the interface; see [START_REF] Bear | Dynamic of Fluids in Porous Media[END_REF][START_REF] Bertsch | Analysis of oil trapping in porous media flow[END_REF][START_REF] Enchéry | Numerical approximation of a two-phase flow problem in a porous medium with discontinuous capillary forces[END_REF] for a detailed discussion. As a result, the usual paradigm in the penalty strategy of DG methods must be revisited. Formulating a suitable penalty strategy to incorporate weakly the nonlinear interface conditions on the pressure and on the saturation in the context of sequential pressure-saturation formulations constitutes the first novel contribution of this work. From a mathematical viewpoint, existence and uniqueness of the weak solution to the coupled system of pressure-saturation equations with general interface conditions was analyzed in [START_REF] Amaziane | Existence of solutions to various rock types model of two-phase flow in porous media[END_REF][START_REF] Bourgeat | Effective model of two-phase flow in a porous medium made of different rock types[END_REF] using vanishing viscosity and homogenization techniques. Moreover, in the past few years, two-phase flows with capillary pressure discontinuities have been studied extensively in the context of Finite Volume methods (FV). An interface problem for a nonlinear degenerate parabolic equation, representative of the saturation equation with discontinuous capillary pressure, was analyzed in [START_REF] Enchéry | Numerical approximation of a two-phase flow problem in a porous medium with discontinuous capillary forces[END_REF]. A FV scheme for the approximation of the solution was proposed and shown to converge to a weak solution. A more general formulation of the interface conditions was considered in [START_REF] Cancès | Nonlinear parabolic equations with spatial discontinuities[END_REF][START_REF] Cancès | Two-phase flows involving capillary barriers in heterogeneous porous media. Interfaces and Free Boundaries[END_REF], where a new definition of a weak solution for such interface conditions was introduced; existence and uniqueness of the weak solution was proven. Interface problems for degenerate parabolic equations including a nonlinear advective term with a discontinuity at the interface for the flux function were analyzed in [START_REF] Cancès | Finite volume scheme for two-phase flow in heterogeneous porous media involving capillary pressure discontinuities[END_REF]. The above FV schemes solely handle the saturation equation with a prescribed total velocity field, whereas the present work will address the coupled system of pressure-saturation equations where the total velocity is also an unknown.

The second contribution of the present work is to propose some possible improvements in the implementation of DG methods for two-phase flows in heterogeneous porous media. These improvements concern the total velocity reconstruction and the use of weighted averages in the formulation of the DG method. Firstly, we consider the total velocity reconstruction introduced in [START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF] for elliptic problems and extended to two-phase flow problems in [START_REF] Ern | Accurate velocity reconstruction for Discontinuous Galerkin approximations of two-phase porous media flows[END_REF]. This velocity reconstruction offers various advantages with respect to the reconstruction proposed in [START_REF] Bastian | Superconvergence and H(div) projection for discontinuous Galerkin methods[END_REF]. While both reconstructions enforce continuity of the normal component of the total velocity at mesh interfaces and hence produce a total velocity field in the space H(div) defined below, the present reconstruction exploits all the degrees of freedom in the DG method and, as such, produces an optimal approximation of the divergence of the total velocity. As a result, it leads to more accurate and stable calculations, as illustrated numerically in [START_REF] Ern | Accurate velocity reconstruction for Discontinuous Galerkin approximations of two-phase porous media flows[END_REF]. A further advantage is that even for piecewise affine discrete pressures, the reconstructed velocity already belongs to the first-order Raviart-Thomas(-Nédélec) finite element space defined below. Therefore, it is possible to consider the natural setting of equal-order interpolation (e.g., piecewise affine) for the pressure and the saturation. In the present work, we extend this reconstruction to the more complex situation where the pressure and the saturation can exhibit nonzero jumps at an interface owing to capillary pressure discontinuity. Secondly, following previous work on DG methods for advection-diffusion equations in heterogeneous media [START_REF] Di Pietro | Discontinuous Galerkin methods for anisotropic semidefinite diffusion with advection[END_REF][START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF], we formulate the consistency terms in the DG method using weighted averages instead of the more usual arithmetic averages, thereby allowing to use the harmonic average of the corresponding diffusivity to penalize interelement jumps. Weighted averages are used in both pressure and saturation equations, and are particularly important in the latter because of possible degeneracies.

The paper is organized as follows. In §2, we briefly restate the governing equations for two-phase flows in porous media in the coupled total pressure/fractional flow formulation and present the interface conditions induced by capillary pressure discontinuities. In §3, we present the sequential DG dis-cretization of the governing equations with weighted averages, the penalty strategy to incorporate the interface conditions, and the accurate total velocity reconstruction. Finally, §4 contains the numerical results. No limiters are used in this work.

Problem formulation

The setting

Let Ω be a bounded, open, polyhedral domain in R d , d ≥ 1, with boundary ∂Ω and outward normal n ∂Ω . The domain Ω is divided into two open, polyhedral subdomains Ω 1 and Ω 2 such that Ω = Ω 1 ∪ Ω 2 and Ω 1 ∩ Ω 2 = ∅. Let Γ := ∂Ω 1 ∩ ∂Ω 2 denote the interface between the two subdomains. More general situations with more than two subdomains can also be considered.

Each subdomain Ω β , β ∈ {1, 2}, represents a porous medium with porosity φ β and intrinsic (absolute) permeability K β ; both quantities are assumed to be constant in each subdomain for simplicity. The mobility of the wetting (α = w) and of the non-wetting phase (α = n) in the subdomain Ω β is denoted by λ αβ , the sum λ β := λ wβ +λ nβ being the total mobility and f β := λ nβ λ β the fractional flux. Let π β be the capillary pressure in the subdomain Ω β . In the sequel, the index β is used to indicate the restriction of a function defined on Ω to the subdomain Ω β .

The capillary pressure and the mobilities are in each subdomain given smooth functions of the non-wetting phase saturation. The non-wetting phase saturation takes values in the interval [S nr , 1 -S wr ], and up to rescaling, we can assume that the same interval is used in both subdomains. Here, S nr and S wr respectively denote the residual saturation of the wetting and non-wetting phase. In each subdomain Ω β , β ∈ {1, 2}, the capillary pressure π β : [S nr , 1 -S wr ) → [0, +∞) is an increasing, continuously differentiable function. We assume that ∀s ∈ [S nr , 1 -S wr ),

π 1 (s) < π 2 (s). (1)
Moreover, we assume that π 2 (S nr ) < lim s→1-Swr π 1 (s). As a result, there is a critical value S * of the non-wetting phase saturation such that

π 1 (S * ) = π 2 (S nr ). ( 2
)
This situation is illustrated in Figure 1. The value π β (S nr ) is called an entry pressure and corresponds to the capillary pressure needed to displace the non-wetting phase from the largest pore in the subdomain Ω β . The difference between entry pressures in Ω 1 and Ω 2 is at the origin of the oil trapping effect [START_REF] Bertsch | Analysis of oil trapping in porous media flow[END_REF]. Finally, regarding the mobilities, we assume that λ nβ : [S nr , 1 -S wr ] → [0, +∞) and λ wβ : [S nr , 1 -S wr ] → [0, +∞) are monotone Lipschitz functions such that λ nβ (S nr ) = 0 and λ wβ (1 -S wr ) = 0. We also assume that the fractional fluxes f β : [S nr , 1 -S wr ] → [0, 1] are increasing functions. Usually, these functions have a so-called S-shape.

Governing equations

The governing equations for two-phase immiscible incompressible flows through the heterogeneous porous medium Ω in the global pressure/fractional flow formulation and in the absence of gravity forces can be written as follows: For a given simulation time T > 0, find (p, s) that satisfy in

Ω β × [0, T] for each β ∈ {1, 2}, -∇• (λ β (s β )K β ∇p β ) = q wβ + q nβ , (3)
u β = -λ β (s β )K β ∇p β , (4)
φ β ∂ t s β + ∇• (u β f β (s β )) -∇• (ǫ β (s β )∇π β (s β )) = q nβ , (5)
where the dependent variables p and s respectively denote the global pressure (in short, the pressure) and the non-wetting phase saturation (in short, the saturation), while u denotes the total velocity. In addition, we have set ǫ β := λ wβ f β K β , and q α , α ∈ {w, n}, denotes volumetric sources or sinks of the corresponding phase in the medium. Observe that the pressure equation does not degenerate, whereas the saturation equation does in the limits s → S nr and s → 1 -S wr .

Interface conditions on Γ are discussed in the following subsection. The initial condition prescribes the saturation over the subdomains Ω β at initial time. We assume that the wetting phase is initially present in the whole domain. Boundary conditions can be of Dirichlet or Neumann type for both the pressure and the saturation. Dirichlet boundary conditions are enforced on the sets ∂Ω D p and ∂Ω D s for the pressure and the saturation, respectively, both sets being assumed to have nonzero (d-1)-dimensional measure. In addition, ∂Ω D s coincides with the inflow boundary defined such that u•n ∂Ω < 0. Neumann boundary conditions enforce the normal component of the total velocity, namely λ(s)K∇p•n ∂Ω , or the normal component of the saturation diffusive flux, namely ǫ(s)∇π(s)•n ∂Ω . For simplicity, we only consider homogeneous Neumann boundary conditions. For a further discussion of practically important boundary conditions in the global pressure/fractional flow formulation, we refer to [START_REF] Chavent | Mathematical Models and Finite Elements for Reservoir Simulation[END_REF][START_REF] Chen | Multiphase flow simulation with various boundary conditions[END_REF].

Interface conditions

We follow here the formulation of interface conditions in [START_REF] Chavent | Mathematical Models and Finite Elements for Reservoir Simulation[END_REF][START_REF] Enchéry | Numerical approximation of a two-phase flow problem in a porous medium with discontinuous capillary forces[END_REF]. The interface conditions for the saturation enforce the continuity of the normal component of the total advective-diffusive flux, prescribe the value of s 2 to S nr if s 1 ∈ [S nr , S * ], and enforce the continuity of the capillary pressure if s 1 ∈ [S * , 1 -S wr ). This yields on Γ,

(u 1 f 1 (s 1 ) -ǫ 1 (s 1 )∇π 1 (s 1 )) •n Γ = (u 2 f 2 (s 2 ) -ǫ 2 (s 2 )∇π 2 (s 2 )) •n Γ , (6)
where n Γ is a unit normal vector to Γ, along with

s 2 = S nr if s 1 ∈ [S nr , S * ], (7) π 1 (s 1 ) = π 2 (s 2 ) if s 1 ∈ [S * , 1 -S wr ). ( 8
)
Introducing the function J : [S nr , 1 -S wr ) → [0, +∞) such that

J(s) = s if s 1 ∈ [S nr , S * ], s -π -1 2 (π 1 (s)) if s 1 ∈ [S * , 1 -S wr ), (9)
where π -1 2 denotes the reciprocal function of π 2 , the interface conditions ( 7)-( 8) yield s 1s 2 = J(s 1 ). [START_REF] Cancès | Nonlinear parabolic equations with spatial discontinuities[END_REF] Observe that J is continuous at S * . The interface conditions for the pressure enforce the continuity of the normal component of the total velocity and prescribe the continuity of the wetting phase pressures p wβ if s 1 ∈ [S nr , S * ] and that of the non-wetting phase pressures p nβ if s 1 ∈ [S * , 1 -S wr ). In the latter case, continuity of the capillary pressures holds owing to [START_REF] Bourgeat | Effective model of two-phase flow in a porous medium made of different rock types[END_REF], thereby implying also the continuity of the wetting phase pressures. Thus, the wetting-phase pressure is always continuous, in agreement with the fact that this phase is always present in both subdomains. Altogether, this yields on Γ,

-λ 1 K 1 ∇p 1 •n Γ = -λ 2 K 2 ∇p 2 •n Γ , (11)
along with

p 1 -p 2 = g(s 1 ) if s 1 ∈ [S nr , S * ], (12)
p 1 -p 2 = ĝ(s 1 , s 2 ) if s 1 ∈ [S * , 1 -S wr ), (13)
where we have introduced the functions g(s 1 ) :=

s 1 Snr f 1 (ξ)π ′ 1 (ξ)dξ + π 1 (S nr ) -π 2 (S nr ), (14)
ĝ(s 1 , s 2 ) :=

s 1 Snr (f 1 (ξ) -1) π ′ 1 (ξ)dξ - s 2 Snr (f 2 (ξ) -1) π ′ 2 (ξ)dξ. ( 15
)
Conditions ( 12) and ( 13) can be readily linked to the continuity of the wetting and non-wetting phase pressures, respectively, using the definition of the global pressure

p β = p wβ + s Snr f β (ξ)π ′ β (ξ)dξ + π β (S nr ) = p nβ + s Snr (f β (ξ) -1)π ′ β (ξ)dξ, ( 16
)
and the fact that π β = p nβp wβ . Finally, observe that the pressure jump (p 1p 2 ) is a continuous function of s 1 since g(S * ) = ĝ(S * , S nr ).

3. The sequential DG method with total velocity reconstruction 3.1. The discrete setting Let {t n } 0≤n≤N be a sequence of discrete times with t 0 = 0 and t N = T and time steps τ n := t n+1t n . Let {T h } h>0 be a family of shape-regular meshes of the domain Ω. For simplicity, meshes are kept fixed in time. A generic element of T h is denoted by T , h T stands for its diameter, and h := max T ∈T h h T represents the mesh size. We say that F is an interior face of the mesh if F has nonzero (d -1)-dimensional measure and if there are distinct T -and T + in T h such that F = ∂T -∩ ∂T + . We say that F is a boundary face if F has nonzero (d -1)-dimensional measure and if there is T ∈ T h such that F = ∂T ∩ ∂Ω. Interior faces are collected in the set F i h and boundary faces are collected in the set F ∂ h . We also set F h := F i h ∪ F ∂ h , and for F ∈ F h , h F denotes the diameter of F . Meshes can possess so-called hanging nodes, resulting for instance from the fact that the subdomains Ω 1 and Ω 2 have been meshed independently; in this case, we make the usual assumption that there is κ > 0 independent of h such that h F ≥ κh T for all F ∈ F h and for all T ∈ T h such that F ⊂ ∂T .

We assume that the meshes are fitted to the partition of Ω into Ω 1 ∪ Ω 2 . In other words, we assume that the interface Γ is exactly covered by a set of faces that is denoted by F Γ h . For further use, we define

F i * h := F i h \ F Γ h .
We also assume that the sets ∂Ω D p and ∂Ω D s where Dirichlet boundary conditions are enforced on the pressure and on the saturation, respectively, are exactly covered by boundary faces, and we denote by F D hp and F D hs the corresponding sets of faces. Since the sets Γ, ∂Ω D p , and ∂Ω D s are known a priori, the above are mild assumptions on the meshes.

For a scalar-or vector-valued function v that is possibly two-valued at F ∈ F i h , we define its jump and average at F as

[[v]] := v --v + , {v} := 1 2 (v -+ v + ), v ± := v| T ± . (17)
These definitions are extended to boundary faces

F ∈ F ∂ h by setting [[v]] = {v} := v| F . ( 18
)
For F ∈ F i h , we define n F as the unit normal vector to F pointing from T - toward T + , whereas for F ∈ F ∂ h , n F coincides with n ∂Ω . The arbitrariness in the orientation of n F is irrelevant in the sequel. For

F ∈ F Γ h , it is convenient to choose n F pointing from Ω 1 toward Ω 2 so that [[v]] = v 1 -v 2 (recall that v β := v| Ω β , β ∈ {1, 2}).

Weighted averages

In the sequel, we consider weighted averages prompted by heterogeneous diffusivity [START_REF] Di Pietro | Discontinuous Galerkin methods for anisotropic semidefinite diffusion with advection[END_REF][START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF]. Let a be a scalar-valued function representing the diffusivity; a can be two-valued on F ∈ F i h with values a T -,F and a T + ,F associated with T -and T + , respectively, whereas a is single-valued on F ∈ F ∂ h . Then, for all F ∈ F i h , introducing the weights

ω T -,F := a T + ,F a T -,F + a T + ,F , ω T + ,F := a T -,F a T -,F + a T + ,F , (19)
such that ω T -,F + ω T + ,F = 1, the weighted average of a function v is defined as

{v} ω := ω T -,F v -+ ω T + ,F v + . ( 20
)
We also denote the harmonic mean of a at F as

a F := a T -,F a T + ,F a T -,F + a T + ,F . (21)
The above definitions are extended to boundary faces by setting {v} ω := v and a F := a. Thus, the definitions ( 19) and ( 21) can still be used for F ∈ F ∂ h upon formally setting a T -,F := a and a T + ,F := +∞. When working with possibly degenerate problems, such as two-phase flow problems, the diffusivity a takes nonnegative (and not necessarily positive) values. The definitions [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] and [START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF] can still be used if only one of the quantities a T ± ,F vanishes. If both quantities vanish, the weights are set to 1 2 ; for further use, the harmonic mean is set to a global reference value for a specified below.

The sequential DG scheme

Let k ≥ 1 be an integer. At each discrete time t n with n ∈ {1, . . . , N }, approximations to the pressure, p n h , and to the saturation, s n h , are sought in the DG finite element space

V k h := {v h ∈ L 2 (Ω); ∀T ∈ T h , v h | T ∈ P k (T )}, (22)
where P k (T ) denotes the vector space of polynomials with total degree ≤ k on T . Observe that approximations for the pressure and the saturation are sought using the same polynomial order. Taking different polynomial orders is also possible, but the present choice appears to be more natural. At the initial time, an approximate saturation s 0 h is constructed by considering, for instance, the L 2 -orthogonal projection onto V k h of the initial datum for the saturation.

Step 1: pressure equation

For n ∈ {0, • • • , N -1}, with s n h ∈ V k h given from the previous time step (n ≥ 1) or by the initial datum (n = 0), we first solve the pressure equation using an interior penalty DG method, that is, we solve for

p n+1 h ∈ V k h such that for all z h ∈ V k h , T ∈T h T λ(s n h )K∇p n+1 h •∇z h - F ∈F i h ∪F D hp F n F •{λ(s n h )K∇p n+1 h } ω [[z h ]] + θn F •{λ(s n h )K∇z h } ω [[p n+1 h ]] + F ∈F i h ∪F D hp γ F σ F h F [[p n+1 h ]][[z h ]] = T ∈T h T (q w + q n ) z h + F ∈F Γ h F -θn F •{λ(s n h )K∇z h } ω + γ F σ F h F [[z h ]] G(s n h ) + F ∈F D hp F -θn F •λ(s n h )K∇z h + γ F σ F h F z h p D , (23)
where p D is the Dirichlet datum for the pressure. Notice that ( 23) is a linear problem for p n+1 h . The interface operator G(s n h ) is defined as

G(s) = g(s n h1 ), if s n h1 ∈ [S nr , S * ], ĝ(s n h1 , s n h2 ), if s n h1 ∈ [S * , 1 -S wr ), (24)
where s n hβ , β ∈ {1, 2}, denotes the restriction of s n h to the subdomain Ω β . The parameter θ ∈ {-1, 0, +1} allows to select the so-called nonsymmetric, incomplete, or symmetric version of the interior penalty DG method, and σ F > 0 is a positive parameter that must be taken larger than a minimal threshold depending on mesh regularity if θ = -1. This dependency on mesh regularity can be avoided by penalizing liftings of jumps instead of jumps as detailed in [START_REF] Bassi | A highorder accurate discontinuous finite element method for inviscid and viscous turbomachinery flows[END_REF][START_REF] Di Pietro | Discontinuous Galerkin methods for anisotropic semidefinite diffusion with advection[END_REF].

The penalty coefficient γ F and the weights ω in the consistency terms are evaluated as described in §3.2.

For F ∈ F i h , F = ∂T -∩ ∂T + , we set a T ± ,F = (λ(s n h )K)| T ± L ∞ (F ) and for F ∈ F ∂ h , F = ∂T ∩ ∂Ω, we set a = (λ(s n h )K)| T L ∞ (F ) . Then, for F ∈ F i * h ∪ F D hp , we set γ F = a F . (25)
Recall that the pressure equation does not degenerate, so that the sole motivation for using weighted averages in this equation is media heterogeneities. It is interesting to recast [START_REF] Ersland | Numerical methods for flow in a porous medium with internal boundaries[END_REF] in a weak residual form by integrating by parts the first term. Observing that for all

F ∈ F i h , [[λ(s n h )K∇p n+1 h z h ]] = {λ(s n h )K∇p n+1 h } ω [[z h ]] + [[λ(s n h )K∇p n+1 h ]]{z h } ω, (26)
with

{z h } ω := ω T + ,F z h | T -+ ω T -,F z h | T + ,

and rearranging terms yields

T ∈T h T -(∇•(λ(s n h )K∇p n+1 h ) + q w + q n )z h + F ∈F i h F n F •[[λ(s n h )K∇p n+1 h ]]{z h } ω + F ∈F ∂ h \F D hp F n F •λ(s n h )K∇p n+1 h z h + F ∈F i h ∪F D hp F -θn F •{λ(s n h )K∇z h } ω + γ F σ F h F [[z h ]] [[p n+1 h ]] ′ = 0, (27)
where for compact notation we have set

[[p n+1 h ]] ′ =      [[p n+1 h ]] if F ∈ F i * h , [[p n+1 h ]] -G(s n h ) if F ∈ F Γ h , p n+1 h -p D if F ∈ F D hp . (28)
Then, it is readily seen that the first term in [START_REF] Van Duijn | The effect of capillary forces on immiscible two-phase flow in heterogeneous porous media[END_REF] weakly enforces the volumetric PDE residual, the second one the flux continuity, the third one the Neumann boundary condition, and the last one the pressure continuity inside Ω 1 ∪ Ω 2 , the interface condition on Γ, and the Dirichlet boundary condition on ∂Ω D p . Notice in particular that the interface condition enforced weakly is

[[p n+1 h ]] = G(s n h ), (29)
which is an O(τ n )-approximation of ( 12)-(13).

Step 2: total velocity reconstruction

Once the approximate pressure p n+1 h has been computed, the total velocity u n+1 h is reconstructed in the Raviart-Thomas(-Nédélec) finite element space of degree l, l ∈ {k -1, k}, where where

RT l (T h ) := {u h ∈ H(div); ∀T ∈ T h , u h | T ∈ [P l (T )] d + xP l (T )}, (30)
H(div) = {v ∈ [L 2 (Ω)] d ; ∇•v ∈ L 2 (Ω)} and x = (x 1 , . . . , x d ) t .
The reconstruction is presented in the case of matching (that is, without hanging nodes), simplicial meshes. In the more general case, the reconstruction is performed on a matching, simplicial refinement of the original mesh T h . The reconstruction prescribes locally on each mesh element T ∈ T h the local degrees of freedom of u n+1 h as follows: For all F ∈ F h with F ⊂ ∂T and for all q ∈ P l (F ),

F (u n+1 h •n F )q = F -n F •{λ(s n h )K∇p n+1 h } ω + γ F σ F h F [[p n+1 h ]] ′ q, (31)
except on Neumann boundary faces where there is no jump contribution, and for all r ∈ [P l-1 (T )] d (there is no condition if k = 1 and l = 0),

T u n+1 h •r = - T λ(s n h )K∇p n+1 h •r + θ F ⊂∂T F ω T,F λ(s n h )n F •K•r[[p n+1 h ]] ′ . ( 32
)
These conditions fully prescribe the degrees of freedom of the vector field u n+1 h [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]. The local degrees of freedom for the spaces RT 0 and RT 1 are illustrated in Figure 2. Furthermore, the computation of u n+1 h is explicit, and thus extremely cheap. Moreover, it is readily seen that the divergence of u n+1 h is optimally convergent in the sense that for all ξ ∈ P l (T ),

T (∇•u n+1 h )ξ = T (q w + q n )ξ. ( 33
)
We refer to [START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF] for more details.

Step 3: saturation equation

To compute s n+1 h , we use a backward Euler scheme in time, together with a spatial discretization of the diffusion term by an interior penalty DG method and a DG method with Godunov flux (or, equivalently, an upwind flux since the fractional flux function is monotone) for the nonlinear advective term. Hence, we solve for

s n+1 h ∈ V k h such that for all v h ∈ V k h , T ∈T h T φτ -1 n s n+1 h v h - T ∈T h T u n+1 h f (s n+1 h )•∇v h + F ∈F h F (u n+1 h •n F )Φ n+1 hF [[v h ]] + T ∈T h T ǫ(s n h )π ′ (s n h )∇s n+1 h •∇v h - F ∈F i h ∪F D hs n F •{ǫ(s n h )π ′ (s n h )∇s n+1 h } ω [[v h ]] + θn F •{ǫ(s n h )π ′ (s n h )∇v h } ω [[s n+1 h ]] + F ∈F i h ∪F D hs δ F σ F h F [[s n+1 h ]][[v h ]] = T ∈T h T φτ -1 n s n h v h + T ∈T h T q n v h + F ∈F Γ h F -θn F •{ǫ(s n h )π ′ (s n h )∇v h } ω + δ F σ F h F [[v h ]] J(s n h1 ) + F ∈F D hs F -θn F •ǫ(s n h )π ′ (s n h )∇v h + δ F σ F h F v h s D , (34)
where s D is the Dirichlet datum for the saturation. Although many terms have been linearized by evaluating them at s n h , (34) is still a nonlinear problem because of the fractional flux being evaluated at s n+1 h . The penalty coefficient δ F and the weights ω in the consistency terms are evaluated as described in §3.2.

For F ∈ F i h , F = ∂T -∩ ∂T + , we set a T ± ,F = (ǫ(s n h )π ′ (s n h ))| T ± L ∞ (F ) and for F ∈ F ∂ h , F = ∂T ∩ ∂Ω, we set a = (ǫ(s n h )π ′ (s n h ))| T L ∞ (F ) .
Then, for F ∈ F i * h ∪ F D hs , we set γ F = a F , and the reference value for a is taken to be a * β := max s∈[Snr,1-Swr) (ǫ β (s)π ′ β (s)) where β ∈ {1, 2} is such that F ⊂ Ω β . For F ∈ F Γ h , the above penalty strategy is revised by setting the weights to 1 2 and δ F = a * 1 a * 2 /(a * 1 + a * 2 ). This way, the interface condition is weakly enforced even in the degenerate case.

The numerical flux Φ n+1 hF is evaluated as

Φ n+1 hF =          f (s n+1,↑ h ) if F ∈ F i * h , {f (s n+1 h )} ω if F ∈ F Γ h , f (s n+1 h ) if F ∈ F ∂ h \ F D hs , f (s D ) if F ∈ F D hs , (35)
where s n+1,↑ h denotes the so-called upwind value of s n+1 h , that is, for

F ∈ F i h , F = ∂T -∩ ∂T + , s n+1,↑ h = s n+1 h | T -if u n+1 h •n F ≥ 0 and s n+1,↑ h = s n+1 h | T + oth- erwise.
The upwind flux is used on all faces except on those at the interface where a weighted centered flux is used. An upwind flux could have been used as well at the interface, but this is unnecessary since the design of the penalty coefficient δ F already ensures that the interface condition is weakly penalized.

It is interesting to recast (34) in a weak residual form by integrating by parts the second and fourth terms, yielding

T ∈T h T φτ -1 n (s n+1 h -s n h ) + ∇•(u n+1 h f (s n+1 h ) -ǫ(s n h )π ′ (s n h )∇s n+1 h ) -q n v h + F ∈F i h F n F •[[-u n+1 h f (s n+1 h ) + ǫ(s n h )π ′ (s n h )∇s n+1 h ]]{v h } ω + F ∈F ∂ h \F D hs F n F •ǫ(s n h )π ′ (s n h )∇s n+1 h v h + F ∈F i h ∪F D hs F -θn F •{ǫ(s n h )π ′ (s n h )∇v h } ω + δ F σ F h F [[v h ]] [[s n+1 h ]] ′ + F ∈F i * h F ω T ↓ ,F |u n+1 h •n F |[[v h ]][[f (s n+1 h )]] + F ∈F D hs F |u n+1 h •n F |[[v h ]][[s n+1 h ]] ′ = 0, (36)
where for compact notation we have set

[[s n+1 h ]] ′ =      [[s n+1 h ]] if F ∈ F i * h , [[s n+1 h ]] -J(s n h1 ) if F ∈ F Γ h , s n+1 h -s D if F ∈ F D hs . (37)
Thus, the first term in (36) weakly enforces the volumetric PDE residual, the second one the flux continuity, the third one the Neumann boundary condition, and the fourth one the saturation continuity inside Ω 1 ∪ Ω 2 , the interface condition on Γ, and the Dirichlet boundary condition on ∂Ω D s . We now clearly see why the penalty coefficient δ F must be nonzero at the interface even in the degenerate case. Furthermore, the interface condition enforced weakly is

[[s n+1 h ]] = J(s n h1 ), (38)
which is an O(τ n )-approximation of [START_REF] Cancès | Nonlinear parabolic equations with spatial discontinuities[END_REF]. Finally, the two terms in the fifth line of (36) result from upwinding. In the first of the two terms, T ↓ denotes the mesh element downwind of F ; moreover, this term is dissipative since f is an increasing function (that is, it yields a nonnegative contribution upon taking v h = s n+1 h ). In the second term, we have used the fact that ∂Ω D s coincides with the inflow boundary to introduce the absolute value of u n+1 h •n F ; moreover, this second term is also dissipative.

Numerical results

The DG method designed in §3 is assessed on three one-dimensional test cases with synthetic data representative of secondary oil recovery problems exhibiting degeneracies (the non-wetting phase saturation vanishes in some parts of the domain) together with pressure and saturation discontinuities at the interface between the two subdomains. The computational domain is Ω = (0, 2) with the interface located at x intf = 1 so that Ω 1 = (0, 1) and Ω 2 = (1, 2). In both subdomains, we set φ = 0.2, S nr = S wr = 0, and choose the mobilities from Brooks-Corey model [START_REF] Bear | Dynamic of Fluids in Porous Media[END_REF] with parameter θ = 2, yielding for s ∈ [0, 1],

λ w (s) = (1 -s) 4 , λ n (s) = s 2 (2s -s 2 ). (39)
The absolute permeability is specified below for each test case. The capillary pressures are given for s ∈ [0, 1] by

π 1 (s) = 5s 2 , π 2 (s) = 4s 2 + 1, (40)
so that the difference in entry pressures between the two subdomains is equal to 1. The critical saturation is S * = 1 √ 5 ≃ 0.4472. There are no sources or sinks (q w = q n = 0). We enforce Dirichlet boundary conditions on the pressure in the form p| x=0 = 1.8 and p| x=2 = 1.0, and mixed boundary conditions on the saturation in the form s| x=0 = 0 and ǫ(s)π ′ (s) ds dx | x=2 = 0. The initial condition is specified below for each test case.

We use a uniform mesh of Ω composed of 80 elements in each subdomain, together with first-order polynomial approximation for both the pressure and the saturation (k = 1 in ( 22)), and the total velocity is reconstructed in the Raviart-Thomas(-Nédélec) finite element space RT 1 (l = 1 in §3.3.2). The time step is set to 0.001. We have verified that the adopted mesh size and time step yield sufficient accuracy. Furthermore, since we are approximating rough solutions with possible degeneracies, spurious oscillations in the saturation profiles triggered by Gibbs-type phenomena can and do occur. As reflected in the results below, these oscillations always remain localized in space to a few cells around shocks, and their amplitude is always very small. Thus, simulations can be performed without employing limiters. To this purpose, the mobilities and capillary pressures are extended by constant values outside the interval [0, 1].

In test case 1, the intrinsic permeability is set to K = 1 over the whole domain. The initial condition for the saturation is s(x, t = 0) = 0.9 if x ∈ (0.1, 0.9), 0 otherwise, (41) so that the non-wetting phase is initially present only in the left subdomain. Results are presented in Figure 3. The saturation profiles are typical of propagating fronts with degeneracies and capillary pressure effects. Some spurious oscillations are visible in the few cells around the propagating fronts at time t = 0.015; these oscillations always remain localized in space and are damped out already at time t = 0.045. When the saturation reaches the critical value S * at the interface, the non-wetting phase penetrates into the right subdomain, and the saturation remains discontinuous at the interface.

Regarding the pressure profiles, the total pressure is discontinuous at the interface, the capillary pressure is first discontinuous and becomes continuous once the non-wetting phase has penetrated inside the right subdomain, in agreement with the interface conditions, and the wetting-phase pressure is continuous before and after penetration, again in agreement with the interface conditions. In addition, it is interesting to observe that owing to the total pressure jump at the interface, the total pressure gradient in each subdomain is larger than the gradient that could be expected simply from the values of the Dirichlet boundary conditions for the pressure. Figure 4 compares the above results, for which the total velocity is accurately reconstructed, to those obtained by simply using the piecewise constant pressure gradient to compute the total velocity. It is seen that this latter procedure leads to substantial spurious oscillations and even to the formation of non-physical shocks at times t = 0.045 and t = 0.1. By comparing the left and right columns of Figure 4, it is observed that these oscillations are localized where the total velocity is not accurately reconstructed. Finally, at longer simulation times (e.g., t = 0.25) when the saturation profiles have become smoother, these spurious oscillations are damped out. These results clearly demonstrate the need for accurate velocity reconstruction in two-phase flow simulations.

Test case 2 extends test case 1 to the situation where the intrinsic permeability is heterogeneous, namely K 1 = 1 and K 2 = 0.1. The same initial condition is used as in test case 1. Results are presented in Figure 5. As for test case 1, the non-wetting phase can penetrate into the subdomain Ω 2 , but because K 2 is now smaller, the non-wetting phase propagates more slowly to the right in Ω 2 , thereby also inducing a different saturation profile in Ω 1 (compare the result at t = 0.25 with that of Figure 3). The saturation exhibits a jump at the interface, and a few spurious oscillations are present in the few cells before the propagating front. The total pressure is discontinuous at the interface, the capillary pressure is first discontinuous and then continuous, while the wetting-phase pressure is always continuous, in agreement with the enforced interface conditions.

Test case 3 uses the same data as test case 2 except that the initial value of the saturation for x ∈ (0.1, 0.9) is now 0.4. Thus, when the propagating front reaches the interface, the saturation is below the critical value so that the non-wetting phase cannot immediately penetrate inside Ω 2 . The nonwetting phase then accumulates upstream of the interface until it reaches eventually the critical value S * and penetrates inside Ω 2 ; this situation is typical of oil trapping effects. Furthermore, the pressure profiles are again in agreement with the the enforced interface conditions.

To conclude, the above results illustrate that the proposed DG method is capable of tackling in a robust and accurate fashion various physical phenomena relevant to two-phase flow problems in heterogeneous porous media with discontinuous capillary pressures. The nonlinear interface conditions, which are weakly enforced in the context of a sequential scheme, are well captured in the numerical simulations. The total velocity reconstruction and the weighted averages are important ingredients to achieve stable computations that can be carried over without using limiters even in degenerate cases. Our aim is now to assess more thoroughly the proposed method in the context of multi-dimensional heterogeneous problems in more complex and realistic test cases related to field studies. From top to bottom: t = 0, t = 0.008, t = 0.045, t = 0.1, and t = 0.25.

Figure 1 :

 1 Figure 1: Capillary pressures in subdomains Ω 1 and Ω 2 and critical non-wetting phase saturation s * .

Figure 2 :

 2 Figure 2: Local degrees of freedom for the spaces RT 0 and RT 1 .

Figure 6 :

 6 Figure6: Test case 3. Left: saturation profiles at various times; right: total pressure (solid), capillary pressure (dashed), and wetting-phase pressure (dotted) at different times. From top to bottom: t = 0, t = 0.008, t = 0.045, t = 0.1, and t = 0.25.

  Figure 4: Test case 1.Comparison of results obtained with (solid lines) and without (dashed lines) velocity reconstruction: saturation (left) and total velocity (right) at different times. From top to bottom: t = 0.008, t = 0.015, t = 0.045, t = 0.1, and t = 0.25.Figure5: Test case 2. Left: saturation profiles at various times; right: total pressure (solid), capillary pressure (dashed), and wetting-phase pressure (dotted) at different times. From top to bottom: t = 0, t = 0.015, t = 0.018, t = 0.03, and t = 0.25.

		0.8 1 1 0.8 0.4 1 1.2											3 4 0.7 3 1.5 0.8 4 2
		0.6 0.8 0.6 0.3											2 0.6 2 1
	s s s s	0.6 0.2											1 0.5 p u 1 p p
		0.4 0.4 0.4 0.1											0 0.4 0 0.5
		0 0.2 0 0.2 0 0.2 0											-2 -1 -2 0.2 0.3 -1 0
	0 0 0 0 -0.1	0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8	1 1 1 1	1.2 1.2 1.2 1.2	1.4 1.4 1.4 1.4	1.6 1.6 1.6 1.6	1.8 1.8 1.8 1.8	2 2 2 2	0 0 0.1 0 0 -0.5	0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8	1 1 1 1	1.2 1.2 1.2 1.2	1.4 1.4 1.4 1.4	1.6 1.6 1.6 1.6	1.8 1.8 1.8 1.8	2 2 2 2
							x x x x											x x x x
		0.8 1 1 0.8 0.4 1 1.2											3 4 0.7 3 1.5 0.8 4 2
		0.6 0.8 0.6 0.3											2 0.6 2 1
	s s s s	0.6 0.2											1 0.5 p u 1 p p
		0.4 0.4 0.4 0.1											0 0.4 0 0.5
		0 0.2 0 0.2 0 0.2 0											-2 -1 -2 0.2 0.3 -1 0
	0 0 0 0 -0.1	0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8	1 1 1 1	1.2 1.2 1.2 1.2	1.4 1.4 1.4 1.4	1.6 1.6 1.6 1.6	1.8 1.8 1.8 1.8	2 2 2 2	0 0 0.1 0 0 -0.5	0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8	1 1 1 1	1.2 1.2 1.2 1.2	1.4 1.4 1.4 1.4	1.6 1.6 1.6 1.6	1.8 1.8 1.8 1.8	2 2 2 2
							x x x x											x x x x
		0.8 1 1 0.8 0.4 1 1.2											3 4 0.7 3 1.5 0.8 4 2
		0.6 0.8 0.6 0.3											2 0.6 2 1
	s s s s	0.6 0.2											1 0.5 p u 1 p p
		0.4 0.4 0.4 0.1											0 0.4 0 0.5
		0 0.2 0 0.2 0 0.2 0											-2 -1 -2 0.2 0.3 -1 0
	0 0 0 0 -0.1	0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8	1 1 1 1	1.2 1.2 1.2 1.2	1.4 1.4 1.4 1.4	1.6 1.6 1.6 1.6	1.8 1.8 1.8 1.8	2 2 2 2	0 0 0.1 0 0 -0.5	0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8	1 1 1 1	1.2 1.2 1.2 1.2	1.4 1.4 1.4 1.4	1.6 1.6 1.6 1.6	1.8 1.8 1.8 1.8	2 2 2 2
							x x x x											x x x x
		0.8 1 1 0.8 0.4 1 1.2											3 4 0.7 3 1.5 0.8 4 2
		0.6 0.8 0.6 0.3											2 0.6 2 1
	s s s s	0.6 0.2											1 0.5 p u 1 p p
		0.4 0.4 0.4 0.1											0 0.4 0 0.5
		0 0.2 0 0.2 0 0.2 0											-2 -1 -2 0.2 0.3 -1 0
	0 0 0 0 -0.1	0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8	1 1 1 1	1.2 1.2 1.2 1.2	1.4 1.4 1.4 1.4	1.6 1.6 1.6 1.6	1.8 1.8 1.8 1.8	2 2 2 2	0 0 0.1 0 0 -0.5	0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8	1 1 1 1	1.2 1.2 1.2 1.2	1.4 1.4 1.4 1.4	1.6 1.6 1.6 1.6	1.8 1.8 1.8 1.8	2 2 2 2
							x x x x											x x x x
		0.8 1 1 1.2 1 0.8 0.4											3 4 0.7 4 0.8 3
		0.6 0.8 0.6 0.3											2 0.6 2
	s s s s	0.6 0.2											1 0.5 p u 1 p
		0.4 0.4 0.4 0.1											0 0.4 0
		0 0.2 0 0.2 0 0.2 0											-2 -1 -2 0.2 0.3 -1
	0 0 0 0 -0.1	0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8	1 1 1 1	1.2 1.2 1.2 1.2	1.4 1.4 1.4 1.4	1.6 1.6 1.6 1.6	1.8 1.8 1.8 1.8	2 2 2 2	0 0 0.1 0 0 -0.	0.2 0.2 0.2 0.2	0.4 0.4 0.4 0.4	0.6 0.6 0.6 0.6	0.8 0.8 0.8 0.8	1 1 1 1	1.2 1.2 1.2 1.2	1.4 1.4 1.4 1.4	1.6 1.6 1.6 1.6	1.8 1.8 1.8 1.8	2 2 2 2
							x x x x											x x x
	Figure 3: Test case 1. Left: saturation profiles at various times; right: total pressure
	(solid), capillary pressure (dashed), and wetting-phase pressure (dotted) at different times.
	From top to bottom: t = 0, t = 0.008, t = 0.015, t = 0.045, and t = 0.25.

Partially supported by the Groupement MoMaS (PACEN/CNRS, ANDRA, BRGM, CEA, EdF, IRSN), France 2 Partially supported by CNPq, Brazil