Solène Bulteau
email: solene.bulteau@univ-nantes.fr

Christophe Berthon
email: christophe.berthon@univ-nantes.fr

Marianne Bessemoulin-Chatard
email: marianne.bessemoulin@univ-nantes.fr

An Asymptotic Preserving scheme for the shallow-water equations with Manning friction using viscous correction of the HLL scheme

Keywords: Finite volume schemes, asymptotic-preserving schemes, conservation laws with source term, shallow-water equations AMS subject classifications. 65M08, 35Q35

The aim of this paper is to propose a numerical scheme for the shallowwater equations with a Manning's friction source term able to preserve the diffusive regime, namely the behavior of the solutions in long time and stiff friction limit. Because of the non-usual rescaling of the friction parameter, a strongly nonlinear derivative operator is involved in the diffusive limit and the development of a numerical scheme possessing this property is not straightforward. The scheme presented here is based on a perturbed HLL scheme with additional viscous terms.

Introduction

This paper is devoted to the development of an asymptotic preserving scheme, based on a perturbed HLL discretization, to approximate the solutions of the shallow-water equations with a Manning's friction source term [START_REF] Manning | On the flow of water in open channels and pipes[END_REF]. The system of interest reads:

     ∂ t h + ∂ x q = 0, ∂ t q + ∂ x q 2 h + gh 2 2 = -k|q|qh -η , (x, t) ∈ R × R + , ( 1
)
where the unknowns, h(x, t) > 0 and q(x, t) ∈ R, both depending on space and time, respectively denote the water height and the discharge. As usual, the water velocity is

given by u = q/h. In the model, g = 9.81 m.s -2 is the gravity constant while k and η denote friction parameters according to the Manning friction model [START_REF] Manning | On the flow of water in open channels and pipes[END_REF]. In general, η is given by 7/3. Concerning the parameter k, called Manning coefficient, it defines the intensity of the friction to be adopted.

Here, we are interested in the behavior of h and q in long time and dominant friction, and faraway from dry areas. Since our aim is to build a numerical scheme able to preserve the correct asymptotic behavior of the solutions, we formally recall the derivation of the asymptotic diffusive regime. A small parameter ε is thus introduced in order to govern the time t and the friction coefficient k as follows:

t ← t/ε and k ← k/ε 2 .
(2)

The friction parameter has to be rescaled differently than the time because of the quadratic term in q in the source term (see [5,[START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF]). With this scaling, the system (1) now reads:

     ε∂ t h + ∂ x q = 0, ε∂ t q + ∂ x q 2 h + gh 2 2 = - k ε 2 |q|qh -η .
(3)

To study the behavior of h and q when ε tends to zero, the following Chapman-Enskog expansions of the unknowns are introduced: h = h 0 + εh 1 + . . . and q = q 0 + εq 1 + . . . .

The here considered water height h is assumed to be positive but, thereafter, we also have to impose the positivity of its first-order term, namely h 0 > 0. Plugging expansions (4) in the rescaled system (3), enforcing ε to zero, we obtain the following limit system:

∂ x q 0 = 0,
-k|q 0 |q 0 (h 0 ) -η = 0.

(5)

Since h 0 > 0, the second expression of the above system necessarily gives q 0 = 0. Then, from (3), we now obtain the following limit problem:

     ∂ t h 0 + ∂ x q 1 = 0, ∂ x g(h 0 ) 2 2 = -k|q 1 |q 1 (h 0 ) -η . ( 6
)
The second expression of the above system leads to the following definition of q 1 , called the local equilibrium (see [20]):

q 1 = -sign(∂ x h 0 ) (h 0 ) η k ∂ x g(h 0 ) 2 2 . ( 7
)
Next, considering the h 0 governing evolution equation in [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF], we obtain the following nonlinear diffusion equation:

∂ t h 0 + ∂ x -sign(∂ x h 0 ) (h 0 ) η k ∂ x g(h 0 ) 2 2 = 0. (8)
In the above equation, we may recognize a non stationary p-laplacian type equation supplemented with a nonlinear flux function. Let us underline that the p-laplacian equation, known to be a degenerate parabolic equation, has been widely studied (see [2,19,23,[START_REF] Segura De León | Regularity for entropy solutions of parabolic p-Laplacian type equations[END_REF] for theoretical aspects and [START_REF] Andreianov | Finite volume schemes for the p-Laplacian on Cartesian meshes[END_REF][START_REF] Barrett | Finite element approximation of the p-Laplacian[END_REF]22,[START_REF] Peton | Étude et simulation d'un modèle stratigraphique advecto-diffusif nonlinéaire avec frontières mobiles[END_REF] for some numerical studies) and appears in several physical problems as, for instance, in non-newtonian fluids [12].

The aim of this work is to derive a scheme for [START_REF] Andreianov | Finite volume schemes for the p-Laplacian on Cartesian meshes[END_REF] such that, in the diffusive regime, the approximate solution is consistent with the limit problem [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF], or equivalently (7)-(8). This property to be satisfied by the numerical scheme, called asymptotic preserving (AP) [5-9, 11, 14, 17, 20], is not straightforward and the development of AP-schemes requires a particular attention. During the two last decades, numerous AP-schemes were proposed in the literature. For instance, in [3,[START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF][START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF], AP numerical strategies have been introduced to accurately approximate the asymptotic diffusive regime issuing from the radiative transfer models. We also refer to [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF] where an AP-scheme is derived to approximate the solutions of the isentropic gas dynamics in the Darcy law regime. Next, in [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF], a generic formulation of AP-schemes is proposed by a suitable extension of the well-known HLL scheme [START_REF] Chalons | Godunovtype schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[END_REF]. At the wide discrepancy with the above mentioned works where the source term rescaling is governed by 1/ε, according to (2), in the present paper we have to deal with a rescaling prescribed by 1/ε 2 . High order source term rescaling in 1/ε m are considered in [5] where the asymptotic diffusive regime is established to be strongly nonlinear. Of course, this nonlinearity is clearly stated in (8). In [5], the authors concluded the asymptotic diffusion study by the derivation of a generic numerical scheme. More recently, in [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF], the shallow-water model with a Manning friction given by ( 1) is adopted and a scheme with interesting asymptotic behavior is derived. However, the diffusion regime satisfied by the introduced method comes from an extension of the numerical technique introduced in [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF] designed for a scaling of the stiffness given by 1/ε. As a consequence, the authors established a correct numerical asymptotic behavior of the water height. But, in [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF], the expected regime verified by the discharge, according to (8), is no longer preserved. In this sense, the numerical scheme proposed in [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF] is at most partially asymptotic preserving. The present work is motivated in the design of a numerical scheme able to fully restore the expected asymptotic diffusive regime given by both relations [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF] and (8).

The paper is organized as follows. Section 2 is dedicated to the presentation of the considered general framework. It allows to introduce several notations useful in the next sections and to present a generic form for all the schemes presented in this paper. Afterwards, in Section 3, the scheme given in [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF] for the system under concern (1) is recalled. The failure of this extension of the scheme presented in [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF] is then highlighted. Indeed, this scheme is able to preserve terms of order zero, namely h 0 given by (8), but not of order one. As a consequence, the relation [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF] to govern q 1 is not preserved and the resulting scheme is said a priori partially asymptotic preserving. In fact, by studying the asymptotic behavior of the discharge, since the scheme adopts a nonphysical scaling in 1/ε, we exhibit an additional relation satisfied once again by h 0 . This new relation combined with (8) make this scheme not asymptotic preserving. This failure is illustrated with several numerical experiments. In Section 4, an extension of the generic method introduced in [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF], respecting the natural rescaling of the source term stiffness, is proposed. This generalization is able to deal with stiffness rescaled in 1/ε m with m ∈ N * . From now on, we emphasize that the here derived scheme turns out to be a revisit of the numerical approach introduced in [5]. Next, this proposed method is applied on the model [START_REF] Andreianov | Finite volume schemes for the p-Laplacian on Cartesian meshes[END_REF], for which m = 2, and the diffusive limit regime, given by ( 7) and (8), is proved to be preserved by the here designed numerical technique. More precisely, the respect of the natural rescaling allows to catch consistent discretizations of the order zero, to get (8), and the order one, to get [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF], of the solutions in the diffusive regime. Finally, Section 5 is devoted to the numerical results to illustrate the interest of the developed scheme.

A Godunov-type scheme with source term

For the sake of completeness in the forthcoming developments, we briefly detail the asymptotic preserving numerical technique introduced in [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF]. To address such an issue, we consider hyperbolic systems with source terms as follows:

∂ t W + ∂ x F (W ) = k (R(x, W ) -W ) , (x, t) ∈ R × R + , ( 9
)
where W ∈ Ω is the vector of the unknowns, with Ω ∈ R N the phase space, F (W ) ∈ R N is the physical flux and R(W ) ∈ R N is a function occurring in the source term. We immediately notice that the system (1) enters in such a formalism with

W = t (h, q) , F (W ) = t q, q 2 h + g h 2 2 and R(x, W ) = t h, q -q|q|h -η . ( 10
)
For the sake of simplicity, we assume k to be a constant and we refer to [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF] for more complex interactions.

Next, we introduce a free parameter, denoted k, such that the system (9) equivalently reformulates as follows:

∂ t W + ∂ x F (W ) = k + k R(x, W ) -W , (x, t) ∈ R × R + , (11)
with

R(x, W ) = k k + k (R(x, W ) -W ) + W. ( 12
)
In fact, from a numerical point of view, k will play the role of a correction to be prescribed in order to recover the expected numerical asymptotic regime.

First, we consider the space discretization made of cells (x i-1/2 , x i+1/2 ) i∈Z with a constant size ∆x such that x i+1/2 = x i-1/2 + ∆x for all i ∈ Z. Next, concerning the time discretization, we adopt t n+1 = t n + ∆t where ∆t > 0 denotes the time increment and where t 0 = 0. As usual, in order to enforce some stability conditions, the time increment will be restricted according to a CFL like condition [START_REF] Courant | On the partial differencial equations of mathematical physics[END_REF]. Now, equipped with (W n i ) i∈Z over each cell (x i-1/2 , x i+1/2 ) at time t n , we search for relevant updated states (W n+1 i ) i∈Z to approximate the solution of (11) at time t n+1 . Following [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF], the required updated states are obtained by adopting a Godunov-type scheme [START_REF] Berthon | Late-time/stiff-relaxation asymptoticpreserving approximations of hyperbolic equations[END_REF]16,[START_REF] Chalons | Godunovtype schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[END_REF]. To address such an issue, we first exhibit an approximate Riemann solver made of four constant states separated by three discontinuities (see Figure 1). According to [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF], the approximate Riemann solver of interest reads as follows:

W x t , W L , W R =                    W L , if x t ≤ -λ LR , α LR W * + (1 -α LR ) R(0 -, W L ), if -λ LR < x t ≤ 0, α LR W * + (1 -α LR ) R(0 + , W R ), if 0 < x t ≤ λ LR , W R , if λ LR < x t , ( 13
)
x 0 | -∆x/2 | ∆x/2 t λ LR -λ LR W L W * L W * R W R
α LR = 2λ LR 2λ LR + k + k LR ∆x . ( 14
)
The main originality in this approximate Riemann solver remains in the introduction of the source term within the two intermediate states separated by a stationary discontinuity. Moreover, for the sake of simplicity in the forthcoming derivations, we here have adopted a symmetric fan made of speed waves -λ LR < 0 < λ LR . Extensions with speed waves λ L < 0 < λ R are straightforward and they are left to the reader. In (13), the vector W * denotes the intermediate state of the classical HLL scheme introduced in [START_REF] Chalons | Godunovtype schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[END_REF], defined by

W * = 1 2 (W R + W L ) - 1 2λ LR (F (W R ) -F (W L )) ,
and the following notation is used in the definition (13) of the Riemann solver:

R(0 ± , W ) = lim x→0 ± R(x, W ).
We underline that, for the specific symetric speeds we consider, the HLL scheme reduces to the Rusanov scheme. In order to obtain the expected updated states, we first put the approximate Riemann solver W (

x-x i+ 1 2 t , W n i , W n i+1 ) at each interface x i+ 1 2 for all i ∈ Z. Next, we define W n+1 i as follows: W n+1 i = 1 ∆x xi x i-1 2 W x -x i-1 2 ∆t , W n i-1 , W n i dx + 1 ∆x x i+ 1 2 xi W x -x i+ 1 2 ∆t , W n i , W n i+1 dx.
With clear notations stated at each interface x i+ 1 2 , after straightforward computations we obtain:

W n+1 i = W n i - ∆t 2∆x α i+ 1 2 F (W n i+1 ) -F (W n i ) + α i-1 2 F (W n i ) -F (W n i-1 ) + ∆t 2∆x α i+ 1 2 λ i+ 1 2 W n i+1 -W n i -α i-1 2 λ i-1 2 W n i -W n i-1 + k∆t 2 α i+ 1 2 R(x - i+ 1 2 , W n i ) -W n i + α i-1 2 R(x + i-1 2 , W n i ) -W n i . ( 15
)
Applying the above generic discretization to the shallow-water equations with Manning friction (1), we obtain

h n+1 i = h n i - ∆t 2∆x α i+ 1 2 q n i+1 -q n i + α i-1 2 q n i -q n i-1 (16) + ∆t 2∆x λ i+ 1 2 α i+ 1 2 h n i+1 -h n i -λ i-1 2 α i-1 2 h n i -h n i-1 , q n+1 i = q n i - ∆t 2∆x α i+ 1 2 (q n i+1 ) 2 h n i+1 + g (h n i+1 ) 2 2 - (q n i ) 2 h n i -g (h n i ) 2 2 ( 17
)
+α i-1 2 (q n i ) 2 h n i + g (h n i ) 2 2 - (q n i-1 ) 2 h n i-1 -g (h n i-1 ) 2 2 + ∆t 2∆x λ i+ 1 2 α i+ 1 2 q n i+1 -q n i -λ i-1 2 α i-1 2 q n i -q n i-1 - k∆t 2 α i+ 1 2 + α i-1 2 q n i |q n i |(h n i ) -η ,
where the wave velocities (λ i+ 1 2 ) i∈Z are fixed as follows:

λ i+ 1 2 = max |q n i | h n i + gh n i , |q n i+1 | h n i+1 + gh n i+1 for all i ∈ Z. ( 18
)
The coefficient α is still defined by [START_REF] Franck | Finite volume scheme with local high order discretization of the hydrostatic equilibrium for the Euler equations with external forces[END_REF], where the parameter k, in order to restore the expected asymptotic preserving property, plays the role of a free parameter and is fixed to catch a discretization of the correct limit in the diffusive regime.

Failure of an asymptotic preserving scheme

In this section, we recall the scheme developed in [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF] and we show that the expected APproperty is not fully satisfied. We adopt the scheme ( 16)-( 17) with the wave velocities [START_REF] Chalons | Godunovtype schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[END_REF]. According to [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF], in order to recover the asymptotic diffusion equation (8), the parameters (k i+ 1 2 ) i∈Z are given by

k i+ 1 2 = 2λ 2 i+ 1 2 k|h n i+1 -h n i | g∆x((h n i+1 ) η + (h n i ) η )(h n i+1 + h n i ) -k for all i ∈ Z. (19)
This correction has been selected in order to get a consistent approximation of the diffusive equation (8) in the asymptotic regime. However, this limit regime is obtained adopting the following rescaling of the time increment ∆t and the friction coefficients k and k LR , using a small parameter ε:

∆t ← ∆t/ε, k ← k/ε and k LR ← k LR /ε.
Once again, we underline that this rescaling is not the natural rescaling of the friction, given by (2). As a consequence, this rescaling is inconsistent with the continuous framework.

After injecting the above rescaling, the scheme ( 16)-( 17) now reads as follows:

h n+1 i = h n i - ∆t 2ε∆x α ε i+ 1 2 q n i+1 -q n i + α ε i-1 2 q n i -q n i-1 (20) + ∆t 2ε∆x λ i+ 1 2 α ε i+ 1 2 h n i+1 -h n i -λ i-1 2 α ε i-1 2 h n i -h n i-1 , q n+1 i = q n i - ∆t 2ε∆x α ε i+ 1 2 (q n i+1 ) 2 h n i+1 + g (h n i+1 ) 2 2 - (q n i ) 2 h n i -g (h n i ) 2 2 ( 21
)
+α ε i-1 2 (q n i ) 2 h n i + g (h n i ) 2 2 - (q n i-1 ) 2 h n i-1 -g (h n i-1 ) 2 2 + ∆t 2ε∆x λ i+ 1 2 α ε i+ 1 2 q n i+1 -q n i -λ i-1 2 α ε i-1 2 q n i -q n i-1 - k∆t 2ε 2 α ε i+ 1 2 + α ε i-1 2 q n i |q n i |(h n i ) -η .
where α, defined by ( 14) is also rescaled to get:

α ε i+ 1 2 = 2λ i+ 1 2 ε 2λ i+ 1 2 ε + (k + k i+ 1 2 )∆x , ( 22
)
and with k i+ 1 2 defined by ( 19). In the following theorem, we exhibit the failure of the above scheme. We only consider the case with a non constant water height for the sake of simplicity. Indeed, a failure in this particular case is enough to invalidate the asymptotic preserving nature of the scheme.

Theorem 3.1. We assume that h n,0 i > 0 for all i ∈ Z and we consider a cell

[x i-1 2 , x i+ 1 2 ] such that h n,0 i ̸ = h n,0 i-1 and h n,0 i ̸ = h n,0
i+1 . Then, in the limit of ε to zero, the first equation of the scheme (20) reads as follows:

h 0,n+1 i -h 0,n i ∆t = 1 ∆x   sign h 0,n i+1 -h 0,n i (h 0,n i+1 ) η + (h 0,n i ) η 2k g 2 (h 0,n i+1 ) 2 -g 2 (h 0,n i ) 2 ∆x (23) -sign h 0,n i -h 0,n i-1 (h 0,n i ) η + (h 0,n i-1 ) η 2k g 2 (h 0,n i ) 2 -g 2 (h 0,n i-1 ) 2 ∆x   ,
while the second equation [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF] reduces to a discretization of

∂ x h 0 = 0. ( 24
)
As expected with the choice of k, the discrete diffusion equation ( 23) is a consistent discretization, when ∆x and ∆t tend to zero, of the diffusion equation (8). The first result (23) thus assesses the asymptotic preserving property of the scheme (20)- [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF] in the sense that the diffusive limit (8) is preserved. This is the meaning of asymptotic-preserving adopted in [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF].

However, here, we also study the behavior of the second equation of the scheme [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF] in the diffusive regime, and the result given by [START_REF] Manning | On the flow of water in open channels and pipes[END_REF] does not coincide with a discretization of the expected local equilibrium [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF]. Moreover, this new statement establishes a constant behavior of (h n i ) i∈Z as ε goes to zero, wich arises in contradiction with the specific non constant water height case we consider.

Proof. The zero-order expansions of (h n i ) n i and (q n i ) n i read

h n i = h 0,n i + O(ε) and q n i = q 0,n i + O(ε).
Moreover, since we consider an interface with a non constant water height, the coefficient α ε is of order ε at interfaces x i-1 2 and x i+ 1 2 :

α ε i-1 2 = O(ε) and α ε i+ 1 2 = O(ε). ( 25
)
Using the above expressions, in the limit of ε to zero, the second equation of the scheme (21) yields to q 0,n i |q 0,n i |(h 0,n i ) -η = 0 for all i ∈ Z. Since h 0,n i > 0, we immediately get q 0,n i = 0. Now, once again with [START_REF] Peton | Étude et simulation d'un modèle stratigraphique advecto-diffusif nonlinéaire avec frontières mobiles[END_REF], the scheme (20) is written as follows in the limit of ε to zero:

h 0,n+1 i = h 0,n i + ∆t 2∆x λ i+ 1 2 lim ε→0 α ε i+ 1 2 ε h 0,n i+1 -h 0,n i -λ i-1 2 lim ε→0 α ε i-1 2 ε h 0,n i -h 0,n i-1
.

(26) Using the definition of α ε i+ 1 2 given by ( 22) and the definition of k i+ 1 2 given by ( 19), we have:

lim ε→0 α ε i+ 1 2 ε = 1 λ i+ 1 2 g((h 0,n i+1 ) η + (h 0,n i ) η )(h 0,n i+1 + h 0,n i ) k∆x|h 0,n i+1 -h 0,n i | . ( 27
)
and the limit scheme (26) coincides with the expected limit (23). Now we prove that the second equation ( 21) converges towards a discretization of (24). Since q 0,n i = 0 and with [START_REF] Peton | Étude et simulation d'un modèle stratigraphique advecto-diffusif nonlinéaire avec frontières mobiles[END_REF], when ε tends to zero, the relation ( 21) reduces as follows:

lim ε→0 α ε i+ 1 2 ε (h 0,n i+1 ) 2 -(h 0,n i ) 2 ∆x + lim ε→0 α ε i-1 2 ε (h 0,n i ) 2 -(h 0,n i-1 ) 2 ∆x = 0.
Moreover, since (27) holds, the above relation can be rewritten

sign h 0,n i+1 -h 0,n i ∆x 1 λ i+ 1 2 ((h 0,n i+1 ) η + (h 0,n i ) η )(h 0,n i+1 + h 0,n i ) 3 h 0,n i+1 -h 0,n i ∆x +sign h 0,n i -h 0,n i-1 ∆x 1 λ i-1 2 ((h 0,n i ) η + (h 0,n i-1 ) η )(h 0,n i + h 0,n i-1 ) 3 h 0,n i -h 0,n i-1 ∆x = 0.
Now, we show that the above discretization is consistent with [START_REF] Manning | On the flow of water in open channels and pipes[END_REF]. To address such an issue, let us consider h 0 (x) a smooth enough function such that h 0,n i = h 0 (x i ) and h 0,n i±1 = h 0 (x i ± ∆x). As a consequence, we easily get:

sign ∂ x h 0 (x i ) + O(∆x) 1 λ i+ 1 2 16(h 0 (x i )) η+3 |∂ x h 0 (x i )| + O(∆x) +sign ∂ x h 0 (x i ) + O(∆x) 1 λ i-1 2 16(h 0 (x i )) η+3 |∂ x h 0 (x i )| + O(∆x) = 0.
The above equality can thus be rewritten as follows:

sign ∂ x h 0 (x i ) + O(∆x) λ i+ 1 2 + sign ∂ x h 0 (x i ) + O(∆x) λ i-1 2 × (h 0 (x i )) η+3 |∂ x h 0 (x i )| + O(∆x) = 0.
Since the wave speeds (λ i+ 1 2 ) i∈Z are positive then, when ∆x tends to zero, the above expression is consistent with the expected limit [START_REF] Manning | On the flow of water in open channels and pipes[END_REF]. The proof is thus achieved.

From now on, let us illustrate the inappropriate behavior of the discretizations given by the scheme (20)- [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF]. To address such an issue, two different initial conditions are considered. They are defined on the spatial domain [-5, 5] and are given as follows.

• Continuous initial condition:

h 0 (x) =          2 if x < -1, 1 2 3 + sin 3πx 2 if -1 ≤ x < 1, 1 else,
and q 0 (x) = 0. (28)

• Discontinuous initial condition:

h 0 (x) = 2 if x < 0, 1
else, and q 0 (x) = 0.

These two initial conditions are displayed on Figure 2. Moreover, the parameter k is set equal to 1, the parameter η is equal to 7/3 and zero-flux boundary conditions are used.

As we can observe in dotted line on Figures 3 and4, the water height in the asymptotic regime is non constant. Thus we are in the framework of the Theorem 3.1 and we expect the scheme (20)-( 21) not to be asymptotic preserving. First, we compare the discretization of the water height given by the equation (20) with the straightforward discretization (23) of the diffusive limit. On Figure 3, we display the approximate water height, at time t = 0.01 and for the continuous initial condition, obtained with the scheme (20)-( 21) for different values of ε, the asymptotic parameter. The dotted line corresponds to the limit water height obtained with (23). The same plots, obtained for the discontinuous initial condition, can be observed on Figure 4. The preservation of the diffusive limit by the scheme (20)-( 21) developed in [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF] cannot be observed, which is in agreement with Theorem 3.1. To complete this illustration, we display on Figures 5 and6 the water height obtained for a fixed value of ε and different number of cells. The water height discretized by the limit scheme ( 23) is still obtained with 800 cells. Moreover, to observe the behavior of the order one of the approximate discharge given by the scheme (20)-( 21), we consider (q n i /ε) n i where (q n i ) n i is the discharge given by [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF]. We compare it with the following discretization of the local equilibrium [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF] :

q n i = -sign h n i+1 -h n i-1 (h n i ) η k g 2 (h n i+1 ) 2 -g 2 (h n i-1 ) 2 2∆x , (30)
where (h n i ) n i is given by the discrete diffusion equation ( 23). This comparison is displayed on Figures 7 and8 for a fixed value of N and different values of ε. We can observe that, for a fixed number of cells, more ε is small and more the order one of the approximate discharge is far from the discrete local equilibrium. This allows us to confirm that the scheme (20)-( 21) is not asymptotic preserving in the sense that it does not preserve both the diffusive limit (8) and the local equilibrium [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF].

Finally, on Figures 9 and 10, we display the order one of the discharge compared to the discrete local equilibrium (30) for a fixed value of ε and different number of cells. Once again, we remark that the scheme (20)-( 21) fails capturing the expected asymptotic behavior. where m ∈ N * . Let us emphasize that m is restricted to 1 in [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF] while we have to fix m = 2 to deal with system (3). To derive the scheme under interest, we extend the formalism (11) in order to correctly mimic the rescaling of the source term in 1/ε m . To address such an issue, we introduce a second free parameter δ > 0 as follows:

∂ t W + ∂ x F (W ) = δ(k + k) R(x, W ) -W , (x, t) ∈ R × R + ,
where R is now given by

R(x, W ) = k δ(k + k) (R(x, W ) -W ) + W.
We notice that δ is a proportional parameter. It will be chosen according to the rescaling in order to restore the correct asymptotic behavior.

The scheme (20)-( 21)-( 22) can then be applied by substituting k +k LR by δ(k +k LR ). The next step consists in injecting into the obtained numerical method the following natural rescaling of the time step ∆t and the parameters k and k LR :

∆t ← ∆t/ε, k ← k/ε m and k LR ← k LR /ε m . (31)
Moreover, the parameter δ is chosen as follows:

δ = ε m-1 .
As usual, the parameter k LR is given at each interface by (19) such that the first equation of the scheme is consistent with the diffusive equation in the limit regime. Now we exhibit the scheme for the shallow-water equations with Manning friction (1). The integer m is equal to 2 and the wave speeds are, once again, defined by [START_REF] Chalons | Godunovtype schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[END_REF]. As mentionned previously, the only modification in the scheme (20)-( 21)-( 22) is the substitution of k + k LR by δ(k + k LR ). This time, the natural rescaling defined by

∆t ← ∆t/ε, k ← k/ε 2 and k LR ← k LR /ε 2 , ( 32
)
can then be adopted and δ fixed to ε to obtain the following rescaled scheme:

h n+1 i = h n i - ∆t 2ε∆x α ε i+ 1 2 q n i+1 -q n i + α ε i-1 2 q n i -q n i-1 (33) + ∆t 2ε∆x λ i+ 1 2 α ε i+ 1 2 h n i+1 -h n i -λ i-1 2 α ε i-1 2 h n i -h n i-1 , q n+1 i = q n i - ∆t 2ε∆x α ε i+ 1 2 (q n i+1 ) 2 h n i+1 + g (h n i+1 ) 2 2 - (q n i ) 2 h n i -g (h n i ) 2 2 (34)
+α ε i-1 2 (q n i ) 2 h n i + g (h n i ) 2 2 - (q n i-1 ) 2 h n i-1 -g (h n i-1 ) 2 2 + ∆t 2ε∆x λ i+ 1 2 α ε i+ 1 2 q n i+1 -q n i -λ i-1 2 α ε i-1 2 q n i -q n i-1 - k∆t 2ε 3 α ε i+ 1 2 + α ε i-1 2 q n i |q n i |(h n i ) -η ,
where

(α ε i+ 1 2
) i∈Z is still given by ( 22). The discrepancy between the above scheme (33)-( 34) and the scheme (20)-( 21), developed in [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF], stays in the introduction of the coefficient δ which modifies the ε-order of the coefficient of the source term discretization. Indeed, we now have a division by ε 3 in (34) and no longer ε 2 as in [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF]. The new parameter δ turns out to be essential to improve the numerical diffusion involved in the scheme in order to deal with the correct scaling (31). We now state the asymptotic behavior satisfied by the scheme (33)-(34). Theorem 4.1. We consider the scheme (33)-( 34)-( 22), where the additional parameter k i+ 1 2 is defined by (19). When ε tends to zero, the two following assertions hold. (i) In the diffusive regime, the equation (34) coincides with a discretization of the local equilibrium [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF].

(ii) In the diffusive regime, the equation (33) coincides with the discrete diffusion equation (23) to discretize (8).

Proof. Let us adopt the following Chapman-Enskog expansions of h n i and q n i for all i ∈ Z, n ∈ N:

h n i = h 0,n i + O(ε) and q n i = q 0,n i + O(ε), ( 35
)
where the order zero of the water height h 0,n i is assumed to be positive for all i ∈ Z and n ∈ N. First, we underline that the following developments easily hold for all i ∈ Z and n ∈ N:

(q n i ) 2 h n i + g (h n i ) 2 2 = (q 0,n i ) 2 h 0,n i + g (h 0,n i ) 2 2 + O(ε).
The discrete discharge evolution law (34) can then be rewritten as follows:

q 0,n+1 i = q 0,n i - ∆t 2ε∆x α ε i+ 1 2 (q 0,n i+1 ) 2 h 0,n i+1 + g (h 0,n i+1 ) 2 2 - (q 0,n i ) 2 h 0,n i -g (h 0,n i ) 2 2 + O(ε) +α ε i-1 2 (q 0,n i ) 2 h 0,n i + g (h 0,n i ) 2 2 - (q 0,n i-1 ) 2 h 0,n i-1 -g (h 0,n i-1 ) 2 2 + O(ε) + ∆t 2ε∆x λ i+ 1 2 α ε i+ 1 2 q 0,n i+1 -q 0,n i + O(ε) -λ i-1 2 α ε i-1 2 q 0,n i -q 0,n i-1 + O(ε) - k∆t 2ε 3 α ε i+ 1 2 + α ε i-1 2 q 0,n i |q 0,n i |(h 0,n i ) -η + O(ε) + O(ε). ( 36
)
Moreover, from the definition of (α ε i+ 1

2

) i∈Z given by ( 22) and the definition of (k i+ 1 2 ) i∈Z given by ( 19), α ε i+ 1 2 can be rewritten as follows:

α ε i+ 1 2 = ε ε + λ i+ 1 2 k∆x|h n i+1 -h n i | g((h n i+1 ) η + (h n i ) η )(h n i+1 + h n i ) -1
.

Using expansions (35), from the above expression we deduce the following limit:

lim ε→0 α ε i+ 1 2 ε =          1 λ i+ 1 2 g((h 0,n i+1 ) η + (h 0,n i ) η )(h 0,n i+1 + h 0,n i ) k∆x|h 0,n i+1 -h 0,n i | , if h 0,n i ̸ = h 0,n i+1 , + ∞, if h 0,n i = h 0,n i+1 .
Even if the limit of α ε i+ 1 2 /ε is not necessarily finite, when ε goes to zero, the dominant term in the equation ( 36) coincides with the source term discretization. As a consequence, we get the following limit for all i ∈ Z:

0 = q 0,n i |q 0,n i |(h 0,n i ) -η . ( 37
)
Since h 0,n i is assumed to be positive, the above expression necessarily gives that the order zero of the discharge q 0,n i vanishes in each cell, and the expansions (35) then rewrite:

h n i = h 0,n i + εh 1,n i + O(ε) and q n i = ε q 1,n i + O(ε) . (38)
Moreover, the following development now holds:

q n i |q n i |(h n i ) -η = ε 2 q 1,n i |q 1,n i |(h 0,n i ) -η + O(ε) . ( 39
)
Now, it is necessary to distinguish the expansion of (α ε i+ 1

2

) i∈Z with respect to the behavior of the water height as follows:

α ε i+ 1 2 =          εα 1 i+ 1 2 + O(ε 2 ) if h 0,n i ̸ = h 0,n i+1 , √ εα 1/2 i+ 1 2 + O(ε) if h 0,n i = h 0,n i+1 and h 1,n i ̸ = h 1,n i+1 , α 0 i+ 1 2 + O(ε) if h 0,n i = h 0,n i+1 and h 1,n i = h 1,n i+1 , (40)
where

α 1 i+ 1 2 = 1 λ i+ 1 2 g((h 0,n i+1 ) η + (h 0,n i ) η )(h 0,n i+1 + h 0,n i ) k∆x|h 0,n i+1 -h 0,n i | , α 1/2 i+ 1 2 = 2 λ i+ 1 2 g(h 0,n i ) η+1 k|h 1,n i+1 -h 1,n i |∆x , α 0 i+ 1 2 =   1 + λ i+ 1 2 2 k|h 2,n i+1 -h 2,n i |∆x g(h 0,n i ) η+1   -1 . ( 41
)
We notice that the quantities α 1 ) i∈Z defined as follows:

β L,ε i+ 1 2 = - ∆t 2ε α ε i+ 1 2 g 2∆x (h n i+1 ) 2 -((h n i ) 2 + k ε 2 q n i |q n i |(h n i ) -η , β R,ε i+ 1 2 = - ∆t 2ε α ε i+ 1 2 g 2∆x (h n i+1 ) 2 -((h n i ) 2 + k ε 2 q n i+1 |q n i+1 |(h n i+1 ) -η .
From (38), ( 39) and (40), β K,ε i+ 1 2 admits the following expansion:

β K,ε i+ 1 2 =              β K,0 i+ 1 2 + O(ε) if h 0,n i ̸ = h 0,n i+1 , 1 √ ε β K,-1/2 i+ 1 2 + O(1) if h 0,n i = h 0,n i+1 and h 1,n i ̸ = h 1,n i+1 , 1 ε β K,-1 i+ 1 2 + O(1) if h 0,n i = h 0,n i+1 and h 1,n i = h 1,n i+1 ,
where we have set

β L,0 i+ 1 2 = - ∆t 2 α 1 i+ 1 2 g 2∆x (h 0,n i+1 ) 2 -((h 0,n i ) 2 + kq 1,n i |q 1,n i |(h 0,n i ) -η , β L,-1/2 i+ 1 2 = - k∆t 2 α 1/2 i+ 1 2 q 1,n i |q 1,n i |(h 0,n i ) -η , β L,-1 i+ 1 2 = - k∆t 2 α 0 i+ 1 2 q 1,n i |q 1,n i |(h 0,n i ) -η , and β R,0 i+ 1 2 = - ∆t 2 α 1 i+ 1 2 g 2∆x (h 0,n i+1 ) 2 -((h 0,n i ) 2 + kq 1,n i+1 |q 1,n i+1 |(h 0,n i+1 ) -η , β R,-1/2 i+ 1 2 = - k∆t 2 α 1/2 i+ 1 2 q 1,n i+1 |q 1,n i+1 |(h 0,n i+1 ) -η , β R,-1 i+ 1 2 = - k∆t 2 α 0 i+ 1 2 q 1,n i+1 |q 1,n i+1 |(h 0,n i+1 ) -η .
Adopting this new notation, we easily rewrite (34) as follows:

q n+1 i = q n i +β L,ε i+ 1 2 + β R,ε i-1 2 + ∆t 2ε∆x λ i+ 1 2 α ε i+ 1 2 q n i+1 -q n i -λ i-1 2 α ε i-1 2 q n i -q n i-1 - ∆t 2ε∆x α ε i+ 1 2 (q n i+1 ) 2 h n i+1 - (q n i ) 2 h n i + α ε i-1 2 (q n i ) 2 h n i - (q n i-1 ) 2 h n i-1
.

(42) Now, let us consider the limit of (42) as ε goes to zero. Since the behavior of α i+ 1 2 and β K,ε i+ 1 2 depends on the values of the water height, two cases must be distinguished to exhibit the limit scheme. First, let us assume that h 0,n i = h 0,n i-1 and/or h 0,n i = h 0,n i+1 , then the limit of (42) as ε goes to zero reads

q 1,n i |q 1,n i |(h 0,n i ) -η = 0. (43)
Next, if we have h 0,n i ̸ = h 0,n i-1 and h 0,n i ̸ = h 0,n i+1 then, in the limit of ε to zero, (42) now tends to

α 1 i+ 1 2 α 1 i+ 1 2 + α 1 i-1 2 g 2 (h 0,n i+1 ) 2 -g 2 (h 0,n i ) 2 ∆x + α 1 i-1 2 α 1 i+ 1 2 + α 1 i-1 2 g 2 (h 0,n i ) 2 -g 2 (h 0,n i-1 ) 2 ∆x = -kq 1,n i |q 1,n i |(h 0,n i ) -η . ( 44
)
To achieve the establishment of (i), we have to show that the above limit scheme, given by ( 43) or (44), is consistent with the local equilibrium [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF]. To address such an issue, we introduce two smooth enough functions, h 0 (x) and q 1 (x), such that (h 0,n i , q 1,n i ) = (h 0 (x i ), q 1 (x i )) and (h 0,n i±1 , q 1,n i±1 ) = (h 0 (x i ± ∆x), q 1 (x i ± ∆x)). Now, in order to deal with the limit scheme, (43) or (44), we have to distinguish the case h 0 (x i ) = h 0 (x i -∆x) and/or h 0 (x i ) = h 0 (x i + ∆x) and the case h 0 (x i ) ̸ = h 0 (x i ± ∆x). Since the expected consistency comes in the limit of ∆x to zero, the case h 0 (x i ) = h 0 (x i -∆x) and/or h 0 (x i ) = h 0 (x i + ∆x), for all ∆x small enough, necessarily implies h 0 (x) equal to a constant in a neighborhood of x i . Then, since h 0 (x i ) > 0, the limit (43) gives q 1 (x i ) = 0, which coincides with the expected local equilibrium [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF] with a constant water height. Next, let us assume that h 0 (x) is not given by a constant. As a consequence we have h 0 (x i ) ̸ = h 0 (x i ± ∆x) as soon as ∆x is small enough. Then, we easily obtain the following limit

lim ∆x→0 α 1 i+ 1 2 α 1 i+ 1 2 + α 1 i-1 2 g 2 h 0 (x i + ∆x) 2 -g 2 h 0 (x i ) 2 ∆x + α 1 i-1 2 α 1 i+ 1 2 + α 1 i-1 2 g 2 h 0 (x i ) 2 -g 2 h 0 (x i -∆x) 2 ∆x = ∂ x gh 0 (x i ) 2 2 ,
to immediately deduce the required consistency of the ε-asymptotic scheme (44) with the local equilibrium [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF] as the water height is not constant. This achieves the establishment of (i).

Now we turn proving (ii). First, we give the behavior of the discretization (33) for ε in a neighborhood of zero. To address such an issue, let us underline that the expansions (38) now write:

h n i = h 0,n i +O(ε) and q n i =      ε q 1,n i + O(ε) if h 0,n i ̸ = h 0,n i-1 and h 0,n i ̸ = h 0,n i+1 , ε 2 q 2,n i + O(ε) otherwise.
The above expansions combined with the expansion (40) of α ε i+ 1 2 imply:

1 2ε∆x α ε i+ 1 2 q n i+1 -q n i = O(ε), ( 45
)
so that the discrete water height evolution law (33) now reads:

h n+1 i = h n i + ∆tO(ε) + ∆t 2ε∆x λ i+ 1 2 α ε i+ 1 2 h n i+1 -h n i -λ i-1 2 α ε i-1 2 h n i -h n i-1 . ( 46
)
To simplify the notations, let us set

γ ε i+ 1 2 = λ i+ 1 2 α ε i+ 1 2 2ε∆x h n i+1 -h n i ,
such that we reformulate (46) as follows:

h n+1 i = h n i + ∆tO(ε) + ∆t γ ε i+ 1 2 -γ ε i-1 2 . ( 47
)
From (40), the following behavior of γ ε i+ 1 2 is obtained:

γ ε i+ 1 2 =                      λ i+ 1 2 α 1 i+ 1 2 2∆x h 0,n i+1 -h 0,n i + O(ε) if h 0,n i ̸ = h 0,n i+1 , √ ε λ i+ 1 2 α 1/2 i+ 1 2 2∆x h 1,n i+1 -h 1,n i + O(ε) if h 0,n i = h 0,n i+1 and h 1,n i ̸ = h 1,n i+1 , ε λ i+ 1 2 α 0 i+ 1 2 h 2,n i+1 -h 2,n i 2∆x + O(ε 2 ) if h 0,n i = h 0,n i+1 and h 1,n i = h 1,n i+1 .
Thus, as soon as h 0,n i = h 0,n i+1 , the quantity γ ε i+ 1 2 is at least of order √ ε. Then, the above expansion rewrites as follows:

γ ε i+ 1 2 = λ i+ 1 2 α 1 i+ 1 2 2∆x h 0,n i+1 -h 0,n i + O( √ ε). ( 48
)
Moreover, by the definition of α 1

i+ 1 2
, given by (41), the following equality holds :

λ i+ 1 2 α 1 i+ 1 2 2∆x h 0,n i+1 -h 0,n i = 1 2∆x g k sign(h 0,n i+1 -h 0,n i ) (h 0,n i+1 ) η + (h 0,n i ) η (h 0,n i+1 ) 2 -(h 0,n i ) 2 ∆x . ( 49
)
As a consequence, the expansion (48) of γ ε i+ 1 2 can thus be rewritten as follows:

γ ε i+ 1 2 = 1 ∆x sign(h 0,n i+1 -h 0,n i ) (h 0,n i+1 ) η + (h 0,n i ) η 2k g 2 (h 0,n i+1 ) 2 -g 2 (h 0,n i ) 2 ∆x + O( √ ε). (50)
Then, equipped with the above expression, the discretization (33), or equivalently (47), thus admits the following behavior:

h n+1 i = h n i + ∆t ∆x     sign h 0,n i+1 -h 0,n i (h 0,n i+1 ) η + (h 0,n i ) η 2k g 2 (h 0,n i+1 ) 2 -g 2 (h 0,n i ) 2 ∆x -sign h 0,n i -h 0,n i-1 (h 0,n i ) η + (h 0,n i-1 ) η 2k g 2 (h 0,n i ) 2 -g 2 (h 0,n i-1 ) 2 ∆x     + ∆tO( √ ε).
In the limit of ε to zero, we immediately obtain the expected consistency with the diffusion equation (8). The proof is thus achieved.

Numerical illustrations

This section is devoted to the numerical representation of the convergence of the discretization given by ( 33)-(34) towards the discretization given by the limit scheme (23). The preservation of the asymptotic behavior of the discharge given by ( 7) is also highlighted. The two initial conditions introduced in Section 3 and given by ( 28) and (29) are considered. Moreover, the parameter k is set equal to 1, the parameter η is equal to 7/3 and zero-flux boundary conditions are used.

We compare the discretization of the water height given by the scheme (33)-(34) with the straightforward discretization (23) of the diffusive limit. As detailed in the proof of Theorem 4.1, the approximate water height given by (33) coincides with the discretization given by this limit scheme when ε tends to zero. On Figures 11 and12, we display the approximate water height, at time t = 0.01 and for both initial conditions, obtained with the scheme (33)-( 34) for different values of ε, the asymptotic parameter. The dotted line denotes the limit water height obtained with (23). We can clearly observe the preservation of the diffusion equation by the scheme (33)-(34). However, we see that when ε goes to zero, the water height discretizations starts by overestimate the limit, then underestimate it, to finish by converging towards the correct limit discretization. This phenomenon is probably due to the transitory discretizations, whose behavior is unknown. To complete this illustration, we display on Figures 13 and 14 the water height obtained for a fixed value of ε and different number of cells. The expected space limit in these regimes stays unknown. Just as an indication, we display the approximation of the diffusive limit as ε tends to zero. The behavior of the space convergence is clearly better than for the scheme developed in [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF], whose results are displayed in Section 3. However, the space convergence seems difficult to be reached. Concerning Figures 13(a) and 14(a), we remark a zig-zag convergence when N grows. In Figure 13(a), for N between 100 and 200, we note that the local maximum decreases, and then increases for N larger than 400. We suppose that this weird convergence is due to the strongly nonlinear differential operator involved in the limit regime. Next concerning Figures 13(b) and 14(b), once again we expect a zig-zag mesh convergence. Nevertheless, to observe this phenomenon, we need very fine meshes involving extremely small values of the time step, so that converging simulations are computationally too expensive. Moreover, to observe the behavior of the order one of the approximate discharge given by the scheme (34), we consider (q n i /ε) n i where (q n i ) n i is the discharge given by (34). We compare it with the discrete local equilibrium (30) on Figures 15 and16 for a fixed value of N and different values of ε. At the discrepancy with results from Section 3, the convergence toward the discrete local equilibrium is illustrated and the scheme (33)-(34) preserves both the diffusion equation (8) and the local equilibrium [START_REF] Buet | Asymptotic preserving and positive schemes for radiation hydrodynamics[END_REF]. Finally, on Figures 17 and18, we display the order one of the discharge compared to the discrete local equilibrium (30) for a fixed value of ε and different number of cells. The spatial convergence of the scheme (33)-( 34) is, here again, better than the spatial convergence of the scheme (20)-( 21), from [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF], displayed in Section 3. Moreover, we remark the same weird convergence behavior as detailed for Figures 13 and14.

Conclusion

In this article, we developed a scheme for the shallow-water equations with Manning friction preserving the correct behavior in long time and dominant friction regime. The proposed method is a generalization of the perturbed HLL discretization introduced in [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF]. In this previous work [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF], the considered framework is linear, in the sense that the source term is linear and leads to a diffusive limit involving a linear operator of Laplacian type. In the case we are interested in here, the shallow-water equations with Manning friction, the source term is quadratic, leading to a diffusive limit involving a nonlinear differential operator of p-Laplacian type. We highlighted the fact that the linear framework could not be directly applied to the considered problem. Indeed, we illustrated the failure of the scheme proposed in [START_REF] Audit | A radiation-hydrodynamics scheme valid from the transport to the diffusion limit[END_REF] to capture correctly both the order zero and one of the solution in the diffusive regime. To overcome this difficulty, we proposed an extension of the generic method introduced in [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF] respecting the natural rescaling of the source term stiffness. This has allowed us to catch consistent discretizations of both the order zero and one of the diffusive limit. This consistency is formally established thanks to Chapman-Enskog expansions of the approximate solutions, and illustrated by several numerical experiments.

Figure 1 :

 1 Figure 1: Structure of the approximate Riemann solver

Figure 2 :

 2 Figure 2: Initial water height.

Figure 3 :

 3 Figure 3: Water height at time t = 0.01 with the continuous initial condition for different values of ε.

Figure 4 :

 4 Figure 4: Water height at time t = 0.01 with the discontinuous initial condition for different values of ε.

Figure 5 :

 5 Figure 5: Water height at time t = 0.01 with the continuous initial condition for different number of cells.

Figure 6 :

 6 Figure 6: Water height at time t = 0.01 with the discontinuous initial condition for different number of cells.

  (a) N = 100 (b) N = 800

Figure 7 :

 7 Figure 7: Order one of the water discharge at time t = 0.01 with the continuous initial condition for different values of ε.

Figure 8 :

 8 Figure 8: Order one of the water discharge at time t = 0.01 with the discontinuous initial condition for different values of ε.

Figure 9 :Figure 10 :

 910 Figure 9: Order one of the water discharge at time t = 0.01 with the continuous initial condition for different number of cells.

1 2 1 2)

 11 are positive and well defined. For the sake of simplicity in the notations, we introduce (β L,ε i+ i∈Z and (β R,ε i+ 1 2

Figure 11 :

 11 Figure 11: Water height at time t = 0.01 with the continuous initial condition for different values of ε.

Figure 12 :

 12 Figure 12: Water height at time t = 0.01 with the discontinuous initial condition for different values of ε.

  (a) ε = 0.01 (b) ε = 0.0001

Figure 13 :

 13 Figure 13: Water height at time t = 0.01 with the continuous initial condition for different number of cells.

  (a) ε = 0.01 (b) ε = 0.0001

Figure 14 :

 14 Figure 14: Water height at time t = 0.01 with the discontinuous initial condition for different number of cells.

  (a) N = 100 (b) N = 800

Figure 15 :

 15 Figure 15: Order one of the water discharge at time t = 0.01 with the continuous initial condition for different values of ε.

Figure 16 :

 16 Figure 16: Order one of the water discharge at time t = 0.01 with the discontinuous initial condition for different values of ε.

Figure 17 :

 17 Figure 17: Order one of the water discharge at time t = 0.01 with the continuous initial condition for different number of cells.

Figure 18 :

 18 Figure 18: Order one of the water discharge at time t = 0.01 with the discontinuous initial condition for different number of cells.

Acknowledgment. The authors are supported by the MoHyCon project (ANR-17-CE40-0027-01), the MUFFIN project (ANR-19-CE46-0004), the Centre Henri Lebesgue (ANR-11-LABX-0020-01) and the Centre de calcul intensif des Pays de la Loire (CCIPL).


Arash Rashidi

Amalya Minasyan

Alexis Cailly

Maher Hamdi

Olivier Redon

Laurent Dussopt

Hussein Yahia

Fast

Fast image deconvolution for enhancement of the resolution in the video rate terahertz imaging

I. INTRODUCTION

ERAHERTZ (THz) imaging has been developing very fast recently. Among various imaging techniques, camera technology has its advantages in terms of speed, bandwidth, and ease of use. I2S TZcam camera is based on CEA-LETI 320240 antenna-coupled bolometer array and has demonstrated state-of-the-art performances for real-time imaging at 0.1-5 THz [START_REF]i2s -TZcam camera[END_REF][START_REF] Meilhan | Performance improvements of THz imagers based on uncooled antenna-coupled bolometer[END_REF]. However, the resolution of the images is often limited by the characteristics of the material traversed as well as by the wavelength of the THz illumination.

Many research works are devoted to the resolution enhancement in the far field of THz imaging through deconvolution techniques [START_REF] Ning | Resolution enhancement in terahertz imaging via deconvolution[END_REF], which is a widely studied and challenging research area in the field of image processing.

In this work, we propose a new fast image deconvolution algorithm using approximated proximal operators. It combines the idea of variable splitting [START_REF] Teodoro | Image restoration and reconstruction using variable splitting and classadapted image priors[END_REF] with new sparsity inducing regularizers. For a THz imaging system, the image processing can be mathematically modeled in the matrix formation by reconstructing an image 𝑥 from an indirect blurred and noisy observation 𝑦: 𝑦 = 𝐻𝑥 + 𝑛 (1) where 𝐻 denotes a known PSF in matrix form and 𝑛 additive Gaussian noise with zero mean and known variance 𝜎 2 . A simple and common approach to solve the problem is using regularization methods i.e., by solving the following optimization problem:

𝑎𝑟𝑔 𝑚𝑖𝑛 𝜇 2 ‖𝐻𝑥 -𝑦‖ 2 + 𝜑(𝐷𝑥) (2)
In our case 𝐷, denotes the discrete gradient operator. The first term of the equation measures the closeness of solution to the observation 𝑦, it is known as data fidelity term. The second term 𝜑(𝑢), the regularizer, enforces a prior knowledge about 𝑥 into the solution.

II. PROPOSED METHOD AND APPROACH

In this section, we will briefly explain the technique to minimize the optimization problem of (2). Here we briefly show the implementation of variable splitting. Variable splitting is a very simple procedure that consists of creating a new variable, say 𝑢, to serve as the argument of 𝜑(. ) under the constraint that 𝑢 = 𝐷𝑥 . The idea is to consider the constrained problem:

𝑎𝑟𝑔𝑚𝑖𝑛 𝑥,𝑢 𝜇 2 ‖𝐻𝑥 -𝑦‖ 2 + 𝜑(𝑢) 𝑠. 𝑡 𝑢 = 𝐷𝑥 (3)
Given the constraint, it must be that 𝑢 = 𝐷𝑥, so the solution to equations ( 3) and ( 2) must be the same or very close. In order to enforce the constraint, we will add an additional term to the cost function that penalizes large differences between 𝑢 and 𝐷𝑥 This then results in the unconstrained optimization problem given by:

𝑎𝑟𝑔𝑚𝑖𝑛 𝑥,𝑢 𝜇 2 ‖𝐻𝑥 -𝑦‖ 2 + 𝜑(𝑢) + 𝛽 2 ‖𝐷𝑥 -𝑢‖ 2 (4)
where 𝛽 is a constant that controls the gain of the penalty term. It also enforces the difference between 𝐷𝑥 and 𝑢 to be as small as possible.

The problem formulated in the form of Eq. ( 3) can be solved by an alternating minimization scheme, i.e., solving two subproblems iteratively:

(𝑆1) 𝑎𝑟𝑔𝑚𝑖𝑛 𝑢 𝜑(𝑢) + 𝛽 2 ‖𝐷𝑥 -𝑢‖ 2 (5) (𝑆2) 𝑎𝑟𝑔𝑚𝑖𝑛 𝑥 𝜇 2 ‖𝐻𝑥 -𝑦‖ 2 + 𝛽 2 ‖𝐷𝑥 -𝑢‖ 2 (6)
An advantage of this formulations is that Fast-Fourier Transform can be used to solve [START_REF] Lorenz | Non-convex variational denoising of images: Interpolation between hard and soft wavelet shrinkage[END_REF] to reduce computational complexity.

III. GENERAL FRAMEWORK AND ALGORITHM

Our algorithm allows the use of multiple regularizers as formulated below:

𝜑 𝑞 𝑝 (𝑥) = { 1 𝑞 (1 -(|𝑥| 𝑝 + 1) -𝑞 ) 𝑞 ≠ 0 𝑙𝑜𝑔 𝑒 (|𝑥| 𝑝 + 1) 𝑞 = 0 (7)
Their approximate proximal operators of (7) are provided as shrinkage/thresholding operators [START_REF] Gholami | A general framework for sparsity-based denoising and inversion[END_REF]: operators work under extremely general conditions and they are fast which is suitable in industry with limited computation capability.

𝑝𝑟𝑜𝑥 1 𝛽 𝜑 𝑞 𝑝 (𝑢) = { (1 - 𝑝 𝛽 ( |𝑢| 𝑝-2 (|𝑢| 𝑝 + 1) 𝑞+1 )) + 𝑢 , |𝑢| > 𝜂 1/𝛽 0 , 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (8)
Depending on the value of 𝑝 and 𝑞 , equation ( 7) defines a potential function, each function has its own characteristic and can be used as prior information in second term of equation ( 2). In Fig. 1 we have depicted various forms of 𝜑 𝑞 𝑝 (𝑥).

We now give the overall algorithm using the general framework for the sub-problem (S1). As outlined in Algorithm 1 below, we minimize Eq. ( 4) by solving the 𝑢 and 𝑥 subproblems separately until the algorithm converges. We keep 𝛽 constant (𝛽 = 256) which makes our algorithm simpler and reduces computational complexity. In step 2 of the algorithm, the solution of subproblem (S1) of equations ( 5) is calculated using the approximate proximal operator of equation ( 8). On step 3 the 𝑢 obtained from step 2 is used to calculate the 𝑥 of Sub problem (S2) of equation ( 6). This process is repeated until the algorithm is converged.

IV. RESULTS

We begin by conducting a numerical experiment to quantify the loss in accuracy by using approximated Proximal Operators for image deconvolution using lena and camera man images. The solution of (4) is computed by two methods: the exact solution, as proposed in [START_REF] Lorenz | Non-convex variational denoising of images: Interpolation between hard and soft wavelet shrinkage[END_REF] and its approximation in (8) for some ℒ 𝑝 = 𝜑 -1 𝑝 (𝑥) quasi-norms. We note that the approach of [START_REF] Lorenz | Non-convex variational denoising of images: Interpolation between hard and soft wavelet shrinkage[END_REF] has very high computational complexity in comparison to approximated Proximal Operators, not suitable for practical applications.

Since the difference cannot be assessed visually, In For a visual demonstration of deconvolution, we used various images at different frequencies acquired with TZcam of I2S. The results of deconvolution on THz images at 0.97 THz and 2.5 THz is shown in Fig. 1. The time of the experiment for 10 iterations is 54 milliseconds on C++ (CPU: Intel(R) Core i7 3.2 GHz). A significant enhancement can be visually assessed. In the enhanced images, shapes and lines can be distinguished easily and the noise level is reduced notably.

V. SUMMARY

A fast deconvolution algorithm based on variable splitting technique with family of sparsity inducing regularizes was proposed and used on THz camera acquired images. The proposed method demonstrated notable enhancement of the resolution and the quality of THz images in various frequencies.

Fast

  image deconvolution for enhancement of the resolution in the video rate

Fig. 1 .

 1 Fig. 1. (a) THz image of a 1951 USAF resolution test chart at 0.97 THz acquired with TZcam and x0.25 magnification lens. (b) Deconvolved image with 𝜑 -1 0.7 . (c) Plots of the Line Pairs pixel in grey values of the group-1 element 2 of the chart of original THz image, (d) deconvolved image. (e) Sample industrial application at 2.5THz. (f) Deconvolved image of (e) with 𝜑 -1 1 .

  Rashidi 1,2 , A. Minasyan 2 , A. Cailly 2 , M. Hamdi 3 , O. Redon 3 , L. Dussopt 3 , H. Yahia 1

	1 INRIA Geostat, Talence, 33405, France
	2 I2S, Pessac, 33608, France
	3 CEA Tech Nouvelle-Aquitaine, Pessac, 33600, France

Here (𝑥) + denotes the positive part of 𝑥. As shown by this equation, it is computed by a simple closed-form expression for values bigger than 𝜂 1/𝛽 and set to zero if smaller. Proximal A.

Table 1

 1 we show the values of Mean Squared Error (MSE) for various values of 𝑝. The MSEs considered here are computed: I) Between the approximated solution and the original image. II) Between the exact solution and the original image.

	The closeness of solutions demonstrates that the use of
	approximated	thresholding	operators	is	almost
	indistinguishable from exact solution.