François Bouchut
email: francois.bouchut@univ-mlv.fr

Sébastien Boyaval
email: sebastien.boyaval@enpc.fr

A NEW MODEL FOR SHALLOW VISCOELASTIC FLUIDS

Keywords: Viscoelastic fluids, Maxwell model, Oldroyd model, Saint Venant model, shallow-water, pseudo-conservative variables, well-balanced scheme

We propose a new reduced model for gravity-driven free-surface flows of shallow viscoelastic fluids. It is obtained by an asymptotic expansion of the upper-convected Maxwell model for viscoelastic fluids. The viscosity is assumed small (of order epsilon, the aspect ratio of the thin layer of fluid), but the relaxation time is kept finite. Additionally to the classical layer depth and velocity in shallow models, our system describes also the evolution of two components of the stress. It has an intrinsic energy equation. The mathematical properties of the model are established, an important feature being the non-convexity of the physically relevant energy with respect to conservative variables, but the convexity with respect to the physically relevant pseudoconservative variables. Numerical illustrations are given, based on a suitable well-balanced finitevolume discretization involving an approximate Riemann solver. Contents 1. Introduction: thin layer approximations of non-Newtonian flows 2. Mathematical setting with the Upper-Convected Maxwell model for viscoelastic fluids 3. Formal derivation of a thin layer approximation 4. The new reduced model and its mathematical properties 5. Finite volume method and numerical results 5.1. Approximate Riemann solver 5.2. Energy inequality 5.3. Numerical fluxes and CFL condition 5.4. Topography treatment 5.5. Numerical results 6. Conclusion 6.1. Physical interpretation from the macroscopic mechanical viewpoint 6.2. Physical interpretation at the microscopic molecular level 6.3. Open questions and perspectives Appendix A.

Introduction: thin layer approximations of non-Newtonian flows

There are many occurences of free-surface non-Newtonian flows over an inclined topography in nature, for instance geophysical flows: mud flows, landslides, debris avalanches . . . . Their mathematical prediction is important, typically for safety reasons in connection with land use planning in the case of geophysical flows. But their modelling is still difficult, as one can conclude from the continuing intense activity in that area (see the reviews [START_REF] Ancey | Plasticity and geophysical flows: A review[END_REF][START_REF] Craster | Dynamics and stability of thin liquid films[END_REF] e.g. plus the numerous references cited therein and below). In this paper, our purpose is (i) to derive a new simple model for a thin layer of viscoelastic non-Newtonian fluid over a given topography at the bottom when the motion is essentially driven by gravity forces and (ii) to numerically investigate the prediction of that simple model in benchmark cases.

Our methodology follows the standard derivation of the Saint-Venant model for gravity-driven shallow water flows, as developped in [START_REF] Gerbeau | Derivation of viscous saint-venant system for laminar shallow water ; numerical validation[END_REF] for instance. For non-Newtonian fluids, there already exist similar projects in the literature. But to our knowledge, they use different models as starting point: power-law and Bingham models [START_REF] Fernández-Nieto | Shallow water equations for non newtonian fluids[END_REF][START_REF] Bresch | Augmented lagrangian method and compressible visco-plastic flows: Applications to shallow dense avalanches[END_REF], or a kinetic model for microscopic FENE dumbbells [START_REF] Narbona-Reina | On a shallow water model for non-newtonian fluids[END_REF] (see also the Section 6.2 for comparison with a kinetic interpretation of our model using Hookean dumbbells). Here, we derive a reduced form of the Upper-Convected Maxwell (UCM) equations, a widely-used differential model for viscoelastic fluids valid in generic geometries, in the specific case of gravity-driven free-surface thin-layer flows over an inclined topography. In particular, the influence of each term in the equations is compared with the aspect ratio (1.1) h/L ≈ ǫ ≪ 1 between the layer depth h and its longitudinal characteristic length L as a function of a small parameter ǫ.

The reduced model obtained is computationally much less expensive to solve numerically than the full UCM model with an unknown free surface (compare for instance with numerical simulations in [START_REF] Pasquali | Free surface flows of polymer solutions with models based on the conformation tensor[END_REF] of a full 3D model). So we can easily investigate its predictions in a number of test cases, to show the capabilities of the model. We note two important aspects from the mathematical viewpoint to discretize our model. It is endowed with a natural energy law (inherited from the UCM model) but has a non-standard hyperbolic structure (the physically relevant energy is not convex with respect to the conservative variables). These features of our model have important consequences on the numerical simulation. Whereas we can only perform numerical simulations in a formal way (because the non-standard hyperbolic structure does not fit in the usual numerical analysis), we can nevertheless confirm that they are physically meaningful (owing to the natural energy law, satisfied at the discrete level).

Regarding the literature, we would like to make two further comments in order to better situate our work and its originality. On the one hand, numerous models for thin layers of non-Newtonian fluids have already been derived in the physics literature. We are aware of only one reduced version of the UCM model which is very close to ours, see [START_REF] Entov | Dynamics of free jets and films of viscous and rheologically complex fluids[END_REF][START_REF] Vm Entov | The dynamics of thin liquid jets in air[END_REF] and a sketch of that work in [START_REF] Renardy | Mathematical analysis of viscoelastic flows[END_REF]. But the reduced model, obtained with another methodology and with a different perspective (ad-hoc model to investigate the break-up and swell of free jets and thin films rather than asymptotic analysis of general fluid equations), finally applies in different conditions (without gravity and topography). The other models we are aware of, typically obtained either with a different methodology or (sometimes and) in different conditions (like spin coating e.g.), are different, see for instance [START_REF] Spaid | Viscoelastic free surface flows: spin coating and dynamic contact lines[END_REF][START_REF] Kang | Nonlinear elastic instability of gravity-driven flow of a thin viscoelastic film down an inclined plane[END_REF][START_REF] Ro | Viscoelastic free surface flows: thin film hydrodynamics of Hele-Shaw and dip coating flows[END_REF][START_REF] Kalliadasis | Steady free-surface thin film flows over topography[END_REF][START_REF] Xin | On non-newtonian lubrication with the upper convected maxwell model[END_REF]. On the other hand, recent works in the mathematical literature also studied reduced models for thin layers of viscoelastic flows. For instance [START_REF] Bayada | Viscoelastic fluids in a thin domain[END_REF][START_REF] Bayada | Viscoelastic fluids in thin domains: a mathematical proof[END_REF] derive reduced models for the Oldroyd-B (OB) system of equations, where a purely viscous component div(η s D(u)) is added to the stress term in the right-hand side of (2.2) in comparison with the UCM equations. But our project is different in essence from the thin layer models obtained for those viscoelastic flows without free surface and essentially driven by viscosity instead of gravity. Recall that here we focus on gravity-driven shallow regimes, and that is why we consider the UCM model in particular rather than the OB model (the viscosity only plays a minor role here).

In Section 2 below, we recall the UCM model for viscoelastic fluids and some of its properties in the mathematical setting that is adequate to our model reduction. Then our new reduced model is derived in Section 3 under a given set of clear mathematical hypotheses (which we classically cannot embed into an existence theory for solutions to the non-reduced UCM system of equations). Section 4 is devoted to the study of some mathematical properties of our new reduced model. In Section 5, we provide numerical simulations in benchmark situations where shallow viscoelastic flows could be advantageously modelled by our new system of equations. Last, in Section 6, a physical interpretation of situations modelled by our system of equations is given in conclusion, along with threads for next studies.

Mathematical setting with the Upper-Convected Maxwell model for viscoelastic fluids

The evolution for times t ∈ [0, +∞) of the flow of a given portion of some viscoelastic fluid confined in a moving domain D t ⊂ R d (d = 2 or 3) with piecewise smooth boundary ∂D t is governed by the following set of equations, the so-called Upper-Convected Maxwell (UCM) model [START_REF] Barnes | An introduction to rheology[END_REF][START_REF] Bird | Dynamics of polymeric liquids[END_REF][START_REF] Renardy | Mathematical analysis of viscoelastic flows[END_REF]:

(2.1) div u = 0 in D t , (2.2)
∂ t u + (u • ∇)u = -∇p + div τ + f in D t , (2.3)
∂ t τ + (u • ∇)τ = (∇u)τ + τ (∇u) T + 1 λ (η p D(u) -τ ) in D t ,
where:

• u : (t, x) ∈ [0, +∞) × D t → u(t, x) ∈ R d is the velocity of the fluid,

• D(u) : (t, x) ∈ [0, +∞) × D t → D(u)(t, x) ∈ R d×d S , where R d×d S denotes symmetric real d× d matrices, is the rate-of-strain tensor linked to the fluid velocity u through the relation (2.4) D(u) = 1 2 (∇u + ∇u T ).

• p : (t, x) ∈ (0, +∞) × D t → p(t, x) ∈ R is the pressure,

• τ : (t, x) ∈ [0, +∞) × D t → τ (t, x) ∈ R d×d S
is the symmetric extra-stress tensor, • η p , λ > 0 are physical parameters, respectively a viscosity only due to the presence of elastically deformable particles in the fluid, and a relaxation time corresponding to the intrinsic dynamics of the deformable particles, • f : (t, x) ∈ [0, +∞) × D t → f (t, x) ∈ R d is a body force.

Notice that we have assumed the fluid homogeneous (with constant mass density, hence normalized to one). We also refer to the Section 6 where more details about the UCM model are given along with a physical interpretation of our results. From now on, we assume translation symmetry (d = 2), we endow R 2 with a cartesian frame (e x , e z ) such that f ≡ -ge z corresponds to gravity and we assume that D t has the following geometry (in particular, surface folding like in the case of breaking waves is not possible):

(2.5) ∀t ∈ [0, +∞) , x = (x, z) ∈ D t ⇔ x ∈ (0, L), 0 < z -b(x) < h(t, x),
where b(x) is the topography elevation and b(x) + h(t, x) is the free surface elevation of our thin layer of fluid. Note that the width h(t, x) is an unknown of the problem (it is a free boundary problem). We shall denote as a x (respectively a z ) the component in direction e x (resp. e z ) of any vector (that is a rank-1 tensor) variable a, and similarly the components of higher-rank tensors : a xx , a xz , . . . We denote by n : x ∈ (0, L) → n(x) the unit vector of the direction normal to the bottom and inward the fluid:

(2.6) 2 .

n x = -∂ x b 1 + (∂ x b) 2 n z = 1 1 + (∂ x b)
We supply the UCM model with boundary conditions for all t ∈ (0, +∞): pure slip at bottom,

(2.7) u • n = 0, for z = b(x), x ∈ (0, L), (2.8) τ n = ((τ n) • n) n, for z = b(x), x ∈ (0, L),
kinematic condition at the free surface N t + N • u = 0 where (N t , N ) is the time-space normal, i.e.

(2.9)

∂ t h + u x ∂ x (b + h) = u z , for z = b(x) + h(t, x), x ∈ (0, L),
no tension at the free surface, (2.10)

(pI -τ ) • (-∂ x (b + h), 1) = 0, for z = b(x) + h(t, x), x ∈ (0, L),
plus (for example) inflow/outflow boundary conditions or periodicity in x. We insist on (2.8) without friction. Adding a friction term in (2.8) would not yield the same result. Finally, the Cauchy problem is supplied with initial conditions

(2.11) u(0, x) = u 0 (x), τ (0, x) = τ 0 (x), h(0, x) = h 0 (x),
assumed sufficiently smooth for a solution to exist. Note indeed that the existence theory for solutions to the UCM system (2.1-2.2-2.3) is still very limited (see e.g. [START_REF] Joseph | Hyperbolicity and change of type in the flow of viscoelastic fluids[END_REF][START_REF] Joseph | Change of type and loss of evolution in the flow of viscoelastic fluids[END_REF]), like for non-Newtonian flows with a free surface (see e.g. [START_REF] Meur | Existence locale de solutions des équations d'un fluide viscoélastique avec frontière libre[END_REF][START_REF] Le | Well-posedness of surface wave equations above a viscoelastic fluid[END_REF] for the so-called Oldroyd-B model with a viscous term in (2.2)).

Last, we recall some essential features of the UCM model (

∂ t σ + (u • ∇)σ = (∇u)σ + σ(∇u) T + 1 λ (I -σ) in D t .
In addition, the following properties are easily derived following the same steps as in [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF] for the Oldroyd-B model (except for the absence of the dissipative viscous term η s |D(u)| 2 ). First, for physical reasons, σ should take only positive definite values (this is easily deduced when σ is interpreted as the Grammian matrix of stochastic processes, see [START_REF] Lelièvre | Problèmes mathématiques et numériques posés par la simulation d'écoulement de fluides polymériques[END_REF] e.g. and Section 6.2). The initial condition (2.11) should thus be chosen so that σ(t = 0) is positive definite. Provided the system (2.1-2.11) has sufficiently smooth initial conditions and the velocity field u remains sufficiently smooth, σ indeed remains positive definite (see [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF] e.g.; where the viscosity η s plays no role in the proof). Second, the system (2.1-2.11) is endowed with an energy (the physical free energy)

(2.14) F (u, τ ) = Dt 1 2 |u| 2 + η p 4λ tr(σ -ln σ -I) -f • x dx
which, following Reynolds transport formula and [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF], is easily shown to decay as

(2.15) d dt F (u, τ ) = - η p 4λ 2 Dt tr(σ + σ -1 + 2I)dx .

Formal derivation of a thin layer approximation

Our goal is to derive a reduced model approximating (2.1-2.11) in the thin layer regime h ≪ L where L is a characteristic length of the flow. We follow the formal approach of [START_REF] Bouchut | A new model of saint venant and savage-hutter type for gravity driven shallow water flows[END_REF][START_REF] Bouchut | Gravity driven shallow water models for arbitrary topography[END_REF]. Our main assumption is thus h/L = O(ǫ), thereby introducing an adimensional parameter ǫ → 0. In the following, we simply write

(H1) h = O(ǫ),
where one should infer the unit (L) and the limit (ǫ → 0). This corresponds to a rescaling of the space coordinates with an aspect ratio ǫ between the vertical and horizontal dimensions. We shall also implicitly use a characteristic time T over which the physical variables vary significantly. Then, our task can be formulated as: find a set of non-negative integers

I = (I ux , I uz , I p , I τxx , I τxz , I τzz )
such that a closed system of equations for variables (ũ, p, τ ) approximating (u, p, τ ) holds and

(3.1) (u -ũ, p -p, τ -τ ) = O(ǫ I )
is a uniform approximation on D t (where the powers are applied componentwise). Note that in (3.1) and all along the paper, O has to be understood componentwise in the unit corresponding to the variable. In particular, the previous assumption (3.1) more explicitly signifies

T L (u -ũ), T 2 L 2 (p -p), T 2 L 2 (τ -τ ) = O(ǫ I ) .
We proceed heuristically, increasing little by little the degree of our assumptions on I. Hopefully, the reduced model found that way corresponds to a physically meaningful regime. We recall that another viewpoint is to find a closed system of equations for depth-averages of the main variables of the system (2.1-2.11). The link with our approach is as follows. The conservation of mass for an incompressible, inviscid fluid governed by (2.1) within a control volume governed by the evolution of the free-surface height as given by the kinematic boundary condition (2.9) and the boundary condition (2.7) at the bottom reads as an evolution equation for the free-surface height h (using the Leibniz rule) where the depth-averaged velocity profile ũx :=

1 h b+h b u x dz enters, (3.2) ∀t, x ∈ [0, +∞) × (0, L) 0 = b+h b (∂ x u x + ∂ z u z )dz = ∂ t h + ∂ x b+h b u x dz = ∂ t h + ∂ x (hũ x ) .
The challenge in the derivation of a reduced model for thin layers is then to find a closure for the evolution of ũx in terms of the variables (ũ, p, τ ). In particular, depth-averaging the equation for u x with the boundary conditions (2.9-2.8-2.10) gives, using again the Leibniz rule,

(3.3) ∂ t b+h b u x dz + ∂ x b+h b u 2 x + p -τ xx dz = [(τ xx -p)∂ x b -τ xz ] | b ,
showing that a typical problem is to write an approximation for b+h b u 2

x and for the source term in the right-hand-side of (3.3) in terms of (ũ, p, τ ).

We now give the detailed system of equations (2.1-2.11) in the 2-d geometry of interest:

∂ x u x + ∂ z u z = 0, (3.4a) ∂ t u x + u x ∂ x u x + u z ∂ z u x = -∂ x p + ∂ x τ xx + ∂ z τ xz , (3.4b) ∂ t u z + u x ∂ x u z + u z ∂ z u z = -∂ z p + ∂ x τ xz + ∂ z τ zz -g, (3.4c) ∂ t τ xx + u x ∂ x τ xx + u z ∂ z τ xx = (2∂ x u x )τ xx + (2∂ z u x )τ xz + η p λ ∂ x u x - 1 λ τ xx , (3.4d) ∂ t τ zz + u x ∂ x τ zz + u z ∂ z τ zz = (2∂ x u z )τ xz + (2∂ z u z )τ zz + η p λ ∂ z u z - 1 λ τ zz , (3.4e) ∂ t τ xz + u x ∂ x τ xz + u z ∂ z τ xz = (∂ x u z )τ xx + (∂ z u x )τ zz + η p 2λ (∂ z u x + ∂ x u z ) - 1 λ τ xz , (3.4f)
where we have used (3.4a) to simplify (3.4f). The boundary conditions (2.7), (2.8) and (2.10) write:

u z = (∂ x b)u x at z = b , (3.5a) -(∂ x b)τ xx + τ xz = -∂ x b -(∂ x b)τ xz + τ zz at z = b , (3.5b) -∂ x (b + h)(p -τ xx ) -τ xz = 0 at z = b + h , (3.5c) ∂ x (b + h)τ xz + (p -τ zz ) = 0 at z = b + h, (3.5d)
while the kinematic condition (2.9), following (3.2), writes

(3.6) ∂ t h + ∂ x b+h b u x dz = 0.
We first simplify the derivation of a thin layer regime by assuming that the tangent of the angle between n and e z is uniformly small (H2) ∂ x b = O(ǫ) as ǫ → 0 , hence only smooth topographies with small slopes are treated here. This restriction could probably be alleviated following the ideas exposed in [START_REF] Bouchut | Gravity driven shallow water models for arbitrary topography[END_REF], though at the price of complications that seem unnecessary for a first presentation of our reduced model. On the contrary, the following assumptions are essential:

(H3) η p = O(ǫ), λ = O(1).
(As explained previously, we recall that the assumptions (H3) hold in the unit of the variable, which is here L 2 /T and T respectively.) As usual in Saint Venant models for avalanche flows, we are looking for solutions without small scale in t and x (thus only with scales T and L), but with scale of order ǫ in z (in fact, ǫL), which can be written formally as

(3.7) ∂ t = O(1), ∂ x = O(1), ∂ z = O(1/ǫ)
in the respective units 1/T, 1/L, 1/L. From now on, for the sake of simplicity, we shall not write explicitly the units as functions of T and L wherever they come into play.

We are looking for solutions with bounded velocity u with bounded gradient ∇u. Thus according to (3.7) and to (3.5a), we are led to the following assumptions on the orders of magnitude (H4)

u x = O(1), u z = O(ǫ), ∂ z u x = O(1),
as ǫ → 0.

(A typical profile for u x reads A(t/T, x/L) + zB(t/T, x/L), with any dimensional functions A and B of the adimensional variables t/T and x/L.) According to (2.3), a typical value for τ is η p D(u). Thus we assume accordingly that (H5) τ = O(ǫ) as ǫ → 0.

We deduce from above that there exists some function u 0 x (t, x) depending only on (t, x) such that

(3.8) u x (t, x, z) = u 0 x (t, x) + O(ǫ).
Then, following the classical procedure [START_REF] Gerbeau | Derivation of viscous saint-venant system for laminar shallow water ; numerical validation[END_REF][START_REF] Bouchut | A new model of saint venant and savage-hutter type for gravity driven shallow water flows[END_REF][START_REF] Bouchut | Gravity driven shallow water models for arbitrary topography[END_REF][START_REF] Marche | Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects[END_REF], we find the following successive implications.

i) From the equation (3.4c) on the vertical velocity u z , we get by neglecting terms in O(ǫ)

(3.9) ∂ z p = ∂ z τ zz -g + O(ǫ) .
Hence ∂ z p = O(1), and the boundary condition (3.5d) gives that p = O(ǫ), indeed (3.10)

p = τ zz + g(b + h -z) + O(ǫ 2 ) .
ii) Next, from the equation (3.4b) on the horizontal velocity u x we get

(3.11) ∂ t u 0 x + u 0 x ∂ x u 0 x = ∂ z τ xz + O(ǫ) .
The

τ xz = (∂ t u 0 x + u 0 x ∂ x u 0 x )(z -b) + O(ǫ 2 ) .
In addition the boundary condition (3.5c) implies that τ xz | z=b+h = O(ǫ 2 ). We conclude therefore that (3.13)

∂ t u 0 x + u 0 x ∂ x u 0 x = O(ǫ), τ xz = O(ǫ 2 ).
iii) The previous result combined with the equation (3.4f) on

τ xz implies ∂ z u x = O(ǫ), hence (3.14) u x (t, x, z) = u 0 x (t, x) + O(ǫ 2 ) .
This "motion by slices" property is stronger than the original one (3.8). iv) Using (3.14) and (3.10) in (3.4b) improves (3.11) to

(3.15) ∂ t u 0 x + u 0 x ∂ x u 0 x = ∂ x (τ xx -τ zz -g(b + h)) + ∂ z τ xz + O(ǫ 2 ) ,
which gives, with the boundary condition (3.5b)

[τ xz -∂ x b(τ xx -τ zz )] | z=b = O(ǫ 3 ), (3.16
)

τ xz = [∂ x b(τ xx -τ zz )] | z=b - z b ∂ x (τ xx -τ zz ) dz + ∂ t u 0 x + u 0 x ∂ x u 0 x + g∂ x (b + h) (z -b) + O(ǫ 3 ). But according to (3.5c) combined with (3.10), one has [τ xz -∂ x (b + h)(τ xx -τ zz )] | z=b+h = O(ǫ 3 ), thus with (3.15) (3.17) τ xz = [∂ x (b + h)(τ xx -τ zz )] | z=b+h - z b+h ∂ x (τ xx -τ zz ) dz + ∂ t u 0 x + u 0 x ∂ x u 0 x + g∂ x (b + h) (z -b -h) + O(ǫ 3
). Therefore, the difference of (3.16) and (3.17) yields

(3.18) ∂ t u 0 x + u 0 x ∂ x u 0 x + g∂ x (b + h) h = ∂ x b+h b (τ xx -τ zz ) dz + O(ǫ 3 ) .
We note that τ xz is then given by (3.16) or (3.17) as a function of u 0 x and (τ xx -τ zz ), and the evolution equation (3.18) for u 0

x is exactly the one that one would have obtained after integrating (3.15) in the e z direction and using the boundary conditions (3.5b) and (3.5c) combined with (3.10). It can also be obtained from (3.3).

v) The result (3.14) with the incompressibility condition (3.4a) and the impermeability condition (3.5a) at the bottom also allows to compute the vertical component of the velocity

(3.19) u z = (∂ x b)u x | z=b - z b ∂ x u x dz = (∂ x b)u 0 x -(z -b)∂ x u 0 x + O(ǫ 3 ) ,
which is of course consistent with our hypotheses about u z = O(ǫ). vi) Collecting all the previous results, (3.4d) and (3.4e)

up to O(ǫ 2 ) give (3.20)        ∂ t τ xx + u 0 x ∂ x τ xx + (∂ x b)u 0 x -(z -b)∂ x u 0 x ∂ z τ xx = 2(∂ x u 0 x )τ xx + η p ∂ x u 0 x -τ xx λ + O(ǫ 2 ), ∂ t τ zz + u 0 x ∂ x τ zz + (∂ x b)u 0 x -(z -b)∂ x u 0 x ∂ z τ zz = -2(∂ x u 0 x )τ zz - η p ∂ x u 0 x + τ zz λ + O(ǫ 2 ),
which closes the system of equations for the reduced model. vii) The previous results which give τ xz at order O(ǫ 3 ), that is (3.16) or (3.17), are consistent with the equation (3.4f) for τ xz at order O(ǫ 3 ), from which one could next obtain an approximation for ∂ z u x up to O(ǫ 2 ), that is

(3.21) ∂ t τ xz + u 0 x ∂ x τ xz + ((∂ x b)u 0 x + (z -b)∂ x u 0 x )∂ z τ xz + 1 λ τ xz = ∂ x (∂ x b)u 0 x (z -b) + ∂ x u 0 x τ xx + η p 2λ + ∂ z u x τ zz + η p 2λ
with τ xx and τ zz given up to order O(ǫ 2 ) by (3.20). This procedure fixes the next term in the expansion (3.14). Note in particular that we do not have u x (t, x, z) = u 0 x (t, x) + O(ǫ 3 ) (dependence on the vertical coordinate subsists at order ǫ 2 ).

To sum up, dropping ǫ, we have obtained a closed system of equations

(3.22)                        ∂ t h + ∂ x (hu 0 x ) = 0, ∂ t (hu 0 x ) + ∂ x h(u 0 x ) 2 + g h 2 2 + b+h b (τ zz -τ xx ) dz = -g(∂ x b)h, ∂ t τ xx + u 0 x ∂ x τ xx + (∂ x b)u 0 x -(z -b)∂ x u 0 x ∂ z τ xx = 2(∂ x u 0 x )τ xx + η p λ ∂ x u 0 x - 1 λ τ xx , ∂ t τ zz + u 0 x ∂ x τ zz + (∂ x b)u 0 x -(z -b)∂ x u 0 x ∂ z τ zz = -2(∂ x u 0 x )τ zz - η p λ ∂ x u 0 x - 1 λ τ zz ,
which allows to compute consistently uniform asymptotic approximations of (u x , u z , p, τ xx , τ zz , τ xz ) as variables of order O(ǫ (0,1,1,1,1,2) ), up to errors in O(ǫ (2,3,2,2,2,3) ). These correspond to approximations of (3.4a)-(3.6) up to O(ǫ (2,2,1,2,2,3,3,3,3,2,3) ).

In (3.22), b depends only on x, h and u 0 x depend on (t, x), while τ xx and τ zz depend on (t, x, z). However, observe that the momentum conservation equation invokes only τ zz dz (integrating those for τ xx and τ zz , see (4.4) below in Section 4) and form a closed system with the equations for the momentum and mass conservation. Another equivalent way to derive the same closed system of equations is to assume that τ xx and τ zz are independent of z (at least at first-order in ǫ). In the rest of this paper, we shall mainly be concerned with that simplified system of equations, whose mathematical properties are easier to study.

The new reduced model and its mathematical properties

The reduced model (3.22) is endowed with an energy equation similar to the one for the full UCM model. Obviously, the whole system of equations for τ in the reduced model rewrite with the entries of the conformation tensor σ = I + 2λ ηp τ . However, since it is diagonal at leading order, we consider only the diagonal part (4.1)

σ 0 = σ xx = 1 + 2λ ηp τ xx 0 0 σ zz = 1 + 2λ ηp τ zz .
The two last equations of (3.22) yield (4.2)

     ∂ t σ xx + u 0 x ∂ x σ xx + (∂ x b)u 0 x -(z -b)∂ x u 0 x ∂ z σ xx = 2(∂ x u 0 x )σ xx - 1 λ (σ xx -1), ∂ t σ zz + u 0 x ∂ x σ zz + (∂ x b)u 0 x -(z -b)∂ x u 0 x ∂ z σ zz = -2(∂ x u 0 x )σ zz - 1 λ (σ zz -1).
These equations imply that σ xx and σ zz remain positive if they are initially. Then, we compute (

∂ t + u 0 x ∂ x + (∂ x b)u 0 x -(z -b)∂ x u 0 x ∂ z 1 2 τ xx - η p 4λ ln 1 + 2λ η p τ xx = (∂ x u 0 x )τ xx - 1 η p τ 2 xx σ xx , ∂ t + u 0 x ∂ x + (∂ x b)u 0 x -(z -b)∂ x u 0 x ∂ z 1 2 τ zz - η p 4λ ln 1 + 2λ η p τ zz = -(∂ x u 0 x )τ zz - 1 η p τ 2 zz σ zz . 4.3)
In order to compute the integral of (4.3) with respect to z, we notice the following formula for any function ϕ(t, x, z) (a combination of the Leibniz rule with boundary conditions at z = b and

z = b + h), (4.4)
b+h b ∂ t + u 0 x ∂ x + (∂ x b)u 0 x -(z -b)∂ x u 0 x ∂ z ϕ dz = b+h b ∂ t ϕ + ∂ x (u 0 x ϕ) + ∂ z (∂ x b)u 0 x -(z -b)∂ x u 0 x ϕ dz = ∂ t b+h b ϕ dz -ϕ b+h ∂ t h + ∂ x b+h b u 0 x ϕ dz -(u 0 x ϕ) b+h ∂ x (b + h) + (u 0 x ϕ) b ∂ x b + (∂ x b)u 0 x -h∂ x u 0 x ϕ b+h -(∂ x b)u 0 x ϕ b = ∂ t b+h b ϕ dz + ∂ x u 0 x b+h b ϕ dz .
Therefore, summing up the two equations of (4.3) and integrating in z gives (4.5)

∂ t b+h b η p 4λ tr(σ 0 -ln σ 0 -I) dz + ∂ x u 0 x b+h b η p 4λ tr(σ 0 -ln σ 0 -I) dz = (∂ x u 0 x ) b+h b τ xx -τ zz dz - 1 η p b+h b τ 2 xx σ xx + τ 2 zz σ zz dz.
Moreover, the classical computation of energy for the Saint Venant model gives

(4.6) ∂ t h (u 0 x ) 2 2 + g h 2 2 + gbh + ∂ x h (u 0 x ) 2 2 + gh 2 + gbh u 0 x + u 0 x ∂ x b+h b (τ zz -τ xx ) dz = 0.
Adding up (4.6) and (4.5) yields (4.7)

∂ t h (u 0 x ) 2 2 + g h 2 2 + gbh + η p 4λ b+h b tr(σ 0 -ln σ 0 -I) dz +∂ x h (u 0 x ) 2 2 + gh 2 + gbh + η p 4λ b+h b tr(σ 0 -ln σ 0 -I) dz + η p 2λ b+h b (σ zz -σ xx ) dz u 0 x = - η p 4λ 2 b+h b tr(σ 0 + [σ 0 ] -1 -2I) dz.
Therefore, we get an exact energy identity for solutions to the reduced model (3.22). Note that to discriminate between possibly many discontinuous solutions (generalized solutions in a sense to be defined, see below the discussion on the conservative formulation), we would naturally require an inequality in (4.7) instead of an equality.

In the case of τ xx and τ zz independent of z, everything becomes more explicit. Using the variables σ xx = 1 + 2λ ηp τ xx and σ zz = 1 + 2λ ηp τ zz (also clearly independent of z), the simplified reduced model then writes

(4.8)                      ∂ t h + ∂ x (hu 0 x ) = 0, ∂ t (hu 0 x ) + ∂ x h(u 0 x ) 2 + g h 2 2 + η p 2λ h(σ zz -σ xx ) = -gh∂ x b, ∂ t σ xx + u 0 x ∂ x σ xx -2σ xx ∂ x u 0 x = 1 -σ xx λ , ∂ t σ zz + u 0 x ∂ x σ zz + 2σ zz ∂ x u 0 x = 1 -σ zz λ ,
while the energy inequality becomes (σ 0 is defined in (4.1)) (4.9)

∂ t h (u 0 x ) 2 2 + g h 2 2 + gbh + η p 4λ h tr(σ 0 -ln σ 0 -I) +∂ x h (u 0 x ) 2 2 + gh 2 + gbh + η p 4λ h tr(σ 0 -ln σ 0 -I) + η p 2λ h(σ zz -σ xx ) u 0 x ≤ - η p 4λ 2 h tr(σ 0 + [σ 0 ] -1 -2I).
In (4.8) and (4.9), b is a function of x and h, u 0 x , σ xx , σ zz depend on (t, x), with h ≥ 0, σ xx ≥ 0, σ zz ≥ 0. From now on, we shall only deal with the simplified reduced model (4.8).

The inequality (4.9) (instead of equality) for possibly discontinuous solutions rules out generalized solutions for which the dissipation -already present in our model ! -is physically not enough (see also [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF] where a similar numerical "entropy" condition is used to build stable finite-element schemes for the viscous UCM model, namely the so-called Oldroyd-B model).

Remark 1 (Limit cases). For the system (4.8), two interesting regimes are important to mention. The first is the standard Saint Venant regime, for which one takes η p /λ = 0. It is obtained in the limit η p → 0 for fixed λ (as opposed to the limit λ → ∞ for fixed η p , some kind of "High-Weissenberg limit" [START_REF] Renardy | Mathematical analysis of viscoelastic flows[END_REF] which is problematic, as we will see in the numerical experiments). The second regime is obtained in the "Low-Weissenberg limit" λ → 0, for fixed η p . Assuming (1 -σ xx )/λ and (1 -σ zz )/λ remain bounded, the system rewritten with τ xx and τ zz gives the viscous Saint Venant system (4.10)

     ∂ t h + ∂ x (hu 0 x ) = 0, ∂ t (hu 0 x ) + ∂ x h(u 0 x ) 2 + g h 2 2 -2η p h ∂ x u 0 x = -gh∂ x b,
with the energy inequality (4.11)

∂ t h (u 0 x ) 2 2 + g h 2 2 + gbh + ∂ x h (u 0 x ) 2 2 + gh 2 + gbh -2η p h ∂ x u 0 x u 0 x ≤ -2η p h(∂ x u 0 x ) 2 .
Remark 2 (Steady states). The source terms (1 -σ xx )/λ and (1 -σ zz )/λ in (4.8) are responsible for the right-hand side that dissipates energy in (4.9). This dissipation has the consequence that steady states are possible only if

(4.12) tr(σ 0 + [σ 0 ] -1 -2I) = 0, i.e. τ = 0 ,
which implies that steady solutions to (4.8) identify with the steady solutions at rest to the standard Saint Venant model:

u 0 x = 0, h + b = cst, σ xx = σ zz = 1.
Remark 3 (Conservativity). The reduced model (4.8) is a first-order quasilinear system with source, but not written in conservative form because of the stress equations on σ xx and σ zz . Indeed, one can put them in conservative form as follows,

       ∂ t (σ xx ) -1/2 + ∂ x (σ xx ) -1/2 u 0 x = -σ -3/2 xx 1 -σ xx 2λ , ∂ t (σ zz ) 1/2 + ∂ x (σ zz ) 1/2 u 0 x = σ -1/2 zz 1 -σ zz 2λ . (4.13)
However, these conservative equations do not help since they are physically irrelevant. Moreover, the physical energy of (4.9) is not convex with respect to these conservative variables σ -1/2 xx and σ 1/2 zz . As a matter of fact, one can show that the energy, that is

E = h (u 0 x ) 2 2 + g h 2 2 + gbh + η p 4λ h (σ xx + σ zz -ln(σ xx σ zz ) -2) , (4.14)
cannot be convex with respect to any set of conservative variables of the form

(4.15) h, hu 0 x , h̟ -1 σ -1/2 xx h , hς -1 σ 1/2 zz h ,
where ̟, ς are smooth functions standing for general changes of variables, see Appendix A.

Nevertheless, the system (4.8) can be written in the following canonical form, strongly reminiscent of the gas dynamics system, (4.16)

         ∂ t h + ∂ x (hu 0 x ) = 0, ∂ t (hu 0 x ) + ∂ x h(u 0 x ) 2 + P (h, s) = -gh∂ x b, ∂ t s + u 0 x ∂ x s = 1 λ S(h, s), with (4.17) s = s xx , s zz = σ -1/2 xx h , σ 1/2 zz h , (4.18) S(h, s) = - σ -3/2 xx 2h (1 -σ xx ), σ -1/2 zz 2h (1 -σ zz ) , (4.19) P (h, s) = g h 2 2 + η p 2λ h(σ zz -σ xx ).
One can compute (4.20)

∂P ∂h |s = gh + η p 2λ (σ zz -σ xx + h 2σ zz h + h 2σ xx h ) = gh + η p 2λ (3σ zz + σ xx ) > 0,
from which we conclude that for smooth b, the system (4.16) is hyperbolic with eigenvalues

(4.21) λ 1 = u 0 x -gh + η p 2λ (3σ zz + σ xx ), λ 2 = u 0 x , λ 3 = u 0 x + gh + η p 2λ (3σ zz + σ xx ),
the second having double multiplicity. One can check that λ 2 is linearly degenerate, while λ 1 and λ 3 are genuinely nonlinear (this follows from computations similar to [27, Example 2.4 p.45] and the first line of (5.33)).

From the particular formulation (4.16), one sees that the jump conditions for a 2-contact discontinuity are that u 0

x and P do not jump (as weak 2-Riemann invariants). However, jump conditions across 1-and 3-shocks need to be chosen in order to determine weak discontinuous solutions in a unique way.

A possible choice of jump conditions is, as explained in Remark 3, to take the conservative formulation (4.13) (or equivalently a conservative formulation related to the variables (4.15), leading to the condition that s does not jump through 1-and 3-shocks). This formulation gives unphysical conservations and nonconvex energy (which could produce numerical under/overshoots), and we shall not make this choice.

Our choice of jump conditions will be rather imposed indirectly by numerical considerations, via the choice of a set of pseudo-conservative variables, i.e. variables for which we shall write discrete flux difference equations. Solving nonconservative systems leads in general to convergence to unexpected solutions, as explained in [START_REF] Manuel | Why many theories of shock waves are necessary: convergence error in formally path-consistent schemes[END_REF]. With a pragmatical point of view, we nevertheless choose the pseudo-conservative variables as (4.22) q ≡ (q 1 , q 2 , q 3 , q 4 ) T := h, hu 0 x , hσ xx , hσ zz T .

In other words, we consider the formal system (4.23)

                     ∂ t h + ∂ x (hu 0 x ) = 0, ∂ t (hu 0 x ) + ∂ x h(u 0 x ) 2 + g h 2 2 + η p 2λ h(σ zz -σ xx ) = -gh∂ x b, ∂ t (hσ xx ) + ∂ x (hσ xx u 0 x ) -2hσ xx ∂ x u 0 x = h -hσ xx λ , ∂ t (hσ zz ) + ∂ x (hσ zz u 0 x ) + 2hσ zz ∂ x u 0 x = h -hσ zz λ .
The choice of these pseudo-conservative variables is good for at least two reasons:

• these variables are physically relevant,

• the energy E in (4.14) is convex with respect to them (see Appendix A).

The second point will make it easier to build a discrete scheme that is energy satisfying (in the sense of the energy inequality (4.9)), while preserving the convex (in the variable q) set (4.24)

U = {h ≥ 0, σ xx ≥ 0, σ zz ≥ 0} ,
which is here the physical invariant domain where the energy inequality (4.9) makes sense. Note that our system is of the form considered in [START_REF] Berthon | Why many theories of shock waves are necessary: kinetic relations for non-conservative systems[END_REF] (see also Remark 4).

Let us mention that for the viscous UCM model, namely the Oldroyd-B model, various numerical techniques are proposed in [START_REF] Lozinski | An energy estimate for the Oldroyd-B model: theory and applications[END_REF][START_REF] Lee | New formulations, positivity preserving discretizations and stability analysis for non-Newtonian flow models[END_REF][START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF][START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF] for the preservation of the positive-definiteness of a non-necessarily diagonal tensor σ in the context of finite-element discretizations.

Finite volume method and numerical results

In this section we describe a finite volume approximation of (4.23). The approximation of the full system is achieved by a fractional step approach, discretizing successively the system (4.23) without the relaxation source terms in 1/λ on the right-hand side of the two stress equations, and these relaxation terms alone. The topographic source term h∂ x b is treated by the hydrostatic reconstruction method of [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] in Subsection 5.4. This approach ensures that the whole scheme is well-balanced with respect to the steady states of Remark 2, because the relaxation terms vanish for these solutions.

The integration of relaxation source terms is performed by a time-implicit cell-centered formula. Note that then the scheme is not asymptotic preserving with respect to the viscous Saint Venant asymptotic regime λ → 0 of Remark 1, for this one would need a more complex treatment of these relaxation terms.

Let us now concentrate on the resolution of the system (4.23) without any source, i.e.

(5.1)

             ∂ t h + ∂ x (hu 0 x ) = 0, ∂ t (hu 0 x ) + ∂ x h(u 0 x ) 2 + P = 0, ∂ t (hσ xx ) + ∂ x (hσ xx u 0 x ) -2hσ xx ∂ x u 0 x = 0, ∂ t (hσ zz ) + ∂ x (hσ zz u 0 x ) + 2hσ zz ∂ x u 0 x = 0, with (5.2) P = g h 2 2 + η p 2λ h(σ zz -σ xx ),
and the energy inequality

(5.3) ∂ t h (u 0 x ) 2 2 + g h 2 2 + η p 4λ h σ xx + σ zz -ln(σ xx σ zz ) -2 +∂ x h (u 0 x ) 2 2 + g h 2 2 + η p 4λ h σ xx + σ zz -ln(σ xx σ zz ) -2 + P u 0 x ≤ 0.
A finite volume scheme for the quasilinear system (5.1)-( 5.2) can be classically built following Godunov's approach, considering piecewise constant approximations of q = (h, hu 0 x , hσ xx , hσ zz ), and invoking an approximate Riemann solver at the interface between two cells. 5.1. Approximate Riemann solver. In order to get an approximate Riemann solver for (5.1), we use the standard relaxation approach, as described in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. It naturally handles the energy inequality (5.3), and also preserves the invariant domain (4.24).

Because of the canonical form of (5.1), which is (4.16) without source, i.e.

(5.4)

       ∂ t h + ∂ x (hu 0 x ) = 0, ∂ t (hu 0 x ) + ∂ x h(u 0 x ) 2 + P = 0, ∂ t s + u 0 x ∂ x s = 0, with (5.5) s = s xx , s zz = σ -1/2 xx h , σ 1/2 zz h ,
we have a formal analogy with the system of full gas dynamics equations. Therefore, we follow the usual Suliciu relaxation approach that is described in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. We introduce a new variable π, the relaxed pressure, and a variable c > 0 intended to parametrize the speeds. Then we solve the system (5.6)

                   ∂ t h + ∂ x (hu 0 x ) = 0, ∂ t (hu 0 x ) + ∂ x (h(u 0 x ) 2 + π) = 0, ∂ t (hπ/c 2 ) + ∂ x (hπu 0 x /c 2 + u 0 x ) = 0, ∂ t c + u 0 x ∂ x c = 0, ∂ t s + u 0 x ∂ x s = 0.
This quasilinear system has the property of having a quasi diagonal form (5.7)

                         ∂ t (π + cu 0 x ) + (u 0 x + c/h)∂ x (π + cu 0 x ) - u 0 x h c∂ x c = 0 , ∂ t (π -cu 0 x ) + (u 0 x -c/h)∂ x (π -cu 0 x ) - u 0 x h c∂ x c = 0 , ∂ t 1/h + π/c 2 + u 0 x ∂ x 1/h + π/c 2 = 0 , ∂ t c + u 0 x ∂ x c = 0 , ∂ t s + u 0 x ∂ x s = 0 .
One deduces its eigenvalues, which are u 0 x -c/h, u 0 x + c/h, and u 0 x with multiplicity 4. One checks easily that the system is hyperbolic, with all eigenvalues linearly degenerate. As a consequence, Rankine-Hugoniot conditions are well-defined (the weak Riemann invariants do not jump through the associated discontinuity), and are equivalent to any conservative formulation. We notice that with the relation (5.5) the equation on s in (5.6) can be transformed back to

(5.8) ∂ t (hσ xx ) + ∂ x (hσ xx u 0 x ) -2hσ xx ∂ x u 0 x = 0, ∂ t (hσ zz ) + ∂ x (hσ zz u 0 x ) + 2hσ zz ∂ x u 0 x = 0.
The approximate Riemann solver can be defined as follows, starting from left and right values of h, hu 0

x , hσ xx , hσ zz at an interface :

• Solve the Riemann problem for (5.6) with initial data completed by the relations (5.9)

π l = P (h l , (σ xx ) l , (σ zz ) l ), π r = P (h r , (σ xx ) r , (σ zz ) r ),
and with suitable values of c l and c r that will be discussed below. • Retain in the solution only the variables h, hu 0

x , hσ xx , hσ zz . The result is a vector called R(x/t, q l , q r ).

Note that this approximate Riemann solver R(x/t, q l , q r ) has the property to give the exact solution for an isolated contact discontinuity (i.e. when the initial data is such that u 0

x and P are constant), because in this case the solution to (5.6) is the solution to (5.1) completed with π = P (h, s).

Then, the numerical scheme is defined as follows. We consider a mesh of cells (x i-1/2 , x i+1/2 ), i ∈ Z, of length ∆x i = x i+1/2 -x i-1/2 , discrete times t n with t n+1 -t n = ∆t, and cell values q n i approximating the average of q over the cell i at time t n . We can then define an approximate solution q appr (t, x) for t n ≤ t < t n+1 and x ∈ R by (5.10)

q appr (t, x) = R x -x i+1/2 t -t n , q n i , q n i+1 for x i < x < x i+1 ,
where

x i = (x i-1/2 + x i+1/2 )/2.
This definition is coherent under a half CFL condition, formulated as

(5.11) x/t < - ∆x i 2∆t ⇒ R(x/t, q i , q i+1 ) = q i , x/t > ∆x i+1 2∆t ⇒ R(x/t, q i , q i+1 ) = q i+1 .
The new values at time t n+1 are finally defined by (5.12)

q n+1 i = 1 ∆x i x i+1/2 x i-1/2
q appr (t n+1 -0, x) dx.

Notice that this is only in this averaging procedure that the choice of the pseudo-conservative variable q is involved. We can follow the computations of Section 2.3 in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF], the only difference being that here the system is nonconservative. We deduce that (5.13)

q n+1 i = q n i - ∆t ∆x i F l (q n i , q n i+1 ) -F r (q n i-1 , q n i ) ,
where

(5.14) F l (q l , q r ) = F (q l ) - 0 -∞ R(ξ, q l , q r ) -q l dξ, F r (q l , q r ) = F (q r ) + ∞ 0 R(ξ, q l , q r ) -q r dξ,
and the pseudo-conservative flux is (5.15)

F (q) = (hu 0 x , h(u 0 x ) 2 + P, hσ xx u 0 x , hσ zz u 0 x ).
In (5.15), the two last components are chosen arbitrarily, since anyway the contributions of F in (5.13) cancel out.

Since the two first components of the system (5.6) are conservative, the classical computations in this context give that for these two components, the left and right numerical fluxes of (5.14) are equal and indeed take the value of the flux of (5.6), i.e. hu 0

x and h(u 0 x ) 2 + π, at x/t = 0. We can notice that while solving the relaxation system (5.6), the variables h, s xx and s zz remain positive if they are initially (indeed this is subordinate to the existence of a solution with positive h, which is seen below via explicit formulas and under suitable choice for c l , c r ). By the relation (5.5) this is also the case for σ xx and σ zz . Therefore, the invariant domain U in (4.24) is preserved by the numerical scheme (5.13), this follows from the average formula (5.12) and the fact that U is convex (in the variable q). Remark 4. The above scheme satisfies the maximum principle on the variable s xx , and the minimum principle on the variable s zz . This means that if initially one has s xx ≤ k for some constant k > 0 (respectively s zz ≥ k), then it remains true for all times.

This can be seen by observing that the set where s xx ≤ k (respectively s zz ≥ k) is convex in the variable q, because according to (5.5), (4.22), it can be written as q 1 q 3 ≥ k -2 (respectively k 2 q 3 1 -q 4 ≤ 0). Then, s is just transported during the resolution of (5.6), while the averaging procedure (5.12) preserves the convex sets. Another proof is to write a discrete entropy inequality for an entropy hφ(s xx ), which is convex if 0 ≤ φ ′ ≤ s xx φ ′′ , take for example φ(s xx ) = max(0, s xxk) 2 /2 (respectively for an entropy hφ(s zz ), which is convex if 0 ≤ -φ ′ ≤ 3s zz φ ′′ , take for example

φ(s zz ) = k -1/3 s zz -3 2 s 2/3 zz + 1 2 k 2/3 for s zz ≤ k, φ(s zz ) = 0 for s zz ≥ k).
We shall not write down the details of this alternative proof.

Energy inequality.

We define in a similar way the left and right numerical energy fluxes (5.16) G l (q l , q r ) = G(q l ) -

0 -∞ E R(ξ, q l , q r ) -E q l dξ, G r (q l , q r ) = G(q r ) + ∞ 0 E R(ξ, q l , q r ) -E q r dξ,
where E is the energy of (4.14) without the topographic term gbh, and

(5.17) G = (E + P )u 0

x is the energy flux. We have from [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] that a sufficient condition for the scheme to be energy satisfying is that (5.18) G r (q l , q r ) -G l (q l , q r ) ≤ 0.

When this is satisfied, because of the convexity of E with respect to q one has the discrete energy inequality

(5.19) E(q n+1 i ) -E(q n i ) + ∆t ∆x i G(q n i , q n i+1 ) -G(q n i-1 , q n i ) ≤ 0,
where the numerical energy flux G(q l , q r ) is any function satisfying G r (q l , q r ) ≤ G(q l , q r ) ≤ G l (q l , q r ).

In order to analyze the condition (5.18), let us introduce the internal energy e(q) ≥ 0 by

(5.20) e = g h 2 + η p 4λ σ xx + σ zz -ln(σ xx σ zz ) -2 , so that (5.21) E = h(u 0 x ) 2 /2 + he,
and (∂ h e) |s = P/h 2 . Then, while solving the relaxation system (5.6), we solve simultaneously the equation for a new variable e,

(5.22) ∂ t ( e -π 2 /2c 2 ) + u 0 x ∂ x ( e -π 2 /2c 2 ) = 0,
where e has left and right initial data e(q l ) and e(q r ). The reason for writing (5.22) is that combining it with (5.6) yields

(5.23)

∂ t h(u 0 x ) 2 /2 + h e + ∂ x h(u 0 x ) 2 /2 + h e + π u 0 x = 0. Define now (5.24) G(q l , q r ) = h(u 0 x ) 2 /2 + h e + π u 0 x x/t=0 . Lemma 1.
If for all values of x/t the solution to (5.6), (5.22) satisfies

(5.25) e ≥ e(q), where here q = R(x/t, q l , q r ), then G r (q l , q r ) ≤ G(q l , q r ) ≤ G l (q l , q r ) and the discrete energy inequality (5.19) holds.

Proof. Since (5.23) is a conservative equation, one has

(5.26) G(q l , q r ) = G(q l ) - 0 -∞ h(u 0 x ) 2 /2 + h e (ξ) -E(q l ) dξ = G(q r ) + ∞ 0 h(u 0 x ) 2 /2 + h e (ξ) -E(q r ) dξ.
Therefore, comparing to (5. [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF]), we see that in order to get the result it is enough that for all ξ

(5.27)

E(R(ξ, q l , q r )) ≤ h(u 0 x ) 2 /2 + h e (ξ),
which is (5.25).

In order to go further, we fix the following notation: in the solution to the Riemann problem for (5.6), there are three waves and two intermediate states, denoted respectively by indices l, * and r, * . Then we have the following sufficient subcharacteristic condition (recall that ∂ h P | s is given by (4.20)).

Lemma 2. If c l , c r are chosen such that the heights h ⋆ l , h ⋆ r are positive and satisfy

∀h ∈ [h l , h ⋆ l ] h 2 ∂ h P | s (h, s l ) ≤ c 2 l , ∀h ∈ [h r , h ⋆ r ] h 2 ∂ h P | s (h, s r ) ≤ c 2 r , (5.28)
then (5.25) holds and thus the discrete energy inequality (5.19) is valid.

Proof. The arguments of decomposition in elementary dissipation terms along the waves used in Lemma 2.20 in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] can be checked to apply without modification.

Lemma 3. Denote (5.29) P l = P (h l , s l ), P r = P (h r , s r ),

a l = ∂ h P | s (h l , s l ), a r = ∂ h P | s (h r , s r ),
and define the relaxation speeds c l , c r by (5.30)

c l h l = a l + 2   max 0, u 0 x,l -u 0 x,r + max 0, P r -P l h l a l + h r a r   , c r h r = a r + 2   max 0, u 0 x,l -u 0 x,r + max 0, P l -P r h l a l + h r a r   .
Then the positivity and subcharacteristic conditions of Lemma 2 are satisfied, and the discrete energy inequality (5.19) holds.

Proof. From (4.20) and (5.5) we have (5.31)

∂ h P | s = gh + η p 2λ 3(hs zz ) 2 + 1 (hs xx ) 2 .
Denoting ϕ(h, s) = h ∂ h P | s , we compute (5.32)

∂ h ϕ| s = ∂ h P | s + h 2 ∂ h P | s g + η p 2λ 6hs 2 zz - 2 h 3 s 2 xx = 1 2 ∂ h P | s 2gh + η p λ 3(hs zz ) 2 + 1 (hs xx ) 2 + gh + η p 2λ 6(hs zz ) 2 - 2 (hs xx ) 2 = 1 2 ∂ h P | s 3gh + 6 η p λ (hs zz ) 2 .
Therefore, we deduce that ϕ satisfies (5.33)

∂ h ϕ| s > 0, ϕ(h, s) → ∞ as h → ∞, ∂ h ϕ| s ≤ 2 ∂ h P | s .
Following [Proposition 3.2] [13] with α = 2, we get the result.

Remark 5 (Bounds on the propagation speeds). Lemma 3 is also valid with the formulas of [Proposition 2.18] [12] instead of (5.30). Here we prefer (5.30) because in the context of possibly negative pressure P these formulas ensure the following estimate on the propagation speeds:

(5.34) max c l h l , c r h r ≤ C |u 0 x,l | + |u 0 x,r | + a l + a r ,
with C an absolute constant. This follows from the property that |P | ≤ h∂ h P | s , which is seen on (4.19)-(4.20).

Numerical fluxes and CFL condition.

The Riemann problem for the relaxation system (5.6), (5.22) has to be solved with initial data q l , q r completed with (5.9), the relation (5.5), e l = e(q l ) ≡ e l , e r = e(q r ) ≡ e r , and (5.29), (5.30). The explicit solution is given, according to [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF], by the following formulae. It has three waves speeds

Σ 1 < Σ 2 < Σ 3 , (5.35) Σ 1 = u 0 x,l -c l /h l , Σ 2 = u 0 x, * , Σ 3 = u 0
x,r + c r /h r , and the variables take the value "l" for x/t < Σ 1 , "l*" for Σ 1 < x/t < Σ 2 , "r*" for Σ 2 < x/t < Σ 3 , "r" for Σ 3 < x/t. The "l*" and "r*" values are given by (5.36)

(u 0 x ) * l = (u 0 x ) * r = u 0 x, * = c l u 0 x,l + c r u 0 x,r + π l -π r c l + c r , π * l = π * r = c r π l + c l π r -c l c r (u 0 x,r -u 0 x,l ) c l + c r , 1 h * l = 1 h l + c r (u 0 x,r -u 0 x,l ) + π l -π r c l (c l + c r ) , 1 h * r = 1 h r + c l (u 0 x,r -u 0 x,l ) + π r -π l c r (c l + c r ) , (5.37) c * l = c l , c * r = c r , s * l = s l , s * r = s r , (5.38) σ * xx,l = σ xx,l h l h * l 2 , σ * xx,r = σ xx,r h r h * r 2 , σ * zz,l = σ zz,l h * l h l 2 , σ * zz,r = σ zz,r h * r h r 2 , (5.39) e * l = e l - (π l ) 2 2c 2 l + (π * l ) 2 2c 2 l , e * r = e r - (π r ) 2 2c 2 r + (π * r ) 2 2c 2 r .
Then we need to compute the left/right numerical fluxes (5.14) that are involved in the update formula (5.13). Since the h and hu 0

x components in (5.6) are conservative, classical computations give the associated numerical fluxes, and we have (5.40)

F l = F h , F hu 0 x , F hσxx l , F hσzz l , F r = F h , F hu 0 x , F hσxx r , F hσzz r ,
where the conservative part involves the Riemann solution evaluated at x/t = 0, (5.41)

F h = (hu 0 x ) x/t=0 , F hu 0 x = (h(u 0 x ) 2 + π) x/t=0 .
More explicitly, (5.41) yields that the quantities between parenthese are evaluated at "l" if

Σ 1 ≥ 0, at "l*" if Σ 1 ≤ 0 ≤ Σ 2 , at "r*" if Σ 2 ≤ 0 ≤ Σ 3 , and at "r" if Σ 3 ≤ 0.
As usual there is no ambiguity in the resulting value when equality occurs in these conditions. The numerical energy flux (5.24) involved in (5. [START_REF] Chupin | The FENE viscoelastic model and thin film flows[END_REF]) can be computed in the same way.

We complete these formulas by computing the left/right numerical fluxes for the variables hσ xx , hσ zz from (5.14), (5.42)

F hσxx l = (hσ xx u 0 x ) l + min(0, Σ 1 ) (hσ xx ) * l -(hσ xx ) l + min(0, Σ 2 ) (hσ xx ) * r -(hσ xx ) * l + min(0, Σ 3 ) (hσ xx ) r -(hσ xx ) * r , (5.43)
F hσxx r = (hσ xx u 0 x ) r -max(0, Σ 1 ) (hσ xx ) * l -(hσ xx ) l -max(0, Σ 2 ) (hσ xx ) * r -(hσ xx ) * l -max(0, Σ 3 ) (hσ xx ) r -(hσ xx ) * r ,
the hσ zz fluxes being computed with the same formulas, replacing "xx" by "zz". The maximal propagation speed is then (5.44)

A(q l , q r ) = max(|Σ 1 |, |Σ 2 |, |Σ 3 |),
and the CFL condition (5.11) becomes

(5.45) ∆tA(q i , q i+1 ) ≤ 1 2 min(∆x i , ∆x i+1 ).

Not that with (5.34) and (5.35) we get (5.46)

A(q l , q r ) ≤ C |u 0 x,l | + |u 0 x,r | + a l + a r
with C an absolute constant, bounding the propagation speeds of the approximate Riemann solver whenever the left and right true speeds remain bounded. This property is more general than the possibility of treating data with vacuum considered in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF].

We have obtained finally the following theorem.

Theorem 1. Consider the system (5.1) with the pressure law (5.2), and denote the pseudoconservative variable by q = (h, hu 0 x , hσ xx , hσ zz ). Under the CFL condition (5.45), the scheme (5.13) with the numerical fluxes F l (q l , q r ), F r (q l , q r ) defined above via (5.40), and with the choice of the speeds (5.29), (5.30), satisfies the following properties. (i) It is consistent with (5.1)-(5.2) for smooth solutions, (ii) It keeps the positivity of h, σ xx , σ zz , (iii) It is conservative in the variables h and hu 0

x , (iv) It satisfies the discrete energy inequality (5. [START_REF] Chupin | The FENE viscoelastic model and thin film flows[END_REF]), (v) It satisfies the maximum principle on the variable s xx , and the minimum principle on the variable s zz , (vi) Steady contact discontinuities where u 0 x = 0, P = cst are exactly resolved, (vii) Data with bounded propagation speeds give finite numerical propagation speed. (viii) The numerical viscosity is sharp, in the sense that the propagation speeds Σ i of the approximate Riemann solver tend to the exact propagation speeds when the left and right states q l , q r tend to a common value.

5.4. Topography treatment. Consider now our system (4.23) with topography, but without the relaxation source terms, i.e.

(5.47)

                 ∂ t h + ∂ x (hu 0 x ) = 0, ∂ t (hu 0 x ) + ∂ x h(u 0 x ) 2 + g h 2 2 + η p 2λ h(σ zz -σ xx ) = -gh∂ x b, ∂ t (hσ xx ) + ∂ x (hσ xx u 0 x ) -2hσ xx ∂ x u 0 x = 0, ∂ t (hσ zz ) + ∂ x (hσ zz u 0 x ) + 2hσ zz ∂ x u 0 x = 0.
With respect to the previous sections, the term -gh∂ x b has been put back, where the topography is a given function b(x). For (5.47), the energy inequality (4.9) is modified only by the fact that there is no right-hand side. Thus it can be written (5.48)

∂ t E + ∂ x G ≤ 0, with (5.49) E(q, b) = E(q) + ghb, G(q, b) = G(q) + ghbu 0 x ,
where E and G are given by (5.21), (5.20), (5.17). Recall that the steady states at rest of Remark 2 are defined by (5.50)

u 0 x = 0, h + b = cst, σ xx = σ zz = 1.
Our scheme for (5.47) is written as (5.51)

q n+1 i = q n i - ∆t ∆x i F l (q n i , q n i+1 , ∆b i+1/2 ) -F r (q n i-1 , q n i , ∆b i-1/2 ) ,
where as before q = (h, hu 0 x , hσ xx , hσ zz ), and ∆b i+1/2 = b i+1 -b i . Thus we need to define the left and right numerical fluxes F l (q l , q r , ∆b), F r (q l , q r , ∆b), for all left and right values q l , q r , b l , b r , with ∆b = b r -b l . We use the hydrostatic reconstruction method of [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] (see also [START_REF] Bouchut | A subsonic-well-balanced reconstruction scheme for shallow water flows[END_REF]), and define

(5.52) h ♯ l = h l -(∆b) + + , h ♯ r = h r -(-∆b) + + , (5.53) q ♯ l = h ♯ l , h ♯ l u 0 x,l , h ♯ l σ xx,l , h ♯ l σ zz,l , q ♯ r = h ♯ r , h ♯ r u 0 x,r , h ♯ r σ xx,r , h ♯ r σ zz,r ,
with the notation x + ≡ max(0, x). Note that we use the notation ♯ instead of * in order to avoid confusions with the intermediate states of the Riemann solver of the previous sections. Then the numerical fuxes are defined by (5.54) F l (q l , q r , ∆b) = F l (q ♯ l , q ♯ r ) + 0, g

h 2 l 2 -g h ♯2 l 2 , 0, 0 , F r (q l , q r , ∆b) = F r (q ♯ l , q ♯ r ) + 0, g h 2 r 2 -g h ♯2 r 2 , 0, 0 ,
where F l and F r are the numerical fluxes (5.40) of the problem without topography.

Theorem 2. The scheme (5.51) with the numerical fluxes F l , F r defined by (5.54), (5.52), (5.53) satisfies the following properties.

(i) It is consistent with (5.47) for smooth solutions, (ii) It keeps the positivity of h, σ xx , σ zz under the CFL condition ∆tA(q ♯ l , q ♯ r ) ≤ 1 2 min(∆x l , ∆x r ) with A defined by (5.44), (iii) It is conservative in the variable h, (iv) It satisfies a semi-discrete energy inequality associated to (5.48), (v) It is well-balanced, i.e. preserves the steady states at rest (5.50).

Proof. We ommit the proof of the points (i) to (iii), which follow the proof of Proposition 4.14 in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. For the proof of (v), consider data q l , q r , b l , b r at rest, i.e. satisfying u 0

x,l = u 0 x,r = 0, h l +b l = h r +b r , σ xx,l = σ xx,r = σ zz,l = σ zz,r = 1. Then from (5.52), (5.53) we get q ♯ l = q ♯ r , the common value q ♯ being q r if ∆b ≥ 0, or q l if ∆b ≤ 0. We observe that then F l (q ♯ l , q ♯ r ) = F r (q ♯ l , q ♯ r ) = F (q ♯ ) with F given in (5.15), and that indeed F (q ♯ ) = (0, gh ♯2 /2, 0, 0). The formulas (5.54) yield F l = (0, gh 2 l /2, 0, 0) = F (q l ), F r = (0, gh 2 r /2, 0, 0) = F (q r ). If this is true at all interfaces, (5.51) gives q n+1 i = q n i , which proves the claim. Let us finally prove (iv). First, the scheme without topography satisfies the discrete energy inequality (5.19). According to [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] section 2.2.2, it implies the semi-discrete energy inequality, characterized by (5.55) G(q r ) + E ′ (q r )(F r (q l , q r ) -F (q r )) ≤ G(q l , q r ), G(q l , q r ) ≤ G(q l ) + E ′ (q l )(F l (q l , q r ) -F (q l )), for all values of q l , q r , and where E ′ is the derivative of E with respect to q. Then, for the scheme with topography, the characterization of the semi-discrete energy inequality writes

(5.56) G(q r , b r ) + E ′ (q r , b r )(F r -F (q r )) ≤ G(q l , q r , b l , b r ), G(q l , q r , b l , b r ) ≤ G(q l , b l ) + E ′ (q l , b l )(F l -F (q l )),
where E and G are defined by (5.49), E ′ denotes the derivative of E with respect to q, and G is an unknown consistent numerical entropy flux. Let us choose

(5.57) G(q l , q r , b l , b r ) = G(q ♯ l , q ♯ r ) + F h (q ♯ l , q ♯ r )gb ♯ ,
where F h is the common h-component of F l and F r , and for some b ♯ that is defined below. Then, noticing that E ′ (q, b) = E ′ (q) + gb(1, 0, 0, 0), we can write the desired inequalities (5.56) as

(5.58) G(q r ) + E ′ (q r )(F r -F (q r )) + F h (q ♯ l , q ♯ r )gb r ≤ G(q ♯ l , q ♯ r ) + F h (q ♯ l , q ♯ r )gb ♯ , G(q ♯ l , q ♯ r ) + F h (q ♯ l , q ♯ r )gb ♯ ≤ G(q l ) + E ′ (q l )(F l -F (q l )) + F h (q ♯ l , q ♯ r )gb l .
But using (5.55) evaluated at q ♯ l , q ♯ r and comparing the result with (5.58), we get the sufficient conditions (5.59)

G(q r ) + E ′ (q r )(F r -F (q r )) + F h (q ♯ l , q ♯ r )gb r ≤ G(q ♯ r ) + E ′ (q ♯ r )(F r (q ♯ l , q ♯ r ) -F (q ♯ r )) + F h (q ♯ l , q ♯ r )gb ♯ , G(q ♯ l ) + E ′ (q ♯ l )(F l (q ♯ l , q ♯ r ) -F (q ♯ l )) + F h (q ♯ l , q ♯ r )gb ♯ ≤ G(q l ) + E ′ (q l )(F l -F (q l )) + F h (q ♯ l , q ♯ r )gb l . We compute now (5.60) E ′ (q) = - (u 0 x ) 2 2 + gh - η p 4λ ln(σ xx σ zz ), u 0 x , η p 4λ (1 -1/σ xx ), η p 4λ (1 -1/σ zz ) ,
and writing

(5.61)

F (q) = hu 0 x , h(u 0 x ) 2 + g h 2 2 + η p 2λ h(σ zz -σ xx ), hσ xx u 0 x , hσ zz u 0 x , G(q) = h (u 0 x ) 2 2 + gh 2 + η p 4λ h(σ xx + σ zz -ln(σ xx σ zz ) -2) + η p 2λ h(σ zz -σ xx ) u 0 x ,
we deduce the identity

(5.62) G(q) -E ′ (q)F (q) = -g h 2 2 u 0 x .
Thus the inequality (5.59) simplifies to

(5.63) -g h 2 r 2 u 0 x,r + E ′ (q r )F r + F h (q ♯ l , q ♯ r )gb r ≤ -g h ♯2 r 2 u 0 x,r + E ′ (q ♯ r )F r (q ♯ l , q ♯ r ) + F h (q ♯ l , q ♯ r )gb ♯ , -g h ♯2 l 2 u 0 x,l + E ′ (q ♯ l )F l (q ♯ l , q ♯ r ) + F h (q ♯ l , q ♯ r )gb ♯ ≤ -g h 2 l 2 u 0 x,l + E ′ (q l )F l + F h (q ♯ l , q ♯ r )gb l .
Now, using (5.54) and the fact that

E ′ (q r ) -E ′ (q ♯ r ) = g(h r -h ♯ r ), 0, 0, 0 , E ′ (q l ) -E ′ (q ♯ l ) = g(h l -h ♯ l
), 0, 0, 0 , the desired inequalities (5.63) rewrite

(5.64)

g h r -h ♯ r + b r -b ♯ F h (q ♯ l , q ♯ r ) ≤ 0, g h l -h ♯ l -b ♯ + b l F h (q ♯ l , q ♯ r ) ≥ 0.
We choose now b ♯ = max(b l , b r ), so that (5.64) can be put in the form

(5.65) h r -h ♯ r -(-∆b) + F h (q ♯ l , q ♯ r ) ≤ 0, h l -h ♯ l -(∆b) + F h (q ♯ l , q ♯ r ) ≥ 0.
Finally, taking into account (5.52), we observe that if h l -(∆b) + ≥ 0 then the second line of (5.65) is an identity, otherwise h ♯ l = 0 and the the second inequality of (5.65) holds because F h (0, q ♯ r ) ≤ 0 by the h-nonnegativity of the numerical flux. The same argument is valid for the first inequality of (5.65), which concludes the proof. Remark 6. The maximum principle property on s xx and minimum principle property on s zz , that hold for the solver without topography, are not valid for the above solver with topography, even if it should hold at the continuous level.

Numerical results.

We now illustrate our model by numerical simulations performed with the scheme described above. Note that the model can be considered independently of its derivation and we explore numerical values beyond the physical regime of the Section 3. We denote by H(x) the Heaviside function with jump +1 at x = 0. For all numerical simulations, we chose Neumann conditions at boundary interfaces.

Test case 1. It is a Riemann problem with initial condition (h, hu 0 x , hσ xx , hσ zz )(t = 0) = (3 -2H(x))(1, 0, 1, 1), without source term (b ≡ 0), that can be interpreted as a "dam" break on a wet floor, with polymeric fluid initially at rest everywhere. We first fix η p = λ = 1 and study the convergence of our scheme with respect to the spatial discretization parameter for 50, 100, 200 and 400 points and a constant CF L = 1/2. The results at final time T = .2 are shown in Fig. 1.

Then, using 400 points and a constant CF L = 1/2, we let η p vary as λ = 1 is fixed. The results are shown in Fig. 2 and3. On the one hand, the limit case η p → 0 coincides with the usual shallow-water model, since then the pressure assumes the same values as in a relaxation scheme for the Saint-Venant equations (independent of hσ xx and hσ zz ) while s xx , s zz become passive tracers and their evolution is only one-way coupled -in fact enslaved -to the autonomous dynamics of the Saint-Venant system of equations. On the other hand, the case η p → ∞ is some kind of "rigid limit". As expected from the formulae (4.21) for the eigenvalues of the Jacobian matrix, the left-going rarefaction wave and the right-going shock wave are all the faster as the viscosity η p increases (we hardly see them at T = .2 for η p = 10 +3 in Fig. 2 and3). This is consistent with the physical notion of rigid limit (sound waves are faster in solids than in liquids). On the contrary, the intermediate wave (a right-going contact discontinuity) is all the slower as η p increases, and the jump of h across it is all the larger. This was not obvious to us at first. It could be explained for instance by the fact that c l + c r (at the denominator of the formulae for u 0

x,⋆ and the two intermediate states 1/h ⋆ l , 1/h ⋆ r in (5.36)) becomes very large when η p increases. Notice that the jumps for σ xx and σ zz are directly related to that for h through (5.38) and are consistently small when that for h is large. In any case, the materials on the left of the intermediate wave (where h is higher) is stretched in direction e x and compressed in direction e z , while it is stretched in direction e z and compressed in direction e x on the right, which we can interpret as the manifestation of a close-to-equilibrium stability property (see Section 6).

We also let λ vary as η = 1 is fixed. The results (still at same given time T = .2) are shown in Fig. 4 and5. One clearly sees here the competition between transport and diffusion of viscoelastic effects (the viscoelastic energy is stored in the new variables, that are transported but also diffused). Indeed, the viscoelastic energy is dissipated all the more rapidly as the relaxation time λ is small (the viscous, "Low-Weissenberg" limit). And the waves are all the more smoothed as λ is small (the source terms, which act as diffusive terms, are all the more important). On the other hand, the jump across the contact discontinuity is all the smaller for σ xx and σ zz as λ is small, but all the larger for h. This is coherent with a reasoning similar to the one above when η p only was varied: a smaller relaxation time λ implies slightly faster rarefaction and shock waves because of (4.21), and a larger coefficient c l + c r in the formulae (5.36).

Test case 2. It is a Riemann problem again without source term b ≡ 0 but with vacuum in the initial condition (h, hu 0

x , hσ xx , hσ zz )(t = 0) = (3 -3H(x))(1, 0, 1, 1), which can be interpreted as a "dam" break on a dry floor. The results in Fig. 6 and Fig. 7 at T = .5 show again that small λ and large η p imply a fast right-going rarefaction wave and a slow contact discontinuity, with a large jump for h and small jumps for σ xx , σ zz at contact discontinuity. On the contrary, large λ and small η p imply a slow right-going rarefaction wave and a fast contact discontinuity, with a small jump for h and large jumps for σ xx , σ zz at contact discontinuity. We nevertheless note that the "High-Weissenberg" limit λ → +∞ (recall Remark 1) is more difficult. In particular, although σ zz remains bounded, σ xx seems to become unbounded close to the wet/dry front. Recalling our Remark 4 (see also the comments below about the occurence of vacuum in Test 3), this is

1 1.5 2 2.5 3 -4 -3 -2 -1 0 1 2 3 4
Water height h for different η p

η p = 10 -1 η p = 10 +1 η p = 10 +3 0 0.5 1 1.5 2 2.5 -4 -3 -2 -1 0 1 2 3 4
Water velocity u 0 x for different η p η p = 10 -1 η p = 10 +1 η p = 10 +3 Longitudinal conformation σ xx for different η p

η p = -1 η p = +1 η p = +3 0 0.5 1 1.5 2 2.5 3 3.5 4 -4 -3 -2 -1 0 1 2 3 4
Transversal conformation σ zz for different η p η p = 10 -1 η p = 10 +1 η p = 10 +3 not surprising: there hold minimum principles σ xx ≥ (kh) -2 , σ zz ≥ (kh) 2 in the High-Weissenberg limit λ → ∞, η p /λ = O(1), -without topography -where the relaxation source terms are neglected (except that here initially k = ∞, but k becomes hopefully finite after some time). Hence a possible blow-up at the front, where h → 0, when λ is too large for the diffusive relaxation source terms to compensate for the transport effects. Note that in any case, the materials on the left of the contact-discontinuity wave (where h is non-zero) is stretched in direction e x and compressed in direction e z , except at the vacuum front where a small region with σ closer to the equilibrium I is constantly seen (may be artificially due to our numerical treatment of the vacuum front). While no actual blow-up occurs, this will also be interpreted as a manifestation of a close-to-equilibrium stability property (see Section 6). Test case 3. It is a benchmark introduced in [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF] to test the treatment of topography, see also [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. For x ∈ (0, 25), we compute until T = .25 the evolution from an initial condition (h, hu 0

x , hσ xx , hσ zz )(t = 0, x) = ((10 -b) + , -350 + 700H(x -50/3), (10 -b) + , (10 -b) + ) over a topography b(x) = H(x -25/3) -H(x -25/2). Two rarefaction waves propagate on the left and right sides of the initial velocity singularity at x = 50/3 so that a vacuum is created in between (in the usual Saint-Venant case). In addition, a couple of rarefaction/shock waves is created at each singular point x = 25/3, 25/2 of the topography, but have much smaller amplitudes than the rarefaction waves at x = 50/3. The results in Fig. 8 and 9 are obtained for various η p , λ at a constant η p /λ = 10 -4 . This particular choice was made because then, at the final time, the system is sufficiently close to the Saint-Venant limit η p /λ → 0 so that the pressure is hardly modified compared with the usual Saint-Venant case.

Compared with the usual Saint-Venant case, the double rarefaction wave centered at x = 50/3 cannot create vacuum but at the single point x = 50/3 where the initial velocity is singular. This can be explained as follows. Assuming that the source terms in the stress equations do not influence much the bounds on σ xx , σ zz , in agreement with our Remark 4, the maximum (resp. minimum) principle holds for s xx (resp. s zz ), and there exists a constant k (depending only on the initial conditions since initially h > 0) such that σ xx ≥ (kh) -2 , σ zz ≥ (kh) +2 . But according to the energy bound, one has that (η p /λ) hσ xx dx remains bounded. We deduce that (η p /λ) h -1 dx remains bounded, and therefore h cannot tend to 0 on a whole interval, but can vanish on a single point. We have then σ xx → +∞ at the singular point, here x = 50/3.

Moreover, another vacuum can be created at x = 25/2 (of course still at a single-point because of the previous reasoning). Of course, this new phenomena could just be a numerical artifact, for instance due to the computation of source terms by the hydrostatic reconstruction method. In any case, assuming it is part of the richer phenomenology of our new model compared with the usual Saint-Venant model (see also Section 6 for possible physical interpretations), we observe that the existence of that phenomena depends on η p , λ. It happens only for η P = 10 -4 , λ = 10 +0 and η P = 10 -3 , λ = 10 +1 in our numerical experiments, not for η P = 10 -5 , λ = 10 -1 . But it occurs for larger λ at η P = 10 -5 (when the stress relaxation terms are less important). Indeed, this phenomena seems triggered mainly by high values of λ (the "High-Weissenberg limit", which by the 0 0.5 way requires very small time steps because of the CFL constraint), and still holds at high values of η p (a "rigid limit" that seems to lead to some kind of a "break-up" of the jet here), even at very large η p when no vacuum occurs in between the rarefaction waves. Notice that the latter phenomenon also induces an additional sign change for u 0 x in between the two vacuum points, whose location seems to depend on λ but not on η p , (another indication that the role of viscoelastic dissipation is dominant here). And compared with the usual Saint-Venant case, velocities also assume much greater value (on the left in particular).

Test case 4. In our last test case, we woud like to assess the treament of another type of topography source terms, with creation of dry/wet fronts, by the hydrostatic reconstruction. A usual test case is Thacker's [START_REF] William | Some exact solutions to the nonlinear shallow-water wave equations[END_REF] e.g., which has analytical solutions in the usual Saint-Venant case. But we could not capture interesting phenomena with Thacker's testcase, in particular because the CFL constraint requires the time step to go to 0 very quickly (on short time ranges) due to the creation of dry fronts where h → 0 and σ xx → +∞ (with λ not too small). Note that this does not necessarily mean that this problem does not have global solutions with finite-energy. A timeimplicit scheme (probably hard to build) might be able to compute finite-energy approximations with a non-vanishing time-step.

Here, we consider the test case proposed by Synolakis [START_REF] Costas | The runup of solitary waves[END_REF] to model the runup of solitary waves. This could be used until interesting final times T = 32.5 after the incidental wave has reflected against the shore and created a dry front, see also [START_REF] Marche | Evaluation of well-balanced bore-capturing schemes for 2D wetting and drying processes[END_REF]. We use (h, hu 0

x , hσ xx , hσ zz )(t = 0, x) = ((1. + h 0 (x) -b(x)) + (1, √ gh 0 (x), 1, 1) as initial condition over a topography b(x) = ((x -40.)/19.85) + , x ∈ (0, 100). The pertubation h 0 (x) = α(cosh( √ .75α(x-acosh( 1/.05)/(.75α)))) -2 models a solitary wave as a function of the parameter with α = .019/.1 according to Synolakis semi-analytical theory.

The results in Fig. 10 and 11 show that it is essentially the variations of η p /λ that influence the water height and velocity among all possible variations of η p , λ. And although the first effect of the variations of η p /λ is on the waves celerity, there is no direct match between variations in η p /λ and a time shift as shown in Fig. 10. On the contrary, the variables σ xx , σ zz depend more on λ alone (recall the importance of relaxation source terms, especially in the case where h → 0), at least for such small values of η p /λ as those tested here (sufficiently close to the Saint-Venant regime for the time step not to vanish, even at high values of λ). The smaller λ is, the stronger the dissipation is and thus balances the large stress values that were induced close to the dry front.

Conclusion

We have proposed a new reduced model for the motion of thin layers of viscoelastic fluids (shallow viscoelastic flows) that are described by the upper-convected Maxwell model and driven by the gravity, under a free surface and above a given topography with small slope (like in the standard Saint Venant model for shallow water). More precisely, we have shown formally that for given boundary conditions and under scaling assumptions (H1-5), the solution to the incompressible Euler-UCM system of equations can be approximated by the solutions to the reduced 0 0.5 model (3.22) in some asymptotic regime. We hope that this asymptotic regime in particular is physically meaningful, that and our new model makes sense, possibly beyond the previous asymptotic regime. (That is why we have studied it mathematically and explored it numerically without constraining ourself to a particular regime, as it is usual in such cases.) Observe that in the end we have obtained a flow model whose dynamics is function of the first normal stress difference only, while the shear part of the stress is negligible and computed as an output of the flow evolution. More specifically, the boundary conditions (2.10-2.8) and the flat velocity profile (consequence of the assumed motion by slices) require compatibility conditions on the bulk behaviour of τ xz inside a thin layer. Before looking in future works for other asymptotic regimes, possibly compatible (under different assumptions) with more general kinematics, we would like to conclude here with a better insight of the physical implications of our reduced model. 6.1. Physical interpretation from the macroscopic mechanical viewpoint. We note that the main differences between our model for shallow (Maxwell) viscoelastic flows and the standard Saint Venant model for shallow water is i) a new hydrostatic pressure (3.10), which is function of the (viscoelastic) internal stresses in addition to the water level h, hence ii) a new hydrodynamic force in the momentum balance (in addition to the external gravity force), which is proportional to the normal stress difference τ xx -τ zz , and iii) variable internal stresses τ xx and τ zz , which have their own dynamics corresponding to a viscoelastic mechanical behaviour (with a finite relaxation time λ = O(1) ; such that one recovers the standard viscous mechanical behaviour only in the limit λ → 0). Moreover, in the asymptotic regime where our non-Newtonian model was derived, with a small viscosity parameter η p = O(ǫ), the strain and stress tensors have the scaling (6.1)

∇u = O(1) O(ǫ) O(ǫ) O(1) , τ = O(ǫ) O(ǫ 2 ) O(ǫ 2 ) O(ǫ) .
One essential rheological feature of our reduced model is thus the ratio ǫ between the shear and elongational components of the stress tensor τ and of the strain tensor ∇u. The fact that our model should mainly describe extensional flows, with small shear (of the same small order as the elongational viscosity), seems to be a strong limitation to the applicability of our model in real situations. Of course, one is likely to need another reduced model (in other asymptotic regimes) to describe flows that are not essentially elongational. Note yet that there are situations where physicists arrive at similar conclusions [START_REF] Entov | Dynamics of free jets and films of viscous and rheologically complex fluids[END_REF][START_REF] Vm Entov | The dynamics of thin liquid jets in air[END_REF] and obtain a very similar one-dimensional model with purely elongational stresses for the description of free axisymmetric jets. By the way, a description of free axisymmetric jets is also well achieved by our model since the pure slip boundary conditions (2.7)-(2.8) is equivalent to assuming a cylindrical symmetry around the symmetry line of the jet, and surface tension effects (neglected in our model) can be included using standard modifications of our no-tension boundary condition (2.10).

Moreover, it seems possible to still include non-negligible shear effects in our model through a parabolic correction of the vertical profile like in [START_REF] Gerbeau | Derivation of viscous saint-venant system for laminar shallow water ; numerical validation[END_REF][START_REF] Marche | Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects[END_REF], as well as surface tension and friction effects of order one at the boundaries. 6.2. Physical interpretation at the microscopic molecular level. A microscopic interpretation of our asymptotic regime can also be achieved using a molecular model of the elastic effects (that is, a model at the molecular level from which the UCM is a coarse-grained version at the macroscopic mechanical level). Following [START_REF]Dynamics of polymeric liquids[END_REF], a typical molecular model that accounts for the elasticity of a fluid invokes the transport of elastically deformable Brownian particles diluted in the fluid (which can often be thought of as large massive molecules like polymers). The simplest model of this kind couples, locally in the physical space, a kinetic theory for "dumbbells" (two point-masses connected by an elastic force idealized as a "spring") with the strain of the fluid.

Let us denote X t (x) the connector vector between the two point-masses of a dumbbell modelling a polymer molecule at position x and time t in the fluid. The collection of vector stochastic processes (X t (x)) t∈(0,+∞) parametrized by x ∈ D t is solution to overdamped Langevin equations (6.2)

dX t + (u • ∇)X t dt = (∇u)X t - 2 ζ F(X t ) dt + 2 k B T ζ dB t
for a given field (B t (x)) t∈(0,+∞) of standard Brownian motions (in Itô sense) where ζ is a friction parameter, k B the Boltzmann constant and T the absolute thermodynamical temperature. The UCM equations can be exactly recovered with the specific choice F(X) = HX. Indeed, the extra-stress τ and the conformation tensor σ are given by Kramers relation

(6.3) τ = η p 2λ (m σ -I) with σ(t, x) = 1 H E Xt ⊗ F( Xt ) = X ⊗ X ψ(t, x, ℓ X)d X
where X = X/ℓ is an adimensional version of X, m = Hℓ 2 kB T is a ratio between the elastic potential energy and the heat of the Brownian bath, and η p = 2λnk B T is the molecular interpretation of the polymer viscosity, with n the number density of polymer chains by unit volume (assumed constant as usual for dilute polymer solutions) and λ a characteristic time for dumbbells. One can always choose ℓ such that m = 1. Then, choosing λ = ζ 4H as a relaxation time, Itô formula allows one to exactly recover the UCM system of equations (2.3) when the solvent is assumed inviscid with a velocity field u(t, x) solution to the Euler equations, on noting that the probability density ψ(t, x, ℓ X) satisfies the following Fokker-Planck equation on the unbounded domain X ∈ R 2 : (6.4)

∂ψ ∂t + u • ∇ψ = -div X [(∇u) X - 1 2λ X]ψ + 1 2λ ∆ X ψ .
We note that in [START_REF] Narbona-Reina | On a shallow water model for non-newtonian fluids[END_REF] a reduced model for shallow viscoelastic flows quite similar to ours has already been derived starting from a coupled micro-macro system like (6.4-6.3-2.2), rather than starting from a coarse-grained system at the macroscopic level like the UCM model. The difference between the Hookean micro-macro system above (equivalent in some sense to the UCM model) and the micro-macro system used in [START_REF] Narbona-Reina | On a shallow water model for non-newtonian fluids[END_REF] is the spring force: it corresponds to FENE dumbbells F(X t ) = X t /(1 -|X t | 2 /b) in [START_REF] Narbona-Reina | On a shallow water model for non-newtonian fluids[END_REF]. The FENE force is more physical because it accounts for a finite extension |X t | < b, but contrary to the Hookean dumbbells, it does not have an exact coarsegrained macroscopic equivalent like the UCM model. Yet, if we follow the same procedure as in [START_REF] Narbona-Reina | On a shallow water model for non-newtonian fluids[END_REF] but for Hookean dumbbells, we can hope to derive a reduced micro-macro model whose coarsegrained version is comparable to our new reduced UCM model. Moreover, if the scaling regimes are the same as in [START_REF] Narbona-Reina | On a shallow water model for non-newtonian fluids[END_REF], then our model should also compare to that in [START_REF] Narbona-Reina | On a shallow water model for non-newtonian fluids[END_REF], for an inviscid solvent, in the infinite extensibility limit b → ∞ (where one formally recovers the Hookean dumbbells from FENE dumbbells). Now, observe that the scaling of our new model implies (6.1) ∇u = γ 0 + O(ǫ) where γ 0 = O(1) is a traceless diagonal matrix with entries ∂ x u 0

x , -∂ x u 0 x . Then (6.4) rewrites (6.5) ∂ψ ∂t

+ u • ∇ψ = 1 2λ div X M ∇ X ψ M + O(ǫ) ,
where M (t, x, X) is a weight function proportional to the Maxwellian e -X T (2λγ 0 -I)X . The approximation (6.5) of (6.4) is consistent with our new reduced model provided it yields a consistent approximation for the stress in (6.3): that is, it suffices to show σ xx , σ zz = O(1) and σ xz = O(ǫ) as ǫ → 0. To this aim, let us define an order-one approximation ψ 0 = ψ + O(ǫ) solution to (6.6)

∂ψ 0 ∂t + u 0 • ∇ψ 0 = 1 2λ div X M ∇ X ψ 0 M .
The point is to estimate the terms (6.7) τ 0 = η p 2λ (σ 0 -I), σ 0 = [X ⊗ X] ψ 0 (X)dX . This is not an easy task because of the coupling between ψ 0 and u 0 . Yet it seems reasonable to assume that ψ 0 remains close to the equilibrium solution M/ M for all times (indeed, the Hookean force is derived from an α-convex potential [START_REF] Arnold | Entropies and equilibria of many-particle systems: an essay on recent research[END_REF]), and in particular the Maxwellian ψ 0 ∝ e -(Ax 2 +Bz 2 +2Cxz) has the scaling A = O(1), B = O(1), C = O(ǫ), which implies that σ (and thus τ ) is diagonal at first order. One then obtains a reduced kinetic model which can be exactly coarse-grained into our new reduced UCM model with Itô formula. A macroscopic consequence of the microscopic assumption above is that the reduced model is well-adapted for elongational flows, which is consistent with our macroscopic intepretation of the model. Indeed, everywhere in the macroscopic physical space, one can only expect a balance of internal elastic energy due to stretching or compressing strains in the directions e x and e z , which is the case in elongational flows. (There is not a high probability of permanently sheared dumbbells.) Moreover, if ψ 0 is actually close to the equilibrium M (the particular case A ≈ 2λ∂ x u 0

x -1, B ≈ -2λ∂ x u 0

x -1, C ≈ 0 of our assumption), then, at first-order, the dumbbells are quite uniformly oriented but stretched in one canonical directione x or e z -and necessarily compressed in the orthogonal one (the level-sets of the distribution function are ellipsoidal with principal axes e x and e z at first order). This was indeed observed in those numerical experiments where no blow-up phenomenon seemed to occur.

The microscopic view is in turn a plausible physical explanation at the molecular level of some macroscopic observations. Recall indeed that one-dimensional simple models similar to our model have already been derived in the past to model axisymmetric free jets of elastic liquids [START_REF] Entov | Dynamics of free jets and films of viscous and rheologically complex fluids[END_REF][START_REF] Vm Entov | The dynamics of thin liquid jets in air[END_REF] with a view to explaining the die swell at the end of an extrusion pipe. Now, a miscroscopic interpretation of the die swell is: the elastic energy stored before the die is released after the die.

The dumbbells, mainly compressed in the radial direction e z before the die, stretch just after the die. This may be responsible for an increase of the jet radius (the free-surface of the jet flow equilibrates with the atmospheric pressure) after a characteristic relaxation time linked to λ, hence the so-called delayed die swell.

Finally, we would like to comment on the results obtained in [START_REF] Narbona-Reina | On a shallow water model for non-newtonian fluids[END_REF] with FENE dumbbells. The main differences with our reduced model (which has the micro-macro interpretation explicited above) are: (i) the relaxation time in [START_REF] Narbona-Reina | On a shallow water model for non-newtonian fluids[END_REF] is assumed small λ = O(ǫ), because then it is possible to compute approximate solutions to the Fokker-Planck equation following the Chapman-Enskog procedure of [START_REF] Degond | Viscoelastic fluid models derived from kinetic equations for polymers[END_REF] ; and (ii) the polymer distribution is mainly radial (ψ 0 is only function of |X|), because the authors claim that this suffices to next imply σ xz = O(ǫ) and, as a consequence, a flat profile for the horizontal velocity like in our model. Then, the scaling regimes are not the same, and the radial assumption is too strong to allow one to recover our ellipsoidal probability distribution. So we cannot directly compare our results though they have a similar flavour. 6.3. Open questions and perspectives. First, regarding the interpretation of our model, one might ask whether the present scaling corresponds to a physical situation actually observed for elastic fluids in nature. In particular, the main questionable assumption is of course the pure-slip and no-friction boundary conditions (2.7-2.8) at the bottom (already unrealistic for Newtonian flows, maybe even more unrealistic for non-Newtonian ones). Second, future works on this topic might consider the following directions:

• derive thin-layer reduced models with other equations modelling non-Newtonian flows, which are believed to better model the rheological properties of real materials (constitutive models like Giesekus, PTT, FENE-P, or other molecular models than the FENE dumbbell model used in [START_REF] Chupin | The FENE viscoelastic model and thin film flows[END_REF][START_REF] Narbona-Reina | On a shallow water model for non-newtonian fluids[END_REF]), and in two-dimensional settings (see [START_REF] Bouchut | Gravity driven shallow water models for arbitrary topography[END_REF][START_REF] Marche | Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects[END_REF] for the standard shallow water model); • derive a reduced model closer to real physical situations, possibly in different regimes, or for instance by using a z-dependent velocity profile u x (possibly a multi-layer model) and different boundary conditions than (2.10) and (2.8) (with surface tension and friction at the bottom), which may lead to find physical regimes where τ xz is not negligible; • give a rigorous mathematical meaning and enhance numerical simulations (well-balanced second-order reconstructions) for non-standard systems of equations like the new one presented here.

We note that multi-layer models are also a path to the modelling of some important physical situations, like a thin layer of polymeric fluids on water to forecast the efficiency of oil slick protection plans.

which is obviously convex with respect to W (at fixed b). We conclude that E is convex with respect to q.

τ

  zz dz, which do not depend on z. Now, using Leibniz rule and the boundary conditions, it is possible to get equations for b+h b τ xx dz and b+h b

Figure 1 .

 1 Figure 1. Convergence of the discretized variables hσ xx in Test case 1

Figure 2 .

 2 Figure 2. Variations of the variables h, u 0 x with η p in Test case 1

Figure 3 .

 3 Figure 3. Variations of the variables σ xx , σ zz with η p in Test case 1

Figure 4 .Transversal conformation σ zz for different λ λ = 10 - 1 λ = 10 +0 λ = 10 +1Figure 5 .

 410110105 Figure 4. Variations of the variables h, u 0 x with λ in Test case 1

  Water height h for different η p , λ η p = 10 -1 λ = 10 +1 η p = 10 +0 λ = 10 +0 η p = 10 +1 λ = 10 -1Water velocity u 0x for different η p , λ η p = 10 -1 λ = 10 +1 η p = 10 +0 λ = 10 +0 η p = 10 +1 λ = 10 -1

Figure 6 . 1 0

 61 Figure 6. Variations of the variables h, u 0 x with η p , λ in Test case 2

Figure 7 .

 7 Figure 7. Variations of the variables σ xx , σ zz with η p , λ in Test case 2

Figure 8 .Figure 9 . 2 - 2 Figure 10 . 2 Figure 11 .

 892210211 Figure 8. Variations of the variable h + b, u 0x with η p , λ in Test case 3. We use different labels for the positive (⊕) and negative (⊖) part of the velocity.

  2.1-2.11). Let σ : (t, x) ∈ [0, +∞) × D t → σ(t, x) ∈ R d×d

	S	be the symmetric conformation tensor linked to the symmetric extra-stress
	tensor τ through the relation
	(2.12)		σ = I +	2λ η p	τ ,
	where I denotes the d-dimensional identity tensor. The UCM model can be written using the
	variable σ instead of τ . Indeed,	ηp 2λ div σ replaces div τ in (2.2), and (2.3) should be replaced with
	(2.13)

Appendix A. Convexity of the energy

In order to check the convexity of E in (4.14) with respect to general variables, we use a Lagrange transformation, see for example Lemma 1.4 in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. Thus E is a convex function of

for given smooth invertible functions ̟, ς, if and only if E/h is a convex function of the Lagrangian variables

Let us denote by V i , i = 1, . . . , 4 the entries of the vector V , then the Lagrangian energy writes

Introduce now the notation

Clearly we only need to look at the convexity with respect to (V 1 , V 3 , V 4 ), and the Hessian matrix H of E/h with respect to these variables (at fixed b) is given by

where prime denotes the derivative with respect to the involved V i . Since V 1 can take any positive value at fixed V 3 or V 4 , the positivity of the diagonal terms give the necessary conditions

Then, writing the positivity of the determinant of the 2 × 2 upper left submatrix of H, and looking at the dominant term when V 1 → ∞ yields the necessary condition

Obviously there is no function Ω(V 3 ) satisfying these conditions, and E is never convex with respect to the considered variables. On the contrary, if we choose the physically natural, but non-conservative, variables q = (h, hu 0

x , hσ xx , hσ zz ), then using the Lagrangian variables