R Émi

Vincent Perrier
email: perrier@math.u-bordeaux1.fr

ASYMPTOTIC EXPANSION OF A MULTISCALE NUMERICAL SCHEME FOR COMPRESSIBLE MULTIPHASE FLOW

The simulation of compressible multiphase problems is a difficult task for modelization and mathematical reasons. Here, thanks to a probabilistic multiscale interpretation of multiphase flows, we construct a numerical scheme that provides a solution to these difficulties. Three types of terms can be identified in the scheme in addition to the temporal term. One is a conservative term, the second one plays the role of a nonconservative term that is related to interfacial quantities, and the last one is a relaxation term that is associated with acoustic phenomena.

The key feature of the scheme is that it is locally conservative, contrarily to many other schemes devoted to compressible multiphase problems. In many physical situations, it is reasonable to assume that the relaxation is instantaneous. We present an asymptotic expansion of the scheme that keeps the local conservation properties of the original scheme. The asymptotic expansion relies on the understanding of an equilibrium variety. Its structure depends, in principle, on the Riemann solver. We show that it is not the case for several standard solvers, and hence this variety is characterized by the local pressure and velocity of the flow. Several numerical test cases are presented in order to demonstrate the potential of this technique.

1. Introduction. In [START_REF] Saurel | Some models and methods for compressible multifluid and multiphase flows[END_REF], we have considered the following model for compressible multiphase flows:

(1.1)

∂α (k) ∂t + u I • ∇α (k) = µ(P (k) -P ( k) ), ∂ α (k) ρ (k) ∂t + div α (k) ρ (k) u (k) = 0, ∂ α (k) ρ (k) u (k) ∂t + div α (k) ρ (k) u (k) ⊗ u (k) + ∇(α (k) P (k) ) = P I ∇α (k) + λ(u ( k) -u (k) ), ∂ α (k) ρ (k) E (k)
∂t + div α (k) (ρ (k) E (k) + P (k) )u (k) = P I u I • ∇α (k) -µP I (P (k) -P ( k) ) + λu I (u ( k) -u (k) ).

In (1.1), P (k) is the pressure of phase Σ k , and the total energy is related to the pressure by E (k) = ε (k) (P (k) , ρ (k) ) + 1 2 u (k) 2 where the internal energy is a concave function of its arguments. We define k = 1 (resp., = 2) when k = 2 (resp., k = 1). The interface velocity and pressure u I and P I are modeled, and this is one of the difficult points; see, e.g., Drew and Passman [START_REF] Drew | Theory of Multicomponent Fluids[END_REF]. The variables µ and λ are relaxation parameters that depend on the interfacial area.

With respect to standard approaches, the seven equation model (1.1) has several original features and properties. First, (1.1) is hyperbolic whatever the choice of the interface parameters. This is in contrast with models that are only conditionally hyperbolic such as in [START_REF] Toumi | A weak formulation of Roe's approximate Riemann solver[END_REF][START_REF] Evje | Hybrid flux-splitting schemes for a two-phase flow model[END_REF]. From the numerical point of view, the interesting consequence is that one may adapt modern standard techniques for the approximation of hyperbolic problems. Another interesting feature is the occurrence of the relaxation parameters µ and λ that can be interpreted as a way to introduce the local flow topology.

However, the model (1.1) leads to several difficult issues:

• A modelization problem. How can we define the interface velocity and pressure? From a rigorous point of view, the situation is rather clear only when λ → +∞ and µ → +∞, in which case we have pressure and velocity equilibrium. This corresponds to interface problems or problems where the phases are intimately mixed. • A mathematical problem. The system cannot be cast into conservative form, which induces difficulties for the definition of a shock and for the derivation of a numerical scheme. This point is in contradiction with the initial goal, the use of standard hyperbolic techniques for the simulation of multiphase flow problems. In the limit λ, µ → +∞, while µ/λ remains bounded, Murrone and Guillard [START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF] have obtained by asymptotic expansion the following model:

∂α (k) ∂t + u • ∇α (k) = ρ ( k) a ( k) 2 -ρ (k) a (k) 2 ρ (k) a (k) 2 α (k) + ρ ( k) a ( k) 2 α ( k) div u, ∂ α (k) ρ (k) ∂t + div α (k) ρ (k) u = 0, ∂ ρu ∂t + div ρu ⊗ u + ∇P = 0, ∂ ρE ∂t + div (ρE + P )u = 0, (1.2)
where ρ and E are the mixture density and the mixture energy defined by ρ = α (1) ρ (1) + α (2) ρ (2) , E = α (1) ρ (1) E (1) + α (2) ρ (2) E (2) ρ and a (k) is the isentropic sound velocity of each phase defined by

a (k) 2 = ∂P (k) ∂ρ (k) s (k)
.

This system was obtained by developing the different pressures and velocities of (1.1) around the equilibrium variables P and u. In what follows, we often denote this model as the "five equations model."

This system can be shown to be unconditionally hyperbolic. Unfortunately, a nonconservative product appears one more time in the volume of fluid equation. This term can be interpreted as follows: when a shock enters a mixture zone, the volume fraction must change because the shock speed is different in the two phases. This term is proportional to the divergence of the velocity and does not play any role in pure interface problems. In contrast, this term does not vanish in general across a shock wave, and it induces a change in volume. However, the difficulty, from a mathematical point of view, is that this term has no precise meaning in the case where u and

ρ ( k) a ( k) 2 -ρ (k) a (k) 2 ρ (k) a (k) 2 α (k) + ρ ( k) a ( k) 2 α ( k)
are simultaneously discontinuous. In that case, it is difficult to derive a numerical scheme.

From Murrone and Guillard analysis, it is clear that something has to be added in the volume of fluid equation to take into account the volume changes across nonlinear waves. In the present contribution, we propose to answer this problem by starting with the scheme of [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] and then to derive an asymptotic expansion from the discrete scheme, instead of the system of PDE (1.1). The scheme of [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF], which is well adapted to approximate multiphase problems, relies on a multiscale representation of the flow. We recall the fact that (1.1) is obtained from an ensemble average of realization; for any realization we have well-localized interfaces between phases, with possibly different scales for bubbles. From the numerical point of view, we write a scheme for each realization, using standard assumptions in the finite volume context. This scheme makes use of one's favorite Riemann solver. Here we consider Godunov, HLL type, acoustic, and relaxation solvers. We take into account the various scales by dividing cells into subcells being filled by any of the two phases. Then we make an ensemble average of the different schemes. This avoids defining precisely interface variables (even though a limit model can be derived), and this also solves the problem of nonconservative products (since the scheme for each realization is obtained from conservation principles). We can identify terms that correspond to relaxation terms as in (1.1); they are the consequence of acoustic effects. The last step, writing the scheme for (1.2), is obtained by doing an asymptotic expansion of the scheme of [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF].

This paper is organized as follows. We first recall in detail the results of [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] with a particular emphasis on the relaxation terms. This is the topic of section 2. The structure of the equilibrium variety for several numerical schemes is studied in section 3. Section 4 is devoted to the derivation of the asymptotic scheme of (1.2). Numerical illustrations are provided in section 5.

2.

A numerical scheme for the seven equations model. The derivation of the scheme uses the ensemble averaging ideas of Drew and Passman [START_REF] Drew | Theory of Multicomponent Fluids[END_REF] combined with the discretization principle introduced by Godunov. In [START_REF] Drew | Theory of Multicomponent Fluids[END_REF], the nonconservative and the modelization terms appear after an ensemble average. The idea of [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] is to begin the derivation of the numerical scheme before the averaging procedure. We describe the scheme in the case of the Godunov solver. All this can easily be adapted to other solvers. We also note that the different terms of (2.3) can be approximated with different solvers; the only constraint is conservation.

At each time step, the flow is described in each computational cell by the average

W j = (α (1)
j , α (1) j ρ
(1)

j , α (1) j ρ
(1)

j u (1) j , α (1)
j ρ (1) j E
(1)

j , α (2)
j , α (2) j ρ
(2)

j , α (2) j ρ
(2)

j u (2) j , α (2)
j ρ (2) j E
(2) j ). We consider a family of random subdivision of the cell

C j = [x j-1/2 , x j+1/2 ] = k [ξ k , ξ k+1 ].
In each of the subcells ]ξ k , ξ k+1 [ , we randomly set the flow variables (see Figure 1)

U (1) j = (ρ (1) j , ρ (1) j u (1) j , ρ (1) j E (1) j ) or U (2) j = (ρ (2) j , ρ (2) j u (2) j , ρ (2) j E (2) j ). ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ £ ¡£ ¡£ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¥ ¡¥ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § § ¡ § ¡ § ¡ § ¡¡¡¡ ¡¡¡ ¡¡¡ ¡¡¡ ¡¡¡ ¡¡¡ ¡¡¡ ¡¡¡ ¡¡¡ ¡¡¡ ¡¡¡ ¡¡¡ ¡¡¡ ¡¡¨© ¡© © ¡© © ¡© © ¡© © ¡© © ¡© © ¡© © ¡© © ¡© © ¡© © ¡© © ¡© © ¡© © ¡© ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ C i D D C C A B x i-1/2 x i+1/2
t + s t Fig. 1. Evolution of each phase after a random subdivision of the cell C i between time t and t + s.

The only constraint on the random subdivision is that it must be consistent with the volume fraction of each phase. Therefore, if X is the indicator function of Σ 1 , the average E(

x j+1/2
x j-1/2 Xdx) must precisely be α where σ is the local interface velocity. We denote by λ max the maximum wave speed in the internal Riemann problems. Under the CFL condition

(j) 1 (x j+1/2 -x j-1/2 ). In the cell C i , the phase Σ 1 satisfies Ci X ∂U ∂t + ∂F ∂x dx dt = 0. Then (2.1) N (ω)-1 j=0 ξj+1 ξj X ∂U ∂t dx + N (ω)-1 j=0
|λ max | s ∆ξ ≤ 1 2 ,
we can integrate (2.1) to get (2.2)

N (ω)-1 l=1 t+s t ξ l+1 ξ l X ∂U ∂t dx dt + N (ω)-1 l=1 t+s t ξ l+1 ξ l X ∂F ∂x dx dt = 0.
Equation (2.2) can be split into three parts: interaction with the left cell (2.3a), interaction with the right cell (2.3c), and internal interactions (2.3b):

(2.3a)

AD D X ∂U ∂t + ∂F ∂x dx dt (2.3b) + N (ω)-1 l=2 t+s t ξ l+1 +s σ(U l i ,U l+1 i ) ξ l +s σ(U l-1 i ,U l i ) X ∂U ∂t + ∂F ∂x dx ds (2.3c) + CBB X ∂U ∂t + ∂F ∂x dx dt = 0.
We study the three terms of (2.3).

• Boundary terms. For the term (2.3a), a straightforward calculation provides

AD D X ∂U ∂t + ∂F ∂x dx dt = AD D ∂XU ∂t + ∂XF ∂x dx dt - AD D U ∂X ∂t + F (U) ∂X ∂x dx dt = x i-1 2 +sσ + (U + i-1 ,U - i ) x i-1 2 X (x, t + s) U(x, t + s)dx -sX x i-1 2 , t + F U i-1 2 + sF lag U + i-1 , U - i [X] j=0 .
Doing the same for the right term (2.3c), we find

BCC X ∂U ∂t + ∂F ∂x dx dt = x i+ 1 2 +sσ -(U + i ,U - i+1 ) x i+ 1 2 X(x, t + s)U(x, t + s)dx -sX x i+ 1 2 , t + F U i+ 1 2 + sF lag U + i , U - i+1 [X] j=N (ω) .
• Internal terms (2.3b). A similar algebra provides

t+s t ξ l+1 +s σ(U l i ,U l+1 i ) ξ l +s σ(U l-1 i ,U l i ) X ∂U ∂t + ∂F ∂x dx ds = s [X] j F lag U j i , U j+1 i + [X] j-1 F lag U j-1 i , U j i + ξj+1+sσ(U j i ,U j+1 i ) ξj +sσ(U j-1 i ,U j i ) X (x, t + s) U(x, t + s)dx - ξj+1 ξj X(x, t)U(x, t)dx.
By summing the different terms and taking the limit when s → 0, we obtain the semidiscrete scheme (2.4)

∂ ∂t   1 ∆x x i+ 1 2 x i-1 2 X(x, t)U(x, t)dx   + 1 ∆x X x i+ 1 2 , t + F U i+ 1 2 -X x i-1 2 , t + F U i-1 2 = 1 ∆x N (ω) j=1 [X] j F lag U j i , U j-1 i -[X] j-1 F lag U j-1 i , U j-2 i + 1 ∆x [X] 0 F lag U + i-1 , U - i + [X] N (ω) F lag U + i , U - i+1 .
We may assume that two adjacent subcells contains different phases, 1 so that the sum in (2.4) can be rewritten as

N (ω)-1 j=1 [X] j F lag U j i , U j-1 i -[X] j-1 F lag U j-1 i , U j-2 i = ±N (ω) i F lag U (2) i , U (1) i
-F lag U (1) i , U (2) i
.

2.1. Averaging procedure. It remains to take the mathematical expectancy of the semidiscrete scheme (2.4). For the first term, we have

E   ∂ ∂t   1 ∆x x i+ 1 2 x i-1 2 X(x, t)U(x, t)dx     = ∂(α (1) i U (1)
i ) ∂t .

Then the scheme can be rewritten as

(2.5) ∂(α (1)
i U

(1)

i ) ∂t + 1 ∆x E X x i+ 1 2 , t + F U i+ 1 2 -E X x i-1 2 , t + F U i-1 2 = λ i F ag U (2) i , U (1) i
-F ag U (1)
i , U

i + 1 ∆x E ([X] 0 ) F ag U + i-1 , U - i + E [X] N (ω) F ag U + i , U - i+1 . (2)
Here we have introduced the notation

λ i = E N (ω) i ∆x ,
which can be interpreted as the number of interfaces per cell. 2

1 Otherwise, if two neighboring cells contain the same phase, we gather them: this is not an extra assumption. 2 In the multidimensional case, more than one relaxation parameter can be introduced. Typically, the same arguments show that there is one relaxation parameter associated with the pressure, two for the velocity in two dimensions, and three in three dimensions. These relaxation parameters depend on thermodynamical parameters, such as the acoustic impedance, but also on the interfacial area between the faces. More details can be found in [START_REF] Abgrall | A scheme for compressible two-phases flows and interface problems[END_REF].

Table 1

Flux indicator for the Eulerian flux where we denote β

(l,p) i+ 1 2 = sign " σ " U l i , U p i+1
"" .

Flow patterns

Left and right states Flux indicator

Σ 1 -Σ 2 U (1) i , U (2) i+1 " β (1,2) i+ 1 2 « + Σ 1 -Σ 1 U (1) i , U (1) i+1 1 Σ 2 -Σ 1 U (2) i , U (1) i+1
" -β

(2,1)

i+ 1 2 « + Σ 2 -Σ 2 U (2) i , U (2)
i U

(2)

i ) ∂t + 1 ∆x E X x i+ 1 2 , t + F U i+ 1 2 -E X x i-1 2 , t + F U i-1 2 = -λ i F ag U (2) i , U (1) i
-F ag U (1) i , U (2)
i - 1 ∆x E ([X] 0 ) F ag U + i-1 , U - i + E [X] N (ω) F ag U + i , U - i+1 .
In (2.5) and (2.6), the expected fluxes E(X(

x i± 1 2 , t + )F(U i+ 1 2 )), E ([X] 0 ) F ag U + i-1 , U - i , and E([X] N (ω) )F ag U + i , U - i+1
are evaluated by a close examination of each possible case.

In the case of E(X(x i± 1 2 , t + )F(U i+ 1 2

)), and specializing on the right interface of the cell C i , four cases may happen:

• either fluid Σ 1 is on the left and fluid Σ 1 is on the right • or fluid Σ 1 is on the left and fluid Σ 2 is on the right • or fluid Σ 2 is on the left and fluid Σ 2 is on the right • or fluid Σ 2 is on the left and fluid Σ 1 is on the right.

Depending on which case occurs, we may or may not take into account the fluxes. Four fluxes are evaluated (and summarized in Table 1); they are flagged when needed following the procedure we describe now. Of course, a similar situation happens on the left interface. Then if we define

P i+ 1 2 (Σ 1 , Σ 1 ) = P X x - i+ 1 2 = 1 and X x + i+ 1 2 = 1 , P i+ 1 2 (Σ 2 , Σ 2 ) = P X x - i+ 1 2 = 0 and X x + i+ 1 2 = 0 , P i+ 1 2 (Σ 1 , Σ 2 ) = P X x - i+ 1 2 = 1 and X x + i+ 1 2 = 0 , P i+ 1 2 (Σ 2 , Σ 1 ) = P X x - i+ 1 2 = 0 and X x + i+ 1 2 = 1 ,
Table 2 Flux indicator for the Lagrangian flux where we denote β

(l,p) i+ 1 2 = sign " σ " U l i , U p i+1 "" . Flow patterns Lagrangian flux Flux indicator Σ 1 -Σ 2 F ag " U (1) i , U (2) i+1 " " β (1,2) i+ 1 2 « - Σ 1 -Σ 1 F ag " U (1)
i , U

" 0

Σ 2 -Σ 1 F ag " U (2)
i , U

" -

" -β (2,1) i+ 1 2 « - Σ 2 -Σ 2 F ag " U (2)
i , U

" 0

we get

E X x i+ 1 2 , t + n F U i+ 1 2 = P i+ 1 2 (Σ 1 , Σ 1 ) F U (1)
i , U

(1) i+1

+ P i+ 1 2 (Σ 1 , Σ 2 ) β (1,2) i+ 1 2 + F U (1)
i , U

(2) i+1

+ P i+ 1 2 (Σ 2 , Σ 1 ) -β (2,1) i+ 1 2 + F U (2) i , U (1)
i+1 .

Similarly, in the case of the Lagrangian term of the right interface, the four cases are summarized in Table 2; a similar situation occurs at the left interface. We get for the right interface

E [X] N (ω) F ag U N (ω) i , U - i+1 = P i+ 1 2 (Σ 1 , Σ 2 ) β (1,2) i+ 1 2 - F ag U (1) i , U (2) i+1
-P i+ 1 2 (Σ 2 , Σ 1 ) β (2,1) i+ 1 2 - F ag U (2) i , U (1) i+1
and for the left one

E ([X] 0 ) F ag U + i-1 , U 0 i = -P i-1 2 (Σ 1 , Σ 2 ) β (1,2) i-1 2 + F ag U (1)
i , U

(2) i+1

+ P i-1 2 (Σ 2 , Σ 1 ) β (2,1) i-1 2 + F ag U (2) i-1 , U (1) i
.

These expressions are given for the phase Σ 1 . Similar expressions can be given for Σ 2 : the variable X is changed into 1 -X, and the phase index for Σ 1 (resp., Σ 2 ) is changed into the one for Σ 2 (resp., Σ 1 ).

A numerical scheme.

As detailed in [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF], a natural choice for the probability is

P i+ 1 2 (Σ 1 , Σ 1 ) = min α (1) i , α (1)
i+1 ,

P i+ 1 2 (Σ 2 , Σ 2 ) = min α (2) i , α (2)
i+1 ,

P i+ 1 2 (Σ 1 , Σ 2 ) = max α (1) i -α (1)
i+1 , 0 ,

P i+ 1 2 (Σ 2 , Σ 1 ) = max α (2) i -α (2) i+1 , 0 .
This leads to the following numerical scheme:

∂ α (1) U (1) ∂t + E (XF) i+ 1 2 -E (XF) i-1 2 ∆x (2.7) = 1 ∆x E F ag ∂X ∂x i,bound + E F ag ∂X ∂x i,relax with E (XF) i+ 1 2 = max α (1) i -α (1) i+1 , 0 β (1,2) i+ 1 2 + F U (1)
i , U

+ min α

i , α

i+1 F U

i , U

-max α

(2)

i -α (2) i+1 , 0 β (2,1) i+ 1 2 - F U (2) i , U (1)
i+1 ,

E F ag ∂X ∂x i,bound = max α (1) i -α (1) i+1 , 0 β (1,2) i+ 1 2 - F ag U (1) i , U (2) i+1
-max α

(2) i -α (2) i+1 , 0 β (2,1) i+ 1 2 - F ag U (2) i , U (1) i+1
-max α

(1)

i-1 -α (1) i , 0 β (1,2) i-1 2 + F ag U (1) i-1 , U (2)
i + max α (2) i-1 -α (2) i , 0 β (2,1) i-1 2 + F ag U (2) i-1 , U (1) i
, and

E F ag ∂X ∂x i,relax = λ i F ag U (2) i , U (1) i
-F ag U (1)
i , U

.

Setting F = 0 and U = 1 leads to the following numerical semidiscrete scheme for the volume fraction:

(2.8)

dα (1) i dt = 1 ∆x max α (1)
i -α

(1)

i+1 , 0 β (1,2) i+ 1 2 - u U (1) i , U (2) i+1
-max α

(2)

i -α (2) i+1 , 0 β (2,1) i+ 1 2 - u U (2) i , U (1) i+1
-max α

i-1 -α

i , 0 β

(1,2) i-1 2 + u U (1) i-1 , U (2)
i + max α (2) i-1 -α (2) i , 0 β (2,1) i-1 2 + u U (2) i-1 , U (1) i
+ λ i u U (2) i , U (1) i
-u U (1) i , U (2) i
.

In practical applications, the discretization is achieved via a splitting method. First, we integrate (with some abuse of notation3 ) (2.9)

W (1) i n+1/2 -W (1) i n ∆t + E (XF ) i+1/2 -E (XF ) i-1/2 ∆x = E F ag ∂X ∂x i,bound
∆x and then a relaxation step (2.10)

W (1) i n+1 -W (1) i n+1/2 ∆t = E F ag ∂X ∂x i,relax .
The flux are computed at time t n . The relaxation step is carried out implicitly.

2.3.

Extension to second order. The scheme (2.9)-(2.10) makes use of several terms (the conservative terms, the boundary terms, and the relaxation terms) that do not play the same role in the derivation. The "conservative" terms are standard, so that any standard high order extension will do. The "relaxation" and "boundary" terms have the same status: they are the agglomeration of fluxes between different phases for which the high order procedure is not immediate.

In order to construct a high order extension of the scheme (2.9)-(2.10), the procedure described in [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] is considered:

1. We extrapolate the physical variables (α ( ) , u ( ) , P ( ) ), = 1, 2, using the MUSCL method. Here, the minmod or van Leer-van Albada limiters are used. 2. We subdivide the cell C i into N regular subcells (N arbitrary) denoted by [y l , y l+1 ] (1 ≤ l ≤ N -1) and average in each subcell the reconstructed variables. 3. Each subcell is randomly subdivided as above. 4. Then we gather the contribution of each subsub-cell and let N → +∞. In the following, U ( ) i+1/2,l (resp., U ( ) i-1/2,r ) represents the variable U ( ) at the cell interface x i+1/2 (resp., x i-1/2 ) on the left (resp., right) after MUSCL extrapolation. If f is any of the extrapolated variables, δ i f represents its limited slope in C i . The jump in α ( ) at y l+1/2 is ∆α ( ) (y l+1/2 ).

The nonconservative terms

E F ag ∂X ∂x i = E F ag ∂X ∂x i,bound + E F ag ∂X ∂x i,relax sum up to (2.11) -β (1,2) i+1/2 + P i+1/2 (Σ 1 , Σ 2 )F ag (U (1) i+1/2,l , U (2)
i+1/2,r ) + β

(2,1) i+1/2

+ P i+1/2 (Σ 2 , Σ 1 )F ag (U (2) i+1/2,l , U (1) i+1/2,r ) -β (1,2) i-1/2 -P i-1/2 (Σ 1 , Σ 2 )F ag (U (1) i-1/2,l , U (2) i-1/2,r ) + β (2,1) i-1/2 -P i-1/2 (Σ 2 , Σ 1 )F ag (U (2) i-1/2,l , U (1) i-1/2,r ) + N -1 l=1 max 0, ∆α (1) (y l+1/2 ) F ag (U (2) 2 (y l+1/2 ), U (1) 1 (y l+1/2 )) - N -1 l=1 max 0, ∆α (2) (y l+1/2 )) F ag (U (2) 2 (y l+1/2 ), U (1) 1 (y l+1/2 )) + N -1 l=1 λ(y l+1/2 ) F ag (U (2) 2 (y l+1/2 ), U (1) 1 (y l+1/2 )) - N -1 l=1 F ag (U (1) 1 (y l+1/2 ), U (2)
2 (y l+1/2 )) , where ∆α (1) (y l+1/2 ) = δ i α (1) (y l+1 -y l ).

We do not need any additional β terms because we account for all the internal cells in ]x i-1/2 , x i+1/2 [ . The third and fourth terms of the above formula converge, when N → +∞, to (2.12)

x i+1/2 x i-1/2
F ag (U (2) (y), U (1) (y))dy max 0, δ i α (1) -

x i+1/2 x i-1/2
F ag (U (1) (y), U (2) (y))dy max 0, δ i α (2) .

A similar formula is obtained for the phase Σ 2 . The integrals in (2.12) are evaluated via the midpoint formula (since second order accuracy is sought).

Similarly, the relaxation terms corresponding to a linear reconstruction of the data can be approximated, thanks to the same interpretation in terms of Riemann sums and to the midpoint rule by (2.13)

λ i F ag (U (2)
i , U

i ) -F ag (U

i , U

i ) .

We finally obtain a second order approximation of the nonconservative terms by

∆x E F ag ∂X ∂x = -β (1,2) i+1/2 + P i+1/2 (Σ 1 , Σ 2 )F ag (U (1) i+1/2,l , U (2) i+1/2,r ) + β (2,1) i+1/2
+ P i+1/2 (Σ 2 , Σ 1 )F ag (U (2) i+1/2,l , U (1)
i+1/2,r ) -β (1,2) i-1/2 -P i-1/2 (Σ 1 , Σ 2 )F ag (U (1) i-1/2,l , U (2) i-1/2,r ) + β (2,1) i-1/2 -P i-1/2 (Σ 2 , Σ 1 )F ag (U (2) i-1/2,l , U (1) i-1/2,r ) + λ i F ag (U (2) i , U (1) i ) -F ag (U (1) i , U (2) i ) + max 0, δα (1) i F ag (U (2) i , U (1) i ) -max 0, δα (2) i F ag (U (1) i , U (2)
i ), (2.14) where δα

(1) i = α (1) i+1/2,l -α (1)
i-1/2,r and δα

(2) i = α (2) i+1/2,l -α (2)
i-1/2,r are the limited slope of α (1) and α (2) in the cell C i .

Hence, the second order extension of the scheme can be geometrically interpreted by adding an additional interface inside the cell C i ; this is the meaning of last two terms of (2.14).

2.4. Extension to other solvers. In (2.7), the numerical flux F that we use at the microscopic level is obtained thanks to an approximate Riemann solver for which it is possible to define a contact speed. The contact speed between the left state U L and the right state U R is denoted by σ(U L , U R ). We also denote by U ± as in Figure 2 for the Godunov solver. This permits us to define the Lagrangian flux

F ag (U L , U R ) = F (U + LR )-σ(U L , U R )U + LR = F (U - LR )-σ(U L , U R )U - LR . Last, we define β (p,q) i+1/2 = sign(σ(U (p) i , U (q) i+1 )).
All the calculations have been performed for the Godunov scheme but can be extended to more general fluxes. The key ingredients of the derivation are, besides the randomization, average procedures, and estimation of the various coefficients, the use of a Riemann solver for which it is possible to define a contact discontinuity speed σ(U, V ). This property is needed because we must define a Lagrangian flux,

F ag (U, V ) = F (U, V ) -σ(U, V )U .
This Lagrangian flux has to be consistent, as well as the base flux F (U, V ).

Even more, the choice of the base flux F and the Lagrangian flux F ag may be independent: we do not really need a relation of the type F ag = F * -σU ; see [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] for more details.

3. Asymptotic model for the numerical scheme. To get the numerical scheme for the five equations model (1.2), we propose following a similar technique as in [START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF], where the continuous model was derived via an asymptotic expansion of (1.1). The very difference is that we will do it at a discrete level. If we set ε i = 1 λi , the discrete scheme for the seven equations model can be formally rewritten as

(3.1) ∂W ∂t + G ∆x = R(W ) ε i ,
and we assume that ε i is large. The first step to get the asymptotic scheme is to determine the equilibrium set

{W such that R(W ) = 0} .
From the previous section, this set is defined as (3.2) V = {W = (α (1) , α (1) U (1) , α (2) , α (2) U (2) ) such that

F ag (W 1 , W 2 ) = F ag (W 2 , W 1 )}.
Then the limit model is obtained by developing the variables around the equilibrium set. We assume that V can be parametrized by a mapping M : u → M (u) (we prove this fact later) and that R(W ) has an expansion of the form

R(W ) = R(M (u)) + εR (M (u))V + o(ε). Thus (3.1) becomes dM u ∂u ∂t + G ∆x = R (M (u))V + o(1).
To get the reduced model, it remains to multiply by the projection P onto the kernel of R (M (u)) in the direction of dM u . This gives the model

(3.3) P dM u ∂u ∂t + P G ∆x = 0.
The aim of this section is to show that for several classical solvers, P does not depend on the solver. We have to show two properties:

1. V is equal to (3.4) {(α (1) , α (1) U (1) , α (2) , α (2) U (2) ) such that u 1 = u 2 and P 1 = P 2 } as in the continuous case. 2. The Jacobian matrix of F ag (U (1) , U (2) ) -F ag (U (2) , U (1) ) in primitive variables is

(3.5) DR(M (u)) =             0 0 0 -λ 0 0 0 λ 0 0 0 0 0 0 0 0 0 0 µ 0 0 0 -µ 0 0 0 µ u λ P 0 0 -µ u -λ P 0 0 0 λ 0 0 0 -λ 0 0 0 0 0 0 0 0 0 0 -µ 0 0 0 µ 0 0 0 -µ u -λ P 0 0 µ u λ P            
as in the continuous case. In the following, we denote by a subscript 12 (resp., 21) the intermediate states obtained by solving an elementary Riemann problem before averaging with left and right states given by the phase Σ 1 (resp., Σ 2 ) on the left and Σ 2 (resp., Σ 1 ) on the right.

3.1. Case of the acoustic solver. Given two states, left and right, described by the physical variables, (ρ L , u L , P L ) for the left state and (ρ R , u R , P R ) for the right state, the Lagrangian flux for the acoustic solver is defined as (0, P , P u ) with the pressure and the velocity of the intermediate state( ) given by (see, e.g., [START_REF] Richtmyer | Difference Methods for Initial-Value Problems[END_REF])

(3.6) P = Z R P L + Z L P R + Z L Z R (u L -u R ) Z R + Z L , u = Z R u R + Z L u L + P L -P R Z R + Z L . In (3.6), Z L,R = ρ L,R a L,R
is the acoustic impedance. Therefore, using the above mentioned conventions, and using the fact that within the cell we consider the left/right states are defined from the conservative variables U (1) for Σ 1 and U (2) for Σ 2 , we have

P 12 = Z 1 P 1 + Z 2 P 2 + Z 1 Z 2 (u 2 -u 1 ) Z 2 + Z 1 , P 21 = Z 1 P 1 + Z 2 P 2 + Z 1 Z 2 (u 1 -u 2 ) Z 1 + Z 2 , u 12 = Z 1 u 1 + Z 2 u 2 + P 2 -P 1 Z 1 + Z 2 , u 21 = Z 1 u 1 + Z 2 u 2 + P 1 -P 2 Z 1 + Z 2 so that u 12 -u 21 = 2 Z 1 + Z 2 (P 1 -P 2 ) , P 12 -P 21 = 2Z 1 Z 2 Z 1 + Z 2 (u 1 -u 2 ) .
We immediately see that if u 12 -u 21 = 0 and P 12 -P 21 = 0, we have u 1 = u 2 and

P 1 = P 2 .
Let us evaluate now the Jacobian matrix of F ag (U (1) , U (2) ) -F ag (U (2) , U (1) ) with respect to the primitive variable on V. We have F ag (U (1) , U (2) ) -F ag (U (2) , U (1) ) = (0, P 12 -P 21 , P 12 u 12 -P 21 u 21 ), and we have to evaluate the partial derivative of this vector with respect to v = ρ, u, and P . Using the fact that the considered state lies in V, a straightforward calculation gives (with u = u 1 = u 2 and P = P 1 = P 2 )

∂ P 12 -P 21 ∂v = 2 Z 1 + Z 2 ∂P 1 ∂v - ∂P 2 ∂v , ∂ P 12 u 12 -P 21 u 21 ∂v = 2Z 1 Z 2 Z 1 + Z 2 u ∂P 1 ∂v - ∂P 2 ∂v + 2 Z 1 + Z 2 P ∂u 1 ∂v - ∂u 2 ∂v
so that we get (3.5) with

λ = 2 Z 1 + Z 2 and µ = 2Z 1 Z 2 Z 1 + Z 2 .
3.2. Case of the exact solver. For the exact solver, the Lagrangian flux can be written as (0, P , P u ), where P and u are the intermediate velocity and pressure when the Riemann problem is solved exactly. Equation (3.4) can be rewritten as (0, P 12 , P 12 u 12 ) = (0, P 21 , P 21 u 21 ) so that u = u 21 = u 12 and P = P 21 = P 12 .

Thus in the (u, P ) variables, the point (u , P ) lies on the 1-wave curve coming from the state (u 1 , P 1 ) and on the 3-wave curve coming from the state (u 1 , P 1 ). Therefore u = u 1 and P = P 1 , because the 1-wave curve is decreasing and the 3-wave curve is increasing: this is a consequence of the convexity of the equation of state [7, page 89].

The same argument holds for the state 2, so that u = u 2 and P = P 2 so that u 1 = u 2 and P 1 = P 2 . Now we estimate the Jacobian matrix. Indeed, we find exactly the same Jacobian matrix as for the acoustic solver because the relations (3.6) are obtained by a linearization in the (u, P ) variables of the true Riemann invariants at (u = u 1 = u 2 , P = P 1 = P 2 ).

3.3. Case of the HLLC solver. The HLLC solver is based on the choice of two wave velocities S L and S R . The intermediate states are calculated by writing the Rankine-Hugoniot relations across these waves; see [11, page 293] for more details. With this solver, the intermediates states are given by

U K = ρ K S K -u K S K -u     1 u E K + (u -u K ) u + P K ρ K (S K -u K )     with K = L or R. The Lagrangian flux is then F ag =   0 ρ K (u -u K ) S L -u L + P L ρ K (u -u K ) S L -u L + P L u   .
Thus F ag 12 and F ag 21 are given by the following expressions:

(3.7a)

F ag 12 =   0 ρ 1 (u 12 -u 1 ) S L 12 -u 1 + P 1 ρ 1 (u 12 -u 1 ) S L 12 -u 1 + P 1 u 12   (3.7b) =   0 ρ 2 (u 12 -u 2 ) S R 12 -u 2 + P 2 ρ 2 (u 12 -u 2 ) S R 12 -u 2 + P 2 u 12   , (3.7c)
F ag 21 =   0 ρ 1 (u 21 -u 1 ) S R 21 -u 1 + P 1 ρ 1 (u 21 -u 1 ) S R 21 -u 1 + P 1 u 21   (3.7d) =   0 ρ 2 (u 21 -u 2 ) S L 21 -u 2 + P 2 ρ 2 (u 21 -u 2 ) S L 21 -u 2 + P 2 u 21   .
The intermediate velocity is given by

(3.8) u = P R -P L + ρ L u L (S L -u L ) -ρ R u R (S R -u R ) ρ L (S L -u L ) -ρ R (S R -u R ) .
For more convenience, we denote by

Z L = ρ L (u L -S L ), Z R = ρ R (S R -u R ).
Then (3.8) becomes

u = P L -P R + Z R u R + Z L u L Z R + Z L .
Here we have used the same notation convention as in section 3.1 for the velocity and pressure. The subscripts 1 and 2 in Z are intended to highlight the fact that the definition of Z depends not only on the left and right state through ρ and u but also on the ordering of the waves through the speed S. For example, Z 12 1 is evaluated with U (1) on the left and U (2) on the right and corresponds to the leftmost wave S. Thus we have

u 12 -u 21 = P 2 -P 1 + Z 21 1 u 1 + Z 21 2 u 2 Z 21 1 + Z 21 2 - P 1 -P 2 + Z 12 1 u 1 + Z 12 2 u 2 Z 12 1 + Z 12 2 = Z 21 1 + Z 21 2 + Z 12 1 + Z 12 2 (Z 21 1 + Z 21 2 )(Z 12 1 + Z 12 2 ) (P 2 -P 1 ) + Z 21 1 Z 12 2 -Z 21 2 Z 12 1 (Z 21 1 + Z 21 2 )(Z 12 1 + Z 12 2 ) u 1 - Z 21 1 Z 12 2 -Z 21 2 Z 12 1 (Z 21 1 + Z 21 2 )(Z 12 1 + Z 12 2 ) u 2 = A(u 1 -u 2 ) + B(P 2 -P 1 ) with (3.9a) A = Z 21 1 Z 12 2 -Z 21 2 Z 12 1 (Z 21 1 + Z 21 2 )(Z 12 1 + Z 12 2 ) , (3.9b) B = Z 21 1 + Z 21 2 + Z 12 1 + Z 12 2 (Z 21 1 + Z 21 2 )(Z 12 1 + Z 12 2 )
.

Let us examine the difference between the second component of (3.7a) and (3.7c), i.e., P 12 -P 21 . We need to evaluate u -u L and u -u R . From (3.8), we have

u -u R = P L -P R + Z L (u L -u R ) Z R + Z L , u -u L = P L -P R + Z R (u R -u L ) Z R + Z L .
This leads to

P 12 -P 21 = Z 12 2 (u 12 -u 2 ) + Z 21 2 (u 21 -u 2 ) = Z 12 2 P 1 -P 2 + Z 12 1 (u 1 -u 2 ) Z 12 1 + Z 12 2 + Z 21 2 P 2 -P 1 + Z 21 1 (u 1 -u 2 ) Z 21 1 + Z 21 2 = Z 21 1 Z 12 2 -Z 21 2 Z 12 1 (Z 21 1 + Z 21 2 )(Z 12 1 + Z 12 2 ) (P 1 -P 2 ) + Z 12 2 Z 12 1 Z 12 1 + Z 12 2 + Z 21 2 Z 21 1 Z 21 1 + Z 21 2 (u 1 -u 2 ) = C(u 1 -u 2 ) -A(P 2 -P 1 ) with (3.10) C = Z 12 2 Z 12 1 Z 12 1 + Z 12 2 + Z 21 2 Z 21 1 Z 21 1 + Z 21 2 .
Clearly, the system u 12 -u 21 = 0, P 12 -P 21 = 0 is equivalent to

A(u 1 -u 2 ) + B(P 2 -P 1 ) = 0, C(u 1 -u 2 ) -A(P 2 -P 1 ) = 0,
which has a determinant -A 2 -BC.

For the system to be invertible, it is sufficient that B and C be strictly positive. This is true in particular if all the Z are strictly positive. We show later in this section and for several choices of speeds S L,R that Z > 0, and thus the system is invertible.

Assume the system is invertible: the equilibrium set is (3.4). It remains to evaluate the Jacobian matrix along the equilibrium set V. This is done as in section 3.1. For example, if we compute the derivative of P 12 -P 21 with respect to v = ρ, u, or P , we have

∂(P 12 -P 21 ) ∂v = ∂A ∂v (P 1 -P 2 ) + ∂B ∂v (u 1 -u 2 ) + A ∂(P 1 -P 2 ) ∂v + B ∂(u 1 -u 2 ) ∂v .
If the variables belong to V, we have

∂(P 12 -P 21 ) ∂v = A ∂(P 1 -P 2 ) ∂v + B ∂(u 1 -u 2 ) ∂v .
Therefore, if A = 0 in the equilibrium set, we get the result for P 12 -P 21 . The same result holds for u 12 -u 21 . Thus, the only thing left is to show that A is zero along (3.4). The remainder of this section is devoted to showing, for several classical choices of speeds S [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF], that 1. Z > 0, 2. A defined in (3.9b) satisfies A = 0 along V.

3.3.1.

Case where S L = u La L and S R = u R + a R . We have

Z 12 1 = ρ 1 a 1 , Z 12 2 = ρ 2 a 2 , Z 21 1 = ρ 1 a 1 , Z 21 2 = ρ 2 a 2 ,
and we see immediately that the Z's are all strictly positive and that, directly from (3.9b), A = 0 along V.

Case where

S L = min (u L -a L , u R -a R ) and S R = max (u R + a R , u L + a L ).
In this case, it is obvious that the Z's are strictly positive so that the equilibrium set is (3.4). Along this equilibrium set, we denote by ā the maximum sound speed and u the velocity of each fluid. We have

Z 12 1 = ρ 1 ā, Z 12 2 = ρ 2 ā, Z 21 1 = ρ 1 ā, Z 21 2 = ρ 2 ā;
therefore, A = 0 along the equilibrium set.

3.3.3. Case where S L = ũã and S R = ũ + ã, where x is the Roe average 4 of x. The equation which gives the intermediate velocity no longer holds because we are not sure that the denominator is nonzero. Nevertheless, we still can use the equations (3.7). The equality of the second and third component of (3.7a) and (3.7b) gives immediately u 12 = u 21 , and we denote by u this velocity. Moreover, by abstracting the second component of (3.7c) and the second component of (3.7a), and by doing the same algebra with (3.7b) and (3.7d), we get

(3.11) (u -u 2 ) S L 12 -S R 21 = 0, (u -u 1 ) S L 21 -S R 12 = 0.
Here the wave velocities are symmetric in u 1 and u 2 , so that S L 12 = S L 21 = S L and S R 21 = S R 12 = S R . We have S L -S R = -2ã = 0, which leads to

u 1 = u 2 = u .
Last, the equality of the second component of (3.7a) and (3.7b) provides P 1 = P 2 . Thus, the equilibrium set is once more (3.4). Along the equilibrium set, we have

Z 12 1 = ρ 1 ã, Z 12 2 = ρ 2 ã, Z 21 1 = ρ 1 ã, Z 21 2 = ρ 2 ã,
where a is the Roe average of a 1 and a 2 , so that Z's are strictly positive, and we can use the formula (3.9a) for A which gives immediately A = 0 on V.

Case of pressure-velocity-based wave speed estimates.

A different approach, that works for perfect gases only, consists in first estimating the pressure P and u and then derivating the estimates of S L and S R . If we suppose that we have estimated P and u , then we choose the following wave speeds:

S L = u L -a L q L , S R = u R + a R q R ,
where (3.12)

q K =    1 if P ≤ P K , 1 + γ + 1 2γ P P K -1 1 2 if P > P K .
Then we find

Z 12 1 = ρ 1 a 1 q 12 L , Z 12 2 = ρ 2 a 2 q 12 R , Z 21 1 = ρ 1 a 1 q 21 R , Z 21 2 = ρ 2 a 2 q 21 L . 4 x := √ ρ L x L + √ ρ R x R √ ρ L + √ ρ R .
These Z's are all strictly nonnegative, so that the equilibrium set is (3.4). On the equilibrium set, P 12 = P 21 so that with (3.12) we have

q 12 L = q 21 R and q 12 R = q 21 L ,
and A is null on the equilibrium set.

3.4. Case of the relaxation solver of [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF]. We approximate the solutions of the Euler equations with those of the following relaxation system:

(3.13)                          ∂ρ ∂t + ∂ (ρu) ∂x = 0, ∂ (ρu) ∂t + ∂ ρu 2 + π ∂x = 0, ∂ (ρE) ∂t + ∂ ((ρE + π) u) ∂x = 0, ∂π ∂t + u ∂π ∂x + a 2 ρ ∂u ∂x = P -π η ,
where a is a parameter chosen to ensure the system is dissipative when the relaxation time η → ∞ (a > max (ρc)). The solution of the Riemann problem for (3.13) can easily be computed because all the fields are linearly degenerate. In particular, we find the following formulas for the intermediate pressure and velocity:

π = π R + π L 2 + a 2 (u L -u R ) , u = π L -π R 2a + u L + u R 2 .
Thus we immediately get

P 12 -P 21 = a 12 + a 21 2 (u 1 -u 2 ) , u 12 -u 21 = 1 2 1 a 12 + 1 a 21 (P 1 -P 2 ) .
As in the case of the acoustic solver, we see immediately that the equilibrium set is (3.4). Last, the study of the Jacobian matrix with the same method as before shows that (3.5) is true with

λ = 1 2 1 a 12 + 1 a 21 and µ = a 12 + a 21 2 .
4. Derivation of the numerical scheme for the five equations model. In the previous section, we have shown that, for several classical solvers, the equilibrium variety is the same as the one in the continuous case, namely (4.1)

V = u 1 = u 2 and P 1 = P 2 .
Therefore, V is of dimension 6 and can be parametrized by the mapping M :

M :         α 1 ρ 1 u P α 2 ρ 2         →             α 1 ρ 1 u P α 2 ρ 2 u P             .
Since this parametrization is naturally written in primitive variables and not in conservative variables, the first thing to do is to provide the form of the numerical scheme for the seven equations model in primitive variables, which is equivalent to the form in conserved variables. The next step is to compute the projector P onto the kernel of R (M (u)) in the direction of dM u and then to evaluate (3.3).

Transformation into primitive form.

In order to simplify the algebra, we rewrite (3.1) as

(4.2) ∂α (k) i ∂t = V F (k) i , ∂ α (k) i ρ (k) i ∂t = M a F (k) i , ∂ α (k) i ρ (k) i u (k)
i ∂t = M o F (k) i , ∂ α (k) i ρ (k) i E (k) i ∂t = EF (k) i . Since (4.3) α (k)
i ρ (k) i ∂u (k) i ∂t = M o F (k) i -u (k) i M a F (k) i
we have

i ∂t = 1 2 u (k) i M o F (k) i + 1 2 u (k) i M o F (k) i -u (k) i M a F (k) i = 1 2 u (k) i M o F (k) i + 1 2 u (k) i M o F (k) i - u (k) i 2 2 M a F (k) i = u (k) i M o F (k) i - u (k) i 2 2 M a F (k) i . 1 2 ∂ α (k) i ρ (k) i u (k) i 2 ∂t = u (k) i 2 ∂ α (k) i ρ (k) i u (k) i ∂t + α (k) i ρ (k) i u (k) i 2 ∂u (k)
Similarly,

EF (k) i = ∂ α (k) i ρ (k) i E (k) i ∂t = ∂ α (k) i ρ (k) i ε (k) i ∂t + 1 2 ∂ α (k) i ρ (k) i u (k) i 2 ∂t and ∂ α (k) i ρ (k) i ε (k) i ∂t = EF (k) i -u (k) i M o F (k) i + u (k) i 2 2 M a F (k) i .
We develop the energy equation to get an equation on the pressure

∂ α (k) i ρ (k) i ε (k) i ∂t = ε (k) i ∂ α (k) i ρ (k) i ∂t + α (k) i ρ (k) i ∂ε (k) i ∂t = ε (k) i M a F (k) i + α (k) i ρ (k) i β (k) i ∂P (k) i ∂t + α (k) i ρ (k) i κ (k) i ∂ρ (k) i ∂t ,
where

β (k) i = ∂ε (k) i ∂P (k) i ρ (k) i and κ (k) i = ∂ε (k) i ∂ρ (k) i P (k) i .
When we combine the mass and the volume fraction equations, we obtain

(4.4) α (k) i ∂ρ (k) i ∂t = M a F (k) i -ρ (k) i V F (k)
i .

Then we have α

i ∂t = EF (k) i -u (k) i M o F (k) i + u (k) i 2 2 M a F (k) i -ε (k) i M a F (k) i -ρ (k) i κ (k) i M a F (k) i -ρ (k) i V F (k) i , that is, ∂P (k) i ∂t = EF (k) i ρ (k) i β (k) i ∂P (k) i ∂t = EF (k) i -u (k) i M o F (k) i + u (k) i 2 2 M a F (k) i -ε (k) i M a F (k) i -α (k) i ρ (k) i κ (k) i ∂ρ (k)
i M o F (k) i α (k) i ρ (k) i β (k) i +   u (k) i 2 2 -ε (k) i -ρ (k) i κ (k) i   α (k) i ρ (k) i β (k) i M a F (k) i + ρ (k) i 2 κ (k) i V F (k) i α (k) i ρ (k) i β (k) i . (k) i α (k) i ρ (k) i β (k) i - u (k)
The scheme (2.7) is equivalent to (4.4), (4.3), and (4.5).

The projection.

Formally, the seven equations numerical scheme in primitive variables is

∂U ∂t + G ∆x = T R(U ) ε j ,
where T is the linear transformation between primitive and conservative variables

T = T 1 0 0 T 2 with (4.6) T i =           1 0 0 0 - ρ i α i 1 α i 0 0 0 - u i α i 1 α i ρ i 0 κ i ρ i      1 0 0 α 1 α 2 d 0 0 - α 1 α 2 d 0 1 0 - ρ 1 α 2 d 0 0 ρ 1 α 2 d 0 0 ρ 1 α 1 α 1 ρ 1 + α 2 ρ 2 0 0 ρ 2 α 2 α 1 ρ 1 + α 2 ρ 2 0 0 0 0 α 1 ρ 2 a 2 2 d 0 0 α 2 ρ 1 a 1 2 d 0 0 0 α 1 ρ 2 d 1 0 - αρ 2 d              with d = α 1 ρ 2 (a 2 ) 2 + α 2 ρ 1 (a 1 ) 2 .
By putting back the system into the reduced conservative variables (α i , α i ρ i , ρu, ρE), 5 Note that

P k -κ k ρ k 2 β k ρ k 2 = a k 2 .
the evaluation of (3.3) provides the following scheme for the five equations model:

∂α 2 ∂t = F V 2 + α 1 α 2 α 2 ρ 1 a 2 1 + α 1 ρ 2 a 2 2     SE 2 α 2 ρ 2 β 2 - u 2 SU 2 α 2 ρ 2 β 2 + u 2 2 2 -ε 2 -ρ 2 κ 2 α 2 ρ 2 β 2 M 2 + ρ 2 2 κ 2 F V 2 α 2 ρ 2 β 2 - SE 1 α 1 ρ 1 β 1 + u 1 SU 1 α 1 ρ 1 β 1 - u 2 1 2 -ε 1 -ρ 1 κ 1 α 1 ρ 1 β 1 M 1 - ρ 2 1 κ 1 F V 1 α 1 ρ 1 β 1     , ∂ (α k ρ k ) ∂t = M k , ∂ (ρu) ∂t = SU 1 + SU 2 , ∂ (ρE) ∂t = SE 1 + SE 2 ,
where we have used the notations of (4.2). We conclude this section by noting that this applies to the first order as well as the second order scheme.

5. Numerical results. In the numerical tests, all the fluids were described by the Stiffened-gas equation of state (EOS):

ε(P, ρ) = P + γP ∞ (γ -1)ρ ,
where γ and P ∞ are constants that depend on the phase under consideration. In this paper, the fluids used in the numerical tests are the water and the air, which have coefficients γ and P ∞ given in Tables 3 and6. In some cases, one side of the tube will be filled with pure fluids, in which case the volume fraction of the fluid that is not here will be initialized with α = 10 -8 . All examples are run with a CFL number of 0.6. The CFL number is computed on the largest wave speed in the flow. The test cases are chosen so that one can see the influence of the div u terms and are ordered by increasing difficulty. The last case is a comparison with experimental data.

We have chosen to use different solvers in the experiments. Several considerations have motivated our choices: the best possible accuracy and the most possible robustness. In interface-like problems, we need to initialize one of the fluids with a very small volume fraction, typically of the order of 10 -8 . In that case, the solver "sees" the fluid as vacuum and might be very sensitive to any numerical error, especially in the case of complex and stiff wave systems. Because of that, some solvers are not suited to some cases. In the following, we have chosen to use the acoustic solver in the tests of section 5.2 (in fact all the solvers are working fine on this case) and the relaxation solver on the other ones. This choice seems to be the best compromise. 5.1. Single fluid test. In that case, the shock tube is filled with water only. At time t = 0, the left part of the tube is at rest, with a pressure of 10 9 Pa, and the right side is at rest too but with a pressure of 10 5 Pa. The separation between both sides is at x = 0.5m, and the profiles are shown at t = 150 µs. The approximated solution is computed with a 1000 point grid mesh and is compared with the analytical solution in Figure 3. The second order scheme has been run. The exact solution and the numerical one are in good agreement.

Pure interface advection.

A one meter long tube is filled on the right with air and on the left with water. The pressure is uniform, equal to 10 5 Pa, and the velocity is uniform too, equal to 1000 m.s -1 . At time t = 0, the discontinuity is localized at x = 0.5m. The exact solution consists in the advection of mass at 1000 m.s -1 , and the pressure and the velocity remain uniform. The results are shown at time t = 200 µs when the exact discontinuity is located at x = 0.7 m. The solution is computed with a 1000 cell mesh and is shown in Figure 4. The results show that the velocity and the pressure remain uniform, and the density and volume of fluids are simply advected, as expected. 5.3. Liquid-gas shock tube. Here we study the evolution of a tube filled on the left with high-pressure water (10 9 Pa) and on the right with atmospheric-pressure air (10 5 Pa). The discontinuity is initially at x = 0, 7 m. The details of the initial conditions are given in Table 4. The computed solutions of the present scheme and of [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] are compared in Figure 5.

Table 4

Details of the conditions for the liquid-gas shock tube of section 5.3. Note that these solutions are computed with the second order scheme. In Figure 6, we display the first and second order accurate solutions for this problem. Clearly, the second order results are much more accurate than the first order ones, as expected.

Two phase flow problem.

In that case, the left and the right part of the shock tube are filled with water and air at the same volume fraction of 0.5. The pressure is of 10 9 Pa on the left and of 10 5 Pa on the right. The discontinuity is initially at x = 0.5 m. The details of the initial conditions are given in Table 5. The computed solutions of the present scheme and of [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] are compared in Figure 7. The results are in good agreement, but we note some discrepancies. The main difference between the present schemes and that of [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF], as it can be seen in Figure 7, consists in the density levels between the contact and the shock wave. The location of the discontinuities as well as the extreme points of the fan coincide exactly up to numerical errors of course. Note, however, that the average densities almost coincide between the contact and the shock and do coincide elsewhere. The pressure and velocity plots superimpose up to numerical errors. In our opinion, the fact that the discontinuities move at the same speed is a good indication that we are computing the same system: our experience indicates that in the case of nonconservative systems, as here, a small modification in the approximation of nonconservative terms implies a large modification of the numerical wave speeds. These are certainly a consequence of the different techniques in time integration: in the original scheme, the relaxation terms are integrated by a splitting technique, but one has to remember that the relaxation terms originate from fluxes. On the contrary, in the present technique, no splitting technique is used. The asymptotic expansion we use has some features of a direct integration of the relaxation terms inside the solver, which is probably better in principle. These algorithmic differences seem here to have a large impact on the volume fraction, but we have no clear explanation why. In Figure 8, we have represented a zoom of the average density between the contact discontinuity and the shock wave for the first and second order versions of the present scheme and the original one. We clearly see the (small) differences. The most interesting phenomenon is that, first, the two first order solutions are very close and, second, the second order solution for the present scheme is much closer from the first order ones than is the second order version of the original scheme. From the algorithmic point of view, the original second order scheme uses a predictor corrector method, and the predicted quantities are projected onto the equilibrium manifold after the predictor and the corrector steps. In the present scheme, this is translated into the choice of the states that are needed in the computation of the projector Π defined in (4.7). In the second order version of the present scheme, Π is evaluated after each substep. Hence, there is no reason why the two sets of results should coincide exactly. It appears, however, that the

Liquid density

Gas density Fig. 5. Liquid-gas shock tube: Comparison between the solutions obtained by the original scheme of [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] (circles) and the present scheme (solid). Both solutions differ only on the densities, in zones where they have few meaning because the corresponding volume fractions are nearly zero.

present technique, based on an asymptotic expansion, is a more stable technique than a simple fractional step method. The lack of exact solution does not enable more accurate comments.

Comparison with experiment.

We make a simulation on an alloy made of epoxy and spinel for which experimental data exist. The tube is impacted with a piston, so that a shock propagates. The experimental setup is sketched in Figure 9. There exists a linear relation between the shock speed u c and the impact velocity u i : u c = a 0 + su i . Epoxy and spinel are described by the stiffened-gas EOS with the coefficient of Table 6.

Initially, both of the solids are at atmospheric pressure, and their densities are ρ 1 = 1185 kg.m -3 and ρ 2 = 3622 kg.m -3 . Their respective volume fractions are α 1 = 0.595 and α 2 = 0, 405. We measure the different variables in x = 0.1 m and in x = 0.6 m. Thanks to these measures (shown in Figure 10) we can compute the shock speed. In Figure 11, we compare our results with the experimental measures taken from [START_REF] Marsh | LASL Shock Hugoniot Data, Los Alamos Series on Dynamic Material Properties[END_REF], those of [START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF], and those obtained by the original scheme and the present scheme.

We see that our results are very close to those of the original scheme and in good agreement with the experimental results. In particular, we see that there is an excellent agreement between the results of the seven equation model with instantaneous relaxation and the present scheme. The agreement is much better than with the Murrone and Guillard scheme for which the div u terms are discretized independently because

• our results are closer from the experimental ones, 0.7 0.75 0.8 0.85 0.9 900 920 940 960 980 1000 Present scheme, first order Original scheme, first order Present scheme, second order Original scheme, second order Fig. 8. Comparison of the first and second order version of the original scheme and the present one for the averaged density. Zoom between the contact discontinuity and the shock wave for the data of Table 5.

¢¡ ¤£

¥ §¦ © § © ! " " $# % '& )( # 10 $2 & 43 "5 4% % 76 2 98 4@ % 2 BA % C6 2 98 4@ % 2 ED

Epoxy

Spinel γ = 2.94 γ = 1.62 P ∞ = 3.2 × 10 9 Pa Pa P ∞ = 141 × 10 8 Pa 6. Conclusions. We have presented a numerical scheme able to compute compressible multiphase schemes that relies on a multiscale description of the flow. Three types of terms can be identified in the scheme in addition to the temporal terms. One is a conservative term, the second one plays the role of a nonconservative term that is related to interfacial quantities, and the last one is a relaxation term that is associated with acoustic phenomena. The key feature of the scheme is that it is locally conservative, contrarily to many other schemes devoted to compressible multiphase problems. The scheme is developed for several Riemann solvers. A second order extension is also described.

In many physical situations, it is reasonable to assume that the relaxation is instantaneous: either the bubbles are very small or, more generally, the interfacial area between phases is large. We present an asymptotic expansion of the scheme that keeps the local conservation properties of the original scheme. The asymptotic expansion relies on the understanding of an equilibrium variety. Its structure depends, in principle, on the Riemann solver. We show that for several standard solvers that this is not the case, and hence this variety is characterized by the local pressure and velocity of the flow.

Several numerical test cases are presented in order to demonstrate the potential of this technique.

Fig. 2 .

 2 Fig. 2. Structure of the Riemann problem.

Fig. 3 .

 3 Fig. 3. Single fluid problem: approximated solution (circles) and exact solution (solid).

Fig. 4 .

 4 Fig. 4. Pure fluid advection: the interface is diffused and centered around x = 0.7 m. The pressure and velocity are left unchanged.

Fluid 1: air α 1 = α α 1 = 1

 111 -α ρ 1 = 50 kg.m -3 ρ 1 = 50 kg.m -3 Fluid 2: water α 2 = 1 -α α 2 = α ρ 2 = 1000 kg.m -3 ρ 2 = 1000 kg.m -3 Global variables P = 10 9 Pa P = 10 5 Pa u = 0 m.s -1 u = 0 m.s -1

Fig. 6 . 5

 65 Fig.6. Comparison between the first and second order accurate schemes.Table 5Details of the conditions for the two phase flow problem of section 5.4.

Fluid 1 :

 1 air α 1 = 0.5 α 1 = 0.5 ρ 1 = 50 kg.m -3 ρ 1 = 50 kg.m -3 Fluid 2: water α 2 = 0.5 α 2 = 0.5 ρ 2 = 1000 kg.m -3 ρ 2 = 1000 kg.m -3 Global variables P = 10 9 Pa P = 10 5 Pa u = 0 m.s -1 u = 0 m.s -1

Fig. 9 .

 9 Fig. 9. Experimental setup.

Fig. 10 .

 10 Fig. 10. Evolution of the variables at x = 0.1 m and x = 0.6 m between the initial time and t = 2. 10 -4 s.

Fig. 11 .

 11 Fig. 11. Shock velocity, comparison between several approaches.

Table 3

 3 EOS coefficients for water and air.

						Air			Water
						γ = 1.4		γ = 4.4
						P ∞ = 0 Pa		P ∞ = 6.10 8 Pa
				1100
						Computed
						exact
				1000
				900
					0	0.2	0.4	0.6	0.8	1
							position (meter)
						Density
								1e+09
													computed
													Exact
	200							8e+08
								6e+08
				Computed
	100			Exact				4e+08
								2e+08
	0	0	0.2	0.4	0.6	0.8	1	0	0	0.2	0.4	0.6	0.8	1
				position (meter)						position (m)
				Velocity					Pressure

Table 6

 6 EOS coefficients for epoxy and spinel.

This is because (2.8) is not written as a difference of flux plus a relaxation term.

Liquid density

Gas density Fig. 7. Two phase flow problem: the present scheme (circles) is compared with the original scheme of [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] (solid lines).

• there is very little difference between the results with the complete seven equation model and the reduced five equations model discretized with our technique.


Christophe Voisin

Jean-Robert Grasso

Éric Larose

François Renard

Evolution of seismic signals and slip patterns along subduction zones: insights from a

Evolution of seismic signals and slip patterns along subduction zones: insights from a friction lab scale experiment.

Christophe Voisin, Jean-Robert Grasso, Eric Larose, François Renard * Shallow dipping subduction zones are the best natural laboratories to observe and measure the whole gamut of slip patterns: Large seismic events break the upper locked portion of the subduction plate into events of magnitude up to 9.3-9.4 so far; M=6.0-7.5 slow slip events (SSE) or silent events occur down-dip the locked portion of the subduction zone 1 ' 2,9 . Deeper along these subduction zones, the oceanic plate slips continuously at the plate tectonics velocity (Figure 1). SSE and silent events are sometimes accompanied by Non Volcanic Tremors (NVT) [START_REF] Dragert | A silent slip event on the deeper Cascadia subduction interface[END_REF] . A tremor represents long duration of high-amplitude seismic signals without clear body wave arrivals [START_REF] Schwartz | Slow slip events and seismic tremor at circumpacific subduction zones[END_REF] .

The physical origin of NVT is poorly understood, being either proposed to be triggered by fluid movements [START_REF] Obara | Nonvolcanic deep tremor associated with subduction in southwest Japan[END_REF] or related to some frictional processes along the subduction interface [START_REF] Ide | Mechanism of deep low frequency earthquakes: Further evidence that deep non-volcanic tremor is generated by shear slip on the plate interface[END_REF][START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF] . Both origins are problematic. If fluids are at the origin of tremors, then we should record tremors updip of the locked zone, where the subducted lithosphere is known to be fluid-saturated. If tremors do find their origin in some frictional process, then we should record NVT not only down-dip the locked zone, but all along the subduction zone. Whatever the process at the origin of NVT, the fact is that they are mostly associated with SSE along subduction zones, suggesting a common origin to both phenomenons. The understanding of SSE is incomplete. Pressure and temperature (P & T) increases with depth are invoked to explain the changes of slip patterns down dip the subduction zone [START_REF] Schwartz | Slow slip events and seismic tremor at circumpacific subduction zones[END_REF] . Aqueous fluids released from subsequent metamorphic dehydration reactions (enhanced by P&T changes) have been proposed to play an important role in causing slip transients and NVT 14 . From the numerical point of view, P&T changes are used to establish the profile of Rate and State friction parameters a and b with depth [START_REF] Liu | Aseismic slip transients emerge spontaneously in threedimensional rate and state modeling of subduction earthquake sequences[END_REF][START_REF] Liu | Spontaneous and triggered aseismic deformation transients in a subduction fault model[END_REF] . b-a is tuned to control the stable/unstable behaviour of the subduction interface, being positive or negative respectively. The depth of the transition zone is set up according to the temperature dependence of b-a. This dependence was experimentally determined for a wet granite gouge [START_REF] Blanpied | Quantitative measure of the variation in fault rheology due to fluid-rock interactions[END_REF] . Other studies impose small but negative (b-a) values at the transition zone to reproduce the change in slip pattern. In this latter case, the conditional stability of slip is insured by either a stiffness that is greater than the critical one [START_REF] Mitsui | Slow slip events controlled by the slab dip and its lateral change along a trench[END_REF] or by a plastic increase of the critical slip distance d c [START_REF] Shibazaki | On the physical mechanism of silent slip events along the deeper part of the seismogenic zone[END_REF] . In all of these models, the subduction interface is imposed to slip at a constant tectonic velocity below the transition zone. Numerical models of friction applied to subduction zones teach us that the transition between earthquakes and SSE occurs down dip the locked zone, where the slip becomes conditionally stable. However, the same situation for stability exists also up-dip of the locked zone [START_REF] Scholz | Earthquakes and friction laws[END_REF] where no SSE has been observed so far. The process(es) at the origin of SSE must explain this difference in the slip pattern above and below the locked zone, respectively. The occurrence of NVT is not investigated by these numerical models.

In the following, we investigate the influence of cumulative slip and ageing of the contact interface on the occurrence of SSE and NVT. The cumulative slip along subduction zones is among the largest on Earth: indeed, a portion of the subducted lithosphere slips over hundreds of kilometres. During this slip, two main processes drive the change in the contact interface: 1/ the sediment input in the subduction factory that constitutes a highly deformable contact interface; 2 / P & T increases, known to drive changes in rock deformation from elastic brittle to plastic creep behaviours. We have devised a friction experiment (Figure 1) that is aimed at reproducing the main features of a natural subduction zone. To recover P & T effects on rock deformation, we allow for plastic deformation to be effective over the time scale of our friction experiment by using a salt slider [START_REF] Voisin | Long term friction: From stick-slip to stable sliding[END_REF][START_REF] Shimamoto | LABORATORY FRICTION EXPERIMENTS AND NATURAL EARTHQUAKES -AN ARGUMENT FOR LONG-TERM TESTS[END_REF] . The slider, an analogue to a portion of the subduction plate, slides against a flat glass surface which represents the overriding continental lithosphere. The slider is pushed at a constant velocity, equivalent to the tectonic velocity. We follow the displacement of the slider through time, as if we were moving along with the buried lithosphere (Figure 1). We record the Acoustic Emission (AE) induced by the slip of the slider.

The slip pattern of the slider is observed to evolve continuously with the cumulative displacement. Acoustics of friction at the scale of 1 square centimetre is generally related to the surface roughness [START_REF] Akay | Acoustics of friction[END_REF] . Recording the AE at high frequency during the friction experiment provides a way to resolve the processes of ageing going on at a lower scale of friction: the contact asperity scale.

In zone NVT appear as a local recollection of the unstable frictional behaviour occurring during the slippage of aseismic slow or silent events down dip the subduction zones [START_REF] Shelly | Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip[END_REF] . The occurrence of slow and silent events, possibly twined to NVT is causally linked to the amount of cumulated slip undergone by the deformable contact interface. This explains why neither SSE nor NVT are recorded up-dip of the locked zone, where the accumulated slip is too small.

The different classes of slip events (earthquakes, slow events, silent events, continuous sliding) are based on observations, either seismic or geodetic. To press further the comparison between our experimental results and the subduction zone events, we estimate slip velocity ratios for the four classes of events (Figure 2). The normalized slip velocity V* (V event / V driving ) is different for subduction zones and for the analogue experiment, as far as fast unstable slips are concerned. Because of a low sampling frequency, we cannot be conclusive on this apparent misfit. Nonetheless, V* scales the same way for subduction and analogue events. Our experimental results suggest that earthquakes, slow and silent events pertain to a continuum of slip patterns driven by the ageing of the subduction interface with cumulated slip. Estimates of the cumulated slip required to observe slow or silent events on actual subduction zones are given by the depth of these events and the dip angles of subduction zones. Compiling data from Kato Correspondence and requests for data should be addressed to CV (cvoisin@obs.ujf-grenoble.fr) Lower line: the normalized slip velocity V* in the four cases (V*=V/V driving ). V driving is cm/year for subduction zones, and m/year for the analogue experiment. In the experiment, V* is about 300, 60, 10, 1, for stick slip, smooth stick slip, smooth oscillations and continuous sliding respectively. For subduction zones, with typical velocities of m/s, cm/s and µm/s, V* is about 10 9 , 10 7 , 10, 1 for earthquakes, slow events, silent events, continuous sliding respectively. The duration of each event is defined by the time period during which the slip velocity is two times greater than the driving velocity. This duration is lower than 0.05 s for the regular stick slip; about 5 s for smooth stick slip; 25 s for smooth oscillations; not defined for continuous sliding. Table 1. A comparison of slip ratio, slip duration ratio and slip velocity ratio (V*) between earthquake and slow or silent slip events for the subduction zone and for the analogue experiment. Slip ratio is slip earthquake / slip event ; Slip duration ratio is duration earthquake / duration event ; Slip velocity ratio is given by V* = V event / V driving . V driving is cm/year for subduction zones, and m/year for the analogue experiment. V event is defined as the ratio of slip over duration. The low temporal resolution of the experiment (20 Hz) precludes to explore high speed velocities, which may explain the discrepancies between V* for fast slips. In spite of this, the concordance of all ratios between the subduction zone and the analogue experiment supports the control of slip pattern by the cumulative displacement and the deformation of the contact interface.

  Laboratoire de Géophysique Interne et Tectonophysique, CNRS, Université Joseph Fourier, Maison des Géosciences, BP 53, 38041 Grenoble cedex 9, FRANCE. * also at Physics of Geological Processes, Oslo, Norway Continuous GPS and broadband seismic monitoring have revealed a variety of disparate slip patterns especially in shallow dipping subduction zones, among which regular earthquakes, slow slip events and silent quakes 1,2 . Slow slip events are sometimes accompanied by Non Volcanic Tremors (NVT), which origin remains unclear 3 , either related to fluid migration or to friction. The present understanding of the whole menagerie of slip patterns is based upon numerical simulations imposing ad hoc values of the rate and state parameters a and b 4-6 derived from the temperature dependence of a and b of a wet granite gouge 7 . Here we investigate the influence of the cumulative slip on the frictional and acoustic patterns of a lab scale subduction zone. Shallow loud earthquakes (stick-slip events), medium depth slow, deeper silent quakes (smooth sliding oscillations) and deepest steady-state creep (continuous sliding) are reproduced by the ageing of contact interface with cumulative displacement 8 . The Acoustic Emission evolves with cumulative displacement and interface ageing, following a trend from strong impulsive events, similar to earthquake seismic signals, to a collection of smaller amplitude and longer duration signals, similar to Non Volcanic Tremors. NVT emerge as the recollection of the local unstable behaviour of the contact interface globally evolving towards the stable sliding regime.

  Figure 2 presents four exerts of this continuous change, taken for different stages of cumulative slip. We observe a change in the frictional behaviour of the slider, from unstable (stick-slip like) to stable slip, passing through quasi stable smooth oscillations. The beginning of the experiment is characterized by large stick-slip events (Figure2, snapshot 1) which amplitude gradually decreases as the slider enters a smooth oscillations sliding regime (snapshots 2-3). The slip pattern of the slider eventually reaches the continuous sliding regime at the imposed driving velocity (snapshot 4). If the amplitude of slip events is decreasing with the cumulative displacement, the duration of slip episodes increases with the cumulative slip. The stick-slip episodes are characterized by a fast stress release (< 0.05 second, limited by the sampling frequency), while the duration of the smooth sliding oscillations is about 30 seconds.In a previous study, we reported on the change in the contact interface geometry with cumulative displacement[START_REF] Voisin | Long term friction: From stick-slip to stable sliding[END_REF] . The change in slip pattern of the slider is accompanied by the development of a striation pattern in the direction of sliding. Here we record the AE associated with the slider slippage and track any change in the acoustic recordings in relation with the change in the contact interface. Figure 3 presents three AE records taken after different increasing cumulative slips. The amplitude of AE signals is decreasing with the cumulative displacement. On the contrary the duration of AE signals increases with the cumulative displacement. The envelope of AE signals is evolving as well. The three examples show a transition from strong impulsive signals to a succession of smaller and longer signals, made of distinct acoustic events.

  our experiments, the progressive evolution of the acoustic signals is related to the contact interface changes and the development of the striation pattern. The records of AE throughout the experiments reveal a progressive change from impulsive events to smaller and longer signals made of numerous acoustic events. It reproduces the in-situ change from earthquake seismic signals to NVT signals observed in subduction zones. Our experiments argue for slow slip events (SSE) and non volcanic tremors (NVT) being two manifestations of a single friction process, seen at two different scales. These experimental results shed a new light on both the slip and seismic patterns along the subduction zones. In our experiments, tremor like signals are emitted by unstable slip at the contact asperity scale (shear failure asperities) occurring during the slider slow slip event. Ageing of the contact interface with cumulated displacement do occur along the subduction interface, driving the change from regular to slow to silent earthquakes. Subduction

[ 2003 ] 4 ;

 20034 we derive a value for cumulated slip of 70 to 230 km, with a mean value of 130 km. Using the appropriated tectonic plate velocities 20 , the necessary time for these slow and silent slip events to emerge is in the 1.5-4 10 6 yr range, with a 3 10 6 yr mean value.Allowing deformation of the contact interface to interfere with friction reveals a continuum of slip patterns. Earthquakes, slow events, silent quakes and continuous sliding appear as different aspects of a holistic process. The ageing of the contact interface with cumulated displacement provides a global framework to capture the occurrence of the different slip patterns and seismic signals along subduction zones. Considering the cumulative displacement as a sine qua non condition for the occurrence of SSE and NVT explains the absence of these latter above the locked zone. Our experimental results are consistent with and rationalize a posteriori: (i) the modelling of aseismic slip transients by a decrease in b-a[START_REF] Liu | Aseismic slip transients emerge spontaneously in threedimensional rate and state modeling of subduction earthquake sequences[END_REF] and an increase in d c (ii) the hypothesis that silent slip and NVT pertain to one and unique phenomenon of friction; (iii) the hypothesis that NVT are local reminiscence of frictional instabilities in these aseismic slip transients[START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF] .

Figure 1 .

 1 Figure 1. (A) Scheme of a subduction zone geometry with the relative location of large earthquakes (EQ), slow and silent events (SSE), and continuous stable sliding (CS). The oceanic plate (pink) is moving at constant imposed velocity (the tectonic loading rate) and is subducted under a continental lithospheric plate (yellow). The experimental analogue to the subduction zone is a salt slider block (light grey) pushed at constant velocity over a planar glass surface. Pressure and temperature effects are mechanically introduced by allowing creep deformation to be efficient at the time scale of the friction experiment. Different frictional slip patterns of the slider are observed: Stick slip (SS) followed by quasi-stable smooth oscillations (SO) and by continuous sliding (CS). (B) Overview of the friction experiment. The salt sample (1 cm 2 surface area) is housed in a nickel iron alloy (Invar©) plate with a low thermal expansion coefficient to minimize thermal perturbations. The salt slider is carefully sliced from a single crystal of NaCl to avoid crack formation. The frictional interface is then roughened with a sand paper (Struers #320 grit) to ensure an initial roughness. A constant continuous velocity of 0.11 µm.s -1 is imposed to the plate. The slider is subjected to a constant normal load of 0.24 MPa and is in contact with a flat glass surface. Five displacement transducers (d 1-5 ) record the plate movement in horizontal (d 1-2 ) and vertical (d 3-5 ) directions. A force gauge records the shear force exerted onto the slider. An accelerometer (reference #4518 from Bruel & Kjaer) is glued on the upper surface of the slider and records the acoustic emission in the 20 Hz-70 kHz frequency range.

Figure 2 .

 2 Figure 2. The continuous variation of the slip pattern of the salt slider. Upper line, from left to right: regular stick slip, smooth stick slip, smooth oscillation, continuous sliding.

Figure 3 .

 3 Figure 3. Three examples of Acoustic Emission records taken after different increasing cumulative slips (0.5 cm; 1.5 cm; 2.5 cm respectively). The horizontal scale is the time in seconds. The time t=0s is the triggering time. The vertical scale represents the vertical acceleration of the slider, given in m.s -2 . The amplitude of AE signals is decreasing with the cumulative displacement (2, 1, 0.05 m.s -2 respectively). On the contrary the duration of AE signals increases with the cumulative displacement (~10 -4 s; ~10 -3 s; ~10 -2 s respectively). The envelope of AE signals is evolving as well. The three examples show a transition from strong impulsive signals to a succession of smaller and longer signals, made of distinct acoustic events. Note that no AE could be recorded during the continuous sliding stage because of the triggered mode of recording.

Figure 4 .

 4 Figure 4. The change in the Acoustic Emission of the slider. Three consecutive experiments increase the cumulative slip by 1 cm each: PA006 (0 to 1cm); PA007 (from 1 to 2cm); PA008 (from 2 to 3cm). For each experiment, all acoustic events are recorded and stacked in order to increase the signal to noise ratio. The results are presented in the form of a spectrogram. The low frequency range of the signals recorded below 10 kHz is associated with the movement of the slider as a block. The high frequency range (around 50 kHz) is associated with the fast slip of the contact asperities. The time t=0s is the triggering time. The time window is focused between -0.02s and +0.02s. Some differences appear in the form of signals recorded. During the first centimetre of slip, the high frequency signal is recorded as a single large acoustic event. During the second centimetre of slip, the acoustic signals have smaller amplitudes but remain in the form of a unique acoustic event. Some weak signals begin to appear in the 0-0.01s time window. During the third centimetre of slip, the amplitude of the records is weaker than previously recorded. The high frequency part of signal is made of a collection of different events occurring at various times, before and after the triggering time.

Figure 1 Figure 2 Figure 3 Figure 4

 1234 Figure 1

Table 1

 1

			subduction zone		analogue experiment
	slip	earthquake ~1 m	slow event ~0.5 m	silent quake ~0.1 m	stick slip ~30 µm	smooth stick slip ~10 µm	smooth oscillation ~5 µm
	slip ratio	1	2	10	1	3	6
	duration	10 s	hours-days	days-years	<0.05 s	~5 s	~25 s
	duration ratio	1	3.10 -3 -10 -4 10 -3 -10 -5	1	< 10 -2	< 2.10 -3
	V *	earthquake ~10 9	slow event ~10 7	silent quake ~2-10	stick-slip 6.10 3	smooth stick-slip ~20	smooth oscillation ~2

Kato, N., A possible model for large preseismic slip on a deeper extension of a seismic rupture plane. Earth and Planetary Science Letters

(1-2), 17 (2003).