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

The aim of this article is to propose a numerical scheme for the shallow-water 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.

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 [24]. 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,13]). With this scaling, the system (1) now reads: 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 + · · · (4) 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: 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: The second expression of the above system leads to the following definition of q 1 , called the local equilibrium (see [20]): [1,4,22,25] 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 (1) such that, in the diffusive regime, the approximate solution is consistent with the limit problem (6), or equivalently (7) and (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,6,7], AP numerical strategies have been introduced to accurately approximate the asymptotic diffusive regime issuing from the radiative transfer models. We also refer to [21] where an AP-scheme is derived to approximate the solutions of the isentropic gas dynamics in the Darcy law regime. Next, in [6], a generic formulation of AP-schemes is proposed by a suitable extension of the well-known HLL scheme [18]. 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 [13], 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 [6] 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 [13], the expected regime verified by the discharge, according to (8), is no longer preserved. In this sense, the numerical scheme proposed in [13] 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 (7) and (8).
The article 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 article. Afterwards, in Section 3, the scheme given in [13] for the system under concern (1) is recalled. The failure of this extension of the scheme presented in [6] 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 (7) 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 [6], 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 (1), 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 (7), 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 [6]. To address such an issue, we consider hyperbolic systems with source terms as follows: 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 For the sake of simplicity, we assume k to be a constant and we refer to [6] for more complex interactions.
Next, we introduce a free parameter, denoted k, such that the system (9) equivalently reformulates as follows: with 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 [10]. 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 [6], the required updated states are obtained by adopting a Godunov-type scheme [15,16,18]. 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 [6], the approximate Riemann solver of interest reads as follows: where the parameter LR is defined by 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 [18], defined by Structure of the approximate Riemann solver. and the following notation is used in the definition (13) of the Riemann solver: . We underline that, for the specific symmetric 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 solverW( for all i ∈ Z. Next, we define W n+1 i as follows: x.
With clear notations stated at each interface x i+ 1 2 , after straightforward computations we obtain: Applying the above generic discretization to the shallow-water equations with Manning friction (1), we obtain where the wave velocities ( i+ 1 2 ) i∈Z are fixed as follows: The coefficient is still defined by (14), 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 [13] and we show that the expected AP-property is not fully satisfied. We adopt the scheme (16) and (17) with the wave velocities (18). According to [13], in order to recover the asymptotic diffusion equation (8), the parameters (k i+ 1 2 ) i∈Z are given by 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 : 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) and (17) now reads as follows: where , defined by (14) is also rescaled to get: 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.
while the second equation (21) reduces to a discretization of 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) and (21) in the sense that the diffusive limit (8) is preserved. This is the meaning of asymptotic-preserving adopted in [13].
However, here, we also study the behavior of the second equation of the scheme (21) in the diffusive regime, and the result given by (24) does not coincide with a discretization of the expected local equilibrium (7). Moreover, this new statement establishes a constant behavior of (h n i ) i∈Z as goes to zero, which 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 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 Using the above expressions, in the limit of to zero, the second equation of the scheme (21) yields to Since h 0,n i > 0, we immediately get q 0,n i = 0. Now, once again with (25), the scheme (20) is written as follows in the limit of to zero: Using the definition of i+ 1 2 given by (22) and the definition of k i+ 1 2 given by (19), we have: 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 (25), when tends to zero, the relation (21) reduces as follows: Moreover, since (27) holds, the above relation can be rewritten Now, we show that the above discretization is consistent with (24). 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: The above equality can thus be rewritten as follows: 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 (24). The proof is thus achieved. ▪ From now on, let us illustrate the inappropriate behavior of the discretizations given by the scheme (20) and (21). 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: • Discontinuous initial condition: 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 and 4, 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) and (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) and (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) and (21) developed in [13] cannot be observed, which is in agreement with Theorem 3.1.
To complete this illustration, we display on Figures 5 and 6 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) and (21), we consider (q n i ∕ ) n i where (q n i ) n i is the discharge given by (21). We compare it with the following discretization of the local equilibrium (7): where (h n i ) n i is given by the discrete diffusion equation (23). This comparison is displayed on Figures 7  and 8 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) and (21) is not asymptotic preserving in the sense that it does not preserve both the diffusive limit (8) and the local equilibrium (7).        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) and (21) fails capturing the expected asymptotic behavior.

AN IMPROVED AP-SCHEME
A generic method is now developed to design asymptotic preserving schemes for the hyperbolic problem (11) with a diffusive regime governed by a rescaling of t and k of the following form: where m ∈ N * . Let us emphasize that m is restricted to 1 in [6] 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: where R is now given by 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)- (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 : Moreover, the parameter is chosen as follows: 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 (18). As mentioned previously, the only modification in the scheme (20)- (22) is the substitution of k + k LR by (k + k LR ). This time, the natural rescaling defined by can then be adopted and fixed to to obtain the following rescaled scheme: where ( ) i∈Z is still given by (22) The discrepancy between the above scheme (33) and (34) and the scheme (20) and (21), developed in [13], 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 (21). 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) and (34). Theorem 4.1. We consider the scheme (33), (34), and (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 (7). (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: 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: The discrete discharge evolution law (34) can then be rewritten as follows: 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), can be rewritten as follows: Using expansions (35), from the above expression we deduce the following limit: Even if the limit of ∕ 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: 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: Moreover, the following development now holds: Now, it is necessary to distinguish the expansion of ( ) i∈Z with respect to the behavior of the water height as follows: We notice that the quantities 1 ) i∈Z defined as follows: From (38)-(40), K, i+ 1 2 admits the following expansion: where we have set Adopting this new notation, we easily rewrite (34) as follows: 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 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 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 (7). 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 (7) 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 to immediately deduce the required consistency of the -asymptotic scheme (44) with the local equilibrium (7) 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: The above expansions combined with the expansion (40) of i+ 1 2 imply: so that the discrete water height evolution law (33) now reads: To simplify the notations, let us set such that we reformulate (46) as follows: From (40), the following behavior of i+ 1 2 is obtained: Thus, as soon as h 0,n i = h 0,n i+1 , the quantity is at least of order √ . Then, the above expansion rewrites as follows: Moreover, by the definition of 1 i+ 1 2 , given by (41), the following equality holds: As a consequence, the expansion (48) of can thus be rewritten as follows: Then, equipped with the above expression, the discretization (33), or equivalently (47), thus admits the following behavior: 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) and (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) and (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 and 12, we display the approximate water height, at time t = 0.01 and for both initial conditions, obtained with the scheme (33) and (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) and (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 [13], whose results are displayed in Section 3. However, the space convergence seems difficult to be reached. Concerning  Figures 13a and 14a, we remark a zig-zag convergence when N grows. In Figure 13a, 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 13b and 14b, 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 and 16 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) and (34) preserves both the diffusion equation (8) and the local equilibrium (7).
Finally, on Figures 17 and 18, 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) and (34) is, here again, better than the spatial convergence of the scheme (20) and (21), from [13], displayed in Section 3. Moreover, we remark the same weird convergence behavior as detailed for Figures 13 and 14.

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 [6]. In this previous work [6], 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 [13] 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 [6] 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.