Tanguy Rivoal

Padé type approximants of Hurwitz zeta function ζ(4, x)

HAL is a

Introduction

The Hurwitz zeta function is defined by

ζ(s, x) = ∞ n=0 1 (n + x) s
for any s, x such that Re(s) > 1 and x ∈ R \ Z ≤0 . For any fixed integer s, ζ(s, x) is meromorphic in C \ Z ≤0 , with poles of order s at each non-negative integer. For any fixed ε > 0, ζ(s, x) has an asymptotic expansion when x → ∞ in the angular sector | arg(x)| < π -ε:

ζ(s, x) ∼ x 1-s s -1 + ∞ k=1 (s) k-1 k! B k x k+s-1 ,
where (B k ) k≥0 is the sequence of Bernoulli numbers.

In this paper, we address the following problem: Given two integers m, n ≥ 0, determine two polynomials A(x) and B(x) ∈ Q[x], of degree ≤ n, such that

A(x)ζ(4, x) + B(x) = O 1 x m+1 (1)
as x → ∞ in any angular sector | arg(x)| < π -ε.

Stated in this form, this is an analytic problem. However, using the asymptotic expansion of ζ(s, x) at x = ∞, Eq. ( 1) can also be interpreted as a Padé type problem at x = ∞ for the formal series

2x -3 + ∞ k=0 (k + 2)(k + 3)B k+1 x -k-4 .
See [START_REF]Simultaneous polynomial approximations of the Lerch function[END_REF]Sec. 2] for details. This Padé problem amounts to solving a linear system with 2n + 2 indeterminates (the polynomial coefficients) and m + n + 1 equations (the vanishing conditions): provided m ≤ n, this system has at least one non-zero solution. The case m = n corresponds to the usual diagonal Padé approximation. The explicit polynomials obtained below are automatically solutions of the associated Padé type problem. Our main result is the explicit determination of a non-zero analytic solution of (1) when m ≤ n/2 (essentially). Unicity of the solution is obviously not guaranteed.

Theorem 1. For any integer n ≥ 0, consider the following Padé type problem: determine two polynomials Q 0,n (x) and

Q 2,n (x) ∈ Q[x], of degree ≤ 4n, such that S n (x) := 3Q 0,n (x)ζ(4, x) + Q 2,n (x) = O 1 x 2n+3 (2)
Problem (2) admits the following solution:

S n (x) = - ∞ k=0 ∂ ∂k k + x + n 2 (k + 1) 2 n (k + 2x) 2 n (k + x) 4 n+1 . (3)
The series converges for all x ∈ C \ {0, -1, -2, -3, . . .}.

Moreover, for the series on the right-hand side of (3), we have

Q 0,n (x) = n j=0 ∂ ∂ε n 2 -j + ε (x + j -n -ε) 2 n (x -j + ε) 2 n (1 -ε) 4 j (1 + ε) 4 n-j |ε=0 and Q 2,n (x) = - 1 6 n j=0 ∂ ∂ε 3 n 2 -j + ε (x + j -n -ε) 2 n (x -j + ε) 2 n (1 -ε) 4 j (1 + ε) 4 n-j j-1 k=0 1 (x + k -ε) 2 |ε=0 .
Diagonal Padé approximants are known for ζ(2, x) and ζ(3, x): the formulas are given in [START_REF] Prévost | A new proof of the irrationality of ζ(2) and ζ(3) using Padé approximants[END_REF] and [START_REF]Simultaneous polynomial approximations of the Lerch function[END_REF]Theorem 2]. However diagonal Padé approximants are not explicitely known for ζ(4, x) and Theorem 1 offers a weaker alternative. The polynomial Q 0,n (x) can also be written more symbolically in the form

Q 0,n (x) = n j=0 ∂ ∂j n 2 -j n j 4 x + j -1 n 2 x + n -j -1 n 2 .
We now let

Q 1,n (x) = - 1 6 n j=0 ∂ ∂ε 3 n 2 -j + ε (x + j -n -ε) 2 n (x -j + ε) 2 n (1 -ε) 4 j (1 + ε) 4 n-j j-1 k=0 1 x + k -ε |ε=0
which is also a polynomial of degree ≤ 4n. Then, our proof will also show that

R n (x) := Q 0,n (x)ζ(3, x) + Q 1,n (x) = ∞ k=0 k + x + n 2 (k + 1) 2 n (k + 2x) 2 n (k + x) 4 n+1 . As x → ∞, k + x + n 2 (k + 1) 2 n (k + 2x) 2 n (k + x) 4 n+1 ∼ 4 n (k + 1) 2 n x 2n+3 (4)
and this suggests that R n (x) = O( 1 x 2n+3 ) as S n (x). However, this is false and in fact one can not prove anything better than R n (x) = O( 1x 2 ). Therefore, the property

S n (x) = O( 1 x 2n+3
) is a non-trivial property, which is not a simple consequence of (4).

We give the proof of Theorem 1 in Section 2 while we make connections with other results in the literature in Section 3. I warmly thank Pierre Bel for his careful reading of a previous version of this paper, and for pointing out that R n (x) is only O( 1x 2 ) and not O( 1 x 2n+3 ).

Proof of Theorem 1

We follow the method used in [START_REF]Simultaneous polynomial approximations of the Lerch function[END_REF] and split the proof in three parts. We also include the case of ζ(3, x) in the first part of the proof. We define the rational function

ρ(t) = t + x + n 2 (t + 1) 2 n (t + 2x) 2 n (t + x) 4 n+1
.

By partial fraction expansion, we have

ρ(t) = 4 s=1 n j=0 E j,n,s (x) (t + x + j) s with E j,n,s (x) = 1 (4 -s)! ∂ ∂t 4-s ρ(t) (t + x + j) 4 |k=-j-x .
Exchanging summations, we thus get

∞ k=0 k + x + n 2 (k + 1) 2 n (k + 2x) 2 n (k + x) 4 n+1 = 4 s=2 n j=0 E j,n,s (x) ζ(s, x) - 4 s=1 n j=0 j-1 k=0 E j,n,s (x) (k + x) s . (5)
Here we must explain how we have disposed of the divergent series ∞ k=0 1 k+x in (5), i.e. why the first sum over s does not start at s = 1. The series on the left-hand side of (5) being convergent, this forces to assign the value

-n j=0 E j,n,1 (x)) j-1 k=0 1 k+x to the divergent expression ( n j=0 E j,n,1 (x)) ∞ k=0 1
k+x+j . Indeed we have n j=0 E j,n,1 (x) = 0 because this is the sum over the residues at all the finite poles of ρ(k), hence also equal to minus its residue at infinity, which is zero. Formally, one should introduce the Lerch series ∞ k=0 z k (k+x) s with |z| < 1 and eventually to let z → 1; see [START_REF]Simultaneous polynomial approximations of the Lerch function[END_REF]Sec. 3.2] for details.

We now observe that ρ

(k) = -ρ(-k -2x -n). Since ρ(-k -2x -n) = 4 s=1 n j=0 E j,n,s (x) (-k -x -n + j) s = 4 s=1 n j=0 (-1) s E n-j,n,s (x) (k + x + j) s ,
we then deduce that E n-j,n,s (x) = (-1) s+1 E j,n,s (x). Therefore

n j=0 E j,n,s (x) = (-1) s+1 n j=0 E j,n,s (x)
which is thus equal to 0 for s = 2 and s = 4, and consequently, the first sum in ( 5) is for s = 3 only:

∞ k=0 k + x + n 2 (k + 1) 2 n (k + 2x) 2 n (k + x) 4 n+1 = n j=0 E j,n,3 (x) ζ(3, x) - 4 s=1 n j=0 j-1 k=0 E j,n,s (x) (k + x) s . (6)
Similarly,

- ∞ k=0 ∂ ∂k k+x + n 2 (k + 1) 2 n (k + 2x) 2 n (k + x) 4 n+1 = - 4 s=1 n j=0 E j,n,s (x) ∞ k=0 ∂ ∂k 1 (k + x + j) s = 4 s=1 n j=0 E j,n,s (x) sζ(s + 1, x) - 4 s=1 n j=0 j-1 k=0 sE j,n,s (x) (k + x) s+1 so that - ∞ k=0 ∂ ∂k k + x + n 2 (k + 1) 2 n (k + 2x) 2 n (k + x) 4 n+1 = n j=0 E j,n,3 (x) 3ζ(4, x) - 4 s=1 n j=0 j-1 k=0 sE j,n,s (x) (k + x) s+1 . (7)

The coefficients are polynomials of degree ≤ 4n

We set Q 0,n (x) := n j=0 E j,n,3 (x) and

Q 1,n (x) := - 4 s=1 n j=0 j-1 k=0 E j,n,s (x) (k + x) s , Q 2,n (x) := - 4 s=1 n j=0 j-1 k=0 sE j,n,s (x) (k + x) s+1
so that the right-hand sides of ( 6) and ( 7) are respectively equal to

Q 0,n (x)ζ(3, x) + Q 1,n (x) and 3Q 0,n (x)ζ(4, x) + Q 2,n (x).
Let us prove that for s ∈ {0, 1, 2}, the Q s,n (x) are in Q[x] and of degree ≤ 4n. We have

Q 0,n (x) = n j=0 ∂ ∂k ρ(k) (k + x + j) 4 |k=-j-x = n j=0 ∂ ∂ 4 ρ( -j -x) | =0 = n j=0 ∂ ∂ n 2 -j + (x + j -n -) 2 n (x -j + ) 2 n ( 1
-) 4 j (1 + ) 4 n-j | =0 . (8)
Eq. ( 8)

shows that Q 0,n (x) ∈ Q[x] and deg(Q 0,n ) ≤ 4n. Furthermore, by Leibniz' formula Q 1,n (x) = - 4 s=1 n j=0 1 (4 -s)! ∂ ∂ 4-s 4 ρ( -j -x) × 1 (s -1)! ∂ ∂ s-1 j-1 k=0 1 k -+ x | =0 = - 1 6 n j=0 ∂ ∂ 3 4 ρ( -j -x) j-1 k=0 1 k -+ x | =0 = - 1 6 n j=0 ∂ ∂ 3 n 2 -j + (x + j -n -) 2 n (x -j + ) 2 n ( 1
-) 4 j ( + 1) 4 n-j j-1 k=0 1 k -+ x | =0
and similarly,

Q 2,n (x) = - 4 s=1 n j=0 (-1) 4-s (4 -s)! ∂ ∂ 4-s 4 R( -j -x) × 1 (s -1)! ∂ ∂ s-1 j-1 k=0 1 (x + k -) 2 | =0 = - 1 6 n j=0 ∂ ∂ 3 4 R( -j -x) j-1 k=0 1 (k -+ x) 2 | =0 = - 1 6 n j=0 ∂ ∂ 3 n 2 -j + (x + j -n -) 2 n (x -j + ) 2 n ( 1
-) 4 j ( + 1) 4 n-j j-1 k=0 1 (x + k -) 2 | =0 . It follows that Q 1,n (x) and Q 2,n (x) are in Q[x]
and of degree ≤ 4n, because for all j ∈ {0, . . . , n}, k ∈ {0, . . . , j -1} and any , we have that

(x + j -n -) n x + k - ∈ Q[x]. 2.3 Proof that S n (x) = O 1 x 2n+3
We shall prove that S n (x) = O 1 x 2n+3 as x → ∞ in any open angular sector that does not contain the negative real axis. The methods of [START_REF]Simultaneous polynomial approximations of the Lerch function[END_REF]Sec. 3.1] can not be used here because they lead to divergent series.

We first observe that it is enough to consider the case x → +∞ on the real axis. Indeed, S n (x) = 3Q 0,n (x)ζ(4, x) + Q 2,n (x) so that we know a priori that S n (x) has an asymptotic expansion in any angular sector that does not contain the negative real axis. Thus the leading power of this expansion can be determined by letting x → +∞ on the real positive axis.

Let N ≥ 0 be an integer. We assume that x ≥ 1. Consider the positively oriented square

C N with sides [-1 2 -iN, N + 1 2 -iN ], [N + 1 2 -iN, N + 1 2 + iN ], [N + 1 2 + iN, -1 2 + iN ], [-1 2 + iN, -1 2 -iN ].
Then by the residues theorem,

1 2iπ C N π 2 sin(πt) 2 ρ(t)dt = N k=0 ∂ ∂k k + x + n 2 (k + 1) 2 n (k + 2x) 2 n (k + x) 4 n+1 .
(The only poles of the integrand inside C N are those of π 2 sin(πt) 2 because x ≥ 1.) For any fixed real number u, π 2 sin(π(u+iv)) 2 = O(e -2π|v| ) as the real number v → ±∞, and ρ(t) = O(1/t 3 ) as t → ∞. Letting N → +∞, it follows that 1 2iπ

-1 2 +i∞ -1 2 -i∞ π 2 sin(πt) 2 ρ(t)dt = - ∞ k=0 ∂ ∂k k + x + n 2 (k + 1) 2 n (k + 2x) 2 n (k + x) 4 n+1 = S n (x).
Then for any t ∈ [-

x 2n+3 S n (x) = 4 n 2iπ -1 2 +i∞ -1 2 -i∞ π 2 sin(πt) 2 (t + 1) 2 n dt.
Similarly, we can prove that 1 2iπ

-1 2 -i∞ -1 2 +i∞ π cot(πt)ρ(t)dt = ∞ k=0 k + x + n 2 (k + 1) 2 n (k + 2x) 2 n (k + x) 4 n+1 = R n (x).
However, we can not deduce from this representation that lim x→+∞ x 2n+3 R n (x) is finite by the method above, because

-1 2 +i∞ -1 2 -i∞ π cot(πt)(t + 1) 2
n dt is divergent. In fact, it turns out that lim x→+∞ x 2n+3 R n (x) is not finite when n ≥ 1, because it can be proved that lim x→+∞ x 2 R n (x) is finite and non-zero.

3 Connections with other works

Cohen's continued fraction for ζ(4, x)

In [START_REF] Cohen | Continued fractions for Gamma products[END_REF], Cohen presented certain continued fractions for values of the Riemann zeta function and the Gamma function. In particular he stated the following one (in his notations):

ζ(4, x + 1) ≈ (2x + 1)/3 1P x (1) + 1 8 2x(2x + 2) 3P x (2) + 2 8 (2x -1)(2x + 3) 5P x (3) + • • • (9)
where

P x ( ) = 2x 4 + 4x 3 + (2 2 -2 + 4)x 2 + (2 2 -2 + 2)x -( -1)( 2 -+ 1).
He wrote that ≈ means "asymptotic expansion as the integer x → ∞", and that it is not an equality.

Maple implementation of Zeilberger's algorithm shows that our sequences (S n (x + 1)) n≥0 , (Q 0,n (x + 1)) n≥0 and (Q 2,n (x + 1)) n≥0 are solutions of the linear recurrence

n 5 U n + (2n -1)P x (n)U n-1 + (n + 1) 3 (n + 2x)(n -2 -2x)U n-2 = 0.
It is then not difficult to prove that Q 2,n (x+1) 3Q 0,n (x+1) are the convergents of Cohen's continued fraction [START_REF] Prévost | Rational approximations to the Hurwitz zeta function[END_REF]. See also Lange's paper [START_REF] Lange | Continued fraction representation for functions related to the Gamma function[END_REF] for many continued fractions related to Hurwitz zeta function, though [START_REF] Prévost | Rational approximations to the Hurwitz zeta function[END_REF] does not seem to be listed.

Cohen then mentioned that Apéry's "continued fraction acceleration" method shows

ζ(4) = 13 C(1) + 2 • 3 • 4 • 1 7 C(2) + 5 • 6 • 7 • 2 7 C(3) + • • • (10)
where

C(n) = 3(2n -1)(45n 4 -90n 3 + 72n 2 -27n + 4). (11)
He also wrote that the convergents an bn of the continued fraction [START_REF]Simultaneous polynomial approximations of the Lerch function[END_REF] are such that

ζ(4) - a n b n ≈ c(-1) n (2 + √ 3) 6n
for some constant c = 0, which is not enough to prove the irationality of ζ(4). The continued fraction [START_REF]Simultaneous polynomial approximations of the Lerch function[END_REF] had been announced before in [START_REF] Cohen | Accélération de la convergence de certaines récurrences linéaires[END_REF], with details given in [START_REF] Cohen | Accélération de la convergence de certaines récurrences linéaires[END_REF].

To conclude this section, we remark that the polynomials Q j,n,r (x) and series stated after Theorem 1 do not seem to satisfy a linear recurrence of order ≤ 2 when r < n. The case r = n is thus very remarkable.

Zudilin's approximations to ζ(4)

In [12, Section 2], Zudilin showed that for any integer n ≥ 0

Z n := - ∞ k=0 ∂ ∂k k + n 2 (k -n) 2 n (k + n + 1) 2 n (k) 4 n+1 = u n 3ζ(4) + v n
where u n and v n are rational numbers. In particular,

u n = n j=0 ∂ ∂j n 2 -j n j 4 n + j n 2 2n -j n 2
and the expression for v n is more complicated. He also proved that (Z n ) n , (u n ) n and (v n ) n are solutions of the linear recurrence

n 5 U n + C(n)U n-1 -3(3n -2)(3n -4)(n -1) 3 U n-2 = 0
where Q(n) is Cohen's polynomial [START_REF] Wilson | Some hypergeometric orthogonal polynomials[END_REF]. It can be verified that un vn coincide with the convergents an bn of ( 10 and thus we recover Zudilin's sequence (v n ) n by the identity

v n = Q 2,n (n + 1) -3u n n j=1 1 j 4 .

Prévost's remainder Padé approximants for ζ(s, x)

In [START_REF] Prévost | A new proof of the irrationality of ζ(2) and ζ(3) using Padé approximants[END_REF], Prévost showed a very original method to prove the irrationality of ζ(2) and ζ(3). We present the slightly modified approach he recently presented in [START_REF] Prévost | Rational approximations to the Hurwitz zeta function[END_REF]. For any integer x ≥ 1, we have For s = 2, the denominators of [n + 1/n](x) are

ζ(s) = x-1 k=1 1 k s + ζ(s, x).
P n (x) = n j=0 n + 1 j + 1 n + j + 2 j + 1 x -1 j , n ≥ 0,
and they satisfy the orthogonality relation

iR P n (x)P m (x) x 2 sin(πx) 2 dx = 0, n = m. For s = 3, the denominators Q n (x) of [n + 2/n](x) are such that Q n (x 2 ) = n j=0 1 j + 1 n + 1 j + 1 n + j + 2 j + 1
x -1 j

x + 1 j , n ≥ 0, and they satisfy the orthogonality relation

iR Q n (x)Q m (x) x 5 cos(πx) sin(πx) 3 dx = 0, n = m.
The two families of orthogonal polynomials (P n ) n and (Q n ) n are specializations of Wilson's orthogonal polynomials [START_REF] Wilson | Some hypergeometric orthogonal polynomials[END_REF].

Recently, Prévost [START_REF] Prévost | Rational approximations to the Hurwitz zeta function[END_REF] proved that the Padé approximants [n + 1/n](x) of ζ(s, x) at x = ∞ converge to ζ(s, x) for any fixed real number s > 1, but convergence is still an open problem when s is a complex number. Moreover, except for s = 2, 3, no expression of these approximants is known, even for s = 4. In this case, the problem is to find explicit expressions for polynomials A n (x) (of degree n) such that iR A n (x)A m (x)

x 8 (2 + cos(2πx)) sin(πx) 4 dx = 0, n = m.

Unfortunately, the weight function x 8 (2+cos(2πx)) sin(πx) 4 is not of the form studied by Wilson. The sequence (Q 0,n (x)) n is not orthogonal for this weight, but it is bi-orthogonal in the following sense: for any n and m such that 0 ≤ m ≤ 2n -1, we have

Bibliography

Tanguy Rivoal, Institut Fourier, CNRS et Université Grenoble Alpes, CS 40700, 38058 Grenoble cedex 9, France tanguy.rivoal@univ-grenoble-alpes.fr

2. 1

 1 Linear forms in 1, ζ(3, x), respectively 1, ζ(4, x)

1 j 4 ,

 14 ); see [12, Section 2, Theorem 2]. We observe that Z n = S n (n+1) and u n = Q 0,n (n+1). Since ζ(4, n+1) = ζ(4)-n j=1 the specialization of Theorem 1 at x = n + 1: S n (n + 1) = 3Q 0,n (n + 1)ζ(4, n + 1) + Q 2,n (n + 1) becomes Z n = u n 3ζ(4) + Q 2,n (n + 1) -

  He then computed explicitly the Padé approximants [n + 1/n](x) at x = ∞ of ζ(2, x), respectively the Padé approximants [n+2/n](x) at x = ∞ of ζ(3, x). After taking x = n+1, he obtained Apéry's famous sequences for ζ(2) and ζ(3).

  iR x m+5 Q 0,n (x) cos(πx) sin(πx) 3 dx = 0 = iR x m+8 Q 0,n (x) 2 + cos(2πx) sin(πx) 4 dx.

3. 4

 4 Beukers and Bel's p-adic irrationality proofsIn[START_REF] Calegari | Irrationality proofs using modular forms[END_REF], Calegari proved the irrationality of the 2-adic numbers ζ 2 (2) and ζ 2 (3), as well as of the 3-adic numbers ζ 2 (3). His proof used overconvergent p-adic modular forms. Later, Beukers[START_REF] Beukers | Irrationality of some p-adic L-values[END_REF] obtained another proof of these facts, of a more classical flavor. In fact, he essentially used Prévost's Padé approximants for ζ(2, x) and ζ(3, x), though his formulas are written differently. The Padé type approximants constructed in[START_REF]Simultaneous polynomial approximations of the Lerch function[END_REF] for ζ(s, x) contain as initial cases Beukers and Prévost approximants; Bel[START_REF] Bel | Fonction Zêta de Hurwitz p-adique et irrationalité[END_REF] used them to prove certain linear independence results for values of p-adic Hurwitz zeta functions. It would be interesting to know if Theorem 1 or its generalization could be used to prove the irrationality of the numbers ζ p (4) for some p. The arithmetic and asymptotic properties of Zudilin's series Z n are not good enough to imply the irrationality of ζ(4), but a modification of Z n conjecturally proves that ζ(4) / ∈ Q.

  for some constant c n > 0 independent of x and t. Since -1/2+iR | π 2 sin(πt) 2 (t + 1) 2 n ||dt| < ∞ for any n, and lim x→+∞ x 2n+3 ρ(t) = 4 n (t + 1) 2

	n , it follows by the dominated convergence
	theorem that
	lim x→+∞

1 2 -i∞, -1 2 + i∞] and any x ≥ 1, we have |x 2n+3 ρ(t)| ≤ c n |(t + 1) 2 n |


Karim Kelfoun

A two-layer depth-averaged model for both the dilute and the concentrated parts of pyroclastic currents

Pyroclastic currents are very destructive and their complex behavior makes the related hazards difficult to predict. A new numerical model has been developed to simulate the emplacement of both the concentrated and the dilute parts of pyroclastic currents using two coupled depthaveraged approaches. Interaction laws allow the concentrated current (pyroclastic flow) to generate a dilute current (pyroclastic surge) and, inversely, the dilute current to form a concentrated current or a deposit. The density of the concentrated current is assumed to be constant during emplacement, whereas the density of the dilute current changes depending on the particle supply from the concentrated current and the mass lost through sedimentation. The model is explored theoretically using simplified geometries as a proxies for natural source conditions and topographies. It reproduces the relationships observed in the field between the surge genesis and the topography: the increase in surge production in constricted valleys, the decoupling between the concentrated and the dilute currents and the formation of surge-derived concentrated flows. The strong non linear link between the surge genesis and the velocity of the concentrated flow beneath it could explain the sudden occurrence of powerful and destructive surges and the difficulty of predicting this occurrence. A companion paper compares the results of the model with the field data for the eruption of Merapi in 2010 [Kelfoun et al.,

submitted] and demonstrates that the approach is able to reproduce the natural emplacement of the concentrated and the dilute pyroclastic currents studied with good accuracy.

Introduction

Pyroclastic currents are fast-moving density currents composed of hot gas and rock fragments.

For small to medium-scale eruptions, generated by column collapse or dome collapse, pyroclastic density currents (PDC) are generally formed of two distinct parts: a concentrated and a dilute part [START_REF] Lacroix | Settling velocity of spheres : a new appaoch[END_REF][START_REF] Sparks | Products of ignimbrite eruptions[END_REF][START_REF] Walker | Lateral variations in the Taupo ignimbrite[END_REF][START_REF] Valentine | Origin of layer 1 deposits in ignimbrites[END_REF]. The concentrated part, termed block-and-ash flow, scoria flow or pumice flow depending on its composition, or more generally called pyroclastic flow [START_REF] Branney | Pyroclastic Density Currents and the Sedimentation of Ignimbrites[END_REF], follows existing valleys and is sensitive to meter-scale variations in the topography. Following the nomenclature used in the volcanological community, the term 'flow' will be used here for the concentrated part of PDCs only [START_REF] Branney | Pyroclastic Density Currents and the Sedimentation of Ignimbrites[END_REF]. Particle concentration in the flow is close to particle concentration at rest, which means that particle interactions probably play a major role in its dynamics. The dilute part of the current is called dilute current, pyroclastic surge or ash-cloud surge. It is formed of fine particles held in suspension by turbulent gas [e.g. [START_REF] Branney | Pyroclastic Density Currents and the Sedimentation of Ignimbrites[END_REF][START_REF] Valentine | Stratified flow in pyroclastic surges[END_REF]. Ash-cloud surges are generally several tens of meters in thickness but they often leave deposits of only a few centimetres-decimetres due to their low particle concentrations (about 1% in volume [START_REF] Druitt | Pyroclastic density currents[END_REF]). Note that this concentration is higher than the maximal concentration of 0.1% that defines dilute mixtures in the multiphase flow literature, i.e. below which collisions can be neglected [e.g. Elgobashi, 1994]. Pyroclastic surges are less confined by the topography, escape easily from the valleys and are very dangerous for the inhabitants. They are believed to be stratified, their density increasing downward [START_REF] Valentine | Stratified flow in pyroclastic surges[END_REF].

The two types of currents interact during their emplacement and can exchange particles and gas. The particles, falling down to the ground through sedimentation of the dilute current, either form a deposit, increase the volume of the concentrated flow or create a surge-derived pyroclastic flow [START_REF] Druitt | A proximal ignimbrite breccia facies on Santorini, Greece[END_REF][START_REF] Calder | Mobility of pyroclastic flows and surges at the Soufrie`re Hills Volcano, Montserrat[END_REF][START_REF] Druitt | Small-volume, highly mobile pyroclastic flows formed by rapid sedimentation from pyroclastic surges at Soufrière Hills Volcano, Montserrat: an important volcanic hazard[END_REF]. Inversely, a dilute current can be generated by the concentrated flow and then can acquire its own dynamics [START_REF] Fisher | Models for pyroclastic surges and pyroclastic flows[END_REF][START_REF] Nakada | Preliminary report on the activity at Unzen Volcano (Japan), November 1990-November 1991: dacite lava domes and pyroclastic flows[END_REF][START_REF] Yamamoto | Pyroclastic flows from the 1991 eruption of Unzen volcano, Japan[END_REF][START_REF] Fisher | Decoupling of pyroclastic currents: hazards assessments[END_REF][START_REF] Kelfoun | A statistical study of trees damaged by the 22 November 1994 eruption of Merapi volcano (Java, Indonesia): relationships between ash-cloud surges and block-and-ash flows[END_REF][START_REF] Bourdier | Decoupling of small volume pyroclastic flows and related hazards at Merapi volcano, Indonesia[END_REF][START_REF] Breard | Coupling of turbulent and non-turbulent flow regimes within pyroclastic density currents[END_REF].

Several approaches have been used to simulate PDCs on digital elevation models (DEMs). For example, the kinematic approach only simulate flow trajectories, considering the flow as a punctual mass [Sheridan andMalin, 1983, Beget and[START_REF] Beget | Two-dimensional kinematic and rheological modeling of the 1912 pyroclastic flow, Katmai, Alaska[END_REF][START_REF] Mcewen | Dynamics of Mount St. Helens' 1980 pyroclastic flows, rockslide-avalanche, lahars, and blast[END_REF][START_REF] Rossano | Numerical simulation of pyroclastic density currents on Campi Flegrei topography: a tool for statistical hazard estimation[END_REF][START_REF] Saucedo | Modeling of pyroclastic flows of Colima Volcano, Mexico: implications for hazard assessment[END_REF]. In the multiphase approach, all the phases (different classes of particles, gases of various compositions) are taken into account and interact together in 3D [e.g. [START_REF] Valentine | Numerical models of plinian eruption columns and pyroclastic flows[END_REF][START_REF] Valentine | Effects of topography on facies and compositional zonations in caldera-related ignimbrites[END_REF][START_REF] Neri | Numerical simulation of collapsing volcanic columns with particles of two sizes[END_REF][START_REF] Neri | Multiparticle simulation of collapsing volcanic columns and pyroclastic flows[END_REF][START_REF] Dartevelle | Numerical modeling of geophysical granular flows: 2. Computer simulations of plinian clouds and pyroclastic flows and surges[END_REF][START_REF] Esposti Ongaro | Transient 3D numerical simulations of column collapse and pyroclastic density current scenarios at Vesuvius[END_REF][START_REF] Esposti Ongaro | Multiphase flow dynamics of pyroclastic density currents during the May 18, 1980 lateral blast of Mount St. Helens[END_REF]. In the depth-averaged method, all the physical properties and the equations are integrated perpendicular to the ground (or sometimes vertically). Compared to the 3D solving of the multiphase approach, this integration allows to speed up the calculation time and to take into account small details of the topography (models can reproduce the natural runout of geophysical flows on high resolution accurate topographies (< 5 m) in a few tens of minutes).

However, a limit of the approach is its accuracy on complex topographies with high curvatures.

For a better understanding of the limits and new developments of the depth-averaged approach on complex topographies, the reader can refer to [START_REF] Hutter | The Savage-Hutter avalanche model: how far can it be pushed?[END_REF], [START_REF] Luca | Two-layer models for shallow avalanche flows over arbitrary variable topography[END_REF], [START_REF] Castro-Orgaz | Nonhydrostatic granular flow over 3-D terrain: New Boussinesq-type gravity waves?[END_REF], [START_REF] Luca | Shallow Geophysical Mass Flows down Arbitrary Topography[END_REF]. The integration of the depth-averaged approach also implies that the flow is wider than it is thick and that vertical acceleration is negligible. The depth-averaged method is often used to simulate granular flows in the laboratory [START_REF] Savage | The motion of a finite mass of granular material down a rough incline[END_REF][START_REF] Pouliquen | Friction law for concentrated granular flows: Application to the motion of a mass down a rough inclined plane[END_REF][START_REF] Gray | Shock waves, dead zones and particle-free regions in rapid granular free-surface flows[END_REF] and geophysical flows like mud flows, pyroclastic flows, landslides, debris avalanches and snow avalanches [e.g. [START_REF] Itoh | Hazard estimation of the possible pyroclastic flow disasters using numerical simulation related to the 1994 activity at Merapi Volcano[END_REF][START_REF] Heinrich | Numerical simulation of the December 1997 debris avalanche in Montserrat, Lesser Antilles[END_REF][START_REF] Turnbull | Mass and Momentum Balance Model of a Mixed Flowing/Powder Snow Avalanche[END_REF][START_REF] Pitman | Computing granular avalanches and landslides[END_REF][START_REF] Kelfoun | Numerical modeling of the emplacement of Socompa rock avalanche, Chile[END_REF][START_REF] Sheridan | Evaluating Titan 2D mass flow model using the 1963 Little Tahoma Peak avalanches[END_REF][START_REF] Patra | Parallel adaptive numerical simulation of dry avalanches over natural terrain[END_REF][START_REF] Lucas | Mobility and topographic effects for large Valles Marineris landslides on Mars[END_REF][START_REF] Kelfoun | Is frictional behaviour suitable for pyroclastic flow simulation: Comparison with a well constrained eruption at Tungurahua volcano (Ecuador)[END_REF]. The depth-averaged method has also been developed to take into account 2 phases, solid and water, for the simulation of debris flows [e.g. [START_REF] Iverson | Flow of variably fluidized granular masses across three-dimensional terrain: 1. Coulomb mixture theory[END_REF][START_REF] Córdoba | TITAN2F: a pseudo-3-D model of 2phase debris flows[END_REF], and to simulate granular flows fluidized by gases [START_REF] Gueugneau | Effects of pore pressure in pyroclastic flows: Numerical simulation and experimental validation[END_REF].

The integration and resolution of the depth-averaged methods is mainly done for currents whose density is constant in time and space, and it has been shown that it can reproduce accurately the velocity, deposit thickness and extension of concentrated pyroclastic flows [START_REF] Kelfoun | Is frictional behaviour suitable for pyroclastic flow simulation: Comparison with a well constrained eruption at Tungurahua volcano (Ecuador)[END_REF]. Theoretically, the depth-averaged approach is less adapted for dilute currents because of their vertical variations in density. However, the displacement of a density current is due to a pressure gradient parallel to the ground and the pressure is more influenced by the total mass above a given point than by its vertical variations. Based on this reasoning, [START_REF] Doyle | Modeling concentrated pyroclastic basal flows from collapsing columns[END_REF] have explored the formation of concentrated flow from a dilute current using a depth-averaged model. Their approach is clearly simplified compared to the multiphase approach but the limited number of variables of the model allows the influence of each one on the global results to be studied.

The aim of the present work is to present a new depth-averaged model that simulates both the concentrated and the dilute part of pyroclastic currents and their interactions. The dilute current can form a concentrated flow, as in [START_REF] Doyle | Modeling concentrated pyroclastic basal flows from collapsing columns[END_REF]. The innovation in our model is that, inversely, the concentrated flow is able to form a dilute current, as often occurs during dome collapses. The model also runs on a pseudo 3D surface (Digital Elevation Model of volcanic topography) and the concentrated part uses a rheology estimated from field observations, which allows the dilute current to form either a deposit or a surge-derived pyroclastic flow depending on the thickness deposited.

The two-layer model

The objective of this work is to develop the simplest model for concentrated and dilute currents that reproduces the first-order features of the natural phenomena. The simulation of the concentrated flow must approximate the velocity of the real flow, the covered area (extent, runout) and the deposit thickness. The model must also address the way the surge is generated by the concentrated flow (Fig. 1), as well as reproducing the velocity, the thickness of the surge and the area it covers. Finally, the surge should lose mass to form a deposit or a surge-derived pyroclastic flow compatible with natural observations. Other phenomena such as the effects of break in slope, the influence of river width changes and the decoupling between the concentrated and the dilute currents [e.g. [START_REF] Ogburn | The effect of topography on ash-cloud surge generation and propagation[END_REF] must also be tested to assess the quality of the model.

Dynamics of the concentrated flow

The dynamics of the concentrated flow is calculated by solving depth-averaged equations of mass [START_REF] Luca | Two-layer models for shallow avalanche flows over arbitrary variable topography[END_REF] and momentum balance, along x in (2) and y in (3), taking into account the mass of the particles exchanged with the dilute current.

    d d s m d d x y h h u h u t x y              [1]         2 d d d 2 d d s d d 1 sin cos 2 x x xy x x x x x d m h u h u h u u t x y T gh gh v u x                       [2]         2 d d d 2 d d s d d 1 sin cos 2 y xy y y y y y y d m h u h u u h u t x y T gh gh v u y                       [3]
where hd is the thickness of the concentrated flow, ( , ) x y u u  u its velocity, x and y the slope of the topography in the x and y directions, t the time and ( , )

x y v v  v
the velocity of the surge (Fig. 1). The concentrated flow is considered to be of constant density, d, in time and space allowing the mass balance to be simplified by a thickness balance in [START_REF] Luca | Two-layer models for shallow avalanche flows over arbitrary variable topography[END_REF].

The variables m and s are the mass fluxes per square meter of the mixture of gas and particles that move from the concentrated to the dilute current (see section 2.3), and of particles that leave the dilute current through sedimentation (see section 2.4), respectively (Fig. 1). The terms at the right-hand side of equations 1 -3, which contain m and s, govern the mass and momentum exchanges between the concentrated and the dilute currents.

The flow is slowed down by the resisting stress exerted by the ground, ( , ) x y

T T  T

, which depends on the rheology of the flow. The rheology of pyroclastic flows is very complex and still not clearly understood. [START_REF] Kelfoun | Is frictional behaviour suitable for pyroclastic flow simulation: Comparison with a well constrained eruption at Tungurahua volcano (Ecuador)[END_REF] and [START_REF] Kelfoun | Suitability of simple rheological laws for the numerical simulation of concentrated pyroclastic flows and long-runout volcanic avalanches[END_REF] have shown that a plastic rheology forms flows and deposits whose characteristics are close to the natural phenomena in terms of extension, thickness, thickness variation, and velocity, although the reason why this is so is not fully understood [START_REF] Kelfoun | Suitability of simple rheological laws for the numerical simulation of concentrated pyroclastic flows and long-runout volcanic avalanches[END_REF]. [START_REF] Kelfoun | Is frictional behaviour suitable for pyroclastic flow simulation: Comparison with a well constrained eruption at Tungurahua volcano (Ecuador)[END_REF] also noted that a velocitydependent term could be added to the plastic rheology to better reproduce the natural emplacement of the Tungurahua pyroclastic flows, formed by the destabilisation of pyroclasts accumulated by fountaining. However, for such eruptions with a relatively small mass flux feeding the flow, the influence of this term on the dynamics was too small to be clearly identified from field observations. [START_REF] Charbonnier | Evaluation of geophysical mass flow models using the 2006 block-and-ash flows of Merapi Volcano, Java, Indonesia: towards a shortterm hazard assessment tool[END_REF] and [START_REF] Charbonnier | Evaluation of the impact of the 2010 pyroclastic density currents at Merapi volcano from high-resolution satellite imagery, field investigations and numerical simulations[END_REF] have reproduced the emplacement of the concentrated pyroclastic flows of the Merapi volcano for the 2006 and 2010 eruptions. They have confirmed that a plastic rheology captures the flow emplacement better than a Coulomb rheology and they have shown that, for the high mass flux related to a dome collapse, the additional term related to the square of the flow velocity must be used to reproduce the natural velocity of the flow. The retarding stress T they have used is expressed by:

0 1 d T a     u T uu u [4]
and its components in x and y are

0 1 x x d x u T T a u     u u and 0 1 y y d y u T T a u     u u .
The best value estimated by [START_REF] Charbonnier | Evaluation of geophysical mass flow models using the 2006 block-and-ash flows of Merapi Volcano, Java, Indonesia: towards a shortterm hazard assessment tool[END_REF] and [START_REF] Charbonnier | Evaluation of the impact of the 2010 pyroclastic density currents at Merapi volcano from high-resolution satellite imagery, field investigations and numerical simulations[END_REF] for the two Merapi eruptions is a1 = 0.01. The velocity related term, a1, influences the flow essentially during the initial phase, when both the mass flux and the flow velocity are high. It might correspond to an additional absorption of energy related to collisions and rock break-up during the lava dome destruction. Neither the temperature and its variations nor the erosion and the incorporation of substratum material in the flow [START_REF] Bernard | Pyroclastic flow erosion and bulking processes: comparing field-based vs. modeling results at Tungurahua volcano, Ecuador[END_REF] are included in the present model. Aggradation is also not included in the simulation of the concentrated part, and the flow stops en masse according to the balance between driving stresses, retarding stresses and inertia.

Dynamics of the dilute current

In the same way as for the concentrated flow, the dilute current (i.e. pyroclastic surge) is simulated by equations of mass and momentum balance. However, a fourth equation must be added to the system because the density varies in time and space. The surge is considered to be incompressible and isothermal following [START_REF] Doyle | Modeling concentrated pyroclastic basal flows from collapsing columns[END_REF][START_REF] Doyle | A two-layer model for the evolution and propagation of concentrated and dilute regions of pyroclastic currents[END_REF]. The turbulent mixing is assumed to be sufficiently intense to maintain a uniform particle concentration perpendicular to the ground. The system describing the surge dynamics is then:

    s s s m m x y h h v h v t x y            [5]     s s s m s s s x s y h h v h v t x y                [6]             1 cos 2 x x x y x x x x x h v h v h v v gh t x y gh u v x                             R [7]                sin 1 cos 2 y x y y y y y y y h v h v v h v gh t x y gh u v y                             R     [8]
Eq. 5 and 6 describe the volume and mass balances (per square meter) for a current whose density varies and whose thickness is not affected by particle settling [START_REF] Doyle | A two-layer model for the evolution and propagation of concentrated and dilute regions of pyroclastic currents[END_REF].

Equations 7 and 8 are the momentum balances. The right-hand terms of Eq. 7 and 8 mean [START_REF] Luca | Two-layer models for shallow avalanche flows over arbitrary variable topography[END_REF] that the surge can move on a slope due to its weight because it is denser than the atmosphere.

It can also move by spreading (2), even on a horizontal surface. The third term (3) includes the stresses that resist to the displacement and the last terms (4) and (5) take into account the momentum gain and lost by exchanges with the concentrated part. Note that no overpressure is taken into account and that blast effects cannot be simulated.

Here, hs is the thickness of the dilute current, ( , )

x y v v  v its velocity and, x and y the slope of the surface formed by the ground plus the concentrated flow in the x and y directions. The density of the surge varies and can be close to the density of the atmosphere, a, which cannot be neglected (Eq. 7-8). The density of the mixture of gas and particles that moves from the concentrated flow to the dilute current is termed m (see section 2.3). As for the concentrated flow, the two terms at the right-hand side of each equation take into account the mass exchange, here between the dilute current and, the concentrated flow or deposit (see section 2.4).

The surge is considered to exhibit turbulent drag stress [e.g. [START_REF] Bonnecaze | Particle-driven gravity currents[END_REF][START_REF] Bonnecaze | Particle-driven gravity currents down planar slopes[END_REF][START_REF] Hogg | The effects of hydraulic resistance on dam-break and other shallow inertial flows[END_REF][START_REF] Doyle | Modeling concentrated pyroclastic basal flows from collapsing columns[END_REF], R, expressed by:

2 s a    R v v [9]
where a2 is a drag coefficient that relates the velocity of the flow to the drag stress exerted by the turbulence. The set of equations ( 5) to ( 8) is valid only where the density of the ash-cloud is higher than the density of the atmosphere, with the ash-cloud behaving as a density current (i.e. an ash-cloud surge).

Air entrainment can have a strong influence on the surge emplacement [START_REF] Bursik | The effects of topography on sedimentation from particleladen turbulent density currents[END_REF][START_REF] Andrews | Dispersal and air entrainment in unconfined dilute pyroclastic density currents[END_REF]: atmospheric gases are ingested by the flow, become heated and thus decrease the surge density. Air ingestion will increase the flow thickness, simultaneously decreasing the flow density [e.g. [START_REF] Hallworth | Entrainment into two-dimensional and axisymetric gravity currents[END_REF]. Thus, the term s s h

 that rules the dynamics of the flow (Eq. 6-8) will not drastically change and, while the density of the surge is significantly higher than the density of the atmosphere, entrainment influences the flow thickness but probably only has a secondary effect on the emplacement (see section 3.7). Air entrainment is not taken into account in the model as in the approach of [START_REF] Doyle | Modeling concentrated pyroclastic basal flows from collapsing columns[END_REF][START_REF] Doyle | A two-layer model for the evolution and propagation of concentrated and dilute regions of pyroclastic currents[END_REF]. This assumption affects the model accuracy for very dilute surges and at the edges of dense surges where air ingestion can decrease the surge density down to the atmosphere density, causing the ash-cloud to become buoyant.

Mass and momentum exchanges: from the concentrated to the dilute current

During the emplacement of the concentrated flow, particles at its surface can be put into suspension to form a mixture of gas and particles that supplies the ash-cloud, changing its thickness and its density (Fig. 1). A number of processes have been proposed for particle transfer from the concentrated flow to the dilute current: (1) turbulence [START_REF] Denlinger | A model for generation of ash clouds by pyroclastic flows, with application to the 1980 eruptions at Mount St. Helens, Washington[END_REF] (2)

elutriation [START_REF] Fisher | Models for pyroclastic surges and pyroclastic flows[END_REF][START_REF] Wilson | The role of fluidization in the emplacement of pyroclastic flows: an experimental approach[END_REF];

(3) decompression of pressurized volcanic gases [START_REF] Woods | The explosive decompression of a pressurized volcanic dome: the 26 December 1997 collapse and explosion of Soufrie`re Hills Volcano, Montserrat[END_REF]; (4) entrainment by heated atmospheric gases [Mc Taggart, 1960;[START_REF] Sparks | Grain size variations in ignimbrites and implications for the transport of pyroclastic flows[END_REF][START_REF] Wilson | Ignimbrite depositional facies: the anatomy of a pyroclastic flow[END_REF]; (5) particle collisions [Spark et al., 1978;[START_REF] Fujii | The 15 September 1991 pyroclastic flows at Unzen Volcano (Japan): a flow model for associated ash-cloud surges[END_REF]. The gas of the mixture can originate either from the concentrated flow or from entrainment and heating of the atmosphere. The density of the mixture, which takes these phenomena implicitly into account, is termed m and the mass flux (per m² of the surface and per s) of gas and particles forming this mixture and supplying the dilute current is termed m. It is impossible to carry out direct measurements of particle exchanges (mass flux rate, the grain size distribution, etc.) between the concentrated and the dilute parts of natural pyroclastic currents. Wind tunnel experiments, in which an air flux is imposed above a granular layer, are useful to understand, at least partially, the mechanisms involved in particle entrainment into a saltation layer. The experiments show that the air flux interacts with the roughness caused by the grains of the granular layer surface and becomes turbulent [Denlinger, 1989]. This turbulence plays an important role in particle mobilization by wind [e.g. [START_REF] Bagnold | The Physics of Blown Sand and Desert Dunes[END_REF]. While various laws can be found in the literature, the mass flux of particles, per square meter, is often considered to be a function of the wind velocity cubed, with no flux below a minimal velocity value [START_REF] Bagnold | The Physics of Blown Sand and Desert Dunes[END_REF]Pye and Tsoar,1990;[START_REF] Tegen | Contribution to the atmospheric mineral aerosol load from land source modification[END_REF][START_REF] Tegen | Contribution of different aerosol species to the global aerosol extinction optical thickness; estimates from model results[END_REF][START_REF] Uno | Trans-Pacific yellow sand transport observed in April 1998: a numerical simulation[END_REF][START_REF] Xuan | Turbulence factors for threshold velocity and emission rate of atmospheric mineral dust[END_REF]. This threshold is generally low (some tens of cm / s) compared to concentrated pyroclastic flow velocities (> 20 m/s) and here it is assumed to be zero for a sake of simplification.

The mass flux of the gas and particle mixture, m (in kg / m² / s), which leaves each square meter of the surface of the concentrated part to supply the dilute part, is thus modelled by: 3 m 3 a   u [START_REF] Patra | Parallel adaptive numerical simulation of dry avalanches over natural terrain[END_REF] where a3 is a parameter that relates the mass flux of particles to the concentrated flow velocity (in kg s 2 m -5 ). Eq. 10 certainly highly simplifies reality, the natural mass flux being potentially complicated by particle collisions, grinding, shocks and hydraulic jumps related to the concentrated flow displacement, by the high temperature and the gas content of the concentrated flow and its degasing, and by the density and velocity of the overlying surge [START_REF] Wilson | Ignimbrite depositional facies: the anatomy of a pyroclastic flow[END_REF]]. Since we do not know the influence of these parameters, we use the simplest equation form deduced from available experiments, while recognizing that more researches are needed to better defined the physical processes involved.

The mixture of particles and gases (atmospheric and lost from the concentrated flow) change both the surge density and its thickness (Fig. 1). The thickness of the surge varies locally with time:

s m m dh dt    [11]
The product of the surge density multiplied by its thickness, s hs, i.e. the mass of the surge per surface unit, also varies according to m:

  s s m d h dt    [12]
The new surge density is deduced from equations ( 11) and ( 12). In the model, the surge velocity is constant perpendicular to the ground and our depth-averaged approach implicitly assumes that vertical mixing is locally instantaneous across the dilute current thickness.

Simultaneously, considering that the mass of the entrained atmosphere is small compared to the mass of gas and particles lost from the concentrated flow, the flow thickness decreases according to m and to its density, d:

d m d dh dt     [13]
The flux of momentum is calculated from the mass flux m  and the velocity of the concentrated part u . The concentrated flow will lose in momentum m  u by per square meter and per second, and the dilute current will gain the same momentum. The velocity of the concentrated flow is not affected by the particle flux because mass and momentum decrease together.

However, this mass transfer can modify the surge velocity.

Mass and momentum exchanges: from the dilute to the concentrated flow

The settling velocity, wsed, of particles in a surge depend on their density, shape and size, and is often given by: [START_REF] Lacroix | Settling velocity of spheres : a new appaoch[END_REF][START_REF] Sparks | Volcanic Plumes[END_REF][START_REF] Dellino | The analysis of the influence of pumice shape on its terminal velocity[END_REF]

  sed 4 3 p s p d s g d w C      [14] [Le
references therein].
where p is the density of particles, dp their diameter and Cd a drag coefficient.

The mass of particles lost by sedimentation (per square meter) depends on this settling velocity wsed, the volumetric fraction of particles cv and the density of the particles p:

s s e d v p c w    [15]
The surge density is given by

  s 1 v p v g c c      
, where g is the density of the hot gas (from volcanic and/or atmospheric origin) of the surge. Note that s decreases progressively to g as the surge loses particles by sedimentation. Because the volumetric fraction of particles is given by:

      s s v g p g g p c             [16]
Eq. 15 becomes:

  s s e d s g w      [17]
The local variation of the mass of the surge per square meter, s hs, through time due to sedimentation is then given by:

    sed s s s s g d h w dt          [18]
To solve equation 18, we follow the assumption of [START_REF] Doyle | A two-layer model for the evolution and propagation of concentrated and dilute regions of pyroclastic currents[END_REF] (their equations 2a-2c) which considers that the sedimentation does not affect the surge thickness but only its density.

Where the density of the surge becomes lighter than the atmosphere, the ash-cloud should lift off. The model then considers that the particles of the ash-cloud leave the system at the point where lift-off should occur and the surge thickness is locally set to zero.

The mass lost by the surge can form a deposit, be added to the existing concentrated flow or form surge-derived concentrated flows [START_REF] Calder | Mobility of pyroclastic flows and surges at the Soufrie`re Hills Volcano, Montserrat[END_REF][START_REF] Druitt | Small-volume, highly mobile pyroclastic flows formed by rapid sedimentation from pyroclastic surges at Soufrière Hills Volcano, Montserrat: an important volcanic hazard[END_REF]. This depends on the thickness of mass deposited, the topographic slope and the rheological parameters chosen for the concentrated part.

The mass of particles lost by the surge equals the mass gained by the concentrated part (concentrated flow or deposit) and is given by:

      sed d d s s s s g d h d h w dt dt           [19]
Because the density of the concentrated part is considered to be constant, equation 12 can be simplified to:

  sed s g d d dh w dt      [20] 2.1.
Numerical resolution of the equations

The system of equations is solved using the VolcFlow code [START_REF] Kelfoun | Numerical modeling of the emplacement of Socompa rock avalanche, Chile[END_REF][START_REF] Kelfoun | Is frictional behaviour suitable for pyroclastic flow simulation: Comparison with a well constrained eruption at Tungurahua volcano (Ecuador)[END_REF], which was developed for the simulation of depth-averaged flows on DEMs. The code has been highly modified to take into account two fluids simultaneously (the concentrated and the dilute currents), the laws simulating the mass and momentum exchanges (Eq. 10-20),

and the density variations of the surge (Eq. 6-8).

Results

The parameters of the following simulations are listed in Tab. 1. If not specified, the parameters used are realistic values either obtained by field measurements or estimated by reproducing natural events in the companion paper [Kelfoun et al., submitted].

Characteristics and evolution of the dilute current

In the model, the surge characteristics depend on the mass, momentum and density of the mixture formed by the concentrated flow, and on the particles lost by sedimentation. The following section illustrates how equations 10 and 17 affect the characteristics of the dilute part.

For this first example, no lateral displacement is taken into account: the mass of the dilute part comes from the underlying concentrated part and subsequently settles on it. Nevertheless, to be consistent with the other sections, the term surge is used to refer to the dilute part. effects on the sedimentation: first, it reduces the settling velocity of a given particle (Eq. 14).

However, because a higher density also means a larger concentration of particles, more particles settle out, thus increasing the mass lost (Eq. 17). Fig 2a shows that the second effect dominates and that a higher surge concentration increases the mass lost. This is true up to a concentration of about 1000 kg/m 3 , which is much higher than the density of a surge. However, except for low density surges (<5 kg/m 3 ), the density has a relatively small influence compared to the velocity of the concentrated flow. The bold curve shows that, with the parameters chosen (Tab.

1), a velocity of the pyroclastic flow of about 15-25 m/s is needed to generate and maintain a surge. Note that another condition for surge generation is that the particles can be held in suspension and can be mixed homogeneously through the whole surge thickness. It is commonly assumed that a particle can be held in suspension if its Rouse number * sed n w k P v  is lower than 2.5, with v* the shear velocity (in the order of magnitude of the surge velocity v ) and k the Von Karman constant, equal to 0.4 [START_REF] Valentine | Stratified flow in pyroclastic surges[END_REF][START_REF] Dellino | A method for the calculation of the impact parameters of dilute pyroclastic density currents based on deposit particle characteristics[END_REF]. The capacity for homogeneous mixing can be estimated using the Stokes number, which is a ratio of the characteristic response times for particles to the fluid turbulence [e.g. [START_REF] Elghobashi | On Predicting Particle-Laden Turbulent Flows[END_REF][START_REF] Sijercic | Simulation of free turbulent particle-laden jet using Reynolds-stress gas turbulence model[END_REF]. The Stokes number is defined by

0 0 0 t v Stk  
, where v0 is the characteristic velocity of the vortices, ℓ0 their characteristic length and t0 is the relaxation time of the particle defined by often found in surge deposits, are easy to be put into suspension (Pn<2.5) and mixed homogeneously (Stk<<1), while particles larger than 1 cm can only be mixed under more exceptional conditions (large eddies, high surge viscosity, low particle density). This should be kept in mind so as not to over interpret the upper part of Fig. 2b.

Surge genesis from a concentrated flow

Equations of mass fluxes are now coupled with Eqs. 1-3 and 5-8 to be studied dynamically. In the first simulation, in 1D (Fig. 3a-b), a mass representing a concentrated flow is released in one go onto a slope in order to approximate a dome collapse. For visibility, a high value of T0 has been chosen so as to create a thick flow. The flow rheology contains a plastic term (Eq. 4), meaning that there is a relation between its thickness and its ability to flow. Where the flow is thick, the driving stresses (related to the slope and the thickness) are higher than the retarding stress (which is constant) and the flow accelerates (Eq. 2). The high flow velocity induces a flux of particles (Eq. 10) that forms a surge when it becomes greater than the mass flux lost by sedimentation (Eq. 17). The mass of the surge increases and the surge acquires its own dynamics (Eq. 7-8). The concentrated flow thins until the driving stress becomes lower than the resisting stress and the velocity decreases until the flow stops. Fewer particles are emitted from the flow and the density of the dilute current decreases rapidly through sedimentation. However, its density can be high enough to enable the dilute current to surpass the front of the concentrated flow (Fig. 3b). The dilute part stops when its density equals the density of the atmosphere and the mass lifts off. Thus, depending on the topography and on the characteristics of the concentrated and the dilute parts (rheological parameters, mass rate, volume, etc.), the dilute current can be slower (Fig. 3a) or faster (Fig. 3b) than the concentrated flow and can stop before or after the runout reached by the concentrated flow.

Concentrated flow genesis from a surge

In the simulation shown in Fig. 3c-d a surge is continuously created from the border of the calculation domain to mimic, for example, a column collapse. The surge moves downstream and progressively loses mass by sedimentation, resulting in a decrease in density. The runout of the surge is determined by the point at which the density of the surge equals the density of the atmosphere (Fig. 3d). Initially, this mass lost by sedimentation of the surge forms a deposit (i.e. a static concentrated part). But, because the thickness of the concentrated part increases with time, the driving stress of the weight increases until the concentrated part begins to flow (Fig. 3d). Concentrated flow genesis occurs in pulses: the mass lost by the surge accumulates on the ground until it is able to flow, then builds up again to the point at which it flows again.

The concentrated part is able to exceed the surge runout (Fig. 3d).

Influence of the slope on surge genesis

The influence of the slope on the surge generated by the concentrated flow in the model can first be analyzed in stationary conditions with a purely plastic rheology. The influence of the turbulent stress (Eq. 4) is generally of secondary importance and is set here to zero.

A concentrated flow reaches a constant velocity when the driving force equals the resisting force. For a purely plastic rheology, and a flow whose thickness does not vary spatially, a stationary state (no variation of h over time) is reached when the driving force equals the plastic basal stress:

0 sin g h     [21]
This gives the thickness needed by the flow to reach a constant velocity:

0 sin h g     [22]
A thicker flow will accelerate to thin and thus reach the equilibrium thickness. Conversely, a thinner flow will decelerate to thicken.

In this simple analysis the valley is considered to be horizontal perpendicular to the valley axis with vertical walls. For a volumetric rate Q of pyroclastic flow, the stationary velocity is given by:

0 sin Q Q g u h w w      [23]
where w is the valley width. Using Eq. 10, the mass flux from the concentrated flow to the dilute current is then:

3 3 m 3 3 0 sin Q g a a w             u
(per m² and per s)

[24]

The mass flux that forms the surge is very sensitive to the concentrated flow velocity and, consequently, to all the parameters that control this velocity, one of which is the slope. For example, increasing the slope from 10° to 30° increases the mass flux by about 25. Decreasing the slope from 10° to 5° decreases the mass flux by about 8.

Another illustration of this slope effect is given by a dynamic analysis of the model (full set of equations). A first simulation is done on a constant slope of 5°. The valley is linear downslope and is defined by a sinusoidal cross section (=50 m, =20 m, Fig. 4a). The concentrated flow is generated at the left side of the domain (x=0). It accelerates to reach ~12 m/s and generates a dilute current that overflows the valley. With the parameters chosen (Tab. 1), the surge reaches a maximal density of ~ 9 kg/m 3 above the flow, and because its density decreases rapidly it moves laterally by about 70 m from the concentrated part (Fig. 4b). Streamlines show that the surge that covers a given area is generated in the proximal area covered by the flow (downstream displacement of less than 100 m).

The simulation in Fig. 4c is done with the same slope and the same valley shape but a steeper slope of 30° is added at a distance of 200 m. The effect of the steeper slope is to accelerate the flow and this occurs even before the break in slope: because the mass on the steeper slope accelerates, the thickness gradient changes. This adds an additional driving stress to the concentrated flow (Eq. 2-3) and the mass accelerates, propagating the thickness change and the acceleration upstream of the break in slope. On the steeper slope the mass accelerates from 12 m/s to about 22 m/s and forms a thick, concentrated surge that spreads out in the downstream area and is able to flow faster than the concentrated flow. The surge formed on the steep slope is able to progress more than 500 m downstream.

The surge deposits are 12% of the initial flow volume for the 5° slope and 30% for the steeper slope. This mass is removed from the flow and, because the flow stop is related to its thickness, a smaller volume means a shorter runout of the flow (not visible for the domain size of Fig. 4).

Influence of river constriction on surge genesis

If a river narrows in a steady state configuration, the volumetric rates of pyroclastic flow are equal before and after the constriction:

w w W W Q u h w u h W   [25]
where w and W are the valley widths and uw, hw, uW and hW are the flow velocity and thickness before and after the constriction, respectively. The thickness and/or the velocity need then to increase to compensate for the reduction in width. With a purely plastic rheology, we have seen that a flow tends to a given thickness for which the driving stress equals the retarding stress and the velocity reaches a steady state. Immediately after the constriction, the flow tends to accelerate to reach the same thickness as before the constriction and hw=hW. The velocity after the constriction is then related to the velocity before the constriction by:

W w u W u w  [26]
The mass flux from the flow to the surge depends on the velocity to the power of three. Then, after the constriction, the mass flux will increase by a factor 3 W w

     
. For example, a reduction of the drainage width by 2 will change the mass flux by 8. Conversely, widening of a valley will decrease the surge genesis.

The effect of a drainage constriction on surge production is simulated by the model in Fig. 5ab. Downstream of the constriction, the valley is sinusoidal (Fig. 4a). Upstream, the central line of the drainage has been enlarged by a 100-m wide incline.

Before the constriction, the velocity calculated by the model with the parameters chosen is about 8 m/s, whereas afterwards it is about 30 m/s. With the parameters chosen (Tab. 1) the mass flux m changes from 5 to 270 kg/m²/s. The mass flux of 5 kg/m²/s is not high enough to compensate the mass lost by particle settling and a very thin surge is created. However, after the constriction, the mass flux of 270 kg/m²/s is larger than the mass lost by particle settling and a surge can form. The surge leaves the drainage to flow over several tens of meters on the interfluves. Thus, the model shows that variations in valley width can have a strong influence of the surge genesis and dynamics.

Decoupling between the currents

Concentrated flows are more sensitive to meter-scale topography than thick, dilute currents.

Thus decoupling between the two flows often occurs at river bends. The following topography was used to check if the model could reproduce this decoupling: a sinusoidal valley cross section, with the first part of the valley linear and parallel to the slope (Fig. 5c), after which the valley curves by 45° before returning to the downslope orientation. The curve follows a circle with a radius of 100 m. Using the chosen parameters (Tab. 1), the concentrated flow formed is about 2.5 m thick with a velocity of ~15 m/s. The mass and volume of the surge increase progressively (max. thickness ~ 25 m). When the current reaches the bend, the surge continues straight forward onto the interfluves, apart from a small part which follows the drainage (Fig. 5d). The concentrated part follows the drainage because its thickness and velocity at the bend, with the parameters chosen, are too low to flow onto the interfluves. It is thus decoupled from the main surge. After the bend, the flow still generates the surge but, due to the lower flow velocity related to its thinning, the new surge formed is thinner and only overflows the interfluves immediately adjacent to the main drainage (Fig. 5e). As shown in figure 5e the thickest surge deposits are found downstream and along the projection of the upstream valley. shows that the densities of the surge and of the mixture play a secondary role in the dynamics of the surge: the shape of the deposit is similar for a surge density of about 10 kg/m 3 and 20 kg/m 3 . This can be explained by Eqs. 7-8 which show that the main parameter that governs the surge dynamics is the local mass (s hs). A thick, light surge has the same mass as a thin, concentrated surge. This conclusion is not valid if surge thickness is close to the depth of the valley: surge can be strongly influenced by the valley shape and can be confined above the concentrated flow. The conclusions are not valid either for very light surges, which rapidly reach the density of the atmosphere. The results justify a posteriori that air entrainment can be considered of second importance as long as the surge density is significantly higher than the atmosphere density.

Influence of the mass rate

Since surge genesis is related, in the model, to the velocity of the concentrated flow, any parameter that causes an increase in the flow velocity also increases the surge production. As explained in the previous sections, a thick flow moves faster than a thin flow. For a given volume of flow, a short duration of formation, and consequently a high volumetric rate, forms thicker flows than low volumetric rates. Thus high volumetric rates are related to faster flows and to more powerful surges.

Simulations carried out with various volumes and volumetric rates show that, close to the source, the principal influence is the mass rate, which controls the velocity of the flow. The areas affected by the surges are similar whatever the volumetric rate and are poorly dependent on the volume. Downstream, however, the velocity reaches an equilibrium. It is no longer influenced by the initial velocity and thus by the mass rate. Here, the velocity and the runout are mainly controlled by the volume of the flow, with a larger volume allowing the concentrated part to flow and to create a surge far from the source.

Change of surge extension with valley infill

Field observation indicate that the valley infill, which decreases the valley's depth, seems to play a role in the destruction of the surrounding interfluves by ash-cloud surges [START_REF] Komorowski | Paroxysmal dome explosion during the Merapi 2010 eruption: Processes and facies relationships of associated high-energy pyroclastic density currents[END_REF][START_REF] Ogburn | The effect of topography on ash-cloud surge generation and propagation[END_REF]. To evaluate the effect of the valley depth on the surge emplacement, simulations were done with the same geometry as for the simulation in Fig. 4b (a 20-m-deep valley, slope 5°) but with the addition of a 50-m high plateau bounding the drainage on half of the calculation domain (closer to the source, Fig. 5f). Remind that such a topography creates high curvatures that lower the model accuracy. This is, however, interesting to explore the model behavior. In the valley surrounded by the plateau, the surge formed with a mixing density m of 10 kg/m 3 is too thin to flow out of the 50-m-deep valley. The surge remains essentially channeled above the flow and particles settle onto it. With more width between the plateau edges and the drainage, the particles can form a thick deposit at the sides of the flow and the model shows that the presence of a surrounding topography favours particle accumulation and, consequently, the formation of surge-derived concentrated flows.

At the downstream end of the deep part, the surge is able to spread out onto the interfluves (Fig. 5f). However, because the surge has settled progressively, the surge thickness and density are not very different to what they would have been without the plateau. With a density m of 2 kg/m 3 , the same mass flux forms a thicker surge that can overflow the 50-m walls of the drainage. However, the edges that bound the flow still play a major role in the surge emplacement, and the distance reached by the surge is strongly reduced (150 m from the valley on the plateau, and 300 m from the valley elsewhere).

Discussion

Compatibility of the model with field observations

Our understanding of the coupling between concentrated and dilute currents as well as the physics of these currents remains very incomplete, as pointed out by [START_REF] Ogburn | The effect of topography on ash-cloud surge generation and propagation[END_REF]. A way to improve our understanding of the physics is to make models, and try to reproduce field observations as closely as possible by changing the model parameters, and then by interpreting the values and the laws obtained.

Ash-cloud surges are known to be more powerful and to overrun the concentrated flow after traversing a break in slope. For example, this has been described at Santiguito, Guatemala [START_REF] Rose | The evolution of Santa Maria Volcano, Guatemala[END_REF], at Ngauruhoe [START_REF] Nairn | Explosive eruptions and pyroclastic avalanches from Ngauruhoe in February 1975[END_REF] and at Merapi, Indonesia [START_REF] Kelfoun | A statistical study of trees damaged by the 22 November 1994 eruption of Merapi volcano (Java, Indonesia): relationships between ash-cloud surges and block-and-ash flows[END_REF][START_REF] Lube | Kinematic characteristics of pyroclastic density currents at Merapi and controls on their avulsion from natural and engineered channels[END_REF][START_REF] Jenkins | The Merapi 2010 eruption: An interdisciplinary impact assessment methodology for studying pyroclastic density current dynamics[END_REF]. The tree directions measured following the 1994 eruption of Merapi [START_REF] Kelfoun | A statistical study of trees damaged by the 22 November 1994 eruption of Merapi volcano (Java, Indonesia): relationships between ash-cloud surges and block-and-ash flows[END_REF] are entirely compatible with the results of Fig. 4.

Surges are also known to detach from the concentrated flow at valley bends [START_REF] Yamamoto | Pyroclastic flows from the 1991 eruption of Unzen volcano, Japan[END_REF][START_REF] Fujii | The 15 September 1991 pyroclastic flows at Unzen Volcano (Japan): a flow model for associated ash-cloud surges[END_REF][START_REF] Kelfoun | A statistical study of trees damaged by the 22 November 1994 eruption of Merapi volcano (Java, Indonesia): relationships between ash-cloud surges and block-and-ash flows[END_REF][START_REF] Charbonnier | Field observations and surface characteristics of pristine block-and-ash flow deposits from the 2006 eruption of Merapi Volcano, Java, Indonesia[END_REF][START_REF] Lube | Kinematic characteristics of pyroclastic density currents at Merapi and controls on their avulsion from natural and engineered channels[END_REF][START_REF] Komorowski | Paroxysmal dome explosion during the Merapi 2010 eruption: Processes and facies relationships of associated high-energy pyroclastic density currents[END_REF]. [START_REF] Loughlin | Eyewitness accounts of the 25 June 1997 pyroclastic flows and surges at Soufrie`re Hills Volcano, Montserrat, and implications for disaster mitigation[END_REF] described the formation of surges after valley constrictions. A strong negative correlation between the valley cross-sectional area and the area covered by detached surges has been demonstrated by [START_REF] Ogburn | The effect of topography on ash-cloud surge generation and propagation[END_REF]. Our model reproduces this effect (Fig. 5b). Their observation that a critical cross-section exists in reality is compatible with the critical velocity needed by the model to form a surge (~20 m/s).

The critical section of [START_REF] Ogburn | The effect of topography on ash-cloud surge generation and propagation[END_REF] might then be explained by the velocity needed to allow a given volumetric rate to flow across a given valley width. [START_REF] Cole | Deposits from dome-collapse and fountain-collapse pyroclastic flows at Soufrière Hills Volcano, Montserrat[END_REF] studied the pyroclastic flows of Montserrat and concluded that there is no direct relationship between the volume of pyroclastic flows and the area covered by surges. This is consistent with the results of the model that the volumetric rate plays a more crucial role in the surge genesis than the volume itself.

Influence of the parameters on surge dynamics

The model has been developed to be as simple as possible. Nevertheless, because it aims to capture the first order behavior of a natural complex phenomenon, about 20 parameters need to be defined. Some of them are known, or can be calculated or estimated in the field (e.g. volume of the flow, gravity, concentrated flow density, topography). But five parameters cannot be directly estimated: the rheological parameters of the concentrated flow (T0 and a1, Eq. 4), the turbulent coefficient of the surge (a2, Eq. 9), the parameter a3 (Eq. 10) that governs the exchange between the concentrated flow and the surge, and the density of the mixture that supplies the surge (m,Eq. 11). This section explains how they can be estimated from field observations. T0 and a1 control the flow emplacement and can be estimated relatively easily. T0 controls the thickness of the concentrated deposit. A high value forms thick, short deposits, a low value forms thin, long deposits. The parameter a1 influences the concentrated flow velocity, and an increase in this value slows down the simulated flows. On a complex topography, the velocity and thus the value of a1 can be estimated depending on whether the flow can cross a given morphology. The results are strongly influenced by the choice of the plastic rheology as a proxy for the complex rheology of pyroclastic flows. The advantage of the plastic rheology compared to the Coulomb law is that is reproduces the correlation between the flow thickness and it flow capacity. This is probably why it successfully reproduces the emplacement of natural pyroclastic flows [START_REF] Kelfoun | Is frictional behaviour suitable for pyroclastic flow simulation: Comparison with a well constrained eruption at Tungurahua volcano (Ecuador)[END_REF][START_REF] Charbonnier | Evaluation of the impact of the 2010 pyroclastic density currents at Merapi volcano from high-resolution satellite imagery, field investigations and numerical simulations[END_REF] as well as debris avalanches [START_REF] Kelfoun | Numerical modeling of the emplacement of Socompa rock avalanche, Chile[END_REF][START_REF] Kelfoun | Suitability of simple rheological laws for the numerical simulation of concentrated pyroclastic flows and long-runout volcanic avalanches[END_REF]. Other rheological laws that would reproduce this correlation would lead to similar results.

The variable a2 controls the surge velocity. It can be estimated using direct observations of the velocity or from the destruction, related to this velocity and to the surge density. The ratio of the length of the surge to its width can be also used to define a2. A low value of a2 tends to increase the length of the surge downstream.

Variable a3 determines the mass flux from the concentrated flow and is estimated from the mass of surge deposits. With a high value of a3, a large mass of particles from the concentrated flow is passed to the surge, and voluminous surge deposits are formed. This parameter implicitly takes into account phenomena that act on the surge genesis (temperature, pore pressure, etc.). Future work must define this dependence more quantitatively, must better explain the mechanism responsible of the particle exchanges, which is certainly more complex than saltation, and must confirm or improve Eq. 10.

Mixture density m determines the thickness of the surge and, consequently, its density. It can be estimated by the elevation reached by the surge along the hills, by the capacity of the surge to flow over the topography or from the estimation of the surge density (using the dynamic pressure for example, [START_REF] Jenkins | The Merapi 2010 eruption: An interdisciplinary impact assessment methodology for studying pyroclastic density current dynamics[END_REF].

Conclusion

A new model has been developed to simulate concentrated pyroclastic flows, pyroclastic surges and their interactions, using two coupled depth-averaged approaches. The analysis of the model behavior shows that the approach and the equations used make it possible to reproduce the behavior of the natural current, at least qualitatively. The model mimics the genesis of surges from concentrated flows and the formation of deposits, and of surge-derived concentrated flows through sedimentation of the surge. It shows that a velocity of about 25 m/s is needed for a concentrated flow to form a surge. All the parameters that increase the concentrated flow velocity can thus be lead to surge formation. This explains the strong influence of the topography on the surge genesis: a break in slope or a river constriction can dramatically influence the surge genesis by accelerating the velocity of the concentrated flow. The model also shows that the volume rate can be more important for the surge genesis than the volume itself. The strong correlation between the volumetric rate and surge genesis can be particularly problematic for hazard assessment. In the case of dome collapse, for example, knowing the dome volume would not be enough in itself to estimate the areas covered by resultant surges.

The volumetric rate is a critical parameter and depends on the way the dome collapses, which is difficult to predict, and could explain the sudden genesis of powerful surges, the difficulty of predicting their subsequent damage, and the potential danger for the population. The model must now be tested against real field cases in order to check its ability to quantitatively reproduce natural phenomenon. In a companion paper [Kelfoun et al., submitted], we show that

Figure

  Figure 2a shows the influence of the surge density and of the flow velocity on the mass balance that controls the surge evolution. The bold curve shows the equilibrium between the mass lost by sedimentation and the mass gained from the concentrated flow. To the right of the bold curve, values are positive and indicate that the mass elutriated from the concentrated flow is higher than the mass lost by sedimentation so that the surge mass increases. To the left, the mass balance is negative: the mass flux from the concentrated flow is insufficient to compensate the mass flux lost by sedimentation and either the concentrated flow cannot form a surge or an existing surge loses mass until it disappears. The density of the surge has two antagonistic

Figure

  Figure 2b shows the same mass transfer rates according to the velocity of the concentrated flow and to different particle sizes. Negative values mean that the surge loses mass and tends to disappear. The null flux, where the mass lost equals the mass gained, is represented in bold for a surge density of 10 kg/m 3 . Two other null-flux curves for a surge density of 1.1 kg/m 3 and 2 kg/m 3 are represented by the dashed and the dotted lines, respectively. Fig. 2b also demonstrates the existence of a threshold in surge formation at around 15-25 m/s. Below this velocity, the

  the surge viscosity. For Stk<<1, particles closely follow fluid streamlines and can be mixed homogeneously. Particles smaller than 0.1 mm in diameter, as

3. 7 .

 7 Influence of the density of the surge ( s ) and of the mixing ( m ) Simulation were carried out using the same parameters as for Fig 5c-e, the only difference being the value of m, the density of the gas-particle mixture that leaves the concentrated flow to form the surge. The value of m controls the value of the surge density s. However, it does not affect the mass flux which remains constant. A low value of m thus forms a thick, dilute surge while a high value forms a thin, dense surge. Comparing the two results (see supplementary material)

Figure 1 :

 1 Figure 1: Flow diagram of the model. The concentrated flow loses a mass flux s of densitym that forms the dilute current (surge). The dilute current is of variable density and loses a mass flux s that can form a concentrated flow or a deposit.

Figure 2 :

 2 Figure 2: (a) Mass fluxes at the interface between the concentrated flow and the surge (in kg per m² and per second, values in italic), based on the velocity of the concentrated flow and the surge density. The mass flux to the surge is given by Eq. 10-12 while the mass is lost according to Eq. 14 and 18. The solid curves are drawn for particle diameters of 0.25 mm. The bold curve indicates that the mass gained balances the mass lost (mass balance = 0). Dashed curves correspond to the equilibrium state for particles of 0.1 mm and 1 mm in diameter. The dark grey zone corresponds to the area where the dilute mixture is lighter than the atmosphere and lifts off. (b) Mass fluxes between the concentrated flow and the dilute current (in kg per m² and per second) based on the velocity of the concentrated flow and the particle diameter (log scale) for a surge density of 10 kg/m 3 . Dashed line: equilibrium state for a density of 2 kg/m 3 and 1.1 kg/m 3 .

Figure 3 :

 3 Figure 3: (a-c) Surge genesis from a concentrated flow. Once generated, the surge has its own dynamics and can overrun the concentrated flow. (a) t = 5 s, (b) t = 15 s. (b-d) Genesis firstly of a concentrated deposit and then of a concentrated flow from the mass lost by the surge. The surge is generated continuously at the left-hand border of the domain with a thickness of 50 m and a density of 20 kg/m 3 . (a) t = 10 s, (b) t = 200 s.

Figure 4 :

 4 Figure 4: Effect of the topographic slope on the area covered by pyroclastic surges. (a) Cross section of the valley. (b) With a slope of 5° the surge extends up to 70 m from the concentrated flow. (c) If a steeper slope is added to the topography, the surge genesis is much stronger. Movie in supplementary material.

Figure 5 :

 5 Figure 5: (a-b) Downstream of a constriction, the higher velocity of the concentrated flow creates more voluminous surges that can leave the main drainage to overflow the interfluves. (c-e) Decoupling between the concentrated and the dilute current. Images show deposit extent for concentrated and dilute currents at (a) t=30 s, (b) 45 s and (c) at rest (150 s). Movie in



  x ,  y topographic slope in x and y  x ,  y slope of the flow surface in x and y T x , T y ] total resisting stress for the flowR = [R x , R y ]total resisting stress for the surge u= [u x , u y ] flow velocity v = [v x , v y ]surge velocity w sed settling velocity of surge particles xx-direction along the slope y y-direction along the slope

Table 1 :

 1 Parameters used in the simulations. For all the simulations, the value of a is set at 1 kg/m 3 . The values of the particle diameters dp are given for Cd = 1. For other values of Cd (e.g.[START_REF] Dellino | A method for the calculation of the impact parameters of dilute pyroclastic density currents based on deposit particle characteristics[END_REF] the particle diameter must be changed to keep the ratio dp/Cd constant. g, gravity; g, density of the gas in the surge; s, surge density; p, particle density; d, concentrated flow density; m, mixing density (section 2.3); a1, turbulent/collisional term for the concentrated flow; a2, turbulent/collisional term for the surge; a3, term for surge genesis; dp, particle diameter; V, volume of the concentrated part, ts, duration of current genesis; _c, the value changes and is calculated by the model, _, parameter not used. All the models described in this article were run in a few minutes on a desktop computer (maximal domain size: 200×500 meshes, calculation speed: ≈1 s of calculation for 1 s in reality on a desktop computer, ≈15 min per simulation).

		g	 g	 s	 p	 d  m T 0	a1	a2	a3	d p	d p	V	t s	min( ) max( )
		m/s²			kg/m 3			kPa	-	-kg s 2 m -5	mm	 unit	m 3	s	-	-
	Fig. 2a - =f ( s , u )	9.81 0.5 0-20 2000	-	-	-	-	-	0.001	1	0	-	-	0	0
	Fig. 2b - =f (d , u )	9.81 0.5	10	-	-	-	-	-	-	0.001 0.01-100 -6.6 -6.6	-	-	0	0
	Fig. 3ab -Surge from flow	9.81 0.64 _c 2400 1600 10 20 0.01 0.2	0.01	10	-3.32	1334 m 3 /m	0	0.38	45
	Fig. 3cd -Flow to surge	9.81 0.64 _c 2400 1600 10 10 0.01 0.2	0.01	10	-3.32	~2600 m 3 /m 200 s	0.38	45
	Fig. 4 -break in slope	9.81 0.64 _c 2400 1600 10	4	0.01 0.2	0.01	0.0625	4	100000	60	5	30
	Fig. 5ab -constriction	9.81 0.64 _c 2400 1600 10	4	0.01 0.2	0.01	0.0625	4	100000	60	10	10
	Fig. 5c-e -bend,  mix 10	9.81 0.64 _c 2400 1680 10	4	0.01 0.1	0.01	0.0625	4	40000	60	20	20
	Fig. 5f -in a deep valley	9.81 0.64 _c 2400 1600 10	4	0.01 0.2	0.01	0.0625	4	200000	30	5	5

Table 1 :

 1 List of variables used

Acknowledgement

This work has be founded by the Domerapi -ANR (French Agence Nationale de la Recherche)

Project (ANR-12-BS06-0012). The is the LabEx ClerVolc contribution n° 249. The data for this paper are available by contacting the corresponding author at k.kelfoun@opgc.fr. I am very grateful to both Dr. Tomaso Esposti Ongaro and an anonymous reviewer for their constructive remarks.

1

Journal of Geophysical Research -Solid Earth

Supporting Information for

A two-layer depth-averaged model for both the dilute and the concentrated parts of pyroclastic currents

Karim Kelfoun

Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, F-63000 Clermont-Ferrand, France.

Contents of this file

Figures S1

Additional Supporting Information (Files uploaded separately)

Captions for Movies S2 to S3

Introduction

The supporting information presents an additional simulation and two movies whose caption are in the manuscript.

2


Yvan Fournier

Emmanuel Jeannot

Benjamin Lorendeau

Experiments with multi-level parallelism runtimes on a CFD code with unstructured meshes

Experiments with multi-level parallelism runtimes on a CFD code with unstructured meshes

Yvan Fournier, EDF R&D (MFEE) Emmanuel Jeannot, INRIA (TADAAM)

Benjamin Lorendeau EDF R&D & INRIA

Abstract

Large applications for parallel computers and more specifically unstructured Computational Fluid Dynamics codes are often based on the bulk synchronous parallelism approach (BSP model) and therefore mostly exploit parallelism using runtime systems like the Message Passing Interface (MPI) which has enabled strong, relatively portable and quite durable performances for these codes for many years. However, as MPI was developed for distributed computing, we do not expect its standalone use to be the best fit for recent many-core architectures.

Indeed, the ever growing performance of high performance machines we are witnessing comes mostly through the aggregation of different computing devices in order to build heterogeneous machines. Those are the combination of traditional computing units (CPUs) with accelerators, namely many-core architectures such as GPGPUs or Intel Xeon-Phi accelerators. This hybridization of HPC clusters requires current scientific code developers to master more and more techniques and programming models in order to harness the quintessence of their machines.

As the post-petascale era has long been foreseen, runtime systems developers have been investigating other parallelism paradigms. Notably, the task based approach has gained a lot of popularity recently as it is expected to deliver portability and performance with a relatively simple programming model. Tasks can be both local or distant, so a single model can handle both internode and intra-node aspects. In addition, computation-communication overlap is straightforward. Opposed to the Coarse Grain Parallelism, performance depends strongly on choosing a good data grain size for each task, which should require performance measuring and tuning but no additional programming effort.

As Code_Saturne [START_REF] Archambeau | Code saturne: A finite volume code for the computation of turbulent incompressible flows-industrial applications[END_REF], our CFD code at EDF R&D is based on unstructured meshes, with a significant part of its code being memory-bound, refined parallelism through the use of MPI + X solutions such as MPI + OpenMP often fails to deliver significant (or any) performance improvements, though it does reduce the memory consumption per thread. Using a simple "loop-local" OpenMP model, as we increase the number of threads per MPI rank, performance drops rapidly, since many secondary loops are not threaded; and avoiding data races often requires specific renumbering strategies, which may not be easily adapted everywhere with a reasonable programming effort. These diminishing returns tend to limit the efforts which are worthwhile to spend in addition to the base MPI model.

We may see HPC current technologies evolution as unfavorable to an unstructured CFD code like Code_Saturne in its current form. This is why we decided to investigate recent HPC techniques and runtime systems for a sustainable, easy to propagate, portable and efficient solution to bring better performance and adaptability to Code_Saturne.

As many teams are already dedicating their work to propose new solvers and solve dense linear algebra, we decided to focus on another part of the puzzle, namely our gradient reconstruction computation. As a significant portion of our main current numerical schemes, it has a high impact over the performance of our code and an intermediate computational intensity. As such, we propose in this article a review of different implementations of our gradient computation towards the implementation of a task based approach through the use of task-based HPC runtime systems, and more specifically the PaRSEC [START_REF] Bosilca | A generic distributed {DAG} engine for high performance computing[END_REF] framework. The Parallel Runtime Scheduling and Execution Control (PaRSEC) framework implements a task-based dataflow-driven programming model aimed at offering high performance while relieving developers of supercomputers' hardware complexity.

We show that our first implementation offers comparable performance while increasing the arithmetic intensity of its computation (see figure 1). Moreover, by removing some data dependencies, our cell based approach paves the way for a more refined grain parallelism approach. We then push this approach to our prime objective -task based approaches -and implement our gradient computation using the PaRSEC runtime. Finally, we propose some insights on the use of We observe roughly same performances while arithmetic intensity of cell based approach is more than the double of the face based approach.

such runtime systems over more common HPC techniques based on their ease of their performance scalability, ease of incremental deployment in a legacy CFD code, and maintainability.

Figure 1 -

 1 Figure1-Computation time comparison between face based and cell based lsq gradient (lower is better). We observe roughly same performances while arithmetic intensity of cell based approach is more than the double of the face based approach.