L Weynans
email: lisl.weynans@inria.fr

Convergence of a cartesian method for elliptic problems with immersed interfaces

Keywords: Finite-differences, cartesian method, elliptic problem, interface discontinuity, interface unknowns, discrete Green's function, convergence Différences finies, problème elliptique, méthode cartésienne, discontinuité au travers de l'interface, inconnues sur l'interface, fonction de Green discrète, convergence

. This method is based on additional unknowns located on the interface, used to express the jump conditions across the interface and discretize the elliptic operator in each subdomain separately. It is numerically second-order accurate in L ∞ -norm. We prove the convergence of the method in two cases: the original second-order method in one dimension, and a first-order version in two dimensions. The proof of convergence takes advantage of a discrete maximum principle to obtain estimates on the coefficients of the inverse matrix. More precisely, we obtain estimates for the sums of the coefficients of several blocks of the inverse matrix. Associated to the consistency error, which has different leading orders throughout the domain, these estimates lead to the convergence results.

Introduction

In this paper we aim to study the convergence of a method developed in a previous paper [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF] to solve on a Cartesian grid the following problem elliptic problem, defined on a domain Ω consisting in the union of two subdomains Ω 1 and Ω 2 , separated by a complex interface Σ (see Figure 1):

-∇.(k∇u) = f on Ω = Ω 1 ∪ Ω 2 (1) u = α on Σ (2) k ∂ u ∂ n = β on Σ (3)
with k constant on each subdomain Ω 1 and Ω 2 , assorted with Dirichlet boundary conditions on δ Ω, defined as the boundary of Ω, and where • means • 2 -• 1 . The domain Ω can have an arbitrary shape.

In the whole paper, we assume that the interface is C 2 and that the solution u of problem ( 1)-( 3) exists and is smooth enough so that our truncation error analyses are valid. We assume, by convention, that the coefficient k is larger in Ω 2 than in Ω 1 (k 2 > k 1 ), and that the vector n is the outside normal for the subdomain Ω 2 . These conventions will be used in the convergence proof. Note that others configurations than the one illustrated on this figure are possible, for instance, Ω 1 separated in several subdomains, and are all covered by our analysis. This elliptic problem with discontinuities across an interface appears in numerous physical or biological models. Among the well-known applications are heat transfer, electrostatics, incompressible flows with discontinuous densities and viscosities [START_REF] Bergmann | A sharp cartesian method for incompressible flows with large density ratios[END_REF], but similar elliptic problems arise for instance in tumour growth modelling, where one has to solve a pressure equation [START_REF] Bresch | Computational modeling of solid tumor growth: the avascular stage[END_REF], or in the modelling of electric potential in biological cells: see for instance [START_REF] Breton | Nonlinear steady-state electrical current modeling for the electropermeabilization of biological tissue[END_REF] or [START_REF] Kavian | Classical electropermeabilization modelling at the cell scale[END_REF] and [START_REF] Leguebe | A second-order cartesian method for the simulation of electropermeabilization cell models[END_REF] where the mentioned numerical method was applied. The method that we study was developed in [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF]. It is based on a finite-difference discretization and a dimension by dimension approach. In order to solve accurately the problem defined by equations ( 1) -(3) near the interface, additional unknowns are defined at the intersections of the interface with the grid, see Figure 2. These interface unknowns are used in the discretization of the elliptic operator near the interface, and prevent to derive specific finite differences formulas containing jump terms, corrective terms, or needing the inversion of a linear system, as in many other second-order Cartesian methods. In Inria order to solve the interface unknowns, the flux jump conditions are discretized and added to the linear system to solve.

In the following we will prove the convergence of the method in two cases: the original second-order method in one dimension, and a first-order version in two dimensions. This first-order version is based on the same ideas as the original method, but the discretization of the normal derivatives across the interface has only a first-order truncation error instead of a second-order for the original method.

The convergence proof is based on a discrete maximum principle, used to provide estimates of the coefficients of the inverse of the discretization matrix. To this purpose we have to prove the monotonicity of the discretization matrix. This monotonicity property is not straightforward for the second-order discretization, since the discretization matrix is not diagonally dominant, due to the discretization of the flux jump conditions across the interface. Then we obtain accurate estimates of the coefficients of the inverse matrix, block by block, in order to account for the different types of truncation errors. Combined to the truncation error expressed block by block, these estimates provide first-or second-order bounds on the numerical error.

In section 2 we describe the numerical schemes, in section 3 and 4 we present the proof for the firstorder version in two dimensions and the second-order version in one-dimension. In section 5 we compare our approach to the literature and in section 6 we present some numerical tests corroborating our analysis.

Description of the numerical schemes 2.1 Interface representation and classification of grid points

In order to describe accurately the geometric configuration in the vicinity of the interface we use the level set method introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. We refer the interested reader to [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF], [START_REF] Sethian | Evolution, implementation, and application of level set and fast marching methods for advancing fronts[END_REF] and [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] for reviews of this method. We recall here some properties that we will use in the following.

• The zero isoline of the level set function, defined here by the signed distance function ϕ:

ϕ(x) =    dist Σ (x)
outside of the interface 0 on the interface, -dist Σ (x) inside of the interface [START_REF] Bergmann | A sharp cartesian method for incompressible flows with large density ratios[END_REF] represents implicitely the interface Σ immersed in the computational domain.

• As recalled in [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF], the level-set, being a distance function, is 1-Lipshitz and almost everywhere differentiable. Moreover, if ϕ is differentiable at a point x, then it satisfies the so-called Eikonal equation at x:

||∇ϕ(x)|| = 1.
• The smoothness of the level-set is in fact strongly related to the smoothness of the interface: as proved in [START_REF] Gilbard | Elliptic Partial Differential Equations of Second Order[END_REF], p 355, if the interface is C 2 , then there exists a real r 0 > 0 such that the level-set is C 2 for all x, y such that |ϕ(x, y)| < r 0 .

• The outward normal vector of the isoline of ϕ passing on x, denoted n(x), can be expressed where ϕ is differentiable as

n(x) = ∇ϕ(x) |∇ϕ(x)| , (5)
In this paper, the level-set is defined so that n is the outside normal for the subdomain Ω 2 . 3) is discretized on a uniform Cartesian grid covering Ω 1 ∪ Ω 2 , see Figure 2. The grid spacing is denoted h. The points on the cartesian grid are named either with letters such as P or Q, or with indices such as M i, j = (x i , y j ) = (i h, j h) if one needs to have informations about the location of the point. We also denote if more convenient x P and y P the coordinates of a point P. We denote by u h i j the approximation of u at the point (x i , y j ). The set of grid points located inside the domain Ω is denoted Ω h .

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
We say that a grid point is irregular if the sign of ϕ changes between this point and at least one of its neighbors, see Figure 2. On the contrary, grid points that are not irregular are called regular grid points. The set of irregular grid nodes is denoted Ω * h . The subset Ω δ h is defined as the set of regular grid points where the stencil for the discrete elliptic operator (described in the following) crosses the isolines ϕ = δ or ϕ = -δ , with δ such that 0 < δ < r 0 . Notably ϕ is C 2 on Ω δ h with bounded derivatives. We define the interface point I i+1/2, j = (x i+1/2, j , y j ) as the intersection of the interface and the segment [M i j M i+1 j ], if it exists. Similarly, the interface point I i, j+1/2 = (x i , y i, j+1/2 ) is defined as the intersection of the interface and the segment [M i j M i j+1 ]. The set of interface points is denoted Σ h , see Figure 2 for an illustration. At each interface point we create two additional unknowns, called interface unknowns, and denoted by u 1,h i+1/2, j and u 2,h i+1/2, j , or u 1,h i, j+1/2 and u 2,h i, j+1/2 . The interface unknowns carry the values of the numerical solution on each side of the interface.

We also denote δ Ω h the set of points defined as the intersection between the grid and δ Ω. They are used to impose boundary conditions.

Second-order discretization in one dimension

We recall here the discretization presented in [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF] applied in one dimension.

• Discrete elliptic operator on a grid point M i

We use the standart three point stencil: M i and its nearest neighbors in each direction, either grid or interface points. We denote u E (resp. u W ) the value of the numerical solution on the nearest point in the east (resp. west) direction, and x E (resp. x w ) its coordinate. The discretization at point M i reads -∇.(k∇u)

h i = -k i 2 x E -x W u E -u i x E -x i - u i -u W x i -x W . (6)

Inria

The truncation error of this discretization is second-order accurate on regular points, and first-order otherwise.

• Discrete jump conditions across the interface On Figure 3, we present a prototypical situation around the interface: the interface point x int = x k+1/2 is located between the grid points x k and x k+1 and we denote dh = x k+1x k+1/2 . We assume for instance that the subdomain Ω 2 is located on the left side of the interface, and Ω 1 on the right side. The normal to the interface is oriented from the left to the right.

The left and right normal derivatives at the interface are computed with second-order formulas using three non-equidistant points:

(∂ n u 1 ) h k+1/2 = 1 + 2d d(d + 1)h (u h k+1 -u 1,h k+1/2 ) - d (1 + d)h (u h k+2 -u h k+1 ), (∂ n u 2 ) h k+1/2 = 3 -2d (1 -d)(2 -d)h (u 2,h k+1/2 -u h k ) - 1 -d (2 -d)h (u h k -u h k-1 ).
We express the jump conditions at point x k+1/2 as

u 2,h k+1/2 -u 1,h k+1/2 = α(x k+1/2 ), (7)
k 2 (∂ n u 2 ) h k+1/2 -k 1 (∂ n u 1 ) h k+1/2 = β (x k+1/2 ). (8)

First-order discretization in two dimensions

We present here the variant of the method of [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF] with the jump on the fluxes discretized with a first-order accuracy.

• Discrete elliptic operator

We use a standard five point stencil with the grid point M i, j and its nearest neighbors, interface or grid points, in each direction. More precisely, we denote u h S the value of the solution on the nearest point in the south direction, with coordinates (x S , y S ). Similarly, we define u h N , u h W and u h E and the associated coordinates (x N , y N ), (x W , y W ) and (x E , y E ). The discretization reads

-∇.(k∇u) h i, j = -∇.(k∇u) h (x i , y j ) = -k i, j u h N -u h i j x N -x i - u h i j -u h S x i -x S 2 x N -x S -k i, j u h E -u h i j y E -y j - u h i j -u h W y j -y W 2 y E -y W . (9)

L. Weynans

The truncation error of this discretization is second-order accurate on regular grid points, and firstorder on irregular grid points.

• Discrete jump conditions across the interface

We discretize the jump conditions (2) and (3) at each interface point I i+1/2, j as:

u 2,h i+1/2, j -u 1,h i+1/2, j = α(I i+1/2, j ), (10)
k 2 (∂ n u 2 ) h i+1/2, j -k 1 (∂ n u 1 ) h i+1/2, j = β (I i+1/2, j ), (11)
and similarly for each interface point I i, j+1/2 . The discretization of the normal derivatives depends of the local geometry of the interface. On Figure 4 one can observe the four possible cases that are met if h is small enough. The first intersection between the normal to the interface and the grid is located on a segment: either

[M i, j , M i, j-1 ], or [M i, j-1 , M i+1, j-1 ], or [M i, j , M i, j+1 ], or [M i, j+1 , M i+1, j+1 ].
The discrete normal derivative is computed as the normal derivative of the linear interpolant of the numerical solution on the triangle composed of the interface point I i+1/2, j and the aforementioned segment. If we denote K this triangle, (x 1 , y 1 ), (x 2 , y 2 ) and (x 3 , y 3 ) its vertices, and u 1 , u 2 and u 3 the associated values, the basis functions on the vertices for the linear interpolation write

λ j (x, y) = α j x + β j y + γ j , j = 1, 2, 3, with α j = y k -y i (x j -x k )(y j -y i ) -(x j -x i )(y i -y k ) , β j = x i -x k (x j -x k )(y j -y i ) -(x j -x i )(y i -y k ) , γ j = x k y i -x i y k (x j -x k )(y j -y i ) -(x j -x i )(y i -y k ) ,
(n x , n y ) being an approximation of the normal at the interface point. With these notations, the approximation of the normal derivative writes for instance for the interface point I i+1/2, j

(∂ n u) h i+1/2, j = (u 1 α 1 + u 2 α 2 + u 3 α 3 )n x + (u 1 β 1 + u 2 β 2 + u 3 β 3 )n y .
This discretization is only first-order accurate because it is based on a linear interpolation.

Elimination of the interface values on the Ω 2 side

We replace the variables u 2,h i+1/2, j and u 2,h i, j+1/2 by u 1,h i+1/2, j + α(I i+1/2, j ) and u 1,h i, j+1/2 + α(I i, j+1/2 ) in the equations ( 8) or [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF], and ( 6) or [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF], in order to eliminate the jump conditions ( 10) or ( 7) from the linear system. Because the jump conditions [START_REF] Brandt | Generalized local maximum principles for finite-difference operators[END_REF] or [START_REF] Ciarlet | Discrete maximum principle for finite-difference operators[END_REF] are expressed exactly, this does not change the truncation errors. In the following we denote A h the matrix of the linear system resulting from this discretization, in one or two dimensions. Peskin The local error array e and the consistency error array τ obey the same linear relationship as the numerical solution and the source terms:

A h e = τ.
For the sake of simplicity we assume that the boundary conditions on δ Ω are Dirichlet boundary conditions imposed exactly. Consequently, the local error e is zero on δ Ω h .

Inria

Convergence proof in two dimensions for the first-order version of the method

First we prove the monotonicity of the discretization matrix, then we use it to apply a discrete maximum principle to the matrix, and obtain estimates on the coefficients of the inverse matrix, block by block.

Monotonicity of the discretization matrix

Here we aim to prove that A h is monotone, that is, that all the coefficients of the inverse matrix of A h are non-negative. Let v be an array of size N + N int , corresponding to N grid points and N int interface unknowns, such that all coefficients of A h v are non-negative, which we denote A h v ≥ 0. We aim to prove that all coefficients of v are non-negative.

To this purpose, we use the following remark:

With the convention used for the normal to the interface, illustrated on Figure 5, if the minimum of v is located on an interface point, then at this interface point the discrete normal derivative in Ω 1 is positive and the discretized normal derivative in Ω 2 is negative.

Indeed, the approximation of the normal derivative is constant, because it computed from a linear interpolation on a triangle. If the minimum of v is located on the considered interface point, then the left normal derivative at this interface point is negative and the right normal derivative at this interface point is positive. Moreover, if the minimum of v is located on an interface point, and if the approximated normal derivative at this point is zero, then the three points values involved in the stencil are equal. Now we consider the minimum of v in the whole domain, interface points included. This value can either be located on a grid point in one of the subdomains Ω 1 or Ω 2 , or on an interface point.

• If the minimum is located on one border of the computational domain

We assume for instance that the minimum of v is v 1, j , located on the grid point M 1, j , such that M 2, j , M 1, j+1 and M 1, j-1 also belong to the computational domain. Therefore, the boundary is located on the left side of M 1, j . The other cases would be treated the same way. The elliptic operator inequality on this grid point yields:

4v 1, j -v 2, j -v 1, j+1 -v 1, j-1 h 2 ≥ 0,
then we have 4v 1, j ≥ v 2, j + v 1, j+1 + v 1, j-1 ≥ 3v 1, j and thus v 1, j ≥ 0. Therefore all values of v are non-negative.

• If the minimum is reached on a grid point in one subdomain sharing at least one point with δ Ω

In this case we denote (i 0 , j 0 ) the indices of the smallest component of v. We assume the grid point is a regular grid point (otherwise the formula would have slightly different weights, but the reasoning would be the same). Using the elliptic operator inequality on this point we can write:

4v i 0 , j 0 -v i 0 +1, j 0 -v i 0 -1, j 0 -v i 0 , j 0 +1 -v i 0 , j 0 -1 ≥ 0, we deduce that v i 0 +1, j 0 = v i 0 -1, j 0 = v i 0 , j 0 +1 = v i 0 , j 0 -1 = v i 0 , j 0 .
Repeating this reasoning on the neighbours of (i 0 , j 0 ), then on the neigbours of the neighbours etc, we deduce that all values in the subdomain, including the boundary values, are equal to v i 0 , j 0 . We use now the reasoning of the last paragraph to conclude that all values of v are non-negative.

• If the minimum is reached on a grid point in one subdomain which does not share any point with δ Ω

We use the notations and geometrical configuration of Figure 5. We assume, without loss of generality, that the subdomain is Ω 1 and we denote (i 0 , j 0 ) the indices of the minimum of v. With the same reasoning as in the previous paragraph, we can prove that all values in the subdomain, including the interface values, are equal to the minimum value v i 0 , j 0 .

Inria

Let us consider one of these interface values, located for instance on point I i+1/2, j . As noticed previously, the fact that all values are equal in the subdomain implies that the normal derivative at the interface point is zero. Due to the fact that A h v ≥ 0, we can write on this interface point:

0 = k 1 (∂ n v 1 ) h i+1/2, j ≤ k 2 (∂ n v 2 ) h i+1/2, j .
On the other side, because the minimum value is reached on this interface point, we have also

(∂ n v 2 ) h i+1/2, j ≤ 0. Consequently, (∂ n v 1 ) h i+1/2, j = (∂ n v 2 ) h i+1/2, j = 0,
and the values of the grid points involved in the stencil for (∂ n v 2 ) h i+1/2, j are equal to the value of the interface point. It means that there are two grid points in the subdomain Ω 2 where the minimum value is reached.

In this paper we have considered so far that there are only two subdomains. Therefore, we know that the subdomain Ω 2 has a non-void intersection with δ Ω, and we use the reasoning of the previous paragraph to conclude. In the case where more subdomains were considered, we would distinguish whether the subdomain Ω 2 has a non-void intersection with δ Ω or not. If not the case, then we would apply again the reasoning of this paragraph, switching from subdomains to subdomains, until finding a subdomain whose intersection with δ Ω is non empty.

• If the minimum is located on one interface point Without loss of generality, we assume that the minimum is located on I i+1/2, j . On this interface point we have the two relationships

(∂ n v 2 ) h i+1/2, j ≤ 0, (∂ n v 1 ) h i+1/2, j ≥ 0.
Furthermore, we know that

k 2 (∂ n v 2 ) h i+1/2, j -k 1 (∂ n v 1 ) h i+1/2, j ≥ 0.
We infer from the previous inequalities that (∂ n v 2 ) h i+1/2, j = (∂ n v 1 ) h i+1/2, j = 0. Therefore there are at least two grid points in each subdomain where the minimum value is reached. We can then follow the reasoning of one of the two last paragraphs.

We have proven that if A h v is non-negative, then v is also non-negative. This property has two implications:

• A h is invertible. Indeed, let us assume that an array v is such that A h v = 0. It means that both A h v and A h (-v) are non-negative. Consequently, we have v ≥ 0 and -v ≥ 0, thus v = 0.

• All values of A -1 h are non-negative.

Notations for discrete Green functions

In the following, the letters P and Q represent either discretization points (on the grid or on the interface) or their indices in the global numerotation of the matrix. For instance, we denote u(P) the coefficient of the row of u with the same index than the point P. Similarly, A h U(P) represents the coefficient of the P-th row of the array A h U, and A h (P, Q) is the coefficient of the P-th row and Q-th column of the matrix A h . We also define respectively by A h (:, Q) and A h (P, :) the Q-th column and the P-th row of the matrix A h .

L. Weynans

For each

Q ∈ Ω h ∪ Σ h , define the discrete Green's function G h (:, Q) = G h (P, Q) P∈Ω h ∪Σ h ∪δ Ω h
as the solution of the discrete problem:

   A h G h (:, Q)(P) = 0, P = Q 1, P = Q P ∈ Ω h ∪ Σ h , G h (P, Q) = 0, P ∈ δ Ω h . (12)
The matrix A h being monotone, all values of G h (:, Q) are positive. For homogeneous Dirichlet boundary conditions we can write the solution of the numerical problem as a sum of the source terms multiplied by the values of the discrete Green function:

u h (P) = ∑ Q∈Ω h ∪Σ h G h (P, Q) (A h u h )(Q), ∀P ∈ Ω h ∪ Σ h .
Now we present the result of Ciarlet in [START_REF] Ciarlet | Discrete maximum principle for finite-difference operators[END_REF], based on a discrete maximum principle, slightly modified to be adapted to our discretization matrix.

Lemma 1. Let S be a subset of points, W a discrete function with W ≡ 0 on δ Ω h , and α > 0 such that:

(A h W )(P) ≥ 0 ∀P ∈ Ω h ∪ Σ h , (A h W )(P) ≥ α -i ∀P ∈ S.
Then

∑ Q∈S G h (P, Q) ≤ α i W (P), ∀P ∈ Ω h ∪ Σ h .

Proof:

Using the definition of the discrete Green function, we can write

A h ∑ Q∈S G h (:, Q) (P) = 1 if P ∈ S, 0 if P ∈ S.
Therefore,

A h W -α -i ∑ Q∈S G h (:, Q) (P) ≥ 0, ∀P ∈ Ω h ∪ Σ h .
As all coefficients of the inverse of A h are non-negative, it leads to

W (P) -α -i ∑ Q∈S G h (P, Q) ≥ 0, ∀P ∈ Ω h ∪ Σ h ,
and finally we obtain an estimate of ∑ Q∈S G h (:, Q) in terms of the coefficients of W :

∑ Q∈S G h (P, Q) ≤ α i W (P), ∀P ∈ Ω h ∪ Σ h .
Our first contribution is to notice that this result can be generalized to several subsets, with both positive and negative lower bounds: Inria Lemma 2. Let S and S be two subsets of points, W a discrete function with W ≡ 0 on δ Ω h , and α > 0, β > 0 such that:

   (A h W )(P) ≥ 0, ∀P ∈ Ω h ∪ Σ h \ S, (A h W )(P) ≥ α -i , ∀P ∈ S, (A h W )(P) ≥ -(β -j ), ∀P ∈ S. Then ∑ Q∈S G h (P, Q) ≤ α i W (P) + α i β -j ∑ Q∈ S G h (P, Q), ∀P ∈ Ω h ∪ Σ h .

Proof:

Using the definition of the discrete Green functions, we can write

A h W (P) ≥ A h α -i ∑ Q∈S G h (:, Q) -β -j ∑ Q∈ S G h (:, Q) (P), ∀P ∈ Ω h ∪ Σ h .

As all coefficients of A -1

h are non-negative, it leads to

W (P) -α -i ∑ Q∈S G h (P, Q) + β -j ∑ Q∈ S G h (P, Q) ≥ 0, ∀P ∈ Ω h ∪ Σ h ,
and finally we obtain the following bound:

∑ Q∈S G h (P, Q) ≤ α i W (P) + α i β -j ∑ Q∈ S G h (P, Q), ∀P ∈ Ω h ∪ Σ h .

Convergence proof

Here we obtain linear relationships between the sums of the coefficients of the blocks of the inverse matrix corresponding to Ω h , Ω * h , Ω δ h and Σ h . Then we combine these relationships to obtain the desired estimates.

Estimates for blocks of grid points and interface points

• Let us define the function f as:

f (x, y) = A -e -ϕ(x,y) ,
with ϕ the signed distance to the interface, negative in Ω 2 and positive in Ω 1 . We choose A such that

f (x, y) ≥ 0, ∀(x, y) ∈ Ω.
Because of ( 5), the function f satisfies:

∂ n f (x, y) = e -ϕ(x,y) (∂ x ϕ n x + ∂ y ϕ n y ) = 1
, ∀(x, y) ∈ Σ.

Therefore:

k 2 ∂ n f (x, y) -k 1 ∂ n f (x, y) ≥ (k 2 -k 1 )e -||ϕ|| ∞ , ∀(x, y) ∈ Σ.

RR n°8872

We denote by B(x M , y M , η) the sphere of center M = (x M , y M ) and radius η. Because the interface Σ is smooth and f is smooth near Σ, there exists a positive real η such that for each point M belonging to Σ, and for all (x, y) ∈ B(x M , y M , η),

k 2 ∂ n f (x, y) -k 1 ∂ n f (x, y) ≥ 1 2 k 2 ∂ n f (x M , y M ) -k 1 ∂ n f (x M , y M ) , ≥ 1 2 (k 2 -k 1 )e -||ϕ|| ∞ . ( 13
)
We consider the following function, positive and infinitely differentiable with compact support:

g(x, y) = exp (-1 1-r 2 ) if r 2 ≤ 1, with r 2 = x 2 +y 2 η 2 , 0 elsewhere. ( 14
)
We define the function F as

F(x, y) = R 2 g(x -s, y -t) f (s,t)ds dt.
If we apply the discrete elliptic operator ( 9) to F we get in the case of a regular point M i, j :

-∇.(k∇F) h i, j = k h 2 4F(x i , y j ) -F(x i + h, y j ) -F(x i -h, y j ) -F(x i , y j + h) -F(x i , y j -h) , = R 2 -∇.(k∇g) h (x i -s, y i -t) f (s,t)ds dt.
The notation -∇.(k∇g)

h (x i -s, y i -t) denotes a discretization similar to -∇.(k∇g) h i, j
, but shifted from M i, j to (x is, y it). A similar formula would be obtained in the case of an irregular grid point. The term -∇.(k∇g) h is consistent with -∇.(k∇g) at first or second-order, depending if we consider a regular or irregular grid point, therefore there exists a positive constant C 1 depending only on f , g, k and Ω, such that

-∇.(k∇F) h i, j ≥ -C 1 , ∀ M i, j ∈ Ω h . (15)
Then with a simple change of variables, let us notice that

F(x + h, y) = R 2 g(x + h -s, y -t) f (s,t)ds dt = R 2 g(x -s, y -t) f (s + h,t)ds dt.
We apply the same treatment on every term appearing in the discretization of the jump of fluxes at the interface and obtain for instance on the interface point I i+1/2, j :

k 2 (∂ n F 2 ) h i+1/2, j -k 1 (∂ n F 1 ) h i+1/2, j = R 2 g(x i+1/2 -s, y j -t) k 2 (∂ n f 2 ) h -k 1 (∂ n f 1 ) h (s,t)ds dt, = B(x i+1/2 ,y j ,η) g(x i+1/2 -s, y j -t) k 2 (∂ n f 2 ) h -k 1 (∂ n f 1 ) h (s,t)ds dt. The notation k 2 (∂ n f 2 ) h -k 1 (∂ n f 1 ) h (s,t) denotes a discretization similar to k 2 (∂ n f 2 ) h i+1/2, j - k 1 (∂ n f 1 ) h
i+1/2, j , but shifted from I i+1/2, j to (s,t).

Inria

Because the discretization of fluxes is consistent and f is smooth across the interface, we can write for h small enough:

k 2 (∂ n f 2 ) h -k 1 (∂ n f 1 ) h (s,t) ≥ 1 2 k 2 (∂ n f 2 ) -k 1 (∂ n f 1 ) (s,t), ∀(s,t) ∈ B(x i+1/2 , y j , η).
Therefore, because of the relationship (13),

k 2 (∂ n F 2 ) h i+1/2, j -k 1 (∂ n F 1 ) h i+1/2, j ≥ 1 4 (k 2 -k 1 )e -||ϕ|| ∞ B(x i+1/2 ,y j ,η) g(x i+1/2 -s, y j -t)ds dt =C 2 > 0. (16)
Now, we denote W the array of the values of F discretized on the grid and interface points, with W ≡ 0 on δ Ω h . The relationship 16 is still applicable for W , as well as the inequality 15 for grid points whose stencil does not contain points belonging to δ Ω h . Let us consider a grid point M i, j whose stencil contains points belonging to δ Ω h . For instance we assume that M i, j-1 belongs to δ Ω h . We can write, because F(x i , y jh) ≥ 0:

-∇.(k∇W ) h i, j = k h 2 4W i, j -W i+1, j -W i-1, j -W i, j+1 -W i, j-1 , = k h 2 4F(x i , y j ) -F(x i + h, y j ) -F(x i -h, y j ) -F(x i , y j + h) , ≥ k h 2 4F(x i , y j ) -F(x i + h, y j ) -F(x i -h, y j ) -F(x i , y j + h) -F(x i , y j -h) , ≥ -C 1 .
Therefore all inequalities obtained for F are still valid for W .

To sum up, there exists two strictly positive constants, C 1 and C 2 depending only on f , g, k and Ω, such that

(A h W )(P) ≥ -C 1 , ∀P ∈ Ω h , (A h W )(P) ≥ C 2 , ∀P ∈ Σ h .

Using lemma (2), it means that

∑ Q∈Σ h G h (:, Q) ≤ W C 2 + C 1 C 2 ∑ Q∈Ω h G h (:, Q). (17)
• Now we consider the exact solution ū of system (1)-( 3), with f = 1, α = β = 0 and u = 0 on δ Ω. We assume that Ω and Γ are smooth enough so that ū exists and is smooth enough for our analysis.

We define the array W as the discretisation of ū on the grid and interface points, with W ≡ 0 on δ Ω h . The discretization of the elliptic operator and the fluxes is consistent at least with first-order accuracy, thus for h small enough, we can write that

-∇.(k∇ W ) h i, j ≥ 1 2 , ∀M i, j ∈ Ω h , k 2 (∂ n W 2 ) h i+1/2, j -k 1 (∂ n W 1 ) h i+1/2, j = O(h), ∀I i+1/2, j ∈ Σ h , k 2 (∂ n W 2 ) h i, j+1/2 -k 1 (∂ n W 1 ) h i, j+1/2 = O(h), ∀I i, j+1/2 ∈ Σ h .

L. Weynans

Therefore, there exists two strictly positive constants, C 3 and C 4 depending only on f , g, k and Ω, such that

(A h W )(P) ≥ C 3 , ∀P ∈ Ω h , (A h W )(P) ≥ -C 4 h, ∀P ∈ Σ h ,
and using lemma (2), it leads to:

∑ Q∈Ω h G h (:, Q) ≤ W C 3 + C 4 C 3 h ∑ Q∈Σ h G h (:, Q). (18)
Combining ( 17) and ( 18) yields

∑ Q∈Σ h G h (:, Q) ≤ W C 2 + C 1 C 2 C 3 W +C 4 h ∑ Q∈Σ h G h (:, Q)
Therefore, for h small enough, we obtain:

∑ Q∈Σ h G h (:, Q) ≤ O(1), (19)
∑ Q∈Ω h G h (:, Q) ≤ O(1). (20)

Estimates for blocks of grid points in Ω *

h

• Let P = M i 0 , j 0 be a grid point belonging to Ω δ h . We consider the fonction FP (x, y) = ln( C r P (x, y)

),

with r P (x, y) = (xx i 0 ) 2 + (yy j 0 ) 2 + h 2 , and C such that FP (x, y) > 0 for all (x, y) ∈ Ω.

Without loss of generality, we assume in the following that x i 0 = y j 0 = 0. We can prove that for every regular grid point M i, j

-∇.(k∇ FP ) h i, j ≥ 0, (21)
and in particular, for the point P itself,

-∇.(k∇ FP ) h i 0 , j 0 ≥ C 5 h 2 , ( 22
)
with C 5 a strictly positive constant.

Proof: On a regular grid point we can write:

-∇.(k∇ FP ) h i, j = k h 2 ln r i-1, j r i+1, j r i, j-1 r i, j+1 r 4 i, j = k 2h 2 ln r 2 i-1, j r 2 i+1, j r 2 i, j-1 r 2 i, j+1 r 8 i, j , = k 2h 2 ln (x i + h) 2 + y 2 j + h 2 (x i -h) 2 + y 2 j + h 2 x 2 i + (y j + h) 2 + h 2 x 2 i + (y j -h) 2 + h 2 (x 2 i + y 2 j + h 2 ) 4
.

Inria

Moreover,

(x i + h) 2 + y 2 j + h 2 (x i -h) 2 + y 2 j + h 2 = (x 2 i + y 2 j + h 2 ) 2 + 3h 4 -2h 2 x 2 i + 2h 2 y 2 j , x 2 i + (y j + h) 2 + h 2 x 2 i + (y j -h) 2 + h 2 = (x 2 i + y 2 j + h 2 ) 2 + 3h 4 -2h 2 y 2 j + 2h 2 x 2 i .
Therefore

r 2 i-1, j r 2 i+1, j r 2 i, j-1 r 2 i, j+1 = [(x 2 i + y 2 j + h 2 ) 2 + 3h 4 ] 2 -4h 4 (x 2 i -y 2 j ) 2 .
We develop the term and remark that

(x 2 i + y 2 j + h 2 ) 2 + 3h 4 2 -4h 4 (x 2 i -y 2 j ) 2 ≥ (x 2 i + y 2 j + h 2 ) 4 ,
which gives us [START_REF] Leguebe | A second-order cartesian method for the simulation of electropermeabilization cell models[END_REF]. The relationship ( 22) is directly obtained by using x i = y j = 0 in the first formula of the proof.

The considered point P belongs to Ω δ h , thus the interface is at a distance bounded independently of h from P, and the function FP is C 2 with derivatives bounded independently of h on irregular grid points and interface points. Thus one can prove with Taylor expansions that there exists strictly positive constants C 6 and C 7 such that

-∇.(k∇ FP ) h i, j ≥ -C 6 , ∀ M i, j ∈ Ω * h , (23)
k 2 (∂ n F2 P ) h i+1/2, j -k 1 (∂ n F1 P ) h i+1/2, j ≥ -C 7 , ∀ I i+1/2, j ∈ Σ h , (24)
k 2 (∂ n F2 P ) h i, j+1/2 -k 1 (∂ n F1 P ) h i, j+1/2 ≥ -C 7 , ∀ I i, j+1/2 ∈ Σ h . (25)
To sum up the previous lines, if we denote WP the array of the values of FP discretized on the grid and interface points, with WP ≡ 0 on δ Ω h , there exists three strictly positive constants, C 5 , C 6 and C 7 , such that

(A h WP )(Q) ≥ 0, ∀Q ∈ Ω h \ Ω * h (A h WP )(P) ≥ C 5 h 2 , (A h WP )(Q) ≥ -C 6 ∀Q ∈ Ω * h , (A h WP )(Q) ≥ -C 7 ∀Q ∈ Σ h ,
Therefore we obtain for each point P belonging to

Ω δ h G h (:, P) ≤ h 2 C 5 WP (:) + h 2 C 7 C 5 ∑ Q∈Σ h G h (:, Q) + h 2 C 6 C 5 ∑ Q∈Ω * h G h (:, Q). ( 26
)
We want to sum this relationship for all points P belonging to Ω δ h .

Let us denote N Ω = diam(Ω) h , which is an upper bound of the number of grid points in each direction. We split Ω δ h into two parts:

the first one denoted Ω δ ,1 h , such that the tangent to the interface at the nearest point of the isoline ϕ = ±δ makes an angle belonging to [-π/4, π/4] ∪ [3π/4, 5π/4],

L. Weynans

and the second one denoted Ω δ ,2 h comprising the remaining points.

We detail the reasoning for Ω δ ,1 h . For any point Q ∈ Ω h and any point P ∈ Ω δ ,1 h , we consider the projection P = (x P , y Q ) of P on the horizontal line going through Q. Such a point satisfies the following relationship:

r P (Q) ≥ r P (Q) = (x P -x Q ) 2 + h 2 and thus ln( C r P (Q) ) ≤ ln( C r P (Q) ) = ln( C (x P -x Q ) 2 + h 2
).

Because the isolines ϕ = ±δ have dimension 1, and because of the definition of Ω δ ,1 h , there exists an integer M > 0 independent of h, such that

∀ 1 ≤ i ≤ N Ω , card P, P ∈ Ω δ ,1 h , ih ≤ x P -x Q ≤ (i + 1)h ≤ M. Consequently for all Q ∈ Ω h , ∑ P∈Ω δ ,1 h WP (Q) = ∑ P∈Ω δ ,1 h ln C r P (Q) ≤ ∑ P∈Ω δ ,1 h ln C r P (Q) ≤ M N Ω ∑ i=1 ln C ih ≤ M h N Ω ∑ i=1 ih (i-1)h ln C x dx = M h N Ω h 0 ln C x dx =0(1) = O( 1 h ).
A similar estimation can be obtained for Ω δ ,2 h , using projections on vertical lines rather than on horizontal lines.

If we sum the inequality (26) on all points in Ω δ

h and use the estimate [START_REF] Huang | Convergence analysis of the immersed interface method[END_REF] for Σ h , we obtain

∑ P∈Ω δ h G h (:, P) ≤ O(h) + O(h) ∑ Q∈Ω * h G h (:, Q). (27)
• We define the function f by :

f (x, y) =    B -1 if |ϕ(x, y)| ≤ h/2, B -e A(|ϕ(x,y)|-h/2) if h/2 ≤ |ϕ(x, y)| ≤ δ , B -e A(δ -h/2) if δ ≤ |ϕ(x, y)|,
with ϕ the signed distance to the interface, negative in Ω 2 and positive in Ω 1 . This function is Lipschitz-continuous on the whole domain. It is also twice differentiable with bounded derivatives, excepted on the isolines of the level-set function ϕ = ±δ , ±h/2. Thus the discrete elliptic operator will be bounded for all grid points, excepted for the grid points in Ω δ h and grid points near the interface, including Ω * h , because the stencil for these points crosses these isolines. The function f satisfies for all (x, y) such that h/2 < |ϕ(x, y)| < δ :

-∇.(k∇ f ) (x, y) = k A 2 (∂ x ϕ) 2 + (∂ y ϕ) 2 = 1 ± A∇.(∇ϕ)(x, y) e A(|ϕ(x,y)|-h/2) ,
Inria because ϕ is the signed distance function. The sign ± in this formula depends on the subdomain to which (x, y) belongs. We choose A and B such that:

kA 2 ± A∇.(k∇ϕ)(x, y) ≥ 1, ∀(x, y) such that |ϕ(x, y)| ≤ δ , f (x, y) ≥ 0, (x, y) ∈ Ω.
For all regular grid points M i, j belonging to Ω h \ Ω δ h , with |ϕ(x i , y j )| < δ , we thus have:

-∇.(k∇ f ) h i, j ≥ 0. ( 28
)
On the other side, for all regular grid points M i, j belonging to

Ω h \ Ω δ h , with |ϕ(x i , y j )| > δ , because the function f is constant, we have -∇.(k∇ f ) h i, j = 0. ( 29
)
The function ϕ is Lipschitz continuous on grid points in Ω δ h , so there exists a strictly positive constant C 8 such that

-∇.(k∇ f ) h i, j ≥ - C 8 h , ∀ M i, j ∈ Ω δ h . (30)
Let us consider a grid point M i, j whose stencil for the elliptic operator crosses one of the isolines

|ϕ| = h/2. If |ϕ(M i, j )| ≤ h/2, then f (M i, j ) = B -1.
Thus the values of the other points involved in the stencil are smaller than f (M i, j ), meaning that -∇.(k∇ f )

h i, j ≥ 0. If |ϕ(M i, j )| > h/2, then
there is at least on point in the stencil satisfying |ϕ| < h/2 and the value of f on such a point is smaller than it would be if f was not truncated at the value B -1. Consequently, we have also in

this case -∇.(k∇ f ) h i, j ≥ 0.
The discontinuity in the first derivative of f on isolines |ϕ| = h/2 is defined such that the discrete elliptic operator applied to f on an irregular point will be computed with at least two first-order derivatives whose value differs by a O(1) amplitude. Consequently the discrete elliptic operator applied to f on an irregular point will scale like 1 h .

As a consequence of the previous lines, there exists a strictly positive constant C 9 such that

-∇.(k∇ f ) h i, j ≥ C 9 h , ∀ M i, j ∈ Ω * h . (31)
Finally, the discrete normal derivative of f computed with the formula involving points in Ω 2 (resp. Ω 1 ) is positive (resp. negative), therefore

k 2 (∂ n f 2 ) h i+1/2, j -k 1 (∂ n f 1 ) h i+1/2, j ≥ 0, ∀ I i+1/2, j ∈ Σ h , (32)
k 2 (∂ n f 2 ) h i, j+1/2 -k 1 (∂ n f 1 ) h i, j+1/2 ≥ 0, ∀ I i, j+1/2 ∈ Σ h . (33)

L. Weynans

Therefore, if we denote W the array of the values of f discretized on the grid and interface points, with W ≡ 0 on δ Ω h , there exists strictly positive constants C 8 and C 9 such that

(A h WP )(Q) ≥ 0, ∀Q ∈ Ω h \ Ω δ h (A h WP )(Q) ≥ - C 8 h , ∀Q ∈ Ω δ h (A h WP )(Q) ≥ C 9 h ∀Q ∈ Ω * h , (A h WP )(Q) ≥ 0 ∀Q ∈ Σ h ,
and we conclude that

∑ Q∈Ω * h G h (:, Q) ≤ C 8 C 9 ∑ Q∈Ω δ h G h (:, Q) + h C 9 W (:). (34)
Combining ( 27) and ( 34) we obtain

∑ Q∈Ω * h G h (:, Q) ≤ C 8 C 9 O(h) + O(h) ∑ Q∈Ω * h G h (:, Q) + h C 9 W (:).
Therefore, for h small enough, we obtain:

∑ Q∈Ω * h G h (:, Q) ≤ O(h). (35)

Convergence result

Finally, we obtain an estimate of the local error on every point P in Ω h ∪ Σ h , with ū the exact solution:

| ū(P) -u h (P)| = | ∑ Q∈Ω h ∪Σ h G h (P, Q)τ(Q)|, ≤ | ∑ Q∈Ω * h G h (P, Q)τ(Q)| + | ∑ Q∈Ω h \Ω * h G h (P, Q)τ(Q)| + | ∑ Q∈Σ h G h (P, Q)τ(Q)|, ≤ O(h)| ∑ Q∈Ω * h G h (P, Q)| + O(h 2 )| ∑ Q∈Ω h \Ω * h G h (P, Q)| + O(h)| ∑ Q∈Σ h G h (P, Q)|, ≤ O(h)O(h) + O(h 2 )O(1) + O(h)O(1) = O(h).
which proves that the numerical solution converges with first-order accuracy to the exact solution in L ∞ -norm.

4 Convergence proof for the one-dimensional case

Monotonicy of the discretization matrix

Here we aim to prove that A h is monotone in spite of the fact that the matrix A h is not diagonally-dominant in the second-order version, due to the discretization terms near the interface. Let v be an array of size N +N int corresponding to N grid points and N int interface unknowns such that A h v ≥ 0. Let us assume that the minimum of v is located on an interface point x int = x k+1/2 . We will prove that, with the notations and orientation of the normal defined on Figure 3, the left normal derivative at this interface point is negative Inria and the right normal derivative at this interface point is positive. Once we have proven this property, the proof of monotonicity of the matrix is exactly the same as in two dimensions, so we do not re-write it.

The left normal derivative at x int is discretized by

(∂ n v 2 ) h k+1/2 = 3 -2d (1 -d)(2 -d)h (v k+1/2 -v k ) - 1 -d (2 -d)h (v k -v k-1 )
.

By hypothesis A h v ≥ 0 hence - v k+1/2 -v k (1 -d)h - v k -v k-1 h ≥ 0, therefore (∂ n v 2 ) h k+1/2 ≤ 2 -d (1 -d)(2 -d)h (v k+1/2 -v k ) ≤ 0.
Moreover, if one can prove that the normal derivative is zero, with the last inequality we can deduce that v k = v k+1/2 . Similarly, the right normal derivative at x k+1/2 is discretized by

(∂ n v 1 ) h k+1/2 = 1 + 2d d(d + 1)h (v k+1 -v k+1/2 ) - d (1 + d)h (v k+2 -v k+1 ). By hypothesis A h v ≥ 0 hence - v k+2 -v k+1 h - v k+1 -v k+1/2 dh ≥ 0, therefore (∂ n v 1 ) h k+1/2 ≥ 1 + d d(d + 1)h (v k+1 -v k+1/2 ) ≥ 0.
Again, if the normal derivative is zero, then v k+1 = v k+1/2 .

Second-order convergence

With exactly the same reasoning than in subsection 3.3 we can prove the estimates [START_REF] Kavian | Classical electropermeabilization modelling at the cell scale[END_REF], [START_REF] Huang | Convergence analysis of the immersed interface method[END_REF] and [START_REF] Breton | Nonlinear steady-state electrical current modeling for the electropermeabilization of biological tissue[END_REF]. We use them to obtain an estimate of the local error on every point P in Ω h ∪ Σ h :

| ū(P) -u h (P)| = | ∑ Q∈Ω h ∪Σ h G h (P, Q)τ(Q)|, ≤ | ∑ Q∈Ω * h G h (P, Q)τ(Q)| + | ∑ Q∈Ω h \Ω * h G h (P, Q)τ(Q)| + | ∑ Q∈Σ h G h (P, Q)τ(Q)|, ≤ O(h)| ∑ Q∈Ω * h G h (P, Q)| + O(h 2 )| ∑ Q∈Ω h \Ω * h G h (P, Q)| + O(h 2 )| ∑ Q∈Σ h G h (P, Q)|, ≤ O(h)O(h) + O(h 2 )O(1) + O(h 2 )O(1) = O(h 2 ).
which proves that the numerical solution converges with second-order accuracy to the exact solution in L ∞ -norm.

Discussion

Numerous numerical methods have been developed for solving the problem (1) -( 3), leading to a secondorder accuracy in maximum norm, among them:

• the pioneering work of Mayo in 1984 [START_REF] Mayo | The fast solution of Poisson's and the biharmonic equations on general regions[END_REF], where an integral equation was derived to solve elliptic interface problems with piecewise coefficients. A second-order and fourth-order Cartesian gridbased boundary integral method for an interface problem of the Laplace equation on closely packed cells was also proposed recently in [START_REF] Ying | A cartesian grid-based boundary integral method for an elliptic interface problem on closely packed cells[END_REF].

• the very well known Immersed Interface Method (IIM) of LeVeque and Li (1994) [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF], and its developments, among them: the fast IIM algorithm of Li [START_REF] Li | A fast iterative algorithm for elliptic interface problems[END_REF] for elliptic problems with piecewise constant coefficients, the Explicit Jump Immersed Interface Method (EJIIM) by Wiegmann and Bube [START_REF] Wiegmann | The explicit jump immersed interface method: finite difference method for PDEs with piecewise smooth solutions[END_REF], the Decomposed Immersed Interface Method (DIIM) by Bethelsen [START_REF] Bethelsen | A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions[END_REF], and the MIIM (maximum principle preserving) by Li and Ito [START_REF] Li | Maximum principle preserving schemes for interface problems with discontinuous coefficients[END_REF].

• the Matched Interface and Boundary (MIB) method [START_REF] Zhou | High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources[END_REF], [START_REF] Zhou | On the fictitious-domain and interpolation formulations of the matched interface and boundary (MIB) method[END_REF], [START_REF] Geng | A two-component matched interface and boundary (mib) regularization for charge singularity in implicit solvation[END_REF], introduced by Zhou et al. : the solution on each side of the interface is extended on fictitious points on the other side. These fictitious values are computed by iteratively enforcing the lowest order interface jump conditions. This method can provide finite-difference schemes of arbitrary high order.

• the Coupling Interface Method, proposed by Chern and Shu [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF], where the discretizations on each subdomain are coupled through a dimension by dimension approach using the jump conditions.

• Recently, Guittet et al. proposed in [START_REF] Guittet | Solving elliptic problems with discontinuities on irregular domains â the Voronoi Interface Method[END_REF] to add degrees of freedom close to the interface and use a Voronoi partition centered at each of these points to discretize the equations in a finite volume approach. Doing so, they obtain a symmetric positive definite linear system and a second-order convergence of the solution.

• In the context of finite element methods, which is quite different from the methodology used here, numerous developments have also been done on cartesian grids, for instance [START_REF] Guo | A group of immersed finite-element spaces for elliptic interface problems[END_REF].

Other classes of Cartesian methods also exist, only first order accurate for interface problems in the case of interface problems, but simpler to implement: Gibou et al. ( [START_REF] Gibou | A second order accurate symmetric discretization of the Poisson equation on irregular domains[END_REF], [START_REF] Gibou | A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem[END_REF]). Let us also mention a new approach to solve a Dirichlet problem by a finite difference analog of the boundary integral equations, presented in [START_REF] Beale | Solution of the dirichlet problem by a finite difference analog of the boundary integral equation[END_REF]. In this paper, the double layer potential is thought as the solution of an interface problem similar to the one considered in this paper.

Concerning the discretization requirements needed to get a second-order spatial convergence, it has been noted since the introduction of Cartesian grid methods that a O(h) truncation error at the points near the interface is enough to get a O(h 2 ) convergence in maximum norm if the discretization is second-order on the regular grid points. However, in the literature, only few works have been devoted to the study of the second-order convergence of Cartesian grid methods for interface problems.

For one-dimensional methods, Huang and Li performed in [START_REF] Huang | Convergence analysis of the immersed interface method[END_REF] a convergence analysis of the IIM, using non-negative comparison functions, and in [START_REF] Wiegmann | The explicit jump immersed interface method: finite difference method for PDEs with piecewise smooth solutions[END_REF] Wiegmann and Bube presented a proof of convergence for one-dimensional problems with piecewise constant coefficients, using a detailed analysis and identification of the coefficients of the matrices involved. In [START_REF] Kavian | Classical electropermeabilization modelling at the cell scale[END_REF], a convergence proof was established in one-dimension for a variant of the method studied in this paper, applied in the context of electropermeabilization models. But this proof was based on a row by row analysis of the discretization matrix, in order to obtain estimates of the coefficients of the inverse matrix. This technique would not be tractable in two dimensions, due to its complexity.

For two-dimensional methods, Beale and Layton [START_REF] Beale | On the accuracy of finite difference methods for elliptic problems with interfaces[END_REF] proved in two-dimensions the second-order convergence for piecewise constant diffusion coefficients, using the fact that a grid function located near the Inria interface can be written as the divergence of a function smaller in norm. In [START_REF] Li | Accurate solution and gradient computation for elliptic interface problems with variable coefficients[END_REF], Li et al. proved the second-order convergence, for the solution and its gradient, in the case of an augmented method, where the jump in the normal derivative of the solution is considered as an additional unknown. The interface problem was rewritten as a new PDE consisting in a leading Laplacian operator plus lower derivatives terms near the interface. With this reformulation it was possible to use the result of [START_REF] Beale | On the accuracy of finite difference methods for elliptic problems with interfaces[END_REF] to prove the convergence. Li and Ito proved in [START_REF] Li | Maximum principle preserving schemes for interface problems with discontinuous coefficients[END_REF] the second-order convergence of their MIIM, using the maximum principle. The proof uses a technical condition related to the location of the interface with respect to the grid point that is not always satisfied. For a slightly different kind of problem, in [START_REF] Beale | Smoothing properties of implicit finite difference methods for a diffusion equation in maximum norm[END_REF] it was proven that the numerical solution of a convection diffusion equation with an interface could allow a O(h) truncation error near the interface and still have a solution with uniform O(h 2 ) accuracy, and first differences of uniform accuracy almost O(h 2 ).

In this paper, the proof is also based on a discrete maximum principle, but differs significantly from the proof in [START_REF] Li | Maximum principle preserving schemes for interface problems with discontinuous coefficients[END_REF] because the discretization is not the same, notably due to the presence of interface unknowns, which makes the monotonicity of the matrix a crucial step in the proof, and lead to a different application of the discrete maximum principle. This result can be considered as a step toward the convergence proof of the second-order method presented in [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF]. In future works, we aim to adapt the ideas presented here to the original method itself. The crucial point being to prove that a discrete maximum principle can be applied to the discretization matrix, two alternatives could be explored:

• One could prove that the discretization matrix for the original second-order method is monotone.

Because the discretization matrix is not diagonally dominant, one would probably need to combine adequatly some elliptic inequalities for the nodes near the interface into the expression of the discrete normal derivative. It may also be necessary to modify the stencil of the flux, but still maintaining its second-order accuracy.

• One could also use the technique presented in [START_REF] Brandt | Generalized local maximum principles for finite-difference operators[END_REF], where non-monotone finite-difference methods are proven to satisfy a generalized local maximum principle, still leading to a convergence result.

Finally, we want to emphasize that the technique that we have used to obtain the bounds on the coefficients of the inverse matrix could also be used to prove the convergence of numerical methods for other numerical methods for elliptic problems, for instance without singular source terms but with discontinuous diffusion coefficients. Indeed, to our knowledge, classical estimates of the discrete Green functions were obtained mainly for smooth diffusion coefficients k, like for instance in [START_REF] Bramble | Pointwise bounds for discrete Green's functions[END_REF].

Numerical study

In this section we provide numerical results only for the first-order method in two dimensions, as the second-order method has already been validated in two-dimensions in [START_REF] Cisternino | A parallel second order cartesian method for elliptic interface problems[END_REF]. This numerical study is not meant to assess performances of the method but simply to corroborate the analysis that we have performed.

In the following we consider a square domain Ω consisting in the union of two subdomains Ω 1 and Ω 2 separated by an interface Σ. We impose exact Dirichlet boundary conditions on the outer boundary of Ω.

Numerical study of the discrete Green functions

Here we study numerically the amplitude in L ∞ -norm, of the different sums of discrete Green functions estimated in subsections 3.3.1 and 3.3.2. We consider an elliptical interface Σ defined as: The amplitude in L ∞ -norm of these sums is presented on Table 1. We observe the same behaviour as the estimates ( 20), ( 19) and ( 35 8.371 ×10 -3 400 1.820 ×10 -1 2.952 ×10 -1

4.164 ×10 -3

Table 1: Numerical amplitude of the different groups of discrete Green functions.

Convergence study: problem 1

It is a test case appearing in [START_REF] Zhou | High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources[END_REF] (MIB method, case 3 of the tests on irregular interfaces) and [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF] (CIM, example 4). We consider an elliptical interface Σ defined as:

( x 18/27
) 2 + ( y 10/27

) 2 = 1.

The exact solution is: u(x, y) =    e x cos(y), inside Σ, 5e -x 2 -y 2 2 otherwise.

We set the diffusion coefficient k = 1 outside the interface, and k = 10 inside the interface. We observe a first-order convergence, as presented in Table 2.

Convergence study: problem 2

It is a test case studied in [START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domains[END_REF]. We consider an spherical interface Σ defined by r 2 = 1/4 with r = x 2 + y 2 . The exact solution is: u(x, y) = e x cos(y) inside Σ 0 otherwise,

The numerical results and orders of convergence are presented in Table 3. We observe again a first-order convergence.

Figure 1 :

 1 Figure 1: Geometry considered: two subdomains Ω 1 and Ω 2 separated by a complex interface Σ.

Figure 2 :

 2 Figure 2: Left: regular nodes described by circles •, irregular nodes (belonging to Ω * h ) described by bullets •, right: nodes belonging to Σ h .

Figure 3 :

 3 Figure 3: Geometrical configuration near the interface in one dimension.

Figure 4 :

 4 Figure 4: All possible stencils for the first-order flux discretization on the left side of the interface, with points involved in the discretization signaled by black circles.

Figure 5 :

 5 Figure 5: Geometrical configuration near the interface in two dimensions.

27 ) 2 = 1 .

 2721 On Figures6-8, we plot respectively the numerical values for ∑Q∈Σ h G h (P, Q), ∑ Q∈Ω h G h (P, Q) and ∑ Q∈Ω * h G h (P, Q).

  ), namely a O(1) behaviour for ∑ Q∈Σ h G h (P, Q) and ∑ Q∈Ω h G h (P, Q), and a O(h) behaviour for ∑ Q∈Ω * h G h (P, Q).

Figure 6 :

 6 Figure 6: Discrete Green functions for Σ h , from top to bottom and left to right, N =50, 100, 200, 400.

Figure 7 :

 7 Figure 7: Discrete Green functions for Ω h , from top to bottom and left to right, N =50, 100, 200, 400.

Figure 8 :Table 2 :

 82 Figure 8: Discrete Green functions for Ω * h , from top to bottom and left to right, N =50, 100, 200, 400.

Inria

RR n°8872

351, Cours de la Libération Bâtiment A 29

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Samuel Coussy

Geoffrey Boissard

Stéphane Belbèze

Dominique Guyonnet

FRENCH FEEDBACK ON THE USE OF INERT WASTE LANDFILL CRITERIA FOR MANAGING EXCAVATED SOILS

Keywords: excavated soils, inert-waste criteria, landfill

The management of excavated soils in the European Union suffers from a deficit of legislation and/or common methodology regarding the definition of soil quality objectives. As the project for a European "Soil Framework Directive" was postponed in 2014, member states that do not have national legislation regarding soil quality criteria generally refer to the inert-waste criteria defined in Annex II of the European Landfill Directive. But one difficulty with this approach is that these criteria were not derived for soil materials, but for waste materials (primarily construction and demolition; C&D waste). Their development relied on a number of working hypotheses; e.g. with respect to leaching behavior, based on the study of waste materials but not of excavated soils. This is a problem because, in France, many geological formations show natural enrichments in inorganic elements, even in areas that are not considered enriched in minerals. This is for example the case around Paris and its suburbs, where several geological formations may release natural concentrations of, e.g., sulfate and selenium. A particular issue is that the earthworks of the "Grand Paris Express" project (a new subway) will generate over 40 million tons of excavated soils from such formations. Managing excavated soils based on the inert-waste criteria of the Landfill Directive would generate excessive costs due to the landfilling of a significant portion of these excavated soils in non-hazardous waste landfills, which would lead to a progressive saturation of landfills in the Paris area. An alternative, more sustainable regarding such earthwork projects, would be to consider excavated soils as a resource rather than as a waste. For this purpose, French administration has developed new criteria for the reuse of excavated soils, based on: (i) the preservation of soil quality, (ii) the preservation of water resources and (iii) compatibility with the future use of the target site, from a public health perspective.

INTRODUCTION

Infrastructure and development projects in France generate each year over 130 million tons of excavated soils. During the next ten years, the construction of a subway for the "Grand Paris Express" (creation of a new circular line and other lines around Paris) will generate over 40 million tons of excavated soils [START_REF] Sgp ; Thierry | 3D geological modelling at urban scale and mapping of ground movement susceptibility from gypsum dissolution: The Paris example (France)[END_REF]. The bulk of these excavated soils is typically landfilled because, in accordance with the Waste Framework Directive; WFD (OJEU, 2008), these soils are considered as waste materials if they are taken away from the site of their excavation. A certain proporition is reused, but very little in comparison with landfilling, mainly because of the liabilities associated with their waste status.

Such practice occurs because the management of excavated soils in the European Union suffers from a deficit of legislation regarding the definition of soil quality objectives. As the project for a European "Soil Framework Directive" was postponed in 2014, member states that do not have national legislation with respect to soil quality criteria, typically refer to the WFD which states that: "The waste status of uncontaminated excavated soils and other naturally occurring material which are used on sites other than the one from which they were excavated should be considered in accordance with the definition of waste and the provisions on by-products or on the end of waste status under this Directive.". For the case of excavated soils that are not used on the site from which they were excavated, the demonstration of absence of contamination typically relies on reference to the inertwaste criteria of Annex II [START_REF] Ojec | Council Decision of 19 December 2002 establishing criteria and procedures for the acceptance of waste at landfills pursuant to Article 16 of and Annex II to[END_REF] of the European Landfill Directive [START_REF] Ojec | Council Directive 1999/31/EC of 26 April 1999 on the landfill of waste[END_REF]. One obvious problem with this approach is that these criteria (leaching criteria for inorganic elements and total content criteria for certain organic compounds) were not derived for natural soil materials but for waste materials (primarily construction and demolition; C&D waste). Another issue is that common anthropic contaminant such as halogenated organic compounds are not included in Annex II.

This paper starts by reviewing the logic and methodology underlying the definition of the inert-waste criteria and then illustrates difficulties encountered in France regarding their application to the management of excavated soils. Finally, the paper proposes an alternative approach for the reuse of excavated soils, which does not rely on the inert-waste criteria. This approach is currently being applied in France.

WORKING HYPOTHESES UNDERLYING THE INERT-WASTE CRITERIA

The landfill acceptance criteria defined by Annex II of the Landfill Directive [START_REF] Ojec | Council Decision of 19 December 2002 establishing criteria and procedures for the acceptance of waste at landfills pursuant to Article 16 of and Annex II to[END_REF][START_REF] Ojec | Council Directive 1999/31/EC of 26 April 1999 on the landfill of waste[END_REF] were developed in the context of the Landfill Directive Technical Adaptation Comitee (TAC; [START_REF] Hjelmar | The EU landfill directive and the Council decision on acceptance criteria for waste to be landfilled: overview and associated R&D needs[END_REF] and were based on a number of working hypotheses. As described in [START_REF] Hjelmar | Development of acceptance criteria for landfilling of waste: An approach based on impact modelling and scenario calculations[END_REF], the overall procedure consisted in assessing, for waste storage scenarios specific to each type of waste (inert, non-hazardous, hazardous), which waste leaching characteristics guaranteed compliance with water quality objectives in an aquifer, at some point located downgradient from the landfill (called a point of compliance; POC). The primary source of water quality objectives was the European drinking water Directive [START_REF] Ojec | Council Directive 98/83/EC of 3 November 1998 on the quality of water intended for human consumption[END_REF].

With respect to waste leaching behavior, it was assumed that leached concentrations followed a decay function:

C = C0 e -(L/S) κ (1)
where: C is leachate concentration as a function of L/S (mg/L); C 0 is the initial peak leachate concentration (mg/L); L/S is the liquid to solid ratio for concentration C (L/Kg);

κ is a kinetic constant describing the rate of concentration decrease (Kg/L). By integrating equation ( 1), the cumulative release of contaminant mass E (in mg/Kg), over the period of time required to reach a specific value of L/S, is obtained from:

E = (C0/κ) (1 -e -(L/S) κ ) (2)
From flow and transport model calculations, maximum values of C 0 were estimated such that water quality objectives at the POCs were always respected. Knowing C 0 and assuming κ is fixed, equation (2) yields the limit cumulative emission value for a given L/S ratio (in France a value of 10 L/Kg is generally used).

Values of the kinetic parameter κ, specific to each potential contaminant (As, Cr, Sb, etc.), were derived from experimental data measured on waste materials. For example, [START_REF] Hjelmar | Development of acceptance criteria for landfilling of waste: An approach based on impact modelling and scenario calculations[END_REF] provide examples of κ values for chloride (0.57 kg/l), antimony (0.11 kg/l) and lead (0.27 kg/l). The larger the value, the faster the release. It was assumed that these κ values controlled leaching at the scale of the landfill. However, the proposed κ values do not come from excavated soils data, whether or not contaminated.

DIFFICULTIES REGARDING THE MANAGEMENT OF EXCAVATED SOILS IN FRANCE

Sulfate and selenium leaching behavior in the excavated soils from Paris and its suburbs

The project of a subway for the "Grand Paris Express" will generate 45 million of tons of excavated soils in 10 years [START_REF] Sgp ; Thierry | 3D geological modelling at urban scale and mapping of ground movement susceptibility from gypsum dissolution: The Paris example (France)[END_REF]. The public contracting authority that manages this project decided in a first step, to base the management of these excavated soils on the landfill acceptance criteria. As a prerequisite, this authority built a detailed database gathering analytical results relative to over 4,000 samples, coming from various geological formations that will be drilled during the subway construction. This database shows that an important proportion of natural (i.e., with no anthropic influence) geological materials exceed the inert-waste criteria for some soluble salts, metals and metalloids.

The example of sulfate and selenium in the natural formation of the Green Clays of the Paris Basin (of early Oligocene age) and of the Sables de Beauchamp [Beauchamp Sands] is shown in Figure 1. The boxplots show how naturally occurring geological formations may present leaching characteristics which do not comply with the inert-waste landfill criteria, thus illustrating difficulties in managing excavated soils on the basis of the Waste Directive and landfill criteria. The red line marks the inert-waste criteria for both sulfate and selenium. data come from the Grand Paris Express database.

It was possible to plot the data of the Beauchamp sands formation on a map of Paris and its suburbs (Figure 2). This was not possible for the data of the Green Clays formation, due to a lack of geographical coordinates. The results indicates that many natural samples exceed the inert-waste criteria for sulfate in the Beauchamp sands formation. For selenium, the samples exceed the inert-waste criteria only in some specific locations in this formation. If we look at the same data in other geological formations of Paris and its suburbs (data not shown here), we can see the same trends, with high sulfate contents widespread all over the area, and local but sometime high enrichments in selenium. The origin of sulfate is due to the presence of gypsum in a lot of geological formations under the city of Paris and its suburbs (e.g. [START_REF] Sgp ; Thierry | 3D geological modelling at urban scale and mapping of ground movement susceptibility from gypsum dissolution: The Paris example (France)[END_REF]. This is particularly the case in the so-called pre-Ludian formations composed of the Bartonian (Late Eocene) Calcaires de Saint Ouen [Saint Ouen Limestone] and Sables de Beauchamp [Beauchamp Sands] and the Lutetian (Middle Eocene) so-called "Marnes et Caillasses" ["Brackish Marl and Limestone"] formation (Figure 3). The origin of the high selenium content in the Paris Basin is natural [START_REF] Gourcy | Origine du sélénium et compréhension des processus dans les eaux du bassin Seine-Normandie[END_REF]. Selenium is related mainly to the presence of Eocene, mainly Ypresian, and Oligocene formations [START_REF] Chabart | Origine des anomalies en sélénium dans les captages d'AEP du département de la Marne. Première approche[END_REF][START_REF] Gourcy | Origine du sélénium et compréhension des processus dans les eaux du bassin Seine-Normandie[END_REF]. The Ypresian formations were deposited during various marine transgressions and regressions and are made of limestones, clays, carbonate sands and fluvial deposits that are rich in organic matter. The Oligocene formations are essentially made of sands.

DISCUSSION

The use of inert-waste criteria for the management of excavated soils in France is not an option since (i) the inert-waste criteria were not derived for the types of materials encountered in this study and (ii) many geological formations show natural enrichments in inorganic elements, even in areas that are not considered as mineralized, as can be seen around Paris for sulfates and selenium. As an alternative, it is proposed to consider excavated soils as a resource and not as a waste destined to be landfilled. Indeed, the landfills in the Paris vicinity are becoming more and more saturated, while excavated soils can be useful in road infrastructures or development projects, with reference to the circular economy.

For this purpose, a national technical committee was implemented in France in 2009, in the framework of the transposition of the European directive on waste into French law, to set up a methodology for excavated soil management. Part of this work was dedicated to the drafting of a guideance document for the reuse of excavated soils coming from contaminated sites [START_REF] Blanc | Guide de réutilisation hors site des terres excavées en technique routière et dans des projets d'aménagement[END_REF]. This guide was updated in 2017 [START_REF] Coussy | Guide de valorisation hors site des terres excavées issues de sites et sols potentiellement pollués dans des projets d'aménagement[END_REF], and there are ongoing efforts to address natural excavated soils (e.g. not contaminated) in a similar fashion. The guide deals with the reuse of these materials for road engineering, city planning and landscaping purposes. The methodology proposed in this guide considers compliance with three important environmental criteria for the safe reuse of excavated soils on a target site:

 Maintaning soil quality of the target site;  Protecting the water resource;  Ensuring compatibility with the future use of the target site from a public health perspective.

To meet these criteria, a tiered 3-level approach (Table 1) was proposed to classify and evaluate excavated soils for reuse purposes, from a country-scale (based on generic guideline values) to a target site-scale (based on the target site background values). This approach is based on both the pseudototal and leachate content of the potential harmful substances in excavated soils. This approach is felt to be a viable alternative to the inert-waste criteria approach because:

 it takes into account all contaminants (natural and anthropic) potentially present in excavated soils. Contaminants are selected on the basis of a minimum list of substances and a historical and documentary study;  it considers the natural background values of the target site, to allow excavated soils coming from geological formations enriched in specific inorganic elements to be reused in a target site of similar geochemical background;  it includes guideline values based on health risk assessment;  the protection of the water resource is ensured by using generic guideline values derived from inert-waste criteria or by a simplified hydrogeological study. This methodology will become mandatory in France because a new end-of-waste ordinance is in preparation for excavated soils. End-of-waste status will be acquired if the above methodology is applied.

CONCLUSIONS

The choice made in France to consider excavated soils as wastes, when they are brought outside their site of origin, have led over the years to absurd (and costly) consequences. The first one is that contractors often prefer to landfill huge amounts of (natural) excavated soils coming from earthwork projects, to be sure that they cannot be accused of causing any sort of environmental degradation (which is generally very unlikely for soils coming from natural geological formations without any sort of mineralization). The second one is that the inert-waste criteria are often used to assess the potential impact of excavated soils on the water resources, while these criteria were derived for waste materials and not for excavated soils. Consequently, many natural formations in France are not "inert" in the sense of Annex II of the European Landfill Directive (1999/31/EC), and are landfilled in non-hazardous waste landfills, which are not designed for receiving large amounts of excavated soils. The only way to tackle the problem and to change for a more sustainable approach regarding the management of earthwork projects is to consider excavated soils as a resource and not as a waste. For this purpose, the new methodology of reuse of excavated soils in France will be an important part of the solution.

Figure 1 .

 1 Figure 1. Boxplots of sulfate (A,C) and selenium (B,D) concentrations in leachates (in mg/kg) from the natural (i.e., not influenced by anthropic activity) formation of the Green Clays (A,B) and Beauchamp Sands (C,D) of the Paris Basin (respectively 122 and 654 samples for each formation leached according to the standard NF EN 12457-2, with a L/S of 10 L/Kg). The blue and green dashes correspond respectively to the upper whisker and the median.The red line marks the inert-waste criteria for both sulfate and selenium. data come from the Grand Paris Express database.

Figure 2 .

 2 Figure 2. Maps of the data from the Beauchamp sands formation for sulfate and selenium. Red triangles indicate values above inert-waste criteria and black dots values below those criteria (1000 mg/kg for sulfate and 0.1 mg/kg for selenium).

Figure 3 .

 3 Figure 3. Simplified stratigraphic column of Paris (from Thierry et al., 2009).

Table 1 .

 1 Principles of the reuse of excavated soils according to a 3-level approach.

	Reuse conditions	A. Maintaining soil	B. Protection of the	C. Compatibility with
		quality of the target site	water resource	the future use of the
				target site (from a
				public health
				perspective)
	Level 1 : Country scale	Generic guideline values	Generic guideline values	Generic guideline values
		(from national natural	(from national natural	(from natural background
		background values)	background values)	values)
	Level

2 : Local urban scale

		Regional or urban	Generic guideline values	Generic guideline values
		guideline values (from	on leachate or simplified	(from generic health risk
		regional or urban	hydrogeological study	assessment)
		background values)
	Level 3 : target site	Background values of the	Generic guideline values	Generic guideline values
	scale	target site	on leachate or simplified	(from generic health risk
			hydrogeological study	assessment) or specific
				health risk assessment

AKNOWLEDGEMENTS

The authors would like to acknowledge Kévin Colard, from Bouygues Travaux Publics, for providing the database from the Grand Paris Express.